Why Do Parts Collide in Nonlinear Static Simulations?

• The penalty method is based on the assumption that two bodies collide (overlap) and an imaginary spring is compressed during the collision.
• The amount of collision is called penetration depth.
• The contact force is then calculated with respect to the penetration depth and stiffness coefficient of the spring.
• The penalty coefficient is the stiffness of the contact (imaginary spring).

  $$F = k \cdot \Delta L $$

  Where, \(k\) is the stiffness of the contact (Penalty coefficient), \(F\) is the contact force, and \(DeltaL\) is the penetration depth.

In the penalty method, it is important to assign a suitable penalty coefficient. A suitable penalty coefficient should keep the penetration depth at a minimum:
If it’s too small, then bodies will overlap.
If it’s too big, it will generate numerical instabilities.
A good practice is to do the following:

• Estimate k (penalty coefficient) and perform an initial simulation
• Check the solver log. Are the Contact penetration residuals small enough? A good threshold might be 1e-6.
• If the contact penetration is small enough, check the solution fields. View displacement results and check the penetration. Is the penetration negligible?
• If penetration residual or penetration depth is large, then increase the k. If the simulation fails (The solution matrix is singular error), then better decrease the k. To tune the k, multiply it with by the factor of 10, 100, or 0.1, 0.01.







Determination of the Penalty Coefficient

The penetration coefficient depends on many factors, such as mesh type, element size on contact surfaces, material properties, load, etc. There is not any particular way to define the penalty coefficient. However, you can use one of the following recommendations. The complexity increases from top to bottom:

1. As a rule of thumb, you can multiply the Young’s modulus by 50 [1]
   $$k = 50 \cdot E $$
   where, \(E\) is the Young’s modulus.
2. G.R. Liu recommends to multiply the Young’s modulus by a factor of 1e5 or 1e8 [2]
   $$k = 1 \cdot 10^{5 \sim 8} \cdot E $$
3. In Code_aster documentation, it is recommended to multiply the largest Young’s modulus of the bodies in contact by a characteristic length [3]
   $$k = d \cdot E_{max} $$
   where, \(d\) is the characteristic length.
4. Zienkiewicz recommends to use the following equation [4]
   $$k = E_{max} \cdot (\frac{1}{d})^p $$
   The characteristic length can be defined as the ratio of the element size to the dimension of the problem domain.
   $$d = \frac{\Delta x}{n} $$
   where, \(d\) is the characteristic length, \(p\) is the order of elements, \(Delta x\) is the element size (nodal spacing) of the smallest element in contact face, and \(n\) is the dimension of the problem domain (SimScale requires 3D domain, therefore n = 3)

In contrast to the Penalty coefficient, penetration is completely forbidden in the Lagrange method.
Thus, the contact force is independent from cell size. This method requires an exact solution (Lagrange equations) instead of iterations, which makes it accurate but unstable in large deformations.

The Augmented Lagrange method can be considered as a combination of Penalty and Lagrange methods. It causes minimal penetration and solution is almost exact. This method is accurate, but not suitable for problems with more than 100,000 nodes. It is also less stable, especially with second order elements.