Geotechnical particle finite element method for modeling of soilstructure interaction under large deformation conditions

Josep Maria Caronell, Lluís Monforte, Matteo O. Ciantia, Maros Arroyo,Antonio Gens

a MECAMAT Group, Department of Engineering, Faculty of Science and Technology, Universitat de Vic-Universitat Central de Catalunya (UVic-UCC), Vic, Spain

b Centre Internacional de Mètodes Numèrics a l′Enginyeria (CIMNE), Barcelona, Spain

c School of Science and Engineering, University of Dundee, Dundee, UK

d Universitat Politècnica de Cataluna - BarcelonaTech, Barcelona, Spain

Keywords:Particle finite element method (PFEM)Structured soils Nonlocal elastoplasticity Contact domain method Soil penetration problems

A B S T R A C T

1. Introduction

Situations in which soil masses are subject to large displacement, frequently involving contact with structures, are ubiquitous in geotechnical engineering: sampling, in situ testing, pile installation, landslides, to name a few. Numerical simulation is clearly useful to advance understanding in these areas. However,the numerical simulation of such problems is a complex task,since the system is full of nonlinearities, contact-related, materialrelated, and also geometrical. The finite element method (FEM)(Zienkiewicz et al.,1999,2005),which has proven successful for the numerical simulation of several multiphysics problems, becomes unreliable once geometric nonlinearities intervene. In Lagrangian formulations,the mesh generally becomes highly distorted,leading to inaccurate results,loss of convergence and calculation stoppage at relatively small displacements (De Borst and Vermeer,1984).

In the past decades, several numerical methods to tackle large strain problems in geotechnical engineering have been proposed to overcome the problems of classic FEM. The application of FEMrelated numerical techniques to model large strain problems in geomechanics can be traced back to the arbitrary Lagrangian Eulerian(ALE)proposal of van den Berg et al.(1996).Subsequently,several ALE frameworks have been developed, such as the remeshing interpolation technique by small strains(RITSS)(Hu and Randolph, 1998) or the efficient ALE (Nazem et al., 2006; Sheng et al., 2009). Other proposals include those based on the material point method (MPM) (Sulsky et al.,1994) or even the point in cell method(Harlow,1964),recently recovered and developed to model problems governed by large movements and deformations with imaginative solutions (de Vaucorbeil et al., 2020). From a noncontinuum perspective, formulations based on the discrete element method (Cundall and Strack, 1979) are also applied to tackling similar problems (e.g. Ciantia et al., 2016; Zhang et al.,2021), although the computational price of DEM is usually steeper, particularly if fluid coupling is included.

During the last decade, the particle FEM (PFEM) has been developed into a viable alternative to deal with large strain problems in geomechanics. PFEM was originally developed for fluid mechanics,using a Lagrangean description of the domain instead of the classical Eulerian approach (Idelsohn et al., 2004; Oñate et al.,2004, 2006). Treating the finite element nodes as particles, a finite element mesh is constructed each time-step to evaluate the solution, thus the mesh is always of good quality. This continuous reconnection of the domain enables rapid boundary detection,thus the method can capture the separation of domains, and in fluid dynamics,detached particles may behave as water drops that may rejoin the domain again. It is worth mentioning that in PFEM, the solution is evaluated with classical finite elements,which facilitates the incorporation of a great deal of FEM know-how.The numerical treatment of the governing equations, boundary conditions and even contact interaction does not differ from what may be employed in FEM,thus any development achieved in classical FEM can be quickly imported into PFEM. Moreover, from an implementation perspective, PFEM algorithms are highly non-intrusive and compatible with FEM solvers, to an extent that some implementation of PFEM even relies on commercial finite element packages (Sabetamal et al., 2021; Yuan et al., 2021). Another important precision is related to the use of the term “particles”,which in PFEM refers to mesh nodes representing “continuum particles” (Tadmor et al., 2012), i.e. the representative volume elements (REVs) that underlie points in a continuum. PFEM is a discretization method to solve equations describing continuum media,fundamentally different in this respect from methods based on direct formulation of element properties and interactions such as DEM. In DEM, elements may be actually designed to represent individual physical particles (e.g. Wu et al., 2021), something excluded by the continuum perspective used in PFEM.

Since its original development, the application of PFEM has spread to problems in solid mechanics(Oliver et al.,2007),thermomechanics (Rodriguez, 2014; Rodriguez et al., 2016), granular materials (Zhang et al., 2013, 2014, 2015; Larsson et al., 2020), fluidstructure interaction (Franci et al., 2016), contact of deformable bodies (Hartmann et al., 2009; Carbonell et al., 2010, 2013), and manufacturing processes(Oñate et al.,2014;Hays,2019),amongst others.In geomechanics,Zhang et al.(2013,2014)proposed a PFEM implementation, using high order elements in conjunction with variational principles, to discretize the governing equations for single-phase continuum media and applied it to studying soil flow problems,including retrogressive landslides in sensitive materials.This formulation was later extended to fluid-saturated bi-phasic porous media(Wang et al.,2021).A different PFEM implementation has been recently proposed using smoothed finite element for the one-phase formulation and applied to the simulation of footings and retrogressive landslides (Zhang et al., 2018, 2021; Yuan et al.,2019).

