Converted document
$\newcommand{\lyxlock}{}$
Section 3.13: Contact Section Up Section 3.13: Contact Section Subsection 3.13.2: Biphasic Contact

### 3.13.1 Sliding Interfaces

A sliding interface can be used to setup a non-penetration constraint between two surfaces. As of version 1.2, three different sliding contact algorithms are available. Although all three are based on the same contact enforcement method, they all differ slightly in their implementation and have been shown to give different performance for different contact scenarios. Each sliding contact implementation is identified by a different type attribute.
sliding-node-on-facet (N2F) This is FEBio's original implementation of sliding contact. It is based on Laursen's contact formulation [39] which poses the contact problem as a nonlinear constrained optimization problem. In FEBio, the Lagrange multipliers that enforce the contact constraints are computed either using a penalty method or the augmented Lagrangian method.
sliding-facet-on-facet (F2F) This implementation is identical to the sliding-node-on-facet implementation but uses a more accurate integration rule: where the former method uses nodal integration, this method uses Gaussian quadrature to integrate the contact equations. This method has been demonstrated to give additional stability and often converges when the former method does not.
sliding-elastic (SE) This sliding contact interface also uses facet-on-facet contact but differs in the linearization of the contact forces, which results in a different contact stiffness matrix compared to the previous two methods. It can be used for frictionless or frictional contact. It may optionally be set to sustain tension to prevent contact surfaces from separating along the direction normal to the interface, while still allowing tangential sliding. This method sometimes performs better than the previous two methods for problems that are dominated by compression. However, the formulation is inherently non-symmetric and therefore will require additional memory and running time. A symmetrized version of this implementation is available (see below), but the symmetric version does not converge as well as the non-symmetric version.
sliding-biphasic (SBP) This method extends the sliding-elastic algorithm to support frictionless biphasic contact (see the next section). This method is described in detail in [10].
sliding-biphasic-solute (SBS) This method is similar to sliding-biphasic. This contact implementation supports biphasic-solute contact (see the next section) [12]. When using biphasic-solute materials, the non-symmetric version must be used.
sliding-multiphasic (SMP) This method is similar to sliding-biphasic-solute. This contact implementation supports multiphasic contact (see the next section) [12]. When using multiphasic materials, the non-symmetric version must be used.
The following table lists the properties that are defined for sliding interfaces. It is important to note that the three different sliding implementations cannot be used interchangeably: not all features are available for each method. The third, fourth and fifth column indicate if a parameter is available for a particular implementation.
 Parameter Description N2F F2F SE SBP SBS SMP Default penalty normal penalty scale factor (1) 1.0 auto_penalty auto-penalty calculation flag (2) 0 update_penalty updated auto-penalty calculation at each step (2) 0 two_pass two-pass flag (3) 0 laugon augmented Lagrangian flag (4) 0 tolerance aug. Lagrangian convergence tolerance (4) 1.0 gaptol tolerance for gap value (4) 0.0 (off) minaug minimum number of augmentations (4) 0 maxaug maximum number of augmentations (4) 10 fric_coeff frictional coefficient (5) 0.0 fric_penalty tangential penalty factor (5) 0.0 ktmult tangential stiffness multiplier (5) 1.0 seg_up maximum number of segment updates (6) 0 (off) smooth_aug smoothed Lagrangian augmentation (7) 0 (off) symmetric_stiffness symmetric stiffness matrix flag (8) 0 search_tol Projection search tolerance (9) 0.01 search_radius search radius (10) 1.0 tension tension flag (11) 0
1. If the augmented Lagrangian flag is turned off (see comment 4), the penalty method is used to enforce the contact constraint. In this method, the contact traction is determined by the gap (i.e. penetration distance) multiplied by the user-defined penalty factor. In the augmented Lagrangian method, the penalty parameter is also used but has a slightly different meaning. In this case, it scales the Lagrange multiplier increment. Due to the different meanings, the user might have to adjust the penalty factor when switching between penalty method and augmented Lagrangian method. In general the penalty method requires a larger penalty factor to reach the same gap than the augmented Lagrangian method. See comment 4 for more information on when to use which method.
2. Choosing a good initial penalty parameter can often be a difficult task, since this parameter depends on material properties as well as on mesh dimensions. For this reason, an algorithm has been implemented in FEBio that attempts to calculate a good initial value for the penalty factor for a particular node/integration piont on the contact interface: Here, is the area of the element the integration point belongs to, is the element volume and is a measure of the elasticity modulus, which is calculated from the elasticity tensor of the element. Although the meaning of depends on the precise material formulation, in general one can regard it as a measure of the small strain Young's modulus of the material, when update_penalty is 0. If update_penalty is set to 1, is recalculated at each iteration to account for the current deformation.
To use this feature, add the following element to the contact section:
<auto_penalty>1</auto_penalty>

