A stable implicit nodal integration-based particle ﬁ nite element method (N-PFEM) for modelling saturated soil dynamics

In this study, we present a novel nodal integration-based particle ﬁ nite element method (N-PFEM) designed for the dynamic analysis of saturated soils. Our approach incorporates the nodal integration technique into a generalised Hellinger-Reissner (HR) variational principle, creating an implicit PFEM formulation. To mitigate the volumetric locking issue in low-order elements, we employ a node-based strain smoothing technique. By discretising ﬁ eld variables at the centre of smoothing cells, we achieve nodal integration over cells, eliminating the need for sophisticated mapping operations after re-meshing in the PFEM. We express the discretised governing equations as a min-max optimisation problem, which is further reformulated as a standard second-order cone programming (SOCP) problem. Stresses, pore water pressure, and displacements are simultaneously determined using the advanced primal-dual interior point method. Consequently, our numerical model offers improved accuracy for stresses and pore water pressure compared to the displacement-based PFEM formulation. Numerical experiments demonstrate that the N-PFEM ef ﬁ ciently captures both transient and long-term hydro-mechanical behaviour of saturated soils with high accuracy, obviating the need for stabilisation or regularisation techniques commonly employed in other nodal integration-based PFEM approaches. This work holds signi ﬁ cant implications for the development of robust and accurate numerical tools for studying saturated soil dynamics.