An implementation called G-PFEM (geotechnical-PFEM) was initially presented by Monforte et al. (2017a). The G-PFEM was based on low-order linear, stabilized finite elements (Monforte et al., 2017b), developed to solve quasi-static hydromechanical problems in fluid saturated porous media (Monforte et al., 2018a)and later extended for dynamic problems using a full Biot formulation(Monforte et al.,2019b).G-PFEM has been mostly applied to studying soil-tool interaction problems, such as cone penetration testing (Monforte et al., 2018b) or soil sampling (Monforte et al.,2021a). These applications have recently emphasized the use of realistic constitutive models for brittle(Monforte et al.,2021b)and structured soils (Monforte et al., 2019a; Hauser and Schweiger,2021).

In Section 2,we present an overview of some key aspects of the numerical technology employed in G-PFEM, starting with the particular ingredients of the method which work in conjunction with the hydro-mechanical problem. The general equations in the strong and weak forms are given and discretized for simplicial finite elements. The numerical examples are developed using an elastoplastic model for structured materials at large strains to describe material behavior. This model is implemented using a nonlocal regularization technique to deal with possible strain localization that is also described.Section 2 includes some aspects of the formulations that have not been presented before, such as the axisymmetric formulation for the field equations or a simplified approach for dynamic problems.

Section 3 is devoted to outline the contact formulations employed in G-PFEM. The first one addresses contacts between a rigid object and deformable media and is the one that has been employed in most G-PFEM applications. The second - newly introduced here - addresses the more general case of contact between two deformable objects and is based on thecontact domainmethod (CDM) (Oliver et al., 2009; Hartmann et al., 2009). The formulations presented are particularized for two-dimensional(2D) problems considering plane-strain and axisymmetric approaches.

Finally,Section 4 presents a number of numerical analysis cases to show the possibilities of G-PFEM. They include novel aspects of problem that have been addressed before,like footing embedment in structured soils or pull-out of screw anchors, but also novel applications of G-PFEM,such as embedded pipeline displacements or free-falling impact of spherical penetrometers in clay.

2. An overview of G-PFEM

2.1. Fundamental aspects of PFEM algorithm

PFEM is a numerical method that combines features of the meshfree particle method and the FEM. As implemented in GPFEM, the method has three main ingredients: the first one is an updated Lagrangian description of the domain, in which the kinematic description uses large displacement and finite strain theory.The governing equations are then solved using different variational principles. Most frequently, as in G-PFEM, using mean weighted residuals (MWRs) like in the Petrov-Galerkin or Ritz-Galerkin methods. Alternatively, the governing equations may also be cast as an optimization problem and solved using second order cone programming(SOCP)(Zhang et al.,2016).The second ingredient is the use of moving particles, i.e. material tied REVs, to define the domain.When particles move,they are reconnected via a Delaunay tessellation (Delaunay, 1934); this supplies a new finite element mesh for the solution of the governing equations in which the particles act as nodes. A downside is that Delaunay reconnection limits the finite element geometry to triangles and tetrahedra for 2D and three-dimensional (3D) models, respectively. The last feature is related to the accuracy of the domain boundary description. After the tessellation, the alpha-shape method(Edelsbrunner and Mucke, 1994) is used to recover the previous domain contour or to identify new boundaries for the computing domain. Note that alpha-shape is not the only technique used to recognize new boundaries, as information on the previous boundary normal is stored in particles and employed to recover contour curvature and improve the quality of the reconnected domain.

In solid mechanics applications,the boundary can be preserved using a constrained Delaunay tessellation, since it guarantees mass conservation and avoids a nonphysical increase of the domain volume. In the examples presented in this work, this option is used. After remeshing, G-PFEM projects the internal variables that describe material strength and plastic straining from old Gauss points to the new ones using a minimum distance criterion. Such criterion, which minimizes numerical diffusion in mesh to mesh mapping, is also an improvement with respect to the original PFEM. Inherited from developments in solid mechanics applications, G-PFEM also includes mesh adaptivity. New particles are inserted in regions where plastic flow takes place in order to increase the spatial resolution of plastic internal variables. In particular, elements whose size is larger than a predefined characteristic size may be divided,if the value of a plastic variable multiplied by the area of the element is larger than a threshold. This improves the mesh shape used for the computation and increases solution accuracy.The combination of all these ingredients allows to model problems in which fluids or plastic solids flow,surpassing the mesh distortion limitations inherent to Lagrangian FEM formulations.

A typical solution algorithm of G-PFEM is conceptually illustrated in Fig. 1. The scheme illustrates acollection of particles or cloud of nodes(C) belonging to the analysis domain, as well as the mesh(M)discretizing the domain with finite elements surrounded by a predefined boundary, representing thevolume(V) of the computing domain.The basic steps of the algorithm are as follows:

Fig.1. Sequence of steps to update in time a “cloud” of nodes representing, in G-PFEM, a soil mass that is progressively penetrated by a rigid structure.

(1) Generate aconstrainedDelaunay tessellation (Mnmesh) using previous mesh nodes as thecollection of particles Cnand the current boundary enclosing volumeVn.

(2) Transfer nodal information to new inserted particles and elemental information from previous finite elementmesh M0to the current finite elementmesh Mn.

(3) Perform the FE solution loop fori:

(i) Contact search for domain-structure interactions using the known boundaryVin.

(ii) Solve the hydro-mechanical FE problem in aLagrangianform using the FE meshMin.

(iii) Update particle displacements and positions and water pressuresVin，Min→Vi+1n，Mi+1n.

(iv) Check convergence of the solution. If not converged,perform another iterationi→i+ 1.

(4) Check mesh accuracy and refine domain inserting new particles in thecollection of particles Cn+1and redefining the boundary of the new volumeVn+1.

(5) If the analysis is not ended, go to Step 1.

2.2. Hydro-mechanical mixed formulation

