Stabilized explicit u -- p w solution in soil dynamic problems near the undrained-incompressible limit

Traditionally, Biot’s formulation is employed to model the behavior of saturated soils. The u (cid:2) p w (solid displacement– pore water pressure) formulation can be considered as the standard one, since involves a good computational performance together with excellent accuracy for slow and moderate speed phenomena. Dynamic processes can be studied even if the acceleration of the water is neglected, what occurs in the undrained limit. It is well-known that u (cid:2) p w formulation might display instabilities in the undrained-incompressible limit. Several techniques have been proposed to overcome this issue, principally within an implicit time integration scheme for small strains. In this paper, a robust implementation of the divergence of the momentum equation technique is presented for an explicit u (cid:2) p w approach within the framework of optimal transportation meshfree scheme at ﬁnite strain. Several examples are provided in order to assess the good performance of the proposed methodology

Left Cauchy-Green tensor b: Body forces vector c: Cohesion (equivalent to the yield stress, r Y ) C (Time integration scheme): damping matrix Material time derivative of h with respect to the solid F ¼ ox oX : Deformation gradient g: Gravity acceleration vector G: Shear modulus h: Nodal spacing H: Hardening modulus, derivative of the cohesion against time.

I:
Second-order unit tensor J ¼ det F: Jacobian determinant k: Intrinsic permeability k: Permeability tensor K: Bulk modulus K s : Bulk modulus of the solid grains K w : Bulk modulus of the fluid M : Mass matrix n: Porosity NðxÞ, rNðxÞ: Shape function and its derivatives p: solid pressure p w : Pore pressure P (Time integration scheme): external forces vector Q: Volumetric compressibility of the mixture R: Internal forces vector s ¼ r dev : Deviatoric stress tensor t: time u: Displacement vector of the solid U: Displacement vector of the water v s ¼ _ u: Velocity vector of the solid v ws : Relative velocity vector of the water with respect to the solid w: Relative displacement vector of the water with respect to the solid Zðx; kÞ: Denominator of the exponential shape function a F , a Q and b: Drucker-Prager parameters b, c: Time integration schemes parameters b, c: LME parameters related with the shape of the neighborhood Dc: Increment in equivalent plastic strain e p : Equivalent plastic strain e: Small strain tensor e 0 : Reference plastic strain j: Hydraulic conductivity k: Lamé constant k: Minimizer of logZðx; kÞ l w : Viscosity of the water m: Poisson's ratio q: Current mixture density q w : Water density q s : Density of the solid particles r: Cauchy stress tensor r 0 : Effective Cauchy stress tensor s: Kirchhoff stress tensor s 0 : Effective Kirchhoff stress tensor U: Plastic yield surface /: Friction angle w: Dilatancy angle Superscripts and subscripts dev : Superscript for deviatoric part e : Superscript for elastic part k : Subscript for the previous step kþ1 : Subscript for the current step p : Superscript for plastic part s : Superscript for the solid part trial : Superscript for trial state in the plastic calculation vol : Superscript for volumetric part w : C superscript for the fluid part relative to the solid one

