Skip to content

Contact tangential behaviour

The frictional stress between contacting surfaces with relative tangential movement is determined by the constitutive contact model, also termed the constitutive interface model. For the modelling of pile driving processes or axial loading of piles in general, the adequate incorporation of frictional forces between pile and soil is a key point for the successful prediction of the overall process. In most geomechanical analyses, a simple Coulomb interface model is used to consider friction. However, the constitutive behaviour of soil-structure interfaces is more complex than can be captured by the simple Coulomb friction model.

The mechanical behaviour of sand-structure interfaces is strongly influenced by the surface roughness \(\kappa\). For relatively smooth surfaces the contact friction angle is often fully mobilised after a short shear path and remains almost constant after mobilisation 1. For rough surfaces, however, the behaviour is similar to the response observed in simple shear or triaxial tests performed on dense sand, where the peak stress is followed by gradual softening with ongoing shearing. Compared to smooth surfaces, the shear band is usually much broader for rough surfaces. These characteristics of the constitutive interface behaviour of sand are schematically shown in this Figurea.

Interface shear tests are performed to determine the mechanical properties of soil-structure interfaces and can, amongst other criteria, categorised according to their boundary conditions. Figureb displays schematics of simple interface shear tests with constant normal stiffness (\(K\) constant, CNS test) and constant normal load (\(t_N\) constant, CNL test). For most BVPs, constant normal stiffness conditions are believed to be more representative, since the shearing of the interface leads to volumetric strain and hence change in the normal contact stress. The influence of the constant normal stiffness and constant normal load boundary conditions on the shear stress vs. tangential displacement curves is also schematically shown in this Figureb. CNS tests do not exhibit a peak in the \(t_T\) vs. \(u_T\) plot and reach higher maximum shear stress with increasing normal stiffness. In case of CNL tests, the curves show a peak while higher normal load results in a higher maximum shear stress. Larger relative displacements between the soil sample and the structure are necessary in case of higher normal load to reach the peak value of \(t_T\) compared to lower normal load .

a) Schematic of the shear stress vs. tangential displacement response of soil-structure interfaces in dependence of the surface roughness (based on ). b) Different boundary conditions used in simple interface shear tests and their influence on the shear stress vs. tangential displacement response.

Coulomb friction

A simple relationship for the sliding condition \(\phi\) is obtained using the well-known Coulomb's law of friction. It is given by

\[\begin{aligned} \phi=||t_{Tj}|| - \tan(\varphi) \langle t_N \rangle - a_T \end{aligned}\label{flow_rule}\]

with \(\langle \bullet \rangle = \mathrm{max}(0,\bullet)\), the friction angle \(\varphi\), and the adhesion \(a_T\) of the contact surface. The tangential adhesion \(a_T\) is only active when the contact interface is closed. This is ensured by linking its activation to the normal contact condition:

\[ a_T = \begin{cases} a_T^0 & \text{for } t_N \geq -a_N \\ 0 , & \text{for } t_N < -a_N \end{cases} \]

where \(a_T^0\) is the prescribed adhesion.

This means that adhesion in the tangential direction can only be mobilised if the interface is either in compression or within the admissible tensile range defined by the normal adhesion \(a_N\). Once the normal traction exceeds this limit (i.e. the interface opens), tangential adhesion vanishes.

The admissible stress domain is illustrated in Figure 1 (light green region). Due to the use of the Macaulay brackets \(\langle t_N \rangle\), frictional resistance is only mobilised under compression. In tension \((t_N<0)\), no frictional contribution is present, and only tangential stresses arising from adhesion \(a_T\) can be transmitted.

It is common in more advanced models to reduce \(a_T\) as a function of tangential slip using a damage or softening law. In the current implementation, however, \(a_T\) is assumed to be constant.

Under compressive normal stress, the maximum transferable shear stress increases linearly with \(t_N\) governed by the friction coefficient \(\mu = \mathrm{tan}(\varphi)\).

Figure 1: Illustration of the admissible state of the Coulomb friction model

The direction of plastic sliding is obtained from the plastic potential \(q\). For the application in geotechnics the form

\[\begin{aligned} q =||t_{Tj}||-\tan(\psi) \langle t_N \rangle \end{aligned} \label{plasticPot}\]