As already pointed out in Section 1, most PFEM proposals for geomechanics rely on a one-phase formulation,where only strictly drained or undrained conditions may be simulated. These simple models have several advantages, amongst them faster computation, but their field of application is limited. A fully coupled formulation modeling bi-phase porous medium extends significantly the range of potential geotechnical problems that may be addressed. Such formulation should also reproduce behavior in both free draining and undrained conditions as limiting cases.

One of the important characteristics of water-saturated soils is that they may have a quasi-incompressible behavior. On the one hand, the application of rapid loads in low-permeability material results in undrained conditions and the fluid saturated porous media behave as an incompressible material. On the other hand,incompressibility may also arise under general drainage conditions for constitutive models that predict null volume change (for instance, critical state models) (Sun et al., 2013). These conditions cause a well-known numerical problem of volumetric locking in low-order finite elements leading to numerical instability.To avoid this numerical pathology, more complex but stable second order finite elements can be used.In G-PFEM,locking is mitigated by the use of a mixed formulation solved with low-order stabilized elements. This approach is well adapted to cases where incompressibility may result from the mechanical behavior of the solid phase itself. It is easily applicable to coupled hydromechanical or singlephase undrained formulations. Mixed formulations introduce extra degrees of freedom per node with respect to the primal formulation,however,in most cases,they offer equal performance at lower computational cost than the use of higher order elements.More details on these stabilization techniques are given elsewhere(Monforte et al., 2017b).

In fluid-saturated porous media, the domain is divided in two parts: the solid skeleton and the porous fraction, which in a saturated scenario will be filled by fluid water. According to the principle of effective stress,the total Cauchy stress tensor,σ,is equal to the sum of the pore pressure,pw, and the effective stress,σ′, i.e.

where 1 stands for the second-order identity tensor.

The volume in the current configuration is given by the determinant of the deformation gradientFalso known as the JacobianJ= det(F). The effective Cauchy stress tensor depends on the strains of the solid skeleton σ′= σ′(F，J).If an approximation of the Jacobian θ is introduced as primary particle variable,in addition to displacements and water pressure, the effective Cauchy stress can be evaluated with the assumed deformation gradient′= σ′(，θ)defined as

where the deviatoric part of the deformation gradientFdis preserved whereas the volumetric partFvis replaced with a nodal approximation (the θ variable). Note that, in this formulation, the Cauchy stress tensor depends on both the displacements and the Jacobian (Wriggers, 2008). Despite that, usual strain-driven stress integration schemes are still suitable for this formulation.Although a 3D version has been also developed(Monforte et al.,2018c),most G-PFEM work to date considers either plane strain or axisymmetric problems.In plane strain conditions,the definition of the assumed deformation gradient is modified to guarantee that the out of plane component of the deformation gradient is equal to unity:

In the axisymmetric case (Fig. 2), the deformation gradient is defined as

Fig. 2. Axisymmetric triangular element.

Consequently, the circumferential effective stress′33must be considered in the definition of the effective stress tensor′for an axisymmetric case:

The balance of mass and linear momentum equations for multiple-phase deformable porous media using a displacement-Jacobian-water pressure (u-θ-pw) formulation is written in the current deformed configuration as

where ρ is the mixture density,astands for the acceleration,andgis the acceleration due to gravity. Because of the inclusion of acceleration,this formulation covers also dynamic cases.Note that this is a formulation only suitable in cases where fluid acceleration is far smaller than that of the solid phase; for different conditions,other formulations including more variables are necessary(Zienkiewicz et al., 1999). A full Biot formulation for G-PFEM,covering all conditions, was presented by Monforte et al. (2019b).Note also that in quasi-static problems (the case for most of the simulations reported later), the first term of Eq. (7) is dropped.

The previous equations require appropriate initial conditions at Ω0, prescribed displacementsand tractionstfor the boundary domain ∂Ωt=and fixed water pressurepw and prescribed water flow

gin the saturated media boundary∂Ωt=. These are written as follows:

for the solid domain boundary and

for the fluid domain boundary. Further information about more complex mixed formulations applied to the mechanical and fluid phase parts can be found in Monforte et al. (2017b).

The finite element discretization of Eq. (8) is presented in Appendix A, employing linear shape functions for all nodal variables. The implementation of a numerical stabilization procedure,i.e.polynomial pressure projection(PPP)(Bochev et al.,2006)in the equations is also illustrated in Appendix A. For quasi-static problems,the system of discretized equations is integrated in time using a completely implicit scheme, whereas a Bossak implicit time integration scheme (Crisfield, 1991) is employed when dynamic effects are considered.

2.3. Constitutive description and integration

2.3.1. Constitutive equations

Different constitutive models for soils have been used in GPFEM. Some of them are classical reference models such as Tresca(Monforte et al., 2017a) or Cam Clay (Monforte et al., 2018a). But more realistic models such as the clay and sand model(CASM)(Yu,1998) with (Hauser and Schweiger, 2021) or without bonding(Monforte et al., 2021b) have been also employed. These elastoplastic models are adapted for large strains using a multiplicative decomposition of the deformation into an elastic and plastic part(Simo and Hughes,1998). This framework ensures frame indifference so that any rigid body motion of the deformable body does not produce spurious stress variations. A hyper-elastic model is typically employed to define the behavior of the material in the reversible deformation region, which ensures that energy is preserved.

In this work, we use a modified Cam Clay model enhanced to consider the effect of soil structure(Nova et al.,2003).This type of model is suitable for the description of the behavior of soft rocks,natural clays or even artificially cemented soils (Gens and Nova,1993; Rouainia and Muir Wood, 2000; Wheeler et al., 2003; Rios et al., 2006; Ciantia, 2018). On the other hand, structure results in an enhanced brittleness and is conductive to strain localization.

