Abstract
In this paper, an efficient and robust methodology to simulate saturated soils subjected to lowmedium frequency dynamic loadings under large deformation regime is presented. The coupling between solid and fluid phases is solved through the dynamic reduced formulation \(up_\mathrm{w}\) (solid displacement – pore water pressure) of the Biot’s equations. The additional novelty lies in the employment of an explicit twosteps Newmark predictorcorrector time integration scheme that enables accurate solutions of related geomechanical problems at large strain without the usually high computational cost associated with the implicit counterparts. Shape functions based on the elegant Local Maximum Entropy approach, through the Optimal Transportation Meshfree framework, are considered to solve numerically different dynamic problems in fluid saturated porous media.
Introduction
Modeling saturated soils under dynamic loads is of crucial importance when researchers deal with fast phenomena, like landslide propagation. Moreover, the more catastrophic 3D problems should also include large deformations. Despite the importance of both aspects, dynamic saturated problems and large deformations, research focused on these topics at once is limited. One of the main reasons why there is a lack of research in this paramount issue of geomechanics might be that a powerful and useful methodology requires complex hydromechanical models including inertial terms coupled with hyperelastoplastic constitutive models where the deformation gradient acts as strain measure. Moreover, as analytical solutions can only be achieved for few idealized configurations, cutting edge numerical techniques must be considered to attain accurate and robust solutions in realworld problems.
The fluid saturated phenomenon has been widely studied in the numerical geotechnical field, where a big range of solutions can be found regarding the formulation considered for the coupled problem (either simplified or complete), the assumptions made with or without accelerations and the numerical techniques used to solve the equations, both in the spatial (mesh or meshfreebased techniques) and temporal dimension (explicit or implicit schemes).
The first formulations aimed to describe the physics behind a saturated porous medium are found in the governing equations introduced by Biot [4], later reviewed by Zienkiewicz and coworkers [53,54,55,56]. Similar equations were obtained by Lewis and Schrefler [20] within the Hybrid Mixture Theory, in this case, starting from the microscopic scale, improving the consistency and robustness of the formulation. Regarding the inertial terms, both acceleration of fluid and solid phases are employed in the complete formulation, covering a wide range of frequencies [18, 39]. This formulation is usually expressed in terms of the relative water displacements, w, which has been proved to be successful [24, 29]. However, other research presents this complete formulation by means of the total displacement of the water, U, as a nodal unknown (Ye et al. [49] and Sladek et al. [46]). Concerning the simplified formulations, the \(uw\) approach is computationally more expensive than the \(up_\mathrm{w}\) since the former employs more degrees of freedom per node. Thus, its utilization is not recommended when the \(up_\mathrm{w}\) formulation is sufficient to capture the complete wave propagation in a saturated soil problem. The \(up_\mathrm{w}\) (solid displacement–fluid pressure) formulation is widely used in dynamics to solve different hydromechanical coupled problems due to its simplicity as well as the high accuracy achieved for a great variety of geomechanical problems (e.g., [10, 53, 54]) .
All these formulations have been usually solved in time through an implicit scheme [1, 6, 8, 15, 34]. Recently, Navas et al. proposed an explicit solution of the \(uw\) formulation with excellent results, see [33]. Explicit schemes are a feasible alternative in which there is no need to compute the tangent stiffness matrix, avoiding the complex linearization of the governing equations. Moreover, the computational effort is minimized as forward values are computed directly from the current one, avoiding the solution of nonlinear equations system when advancing in time. Finally, explicit schemes allow a more efficient use of multicore processor, thus facilitating a parallel programming paradigm.
Regarding the application of the Biot’s equations under large deformation regime, the first works were carried out by Diebels and Ehlers [14], Borja et al. [6, 8] and Armero [1] who tested their models by simulating the constitutive behavior of the the solid phases with an hyperelastic, CamClay and Drucker–Prager theories, respectively. Around the same period, Ehlers and Eipper [15] applied a new NeoHookean constitutive model to represent the compaction of the soil up to the solid compaction point. An interesting extension was made by Sanavia et al. to unsaturated soils [42,43,44]. Again, most of these models were solved employing implicit schemes where the linearization of the \(up_\mathrm{w}\) equations was necessary. There is a scarcity of examples in the specialized literature of explicit solutions in time for the \(up_\mathrm{w}\) formulation under large deformation regime. However, for saturated porous media undergoing a fast deformation process, this type of time integration schemes is a feasible alternative, as the usual restriction required for the time step to attain a stable solution can be assumed, as far as numerical efficiency is concerned. The present research aims to cover the lack of explicit time integration solutions for the \(up_\mathrm{w}\) formulation undergoing large deformations.
In recent years, in the computational mechanics field, large strain approaches go hand in hand with meshfree methods due to their numerous advantage to reproduce large relative displacements. In the geotechnical field, this combination of tools shows excellent results in problems such as landslides, liquefaction or other natural disasters. Saturated soils are also modeled through these approaches. Recent promising works can be found in the literature like Pastor et al. [38] with the Smooth Particle Hydrodynamics (SPH) and the works of Bandara and Soga [3], Ceccato and Simonini [11] or Zhao and Choo [51] with the Material Point Method (MPM). Precisely, with this meshfree scheme, we find excellent contributions to the explicit \(up_\mathrm{w}\) approach (see [50, 52]). The small strain approach is employed within this research.
The goal of the present research is the proposal of a robust predictorcorrector explicit algorithm for the \(up_\mathrm{w}\) formulation at large strain where the spatial domain has been discretized into nodes and material points following the Optimal Transportation Meshfree (OTM) scheme of Li et al. [21]. The shape functions developed by Arroyo and Ortiz [2] based on the principle of maximum entropy [37] are also employed.
The rest of the paper is organized as follows. The Biot’s equations are presented in Sect. 2, with emphasis within the \(up_\mathrm{w}\) formulation. The constitutive models employed to model the solid behavior are summarized in Sect. 3. The discretization techniques, highlighting the meshfree solution and the explicit methodology, are provided in Sect. 4. Applications to various problems are illustrated in Sect. 5. Relevant conclusions and future lines are drawn in Sect. 6. The definitions of all symbols used in the equations are provided in the nomenclature appendix.
Biot’s equations: \(up_\mathrm{w}\) formulation
The Biot’s equations [4] are based on formulating the mechanical behavior of a solidfluid mixture, the coupling between different phases, and the continuity of flux through a differential domain of saturated porous media. Hereinafter, the balance equations will be derived from Lewis and Schrefler [20] in the spatial setting (see [20] or [43, 44] for the kinematic equations), departing from the more general equation, and, in order to reach the compact \(up_\mathrm{w}\) form, making the necessary hypotheses.
Concerning the notation, bold symbols are employed herein for vectors and matrices as well as regular letters for scalar variables. Let \(\varvec{u}\) and \(\varvec{U}\) represent the displacement vector of the solid skeleton and the absolute displacement of the fluid phase, respectively. Since in porous media theory is common to describe the fluid motion with respect to the solid, the relative displacement of the fluid phase with respect to the solid one, \(\varvec{w}\), is introduced and expressed as [25]
where \(S_\mathrm{w}\) is the degree of water saturation and n the soil porosity. Note that \(\varvec{ \left( Uu\right) }\) is usually termed as \(\varvec{u}^{ws}\) in the literature [20].
Let \(\rho \), \(\rho _{w}\) and \(\rho _{s}\), respectively, represent the mixture, fluid phase, and solid particle densities, the mixture density can be defined as function of the porosity:
In the above equations, the porosity, n, is the ratio between the voids volume, \(V_\mathrm{v}\), and the total volume, \(V_\mathrm{T}\):
where \(V_s\) is the volume of the solid grains.
Next, we first explain in detail the derivation of mass balance and linear momentum equations for a fluid saturated multiphase media. Then, the final \(up_\mathrm{w}\) formulation is presented. The following equations are first given by Lewis and Schrefler [20]. In this research, \(D^s/Dt\) denotes the material time derivative with respect to the solid, considering:
where \(\varvec{a}^s\) and \(\varvec{a}^{ws}\) are the solid acceleration and the relative water acceleration with respect to the solid, respectively, being the proposed expressions based on the relationships \(\varvec{\dot{u}} \equiv \varvec{v}^{s}\) and \(\varvec{\dot{w}} \approx nS_\mathrm{w} \varvec{v}^{ws}\).
Derivation of the mass balance equation
The general mass balance equation in a multiphase media for compressible grains is presented next. Let \(p_\mathrm{w}\), \(p_\mathrm{g}\) represent the water and gas pressures, respectively, T, the temperature, then this general mass balance equation is written as follows,
where the right hand side term represents the quantity of water lost through evaporation for unit time and volume. The thermal expansion coefficient of the solidfluid mixture, \(\beta _\mathrm{sw}\), is a combination of that of the solid, \(\beta _\mathrm{s}\), and the fluid, \(\beta _\mathrm{w}\):
In addition, \(\alpha \) is the Biot’s coefficient:
where K denotes the bulk modulus of the solid skeleton. Biot’s coefficient may be usually assumed equal to one in soils as the grains are less deformable than the mixture.
In the current work, the soil is assumed to be totally saturated, i.e., \(V_\mathrm{v}\) coincides with the water volume, which results \(S_\mathrm{w}\) equals to one and \(S_g=0\). As we also consider isothermal multiphase media, \(D^s T/D t =0\), \(e^w=0\), consequently, \(D^s S_w/D t =0\). Taking into account all these assumptions, the volumetric compressibility of the mixture, Q [54], can be calculated as
where \(K_s\) is the bulk modulus of the solid grains, whereas \(K_w\) is the compressive modulus of the fluid phase (usually water).
If additionally the spatial variation of the fluid density is neglected and we take into consideration Eq. (7), (4) is simplified as,
Linear momentum balance equations
On the one hand, the relative velocity of the fluid, \(\dot{\varvec{w}}\), in Eq. (4) is defined through the generalized Darcy law as [20]
where \(\varvec{g}\) represents the gravity acceleration vector, \(\varvec{k}\), the intrinsic permeability tensor of the porous matrix in water saturated condition, considered isotropic in this research (\(\varvec{k}=k \varvec{I}\)), \(k^{rw}\) is the water relative permeability parameter (a dimensionless parameter varying from zero to one) and \(\mu _\mathrm{w}\) is the dynamic viscosity of the water [Pa \(\cdot \) s]. The intrinsic permeability k, expressed in [m\(^{2}\)], is related with the notion of hydraulic conductivity, \(\kappa \) [m/s], by the following equation
On the other hand, according to Lewis and Schrefler [20], the linear momentum balance equation for the multiphase system can also be expressed as the summation of the dynamic equations for the individual constituents relative to the solid as, i.e.,
where the convective terms, related to the acceleration terms, have been neglected, which is normal in soils.
In the calculation of the internal forces of the soil, the Terzaghi’s effective stress theory [47] will be followed, which is defined as follows:
where \( \varvec{ \sigma '} \) and \(\varvec{ \sigma }\) are the respective effective and total Cauchy stress tensors (positive in tension), whereas \(\mathbf{I} \) is the secondorder unit tensor. Contrary, pore pressure \(p_\mathrm{w}\) is assumed positive for compression.
Plugging Eq. (12) into Eq. (11), the linear momentum equation can be written as follows
The \(up_\mathrm{w}\) formulation
Considering the three Biot’s equations, the \(\varvec{u}p_\mathrm{w}\) assumes that accelerations of the fluid phase are negligible. Thus, Eq. (13) yields:
Moreover, in order to avoid the employment of \(\varvec{w}\) as a degree of freedom of our problem, Eqs. (8) and (9) can be combined and the mass equation can be expressed as
Constitutive models for the solid phase
In this Section, the hyperelastic and hyperelastoplastic models, employed within this research, are outlined. Further information of both constitutive laws can be found in [30, 33, 34].
NeoHookean material model extended to compressible range
In this research, the NeoHookean constitutive behavior has been considered as a extension of the elastic one in the large strain regime. Moreover, among several variants, the one proposed by Ehlers and Eipper [15] has been chosen. This law takes into consideration the compaction point of the soil, from the influence of the initial porosity \(n_0\) and the Jacobian, calculated as the determinant of the deformation gradient \(\varvec{F}\), in the following manner:
where \(\varvec{\tau }'\) and \(\varvec{b}\) are the effective Kirchhoff stress tensor and the left Cauchy–Green tensor, respectively, whereas G and \(\lambda \) are the Lamé constants.
Drucker–Prager yield criterion
In order to reproduce frictionalcohesive behavior at large strain, the traditional Drucker–Prager yield criterion [41, 44] has been extended to large strain procedure. This methodology follows the work of Ortiz, Simo and coworkers [12, 36, 45] to relate the left Cauchy–Green strain tensor \(\mathbf {b}\), calculated at the current configuration, and the small strain tensor \(\varvec{\varepsilon }\). Indeed, for the current loading step, \(k+1\), the trial elastic deformations, pressure (\(p^{trial}_{k+1}\)), and the deviatoric stress tensor (\(\mathbf {s}^{trial}_{k+1} \)) are computed as the elastic deformations, pressure and the deviatoric stress tensor are computed as:
where K and G represent the bulk and shear moduli of the solid, respectively. \(\Delta \mathbf {F}_{k+1} \) is the incremental form of the deformation gradient, calculated as:
where \(\nabla N^a\) is the gradient of the shape function, in this case, the Local MaxEnt, defined in Sect. 4.1.
Regarding the Drucker–Prager yield criterion, the employed methodology allows to distinguish if the location of the stress state is on the cone or apex before calculating the plastic strain. The yield conditions for the classical and apex regions, respectively, are:
where \(\Delta \gamma _1=\frac{\Vert \mathbf{s} ^\mathrm{trial}_{k+1}\Vert }{2G}\), \(\Delta \gamma \) and \(\Delta \gamma _2\) are the objective functions to be calculated in the Newton–Raphson scheme for the classical or apex regions accordingly. \(c_k\) is the cohesion of the material, H the hardening parameter and \(\alpha _F\),\(\alpha _Q\) and \(\beta \) are material parameters that depend on friction and dilatancy angles as well as the shape of the yield surface, taking into account that the Drucker–Prager criterion employs a cone to approximate the Mohr–Coulomb surface and this cone can be outer or inner to the aforementioned surface (more information is found in [33]).
A limit value for the pressure, \(p_\mathrm{lim}\), is necessary to know which algorithm is to be employed. If the trial pressure is lower than this limit, classical returnmapping algorithm is employed, being this limit written as:
The equivalent plastic strain, \({\overline{\varepsilon }}^p_{k+1}\), is calculated in different ways depending on the stress state, whether it is in the classical or in the apex region:
Discretization of the solution: explicit scheme
To solve the aforementioned coupled problem in the time domain, the standard central difference explicit Newmark time integration scheme is employed. If the current time step is numbered as \(k+1\), and assuming the solution in the previous step k has been already obtained (hence it is known), a relationship between \(\mathbf{u} _{k+1}\), \(\dot{\mathbf{u }}_{k+1}\) and \(\ddot{\mathbf{u }}_{k+1}\) is established according to a finite difference scheme, as follows:
Similarly, the pore pressure, evaluated at material point level, can be expressed in terms of its derivative.
When the Newmark scheme parameters, \(\gamma \) and \(\beta \), are set to 0.5 and 0, respectively, the central difference scheme is obtained. In the present research, \(\theta =\gamma =0.5\). Rearranging terms, Predictor and Corrector terms can be obtained:
being the underlined terms the ones of the predictor step, which will be called \(\dot{u}_{k+*}\) and \(p_{\mathrm{w}_{k+*}}\). For further details, the reader is referred to Sect. 4.2.1.
About the numerical stability of the proposed methodology, it is guaranteed when the Courant–Friedrichs–Lewy (CFL) condition is satisfied. In particular, the time step, \(\Delta t\), should be small enough to ensure that the compressive wave can travel between nodes, i.e.,
where \(\Delta x\) represents the minimum mesh size, further details about the procurement of this parameter will be given in the following subsection. Finally, \(V_\mathrm{c}\) is the pwave velocity (see [55]), which is defined by
Spatial discretization
The Optimal Transportation Meshfree [17, 21, 22] has been demonstrated to perform reasonably well in geotechnical problems and, specifically, in multiphase problems [31]. It is based in the conjunction of material points and nodes. As mentioned before, the shape functions are based on the work of Arroyo and Ortiz [2], who defined the Local MaxEnt shape function (LME) of the material point \(\varvec{(}x)\) with respect to the neighborhood \(\varvec{(}x_a)\) as follows:
where the computation is done along a neighborhood \(N_b\) and
The first derivatives of the shape function can be obtained from the own shape function and its Hessian matrix J by employing the following expression:
The parameter \(\beta _{_\mathrm{LME}}\) defines the shape of the neighborhood, and it is related with the discretization size (or nodal spacing), h, and the constant, \(\gamma _{_\mathrm{LME}}\), which controls the locality of the shape functions, as follows,
It bears emphasis that \(\varvec{\lambda }^*(\mathbf{x})\) comes from the minimization of the function \(g(\varvec{\lambda })=\log Z(\mathbf{x}, \varvec{\lambda })\) to guarantee the maximum entropy. Moreover, in the remapping of the shape function, before recomputing the aforementioned minimization process, it is necessary to update the neighborhood and the parameter \(\beta _{_{\mathrm{LME},k+1}}^p < \beta _{_{\mathrm{LME},k}}^p\) in order to improve the stability.
The discretization size, h, is an interesting topic when dealing with explicit schemes in the OTM methodology. Although the neighborhood or influence radius is larger than in the traditional FEM [21], we should take h as the distance between the current node and the closest one since it is the more limiting one. Furthermore, the value of \(\Delta x\) expressed in Eq. 29 will be obtained as the minimum value of h for each node in the whole nodal set.
By employing the outlined shape functions and applying Galerkin procedure to the weak form of Eqs. (11) and (13) (See [41, 44] for details), the following matrix equations appear:
where the internal and external forces are defined as:
and the mass and damping matrices, constructed as lumped matrices in order to alleviate the computational effort of the explicit scheme, are written as follows:
being \(V_p\) and \(N_p\) the volume and the neighborhood of a material point P, respectively, \(\mathbf {B}\) the symmetric shape function gradient operator and \(\mathbf {m}\) the identity matrix in Voigt notation. Thus, \(\mathbf {B}\mathbf {m}\) reproduces the divergence operation.
Explicit integration
The proposed scheme seeks the value of the solid acceleration, \(\ddot{\varvec{u}}\), calculated from Eq. (35). It is worth mentioning that the subscript \(k+1\) is employed for the current step and k in the previous one. Furthermore, in this calculation, it is necessary to predict the internal forces from the values of the predicted solid displacement, \(\varvec{u}_{k+*}\), and the predicted pore pressure, \(p_{w_{k+*}}\). The stress has to be calculated in this predicted step as well:
Moreover, the approximation of the logarithmic strain as the measure to be employed in the deformed configuration has been demonstrated to provide good performance when large deformations are modeled (see [7, 9, 44]). In the present research, the tensor \(\mathbf {b}\), the Left Cauchy–Green strain tensor (\(\mathbf {b}=\mathbf {F}\mathbf {F}^T\)) depends on the displacement on the predicted step as follows:
Once the solid acceleration is reached, the pore pressure velocity can be calculated from Eq. (36). Also, in this equation, water internal forces and solid velocities have to be evaluated in the predicted step, \(k+*\).
All these ingredients are those which integrate the Newmark PredictorCorrector explicit algorithm for the \(\varvec{u}p_\mathrm{w}\) formulation at large strain. Its numerical implementation is explained in the following section.
Explicit algorithm within the OTM framework
The pseudoalgorithm of the whole model can be written as follows. The employment of the superscript p for material point calculations has to be pointed out.

1.
Explicit Newmark Predictor (\(\gamma =0.5\), \(\beta =0\))
$$\begin{aligned} u_{k+1}= & {} u_k+\Delta t \dot{u}_{k}+0.5\Delta t^2 \, \ddot{u}_k, \\ \dot{u}_{k+*}= & {} \dot{u}_{k}+(1\gamma )\Delta t \, \ddot{u}_{k}, \\ p_{w_{k+*}}= & {} p_{w_k}+(1\gamma )\Delta t \, \dot{p}_{w_k}. \end{aligned}$$ 
2.
Nodes and Material points position update
$$\begin{aligned}&x_{k+1}=x_{k}+\Delta u_{k+1},\\&x_{k+1}^p=x_{k}^p+\sum _{a=1}^{Nb}\Delta u_{k+1}^a \otimes N^a(x^p_{k}). \end{aligned}$$ 
3.
Deformation gradient calculation and related parameters
$$\begin{aligned} \mathbf {F}_{k+1}= & {} \Delta \mathbf {F}_{k+1} \mathbf {F}_{k}, \\ V= & {} JV_0=\det \mathbf {F} \, V_0,\\ n= & {} 1\frac{1n_0}{J}. \end{aligned}$$ 
4.
Update density and recompute lumped mass
$$\begin{aligned} \rho _{k+1}=n_{k+1}\rho _\mathrm{w}+(1n_{k+1})\rho _s. \end{aligned}$$ 
5.
Remapping loop, reconnect the nodes with their new material neighbors.

6.
Constitutive relations from the ElastoPlastic model, \(\varvec{\sigma '}_{k+*}\) and internal forces \(\varvec{R}^s_{k+*}\) and \(\varvec{R}^w_{k+*}\).

7.
Calculate \(\varvec{\ddot{u}}_{k+1}\) from Eq. (35):
$$\begin{aligned} \ddot{\varvec{u}}_{k+1} = \left[ \varvec{M}^s \right] ^{1}\left[ \varvec{R}^s_{k+*} \varvec{R}^w_{k+*} + \varvec{f}^{ext, \, s}_{k+1}\right] \end{aligned}$$ 
8.
Calculate \(\dot{p_\mathrm{w}}_{k+1}\) from Eq. (36):
$$\begin{aligned} \dot{p_\mathrm{w}}_{k+1}= \varvec{C} \dot{\varvec{u}}_{k+*} + \varvec{M}^w \ddot{\varvec{u}}_{k+1} + \varvec{f}^{\mathrm{ext}, \, w}_{k+1}\varvec{R}^\mathrm{w}_{k+*} \end{aligned}$$ 
9.
Explicit Newmark Corrector
$$\begin{aligned}&\dot{u}_{k+1}=\dot{u}_{k+*}+\gamma \Delta t \, \ddot{u}_{k+1},\\&{p}_{\mathrm{w}_{k+1}}={p}_{\mathrm{w}_{k+*}}+\gamma \Delta t \, \dot{p}_{\mathrm{w}_{k+1}}. \end{aligned}$$
Verification examples
This section is composed by two different problems. The first one deals with a consolidation, either pseudostatic or cyclic one, in order to validate the model in typical porous media applications. The second one, seeking the assessment of the performance of the proposed algorithm in a real geotechnical problem, studies the failure of a vertical wall of saturated soil.
Consolidation of a column of soil
In the following two examples, an idealization of a semiinfinite stratum of soil through a 2D column is employed, which is a traditional procedure seen in the literature. This column has a height \(H_\mathrm{T}=10\) m and a width \(L=1\) m. Lateral movements are prevented as well as the vertical movement of the rigid base. On the top, the drainage is allowed (\(p_\mathrm{w}=0\)). This geometry and boundary conditions are depicted in Fig. 1. Also, shape of both loads is depicted for the following problems, large deformation and dynamic consolidations, Sects. 5.1.1 and 5.1.2, respectively.
A regular nodal discretization of 0.5 m. size is employed, taking into account that the last top meter of the stratum is discretized with a 0.25 m. size in order to capture properly the wave provoked by the load. A similar mesh was proposed by Sabetamal et al. [40].
Large strain consolidation
Our goal in this section is the verification of the presented methodology when large deformation occurs. Considering this, the consolidation problem solved by Li et al. [23] is performed as a reference. The aforementioned geometry, seen in Fig. 1a, is adopted. The column of soil is loaded following the curve of Fig. 1b; increasing to reach \(P_\mathrm{max}\) at \(t_0=0.05\) s, when the pressure is kept constant until the end of the simulation (0.5 s). The soild and water parameters are listed in Table 1, being the NeoHookean material of Eq. (16) assumed in this case.
The verification is made against the solution proposed by Li [23]. The settlement of the top surface along time is checked for two different values of \(P_\mathrm{max}\), namely 2 and 8 MPa, that provide two different scenarios, small and large deformation regimes. The obtained solutions are seen in Fig. 2 for the two cases. Three different solutions are depicted: Static \(up_\mathrm{w}\) (Li [23]), Dynamic \(uw\) (Navas [33]) and Dynamic \(up_\mathrm{w}\) (present research). At the final of the consolidation, similar values of the settlement are achieved. Since inertial terms are included in the proposed methodology, the comparison along the entire process described by Li [23] is not possible, since in that research the quasi static \(up_\mathrm{w}\) formulation is assumed. Consequently, a ramped loading, contrary to the stepwise one employed in [23], is necessary in our case to avoid nonphysical sudden loading. Similarly, the results are not comparable against the \(uw\) formulation since fluid acceleration, neglected in the present research, was considered. In addition, the existence of displacements larger than the final settlement between 0.18 and 0.3 s. can be attributed to the above observation regarding the uw formulation. Since the application of the load is very quick, an impact fluid wave appears, whose propagation is neglected in the proposed formulation.
Additionally, the obtained settlement is compared, for both loading states (2 and 8 MPa), against the small strain solution, which was provided also by Li [23]. In Fig. 3, this comparison is plotted. As much is the deformation, as more important is the employment of the Finite Deformation regime since, as it is seen in this application, spurious results can be obtained.
Dynamic consolidation
Since soil inertial terms are considered in the proposed \(up_\mathrm{w}\) formulation, a dynamic problem has been proposed in order to see the performance of the proposed methodology. An interesting test was firstly studied by Sabetamal et al. [40] and later by Monforte et al. [28] and Navas [32, 35]. Also, the NeoHookean material is utilized. The material properties provided in Table 2, and the sinusoidal load, shown in Fig. 1c, are employed. In the aforementioned research, complete formulation (\(uwp_\mathrm{w}\) and \(uw\) ) results were provided. In this case, \(up_\mathrm{w}\) solutions of the pore pressure at different locations are presented against the stabilized \(uw\) one in Fig. 4. Slightly differences are encountered. Following, possible reasons are detailed.
On the one hand, the differences between the \(uw\) and \(up_\mathrm{w}\) solutions are small. This is due to the frequency of the load, which is not high enough to provoke water waves and, thus, the acceleration of the water phase can be neglected. Thus, we have to take into account that, following the research of Zienkiewicz and coworkers [53], the configuration of this model lies on the denominated Zone I, where dynamic terms can be neglected (See point 1 of Fig. 5). This is the reason to have similar results for both \(uw\) and \(up_\mathrm{w}\) formulations.
Zones of Fig. 5 depend on the geometry, elastic properties, frequency of the load and permeability. By fixing the rest of the parameters and tuning the frequency from 25 to 200 (Point 2 in Fig. 5) and 500 Hz (Point 2 in Fig. 5), our problem becomes Zone II and III, respectively, where dynamic terms are important. Thus, in Figs. 6 and 7, the pore pressure evolution for both approaches is presented for 200 and 500 Hz. It is noticeable the difference, since the \(up_\mathrm{w}\) is not able to capture several peaks that the \(uw\) does, more displayed for 500 Hz. Indeed, differences are more severe when it is measured deeper in the column, possibly for the undrained behavior. It must be pointed out that, for 200 Hz, no differences should be found. However, the \(up_\mathrm{w}\) solution is not able to reach \(uw\). Although point 2 is close to the border of Zone III, the figure proposed by Zienkiewicz and coworkers [53] may be updated, at least for the finite strain theory.
On the other hand, the second comparison is made in the settlement. In Sabetamal et al. [40], we find also the comparison against the analytical solution proposed by De Boer [13], corresponding to incompressible constituents. In Fig. 8, the settlement is plotted for the first 6 m from the top in two instants: 0.135 s. and 0.155 s. There is a slight difference between the peaks.
Vertical cut
In this Section, the current methodology is applied to the drainage of a square domain of saturated soil loaded on the top right half by a rigid footing. This load provokes the failure of the material in a typical vertical cut, whose shear band varies depending on the material properties, described in Sect. 3 for a hyperelastoplastic material. Precisely, the importance of this example lies in fact that, depending on the dilatancy angle, the formation of the shear band and the deformation pattern as well as the pore pressure may vary. For all the cases, the friction angle is kept at 20\(^\circ \).
The same problem was previously studied by Sanavia et al. [42, 43] and Navas et al. [33, 34] for both quasi static and dynamic regimes, respectively. The geometry and material properties are shown in Fig. 9. A displacement of 1 m on the loaded boundary, \(\Gamma _4\), is imposed gradually during the simulation, which has been fixed of 50 s. A regular 12 \(\times \) 12 nodal discretization is employed, which corresponds to a nodal spacing of 0.833 m. The time step is of 5 ms.
Results are depicted, at the final stage, in Fig. 10. In the referred bibliography, we found similar distributions of pore pressure and plastic strain for dilatant, contractive or neutral soils. However, it is worth mentioning that those results were obtained with different coupled formulation, what leads to small difference of the obtained values. Despite this fact, the trend of the behavior of the soil is well captured.
About the shear band, it can be observed that there are no big variations on the obtained peak values of the equivalent plastic strain when the dilatancy angle changes, being slightly bigger when the dilatancy angle decreases. However, an important decrease in the shear band slope is noticed when dilatancy decreases. For associate plasticity, \(\psi =20^\circ \), the shear band almost reaches the toe of the lateral wall. It should be noticed that the formation of shear bands induces to lockingbased instabilities. Those should be overcome with the appropriate techniques.
In addition, the effect of the plastic dilatancy (contractancy) is evidenced by the negative (positive) pore pressure within the shear band zone (see Fig. 10). Moreover, in the case of zero dilatancy angle, no marked pore pressure variation is observed within the shear band zone. Similar phenomena were obtained in the cited research. In order to study the evolution of the principal results of the problem, the history of the pore pressure in a material point close to the shear band (Point P, see Fig. 9) is depicted in Fig. 11.
For positive dilatancy values, smooth pore pressure evolution is observed. In addition, the peak pressure signals the initiation of plastic strain localization or shear band. The further extension of the shear band is accompanied by the continuous decreasing of the pore pressure. The material with dilatancy equal to 0\(^\circ \) experiences a softer decreasing (close to a 0\(^\circ \) slope), in this case, due to the dissipation of the pore pressure in the permeable boundary, not because of the shear band. From the same figure, it can be seen a very unstable behavior of the soil of contractive angle. This happens since the soil does not admit more load: it is completely failed. In the case of \(\psi =20^\circ \), it can be noted in Fig. 11 that the negative water pressure within the shear band is smaller than zero, indicating the possible occurrence of cavitation when lower values of − 98986 Pa (at ambient temperature) are reached. This phenomenon should be modeled by extending the formulation of this work to unsaturated conditions in order to properly capture this phenomenon.
It must be pointed out that, in this research, the sought goal is the assessment of the performance of the proposed algorithm within this geotechnical problem. Other interesting studies of the performance of the Optimal Transportation Method were carried out in [33]. Also, regular distributions provide better results. Finally, the importance of the neighborhood size was assessed, concluding that larger values of \(\gamma _{_\mathrm{LME}}\) (which corresponds to smaller neighborhood) reduce the spurious smoothing out of the shear band, being the best results obtained for \(\gamma _{_\mathrm{LME}}\)=1.4. It is worth mentioning that local processes such as plastic shear band localization are still influenced by the nodal spacing, as well as the achievement of a smoother pore pressure distribution. Following, the influence of the mesh size together with the time step is provided.
Study of the discretization size versus time step
Three different meshes, with different discretization size, have been employed: 12 \(\times \) 12 nodes (\(h=0.833\) m), 14 \(\times \) 14 nodes (\(h=0.714\) m), and 16 \(\times \) 16 nodes (\(h=0.625\) m). For this study, the case with \(\psi =0^\circ \) is employed. About the results, according as the mesh gets finer, the maximum percentage of the CFL condition able to be employed becomes smaller. Up to 25% for 12 \(\times \) 12, 20% for 14 \(\times \) 14 and 15% for 16 \(\times \) 16, the simulation becomes unstable. In Fig. 12, the computational costs of all the simulations from 5 of CFL to 25% (when it was possible) are given. It can be observed how the computational time increases exponentially as the time step decreases. Moreover, whether the discretization size gets finer the computational time increases dramatically as well.
Thus, the question of how good is the discretization size for the present problem arises. In Fig. 13, the pore pressure (in Pa) and equivalent plastic strain at 50 s of the square domain for the cases 16x16 and \(\Delta t=5\%\) CFL and 12x12 and \(\Delta t=25\%\) CFL are depicted, i.e., the largest and the shorter simulations of this study. This comparison gives us the idea of the improvement for a finer spatial discretization and smaller time steps. On the one hand, the pore pressure presents very similar distributions, what let us think that an accurate pore pressure distribution can be captured for relatively coarse discretization nodes at a relatively large time step. On the other hand, a clearer shear band is captured as the nodal size becomes finer. Moreover, since it is more localized, for the same energy, the values the equivalent plastic strain of the finer shear bands are bigger.
As it was mentioned before, up to 25% for 12 \(\times \) 12, 20% for 14 \(\times \) 14 and 15% for 16 \(\times \) 16, the simulation becomes unstable. It is important to remark this issue as one drawback of the proposed methodology when the material reaches its failure. The simulation behaves stable until the material plastifies. In Fig. 14, the distribution of pore pressure and equivalent plastic strain is seen for the final of the simulation (approximately 10 s.). The pore pressure distribution affects the constitutive behavior, leading to a spurious plastic zone. This pore pressure distribution is the typical one for undrainedincompressible materials: alternate negativepositive values. Although the problem is not in the undrainedincompressible limit, it is necessary the stabilization due to the behavior of the material once the plastification is reached. According with the obtained results of this study, future work will be lead to the stabilization of the proposed formulation.
A remark on the shear band width against the discretization size (Fig. 14) has to be pointed out. These numerical results seem to show a dependency of the shear band width on the discretization size, h. Possible reasons could be due to: (i) the shear band width of the problem is smaller than h, or (ii) the Laplacian term in the mass balance equation is not enough to regularize the strain localization problem upon h refinement. In this case, a regularization procedure could be adopted by modifying the Drucker–Prager criterion adopted in this work (e.g., Perzyna viscoplasticity) as shown in [19].
Finally, as a conclusion of this study, it must be pointed out that, for the appropriate capture of the pore pressure distribution is not necessary to refine the mesh. However, \(\Delta t\) in terms of the % of the CFL required to reach stable solutions seems to be relatively low. Moreover, if the problem needs to refine the mesh for any other reason (such as the capture of the shear band), even smaller time steps are required to reach stable results. Thus, parallelization techniques should be employed in order to get accurate results in an adequate computational time. Indeed, further work is required to reach explicit stable solutions for dynamic saturated problems at large deformations within a MaxEnt OTM framework and will be the topic for future work.
Conclusion
A new methodology to model and compute biphase saturated soils at large strains under low/medium frequency loads, by means of an Optimal Transportation Meshfree scheme with an explicit predictorcorrector time integration approach, is proposed.
The robustness of the proposed formulation is assessed by applying it to different wellknown geomechanical initial boundary value problems, with both elastic and plastic soil behavior, achieving excellent results. The first example carried out is a consolidation at large strains that was proposed firstly by Li et al. [23]. The behavior of the soil when the range of deformation is big is perfectly captured. In the second example, the model is employed under high frequency loading conditions with a hyperelastic medium. The \(up_\mathrm{w}\) formulation provides a good performance under low/medium frequency loads, but it is not well suited for high frequency loads. The model is robust and captures both displacement and pore water pressure. Zones of applicability, proposed by Zienkiewicz [55], may be revised in accordance with the results provided in this manuscript. Indeed, the validity when finite strains are employed should be assessed.
Finally, in the last case of analysis, a vertical cut is conducted for a hyperelastoplastic saturated porous material under a Drucker–Prager flow rule. The proposed model is capable of capturing the complex pore water pressure evolution within the highly distorted plastic shear band in accordance with the dilatancy of the material. Furthermore, the results obtained in the present manuscript are in agreement with the work of Sanavia [43], i.e., contractive materials accumulate pore pressure within the shear band while in the dilatant shear band a reduction in pore pressure is observed.
One of the main conclusions driven by the good performance of the proposed methodology is its extension to other particlebased numerical techniques. Previous research of the same authors, regarding the Material Point Method, shows the excellent fulfillment with Local MaxEnt shape function and an explicit predictorcorrector scheme (see [26, 27]), being both numerical techniques employed within this research. The robustness of the explicit scheme here presented encourages the authors to study other coupled formulations as well as the possibility of making dynamic relaxation techniques in order to extend the range of applicability to long simulations. Even though, since the explicit predictorcorrector time integration approach, considered in the present research, seems to capture adequately the complex hydromechanical behavior at large strains, there are other explicit time integration strategies (Runge–Kutta schemes, embedded Runge–Kutta schemes, symplectic algorithms, Taylor–Galerkinbased techniques, etc. [5, 16, 48]) with already proven capabilities in other scientific fields that should be also considered.
References
 1.
Armero F (1999) Formulation and finite element implementation of a multiplicative model of coupled poroplasticity at finite strains under fully saturated conditions. Comput Methods Appl Mech Eng 171(3–4):205–241. https://doi.org/10.1016/S00457825(98)002114
 2.
Arroyo M, Ortiz M (2006) Local maximumentropy approximation schemes: a seamless bridge between finite elements and meshfree methods. Int J Numer Meth Eng 65(13):2167–2202. https://doi.org/10.1002/nme.1534
 3.
Bandara S, Ferrari A, Laloui L (2016) Modelling landslides in unsaturated slopes subjected to rainfall infiltration using material point method. Int J Numer Anal Meth Geomech 40(9):1358–1380. https://doi.org/10.1002/nag.2499
 4.
Biot MA (1956) Theory of propagation of elastic waves in a fluidsaturated porous solid. I. Lowfrequency range. J Acoust Soc Am 28(2):168–178. https://doi.org/10.1121/1.1908239
 5.
Blanc T, Pastor M (2012) A stablized RungeKutta Taylor smoothed particle hydrodynamics algorithm for large deformation problems in dynamics. Int J Numer Methods Eng 91(June):1427–1458. https://doi.org/10.1002/nme
 6.
Borja RI, Alarcón E (1995) A mathematical framework for finite strain elastoplastic consolidation. Part1: balance laws, variational formulation, and linearization. Comput Methods Appl Mech Eng 122(94):145–171
 7.
Borja RI, Tamagnini C (1998) CamClay plasticity part III: extension of the infinitesimal model to include finite strains. Comput Methods Appl Mech Eng 155(1–2):73–95. https://doi.org/10.1016/S00457825(97)001412
 8.
Borja RI, Tamagnini C, Alarcón E (1998) Elastoplastic consolidation at finite strain. Part 2: finite element implementation and numerical examples. Comput Methods Appl Mech Eng 159:103–122
 9.
Borja RI, Tamagnini C, Amorosi A (1997) Coupling plasticity and energyconserving elasticity models for clays. J Geotech Geoenviron Eng 123(October):948–957
 10.
Cao TD, Sanavia L, Schrefler BA (2016) A thermohydromechanical model for multiphase geomaterials in dynamics with application to strain localization simulation. Int J Numer Methods Eng 107(January):312–337. https://doi.org/10.1002/nme
 11.
Ceccato F, Simonini P (2016) Numerical study of partially drained penetration and pore pressure dissipation in piezocone test. Acta Geotech 12:195–209
 12.
Cuitiño A, Ortiz M (1992) A materialindependent method for extending stress update algotithms from smallstrain plasticity to finite plasticity with multiplicative kinematics. Eng Comput 9:437–451
 13.
De Boer R, Ehlers W, Liu Z (1993) Onedimensional transient wave propagation in fluidsaturated incompressible porous media. Arch Appl Mech 63(1):59–72. https://doi.org/10.1007/BF00787910
 14.
Diebels S, Ehlers W (1996) Dynamic analysis of fully saturated porous medium accounting for geometrical and material nonlinearities. Int J Numer Methods Eng 39(1):81–97
 15.
Ehlers W, Eipper G (1999) Finite Elastic Deformations in LiquidSaturated and Empty Porous Solids. Transp Porous Media 34(1986):179–191. https://doi.org/10.1007/9789401145794_11
 16.
Hairer E, Norsett SP, Wanner G (1993) Solving ordinary differential equations I. Nonstiff problems. Springer, Berlin
 17.
Huang D, Weißenfels C, Wriggers P (2019) Modelling of serrated chip formation processes using the stabilized optimal transportation meshfree method. Int J Mech Sci 155(March):323–333. https://doi.org/10.1016/j.ijmecsci.2019.03.005
 18.
Jeremić B, Cheng Z, Taiebat M, Dafalias YF (2008) Numerical simulation of fully saturated porous materials. Int J Numer Anal Meth Geomech 32:1635–1660. https://doi.org/10.1002/nag.2347
 19.
Lazari M, Sanavia L, Schrefler BA (2015) Local and nonlocal elastoviscoplasticity in strain localization analysis of multiphase geomaterials. Int J Numer Anal Methods Geomech 39(March 2007):1570–1592. https://doi.org/10.1002/nag
 20.
Lewis RW, Schrefler BA (1998) The finite element method in the static and dynamic deformation and consolidation of porous media. Wiley, New Jersey
 21.
Li B, Habbal F, Ortiz M (2010) Optimal transportation meshfree approximation schemes for fluid and plastic flows. Int J Numer Methods Eng 83(June):1541–1579. https://doi.org/10.1002/nme
 22.
Li B, Stalzer M, Ortiz M (2014) A massively parallel implementation of the optimal transportation meshfree (pOTM) method for explicit solid dynamics. Int J Numer Meth Eng 100:40–61
 23.
Li C, Borja RI, Regueiro RA (2004) Dynamics of porous media at finite strain. Comput Methods Appl Mech Eng 193(36–38):3837–3870. https://doi.org/10.1016/j.cma.2004.02.014
 24.
LópezQuerol S, Blázquez R (2006) Liquefaction and cyclic mobility model in saturated granular media. Int J Numer Anal Meth Geomech 30(5):413–439. https://doi.org/10.1002/nag.488
 25.
LópezQuerol S, Fernández Merodo JA, Mira P, Pastor M (2008) Numerical modelling of dynamic consolidation on granular soils. Int J Numer Anal Meth Geomech 32:1431–1457. https://doi.org/10.1002/nag
 26.
Molinos M, Navas P, Manzanal D, Pastor M (2021) Local maximum entropy material point method applied to quasibrittle fracture. Eng Fract Mech 241:107394. https://doi.org/10.1016/j.engfracmech.2020.107394
 27.
Molinos M, Navas P, Pastor M, Stickle MM (2021) On the dynamic assessment of the localmaximum entropy material point method through an explicit predictorcorrector scheme. Comput Methods Appl Mech Eng 374:113512. https://doi.org/10.1016/j.cma.2020.113512
 28.
Monforte L, Navas P, Carbonell JM, Arroyo M, Gens A (2019) Low order stabilized finite element for the full Biot formulation in soil mechanics at finite strain. Int J Numer Anal Meth Geomech 43:1488–1515. https://doi.org/10.1002/nag.2923
 29.
Navas P, LópezQuerol S, Yu RC, Li B (2016) Bbar based algorithm applied to meshfree numerical schemes to solve unconfined seepage problems through porous media. Int J Numer Anal Meth Geomech 40(6):962–984. https://doi.org/10.1002/nag.2472
 30.
Navas P, LópezQuerol S, Yu RC, Pastor M (2018) Optimal transportation meshfree method in geotechnical engineering problems under large deformation regime. Int J Numer Meth Eng 115(10):1217–1240. https://doi.org/10.1002/nme.5841
 31.
Navas P, Manzanal D, Martín Stickle M, Pastor M, Molinos M (2020) Meshfree modeling of cyclic behavior of sands within large strain generalized plasticity framework. Comput Geotech 122:103538. https://doi.org/10.1016/j.compgeo.2020.103538
 32.
Navas P, Pastor M, Yagüe A, Stickle MM, Manzanal D, Molinos M (2021) Fluid stabilization of the uw biot’s formulation at large strain. Int J Numer Anal Methods Geomech 45(3):336–352. https://doi.org/10.1002/nag.3158
 33.
Navas P, Sanavia L, LópezQuerol S, Yu RC (2018) Explicit meshfree solution for large deformation dynamic problems in saturated porous media. Acta Geotech 13:227–242. https://doi.org/10.1007/s1144001706127
 34.
Navas P, Sanavia L, LópezQuerol S, Yu RC (2018) uw formulation for dynamic problems in large deformation regime solved through an implicit meshfree scheme. Comput Mech 62:745–760. https://doi.org/10.1007/s004660171524y
 35.
Navas P, Yu RC, LópezQuerol S, Li B (2016) Dynamic consolidation problems in saturated soils solved through uw formulation in a LME meshfree framework. Comput Geotech 79:55–72. https://doi.org/10.1016/j.compgeo.2016.05.021
 36.
Ortiz M, Simo JC (1985) A unified approach to finite deformation elastoplastic analysis based on the use of hyperelastic constitutive equations. Comput Methods Appl Mech Eng 49(2):221–245. https://doi.org/10.1016/00457825(85)900611
 37.
OrtizBernardin A, Puso MA, Sukumar N (2004) Construction of polygonal interpolants: a maximum entropy approach. Int J Numer Meth Eng 61(12):2159–2181. https://doi.org/10.1002/nme.1193
 38.
Pastor M, Tayyebi SM, Stickle MM, Yagüe A, Molinos M, Navas P, Manzanal D (2021) A depth integrated, coupled, twophase model for debris flow propagation. Acta Geotech Online. https://doi.org/10.1007/s11440020011144
 39.
Ravichandran N, Muraleetharan KK (2009) Dynamics of unsaturated soils using various finite element formulations. Int J Numer Anal Meth Geomech 33:611–631. https://doi.org/10.1002/nag
 40.
Sabetamal H, Nazem M, Sloan SW, Carter JP (2016) Frictionless contact formulation for dynamic analysis of nonlinear saturated porous media based on the mortar method. Int J Numer Anal Meth Geomech 40(1):25–61. https://doi.org/10.1002/nag.2347
 41.
Sanavia L, Pesavento F, Schrefler BA (2006) Finite element analysis of nonisothermal multiphase geomaterials with application to strain localization simulation. Comput Mech 37(4)(4):331–348. https://doi.org/10.1007/s0046600506736
 42.
Sanavia L, Schrefler BA, Stein E, Steinmann P (2001). In: Ehlers W (ed) Modelling of localisation at finite inelastic strain in fluid saturated porous media. Proc. Kluwer Academic Publishers, pp 239–244
 43.
Sanavia L, Schrefler BA, Steinmann P (2001) A Mathematical and numerical model for finite elastoplastic deformations in fluid saturated porous media. In: Capriz G, Ghionna V, Giovine P (eds) Modeling and mechanics of granular and porous materials. Engineering and technology, series of modeling and simulation in science. pp 297–346
 44.
Sanavia L, Schrefler BA, Steinmann P (2002) A formulation for an unsaturated porous medium undergoing large inelastic strains. Comput Mech 28(2):137–151. https://doi.org/10.1007/s0046600102778
 45.
Simo JC, Hughes TJR (2004) Interdisciplinary applied mathematics, vol 7. Computational inelasticity 79. https://doi.org/10.1086/425848
 46.
Sladek J, Sladek V, Schanz M (2014) A meshless method for axisymmetric problems in continuously nonhomogeneous saturated porous media. Comput Geotech 62:100–109. https://doi.org/10.1016/j.compgeo.2014.07.006
 47.
Terzaghi KV (1925) Engineering newsrecord. Princ Soil Mech 95:19–27
 48.
Torabi M, Pastor M, Martín Stickle M (2020) Threestep predictorcorrector finite element schemes for consolidation equation. Math Probl Eng. https://doi.org/10.1155/2020/2873869
 49.
Ye F, Goh SH, Lee FH (2014) Dualphase coupled uU analysis of wave propagation in saturated porous media using a commercial code. Comput Geotech 55:316–329. https://doi.org/10.1016/j.compgeo.2013.09.002
 50.
Zhang H, Wang K, Chen Z (2009) Material point method for dynamic analysis of saturated porous media under external contact/impact of solid bodies. Comput Methods Appl Mech Eng 198(17):1456–1472. https://doi.org/10.1016/j.cma.2008.12.006
 51.
Zhao Y, Choo J (2020) Stabilized material point methods for coupled large deformation and fluid flow in porous materials. Comput Methods Appl Mech Eng 362:112742. https://doi.org/10.1016/j.cma.2019.112742
 52.
Zheng Y, Gao F, Zhang H, Lu M (2013) Improved convected particle domain interpolation method for coupled dynamic analysis of fully saturated porous media involving large deformation. Comput Methods Appl Mech Eng 257:150–163. https://doi.org/10.1016/j.cma.2013.02.001
 53.
Zienkiewicz OC, Chan AHC, Pastor M, Paul DK, Shiomi T (1990) Static and dynamic behaviour of geomaterials: a rational approach to quantitative solutions. Part I: fully saturated problems. Proc. R. Soc. Lond. A429:285–309
 54.
Zienkiewicz OC, Chan AHC, Pastor M, Schrefler BA, Shiomi T (1999) Computational geomechanics with special reference to earthquake engineering, vol 1999. John Wiley, UK
 55.
Zienkiewicz OC, Chang CT, Bettes P (1980) Drained, undrained, consolidating and dynamic behaviour assumptions in soils. Géotechnique 30(4):385–395. https://doi.org/10.1016/j.ocecoaman.2012.02.008
 56.
Zienkiewicz OC, Shiomi T (1984) Dynamic behaviour of saturated porous media: the generalized Biot formulation and its numerical solution. Int J Numer Anal Meth Geomech 8(1):71–96. https://doi.org/10.1002/nag.1610080106
Acknowledgements
This research was funded by the Ministerio de Ciencia e Innovación, under Grant Number, PID2019105630GBI00; and the European Research CouncilH2020 MSCARISE, Grant Agreement No 101007851 (DISCO2STORE), being both greatly appreciated. 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 second author really appreciates the Entrecanales Ibarra Foundation for his undergraduate scholarship.
Funding
Open Access funding provided thanks to the CRUECSIC agreement with Springer Nature.
Author information
Affiliations
Contributions
Conceptualization and mathematical methodology were contributed by PN and MMS. Implementation was contributed by MM and PN. Validation was contributed 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.
Corresponding author
Ethics declarations
Conflict of interest
The authors declare that they have no conflict of interest.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendices
Appendix
Nomenclature

\(\varvec{a}^s \equiv \varvec{\ddot{u}}\): acceleration vector of the solid = material time derivative of \(\varvec{v}^s\)

\(\varvec{a}^{ws}\): relative water acceleration vector with respect to the solid = material time derivative of \(\varvec{v}^{ws}\) with respect to the solid

\(\varvec{b}=\varvec{F}\varvec{F}^{T}\): left Cauchy–Green tensor

\(\varvec{\overline{b}}\): body forces vector

c: cohesion (equivalent to the yield stress, \(\sigma _Y\))

\(\varvec{C}\) (time integration scheme): damping matrix

\(\frac{D^s\Box }{Dt}\) \(\equiv \dot{\Box }\): material time derivative of \(\square \) with respect to the solid

\(\varvec{F}=\frac{\partial \varvec{x}}{\partial \varvec{X}}\): deformation gradient

\(\varvec{g}\): gravity acceleration vector

G: shear modulus

h: nodal spacing

H: hardening modulus, derivative of the cohesion against time.

\(\varvec{I}\): secondorder unit tensor

\(J=\text{ det }\varvec{F}\): Jacobian determinant

k: intrinsic permeability

\(\varvec{k}\): permeability tensor

K: bulk modulus of the solid skeleton

\(K_s\): bulk modulus of the solid grains

\(K_\mathrm{w}\): bulk modulus of the fluid

\(\varvec{M}\) : mass matrix

n: porosity

\(N(\varvec{x})\), \(\nabla N(\varvec{x})\): shape function and its derivatives

p: solid pressure

\(p_\mathrm{w}\): pore pressure

\(\varvec{P}\) (time integration scheme): external forces vector

Q: volumetric compressibility of the mixture

\(\varvec{R}\): internal forces vector

\(\varvec{s}=\varvec{\sigma }^\mathrm{dev}\): deviatoric stress tensor

t: time

\(\varvec{u}\): displacement vector of the solid

\(\varvec{U}\): displacement vector of the water

\(\varvec{v}^s=\varvec{\dot{u}}\): velocity vector of the solid

\(\varvec{v}^{ws}\): relative velocity vector of the water with respect to the solid

\(\varvec{w}\): relative displacement vector of the water with respect to the solid

\(Z(\varvec{x},\varvec{\lambda })\): denominator of the exponential shape function

\(\alpha _{F}\), \(\alpha _{Q}\) and \(\beta \): Drucker–Prager parameters

\(\beta \), \(\gamma \): time integration schemes parameters

\(\beta _{_\mathrm{LME}}\), \(\gamma _{_\mathrm{LME}}\): LME parameters related with the shape of the neighborhood

\(\Delta \gamma \): increment of equivalent plastic strain

\({\overline{\varepsilon }}^p\): equivalent plastic strain

\(\varvec{\varepsilon }\): small strain tensor

\(\varepsilon _0\): reference plastic strain

\(\kappa \): hydraulic conductivity

\(\lambda \): Lamé constant

\(\varvec{\lambda }\): minimizer of log\(Z(\varvec{x},\varvec{\lambda })\)

\(\mu _\mathrm{w} \): viscosity of the water

\(\nu \): Poisson’s ratio

\(\rho \): current mixture density

\(\rho _\mathrm{w}\): water density

\(\rho _s\): density of the solid particles

\(\varvec{\sigma }\): Cauchy stress tensor

\(\varvec{\sigma '}\): effective Cauchy stress tensor

\(\varvec{\tau }\): Kirchhoff stress tensor

\(\varvec{\tau '}\): effective Kirchhoff stress tensor

\(\Phi \): plastic yield surface

\(\phi \): friction angle

\(\psi \): dilatancy angle
Superscripts and subscripts

\(^\mathrm{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

\(^\mathrm{trial}\): superscript for trial state in the plastic calculation

\(^\mathrm{vol}\): superscript for volumetric part

\(^{w}\): superscript for the fluid part relative to the solid one
Rights and permissions
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/.
About this article
Cite this article
Navas, P., Molinos, M., Stickle, M.M. et al. Explicit meshfree \({{\varvec{u}}}{{\varvec{p}}}_\mathbf{\mathrm{w}}\) solution of the dynamic Biot formulation at large strain. Comp. Part. Mech. (2021). https://doi.org/10.1007/s40571021004368
Received:
Revised:
Accepted:
Published:
Keywords
 Biot’s equations
 Meshfree
 Newmark predictorcorrector
 Explicit approach
 Large strains