Introduction
Saturated soil represents a multiphase material comprising a solid skeleton and fluid occupying voids.The study of saturated soil dynamics has long been recognised as a challenging problem (Zienkiewicz and Shiomi, 1984), necessitating the simultaneous analysis of both solid and fluid phases and their interactions (Zienkiewicz et al., 1980).Time-dependent processes in saturated soils (Bjerrum, 1967;Olson, 1998), such as soil consolidation and the dissipation of pore water pressure, introduce additional complexity for numerical modelling.This complexity arises from the need for numerical schemes with appropriate time-stepping and convergence criteria (Prevost, 1983;Zienkiewicz et al., 1999).
In situations where soils undergo significant deformation, as observed in landslides and foundation penetration processes, numerical modelling becomes even more challenging.This is due to potential issues related to mesh distortion and the evolution of free-surfaces.
Over the past decades, significant advances have been made in the development of numerical methods for modelling saturated soil dynamics under large deformation.These methods include the Arbitrary Lagrangian-Eulerian method (Sabetamal et al., 2014), Smoothed Particle Hydrodynamics method (Pastor et al., 2009;Bui and Nguyen, 2017;Peng et al., 2017), Material Point Method (Zhang et al., 2009;Bandara and Soga, 2015;Yamaguchi et al., 2020;Kularathna et al., 2021), and the Particle Finite Element Method (PFEM) (Monforte et al., 2019;Wang et al., 2021a;Carbonell et al., 2022;Jin and Yin, 2022), among others.In these studies, the coupling between the solid skeleton and the pore fluid is described through either the Biot formulation in a set of governing equations for the solid-fluid mixture (Zienkiewicz and Shiomi, 1984) or through models of interaction forces added to two separate sets of governing equations for the two phases (Drew, 1983).Due to its effectiveness and efficiency, the Biot formulation has been widely adopted in the development of these methods for various practical applications (Pastor et al., 2009;Zhang et al., 2009;Monforte et al., 2019;Wang et al., 2021a;Jin and Yin, 2022).While interaction models can capture more complex coupling processes between two-phase materials, the computational cost is higher since the solid and fluid phases are separately simulated (Bandara and Soga, 2015;Bui and Nguyen, 2017).Among these developments, explicit time integration schemes with conditional stability are most commonly used, although they are not computationally efficient for analysing geotechnical problems characterised by low to medium frequency.Recently, semi-implicit and implicit time integration schemes have been implemented in these methods to improve the computational efficiency (Kularathna et al., 2021;Wang et al., 2021a;Carbonell et al., 2022;Lian et al., 2023).
The PFEM is a hybrid method that combines finite element analysis with the feature of particle approaches.Originally developed in the fluid dynamics community to model free-surface flow and fluid-structure interaction problems (Oñate et al., 2004), the PFEM has gradually gained recognition in geotechnical community as well.Various versions of the PFEM have been developed to tackle large deformation problems in geotechnical engineering applications.The PFEM has become popular due to its two key features: (i) it naturally inherits the demonstrated capability of the FEM in describing the complex nonlinear behaviour of geomaterials, and (ii) pre-existing FEM codes can be easily extended to incorporate the particle technique (Oñate et al., 2004;Cremonesi et al., 2011;Zhang et al., 2013) for large deformation analysis.
Among recent PFEM developments, one important improvement is the adoption of a nodal integration scheme (Zhang et al., 2018).This scheme eliminates the mapping requirement after the re-meshing operation in modelling history-dependent materials, as all state variables are stored at mesh nodes.Remarkably, in most developed nodal integration-based PFEM versions, it has been shown that numerical instabilities may arise in the dynamic analysis of both single-phase (Jin et al., 2021a;Yuan et al., 2023a) and two-phase materials (Jin and Yin, 2022;Wang et al., 2023a;Yuan et al., 2023b) if no stabilisation technique is imposed.For singlephase materials, instability arises from low-order integration, known as temporal instability, which has also been observed in nodal integration-based smoothed FEM (Zhang and Liu, 2010) and various other mesh-free methods (Belytschko et al., 2000;Chen et al., 2001).Concerning two-phase materials, the additional infsup condition (Franco and Bathe, 1990) should be fulfilled since equal order shape functions are used.The introduction of stabilisation may involve different techniques or tunable parameters for multiple problems, for which no universal criteria are available (Belytschko et al., 2000;Puso et al., 2008;Silva-Valenzuela et al., 2020;Wei et al., 2020).Furthermore, in some developed nodal integration-based PFEM versions (Jin and Yin, 2022;Yuan et al., 2022), explicit or semi-implicit schemes were adopted, which may be less efficient in simulating long-term geotechnical problems, such as consolidation settlement in clay.
To circumvent the use of stabilisation techniques and develop an implicit PFEM version, the second-order cone programming (SOCP) is considered as a promising solution.The adoption of mixed variational principles, consolidating all the governing equations into a single functional, can achieve better accuracy in stress calculation, compared to displacement-based FEM (Reddy, 2002).Previous studies have demonstrated that an accurate and efficient implicit time integration scheme is compatible with the SOCPbased computational framework (Zhang et al., 2013(Zhang et al., , 2019(Zhang et al., , 2023;;Wang et al., 2021aWang et al., , 2023b;;Wang and Lei, 2023).Recently, the temporal stability of an implicit nodal integration-based PFEM (N-PFEM) using SOCP has been validated for dynamic total stress analysis of single-phase materials (Zhang et al., 2023), and in the dynamic effective stress analysis of saturated media (Wang et al., 2023c) without the need for ad-hoc stabilisation techniques.However, the N-PFEM proposed in the work of Wang et al. (2023c) is limited to simulating saturated soils in two extreme conditions, i.e., fully undrained and fully drained circumstances.In other words, it can only be applied to the limiting scenarios, i.e., shortterm analyses of saturated soils with very low permeability and long-term analyses where excess pore water pressure is fully dissipated.It cannot capture water transport in porous media and its sequential effects on soil deformation.Water transport plays an important role in many geotechnical problems.One example is soil consolidation settlement.Another noteworthy scenario involves the initiation of catastrophic failure, exemplified by the determination of the timing and mechanisms of landslides (Iverson et al., 1997;Torres et al., 1998;Berti and Simoni, 2005).
In this study, we further extend the N-PFEM for the coupled dynamic analysis of saturated media, considering water seepage.This development is based on the generalised Hellinger-Reissner (HR) variational principle as proposed in the work of Wang et al. (2021a).The coupled N-PFEM developed in this study is assessed by studying a range of challenging geotechnical problems, such as wave propagation, strip footing, and consolidation, which require the fulfilments of temporal stability and the inf-sup condition.It is demonstrated that, while inheriting the advantages of the conventional nodal integration-based PFEM for modelling large deformation problems, the developed N-PFEM also successfully reproduces smooth effective stress and pore water pressure fields without the need for any stabilisation technique.

Model description
In this section, the governing equations and solution schemes of the current numerical model are briefly introduced.More theoretical and technical details can be found in our previous studies (Zhang et al., 2013(Zhang et al., , 2019(Zhang et al., , 2022;;Wang et al., 2021a,b).

Governing equations
A two-dimensional (2D) domain, V, comprising saturated porous media and bounded by a surface S, is considered, as shown in Fig. 1.
The governing equations for dynamic saturated media are below (Zienkiewicz et al., 1999;De Boer, 2012): (a) The equilibrium equation for the mixture Fig. 1.Domain of a saturated medium and its boundary partition.S t , S u , S p and S q represent the associated surfaces subjected to traction, displacement, pore water pressure, and fluid flux, respectively.
where s 0 is the effective stress; m is the vector, m ¼ [1,1,0]; p is the pore water pressure; b is the body force; r is the density; V is the linear strain-displacement operator; and v is the velocity.
(b) Displacement-strain relation where 3 is the strain and u is the displacement.
where r f , g f , and b f are the density, the unit weight, and the body force of a fluid, respectively; k is the Darcy hydraulic conductivity; and w is the superficial velocity.
(d) Mass balance equation where w can be eliminated by substituting Eq. (3) into Eq.( 4): where N is a matrix consisting of components of the outward normal vector of the boundary; and t, u p , p p and q p are the prescribed traction, displacement, pore water pressure and fluid flux, respectively.
(f) Constitutive equations of elastoplasticity where F is the yield function; 3 e and 3 p are the elastic and plastic strains, respectively; C is the elastic compliance matrix; l is the plastic multiplier; and G is the plastic potential.The associated flow rule assumes G ¼ F. The incremental forms of Eqs. ( 10) and ( 11) with their complementary conditions can be expressed as The above governing equations are the so-called u-p form suitable for dynamic analysis of geotechnical problems of low to medium frequency (Zienkiewicz et al., 1999).The porosity n can be updated by _ n ¼ ð1 À nÞV$v, assuming the solid skeleton is incompressible (De Boer, 2012).The density of the mixture is where r s is the density of the solid.

Time discretisation
Using the q-method (Wood, 1990), the effective stress, velocity, and pore water pressure are discretised in time as where the subscripts n and nþ1 denote the known and unknown states; Dt is the time step; and q 1 , q 2 and q 3 are the numerical parameters ranging from 0 to 1. Substituting Eqs (13)e(15) into the equilibrium Eq. ( 1) and boundary conditions Eqs. ( 6) and ( 9) gives Wang et al. (2021a) demonstrated that the above time marching scheme is unconditionally stable for dynamic analysis of saturated media if q 1,2,3 !0.5.In this study, we use q 1 ¼ q 2 ¼ q 3 ¼ 1 unless otherwise specified.

Spatial domain discretisation with nodal integration
Three-node triangular elements are adopted for finite element discretisation, with a nodal integration scheme implemented over smoothing cells.This scheme has been demonstrated to be effective in addressing volumetric locking issues associated with linear elements (Zeng and Liu, 2018;Zhang et al., 2018;Meng et al., 2020).Following standard finite element notations, the displacement field u is approximated by where the hat symbol ' ^' represents variables at mesh nodes, and N u is the matrix containing the linear interpolation shape functions.
The strain tensor 3 over finite elements is then expressed as Due to the linear shape function, the strain tensor ε is uniform within each three-node triangular element.Hence, the smoothed strain at each mesh node can be evaluated as a weighted average of strain over the smoothing cell (Chen et al., 2001;Liu et al., 2007;Zeng and Liu, 2018;Zhang et al., 2018;Meng et al., 2020).As illustrated in Fig. 2, cells are constructed by connecting the centroid of each triangle to its three mid-edge nodes.The smoothed strain at the kth node (or the strain across the kth smoothing cell) is where A s k is the area of the kth smoothing cell U s k , N s is the number of elements contributing to the kth smoothing cell, and A e i , B e u;i and b u e i are respectively the area, strain gradient matrix and nodal displacement of the ith linear triangular element involved in the kth smoothing cell.The relationship between the smoothed strain and nodal displacement at the kth node can be further expressed as (Meng et al., 2020): u;i being the smoothed strain gradient matrix.
Other master fields, such as the effective stress s 0 , the pore water pressure p, and the dynamic force r, are also assumed uniform across each smoothing cell s 0 z s 0 ; rzr; pzp; czc (24) where the overbar represents the value over the smoothing cell (also the value at the mesh node associated with this smoothing cell).For brevity, an intermediate variable c (Zhang et al., 2016;Wang et al., 2021a).

Optimisation problem
According to the generalised Hellinger-Reissner variational principle proposed by Wang et al. (2021a), the discretised governing equations for saturated soil dynamics are equivalent to the following min-max problem: The equivalence between this min-max problem and the governing equations can be verified by taking functional derivatives of the associated Lagrangian functional (Wang et al., 2021a), as shown in Appendix A. However, it is noteworthy that, in this study, all the matrices in the min-max problem (25) are calculated over smoothing cells, which differs from those in the work of Wang et al. (2021a), where they were estimated on finite elements.This difference arises due to the utilisation of nodal integration in this study.Specifically, the matrices and vectors in the min-max problem (25) are where I is an identity matrix, and U s is the smoothing cell. min The min-max problem ( 25) can be further reformulated as a standard SOCP problem with different yield criteria, such as Tresca, MohreCoulomb and Drucker-Prager functions (Zhang et al., 2013(Zhang et al., , 2019)), and the implementation details have been documented by Wang et al. (2021b).In Appendix B, the reformulated SOCP problem is presented.

N-PFEM
The PFEM treats mesh nodes as free particles within the framework of Lagrangian FEM.The evolving computational domain is identified using the alpha-shape method (Edelsbrunner and Mücke, 1994;Oñate et al., 2004;Cremonesi et al., 2010): for a cloud of points with a characteristic length h, an examination is performed on each point to check whether it is possible to place a sphere with a radius ah (a is a pre-defined factor usually varying from 1.2 to 1.6) such that it only contains that point; if possible, the point is a boundary point; otherwise, it is an internal point.During the PFEM simulation, the boundary of the computational domain associated with a new mesh is obtained by removing unnecessary elements from the mesh generated through Delaunay triangulation of updated mesh nodes (Cremonesi et al., 2010).However, the quality of the resulting mesh is usually insufficient for classical FEM analysis to handle large deformation.Therefore, a re-meshing or mesh smoothing operation on the identified domain is needed to ensure mesh quality (Zhang et al., 2013;Cante et al., 2014;Meduri et al., 2019;Wang et al., 2022).Alternatively, the smoothing technique employed in smoothed finite element methods can alleviate this issue by using cells rather than finite elements for spatial domain discretisation (Yuan et al., 2019;Zhang et al., 2018Zhang et al., , 2021Zhang et al., , 2022)).
In this paper, a nodal integration scheme is introduced into PFEM for poro-elastoplastic analysis of saturated porous media.The basic calculation procedures of the developed N-PFEM within each time step are as follows: (i) A cloud of particles (vertices) is used to represent the computational domain (Fig. 3a); (ii) The alpha-shape method is employed to identify the computational domain, which also generates a new mesh through Delaunay triangulation (Fig. 3b); (iii) Smoothing cells are constructed for each node of triangular elements (Fig. 3c); (iv) Incremental finite element analysis is performed; (v) Field variables at all nodes are updated (e.g. the coordinates of particles are updated using the resolved incremental displacements, as shown in Fig. 3d).
In procedure (ii), mesh refinements, involving actions such as node deletion and insertion, are necessary to ensure a high-quality mesh, as indicated by Meduri et al. (2019).These measures are particularly essential when addressing situations with profound deformations, such as those observed in flow-like landslides, breaking waves, and splashes.Additionally, there is the possibility of introducing additional nodes to enhance the precision of numerical simulation.It is worth noting that these refinements could potentially lead to a modification in the total number of nodes, necessitating variable mapping procedures for the newly added nodes.
In the proposed N-PFEM, we adopt the infinitesimal strain (Eq.( 2)) with an updated geometry for large deformation analysis.The use of infinitesimal strain for large deformation analysis has proven to be effective in geotechnical applications (Hu and Randolph, 1998;Tian et al., 2014;Wang et al., 2015;Kong et al., 2018).Notably, the PFEM with this strategy has demonstrated its effectiveness in reproducing large deformational phenomena observed in laboratory experiments.Examples include water dam break flow (Zhang et al., 2019), dry (Zhang et al., 2014) and saturated (Wang et al., 2021a) granular collapses, and submarine landslide-generated waves (Zhang et al., 2019(Zhang et al., , 2022)).In these cases, good agreements were achieved between numerical simulations and experimental measurements.Its effectiveness will be further demonstrated for dealing with saturated soil under large deformation through several benchmarks to be presented in Section 3.

Numerical examples
To validate and demonstrate the proposed model for dynamic analysis of saturated media, numerical simulations are conducted for four benchmarks.These benchmarks cover both elastic and elastoplastic scenarios as well as both small and large deformations.The first two examples involve simulations where node positions are not updated, representing small deformation cases.In contrast, the third and fourth examples utilise the complete N-PFEM procedures, as outlined in Section 2.4.2, to perform the large deformation analysis.

Wave propagation
In this example, we study a typical benchmark for the dynamic analysis of saturated media (Markert et al., 2010): the wave propagation in a 2D rectangular saturated weightless poroelastic medium.The model setup is shown in Fig. 4a, where a surface load  with a width W f ¼ 1 m is applied at the top surface of the domain.The load obeys the function f(t) ¼ 100sin(25pt)(1ÀH(tÀ0.04)) kPa, where H denotes the Heaviside step function, and lasts for 0.04 s (Fig. 4b).Material parameters used in this study are the same as those in the work of Markert et al. (2010): the density of solid r s ¼ 2000 kg/m 3 , the density of fluid r f ¼ 1000 kg/m 3 , porosity n ¼ 0.33, Young's modulus E ¼ 14.5 MPa, Poisson's ratio n ¼ 0.3, and the Darcy hydraulic conductivity k ¼ 10 À2 m/s.The computational domain is discretised using 16,586 nodes and 32,643 three-node triangular elements, based on which 16,586 smoothing cells are constructed.The time step adopted in the simulation is 0.005 s, facilitating an accurate representation of wave propagation over a duration of 0.2 s.Consequently, a total of 400 simulation steps are executed.
Fig. 4c and d presents the evolution history of the displacement at point A and the pressure at point B obtained from our simulation (SOCP-L).These results are also compared with the reference solutions from the work of Markert et al. (2010), the standard nodal integration-based smoothed PFEM without stabilisation techniques from the work of Wang et al. (2023a) (denoted as NS-PFEM), and the simulation results obtained from the same SOCP-based FE computational framework reported by Wang et al. (2021a) using quadratic displacement/linear stress mixed elements (SOCP-Q).Satisfactory agreements are observed between SOCP-Q, SOCP-L, and the results of Markert et al. (2010), demonstrating the accuracy of the proposed nodal integration scheme for the dynamic analysis of saturated media.As for the NS-PFEM, incorrect results (i.e.wrong shape of the displacement at A and the oscillation of pore water pressure at B) are obtained due to the temporal instability-induced non-physical oscillations, as reported by Wang et al. (2023a).The slight discrepancy in the displacement obtained from SOCP-L, as compared to SOCP-Q and the results from the work of Markert et al. (2010), can be attributed to the utilisation of shape functions one order lower in the current method.It is worth noting that SOCP-Q (Wang et al., 2021a) employs quadratic/linear mixed triangular elements, resulting in a higher number of mesh nodes (i.e.65,814).Consequently, SOCP-Q exhibits lower computational efficiency compared to the nodal integration-based method developed in our study (i.e.SOCP-L).This is evident from the computational times of 240 min for SOCP-L and 759 min for SOCP-Q in solving this benchmark problem, using a Windows 10 Desktop with an Intel i7 CPU @ 3.00 GHz.
To further demonstrate the numerical stability of our proposed method for dynamic analysis of saturated media, wavefields of pore water pressure at t ¼ 0.05 s are illustrated in Fig. 5, where the displacement is scaled by a factor of 250.The simulation results obtained from the developed SOCP-based method using linear elements and nodal integration (referred to as SOCP-L) are compared with those from the conventional FE simulation using linear elements, both with and without stabilisation as described by Monforte et al. (2019).As can be seen, unphysical oscillation is present in the conventional FE simulation results when linear elements are used without stabilisation techniques (Fig. 5c).In contrast, our method demonstrates the ability to generate a smooth wavefield, as depicted in Fig. 5a, even in the absence of any stabilisation techniques on linear elements.While for the simulation of this problem using the NS-PFEM model (Wang et al., 2023a), severe numerical noises related to temporal instability are observed if stabilisation techniques are not applied.Notably, the results obtained from our method exhibit good agreement with those using the same SOCP-based formulation with quadratic/linear mixed elements (Fig. 5b) and conventional FE simulation with stabilisation techniques (Fig. 5d).

Strip footing and long-term consolidation
In this section, the long-term consolidation of both elastic and elastoplastic saturated soil subjected to footing under small deformations is considered.This problem has also been examined in various previous studies (Jin et al., 2021b;Yuan et al., 2023b).The model setup is as shown in Fig. 6a, consistent with those in the work of Manoharan and Dasgupta (1995) and Yuan et al. (2023b).Material parameters are as follows: the density of solid r s ¼ 2000 kg/m 3 , the density of fluid r f ¼ 1000 kg/m 3 , porosity To facilitate the analysis, a scaled value of time, T v , is introduced, where T v ¼ 1 corresponds to 4680 d.The load with a width of 3 m is applied at the top surface and described by a ramp function over a period of 0.01T v and then is kept constant, as depicted in Fig. 6b.The simulation is conducted for a duration of 1000T v , and it consists of two stages (Yuan et al., 2023b): the first loading stage (from t ¼ 0 to t ¼ 0.01T v ) is simulated using 100 constant time steps, while the long-term consolidation stage extends until t ¼ 1000T v with a constantly increased time step Dt by a factor of 1.1 at each step.
Following the work of Manoharan and Dasgupta (1995), the displacements are normalised concerning the width of the load W f , and pore water pressure is normalised concerning the maximum value of the load, i.e. 100 kPa.Fig. 7a illustrates the simulation results for the elastic behaviour of the saturated porous media.It shows that Points A and B experience rapid settlement during the loading phase (t 0.01T v ), followed by gradual deformation due to consolidation (0.01T v t 100T v ), and finally reach a steady state after complete dissipation of pore water pressure (t !100T v ).The normalised displacements at points A and B agree well with the independent numerical results reported by Manoharan and Dasgupta (1995).Additionally, the normalised pore water pressure recorded at point B is nearly consistent with the results reported by Manoharan and Dasgupta (1995), as shown in Fig. 7b.The results of Yuan et al. (2023b) are not presented in this figure since they exhibit a high degree of similarity to the curves depicted in Fig. 7a and b.The same numerical configuration with an elastoplastic medium is further examined.The associated Mohr-Coulomb yield criterion, with cohesion c ¼ 10 kPa and friction angle 4 ¼ 20 , is adopted in this simulation which is in line with that reported by Sabetamal et al. (2016).The normalised displacement records at points A and B and normalised pore water pressure records at point A are presented in Fig. 8a and b, respectively.These results exhibit a high level of agreement with the results of Manoharan and Dasgupta (1995) and Sabetamal et al. (2016).
To demonstrate the accuracy of our model, we present the pore water pressure fields at two specific instants: t ¼ 0.001T v and 0.01T v , which have been previously studied in two separate versions of smoothed PFEM (Jin et al., 2021b;Yuan et al., 2023b).In Fig. 9, numerical results are compared with those obtained from the smoothed PFEM without and with stabilisation techniques (Jin et al., 2021b;Yuan et al., 2023b).The numerical solutions from the present model at the two instants for both elastic and elastoplastic analyses are shown in Fig. 9aed and g.These results exhibit smoother distributions compared to jagged patterns obtained by the smoothed PFEM without stabilisation techniques (i.e.Fig. 9bee  and h) (Jin et al., 2021b;Yuan et al., 2023b).Notably, when comparing Fig. 9b with Fig. 9e and h, spurious oscillations are more pronounced in the elastoplastic analysis for the smoothed PFEM without stabilisation techniques.In contrast, the present model accurately captures the realistic pore water pressure field for both elastic and elastoplastic analyses, similar to the results obtained in the smoothed PFEM with stabilisation techniques (Jin et al., 2021b;Yuan et al., 2023b) as shown in Fig. 9cef and i, demonstrating the numerical stability of our model for elastic and elastoplastic analyses of saturated media.

Self-weight consolidation
In this section, the elastic behaviour of a saturated medium under large deformation is investigated by considering a 2D selfweight consolidation of a block.This configuration has been previously studied as a standard benchmark (Zhao and Choo, 2020;Zheng et al., 2021;Yuan et al., 2022).The elastic block is a rectangle with a width of 4 m and a height of 2 m.Due to symmetry, a square with dimensions of 2 m is adopted in the simulation.The model setup, depicted in Fig. 10a, illustrates the mesh configuration and boundary conditions.The material parameters used in the work of Zheng et al. (2021)   constant time step of 0.005 s is used in this simulation for a duration of 0.5 s.
The numerical results of pore water pressure at three time instants, i.e. t ¼ 0.1 s, 0.3 s and 0.5 s, along with the displacement at t ¼ 0.5 s, are presented in Fig. 10c.Similar results have been observed in previous studies (Zheng et al., 2021(Zheng et al., , 2022;;Yuan et al., 2022), where stabilisation techniques or a fractional step scheme were employed to ensure accuracy.The pore water pressure at three monitoring points, denoted as A, B, and C (as illustrated in Fig. 10a), is compared with the results of Yuan et al. (2022).Remarkably, these comparisons give good agreements, which demonstrate the validity and accuracy of our model for elastic analysis of saturated media under large deformation.

Vertical cut of a saturated soil block
In this section, a typical example of a poro-elastoplastic medium involving large deformation (Sanavia et al., 2002;Navas et al., 2018) is analysed to validate the capability of the proposed model in handling practical geotechnical problems.The model setup of a rigid footing on a square domain of saturated soil, same as that studied by Navas et al. (2018), is shown in Fig. 11.The model parameters are shown in Table 1.
The computational domain is discretised into 31,008 linear triangular elements with 15,758 nodes, such that the total number of cells is 15,758.In the study of Navas et al. (2018), a total vertical displacement of 2 m is imposed on the loading area at a constant rate of 0.02 m/s.In this study, the load is enforced through an incremental vertical displacement of Du y ¼ À0.02 m with a time step of Dt ¼ 1 s, resulting in 100 simulation steps.
The computed pore water pressure field at two instants, t ¼ 50 s and t ¼ 100 s are first compared with available results.Yuan et al. (2023b) employed different values of the stabilisation parameter s to define a stabilising term that is added to the equilibrium equation of fluid mass, to investigate the performance of a developed Polynomial-Pressure-Projection technique in stabilising the pore water pressure field in this challenging example.Since this problem involves large deformation elastoplastic analysis, the spurious oscillations in pore water pressure field are more pronounced compared to both large deformation elastic analysis and small deformation elastoplastic analysis.To alleviate this issue, a large stabilisation parameter s ¼ 100 was used by Yuan et al. (2023b).It should be noted that a different value of stabilisation parameter 10 is used by Yuan et al. (2023b) to stabilise the numerical solution to the problem shown in Section 3.2.While the present N-PFEM model can capture the smooth distributions of field variables without using any stabilisation techniques.The distributions of pore water pressure calculated from the N-PFEM at t ¼ 50 s and t ¼ 100 s are on a deformed mesh in Fig. 12a  and b, which agree well with the solution from a smoothed PFEM model using a stabilisation parameter of 100 shown in Fig. 12c and d (Yuan et al., 2023b).However, due to the involvement of large deformations, a slightly jagged distribution can still be observed.This can be further improved by increasing the level of discretisation as suggested by Navas et al. (2018).To illustrate the challenges in capturing stabilised pore water pressure, the computed pore water pressure fields at t ¼ 100 s obtained using the smoothed PFEM without stabilisation (i.e. the stabilisation parameter is set at 0) and with a stabilisation parameter of 10 are plotted in Fig. 12e and f, respectively (Yuan et al., 2023b).These results highlight the advantage of the proposed nodal integration-based method, since no stabilisation techniques are required when using linear triangular elements.
The calculated equivalent plastic strains at t ¼ 50 s and t ¼ 100 s are also illustrated in Fig. 13a and b, respectively.The captured shear bands are similar to the ones (Fig. 13c and d) observed by Yuan et al. (2023b) using the smoothed PFEM with a stabilisation parameter of 100.The deformed configurations in Fig. 13 depict the accumulation of plastic strain in narrow zones due to the slip of one part of soil on the other, which has also been reported in the work of Navas et al. (2018) and Sanavia et al. (2002).The present results support the validity of our model for large deformation elastoplastic analysis of typical geotechnical problems.

Conclusions
To develop an efficient and robust PFEM framework for dynamic analysis of saturated media in the context of geotechnical engineering applications, a novel nodal integration-based approach called N-PFEM is proposed.This framework incorporates a nodal integration technique into a generalised HR variational principlederived implicit PFEM formulation.By utilising low-order elements with a strain smoothing technique, the model can significantly reduce the computational cost and eliminate the volumetric locking issue.Additionally, tedious mapping operations after remeshing in the original PFEM are avoided, as all field variables are evaluated at mesh nodes and nodal integrations are performed over cells.
It is demonstrated that the proposed HR variational principlederived formulation resolves the instability issues observed in nodal-based mesh-free methods and other nodal integration-based PFEM formulations for dynamic saturated soils without the need for ad-hoc stabilisation techniques.The model is validated against a series of benchmark problems, exhibiting excellent agreement with reference solutions.Our work provides a valuable reference for the development of advanced numerical tools capable of simulating geotechnical applications involving large deformation and saturated conditions.

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.(Yuan et al., 2023b); and (d) t ¼ 100 s from the smoothed PFEM with a stabilisation parameter of 100 (Yuan et al., 2023b).

Fig. 3 .
Fig. 3. Schematic of the N-PFEM calculation procedures from time step n to nþ1: (a) a cloud of particles at t n ; (b) triangular elements obtained from the alpha-shape method at t n ; (c) constructed node-based smoothing cells at t n ; and (d) updated particle positions at t nþ1 .

Fig. 2 .
Fig. 2. Schematic of node-based smoothing cells created based on adjacent triangular elements.Each smoothing cell consists of N s sub-cells and each sub-cell is one-third of the corresponding linear triangular element.

Fig. 5 .
Fig.5.Wavefields of pore water pressure at t ¼ 0.05 s calculated by: (a) our nodal integration-based SOCP-L and (b) SOCP-Q methods, in comparison to (c) the conventional FE with linear elements without stabilisation(Monforte et al., 2019) and (d) with stabilisation(Monforte et al., 2019).The displacement field is scaled by a factor of 250.

Fig. 7 .
Fig. 7. Numerical results of strip footing on poroelastic medium compared with those of Manoharan and Dasgupta (1995): (a) normalised displacement at points A and B; and (b) normalised pressure at point A.
Fig. 8. Strip footing: (a) normalised displacements at points A and B compared with the results of Manoharan and Dasgupta (1995) and Sabetamal et al. (2016) for elastoplastic case; and (b) normalised pore water pressure at point B compared with the results of Sabetamal et al. (2016) for elastoplastic case.

Fig. 9 .
Fig. 9. Numerical results of pore water pressure obtained in this study and those by two versions of smoothed PFEM (Jin et al., 2021b; Yuan et al., 2023b): (a) elastic analysis at t ¼ 0.01T v from this study; (b) elastic analysis without stabilisation at t ¼ 0.01T v from the study of Yuan et al. (2023b); (c) elastic analysis with stabilisation at t ¼ 0.01T v from the study of Yuan et al. (2023b); (d) elastoplastic analysis at t ¼ 0.001T v from this study; (e) elastoplastic analysis without stabilisation at t ¼ 0.001T v conducted by Jin et al. (2021b); (f) elastoplastic analysis with stabilisation at t ¼ 0.001T v conducted by Jin et al. (2021b); (g) elastoplastic analysis at t ¼ 0.01T v from this study; (h) elastoplastic analysis without stabilisation at t ¼ 0.01T v conducted by Yuan et al. (2023b); and (i) elastoplastic analysis with stabilisation at t ¼ 0.01T v conducted by Yuan et al. (2023b).

Fig. 10 .
Fig. 10.2D self-weight consolidation: (a) model setup; (b) evolution of gravity acceleration; (c) contours of pore water pressure p and the magnitude of displacement u at different instants; and (d) time history of pore water pressure at points A, B and C, compared with those in the work of Yuan et al. (2022).

Fig. 11 .
Fig. 11.Model setup for the vertical cut of a saturated soil block.

Fig. 12 .
Fig. 12. Numerical solutions of pore water pressure: (a) at t ¼ 50 s from this study; (b) at t ¼ 100 s from this study; (c) at t ¼ 50 s from the smoothed PFEM with a stabilisation parameter of 100; (d) at t ¼ 100 s from the smoothed PFEM with a stabilisation parameter of 100; (e) at t ¼ 100 s from the smoothed PFEM without stabilisation (the stabilisation parameter is set at 0); and (f) at t ¼ 100 s from the smoothed PFEM with a stabilisation parameter of 10. (cef) are from the work of Yuan et al. (2023b).

Fig. 13 .
Fig. 13.Numerical solutions of equivalent plastic strain at t ¼ 50 s and t ¼ 100 s: (a) t ¼ 50 s from this study; (b) t ¼ 100 s from this study; (c) t ¼ 50 s from the smoothed PFEM with a stabilisation parameter of 100(Yuan et al., 2023b); and (d) t ¼ 100 s from the smoothed PFEM with a stabilisation parameter of 100(Yuan et al., 2023b).

Table 1
Model parameters adopted in this study.
Parameter Value Density of solid, r s (kg/m 3 ) 2700 Density of fluid, r f (kg/m 3 )