The yield surface of the model, which is based on that of modified Cam Clay, is graphically illustrated in Fig. 3. On the one hand, the reference unstructured soil has a preconsolidation pressurepsand null tensile strength. On the other hand,pmandptaccount for the effect of structure; the former corresponds to the increase in the yield stress along isotropic compression paths whereas the latter is the tensile strength. Although other hypotheses are possible, generally, these two hardening variables are considered proportional (Ciantia and di Prisco, 2016), beingcthe proportionality factor.Mathematically,the yield locus is defined as

Fig. 3. Yield surface in the p′-q plane for axisymmetric triaxial compression.

whereq=，in whichJ2is the second invariant of the effective deviatoric Kirchhoff stress tensor, τ′;Mis the slope of the critical state line in thep′-qplane, in whichp′= tr(τ′)/3 is the first invariant of the Kirchoff stress tensor;andp*andp*ccan be written as

In correspondence of the modified Cam Clay, the preconsolidation pressure of the unstructured soil depends on the volumetric plastic strains. Due to plastic straining, structure degrades; thus,ptdecreases due to both volumetric and deviatoric plastic strains.

where ρs，ρtand χtare the constitutive parameters; andlpis the spatial plastic velocity gradient.

The elastic response is characterized by means of a hyper-elastic model incorporating a tensile range(Tamagnini et al.,2002),which is formulated in terms of the Hencky strain and the Kirchhoff stress tensor. Finally, associated plasticity is assumed.

2.3.2. Stress integration

Elastoplastic models are usually integrated in G-PFEM using an explicit stress integration technique (Monforte et al., 2014) based on Sloan et al. (2001). The algorithm includes adaptive substepping and a yield surface drift correction algorithm. An alternative stress integration technique is also employed in this work,the so-called IMPLEX technique(Oliver et al.,2008),which provides extra robustness and computability with respect to usual methods.The interested reader is referred to Monforte et al. (2019a) for further details on the application of the IMPLEX technique to this constitutive model.

2.3.3. Nonlocal regularization

The simulation of brittle materials is conductive to strain localization (Zienkiewicz et al., 1995) from which mesh dependency of the solution may result(Galavi and Schweiger, 2010).In the context of finite elements, strain localization is highly influenced by the mesh: the width of the shear band is typically related to the element size whereas its direction is sometimes controlled by preferential alignment of the elements.Therefore,as the mesh is further refined, the thickness of the shear band decreases,and the energy dissipated in the shear band tends to zero.This pathological mesh-dependence may be mitigated using regularization techniques, which incorporate a length scale to the constitutive model thereby enforcing the width of the localized region.

Among regularization techniques,the nonlocal integral type has the advantage of not changing the field equations, which in turn results in a quite straight-forward implementation (Galavi and Schweiger, 2010; Mánica et al., 2018).

A nonlocal integral type regularization technique is used here to mitigate the pathological mesh-dependence that exhibit numerical simulations where softening is encountered.

In this approach, the constitutive model is evaluated replacing some variables with its nonlocal counterpart, which is a spatial average in a neighborhood.Therefore,the constitutive response of a Gauss point is influenced by all the other integration points within a neighborhood, the size of which is determined with a characteristic length. The expression of a nonlocal variable ~β is

wherew(x，||x-y||)is the weighting function for pointxcontrolling the influence of its neighbors in terms of their relative distance.

In this work,plastic strains are considered as nonlocal variables;from these values,psandptmay be obtained by integrating analytically Eqs.(13)and(14).The weighting function proposed by Galavi and Schweiger(2010)is employed here as it has been found to outperform other weighting functions in removing mesh bias.It is shown (Mánica et al., 2018) how, using this approach, the thickness of the shear band is related to a characteristic length scale enforced through the function weighting the neighbors′relative distance.

3. Treatment of the contact interaction

3.1. General aspects

The object of“contact interactions”in G-PFEM are two separate continuum bodies. As a basic pre-requisite, accurate contact interactions require a precise definition of the domain contour; in PFEM, this involves continuous tracking of free surfaces of the domain. In most geotechnical applications, the fragmentation and re-merging of the domain that is characteristic of PFEM applications in fluids is not usually necessary.Taking advantage of that,GPFEM constrains remeshing to preserve the external boundaries of the domain during the analysis, thus obtaining a more precise resolution. This appears somewhat simpler than alternatives employed with other particle techniques, such as smoothed particle hydrodynamics (SPH) or MPM (Bardenhagen et al., 2000;Gonzalez Acosta et al., 2021).

A contact interaction model can be used to capture the mechanical forces emerging from the contact of two material subdomains.The interaction can also consider not only the mechanical contact,but also the flux of heat or water between one domain and another. Usually, there are two main steps to include contact in a finite element model: (i) the geometrical detection and (ii) the numerical formulation of the contact constraint that will apply to the governing equation.A contact model can be considered an extra challenge for any formulation, in this case for the G-PFEM, as it introduces the need for the numerical treatment of a mathematical discontinuity in the numerical analysis. It also allows for the introduction of nonlinear friction models between contacting surfaces,increasing the modeling capabilities and giving more realistic simulations. Penetration problems that will be presented in this paper are not going to require the flux of water between the contacting domains, because one of the domains (the structure) is considered impervious.However,the novel implementation of the contact domain in G-PFEM,presented later,allows for the modeling of this interaction for all primary variable fields of the problem,including the flux of water if needed. In that case, a scalar term related to the Cauchy water pressure must be considered for the Darcy flow equation.

