Contact Analysis Theory
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 penalty method, no extra degrees of freedom are needed for the contact pressure or the friction traction vector. These results are just computed from the displacements and the penalty stiffness. This means that no special solver strategies are necessary.
Identity and Contact Pairs in the COMSOL Multiphysics Reference Manual
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.
Before the boundaries come in contact, the source point found is not necessarily the point on the source boundary closest to the destination point. However, as the boundaries approach one another, the source point converges to the closest point as the gap distance approaches zero.
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
The correction of the source side offset is necessary since the normal on the source side is not necessarily pointing toward the current point on the destination. When in contact, the normals will be exactly opposite.
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
where the index ‘m’ indicates that the coordinates are taken as the mapped coordinates from the source side,
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.
Augmented Lagrangian Method
Using the special gap distance variable (solid.gap), the penalized contact pressure Tnp is defined on the destination boundary as
(3-90)
pn is the user defined normal penalty factor.
The penalized friction traction Ttp is defined on the destination boundary as:
(3-91)
where Tt,trial is defined as
(3-92)
In Equation 3-92, pt is the user defined friction traction penalty factor. The critical sliding resistance, Tt,crit, is defined as
(3-93)
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.
Penalty Method
The penalty algorithm is essentially based on the insertion of a stiff spring between the contacting boundaries. The penalty factor can be interpreted as the stiffness of that spring. A high value of the stiffness fulfills the contact conditions more accurately, but if it is too high, it makes the problem ill-conditioned and unstable. The contact pressure in the normal direction is computed as
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
where
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 same notation as for the description of the augmented Lagrangian method above is used.
The springs between the contacting boundaries are added using weak expressions on the destination boundary.
Friction Models
The friction model is either no friction or Coulomb friction.
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.
Directions of the Contact Forces
Since the contact model always makes use of the geometric nonlinearity assumption, the contact pressure and the friction force represent, respectively, the normal and tangential components of the nominal traction. This is the force density with components in the current configuration (that is, on the spatial frame) but related to the undeformed area of the corresponding surface element.
Adhesion
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.
The adhesion starts acting when the adhesion criterion is met in the previous time or parameter step. An auxiliary degree of freedom located at Gauss points is used as an indicator of whether the adhesion criterion has been met or not.
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.
The shearing deformation is resisted using a shear stiffness Kτ, related to the normal stiffness by
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
The shear traction is proportional to the relative deformation between the source and destination boundaries measured from the position when the adhesion criterion was first fulfilled, so that
where s is the slip vector and Ft is the deformation gradient on the boundary.
Decohesion
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.
Up to the maximum stress, the adhesive layer is linearly elastic. At higher deformation, the stress decreases, and a damage is assumed to have occurred. The damage function is a scalar function of the maximum relative displacement which has occurred between the two boundaries. The stiffness in the damaged state is permanently reduced, even on unloading.
In general, the stress state is a combination of tension and shear, so the data for the two modes must be combined into a multiaxial decohesion rule. For two of the decohesion laws you can specify details of the multiaxial behavior.
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.
In this expression, the notation
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.
Linear Separation Law
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.
Figure 3-20: Stress vs. displacement curve for the linear separation law
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 mixed mode failure initiation displacement is then defined as
The implication is that compressive displacements are ignored, and only tensile and tangential displacements contribute.
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.
From these relations, the total decohesion failure displacement in the mixed mode can be computed. For a power law criterion, the expression is
For the Benzeggagh-Kenane criterion, the corresponding expression is
The damage function can be described as
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
Polynomial Separation Law
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.
Figure 3-21: Stress vs. displacement curve for the polynomial separation law
All expressions presented under Linear Separation Law are still valid except the damage function, which now is
Multilinear Separation Law
The multilinear separation law introduces a region of constant stress, similar to a perfectly plastic material model, as shown in Figure 3-22.
Figure 3-22: Stress vs. displacement curve for the multilinear separation law
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 decohesion displacement is
In this model, a linear mixed mode failure law is assumed:
The mixed mode plastic failure displacement up is defined as
The mixed mode failure displacement is
The damage function can be written as