In COMSOL Multiphysics you can model contact between a group of boundaries in 2D or 3D. There are two algorithms available: an augmented Lagrangian method and a penalty method.

In the augmented Lagrangian method, the system of equations is solved in a segregated way. Augmentation components are introduced for the contact pressure Tn and the components Tti of the friction traction vector Tt. An additional iteration level is added where the usual displacement variables are solved separated from the contact pressure and traction variables. The algorithm repeats this procedure until it fulfills a convergence criterion.

In the following equations F is the deformation gradient matrix. When looking at expressions evaluated on the destination boundaries, the expression map (E) denotes the value of the expression E evaluated at a corresponding source point, and g is the current gap distance between the destination and source boundary.

Both the contact map operator map (E) and the gap distance variable are defined by the contact element elcontact. For each destination point where the operator or gap is evaluated, a corresponding source point is sought by searching in the direction normal to the destination boundary.

It is possible to add an offset value both to the source (osrc), and to the destination (odst). If an offset is used, the gap will be computed with the geometry treated as larger (or smaller, in the case of a negative offset) than what is actually modeled. The offset is given in a direction normal to the boundary. It can vary along the geometry, and may also vary in time (or by parameter value when using the continuation solver).

The effective gap distance dg used in the analysis is thus

In some cases, the discretization caused by the meshing of curved boundaries will give small irregularities in the gap values. You can then choose to adjust the gap by the computed gap value in the initial configuration, ginit. In this case, the effective gap distance is

When including friction in the contact problem, it is important to keep track of not only the gap in the normal direction, but also the slip between the boundaries in the tangential direction. The slip s since the last converged step is defined as

and Xm,old is the value of Xm in the last converged time or parameter step. The coordinates are material (undeformed) coordinates. The slip vector is thus approximated using a backward Euler step. The deformation gradient F contains information about the local rotation and stretching.

Using the special gap distance variable (solid.gap), the penalized contact pressure Tnp is defined on the destination boundary as

pn is the user defined normal penalty factor.

The penalized friction traction Ttp is defined on the destination boundary as:

where Tt,trial is defined as

In Equation 3-92, pt is the user defined friction traction penalty factor. The critical sliding resistance, Tt,crit, is defined as

In Equation 3-93, μ is the friction coefficient, Tcohe is the user defined cohesion sliding resistance, and Tt,max is the user defined maximum friction traction.

In the following equation δ is the variation (represented by the test() operator in COMSOL Multiphysics). The contact interaction gives the following contribution to the virtual work on the destination boundary:

where wcn and wct are contact help variables defined as:

where i is the augmented solver iteration number and γfric is a Boolean variable stating if the parts are in contact or not.

where dg is the effective gap distance, pn is the penalty factor, and p0 is the pressure at zero gap. The latter parameter can be used to reduce the overclosure between the contacting boundaries if an estimate to the contact pressure is known. In the default case, when p0 = 0, there must always be some overclosure (a negative gap) in order for a contact pressure to develop.

In case of friction, the tangential force Tt is computed as

This can be viewed as a tangential spring giving a force proportional to the sliding distance. The penalty factor pt can be identified as the spring constant. The sticking condition is thus replaced by a stiff spring, so that there is a small relative movement even if the force required for sliding is not exceeded.

The definition of Tt,crit depends on the algorithm used for computing the normal contact.

The frictional coefficient μ is defined as

where μs is the static frictional coefficient, μd is the dynamic frictional coefficient, vs is the slip velocity, and αdcf is a decay coefficient.

You can model a situation where two boundaries stick together once they get into contact by adding an Adhesion subnode to Contact. Adhesion can only be modeled when the penalty contact method is used. The adhesion formulation can be viewed as if a thin elastic layer is placed between the source and destination boundaries when adhesion is activated.

When adhesion is active, the tensile stress σn in the normal direction is computed from the effective gap dg as

The default value for the stiffness in the normal direction Kp is the penalty stiffness of the contact condition, but you can also enter it explicitly. In compression, the penalty stiffness is always used. This leads to a bilinear stiffness in the normal direction when Normal stiffness is User Defined.

where nτ is a coefficient with the default value 0.17. This coefficient can either be input explicitly, or be computed from a Poisson’s ratio. A plane strain assumption is used for this conversion, giving

where s is the slip vector and Ft is the deformation gradient on the boundary.

When adhesion is active, it is possible to break the bond between the source and destination boundaries by specifying a Traction separation law. Three different such decohesion laws are available.

The decohesion behavior in both tension or shear is based on two quantities, the maximum allowed stress and the energy release rate. Tearing in the normal direction is called mode I, and shear is called mode II. In the theory below, a subscript i can take the values I and II.

As a measure of the displacement in the adhesive layer, the mixed mode relative displacement um, is defined as weighted combination of the normal direction displacement uI and the tangential displacement uII.

has been used. Thus, um is insensitive to compressive displacement. It is um which is compared with the various mixed mode failure displacements described below.

In the linear separation law, the decreasing part of the stress vs. displacement curve is linear. The graph in Figure 3-20 is generic, and valid both for mode I and mode II.

The failure initiation displacement ui0 for each mode is determined from the maximum stress and the elastic stiffness, as

The decohesion displacement uif for each mode is implicitly determined from the given energy release rate (which is the area under the stress vs. displacement curve):

The final mixed mode failure displacement depends on the selected Failure criterion. The Power Law failure criterion is defined as

whereas the Benzeggagh-Kenane failure criterion is defined as

The exponent η is called the Mode mixity exponent.

The secant stiffness of the damaged adhesive layer (indicated by the red line in Figure 3-20) is

When D > 0, this stiffness is relevant only when the displacement is below the maximum experienced displacement, that it when moving on the red line in Figure 3-20, Figure 3-21, and Figure 3-22.

In an onloading state with D > 0, the stiffness is negative. The negative slope of the curve

The polynomial separation law differs from the linear separation law in that the decreasing part of the stress curve is represented by a cubic polynomial. This polynomial is chosen so that it has zero slope in both ends, which means that the secant stiffness does not have any discontinuities. The area under the cubic function is the same as that of the linear function when the failure displacement uif is the same. That is, for a given energy release rate Gic, the same uif is obtained.

All expressions presented under Linear Separation Law are still valid except the damage function, which now is

The multilinear separation law introduces a region of constant stress, similar to a perfectly plastic material model, as shown in Figure 3-22.

This law requires one more material parameter λ, describing the width of the “plastic” region. The Shape factor λ is the ratio between the plastic (constant stress) part of Gic and the total “inelastic” part of Gic:

Note that the shape factor is assumed to be the same for both tension and shear. The value λ = 0 corresponds to the linear separation law.

The plastic displacement uip can thus be expressed as

The mixed mode plastic failure displacement up is defined as