3.2. The penalty method

Two different strategies have been developed in G-PFEM to deal with contact constraints. The first one addresses the case - frequent in geotechnical applications - of an object that is much stiffer than the ground with which interacts.For such cases,it is advantageous to treat that object as a rigid wall applying forces into a ground domain(Fig.4). Some objects of interest,e.g.a pipe,have geometries that allow defining a smooth contour by mathematical parametrization of the wall geometry; in other cases,several walls may be needed to define the object,e.g.a screw pile.However,the object is defined,and the mathematical model of the normal contact constraint is added to the linear momentum balance equations, as follows.

Fig. 4. Contact forces (red vectors) of deformable domain particles contacting with a rigid wall.

At each rigid wall surfaces,an outward normal direction(Fig.5)for any contacting particleIis defined asnI.When a particleIis in contact with one of the surfaces,the normal gapgNis calculated by the projection of the exceeding distance vector into the corresponding normal direction:

Fig. 5. Kinematics of a particle contacting with a rigid wall.

whereC= 2π, andx1is the radial coordinate in axisymmetric conditions, whereasC= 1 andx1= 1 in other conditions.

The contact constraint adds a nonlinear term to the governing equations,which requires linearization for integration purposes.A simplified general expression is used in this case, giving the definition of the contact stiffness matrix as

The tangential contact constraint is generally modeled as cohesive-frictional in G-PFEM. For a purely frictional case, the tangent vector,the tangential gap and contact force can be written as follows:

An elastoplastic analogy can be applied to defining the tangential counterpart of the contact force. Further details on this approach can be found in Monforte (2018).

3.3. The CDM

The previous approach minimizes the complexity of the contact mechanics problem and gives a proper solution for structures in which deformation can be neglected. A different, more general,approach to define contact interaction in G-PFEM starts from the premise that all contacting domains are deformable. Traditional contact methods can be applied in this case,considering that,at the computing stage,a finite element mesh is used.However,there is a technique that was specially designed for PFEM, i.e. the CDM(Hartmann et al.,2009;Oliver et al.,2009).This technique embeds contact detection within the remeshing stage of the general PFEM algorithm. In the CDM, an ancillary contact mesh is created upon remeshing (Fig. 6). That contact mesh serves two purposes: (i)detection and definition of the contact surface and(ii)computation of the contact constraint term.

Fig. 6. Ancillary triangular mesh used for contact detection and computation in the CDM.

The “contact-domain” ancillary mesh is generated through a constrained Delaunay tessellation of the exterior of the deformable domains after a boundary shrinkage operation.For each one of the resulting triangular ancillary elements,the normal vectornand the tangent vectortare defined by its base.As illustrated by the contact triangle of Fig. 6, contact gaps are defined using the relative position and movement of the vertex not in the base (node 3) of the ancillary contact-domain element. Knowing the magnitudes at configurationn, the updated definition of the normal and tangential gaps is given by

where∇·uis the material gradient of the incremental displacement field of the two contacting bodies and the projection of the gradient to a given directionnis (∇·u)·n= ∂u/∂n. The evaluation of this term can be reduced to a pure geometrical problem (Hartmann et al., 2009). From this definition, a local element-wise constraint enforcement is applied to defining normal and tangential discrete Lagrange multipliers as

wheretNandtTare respectively the normal and tangential components of the traction vectort= σ·nresulting from the stresses of the base-side domain element; and τ is the element-wise constant stability parameter obtained as

whereEminis the minimal Young’s modulus of the contacting bodies;Lis the base-side length of the contact domain element(Fig.6);and αstab∊[1，2]is a dimensionless user-defined parameter,independent of the mesh size and introduced to reduce the over constraint derived from the non-smoothed definition of boundaries with linear finite elements.

After enforcing the contact constraints via the determination of the discrete Lagrange multipliers, the resulting contact contributions can be computed. Therefore, the discretized contact virtual work can be expressed using the introduced approximations for each four nodded contact patch as

corresponding to the normal and tangential counterparts.Using the of the discretization of the displacement incrementuh=Nu· ~u，vectorsNn，Tn,NtandTtare defined as

The directional derivatives appearing in the last expressions can be evaluated using the geometrical characteristics of the contact domain. Complete expressions for the linearization of the contact constraints are developed in detail in Hartmann et al. (2009).

4. Representative numerical simulations

In this section, four examples are presented to illustrate the capacities of G-PFEM. The first set of simulations, i.e. the indentation of a deformable footing in a porous, soft rock, is used to showcase the mesh-independence properties of the formulation for both drained and undrained conditions. The two contact strategies outlined above are compared in the second analysis which describes a pipeline seabed interaction problem. In the third example, i.e. the uplift of a screw pile, the effect of structural deformability on model response is examined. Finally, the last simulation illustrates a dynamic application by examining the free fall of a sphere into a clayey material. Except where otherwise indicated, the soil-structure contact is assumed smooth in all the simulations presented.

4.1. Indentation of a rigid strip footing

The first numerical analysis corresponds to the indentation of an elastic,strip footing in a porous,soft rock.A parametric study of this problem,focusing on the effect of permeability and using different numerical settings,was presented in Monforte et al.(2019a).It was shown there how the failure mode was switched from diffuse contractive failure to shear-localized bearing capacity failure as permeability decreased and undrained behavior emerged.