When the auto-penalty flag is on, the value of the penalty parameter serves as a scale factor for the automatically-calculated penalty factor.
3. Each sliding interface consists of a primary surface and a secondary surface. The primary surface is the surface over which the contact equations are integrated and on which the contact tractions are calculated. The secondar surface is used to measure the gap function and to define the necessary kinematic quantities such as surface normals and tangents. This approach is usually referred to as the single-pass method. When using the single-pass algorithm, the results can be influenced by the choice of primary and secondary surfaces. It is best to use the most tessellated surface as the primary and the coarsest as the secondary surface. To resolve the bias issue, one can also use a two-pass algorithm for enforcement of the contact constraint. In this case, a single pass is performed first, after which the primary and secondary surfaces are swapped and another pass is performed. When using the two-pass method, the definition of primary and secondary surfaces is arbitrary. In most problems, the single pass is sufficient to enforce contact; with a judicious choice of primary-secondary pair and contact parameters, good results can be obtained. If however, the single pass does not give good answers, for example, when due to the geometry's curvature the gap cannot be small enough with a single pass, the two-pass method can be used, although at the expense of more calculations.
If one of the contacting surfaces is rigid, a slightly different approach is recommended. In this case, it is best to pick the rigid surface as the secondary surface and to use a single pass algorithm. The reason is that the nodal degrees of freedom on the rigid surface are condensed to the rigid degrees of freedom and if the rigid surface is the primary surface, the reaction forces may not propagate correctly to the secondary surface. This is especially true if the rigid degrees of freedom are fixed.
4. In the presence of a sliding interface (and other contact interfaces), FEBio needs to calculate the contact tractions that prevent the two participating surfaces from penetrating. In general these tractions can be found using the method of Lagrange multipliers. However, the direct calculation of these multipliers has several computational disadvantages and therefore FEBio approximates the multipliers using one of two alternative methods: the penalty method and the augmented Lagrangian method. In the former method, the multipliers are approximated by the gap (i.e. penetration distance) scaled by a suitably chosen penalty factor. In many cases, this method is sufficient to get good results. Since the correctness of a contact solution is directly determined by the amount of penetration at the converged state, the user has direct control over the quality of the solution. By increasing the penalty factor, the penetration is reduced. However, in some cases, especially in large compression problems, the penalty factor required to achieve an acceptable amount of penetration has to be so large that it causes numerical instabilities in the non-linear solution algorithm due to ill-conditioning of the stiffness matrix. In these cases, the augmented Lagrangian method might be a better choice. In this method, the multipliers are determined iteratively where, in each iteration, the multiplier's increments are determined with a penalty-like method. The advantage of this method is twofold: due to the iterative nature, the method will work with a smaller penalty factor, and in the limit, the exact Lagrange multipliers can be recovered.
To turn on the augmented Lagrangian method, simply add the following line to the contact section:
<laugon>1</laugon>

To turn off the augmented Lagrangian method, either set the value to 0 or remove the parameter altogether. The convergence tolerance is set as follows:
<tolerance>0.01</tolerance>

With this parameter set, the augmented Lagrangian method will iterate until the relative increment in the multipliers is less than the tolerance. For instance, setting the tolerance parameter to 0.01 implies that the augmented Lagrangian method will converge when there is less than a 1% change in the L2 norm of the augmented Lagrangian multiplier vector between successive augmentations. Alternatively, the user can also specify a tolerance for the gap value. In this case, the iterations will terminate when the gap norm, which is defined as the averaged L2 norm, (, the Macauley Bracket) is less than the user-specified value:
<gaptol>0.001</gaptol>

