Cone penetration test: Mohr–Coulomb model with updated-Lagrangian kinematics (numgeo vs. Abaqus)
This benchmark compares axisymmetric cone penetration simulations performed with numgeo and Abaqus. The geometry, mesh topology, soil parameters, initial stress state, surface surcharge, and penetration history are matched as closely as the program-specific implementations allow. The purpose is to compare the total vertical cone reaction and the stress field during frictional large-deformation contact.
The dry soil is modelled with an elastic–perfectly plastic Mohr–Coulomb model. The cone is linear elastic and much stiffer than the soil. Penetration is prescribed in an implicit dynamic step, and geometric nonlinearity is enabled so that the mesh and contact configuration are updated throughout penetration.
Figure 1: Horizontal stress during cone penetration in Abaqus and **numgeo**.
Tip
The numgeo and Abaqus input files, the numgeo contact-property subroutine, reference data, and post-processing scripts can be downloaded here.
Material model
The numgeo soil definition is:
*Material, name=Kali3, phases=1
*Mechanical = Mohr-Coulomb-2
50d3, 0.3, 5, 29.8, 9.74, 0.
*Density
1.42
The parameters are:
| Parameter | Value |
|---|---|
| Young's modulus \(E\) | 50 MPa |
| Poisson's ratio \(\nu\) | 0.3 |
| Cohesion \(c\) | 5 kPa |
| Friction angle \(\varphi\) | \(29.8^\circ\) |
| Dilation angle \(\psi\) | \(9.74^\circ\) |
| Tensile cut-off \(p_t\) | 0 kPa |
| Density \(\rho\) | 1.42 t/m\(^3\) |
Mohr-Coulomb-2 expects the friction and dilation angles in degrees. It uses the exact Mohr–Coulomb yield surface and returns an elastic material Jacobian; see the material reference.
The corresponding Abaqus model uses its built-in elasticity and Mohr–Coulomb plasticity definitions; no material user subroutine is required:
*Elastic
50d3, 0.3
*Mohr Coulomb
29.8, 9.74
*Mohr Coulomb Hardening
5, 0
5, 0.1
*Density
1.42
Initial state, boundary conditions, and penetration
The soil is initialised with a geostatic stress field and subjected to a 107 kPa surcharge at its top boundary. The base is fixed vertically, the outer radial boundary is fixed horizontally, and the cone is constrained horizontally.
*Initial conditions, type=stress, geostatic
soil_all, 0, -107, -1.53, -128.7, 0.5, 0.5
*Dsload, instant
surf_soil_top, P, -107
*Boundary
cpt_all, u1, 0.0
soil_bottom, u2, 0.0
soil_right, u1, 0.0
The penetration amplitude in the supplied numgeo input is:
*Amplitude, name=LoadingRamp20mms, type=Ramp
0.0, 0.0, 20.0, 1.0
*Boundary, amplitude=LoadingRamp20mms
cpt_top_left, u2, -1.0
Despite the historical amplitude name, these values prescribe a constant penetration rate of \(1/20\) m/s, or 50 mm/s. The penetration step lasts 1.5 s, giving a final prescribed penetration of 75 mm. Consequently, the time coordinate in Figures 2 and 3 can be converted to penetration depth using \(z=0.05t\) m.
Contact and friction
The numgeo model uses element-based mortar (EBM) contact with a normal penalty formulation and Coulomb friction. Separation is permitted. The supplied archive uses the legacy compact interaction syntax, which remains valid for this benchmark:
*Interaction, name=penalty, Mechanical=penalty, stiffness_factor=40, user
*Friction, model=Coulomb
*Contact Pair, interaction=penalty, discretisation=element mortar
surf_soil_CPT, surf_CPT_soil
The normal stiffness factor is set to 40, twice the default value of 20. It scales the representative stiffness derived from the adjacent continuum element. Because the non constant option is not active, this normal penalty value is not recalculated during penetration. The current syntax and parameter definitions are documented under contact interaction and normal contact behaviour.
The user option activates the supplied user_contact_properties routine. The current archive sets the tangential penalty stiffness to \(\varepsilon_T=30\times10^4\) in the consistent model-unit system and assigns \(\mu=0.381\)—equivalent to an interface friction angle of approximately \(20.8^\circ\)—only to the physical cone and shaft contact region. Friction is set to zero on the artificial zipper extension beneath the physical cone:
props(:) = 0.0d0
props(1) = 30d4
if (physical_contact_region) then
props(2) = 0.381
endif
The pseudocode above summarises the spatial condition; the complete coordinate-dependent implementation is included in the downloadable archive. User-defined contact properties are needed here because the friction coefficient changes along the contacting surface. The underlying Coulomb law is described in the contact theory.
In Abaqus, the corresponding spatial distinction is implemented with a field-dependent friction coefficient:
*Surface Interaction, name=IntProp-1
1.
*Friction, dependencies=1
0.0, , , 0.0, 0.0, 0.0
0.381, , , 0.5, 0.5, 0.5
Comparison of results
Figures 2 and 3 plot the total vertical reaction force of the cone against analysis time. The ordinate is therefore not the area-normalised cone resistance \(q_c\). If \(q_c\) is required, the reaction force must be divided by the projected cone area; the time coordinate can additionally be converted to penetration depth with the relation given above.
Figure 2 compares three tangential penalty stiffnesses, \(\varepsilon_T=56\times10^4\), \(40\times10^4\), and \(30\times10^4\), with the Abaqus result. The overall force level is similar in both programs. Increasing \(\varepsilon_T\) enforces the elastic sticking range more strongly but also increases the high-frequency oscillations in the numgeo reaction force. Of the three investigated values, \(30\times10^4\) gives the smoothest response and the closest visual agreement with Abaqus.
Figure 2: Total cone reaction for three tangential penalty stiffnesses in **numgeo**, compared with Abaqus.
Figure 3 uses the default normal stiffness factor of 20 together with \(\varepsilon_T=30\times10^4\). Compared with the factor of 40 used in Figure 2, the lower normal penalty permits more contact penetration and produces markedly larger reaction-force oscillations. This result illustrates the sensitivity of zipper-based penetration analyses to the normal penalty parameter. A very large penalty is not automatically preferable, however, because it can cause ill-conditioning and convergence problems; the factor should therefore be checked in a sensitivity study.
Figure 3: Total cone reaction using the default normal stiffness factor in **numgeo**, compared with Abaqus.