Here only two permeabilities are considered,k= 10-10m/d andk= 0.1 m/d, corresponding to undrained and drained conditions,respectively.The footing has a height of 0.5 m and a width ofB= 1 m. Because of the problem symmetry, only half of the footing is represented. The domain expands 5Bin the horizontal direction and 3:8Bin the vertical direction. Null excess water pressure is prescribed on the upper free surface and at the bottom of the domain. Null displacements on both directions are prescribed at the bottom of the domain, whereas null horizontal displacements are prescribed on the left and right boundaries. A constant velocity, equal toB/2 a day, is prescribed on the upper boundary of the footing. A sketch of the geometry and the considered boundary conditions are depicted in Fig. 7. The constitutive parameters employed for both footing and soft rock are reported in Table 1.The characteristic length of the non-local model is set tolc= 0:025B.

Table 1 Soil parameters in the example simulations.

Fig. 7. Sketch of the geometry for indentation of a strip footing.

During the computation of this problem,h-adaptive techniques are used: new nodes are inserted in regions where plastic flow is large. In particular, new nodes are inserted in regions where the nonlocal plastic deviatoric strain is higher than a threshold, until the minimum element size attains a predefined length,lmin. As a result, coarse elements are used in regions still in elastic regime

whereas a fine discretization is applied to zones of intense plastic straining. The problem is computed repeatedly, setting three different minimum values of element size as different fractions of the non-local characteristic length, i.e.lmin= 0.4lc,lmin= 0.1lcandlmin=0.075lc.These different mesh refinement controls ensure that the number of neighboring Gauss points from which non-local variables are computed is different in the three cases and therefore represent a strong test on the mesh-insensitive nature of the algorithm. Monforte et al. (2019) illustrated the meshindependence of this formulation in simplified cases (biaxial shearing) where the number of elements and nodes in each simulation was kept constant and all the elements of each mesh had a similar size.In the more complex problem analyzed here,the focus is on the effect ofh-adaptive techniques and different element-sizes on the mesh-dependence of the solution.

Fig. 8 presents the footing penetration curves for drained and undrained conditions for the three levels of initial mesh refinement. It is noticeable how undrained conditions result in a somewhat stiffer elastic response. As noted in Monforte et al. (2019a),although structural yielding takes place at about the same load level, the response afterwards is very different for drained and undrained cases. For the drained condition of diffuse failure, a smoothly hardening structural response is predicted and the three curves practically overlap,showing minimal effect of variable mesh refinement.In undrained conditions,the footing pressure plateaus after yielding, but then a sudden drop of resistance appears at 0:14B.Monforte et al.(2019a)noted that this fall corresponds to the moment in which the shear bands define a fully formed plastic mechanism.There is a slight effect of the different mesh refinement settings on that response, with the coarsely-limited mesh predicting a less abrupt post-plateau capacity fall and mobilizing a somewhat larger residual capacity.

Fig. 9 reports the accumulated plastic shear strain at the final simulation stage for undrained conditions. The localized shear bands define a bearing capacity failure mechanism with minimal effect of the degree of mesh refinement. Fig. 9 also presents the final finite element mesh. As expected, the mesh adaptation routines maintain a coarse mesh in areas in elastic regime whereas a fine discretization appears in the regions where strain localization takes place.

Fig. 8. Rigid strip footing load-displacement curve for drained and undrained conditions for the three levels of discretization.

4.2. Vertical and horizontal displacement of a pipeline

The second numerical example compares the two contact formulations implemented in G-PFEM, i.e. the one using a penalty method for rigid walls and the more general one based on the CDM.A plane strain simulation of vertical embedment and subsequent lateral displacement of a pipeline is employed to illustrate their comparative performance.The pipeline has a diameter of 1 m and is initially located just in contact with the seabed.It is then advanced vertically at a rate of 1 m/d and then fixed vertically and displaced horizontally, again at 1 m/d rate. The simulation domain is a rectangle of 14 m wide and 7 m deep. The pipe initially contacts the midpoint of the upper boundary.Null displacements are prescribed at the bottom of the domain and null horizontal displacements are prescribed at the vertical boundaries (Fig.10).

The constitutive parameters of both the soil and the pipeline are presented in Table 1. There is no structure in this material (pt= 0)but the isotropic preconsolidation pressure is relatively high(ps= 400 kPa), resulting in overconsolidation ratios of around 20.For these conditions, the model predicts - for shearing in triaxial compression - a slightly brittle undrained strength (peak/critical=1.5)that is practically uniform in depth.The permeability and loading rate selected ensure undrained behavior. The comparison here is focused on the relative performance of the contact algorithms and not on the structure deformability,thus in the CDM simulation, the pipeline is represented by a full section of a high stiffness elastic material.

Fig.11a presents the normalized penetration resistance during the vertical motion of the pipe. Almost coincident results are obtained for the two different contact formulations. The penetration curve (vertical force normalized by the residual undrained strength) is very similar to the backbone curve presented by Chaterjee et al. (2012a) to summarize a series of ALE simulations with a softening rate-dependent Tresca model.The main difference lies in the initial slope, which reflects the lower soil stiffness employed in this simulation(the rigidity indexE/suis around 50 in this simulation, whereas it was 500 in Chaterjee et al. (2012a),beingsuthe undrained strength of the material).

Fig. 9. Final finite element mesh (upper) and accumulated plastic deviatoric strain (lower) for the three levels of discretization: (a, d) Coarse, (b, e) medium, and (c, f) fine.

Fig.10. Sketch of the geometry for displacement of a pipeline.

Fig. 11b presents the evolution of the horizontal and vertical forces acting on the pipeline during its motion. Because of the prescribed purely horizontal motion, the vertical force drops slightly more than 50%once the lateral displacement starts.At large horizontal displacements,the so-called“friction ratio”of horizontal to vertical force acting on the pipe rises to 1.1.These results are well aligned with those obtained in previous ALE simulations of the same problem (Chatterjee et al., 2012b). Nevertheless, the main point of interest here is that they seem independent of the contact algorithm employed in the simulation.