However, one must be careful when specifying a gap tolerance. First note that the gap tolerance is an absolute value (unlike the tolerance which is a relative value), so this parameter depends on the dimensions of the model. Also, there are cases when a gap tolerance simply cannot be reached due to the geometry of the model in which case the augmentations may never converge.
Note that both convergence parameters may be defined at once. In that case, FEBio will try to satisfy both convergence criteria. On the other hand, setting a value of zero will turn off the convergence criteria. For example, the default value for gaptol is zero, so that FEBio will not check the gap norm by default.
Finally, the minaug and maxaug can be used to set a minimum and maximum number of augmentations. When the maxuag value is reached, FEBio will move on to the next timestep, regardless of whether the force and gap tolerances have been met. When specifying a value for minaug, FEBio will perform at least minaug augmentations.
5. The sliding-node-on-facet and sliding-elastic contact implementations are currently the only contact algorithms that supports friction. For node-on-facet, three parameters control the frictional response: fric_coeff is the material's friction coefficient and its value must be in the range from 0.0 to 1.0; fric_penalty is the penalty factor that regulates the tangential traction forces, much like the penalty parameter regulates the normal traction force component; the parameter ktmult is a scale factor for the tangential stiffness matrix. It is default to 1.0, but it is observed that reducing this value might sometimes improve convergence. For sliding-elastic, only fric_coeff is needed to use frictional contact. Frictional contact is inherently dissipative and thus history-dependent. The solution may vary with the size of time increments, the cycle of loading, or Lagrangian augmentations (when laugon=1).
6. In a contact problem, FEBio calculates the projection of each node of the primary surface onto the secondary surface. As a node slides across the secondary surface, the corresponding surface facet can change. However, in some cases, switching facets is undesirable since it might cause instabilities in the solution process or a state in which the node oscillates continuously between two adjacent facets and thus prevents FEBio from meeting the displacement convergence tolerance. The parameter seg_up allows the user to control the number of segment updates FEBio will perform during each time step. For example, a value of 4 tells FEBio it can do the segment updates during the first four iterations. After that, nodes will not be allowed to switch to new segments. The default value is 0, which means that FEBio will do a segment update each iteration of each timestep.
7. The augmented Lagrangian method may produce a non-smooth contact pressure distribution at the contact surface, even though stresses within the underlying elements may remain relatively smoother. Turning this flag on changes the method of updating the Lagrange multiplier to use the projection of the element stress to the corresponding contact face. This feature only works with some element types (HEX8G8, HEX8G1, TET4G4, TET4G1, PENTA6G6, TET10G1, TET10G4, TET10G8, TET10GL11, TET15G4, TET15G8, TET15G11, TET15G15, PYRA5G8).
8. The sliding-elastic, sliding-biphasic, sliding-biphasic-solute and sliding-multiphasic contact implementations for sliding contact are inherently non-symmetric formulations. Symmetrized versions of these algorithms do not perform as well as the nonsymmetric version so it is recommended to use the latter. The following line in the contact element controls which version of the algorithm is used.
<symmetric_stiffness>0</symmetric_stiffness>

A value of 1 uses the symmetric version, whereas a value of 0 uses the non-symmetric version. A similar line may be included in the Control element to use a non-symmetric global stiffness matrix. In some cases, a non-symmetric contact stiffness may converge well even when the global stiffness matrix is symmetric.
9. The search_tol parameter defines the search tolerance of the algorithm that projects nodes of the primary surface onto facets of the secondary syrface. A node that falls outside an element, but whose distance to the closest element's edge is less than the search tolerance is still considered inside. This can alleviate convergence problems when nodes are projected onto edges of elements and due to numerical error may be projected outside the surface.
10. The search_radius is a scale factor which multiplies the diagonal length of the model's bounding box to produce a dimensional search radius used by the algorithm that projects points onto facets. When the distance between the point and the facet exceeds the dimensional search radius, that projection is ignored. This can alleviate convergence problems when surfaces have multiple folds and the projection produces multiple solutions, only one of which (the closest distance) is valid.
11. The tension flag determines whether the contact interface can sustain tension and compression (tension=1) or only compression (tension=0).
Section 3.13: Contact Section Up Section 3.13: Contact Section Subsection 3.13.2: Biphasic Contact