with the dilatancy angle \(\psi\) is common. The plastic potential is only an auxiliary construct in order to be able to use a non-associated flow rule (nAFR) in addition to the so-called associated flow rule (AFR) (\(m_i=n_i\),\(\psi=\varphi\)). In the nAFR, the flow rule \(m_i\) is determined from Eq. (\(\ref{flow_rule}\)) vertically (\(\psi\)=0) standing on the flow surface or inclined (\(\psi \neq 0\)) using the dilatancy angle \(\psi\). Only the nAFR is currently available in numgeo.

As is shown in this Figure 2, the elastic (sticking) range of the Coulomb friction model either requires the definition of a relative tangential displacement \(u_{T, crit}\) required to fully mobilise the maximum shear stresses \(t_{T, crit}\), indirectly defining the tangential penalty factor of the frictional contact \(G\), or the direct definition of \(\epsilon_T\) as model parameter. The latter approach is implemented in numgeo.


Figure 2: Illustration of the Coulomb friction model using the tangential penalty factor \(\epsilon_T\)

Typically, the tangential penalty factor \(\epsilon_T\) (or tangential stiffness) is given by

\[ \epsilon_T = s_T G^e \label{eq:tangential-contact-stiffness} \]

where \(G^e\) is a representative underlying element shear stiffness and \(s_T\) is a scaling factor.

As for the normal stiffness \(K^e\) in Contact normal stiffness, only few open literature present formulas to estimate \(G^e\). In these, \(G^e\) is often defined as a fraction of the contact normal stiffness \(K^e\) or as a scaled quantity based on the shear stiffness of the adjacent continuum element.

The latter option is attractive for the analysis of soils, for which the shear modulus is a function of stress, density and additional state variables. Using the stiffness-based approach, \(G^e\) essentially can consider all the effects also considered by the adopted continuum model (incrementally non-linearity, effects resulting from load reversals). To do so, at each contact integration point (see here), the constitutive Jacobian \(\textbf{J}\) of the adjacent element is determined by calling the material routine. Then, \(G^e\) is defined by 2

\[ G^e = \dfrac{2}{d_s \cdot (ndim\cdot2-3)} \sum^{2\cdot ndim}_{i=4} J_{ii} \label{eq:rep-shear-modulus} \]

\(d_s\) is the virtual thickness of the shear zone and has the unit of a length. \(\textbf{J}\) has the unit of stress and \(G^e\) is force per volume. Multiplying \(\epsilon_T\) by \(u_T\) then gives \(t_T\), which is the shear stress. From a soil-mechanical perspective, \(d_s\) could be approximated as a multiple of the mean grain diameter of the soil 2. However, such small values for \(d_s\) might lead to a significant increase in \(G^e\) which could lead to ill-conditioning of the global stiffness matrix and to convergence issues. To avoid such numerical difficulties numgeo uses \(d_s=0.5\) m per default. This value was determined based on comparative studies with the commercial finite element code Abaqus. With \(d_s=0.5\) m, Eq. (\(\ref{eq:tangential-contact-stiffness}\)) reduces to:

\[ \epsilon_T = s_T G^{mat} ~~~\text{with:}~~~ G^{mat}=\dfrac{4}{(ndim\cdot2-3)}\sum^{2\cdot ndim}_{i=4} J_{ii} \]

\(G^{mat}\) is being determined automatically by numgeo. The scaling factor \(s_T\) has a default value of 1 and can be modified by the user.

Prescribing \(\epsilon_T\)

Instead of using above approach to calculate \(\epsilon_T\) based on the shear modulus of the adjacent element, the user can also supply \(\epsilon_T\) using a user subroutine. In this case, \(\epsilon_T\) can as well be spatially changing as well as change with time depending on the users' coding.

The representative tangential stiffness \(G^e\) is evaluated for the two contacting elements of the interface section. The minimum of both is used to calculate \(\epsilon_T\)

Mohr-Coulomb friction

The Mohr-Coulomb friction model extends the Coulomb friction model by incorporating dilatancy through a non-associated flow rule. While the Coulomb model treats the normal and tangential contact as uncoupled during plastic sliding (the normal stress is unaffected by tangential slip), the Mohr-Coulomb model captures the physically observed phenomenon that rough interfaces dilate upon shearing, which in turn modifies the normal contact stress through the penalty constraint.

The yield function is identical to the Coulomb model:

\[ \phi = \| t_{T} \| - \mu \, \langle t_N \rangle - a_T \]

with \(\mu = \tan(\varphi)\). The plastic potential, however, now plays an active role:

\[ q = \| t_{T} \| - \tan(\psi) \, \langle t_N \rangle \]

In the Coulomb model the dilatancy angle \(\psi\) does not influence the stress return because normal and tangential directions are treated independently. In the Mohr-Coulomb model, the flow direction derived from \(q\) couples the tangential slip to the normal direction, so that plastic sliding produces a normal displacement increment proportional to \(\tan(\psi)\).

Elastic trial state

Given the stress state \(\boldsymbol{t}^{n}\) from the previous converged increment and the displacement increment \(\Delta \boldsymbol{u}\), the trial state is computed as:

\[ t_N^{\mathrm{trial}} = -\epsilon_N \left( g_N - c \right) - p^w \]
\[ t_{T,j}^{\mathrm{trial}} = t_{T,j}^{n} + \epsilon_T \, \Delta g_{T,j} \qquad j = 1, \ldots, n_{\mathrm{dim}}-1 \]

where \(\epsilon_N\) is the normal penalty factor, \(\epsilon_T\) is the tangential penalty factor, \(g_N\) is the normal gap, \(c\) is the clearance, \(p^w\) is the pore water pressure and \(\Delta g_{T,j}\) is the tangential gap increment. The normal contact uses a total formulation (via the penalty method), while the tangential contact is incremental.

Yield check

The friction capacity at the trial normal stress is:

\[ \tau_{\max} = \mu \, t_N^{\mathrm{trial}} + a_T \]

The yield function is evaluated at the trial state. In 2D:

\[ \phi^{\mathrm{trial}} = | t_T^{\mathrm{trial}} | - \tau_{\max} \]

In 3D, the resultant trial shear stress is used:

\[ \phi^{\mathrm{trial}} = \| \boldsymbol{t}_T^{\mathrm{trial}} \| - \tau_{\max} \qquad \text{with} \quad \| \boldsymbol{t}_T^{\mathrm{trial}} \| = \sqrt{ (t_{T,1}^{\mathrm{trial}})^2 + (t_{T,2}^{\mathrm{trial}})^2 } \]

If \(\phi^{\mathrm{trial}} \leq 0\), the trial state is admissible (elastic stick) and accepted without modification. If \(\tau_{\max} \leq 0\) (tensile normal stress exceeding the adhesion), no friction can be mobilised and the shear stress is set to zero.

Return mapping (plastic slip)

When \(\phi^{\mathrm{trial}} > 0\) and \(\tau_{\max} > 0\), the stress must be returned to the yield surface. The return mapping is performed in a single step. The key quantity is the plastic denominator:

\[ h = \epsilon_T + \epsilon_N \, \mu \, \tan(\psi) \]

Note that for \(\psi = 0\) this reduces to \(h = \epsilon_T\). The plastic multiplier is:

\[ \Delta \lambda = \frac{\phi^{\mathrm{trial}}}{h} \]

The returned stresses are:

\[ t_N = t_N^{\mathrm{trial}} + \Delta \lambda \, \epsilon_N \, \tan(\psi) \]
\[ t_{T,j} = d_j \left( \| \boldsymbol{t}_T^{\mathrm{trial}} \| - \Delta \lambda \, \epsilon_T \right) \]

where \(d_j = t_{T,j}^{\mathrm{trial}} / \| \boldsymbol{t}_T^{\mathrm{trial}} \|\) is the unit direction of the trial shear traction (in 2D: \(d = \mathrm{sign}(t_T^{\mathrm{trial}})\)). The returned shear stress magnitude equals \(\mu \, t_N + a_T\) exactly, which can be verified by substitution.

The dilatancy effect is visible in the expression for \(t_N\): during plastic slip, the normal stress increases by \(\Delta \lambda \, \epsilon_N \, \tan(\psi)\). Physically, the interface dilates upon shearing, and the penalty constraint converts this dilation into additional compressive normal stress, which in turn increases the friction capacity. For \(\psi = 0\), the normal stress remains unmodified and the model reduces to the Coulomb friction model.

Consistent elastoplastic tangent

The consistent (algorithmic) tangent \(\mathbf{D}^{ep}\) is derived from the standard formula for single-surface perfect plasticity:

\[ \mathbf{D}^{ep} = \mathbf{D}^{e} - \frac{(\mathbf{D}^{e} \, \boldsymbol{m})(\boldsymbol{n}^T \mathbf{D}^{e})}{\boldsymbol{n}^T \mathbf{D}^{e} \, \boldsymbol{m}} \]

where \(\boldsymbol{n} = \partial \phi / \partial \boldsymbol{t}\) is the yield surface gradient and \(\boldsymbol{m} = \partial q / \partial \boldsymbol{t}\) is the flow direction, both evaluated in the local contact stress space \(\boldsymbol{t} = [t_N, t_{T,1}, \ldots]^T\).

2D

The elastic stiffness, yield gradient and flow direction are:

\[ \mathbf{D}^{e} = \begin{bmatrix} \epsilon_N & 0 \\ 0 & \epsilon_T \end{bmatrix}, \quad \boldsymbol{n} = \begin{bmatrix} -\mu \\ s \end{bmatrix}, \quad \boldsymbol{m} = \begin{bmatrix} -\tan(\psi) \\ s \end{bmatrix} \]

with \(s = \mathrm{sign}(t_T^{\mathrm{trial}})\). The denominator is \(\boldsymbol{n}^T \mathbf{D}^{e} \, \boldsymbol{m} = \epsilon_N \, \mu \, \tan(\psi) + \epsilon_T = h\). The resulting tangent is:

\[ \mathbf{D}^{ep} = \frac{1}{h} \begin{bmatrix} \epsilon_N \, \epsilon_T & s \, \epsilon_N \, \epsilon_T \, \tan(\psi) \\[4pt] s \, \epsilon_N \, \epsilon_T \, \mu & \epsilon_N \, \epsilon_T \, \mu \, \tan(\psi) \end{bmatrix} \]

For \(\psi = 0\): \(h = \epsilon_T\) and

\[ \mathbf{D}^{ep}\big|_{\psi=0} = \begin{bmatrix} \epsilon_N & 0 \\ s \, \mu \, \epsilon_N & 0 \end{bmatrix} \]

which is identical to the tangent of the Coulomb friction model during slip.

3D

In 3D, the return is radial in the \((t_{T,1}, t_{T,2})\) plane. The consistent tangent has the structure:

\[ D^{ep}_{11} = \frac{\epsilon_N \, \epsilon_T}{h} \]
\[ D^{ep}_{1,j+1} = d_j \, \frac{\epsilon_N \, \epsilon_T \, \tan(\psi)}{h} \]
\[ D^{ep}_{j+1,1} = d_j \, \frac{\epsilon_N \, \epsilon_T \, \mu}{h} \]
\[ D^{ep}_{j+1,k+1} = r \, \epsilon_T \left( \delta_{jk} - d_j \, d_k \right) + d_j \, d_k \, \frac{\epsilon_N \, \epsilon_T \, \mu \, \tan(\psi)}{h} \]

where \(r = \| \boldsymbol{t}_T \| / \| \boldsymbol{t}_T^{\mathrm{trial}} \|\) is the ratio of returned to trial shear stress magnitude, and \(\delta_{jk}\) is the Kronecker delta. The first term in \(D^{ep}_{j+1,k+1}\) represents the secant stiffness perpendicular to the sliding direction (arising from the change in the return direction \(d_j\) upon perturbation of the trial stress), while the second term captures the stiffness in the sliding direction due to the dilatancy-friction coupling.

Continuum models as contact models

Hypoplastic friction model

Please refer to 2.

Sanisand friction model

Please refer to 2.

  1. Jakob Vogelsang. Untersuchungen zu den Mechanismen der Pfahlrammung. PhD thesis, Veröffentlichung des Instituts für Bodenmechanik und Felsmechanik am Karlsruher Institut für Technologie (KIT), Heft 182, 2016. ↩

  2. Patrick Staubach, Jan Machaček, and Torsten Wichtmann. Novel approach to apply existing constitutive soil models to the modelling of interfaces. International Journal for Numerical and Analytical Methods in Geomechanics, 46(7):1241–1271, may 2022. URL: https://onlinelibrary.wiley.com/doi/10.1002/nag.3344, doi:10.1002/nag.3344. ↩↩↩↩