Skip to content

Constitutive interface formulation (friction)

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) t_{N}-c \end{aligned}\label{flow_rule}\]

with the friction angle \(\varphi\) and the adhesion \(c\) of the contact surface. 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) t_{N} \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 AFR is currently available in numgeo.

As is shown in this Figure, 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.

Illustration of the Coulomb friction model using the tangential penalty factor G ($\epsilon_T$ in the text)

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 view 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\), Eq. (\(\ref{eq:tangential-contact-stiffness}\)) reduces to:

\[ \epsilon_T = s_T G^{mat} ~~~\text{with:}~~~ G^{mat}=\dfrac{1}{(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\)

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