Introduction
Porous media dynamics at large strain is a cutting-edge issue in the field of computational geomechanics. In this respect, the u À w À p w and u À w Biot formulations, where no dynamic terms are suppressed, are of paramount significant [27]. However, these formulations are computationally expensive compared with the u À p w formulation, since the latter requires, for each node, fewer degrees of freedom than the former ones. Although there are several problems where it is necessary to employ a complete formulation, since every dynamic term becomes important, the majority of the dynamic problems in saturated soils can be covered by the u À p w formulation without loss of accuracy, as most of them fall within its range of application [64]. Nevertheless, these formulations are not exempted from numerical drawbacks depending on the spatial discretization of the numerical method, even for the meshfree approach herein presented. An unstable behavior in numerical solutions of saturated porous media is usually manifested by an over-stiffening of the mechanical solution and a great and fictitious space variability of the water pressure field, occasionally leading to non-uniqueness and mesh dependence of the solutions as well as some instabilities. The occurrence of this kind of pathologies is connected with the undrained-incompressible limit of the coupled problem. While this type of pathology in the u À w À p w and the u À p w formulations is connected with the undrained-incompressible limit, that is, analogous to the incompressible limit in solid mechanics problems [46], in the u À w, the pathology appears when the compressibility of the fluid is some orders of magnitude higher than the compressibility of the mixture [40]. In order to overcome this numerical issue in the hydromechanical response of a saturated porous media, several alternatives have been proposed in the specialized literature.
The first approach aims to stabilize the pore pressure term by the inclusion of new terms in the formulation. Monforte et al. [34] assessed two of these techniques: Divergence of the momentum equation (DME) and polynomial projection technique (PPP). The DME was proposed for water problems by Hafez and Soliman [18] and later improved for the u À p w formulation by Pastor [45]. On the other hand, the PPP [14] was applied to the u À p w first by White and Borja [59] and extended by Sun et al. [57] and by Gavagnin et al. [17], being the first one for small strain and the latter for finite strain.
A different approach relies on an average pressure projection. The intrinsic idea is inspired by the work of Hughes [22] and the recent developments of the B-bar method [16]. Furthermore, the strategy behind is analogous to the diamond elements one (Hauret,Kuhl and Ortiz [19]). This approach has been extended to large strain formulations by Simo et al. [54,55] and De Souza Neto and coworkers [11,12], considering the F-Bar technique as a volume averaging procedure, changing the volumetric part of the strain tensor by a non-local one. More recently, Sun and coworkers [56,57] have proposed a novel alternative, starting from the work of Moran et al. [35]. Related techniques have been extended to multiphase formulations for small and large strains in the context of the Optimal Transportation Meshfree method [36,37,[40][41][42][43]. Also in the finite volume field exist several approaches that stabilize the poromechanical instabilities through macroelement techniques (see the work of Camargo et al. [8,9]).
Furthermore, stable formulations can be also achieved by tuning the time discretization, as the split-operator [23,68]. Some of these implementations employ fractional steps techniques that provides a stabilization which allows the employment of the same order of interpolation in mixed formulations. This was first noticed by Schneider et al. [52] and Kawahara and Ohmiya [25]. Zienkiewicz and coworkers [66] explained the reasons why the fractional step technique provides stabilization for computational fluid dynamic problems. This technique was also applied by Pastor et al. [46] for coupled analysis in saturated soil problems, and extended later by Li et al. [31], being worth mentioning the works of White and Borja [59] and Mira et al. [33]. In the geotechnical field, fractional step techniques have been proposed within several meshfree techniques such as the SPH formulations proposed by Blanc and Pastor [4] and the MPM ones proposed by Kularathna et al. [26].
Finally, different order of interpolation for pore pressure and displacement has been the main technique to obtain smooth results, although it requires a stronger computational effort as requires larger number of nodes for the solid phase than for the pore water pressure. Mira et al. [33] in the u À p w formulation and Jeremić et al [24] within u À U À p w formulation employed different order of interpolation for the different variables.
A summary of the different stabilization techniques in the context of small strain regime can be found in Markert et al. [32].
In this research, the explicit u À p w is employed due to its computational benefits, which come from both the explicit Newmark predictor-corrector scheme, which was first developed by Navas et al. [40] for the explicit u À p w formulation, and the intrinsic computational saving benefits of the employment of the u À p w against the complete formulations. Regarding the application of this formulation in a finite deformation approach, the first works were tested simulating the constitutive behavior of the solid phases with linear elastic (Diebels and Ehlers [13]), Cam-Clay (Borja et al. [5,6]) and Drucker-Prager theories (Armero [1]) as well as with anisotropic materials [62]. At the same time, Ehlers and Eipper [15] applied a new elastic nonlinear constitutive model to reproduce the compaction of the soil. Implicit schemes were employed in these works, being necessary the linearization of the derivatives of the u À p w equations. Furthermore, Sanavia et al. [49] considered several neglected terms in the linearization of the previous works and extended the methodology to unsaturated soils [51]. Moreover, with meshfree schemes, we find excellent contributions to the explicit u À p w approach (see [60,63]) in the small strain approach. The recent work of Navas et al. [40] shows excellent results within an explicit scheme at large strain.
In the present research, a novel robust stabilized predictor-corrector explicit algorithm for the u À p w formulation at large strain is proposed. The numerical framework considered in order to implement this new stabilized formulation is the optimal transportation meshfree (OTM) scheme [2,28]. The door is opened to any extension to traditional or novel computational techniques, which can be easily made by adapting the spatial discretization.
The rest of the paper is structured as follows. The problem formulation is summarized in Sect. 2. The explicit methodology, as well as the constitutive models of the solid response, is drawn in Sect. 3. The stabilization technique is detailed in Sect. 4. Applications to several problems are illustrated in Sect. 5. Relevant conclusions are depicted in Sect. 6. The definitions of all symbols used in the equations are provided in the nomenclature appendix.

Biot's equations: u-p w formulation
The u À p w formulation is one of the most common forms of the Biot's equations [3]. This theory is based on three statements: the mechanical response of a solid-fluid mixture, the continuity of flux through a differential domain of saturated porous media and the coupling between the aforementioned phases. The derivation of the formulation is widely encountered in the literature. All the equations presented hereinafter are derived from those presented by Lewis and Schrefler [27], later proposed for the finite strain setting by Sanavia and Navas [42,50,51].
Respecting the notation, bold symbols has been employed herein for vectors and matrices as well as letters of the Latin alphabet for scalar variables. Let u represents the displacement vector of the solid skeleton, w the displacement vector of the fluid phase with respect to the solid and p w the pore water pressure, which is assumed positive for compression. Terzaghi's effective stress [58] is considered being defined for incompressible solid constituents as: where r 0 is the effective Cauchy stress tensor, r is the total Cauchy stress tensor (both positive in tension) and I is the second order unit tensor. Furthermore, D s =Dt is the material time derivative following the movement of the solid particles, also defined as _ h. The main assumption with respect to the full Biot's formulation lies on neglecting the accelerations of the fluid phase. Thus, the linear balance equation of the mixture reads: The continuity equation, also known as mass balance equation, is employed as well: Moreover, both linear momentum balance equation of the fluid phase, defined through the Darcy's theory, and the equilibrium of the mass, can be written together in order to reduce the number of degrees of freedom of our problem: In equations (2) and (4), q is the density of the mixture, while q w and q s are the fluid and solid particle densities, respectively. These three densities are related to each other through the expression: where n is the porosity defined by the ratio of volume of voids, V v , with respect to the total volume, where V s is the volume of the solid grains: The soil is considered totally saturated. Thus, V v is equal to the water volume. Furthermore, the volumetric compressibility of the mixture, Q [64], has been calculated considering incompressible solid constituents as where K w is the bulk modulus of the fluid phase, usually considered as water.
On the other hand, g represents the gravity acceleration vector and , the permeability matrix. This matrix can be written in terms of k rw , the hydraulic conductivity, j [m/s], and the specific weight of the water through the following equation: k rw is set as 1 when full saturation is to be considered.

Implementation tools
In this section, spatial and time discretization, as well as the constitutive models employed in this research, are briefly explained.

Spatial discretization
The optimal transportation meshfree [21,28,29], as one of the most trending meshfree methods, is employed in the present research. Its fulfillment within geotechnical multiphase problems [38] makes this methodology a robust alternative to material point method or Smoothed Particle Hydrodynamics, among others. Similarly to the first one, OTM is based on a particle-to-node interpolation, where the shape functions that relates both sets are based on the work of Arroyo and Ortiz [2], who defined the Local Max-Ent shape function (LME) of the material point with respect to the nodal neighborhood. For the readers' knowledge, the equation of the shape function and the achievement of the first derivatives of the shape function, as well as the details of every parameter, are depicted in [2,28].
Applying the standard Galerkin procedure to the weak counterpart of Eqs. (2) and (4) (See [48,51] for details), the following equations are obtained: where: while the mass and damping matrices, are written as follows: In these expressions V p is the volume associated with the material point P while N p range from the neighborhood points associated with P. The matrix B is the symmetric shape function gradient operator while m the identity matrix in Voigt notation. Thus, Bm reproduces the divergence operation.

Explicit integration
The coupled problem in the time domain is solved in this research through the explicit central difference Newmark time integration scheme. Let the current step k þ 1, and assuming the solution in the previous step, k, has been already obtained, the relationship between displacement, velocity and acceleration, u kþ1 , _ u kþ1 and € u kþ1 , follows the Newmark scheme for each phase: By considering h ¼ c ¼ 0:5 and b ¼ 0, the Predictor and Corrector terms are read as follows: where the underlined terms are the predictor ones, denoted as _ u kþÃ and p w kþÃ , respectively. On the other hand, the stability of this method is assured through the Courant-Friedrichs-Lewy (CFL) condition, being the time step, Dt, smaller than h V c , where h represents the discretization size and V c is the velocity of the p-wave (see [65]), defined by: It bears to emphasize several aspects about the predicted step. About the stress, it has to be calculated in this predicted step through the predicted displacement as follows: The logarithmic strain, as the strain measure to be employed in the deformed configuration, has to be calculated from the tensor b, the Left Cauchy-Green strain tensor (b ¼ FF T ), which depends on the displacement on the predicted step as well: Moreover, solid velocities and water internal forces must be evaluated in the predicted step, k þ Ã. In Section 4, all these premises will be considered in a pseudo-algorithm, written once the stabilization is taken into account.

Constitutive models for the solid phase
A brief note of the employed constitutive models is outlined in this Section: a hyper-elastic model for soils and another one that involves plastic deformation which follows the Drucker-Prager failure criterion, both developed in a finite strain approach. As mentioned before, for elastic materials, a nonlinear law, considering the influence of the Jacobian calculation, the initial porosity n 0 and the compaction point of the soil, is employed (see Ehlers and Eipper [15]). This law follows the subsequent definition of the Kirkchoff stress tensor: In order to reproduce elasto-plastic behavior at large strain, the Drucker-Prager yield criterion proposed by Sanavia and coworkers for large strain is employed [48,51]. The relationship between left Cauchy-Green strain tensor b, calculated at the current configuration, and the small strain tensor e is made through Ortiz, Simo and coworkers' research [10,44,53]: e e trial kþ1 ¼ Trial stress measures are also computed from the elastic trial strain as: where G and K represent the shear and bulk moduli of the solid, respectively. The yield conditions for the classical and apex regions are calculated in a different manner for the proposed yield criterion. Once known which algorithm to employ, one of the following expressions applies: where Dc 1 ¼ ks trial kþ1 k 2G , Dc and Dc 2 are the objective functions in the Newton-Raphson scheme for the classical or apex regions accordingly. In Eqs. 20) and (21), H is the hardening parameter of the material, c k represents the cohesion, and a F ,a Q and b depend on the friction and dilatancy angles as well as the shape of the yield surface. For details of the calculation of the equivalent plastic strain, see [48,51].

Stabilization of the water phase
The necessity of stabilization in the developed u À p w formulation lies on the work of Zienkiewicz and Taylor [67]. Assuming incompressible constituents, the Biot modulus Q is larger than the compressibility of the solid skeleton. Traditionally, the usual condition required for the u À p w problem is given by n u ! n p w (see, for instance, Pastor et al [45]). However, this condition, although necessary, must be fulfilled in any assembly of elements of the mesh [45,67]. See [20] for a detailed study of the stability conditions of the three-field formulation.
The pore pressure term in the algebraic equation has been traditionally the one which suffers the aforementioned oscillations. Thus, the majority of the available techniques have this term as the target one, departing from the stabilization techniques employed in the context of Computational Fluid Dynamics. Brezzi and Pitkaranta [7] proposed to add to the pressure term, the so-called K pp position in the implicit mixed solid u À p formulation, the stabilizing term: which is equivalent to the one proposed by Hughes et al [23]. This term, when working with the u À p w implicit saturated formulation, is analogous to the K p w p w term, i.e., the pore pressure position of the continuity equation, Eq. (4) (see Pastor et. al. [45] for details). Based on Petrov-Galerkin theory; the added term is proportional to the product of the linear momentum equation by h 2 2 l rN. Indeed, both techniques are proportional to h, the mesh size in Finite Element approaches or the nodal distance in meshfree methodologies. This fact leads to an interesting issue when the mesh is refined: the proposed term tends to zero and, thus, the consistency is preserved. The methodology to be explained hereinafter is inheritor of the work of Hafez and Soliman [18] proposed to stabilize the Navier-Stokes equations and the one proposed by Pastor et. al. [45], which adapted the former one to the hydromechanical problem in quasi-static conditions (u À p w formulation). That work was later extended to a formulation valid for the full Biot equations by Monforte et. al. [34]. All these works are based on adding the divergence of the momentum equation to the mass balance equation.
In our case, in the explicit scheme, Eq. (10) may suffer numerical instabilities in the undrained limit. If , tends to zero, M w , f ext; w , and R wÃ tend to be much smaller than C _ u. It leads to the situation that the pore pressure velocity depends only on the solid velocity, what means that no influence of the pore water phase affects the results. Due to this fact, spurious pore pressure oscillations are obtained. Moreover, in the incompressible limit, the increase in the compressibility of the mixture, Q, aggravates this situation since tends to increase the obtained result.
Thus, the objective is to add a new term, which is equal to zero in order to avoid nonphysical results, but avoids leaving the pore pressure value just in the hands of the velocity of the solid phase. Following the aforementioned works ( [34,45]), in our proposal, the divergence of the linear momentum balance equation of the mixture, Eq. (2), will be employed. This equation, multiplied by the parameter a= a(h), yields: taking into account the divergence of the volumetric forces is equal to zero. If the solid skeleton behavior is elastic, which can be perfectly assumed even if the material plasticizes, the divergence of the solid displacement can be expressed in terms of the volumetric stress invariant, p 0 , and the bulk modulus, K, as: From the same assumption, the divergence of the internal forces can be expressed as (see Pastor et. al. [45] for details): where k and G are the Lamé constant and shear modulus of the soil skeleton.
On the other hand, the displacement of the water is negligible if small permeability is employed, which is the case of the undrained limit, the goal of the proposed stabilization technique. Thus, from the integration on time of the equation (3) and plugging definition (24), we can obtain: Thus, considering Eq. (25), the divergence of the internal forces of both fluid and solid phases yields: Taking into account Eqs. (23) and (27), the divergence of momentum of the mixture can be rewritten as: We have considered the incompressible limit, where k and G are much smaller than the bulk modulus of the water phase, Q.
Pastor and coworkers proposed that the parameter a may depend on the density of the mixture and the critical time step, which is the time step that is able to capture the P wave of the saturated media. The main components of this reference time step are the characteristic P-wave velocity, V c , and the element size, h: Finally, as Eq. (28) is equal to zero, we can add directly to Eq. (4) without losing physical meaning: Rewriting Eq. (30) in a matrix form, the stabilized system of equations yields: where the Ã terms are the stabilized terms that are built as:

Explicit stabilized pseudo-algorithm
The explicit scheme is similar to the one proposed by Navas et al. [39] with the introduction of the stabilization terms. Following, the pseudo-algorithm of the whole model is proposed, being the superscript p employed for material point calculations.

Explicit Newmark Predictor
2. Update position of Nodes and Material points Du a kþ1 N a ðx p k Þ:

Calculation of the deformation gradient and related parameters
Du a kþ1 rN a ðx p k Þ; 4. Update density and mass. Construct lumped mass q kþ1 ¼ n kþ1 q w þ ð1 À n kþ1 Þq s :

Verification examples
Three different verification examples are considered in this section. The first verification case deals with a consolidation process in a saturated column of soil subjected to a bilinear loading at the surface The second one allows us to verify also the 2D component behavior, analyzing the behavior of the pore water pressure distribution along a 2D stratum when it is loaded by a strip footing. Finally, the third verification case assesses the failure of a vertical wall of saturated soil, aiming to analyze the performance of the proposed stabilization technique in a traditional geotechnical problem.

Consolidation of a column of soil
In this problem, an idealization of a semi-infinite stratum of soil through a 2D column is employed. The height and width of the column are H T ¼ 10m and L ¼ 1m, respectively. Horizontal displacement is forbidden in the vertical boundaries, while vertical displacement is blocked in the rigid base. In addition, the top boundary is considered perfectly drained (p w ¼ 0). This geometry and boundary conditions are depicted in Fig. 1. Also, the load ramp considered is presented. Sabetamal et al. [47] proposed a particular mesh in order to be able of a proper capturing of the wave provoked by the load. The upper meter of the stratum is discretized with a 0.25 m.; meanwhile, a discretization of 0.5 m. size is employed for the rest of the stratum. A similar discretization is employed in this research within the OTM configuration.
The objective of the present section is the verification of the methodology in the undrained-incompressible limit. A stratum of soil suffers the application of a pressure which follows the curve of Fig. 1B. The load is applied gradually until reaching P max at t 0 ¼ 7:5 s. Following, the l when the load is kept constant until 7.5 s. The water and solid parameters are presented in Tab. 1, being the nonlinear elastic material of Eq. (15) employed for the solid model of the problem.
The solutions obtained, with and without the stabilization technique, are shown in Fig. 2 for the two cases. In this figure, the pore pressure distribution, when the whole load is applied (7.5 s), is plotted. The typical positive-negative pattern of an unstable pore pressure distribution is observed for the unstabilized case. Also, spurious maximum values, close to 500 Pa, are observed in this case; meanwhile, the maximum value of the stabilized case is close to the total amount of load, 100 Pa. Also, in Fig. 2, the pore pressure distribution obtained with a quadratic finite element scheme can be observed. A similar pattern can be appreciated between the proposed approach and the quadratic finite element one. However, oscillations close to the top of the column appear in the distribution obtained with the finite element approach. These oscillations are not present in the results obtained with the proposed approach.
If this pore pressure is plotted along the depth, Fig. 3, the oscillations are seen for the unstabilized case in the 3 first meters of the stratum of the column. The stabilization technique allows to mitigate this unrealistic behavior, smoothing out the solution without losing physical information. Although the pore pressure distribution obtained with the quadratic finite element approach improves the unstabilized result, it still presents some oscillations close to the top of the column.

Strip footing load
The following example is based on a saturated stratum loaded by a vibrating strip footing. Fig. 4 presents the geometry, boundary conditions as well as the history load of the aforementioned footing. This verification case was firstly studied by Li et al. [30] through the finite element method and later with the Material Point Method by Zhao and Choo [61]. The permeability of the soil stratum in the present work is similar to the one considered by Zhao an (A) (B) Fig. 1 Geometry, load and boundary conditions of the column of soil problem Choo, which is lower than the original one. The traditional Neo-Hookean model is employed. The material parameters are depicted in Fig. 4. It is worth mention that comparison against these works is not possible since no damping parameter is employed in the present research. However, it is possible to capture similar trends that were observed in the previous works. Moreover, they considered an incompressible pore water phase, whereas compressibility of the pore water is taken as 1e10 Pa in the present work. First, the spatial distribution of the pore pressure shows oscillations, since it lies on the incompressible-undrained limit (Fig. 5). The values are not so severer than in the previous example, since no negative values are achieved. It can be due to the 2D configuration of the problem or the dynamic behavior of the load. However, unrealistic peak values are obtained.
The phenomenon is studied at t=0.1 s. In Fig. 5, the distribution of the pore pressure is depicted for the unstabilized and stabilized solutions and also the quadratic finite element scheme. The unstabilized profile shows a spurious oscillation, which is alleviated in the stabilized and the finite element solutions. The difference between the stabilized and the finite element solutions lies in the peak, which is more pronounced for the finite element one. This peak is quickly alleviated along the depth in both cases, being quicker in the finite element approach than in the stabilized solution. The value of the peak can be observed in Fig. 6. In this figure, the distribution of the pore pressure along the depth is depicted, located under the edge of the strip footing. While the stabilized solution and the quadratic one show a smooth transition, with a peak of around 3 MPa, the unstabilized solution oscillates around the correct solution providing unreal values, with a peak of 4 MPa.
Next, the time evolution of the pore water pressure at two different locations (A and B in Fig. 4, i.e., a shallow and a deep location) is analyzed.
As it can be appreciated in Fig. 7, the main difference can be encountered in the location at point A, shallow zone. Effectively, in Fig. 7A, there is an important overestimation of the peak values that is mitigated along the depth, as it can be observed in Fig. 7B. Although the overall trend is not greatly affected with depth, several peaks appear. In the   present case, as the material is Neo-Hookean, these peaks scarcely affect the final result. However, with different constitutive materials, as in the following application case, the found peaks would alter the final result.
The differences between the quadratic and the stabilized solutions are depicted clearer in point B, where the quadratic one is able to dissipate pore pressure more efficiently than the proposed solution.

Vertical cut
In this final application, the proposed stabilization technique is applied to model a drainage process induced by a strip rigid footing over a square domain of a hyperelastoplastic saturated soil. The failure of the material, leaded by the applied load, takes the shape of a typical vertical cut with an inclined shear band. A similar problem was previously assessed by Sanavia et al. [49,50] in a quasi-static regime and by Navas et al. [41,42] for a dynamic approach, considering different dilatancy angles in order to see the performance of the soil depending on the properties. In this case, only the case of friction angle 20 and dilatancy 0 is studied. The geometry and material properties are shown in Fig. 8A gradual displacement of 1 m at the boundary C 4 is applied along the 10 s of the simulation, being the employed time step 1 ms. A nodal spacing of 0.833 m., forming a regular 12Â12 nodal discretization, is considered.
Pore pressure results are depicted, at the final stage, in Fig. 9. For the stabilized simulation, it was reached the final time of 10 s. However, for the non-stabilized case, the simulation breaks out at 1 s due to the influence of the spurious pore water pressure distribution over the solution process. In the third row, the solution obtained with a T6-P3 finite element technology, without stabilizing, is also shown. The results are similar to the stabilized one. In the referred bibliography, we found distributions of pore pressure similar to the stabilized one, taking into account that, in this case, the permeability was lowered to reach undrained conditions. Despite this fact, the trend of the behavior of the soil is well captured. Contrary, the unstabilized one presents the typical oscillations that become bigger due to the plastic behavior of the soil.   Fig. 10. It can be observed how the stabilized one and the quadratic one are similar to those reached by Sanavia and coworkers for the case / ¼ 20 and w ¼ 0 . Nevertheless, as we mentioned before, the spurious pore water pressure spatial distribution affects the spatial distribution of plastic zones and the subsequent development of the shear band mechanism.

Conclusions
The stabilization of the explicit Newmark predictor-corrector solution of the Biot's u À p w formulation is provided in this manuscript. This formulation has been at large strain within the optimal transportation meshfree method, although it can be extended to some other well-known meshfree schemes, such as smooth particle hydrodynamics or material point method. The stabilization is achieved through the addition of a stabilization term, which comes from the divergence of the momentum equation.
The proposed methodology was previously demonstrated to work in a robust way in dynamic problems of saturated soils, mainly under drained conditions, but in this research, it has been extended to the undrained-incompressible limit. Furthermore, the methodology accounts with the potential strength of an explicit scheme. Elastic and plastic materials have been tested, showing excellent results, mainly with the second one, since the spurious peak values are catastrophic by means of unreal plasticizing of the material.
The first example has demonstrated the performance of the proposed methodology in the traditional Terzaghi's scheme. The second one extends this elastic scheme to 2D conditions in order to visualize the propagation of the waves along the domain. Harmonic loading has been employed to assess the behavior of the pore pressure along the time. Finally, a typical plastic problem has been analyzed. Under the light of the proposed stabilization technique, the expected shear band is attained. However, without stabilizing the problem, the plastic zone dissipates creating spurious plastic zones along the domain, being impossible the formation of the shear band. The wide range of applications allows us to validate the proposed methodology under dynamic conditions. The present research opens the door of the explicit simulations of saturated soils by some other meshfree methods. Next steps will present this formulation in the material point method. Moreover, this methodology is feasibly transferable to unsaturated conditions by adding the influence of the degree of saturation and the air phase. Because of the nature of the explicit scheme, this formulation arises easily from the proposed in this paper. Authors would like to thank the administrative and technical support of the ETSI Caminos, Canales y Puertos, from the Universidad Politécnica de Madrid, as well. Additionally, the fifth author really appreciates the Entrecanales Ibarra Foundation for his undergraduate scholarship.
Author contributions Conceptualization and mathematical methodology were carried out by PN and MMS Implementation was done by MM and PN Validation was done by AY and DM Meanwhile, supervision, project administration and funding acquisition belong to MP. All authors have contributed to the writing and original draft preparation, and they read and agreed to the published version of the manuscript.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.

Declarations
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.