Contact normal behaviour
The penalty regularisation approach used to enforce the normal contact constraints requires the penalty factor \(\epsilon\) between the two contact surfaces. The amount of penetration between the two surfaces depends on the magnitude of \(\epsilon\). Higher values decrease the amount of penetration but can lead to ill-conditioning of the global stiffness matrix and to convergence issues. Typically, the stiffness is given by
\[
\epsilon = s K^e
\]
where \(K^e\) is a representative underlying element stiffness and \(s\) is a scaling factor. Only few open literature present formulas to estimate \(K^e\). In numgeo we approximate \(K^e\) based on the constitutive Jacobian \(\textbf{J}\) of the adjacent element:
\[
K^e = \boldsymbol{J}_{ii}/3.
\]
The default value of the stiffness factor is \(s=20\).