Fig.12 depicts the excess water pressure and accumulated plastic shear strain at the end of the problem.Plastic strain delineates shear localization features.Shear localization takes place at the soil-pipe contact,but also indicates a vertical bearing capacity failure mechanism at the initial pipe location and a series of passive wedges crossing the soil that accumulates in front of the pipe. There are some differences in the number and development of the lateral wedge localization features between the two simulations; these may be attributed to the slightly different contact forces due to the different formulations.The excess pore water pressure is negative in the area vacated by the pipe as a result of unloading,but also in the displaced passive wedge as a result of dilatant shearing.Again,some subtle differences between the two contact formulations appear in the wedge area.These small differences in the development of the localization patterns do not affect the structural response.

Finally, Fig. 13 illustrates the effect of contact friction. The increased effect on pipeline vertical resistance of contact friction agrees with previous work on this topic. Importantly, the example shows how the very similar performance of the different contact algorithms is maintained also in cases inwhich the contact is frictional.

4.3. Uplift capacity of a screw pile

Screw piles, generally made of high strength steel, consist of a shaft with a number of helices attached to it. Screw piles require less installation effort as they are penetrated into the ground by applying a turning moment at the head of the pile (Mohajejarni et al., 2016), and are frequently employed as anchors. Monforte et al. (2019c) investigated the pull-out capacity of single and double-helix screw piles using a Tresca single-phase model and rigid-wall model for the pile. Here the problem is revisited in a coupled simulation using CDM and the structured soil constitutive model.

The problem is simplified by assuming axisymmetric conditions.The pile shaft has 1.7 m length and 0.05 m radius and helices have 0.2 m radius. Two cases are examined with a single helix or two helices, spaced 0.2 m apart (Fig. 14). Linear elastic behavior is assumed for the pile (Table 2). The constitutive parameters of the soil are reported in Table 1.The base case corresponds to a slightly over-consolidated clay,with no structure(pt= 0).The installation phase is not modeled and the pull-out is represented by applying a vertical velocity of 1 radius per hour at the head of the pile.Due to the low permeability considered, undrained conditions prevail.

Table 2 Structure parameters in the example simulations.

Fig. 15a shows the evolution of the pile resistance in terms of normalized pile uplift for different values of pile stiffness.Although all curves converge at large displacements, stiffer piles mobilize higher soil resistance faster as they displace. This effect appears almost identical for the case with one helix and for the case with two helices (Fig.15b). The failure mechanism for a single helix is one of bearing capacity ahead of the pile to which,in the case of a double helix,a cylindrical shear surface connecting the helix tips is added(Fig.16).This cylindrical shear surface is less sensitive to the helix wing deflection.

The study of stresses in the structure is also possible thanks to the use of the contact domain.To showcase this,Fig.17 reports the vertical stress on the helices at an uplift of 0.15 m.The stiffer screws present larger vertical stress as they mobilize a larger resistance.Also, the vertical stress is much larger in the upper helix with respect to the lower one.

Fig.11. Effect of contact formulation(rigid:object as moving wall;def:CDM)on a pipeline insertion and lateral displacement simulation:(a)Vertical resistance versus normalized depth, and (b) Evolution of the vertical and horizontal soil-pipeline contact forces.

Fig.12. Pipeline insertion and lateral displacement simulation:(a,b)Accumulated plastic shear strain,and(c,d)Excess water pressure(kPa)for pipeline as prescribed rigid wall(a,c), and pipeline as deformable body (b, d).

4.4. Free falling spherical penetrometer

Free falling spherical penetrometers have been recently developed as a rapid low-cost tool for offshore site investigation(Morton and O’Loughlin, 2012). The instrument is dropped overboard and penetrates the seabed employing the kinetic energy acquired by free falling through the water column. The penetrometer motion during deployment is recorded by inertial sensors and later interpreted to infer soil undrained strength (Morton et al., 2016). A comprehensive ALE based numerical parametric study of the problem has been recently presented (Mana et al., 2018). The problem is here used to illustrate the capabilities of G-PFEM in a dynamic setting.

Fig.13. Effect of contact friction on pipeline insertion curve.

Fig.14. Sketch of the geometry for screw pile.

Fig.15. Effect of structural stiffness on the pull-out curve of single and double helix screw piles:(a)Pull-out forces,and(b)Pull-out forces at reduced stiffness normalized by that at higher stiffness.

The sphere diameter is 20 cm and impacts vertically on the soil.An axisymmetric model is adopted, the simulation domain is a square measuring 2 m by 2 m in size. Null displacements are prescribed on the bottom boundary of the soil, whereas null radial displacements are prescribed at the outer radius. Drainage is allowed at the top and bottom boundaries.

Fig. 18 presents a sketch of the geometry of the problem. The constitutive parameters of the soil are reported in Table 1: they describe a rate-independent, slightly over-consolidated unstructured soft clay. Because of the low operative stiffness of the soil skeleton,the accelerations of fluid and solid skeleton are presumed equal and the use of a full Biot formulation(Monforte et al.,2019b)is not necessary. The properties of the penetrometer are listed in Table 2; they correspond to an elastic, relatively light-weight and low-stiffness object,different from the steel penetrometers that are more frequently employed.

The dynamics of the impacting sphere are illustrated in Fig.19 for a series of analyses in which the impact velocity is varied.The sphere suffers a very intense deceleration (up to 20g) that, almost independently of initial velocity, stops penetration in about 30 ms.Because the impacting sphere is elastic no remeshing takes place on it and the curves are recorded in a node at the top of it;acceleration curves were slightly smoothed to avoid the clutter caused by elastic waves on the sphere. A small rebound is visible at the end of the impact.This rebound reflects elastic energy stored in the soil(Zhang et al.,2021)and,as observed also by others(Zambrano-Cruzatti and Yerro,2020),can be reduced by increasing the elastic modulus.

Penetration curves for two different series of impacts are presented in Fig.20 showing how the effect of increased impact velocity or increased sphere mass increases penetration depths. Both variables (sphere mass or impact velocity) have a direct effect on the kinetic energy at impact: plotting this variable versus penetration depth, a clear relationship appears (Fig. 21a). This agrees with experimental (O’Loughlin et al., 2013) and numerical (Zhang and Evans, 2019) studies of dynamic offshore anchors showing a direct relation between total available energy at impact and penetration depth. A similar relation also holds for other procedures involving soil impact,such as the standard penetration test(Zhang et al.,2021).In general,total available energy at impact includes kinetic and potential terms,but the small penetration achieved in the ball impacts simulated here makes the potential term negligible.

Mana et al. (2018), in their parametric numerical study of spherical penetrometers, elaborated on the impact energy versus depth relation.They defined a dynamic bearing capacity factor,Nd,to transformEs(z), the fraction of impact energy spent by the impactor at depthz, into plastic work as

Fig.16. Effect of structural stiffness on screw pile pull-out.Case with two helices(a,b)and case with a single helix(c,d).Plastic strain at normalized displacement u/D=0.15 under(a, c) Young’s modulus of the pile E = 0:5 GPa, and (b, d) E = 500 GPa.

Fig.17. Effect of structural stiffness on screw vertical stress(kPa)at uz=0.15 m:(a,c)Young’s modulus of the pile E=0.5 GPa,and(b,d)E=500 GPa.Case with a single helix(a,b)and case with two helices (c, d).

Fig.18. Sketch of the geometry for free falling sphere.

whereApis the projected area of the sphere in contact with the soil.They went on to obtain the following envelope for the dynamic bearing capacity factor at maximum penetration depth:

wheredpis the penetration depth,Dis the sphere diameter,andAis a coefficient that depends on soil viscous properties and is 5.5 for a non-viscous case. Fig. 21b compares the results of the G-PFEM simulations with the Mana et al. (2018) envelope. The values ofsuemployed correspond to those obtained in triaxial compression after isotropic compression to confining stress of 10 kPa.It appears that,despite using different numerical method,constitutive model and sphere characteristics,the G-PFEM results are very close to the previously obtained envelope.

Fig.19. Time evolution of (a) vertical displacement, (b) velocity, and (c) acceleration for free falling sphere.

Fig. 20. Penetration curves of free falling sphere for (a) different impact velocities at fixed mass (m0) and (b) for different sphere masses at fixed impact velocity (3 m/s).

5. Conclusions

This work has presented some characteristics of G-PFEM, a PFEM implementation specially designed for geomechanical problems involving large deformations and soil-structure interaction. PFEM adopts a particle description of the domain, relying on the FEM to solve the balance equations and Delaunay tessellation for remeshing. In G-PFEM, an updated Lagrangian description of motion, using large displacement and finite strain theory, is employed. The governing equations are solved using mixed, stabilized finite elements to bypass volumetric locking, resulting from undrained conditions or from quasi-incompressible response at critical state. Remeshing minimizes the distortion of the mesh while preserving the boundary of the domain.

Fig.21. (a)Relation between impact energy and depth of penetration,and(b)Evolution of dynamic bearing capacity factor for impacts(red curve:limit envelope from Mana et al.(2018)).

Although not strictly necessary, to use all the potential of GPFEM,constitutive models need to be adapted to large strains.This step is illustrated here for a classic structured soil model, reformulated using a multiplicative split of elastoplastic deformation.This model is well adapted to represent brittle responses of materials such as porous rocks or structured clays. Material brittle responses, conducive to strain localization, are frequent in geomechanics.This may lead to pathological mesh-dependence.To deal with this problem, G-PFEM has been equipped with a nonlocal regularization technique.

Two contact formulations are presented, one for the contact between a porous media and a rigid object and another, only recently included in G-PFEM, to simulate the interaction between two deformable objects.In the former,a parametric representation of the rigid object is adopted and contact constraints are imposed by means of a penalty method.The latter relies on a virtual mesh for the interface and Lagrange multipliers to introduce contact constraints.

A set of simulations has been presented to showcase G-PFEM performance. We first demonstrate the mesh insensitive nature of the solution in drained and undrained conditions by examining the insertion of a deformable footing into a soft,porous rock.When the object that interacts with the soil is very rigid, both contact formulations give very similar results: this is illustrated using the interaction of a rigid pipeline with the seabed.However,sometimes deformation of the object is itself of interest or may affect the outcomes of its interactions with the soil,as shown in the case of a screw-pile pull-out.Finally,the dynamic capabilities of the code are illustrated by analyzing the impact of a free-falling ball into a clayey seabed.

To facilitate a wide application of the method, G-PFEM is currently developed as open-source and the code is available in a public repository(https://gitlab.com/pfem-research/kratos/).

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

The authors acknowledge financial support by Severo Ochoa Centre of Excellence (2019-2023) Grant No. CEX2018-000797-S funded by MCIN/AEI/10.13039/501100011033 and research projects BIA2017-84752-R and PID2020-119598RB-I00.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jrmge.2021.12.006.