A computational periporomechanics model for localized failure in unsaturated porous media

10/29/2020
by   Shashank Menon, et al.
University of Florida
0

We implement a computational periporomechanics model for simulating localized failure in unsaturated porous media. The coupled periporomechanics model is based on the peridynamic state concept and the effective force state concept. The coupled governing equations are integral-differential equations without assuming the continuity of solid displacement and fluid pressures. The fluid flow and effective force states are determined by nonlocal fluid pressure and deformation gradients through the recently formulated multiphase constitutive correspondence principle. The coupled peri-poromechanics is implemented numerically for high-performance computing by an implicit multiphase meshfree method utilizing the message passing interface. The numerical implementation is validated by simulating classical poromechanics problems and comparing the numerical results with analytical solutions and experimental data. Numerical examples are presented to demonstrate the robustness of the fully coupled peri-poromechanics in modeling localized failures in unsaturated porous media.

READ FULL TEXT VIEW PDF

Authors

page 15

page 16

page 18

page 20

page 23

page 24

page 26

page 29

03/23/2021

A stabilized computational nonlocal poromechanics model for dynamic analysis of saturated porous media

In this article we formulate a stable computational nonlocal poromechani...
08/23/2021

Computational multiphase periporomechanics for unguided cracking in unsaturated porous media

In this article we formulate and implement a computational multiphase pe...
08/25/2021

Coupling of IGA and Peridynamics for Air-Blast Fluid-Structure Interaction Using an Immersed Approach

We present a novel formulation based on an immersed coupling of Isogeome...
01/02/2020

On the modelling, linear stability, and numerical simulation for advection-diffusion-reaction in poroelastic media

We perform the linear growth analysis for a new PDE-based model for poro...
11/30/2020

An accelerated hybrid data-driven/model-based approach for poroelasticity problems with multi-fidelity multi-physics data

We present a hybrid model/model-free data-driven approach to solve poroe...
08/07/2020

An integrated numerical model for coupled poro-hydro-mechanics and fracture propagation using embedded meshes

Integrated models for fluid-driven fracture propagation and general mult...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Strain localization of geomaterials is a significant problem because it is closely related to the onset of failure of geomaterials (e.g., [1, 2, 3, 4, 5]). For example, the failure of soil slopes is a typical strain localization problem in which deformation usually occurs in a narrow banded zone with a finite thickness [6, 7]. Physically, strain localization in solid materials is characterized by the concentration of strain or displacement discontinuity in a narrow banded zone. Mathematically, strain localization is a bifurcation problem from a continuous deformation to a discontinuous one. Over the past decades, the major studies of shear banding are focused on dry or fully saturated geomaterials [2]. However, in recent years, strain localization in unsaturated soils has been the subject of tremendous interests because of the importance of practical engineering applications (e.g., unsaturated soil slope failure and geo-energy production and storage) [8, 9] (see [10] for a recent review of this subject). Unsaturated soils are three-phase porous media consisting of a deformable solid skeleton, pore water, and pore air [11, 12, 13]. Strain localization in unsaturated soils involves coupled multiphysics processes such as the coupled solid deformation and unsaturated fluid flow process (e.g., [14, 15, 16, 17]). The numerical simulations via the mixed finite element method have demonstrated that material heterogeneities (e.g., density and the degree of saturation) have a first-order triggering role in shear banding in geomaterials [5, 18]. Multiphysics processes are known triggering agents in the failure of multiphase materials such as soils [14, 15, 19]. The strain localization process usually presents multiscale and heterogeneous characteristics. Numerical methods such as the finite element method [20, 21] are useful tools to study the onset of strain localization in geomaterials [22, 23].

The finite element method is formulated in terms of the weak or variational form of partial differential equations assuming continuous field variables such as displacement and fluid pressures

[24, 25]. For this reason, the finite element method has difficulties in modeling localized failures/cracks because of the singularities across discontinuities [26]. From the mathematical viewpoint, one underlying reason is the loss of ellipticity of the partial differential equation at the onset of extreme material failures (e.g., shear bands) [27, 28, 21]. Typical regularization strategies or techniques such as rate-dependent plastic models, higher-order gradient models, nonlocal models, the strong discontinuity approach, phase-field models, the extended finite element method have been formulated for single-phase solid materials to overcome the ill-posedness problem[29, 30, 31, 32, 33, 34, 35, 36, 37]. Recently these remedies for single-phase materials have been applied to simulate shear bands or fractures in multiphase materials [38, 39, 40, 41, 42, 43, 44, 45, 46]. For instance, Réthoré et al [47, 48] formulated a finite element model for simulating fluid flow in a fracturing unsaturated porous medium in which the discontinuity is described by exploiting the partition-of-unity property of finite element shape functions. Ehlers and Luo [45] proposed a phase-field approach embedded in the theory of porous media for modeling dynamic hydraulic fracturing. Oka et al [40] developed an elasto-viscoplastic model for strain localization analysis of unsaturated soils under dynamic conditions. Roth et al [46] presented a three-dimensional extended finite element model for application to complex flow in discontinuities.

Those aforementioned enhanced finite element methods (e.g., strong discontinuities) for modeling shear bands except the nonlocal models usually assume a zero thickness of shear bands. However, the thickness of shear bands in geomaterials can be in the order of several grain sizes as observed in the laboratory testing [49, 50, 51, 52] or be of the scale of meters in the field [53, 54, 55]. Thus, in general, a length scale is required to realistically model the formation of shear bands at the continuum scale. Similarly, fluid flow in natural porous media implies that the flux at a given location may depend on a possibly large neighborhood that includes complex multiscale pathways of varying length and hydraulic conductivity. In general, there are both physical and numerical reasons for formulating a nonlocal framework for modeling variably saturated geomaterials. The physical reasons include: (i) heterogeneity of microstructures and its homogenization at a small scale on which the smoothed strain field cannot be considered as uniform [56]; (ii) grain or pore size and geometrical effects [57]; (iii) length scales of modeling fluid flow in porous media at the continuum scale through the thermodynamically constrained averaging theory and upscaling methods [58, 59]; and (iv) solid-fluid interaction at the nanoscale considering long-range forces [60, 61]. For example, de Borst et al [27]

presented a nonlocal mathematical model for strain localization using the Cosserrat continuum model of solid mechanics, wherein individual particles are enriched with three rotational degrees of freedom alongside the pre-existing three translational degrees of freedom. Zhang et al

[53] explored an internal length scale in multiphase materials such as variably saturated porous media and the relation of this length scale to the strain localization behavior under dynamic loading conditions. The nonlocal continuum concept has emerged as an effective means for regularizing the boundary value problems with strain softening, capturing the size effects, and avoiding spurious localization that gives rise to pathological mesh sensitivity in numerical computations [28]. For instance, Lazari et al [19] employed nonlocal rate-dependent plasticity to regularize the numerical solution in the strain localization process of unsaturated geomaterials.

To model localized failure in unsaturated geomaterials, a promising novel approach is peridynamics. Peridynamics was proposed to unite the mathematical modeling of discontinuities, porous media, and discrete particles within the same mathematical framework [62, 63]. Classical theory assumes that all bodies have a continuous, even distribution of mass and that all internal forces act at zero distance. In peridynamics, the governing equations are based on the concept that forces act on the material points over a finite distance. The classical theory falls apart at smaller length scales where long-range forces begin affecting material behavior. The field equations in peridynamics are integral-differential equations and are valid at discontinuities such as shear bands. Peridynamics has been applied to study the fluid flow and deformation problems in solid materials [64, 65, 66, 67, 68, 69, 70]. For instance, Turner [67] proposed a peridynamic model accounting for fluid-structure interaction for modeling fractures. Jabakhanji et al [64] proposed a bond-based peridynamic formulation for fluid flow in porous media. Katiyar et al [71] formulated a state-based peridynamic model for fluid flow in geomaterials. Ouchi et al [66] presented a coupled state-based peridynamics formulation to simulate fluid-driven fractures in a poroelastic medium. In [69], a peridynamic elastoplastic model was proposed for unsaturated geomaterials through the original correspondence principle [63]. It is noted that the formulation in [69] is for the single-phase geomaterials because the unsaturated water flow equation is missing in the mathematical model.

In this article, we implement a computational periporomechanics model for unsaturated porous media in which the unsaturated fluid flow equation is coupled with the force balance equation. The mathematical model is formulated based on the peridynamic state concept [63] and the peridynamic effective force concept for unsaturated porous media [72]. The coupled field equations are integral-differential equations (integration in space and differentiation in time) without assuming the continuity of solid skeleton displacement and pore fluid pressures. The intrinsic multiphase length scale in the field equations can be utilized to capture the length scale characteristics associated with the multiphysics behavior of unsaturated porous media (e.g., thickness of deformation band and pore pressure band). Here the intrinsic multiphase length scale means that a length scale exists in both the solid skeleton deformation and fluid flow field equations. For numerical implementation, the governing equations are discretized in space by a multiphase Lagrangian meshfree method and in time by a fully implicit method [20]. Due to its nonlocal formulation, the periporomechanics model is computationally more demanding than numerical methods based on classical local theory. To circumvent this problem, the implemented numerical model is enhanced for high-performance computing through open MPI (message passing interface) [73, 74]. To validate the numerical implementation, numerical simulations of Terzaghi one-dimensional consolidation problem [75], drainage of a soil column [76], and Mandel’s slab problem [77] are presented. Numerical simulations of localized failure in unsaturated porous media are performed through the implemented computational periporomechanics model. It is demonstrated by the detailed numerical results in Section 5 that the computational periporomechanics model is robust in simulating the formation of deformation and fluid flow bands in unsaturated porous media.

The remainder of the article is structured as follows. Section 2 introduces the unsaturated periporomechanics under quasi-static condition assuming a passive atmospheric pressure. Section 3 details the numerical implementation via an implicit multiphase Lagrangian meshfree method. Section 4 deals with validation of the numerical implementation. Section 5 presents numerical simulations of localized failure in unsaturated soils via the implemented computational periporomechanics model, followed by a summary in Section 6. Sign convention in continuum mechanics is adopted, i.e., the tensile stress in solid is positive, and pore fluid pressure in compression is positive.

2 Mathematical model

This section introduces the unsaturated periporomechanics model under quasi-static condition assuming passive air pressure. For the sake of notation and nomenclature, we briefly introduce the effective stress concept of classical unsaturated poromechanics (e.g., [12, 13]). Unsaturated soil is a three-phase porous medium consisting of a solid skeleton, liquid water, and gaseous phases. The difference between pore air and pore water pressures is usually termed as matric suction. Here, we assume passive air pressure in unsaturated soils to be at equilibrium with external atmospheric pressure and thus can neglect its effect [24]. In this case, matrix suction is merely negative pore water pressure. Here the generalized form of effective stress for the solid skeleton (e.g., [78, 79]) is adopted.

(1)

where

is the effective Cauchy stress tensor,

is the total Cauchy stress tensor, and are the pore air and pore water pressures respectively, is the Kronecker delta, and is the effective degree of saturation [79] defined as

(2)

where is the degree of saturation, and is the residual and maximum degree of saturation, respectively. Under the assumption of passive pore air pressure, equation (1) can be written as

(3)

2.1 Coupled unsaturated periporomechanics model

Peridynamics theory assumes a material is composed of material points, each interacting with its neighboring points within a nonlocal region [62, 63]. Each interaction pair of a material point with its neighboring material point in the horizon is referred to as a “bond”. Peridynamic states are mathematical objects that can apply to continuous or discontinuous functions [63]

. Peridynamic states depend upon position and time and operate on a vector connecting any two material points. Depending on whether the value of this operation is a scalar or vector, a peridynamic state is called a scalar-state or a vector-state. For instance, the deformation state is a function that maps any mechanical bond onto its image under the deformation. A force state contains the forces within bonds of all lengths and orientations. In general, a peridynamic state can be viewed as a data set to extract information about the mechanical or physical state of material points.

Figure 1: Schematic of kinematics of the coupled periporomechanics assuming atmospheric air pressure.

Figure 1 sketches the kinematics of unsaturated periporomechanics in which pore air material points are not shown under the assumption of passive air pressure. The relative position of the solid material points and in the reference configuration is given by the reference position vector denoted as the mechanical bond vector

(4)

The relative position of the same solid material points in the deformed configuration is given by the deformation position vector state

(5)

where and are the position vectors of the solid material points at positions and , respectively, in the deformed configuration. Let and are the displacements of material points at and , respectively. We have

(6)

The displacement vector states for bonds and are defined as

(7)

Substituting the displacement vector state and the reference vector state at material point into Equation (5) generates

(8)

In the unsaturated periporomechanics model, the Lagrangian coordinate system is adopted for the solid material points and the relative Eulerian coordinate system moving with respect to the solid phase is employed for the pore fluid material points (see Figure 2). Thus, the relative position of the same pore water material points in the deformed configuration is given by the deformation position vector state

(9)

where and are the position vectors of the pore water material points at positions and , respectively, in the deformed configuration. The relative deformed position vector state of the pore fluid material points with respect to the solid material points is

(10)

Furthermore, we define the pore water pressure states at material points and operating on the bonds, and , respectively, in the current configuration as

(11)

where and are pore water pressures at material points and , respectively.

Figure 2: Schematics of a solid material point (in gray and red) and a pore water material point (in blue) and their corresponding neighboring material points within the horizon: Initial state (left) and deformed state (right).

Assuming a passive pore air pressure the effective force state can be defined as

(12)

where is the total force state and is the force state of pore water, which is related to pore water pressure. It has been demonstrated that the effective force state is an energy conjugate to the deformation state of the solid skeleton [72]. It follows from (12) that the peridynamic linear momentum balance of unsaturated porous media under the quasi-static condition at can be written as

(13)

where and are the total force states at material points and , respectively, is the neighborhood of with a characteristic length , is the intrinsic density of the soil solid, is the intrinsic density of pore water, is the porosity of the mixture, and is the gravity acceleration vector.

Let and be the rates of volume change scalar state of the solid skeleton at material points and , respectively. Let and be the mass flow scalar states of pore water relative to the solid skeleton at material points and , respectively. Assuming that the soil solid is incompressible, the peridynamic mass balance equation of the mixture can be written as

(14)

where is the bulk modulus of pore water.

2.2 Constitutive and physical models

In the proposed periporomechanics model, the constitutive and physical models for the solid deformation and unsaturated fluid flow are cast through the multiphase constitutive correspondence principle [72]. The nonlocal deformation gradient can be determined from the deformation state as

(15)

where is a weighting function and is the peridynamic shape tensor [63] that is defined as

(16)

The symbol ‘’ is the peridynamic equivalent of the tensor-product and represents a summation of the product of the two quantities over the entire horizon. The nonlocal deformation gradient represents a relative motion between a material point and its neighbors in the horizon.

Next, through the multiphase correspondence principle that the effective force state at material point can be written as

(17)

where the weighting function is expressed as for brevity, is the second-order identity tensor, and is Jacobian (determinant) of . From (17) the total force state of the mixture reads

(18)

Given (18), the linear momentum balance equation of unsaturated porous media can be rewritten as

(19)

where the second term in the integrand is evaluated at material point .

The effective stress in equation (19) can be determined by a classical elastoplastic constitutive model for unsaturated soils through the multiphase constitutive correspondence principle. In unsaturated poromechanics, the degree of saturation can be determined through the soil water retention curve [80] given matric suction (i.e., negative pore water pressure). Without loss of generality, the soil water retention curve is represented by van Genuchten equation [81].

(20)

where is a scaling factor, and and are fitting parameters for the experimental data.

For self-completeness, the key components of a classical elastoplastic constitutive model are summarized as follows. The total strain in the rate form reads,

(21)

where is the total strain tensor, is the elastic strain tensor, is the plastic strain tensor, and . The effective stress is determined by the elastic model as

(22)

where is the elastic stiffness tensor and .

(23)

where and are the bulk and shear moduli of the solid skeleton, respectively. The plastic strain is determined by the plastic model. The yield surface is a function of the effective stress and the apparent pre-consolidation pressure and is expressed as

(24)

where is the mean effective stress, is the apparent preconsolidation pressure, is the deviator stress, and is the slope of the critical state line. Here and are determined by the following equations

(25)

The size of the yield function is dependent on the apparent preconsolidation pressure, which is a function of the plastic volumetric strain , matric suction, and the degree of saturation [82, 14]. The plastic strain rate is determined by

(26)

where is the so-called plastic multiplier, and .

Similarly, the non-local gradient of pore water pressure at the material point in the current configuration can be determined from the pore water pressure state as

(27)

Given the nonlocal gradient of pore water pressure, the volume flux of the pore water and at material point and point can be determined by the generalized Darcy’s law for unsaturated fluid flow [83] as

(28)

where is the intrinsic hydraulic conductivity for pore water, is the relative permeability of pore water under unsaturated condition, and is determined at material point . Given the soil water retention curve, the relative permeability of the pore water phase [84] is expressed as

(29)

where is the same parameter in Equation (20). Based on the multiphase correspondence principle for porous media [72], the mass flow state for the pore water at material points and can be defined as

(30)

It follows that the peridynamic mass balance of the fluid phase under unsaturated condition can be written as

(31)

Th first two rate terms in equation (31) are the same as in classical poromechanics and only the divergence terms are replaced by peridynamic formulations in line with a nonlocal vector calculus [63, 85]. All the input material parameters in the constitutive and physical models for the solid skeleton and fluid phases are measurable variables.

3 Numerical implementation

This section concerns the numerical implementation of the aforementioned unsaturated periporomechanics model. It is assumed in the mathematical formulation of periporomechanics that a body of unsaturated porous media consists of equal-sized mixed material points. In line with this hypothesis, the integral-differential equations of periporomechanics can be discretized by a Lagrangian meshless method in spatial domain [86] and a fully implicit method in temporal domain [20]. In the proposed method, the problem domain is discretized into a finite number of equal-sized mixed material points (i.e., pore fluid and solid skeleton). Each material point has two kinds of degrees of freedom, i.e., displacement and pore fluid pressure. Thus, the meshless method is different from the one in [87] or [72], which is for modeling single-phase solid materials. We refer to the literature for an in-depth report on meshfree methods [88, 89]. The message passing interface is utilized to enhance the implemented computational periporomechanics model for high-performance computing [73, 74].

3.1 Implicit time integration

An implicit time integration for is employed for the unsaturated periporomechanics framework. The fundamental unknowns at each mixed material point are displacement vector and pore water pressure . Given the values of displacement vector and pore water pressure at time step , the implicit time integration [20, 21] is used to determine the displacement vector and pore water pressure at time step . At time step , we can write out the momentum balance as

(32)

where the variables at mixed material points and are evaluated at time step , and the subscript is omitted for brevity. The effective force state at time step is determined by the effective stress through Equation (12). Given the nonlocal deformation gradient, the effective stress is computed through the elastoplastic constitutive model introduced in the previous section. In line with the global implicit time integration, this elastoplastic constitutive model is numerically implemented through an implicit algorithm (e.g., [90, 91, 92, 15]). We refer to [93] for the implementation of the numerical algorithm.

The backward finite difference is adopted to approximate the rate terms of the mass balance equation in time domain. Then we have

(33)

where a subscript for variables at time step is omitted for brevity, and variables at time step are noted with a subscript . Multiplying on both sides of Equation (33) generates

(34)

3.2 Linearization

The tangent stiffness matrix is required when solving nonlinear equations by Newton’s method. In general, the tangent stiffness matrix is a function of the chosen linearization point and must be updated every time step when the nodal degrees of freedom (displacement and pore water pressure) are updated. The tangent stiffness matrix relates an infinitesimal disturbance in the degree of freedom to change in the global force state and the global fluid flux state. For the state-based periporomechanics the tangent stiffness matrix contains non-zero entries for each node that interacts directly as well as nodes that share a neighbor material point, which makes the Jacobian inherently denser than that of a local poromechanics model. Figure 3 sketches the incremental kinematics for the solid phase and the pore fluid phase, respectively.

Figure 3: Incremental kinematics for the solid phase and the pore water pressure of the fluid phase. Note: and are iteration numbers at time step .

The linearization of the momentum balance equation with respect to the displacement vector states and , and pore water pressures and at time step assuming an infinitesimal deformation is

(35)

Similarly, the linearization of the mass balance equation at time step is

(36)

3.3 Multiphase meshfree spatial discretization

Figure 4 sketches the two-phase meshfree discretization of the problem domain.

Figure 4: (a) Schematic of the multiphase meshfree spatial discretization of an unsaturated porous material body, and (b) zoom-in view of material point and its neighbors within the horizon . Note: the gray dot represents the solid skeleton and the blue circle represents the pore water.

For a periporomechanics material point , there are totally neighbor material points within its horizon. Here by the periporomechanics material point we mean that the material point has two states - the solid skeleton state and the pore fluid state. Therefore, each periporomechanics material point is endowed with two kinds of degree of freedom, the displacement of the solid skeleton and the pressure of the pore water. The momentum balance equation at material point becomes

(37)

where is the volume of material point within the horizon of material point . Similarly, the mass balance equation at material point is

(38)

At time step , the linearization of the momentum balance is

(39)

Similarly, the spatial discretization of the linearized mass balance equation is

(40)

The tangent matrix must be assembled by looping over all the individual particles and neighbors and account for the effect of perturbations in all degrees of freedom in every single bond on the force state and the fluid flux state at the point. The linearization must be considered for the periporomechanics material point itself as well as all its neighbors for each degree of freedom such as the displacement and the pore water pressure. Parallel computing capability of the implemented numerical model is realized through Open MPI, an open-source MPI (message passing interface) library

[87, 86]. Here MPI is a communication protocol for programming parallel computers. We refer to [73, 74] for the technical details regarding the implementation of Open MPI for high-performance computing.

4 Validation of the numerical implementation

In this section we validate the implemented computational periporomechanics model by conducting numerical simulations of classical poromechanics problems including Terzaghi’s one-dimensional consolidation problem, the drainage of a soil column, and Mandel’s slab problem [75, 76, 77].

4.1 Terzaghi’s one-dimensional consolidation

The Terzaghi’s one-dimensional consolidation problem is simulated, and the numerical result is compared against the closed-form solution. For the specimen, the height is 20 cm, the width is 5 cm, and the out-of-plane thickness is 1 cm. The soil skeleton is assumed as a homogeneous elastic material with bulk modulus kPa and shear modulus = 1.8 kPa. The fluid density is 1000 kg/m and is assumed incompressible. Hydraulic conductivity of the soil skeleton is m/s and the initial porosity is = 0.5. Figure 5 sketches the model set-up and the boundary conditions of the fully saturated soil specimen. Both lateral boundaries are constrained against horizontal deformation and are free to deform vertically. The bottom boundary is fixed in all three directions. The problem domain is discretised into 400 material points with a characteristic dimension of = 0.5 cm. The horizon is assumed as .

Figure 5: Schematic of the geometry and boundary conditions of the Terzaghi consolidation problem.

The initial and boundary conditions for the fluid phase are as follows. At the initial state, all boundaries for the fluid phase are impervious. The top boundary is prescribed with a zero fluid pressure (i.e., the atmospheric pressure) during the consolidation process. A surcharge in the form of body force, equivalent to 50 kPa, is applied to the top boundary at t = 0 s and is held constant during the simulation. The total simulation time is 2000 s with an equal time step = 1 s. Figure 6(a) plots the excess pore water pressure from the numerical simulation against the closed-form solution [75] at different consolidation stages. Figure 6(b) shows the convergence of the numerical simulation to the analytical one with an increasing nonlocal parameter ‘m’: the ratio of the horizon to the characteristic length of the material point (i.e., ). It is demonstrated from Figure 6 that the numerical solution is in good agreement with the closed-form solution.

Figure 6: (a) Variations of the excess water pressure versus the height of the soil specimen from the numerical and closed-form solutions at different consolidation times. (b) Convergence of the periporomechanics solution to the closed-form solution with increasing nonlocal parameter. Note: T is a dimensionless time factor [75].

4.2 Drainage of a soil column

In this example, the experimental test of drainage of a sand column [76] is simulated numerically. The numerical results are compared with the experimental data. The sand column is 1 m in height and 0.25 m in width. The sample is discretized into uniform material points with a characteristic length of = 0.2 m. The horizon is assumed as . For the material parameters, Young’s modulus is 1.3 MPa, Poisson’s ratio is 0.4, soil density is 2000 kg/m, water density is 1000 kg/m, porosity is 0.2975, intrinsic permeability is 4.5 10m. From the data provided by Liakopoulos, the degree of saturation and the relative permeability of water are determined by the following equations [24].

(41)
(42)

For the solid phase, the lateral surfaces are constrained against lateral displacement but free to deform vertically, and the bottom is constrained for vertical displacement. For the fluid phase, both the upper and lower ends of the column are exposed to atmospheric pressure, and zero flux conditions are imposed on the lateral surfaces. The sand column is initially fully saturated with water. The water drainage is solely driven by the gravity load. The results are plotted in Figures 7 and 8. In Figure 7, the pore water pressure from numerical simulations at different times is compared to the experimental results. The numerical results are in good agreement with the experimental data except for the early-stage results. In Figure 8, the fluid flow velocity at the drainage boundary predicted by the numerical simulation is compared against the experimental result. It is shown that the numerical result is close to the experimental data. It shall be noted that because the input parameters adopted from the literature are approximate ones the numerical results can not match the experimental data exactly.

Figure 7: Comparison of the pore water pressures from the numerical simulation and experimental data [76] in the soil column at different times.
Figure 8: Comparison of the fluid flow velocities from the numerical simulation and experimental data [76] at the drainage boundary over time.

4.3 Mandel’s slab problem

In this example, Mandel’s slab problem is simulated to examine if the implemented computational periporomechanics model can capture the so-called Mandel-Cryer effect [94, 84] under saturated condition. Here the Mandel-Cryer is a non-monotonic pore fluid pressure response that can only be captured by a fully coupled poromechanics solution. This phenomenon was first studied by Mandel with a saturated isotropic poroelastic specimen [77] and later by Cryer with an isotropic poroelastic sphere under hydrostatic pressure [95]. Both studies showed that the induced pore fluid pressure at the center of the specimen is higher than the applied pressure during the initial loading stage. The physical reason behind this phenomenon is that fluid flow and pressure dissipation are retarded by the permeability and the distance to the drainage boundary. In contrast, the applied load and pore pressure response may be immediate and constant. Thus, when the pore fluid pressure near the drainage boundary starts to dissipate, the consequent rapid and localized deformation leads to an instantaneous load transfer to other parts of the porous body, which may lead to the development of additional excess pore pressure.

Figure 9: Sketch of the model setup, loading protocol and the boundary conditions. Note: the rollers are free to move horizontally and stands for the external load.

For Mandel’s slab problem, an infinitely long rectangular saturated soil specimen with free drainage conditions at both lateral boundaries is suddenly subjected to a constant vertical surcharge through a rigid and friction-free plate of the same width as the soil sample [77] (see Figure 9). The same material parameters for the Terzaghi’s consolidation in 4.1 are adopted. The surcharge is equal to the initial excess pore pressure in the specimen. The lateral boundaries are stress-free and can deform freely. The deformation of the specimen is under two-dimensional conditions. The dimension of the specimen is 10 m in width, 3 m in height, and 0.4 m in out-of-plane thickness. Figure 9 sketches the problem domain. The specimen is discretized into 1500 equal-sized material points with a characteristic length scale of = 0.2 m. The horizon is assumed as = 0.4 m. The initial excess pore fluid pressure is 1 kPa throughout the specimen. At t = 0, both lateral boundaries are forced to zero pressure, and a constant vertical surcharge equivalent to 1 kPa is applied to the top boundary. The numerical result is compared against the closed-form solution [77, 96]. Figure 10 plots the variation of the water pressure at the center of the Mandel’s slab obtained by the numerical simulation and the closed-form solution. It is shown that the numerical result generated by the proposed computational periporomechanics model are in good agreement with the closed-form solution. Thus, it can be concluded that the fully coupled computational periporomechanics model can capture the Mandel-Cryer effect in the classical coupled poromechanics problem.

Figure 10: Comparison of the numerical and closed-form solutions of the water pressure at point A at the center of Mandel’s slab.

5 Numerical examples

Numerical examples of localized failure in unsaturated soils under different scenarios are presented in this section. Example 1 is a shear and compression test of a two-dimensional unsaturated soil sample. Example 2 concerns strain localization in a simple compression test under two-dimensional condition. Examples 1 and 2 study the influence of the horizon on shear bands. Example 2 also analyzes the impact of the intrinsic hydraulic conductivity and spatial discretization on the numerical results. Example 3 simulates the bearing capacity of a footing on unsaturated soils. For all examples, the same intrinsic length scales (i.e., horizon) are adopted for both the fluid and solid phases. For simplicity, the unit weighting function is used for all terms in the mathematical model.

5.1 Example 1: Shear and compression of a two-dimensional specimen

In this example, we simulate a two-dimensional unsaturated soil specimen under combined shear and compression loading. As shown in Figure 11 (a), the soil specimen is 1 cm by 1 cm in the plane and the out-of-plane thickness is 0.025 cm. The sample is discretized into 12,800 material points with a characteristic length = 0.125 mm. Figure 11(b) sketches the spatial discretization. The horizon is 0.25 mm and the thickness of boundary layers are 0.25 mm. The material parameters for the solid skeleton are assumed as follows: bulk modulus is 38.9 MPa, shear modulus is 13 MPa, swelling index is 0.03, slope of critical state line is 1.0, compression index is 0.12, and specific volume at unit preconsolidation pressure is 2.76, and fitting parameters for preconsolidation pressure [69] are 0.185 and 1.42. The parameters for soil water characteristic curve are: = 0.0, = 1.0, = 20 kPa, = 2.0, = 0.5. The soil density is 2000 kg/, water density is 1000 kg/, initial porosity is 0.5, and hydraulic conductivity is 1 m/s. The same material parameters are also adopted for the numerical examples presented in Section 5.2.

The bottom of the specimen is fixed. All fluid phase boundaries are no-flux conditions. For the initial state, a uniform isotropic stress state with a magnitude of -100 kPa is prescribed on the skeleton, and apparent pre-consolidation pressure is -300 kPa. In the sample, the initial suction is assumed 20 kPa, and the degree of saturation is 0.8. The constant confining pressure on the lateral boundaries is 84 kPa, which is imposed through an equivalent total force density [93]. For the loading protocol, the soil specimen is subjected to a vertical compression load and an in-plane shear load on the top boundary. Both loads are prescribed by constant displacement rates of m/s in the direction and m/s in the direction, respectively (see Figure 11). The confining pressure is applied on both lateral boundaries by a linear ramp function in 0.04 s. Then, the external loads are imposed on the top boundary with the prescribed constant rates. The total simulation time is 24 s.

Figure 11: (a) Model setup, and (b) spatial discretization (the problem domain in blue and the boundary layers in gray).
Figure 12: Reaction force and its direction versus the vertical and horizontal displacements at the top boundary. Note the markers at displacements: = 0.13 mm and = 0.24 mm
Figure 13: Contours of equivalent plastic shear strain at displacements: (a) = 0.13 mm and (b) = 0.24 mm.
Figure 14: Contours of plastic volumetric strain at displacements: (a) = 0.13 mm and (b) = 0.24 mm.

Figure 12 plots the reaction force versus the vertical and horizontal displacements at the top boundary, as well as the direction of the reaction force and its eventual variation with increasing displacement. The markers in Figure 12 denote the vertical displacements: = 0.14 mm and = 0.24 mm, respectively. The results of plastic strains, pore water pressure and fluid flux vector for the load steps are shown in Figures 13, 14, 15 and 16. Figure 13 shows the contours of equivalent plastic shear strain (i.e., the second invariant of the plastic deviator strain tensor [14]). It is noted that the boundary layers are omitted here. Figure 14 plots the contours of plastic volume strain. The results in both figures demonstrate that a single shear band with a finite thickness has been initiated and developed in the specimen. The soil skeleton experiences positive plastic volumetric strain in the localized shear zone. Thus, it can be concluded that the shear band is a dilative one. Figure 15 plots the pore water pressure at the selected load steps. The results in Figure 15 show that a single banded zone has formed in the contour of pore water pressure. The banded zone of pore water pressure is consistent with the shear band of the solid skeleton shown in Figures 13 and 14. The decrease in the pore water pressure inside the shear band (equivalently, suction increases) can be in part due to the plastic dilation in the banded deformation zone.

Figure 15: Contours of pore water pressure (kPa) at displacements: (a) = 0.13 mm and (b) = 0.24 mm
Figure 16: Contours of fluid flux vector and the magnitude at displacements: (a) = 0.13 mm and (b) = 0.24 mm

Figure 16 plots the snapshots of the fluid flux vector on the undeformed configuration at the selected load steps. The results demonstrate that the water flows into the center of the specimen at the inception of shear banding. After the shear band occurs, more water flows into the banded zone. This can be due to the fact that the pore water pressure in the banded zone is lower than the pore water pressure outside the shear band (refer to Figure 15). As the plastic volumetric strain (see Figure 14) increases in the banded zone, the pore water pressure decreases and the magnitude of the water flux rises because of the increased difference (or discontinuity) between the pore water pressures within and outside the shear band. It can be concluded from the example that the proposed periporomechanics model can capture the finite thickness of the shear band and the localized zone of the water pressure in the specimen.

5.1.1 Influence of horizon

The horizon in the proposed periporomechanics is an intrinsic length scale for both the fluid and solid phases. The impact of the horizon on the formation of shear banding under shear and compression conditions is examined here. In so doing, we re-run the simulation in Section 5.1 with two more horizons, = 0.375 mm and =0.5 mm while keeping all other simulation parameters unchanged. Figure 17 compares the loading capacity curves (the reaction force versus displacement) generated by the simulations with the three horizons, respectively. The results show that the peak reaction forces (or applied loads) are different for the simulations with these three horizons. It is evident that the force responses are almost identical at the early loading stages. After the peak load, the loading curves diverge, which demonstrates the load capacity of the specimen in the post-localized regime is sensitive to the intrinsic length scale (horizon). Specifically, the simulation with a larger horizon (e.g., =0.5 mm) generates a relatively larger loading capacity in the post-localization regime. Figures 18 and 19 plot the contours of equivalent plastic shear strain and pore water pressure at the vertical displacement = 0.24 mm, respectively. The results obviously show that the internal length scale impacts both the formation of shear deformation band and pore pressure band. The simulation with = 0.25 mm (the smallest horizon in all three cases) produces the narrowest shear band and the maximum intensity of strain and pore water pressure within the shear band in all three cases. From the results in Figures 17 - 19, we may conclude that both the deformation pattern, pore water pressure distribution, the load capacity of the specimen are dependent on the horizon size in the numerical model. For instance, under the same other conditions, the simulation with a larger horizon generates less intense plastic deformation in the localized zone than the simulation with a smaller horizon. Furthermore, the results show that the horizon may be correlated to the width of the shear band. It may imply that the proposed periporomechanics model can better capture the thickness of shear bands in unsaturated soils by adjusting the horizon.

Figure 17: Comparison of reaction force for the cases of = 0.25 mm, = 0.375 mm and = 0.5 mm.
Figure 18: Contours of equivalent plastic shear strain at the applied displacement = 0.24 mm for (a) = 0.25 mm, (b) = 0.375 mm, and (c) = 0.5 mm.
Figure 19: Contours of pore water pressure (kPa) at the applied displacement = 0.24 mm for (a) = 0.25 mm ,(b) = 0.375 mm, and (c) = 0.5 mm.

5.2 Example 2: Compression of a two-dimensional specimen

In this example, we conduct a compression test of a two-dimensional unsaturated soil specimen with a constant confining pressure. Figure 11 (a) sketches the problem domain and the displacement loading condition. The specimen has the same geometrical size as the one for Example 1 (Figure 11 (a)) in 5.1. The bottom is fixed in all directions. Both lateral boundaries are prescribed with constant confining pressure. All fluid phase boundaries are no-flux. The specimen is discretized into 12,800 equal-sized material points with a characteristic length = 0.125 mm. The horizon is 0.25 mm. For the initial state, a uniform isotropic stress state of -100 kPa is prescribed on the skeleton, and apparent pre-consolidation pressure is -300 kPa. The initial suction is 20 kPa. Equivalently, the pore water pressure is -20 kPa. The constant confining pressure on the lateral boundaries is 84 kPa, which is imposed through an equivalent total force density [93]. For the loading protocol, the confining pressure is applied by a linear ramp function in 0.04 s. Then, the top boundary is subjected to a vertical compression with a rate of m/s. The total simulation time is 24 s.

Figure 20: Reaction force versus the vertical displacement at the top boundary.

Figure 20 plots the reaction force versus the vertical displacement of the top boundary. The markers in Figure 20 correspond to two applied displacements: (a) = 0.18 mm and (b) = 0.72 mm, respectively. The peak reaction force occurs at = 0.18 mm. Simulation results at both load stages are shown in Figures 21, 22, 23 and 24. All contours depicted here are drawn on the deformed configuration of the specimen. Figure 21 presents the snapshots of equivalent plastic shear strain. The results show the formation of two conjugate shear bands in the specimen. Figure 22 presents the snapshots of the plastic volumetric strains. The results in Figure 22 show that positive plastic volumetric strains occur in localized shear zones. Figure 23 presents the snapshots of pore water pressure. Similar to the results for the deformation field in Figure 21 and 22, two conjugate banded zones have formed in the pore water pressure field. Within the shear band, the pore water pressure is smaller than that outside the shear band. It is also evident that both the deformation band and the fluid pressure band have a similar finite thickness in that the same horizon is adopted for both the solid and the fluid. Figure 24 shows the snapshots of the fluid flux vector and the magnitude at the different load stages. The results show that the pore water flows into the shear band. This is due to the fact that the pore water within the shear band decreases (or equivalently suction increases) in the formation of dilative shear bands.

Figure 21: Contours of equivalent plastic shear strain at displacements:(a) = 0.18 mm and (b) = 0.72 mm.
Figure 22: Contours of plastic volume strain at displacements: (a) = 0.18 mm and (b) = 0.72 mm.
Figure 23: Contours of pore water pressure (kPa) at displacements: (a) = 0.18 mm and (b) = 0.72 mm.
Figure 24: Contours of fluid flow vector and the magnitude at displacements: (a) = 0.18 mm and (b) = 0.72 mm.

5.2.1 Influence of horizon

We study the influence of the horizon on the formation of shear bands under the symmetrical loading condition. In so doing, the previous simulation is repeated with two more horizons, = 0.375 mm and = 0.5 mm while keeping all other parameters unchanged. Figure 25 compares the loading capacity curves (i.e., the applied loads versus the vertical displacement the top boundary) generated by the simulations with the three horizons, respectively. The results show that the loading response curves are slightly sensitive to the horizon size in the post-localization regime. The simulation with a larger horizon size generates a higher loading response at the same load step.

Figure 25: Comparison of force versus applied vertical displacement for = 0.25 mm, (b) = 0.375 mm, and (c) = 0.5 mm.
Figure 26: Contours of pore water pressure at = 0.72 mm for for (a) = 0.25 mm, (b) = 0.375 mm, and (c) = 0.5 mm.
Figure 27: Contours of equivalent plastic shear strain at = 0.72 mm for (a) = 0.25 mm, (b) = 0.375 mm, and (c) = 0.5 mm.

Figures 26 and 27 show the contours of pore water pressure and plastic shear strain at the last load step of the simulations with the three horizons. Similar to the results of the shear example, the horizon size can influence the deformation and pore water pressure in the post-localization regime (e.g., the thickness of shear bands and the magnitude of the pore water pressure in the shear band). For instance, the simulation with the smallest horizon size generates the most intense plastic shear strain within the shear bands.

5.2.2 Influence of hydraulic conductivity

We analyze the influence of the intrinsic hydraulic conductivity on the numerical results. Assuming a constant viscosity (i.e., isothermal condition) and density (i.e., low water pressure) of pore water [83], this analysis is equivalent to studying the influence of intrinsic permeability on the numerical results. For this purpose, we re-run the simulation in Section 5.2 with three different values of hydraulic conductivity: (a) = 3 m/s, (b) = 3 m/s, and (c) = 3 m/s. Figure 28 plots the reaction force versus the vertical displacement for the three simulations. Figure 29 portrays the contour of pore water pressure at = 0.72 mm for the three simulations. Figure 30 presents the contours of equivalent plastic shear strains at = 0.72 mm.

Figure 28: Plot of reaction force versus vertical displacement for simulations with three intrinsic hydraulic conductivities: (a) = 3 m/s, (b) = 3 m/s, and (c) = 3 m/s.
Figure 29: Contours of pore water pressure (kPa) at = 0.72 mm for simulations with three intrinsic hydraulic conductivities: (a) = 3 m/s, (b) = 3 m/s, and (c) = 3 m/s.
Figure 30: Contours of equivalent plastic shear strain at = 0.72 mm for simulations with three intrinsic hydraulic conductivities: (a) = 3 m/s, (b) = 3 m/s, and (c) = 3 m/s.

The contours of water pressure show that the hydraulic conductivity affects the magnitude of the water pressure in the post-localization regime. For example, the results in Figure 29 the simulation with a smaller value of intrinsic hydraulic conductivity generates a more considerable pressure change inside the shear band. Furthermore, the results in Figure 30 show that the intrinsic permeability has a negligible impact on the thickness of shear bands as well as the magnitude of equivalent plastic shear strain within the shear band. This observation may be due to the strong nonlocality for the fluid and solid phase in the coupled periporomechanics model. Note further study on the aspect shall be conducted in the future.

5.2.3 Influence of spatial discretization

To examine the influence of spatial discretization on the numerical results, we rerun the simulation in Section 5.2 with two spatial discretization schemes: (a) the coarse one with 3200 material points ( = 0.25 mm) and (b) the fine one with 12,800 material points ( = 0.125 mm), respectively. To compare the results, the same horizon = 0.375 mm is adopted for both simulations. Figure 31 plots the reaction force over the vertical displacement at the top boundary for the simulations with two spatial discretization schemes. It is apparent that the force predicted in both simulations is almost identical.

Figure 31: Plot of reaction force versus vertical displacement with two different spatial discretization schemes.
Figure 32: Contours of pore water pressure (kPa) at = 0.72 mm using (a) coarse discretization and (b) fine discretization in spatial domain.
Figure 33: Contours of equivalent plastic shear strains at = 0.72 mm using (a) coarse discretization and (b) fine discretization in spatial domain.

Figure 32 depicts the contours of water pressure superimposed on the deformed configuration at = 0.72 mm. Figure 33 plots the contours of equivalent plastic shear strains at = 0.72 mm. The results in both figures evidence that the equivalent plastic shear strain and water pressure are insensitive to spatial discretizations. For both cases, the plastic shear deformation and the water pressure localize in a well-defined zone of finite thickness that is correlated to the horizon. It may be concluded that the proposed coupled periporomechanics is not sensitive to the spatial discretization, which is different from the standard finite element method. Thus, the proposed periporomechanics can be a valuable tool in the analysis of localization phenomena in unsaturated porous media.

5.3 Example 3: Bearing capacity of a footing on unsaturated soils

This example deals with the bearing capacity of footing on unsaturated soils. Specifically, this example aims to show the proposed computational periporomechanics model can qualitatively simulate the localized failure mode of unsaturated soils under a footing. Figure 34 sketches the geometry of the problem domain and the applied boundary conditions. A footing on which the displacement load is imposed is placed on the top of the soil domain. The dimension of the soil domain is 4 m 1 m 0.02 m. The footing’s dimension is 0.8 0.01 0.02 m. As shown in Figure 34, the lateral boundaries are constrained from horizontal deformation and are free to deform vertically. The bottom is constrained from vertical deformation. All fluid phase boundaries are no-flux conditions. The problem domain is discretized into 60,000 mixed material points. The multiphase length scale of is assumed 0.022 m.

Figure 34: Problem setup.

The material parameters for the skeleton are as follows. Bulk modulus is 55.56 kPa, the shear modulus is 18.52 kPa, swelling index is 0.05, the slope of critical state line slope is 1.0, compression index is 0.15, and specific volume at unit pre-consolidation pressure is 2.0, initial apparent pre-consolidation pressure is -200 kPa, and other parameters are the same adopted in examples 1 and 2. The parameters for soil-water characteristic curve are: = 0.0, = 1.0, = 10 kPa, = 1.25, and = 0.2. The soil density is 2500 kg/, water density is 1000 kg/, initial porosity is 0.25, and hydraulic conductivity is 1 m/s. The initial suction in the soil is assumed 20 kPa. The corresponding degree of saturation and relative permeability are 0.8 and 0.01, respectively. The loading protocol is divided into two steps. In step 1, the soil’s initial geostatic stress state is prescribed by imposing a gravity load on the soil. Then, a null deformation condition is assigned to the soil before a load is imposed on the footing. In step 2, the footing is given a vertical displacement = 0.04 m at the constant rate of 1.67 mm/s in 15 s.

Figure 35: Plot of the reaction force on the footing versus the applied vertical displacement.
Figure 36: Contour of (a) the magnitude and (b) direction of the displacement vector of the soil at = 2.5 cm.
Figure 37: Contours of (a) equivalent plastic shear strain, (b) plastic volumetric strain, and (c) pore water pressure (kPa) at = 2.5 cm.

Figure 35 plots the reaction force of the footing versus the vertical displacement. As shown in Figure 35, the bearing capacity of the footing is around 7.5 kN. The contour of the magnitude of displacement is presented in Figure 36 (a). The displacement vector is plotted in Figure 36 (b) in which the arrowhead denotes the direction, and the length of the arrow tail shows the scaled magnitude. The contours of equivalent plastic shear strain, plastic volumetric strain (or plastic dilatation), pore water pressure in the soil at = 2.5 cm are depicted in Figure 37 (a), (b) and (c). The results show that a typical localized shear failure mode of soil under shallow foundations [75]. For instance, the settlement of the footing has caused the formation of symmetric slip surfaces in the soil and soil upheaval on both sides of the footing. It is apparent from Figure 37 (b) that the shear bands under the footing are compactive while those associated with soil upheaval are dilatant. It is also revealed from Figure 37 (c) that matric suction decreases (pore pressure rises) in the compactive shear bands beneath the footing while matric suction increases in the dilative shear bands. The results are corroborated by the data shown in Figures 38 and 39 which present the variation of the equivalent plastic strain and pore water pressure versus the footing’s vertical displacement at points inside and outside the shear bands denoted in Figure 36 (a).

Figure 38: Variation of (a) equivalent plastic shear strain and (b) pore water pressure versus the footing’s vertical displacement at the points A and B (inside and outside banded compaction zones, respectively).
Figure 39: Variation of (a) equivalent plastic shear strain and (b) pore water pressure at the points C and D versus the footing’s vertical displacement (inside and outside banded dilation zones, respectively).
Figure 40: Plot of reaction force on the footing over the applied vertical displacement for Simulation I and Simulation II.
Figure 41: Contours of equivalent plastic shear strain for (a) Simulation I and (b) Simulation II at = 2.5 cm.
Figure 42: Contours of pore water pressure for (a) Simulation I and (b) Simulation II at = 2.5 cm.

To test the solution’s uniqueness, we rerun the simulation with a relatively finer spatial discretization consisting of 100,000 mixed material points. The same length scale and input parameters are adopted. The comparison between the two simulations is presented in Figures 40, 41 and 42. Here simulations I and II stand for the simulations with 60,000 and 100,000 material points respectively. Figure 40 shows that the reaction force curves from both simulations are almost identical. Moreover, as shown in Figures 41 and 42 respectively the contours of equivalent plastic shear strain and pore water pressure at = 2.5 cm from both simulations are in good agreement. It is noted that all the numerical simulations are run in parallel on a supercomputer. To show the performance of the Open MPI implementation and scalability of the numerical code, the total wall-clock time for both simulations run in parallel with different number of CPU (central processing unit) processors is summarized in Table 1. With the same number of CPU processors, the wall-clock time of simulation I with 60,000 multiphase material points is much shorter than that of simulation II with 100,000 multiphase material points. Doubling the number of CPU processors reduces the wall-clock time dramatically for both simulations, although the reduction of time is less than half. The latter can be because more message passing occurs between processors for the simulations with 64 CPUs than those with 32 CPUs.

Number of CPU Processors Wall-clock Time
Simulation I (60,000 material points) Simulation II (100,000 material points)
32 15,600 s 44,500 s
64 8,900 s 33,000 s
Table 1: Comparison of wall-clock time for simulations I and II

6 Closure

In this article, we have implemented a computational periporomechanics model for simulating localized failure in unsaturated porous media. The mathematical model is based on the peridynamic state concept and the effective force state concept. It is assumed that the pore air pressure in unsaturated porous media is passive (i.e., atmospheric pressure). The coupled governing equations are integral-differential equations without assuming the continuity of solid displacement and pore water pressure (or matric suction). Both equations have an intrinsic two-phase length scale - the horizon. In the coupled periporomechanics model, it is assumed that unsaturated porous media consist of equal-sized mixed material points, which have two kinds of degrees of freedom, i.e., the displacement and pore water pressure. In line with this hypothesis, we have proposed an implicit multiphase Lagrangian meshless method to numerically implement the periporomechanics model. The numerical implementation is validated by simulating Terzaghi’s consolidation problem, drainage of a soil column, and Mandel’s slab problem. Numerical examples under different scenarios are conducted to demonstrate that the computational periporomechanics model is robust in modeling the formation of strain localization in unsaturated porous media. It is shown that the intrinsic multiphase length scale can be utilized to characterize the finite thickness of the localization of deformation and pore water pressure. Length scale sensitivity analysis has shown that the multiphase length scale can impact the formation of shear bands and the loading capacity of unsaturated soils in the post-localization regime. It was found from the numerical analysis that spatial discretization and intrinsic hydraulic conductivity have mild effects on localized failure of unsaturated soils. This may be due to the fact that the identical length scale (i.e., the horizon) is chosen for both the pore fluid and solid phases. It is shown from the footing example that the computational periporomechanics model can replicate the shear failure modes of soils under a strip footing at the field scale.

Acknowledgment

Support for this work was provided by the US National Science Foundation under Contract numbers 1659932 and 1944009. The authors are grateful to the anonymous reviewers for their constructive reviews.

References

  • Rudnicki and Rice [1975] Rudnicki, J.W., Rice, J.. Conditions for the localization of deformation in pressure-sensitive dilatant materials. Journal of the Mechanics and Physics of Solids 1975;23(6):371–394.
  • Vardoulakis and Sulem [1995] Vardoulakis, I., Sulem, J.. Bifurcation analysis in geomechanics. Tech. Rep.; Blackie academic & professional; 1995.
  • Borja [2002] Borja, R.I.. Bifurcation of elastoplastic solids to shear band mode at finite strain. Computer Methods in Applied Mechanics and Engineering 2002;191(46):5287–5314.
  • Borja et al. [2013a] Borja, R.I., Song, X., Rechenmacher, A.L., Abedi, S., Wu, W.. Shear band in sand with spatially varying density. Journal of the Mechanics and Physics of Solids 2013a;61(1):219–234.
  • Song [2017] Song, X.. Transient bifurcation condition of partially saturated porous media at finite strain. International Journal for Numerical and Analytical Methods in Geomechanics 2017;41(1):135–156.
  • Borja and Aydin [2004] Borja, R.I., Aydin, A.. Computational modeling of deformation bands in granular media. i. geological and mathematical framework. Computer Methods in Applied Mechanics and Engineering 2004;193(27-29):2667–2698.
  • Lu and Godt [2013] Lu, N., Godt, J.W.. Hillslope hydrology and stability. Cambridge University Press; 2013.
  • Gens [2010] Gens, A.. Soil–environment interactions in geotechnical engineering. Géotechnique 2010;60(1):3–74.
  • Morse et al. [2017] Morse, M.S., Lu, N., Wayllace, A., Godt, J.W.. Evolution of strain localization in variable-width three-dimensional unsaturated laboratory-scale cut slopes. Journal of Engineering Mechanics 2017;143(9):04017085.
  • Likos et al. [2019] Likos, W.J., Song, X., Xiao, M., Cerato, A., Lu, N.. Fundamental challenges in unsaturated soil mechanics. In: Geotechnical Fundamentals for Addressing New World Challenges. Springer; 2019:209–236.
  • Alonso et al. [1990] Alonso, E.E., Gens, A., Josa, A.. A constitutive model for partially saturated soils. Géotechnique 1990;40(3):405–430.
  • Fredlund and Rahardjo [1993] Fredlund, D.G., Rahardjo, H.. Soil mechanics for unsaturated soils. John Wiley & Sons; 1993.
  • Lu and Likos [2004] Lu, N., Likos, W.J.. Unsaturated soil mechanics. Wiley; 2004.
  • Borja et al. [2013b] Borja, R.I., Song, X., Wu, W.. Critical state plasticity. part vii: Triggering a shear band in variably saturated porous media. Computer Methods in Applied Mechanics and Engineering 2013b;261:66–82.
  • Song and Borja [2014a] Song, X., Borja, R.I.. Mathematical framework for unsaturated flow in the finite deformation range. International Journal for Numerical Methods in Engineering 2014a;97(9):658–682.
  • Song and Borja [2014b] Song, X., Borja, R.I.. Finite deformation and fluid flow in unsaturated soils with random heterogeneity. Vadose Zone Journal 2014b;13(5).
  • Song et al. [2018] Song, X., Wang, K., Ye, M.. Localized failure in unsaturated soils under non-isothermal conditions. Acta Geotechnica 2018;13(1):73–85.
  • Song et al. [2017] Song, X., Ye, M., Wang, K.. Strain localization in a solid-water-air system with random heterogeneity via stabilized mixed finite elements. International Journal for Numerical Methods in Engineering 2017;.
  • Lazari et al. [2015] Lazari, M., Sanavia, L., Schrefler, B.. Local and non-local elasto-viscoplasticity in strain localization analysis of multiphase geomaterials. International Journal for Numerical and Analytical Methods in Geomechanics 2015;39(14):1570–1592.
  • Hughes [2000] Hughes, T.J.. The finite element method: linear static and dynamic finite element analysis. Courier Corporation; 2000.
  • De Borst et al. [2012] De Borst, R., Crisfield, M.A., Remmers, J.J., Verhoosel, C.V.. Nonlinear finite element analysis of solids and structures. John Wiley & Sons; 2012.
  • Loret and Prevost [1991] Loret, B., Prevost, J.. Dynamic strain localization in fluid-saturated porous media. Journal of Engineering Mechanics 1991;117(4):907–922.
  • Lu et al. [2017] Lu, M., Zhang, H., Zheng, Y., Zhang, L.. A multiscale finite element method for the localization analysis of homogeneous and heterogeneous saturated porous media with embedded strong discontinuity model. International Journal for Numerical Methods in Engineering 2017;112(10):1439–1472.
  • Lewis and Schrefler [1998] Lewis, R.W., Schrefler, B.A.. The finite element method in the static and dynamic deformation and consolidation of porous media. BOOK; John Wiley; 1998.
  • Zienkiewicz et al. [1999] Zienkiewicz, O.C., Chan, A., Pastor, M., Schrefler, B., Shiomi, T.. Computational geomechanics. Citeseer; 1999.
  • Stakgold and Holst [2011] Stakgold, I., Holst, M.J.. Green’s functions and boundary value problems; vol. 99. John Wiley & Sons; 2011.
  • De Borst [1991] De Borst, R.. Simulation of strain localization: a reappraisal of the cosserat continuum. Engineering Computations 1991;8(4):317–332.
  • De Borst et al. [1993] De Borst, R., Sluys, L., Muhlhaus, H.B., Pamin, J.. Fundamental issues in finite element analyses of localization of deformation. Engineering Computations 1993;10(2):99–121.
  • Needleman [1988] Needleman, A.. Material rate dependence and mesh sensitivity in localization problems. Computer Methods in Applied Mechanics and Engineering 1988;67(1):69–85.
  • De Borst and Duretz [2020] De Borst, R., Duretz, T.. On viscoplastic regularisation of strain-softening rocks and soils. International Journal for Numerical and Analytical Methods in Geomechanics 2020;44(6):890–903.
  • De Borst and Mühlhaus [1992] De Borst, R., Mühlhaus, H.B.. Gradient-dependent plasticity: formulation and algorithmic aspects. International Journal for Numerical Methods in Engineering 1992;35(3):521–539.
  • Pijaudier-Cabot and Bažant [1987] Pijaudier-Cabot, G., Bažant, Z.P.. Nonlocal damage theory. Journal of Engineering Mechanics 1987;113(10):1512–1533.
  • Simo et al. [1993] Simo, J.C., Oliver, J., Armero, F.. An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids. Computational Mechanics 1993;12(5):277–296.
  • Borja and Regueiro [2001] Borja, R.I., Regueiro, R.A.. Strain localization in frictional materials exhibiting displacement jumps. Computer Methods in Applied Mechanics and Engineering 2001;190(20-21):2555–2580.
  • Miehe et al. [2010] Miehe, C., Hofacker, M., Welschinger, F.. A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. Computer Methods in Applied Mechanics and Engineering 2010;199(45-48):2765–2778.
  • Moës and Belytschko [2002] Moës, N., Belytschko, T.. Extended finite element method for cohesive crack growth. Engineering fracture mechanics 2002;69(7):813–833.
  • Belytschko et al. [2013] Belytschko, T., Liu, W.K., Moran, B., Elkhodary, K.. Nonlinear finite elements for continua and structures. John wiley & sons; 2013.
  • De Borst et al. [1999] De Borst, R., Pamin, J., Geers, M.G.. On coupled gradient-dependent plasticity and damage theories with a view to localization analysis. European Journal of Mechanics-A/Solids 1999;18(6):939–962.
  • Ehlers et al. [2004] Ehlers, W., Graf, T., Ammann, M.. Deformation and localization analysis of partially saturated soil. Computer Methods in Applied Mechanics and Engineering 2004;193(27-29):2885–2910.
  • Oka et al. [2019] Oka, F., Shahbodagh, B., Kimoto, S.. A computational model for dynamic strain localization in unsaturated elasto-viscoplastic soils. International Journal for Numerical and Analytical Methods in Geomechanics 2019;43(1):138–165.
  • De Borst et al. [2006] De Borst, R., Réthoré, J., Abellan, M.A.. A numerical approach for arbitrary cracks in a fluid-saturated medium. Archive of Applied Mechanics 2006;75(10-12):595–606.
  • Callari et al. [2010] Callari, C., Armero, F., Abati, A.. Strong discontinuities in partially saturated poroplastic solids. Computer Methods in Applied Mechanics and Engineering 2010;199(23-24):1513–1535.
  • Mroginski and Etse [2014] Mroginski, J.L., Etse, G.. Discontinuous bifurcation analysis of thermodynamically consistent gradient poroplastic materials. International Journal of Solids and Structures 2014;51(9):1834–1846.
  • Mohammadnejad and Khoei [2013] Mohammadnejad, T., Khoei, A.. Hydro-mechanical modeling of cohesive crack propagation in multiphase porous media using the extended finite element method. International Journal for Numerical and Analytical Methods in Geomechanics 2013;37(10):1247–1279.
  • Ehlers and Luo [2017] Ehlers, W., Luo, C.. A phase-field approach embedded in the theory of porous media for the description of dynamic hydraulic fracturing. Computer methods in applied mechanics and engineering 2017;315:348–368.
  • Roth et al. [2020] Roth, S.N., Léger, P., Soulaïmani, A.. Fully-coupled hydro-mechanical cracking using xfem in 3d for application to complex flow in discontinuities including drainage system. Computer Methods in Applied Mechanics and Engineering 2020;370:113282.
  • Réthoré et al. [2007a] Réthoré, J., de Borst, R., Abellan, M.A.. A two-scale approach for fluid flow in fractured porous media. International Journal for Numerical Methods in Engineering 2007a;71(7):780–800.
  • Réthoré et al. [2007b] Réthoré, J., de Borst, R., Abellan, M.A.. A discrete model for the dynamic propagation of shear bands in a fluid-saturated medium. International Journal for Numerical and Analytical Methods in Geomechanics 2007b;31(2):347–370.
  • Mühlhaus and Vardoulakis [1987] Mühlhaus, H., Vardoulakis, I.. The thickness of shear bands in granular materials. Geotechnique 1987;37(3):271–283.
  • Nemat-Nasser and Hori [1999] Nemat-Nasser, S., Hori, M.. Micromechanics: overall properties of heterogeneous materials; vol. 37. Elsevier; 1999.
  • Rechenmacher and Finno [2003] Rechenmacher, A.L., Finno, R.J.. Digital image correlation to evaluate shear banding in dilative sands. Geotechnical Testing Journal 2003;27(1):13–22.
  • Alshibli and Hasan [2008] Alshibli, K., Hasan, A.. Spatial variation of void ratio and shear band thickness in sand using x-ray computed tomography. Géotechnique 2008;58(4):249–257.
  • Zhang et al. [1999] Zhang, H., Sanavia, L., Schrefler, B.. An interal length scale in dynamic strain localization of multiphase porous media. Mechanics of Cohesive-frictional Materials: An International Journal on Experiments, Modelling and Computation of Materials and Structures 1999;4(5):443–460.
  • Zhang and Schrefler [2004] Zhang, H., Schrefler, B.. Particular aspects of internal length scales in strain localization analysis of multiphase porous materials. Computer methods in applied mechanics and engineering 2004;193(27-29):2867–2884.
  • Schrefler et al. [2006] Schrefler, B., Zhang, H., Sanavia, L.. Interaction between different internal length scales in strain localization analysis of fully and partially saturated porous media—the 1-d case. International journal for numerical and analytical methods in geomechanics 2006;30(1):45–70.
  • Bažant and Jirásek [2002] Bažant, Z.P., Jirásek, M.. Nonlocal integral formulations of plasticity and damage: survey of progress. Journal of Engineering Mechanics 2002;128(11):1119–1149.
  • Dormieux et al. [2006] Dormieux, L., Kondo, D., Ulm, F.J.. Microporomechanics. John Wiley & Sons; 2006.
  • Gray and Miller [2014] Gray, W.G., Miller, C.T.. Introduction to the thermodynamically constrained averaging theory for porous medium systems. Springer; 2014.
  • Cushman [2013] Cushman, J.H.. The physics of fluids in hierarchical porous media: Angstroms to miles; vol. 10. Springer Science & Business Media; 2013.
  • Israelachvili [2011] Israelachvili, J.N.. Intermolecular and surface forces. Academic press; 2011.
  • Song and Wang [2019] Song, X., Wang, M.C.. Molecular dynamics modeling of a partially saturated clay-water system at finite temperature. International Journal for Numerical and Analytical Methods in Geomechanics 2019;43(13):2129–2146.
  • Silling [2000] Silling, S.A.. Reformulation of elasticity theory for discontinuities and long-range forces. Journal of the Mechanics and Physics of Solids 2000;48(1):175–209.
  • Silling et al. [2007] Silling, S.A., Epton, M., Weckner, O., Xu, J., Askari, E.. Peridynamic states and constitutive modeling. Journal of Elasticity 2007;88(2):151–184.
  • Jabakhanji and Mohtar [2015] Jabakhanji, R., Mohtar, R.H.. A peridynamic model of flow in porous media. Advances in Water Resources 2015;78:22–35.
  • Oterkus et al. [2017] Oterkus, S., Madenci, E., Oterkus, E.. Fully coupled poroelastic peridynamic formulation for fluid-filled fractures. Engineering Geology 2017;225:19–28.
  • Ouchi et al. [2015] Ouchi, H., Katiyar, A., York, J., Foster, J.T., Sharma, M.M.. A fully coupled porous flow and geomechanics model for fluid driven cracks: a peridynamics approach. Computational Mechanics 2015;55(3):561–576.
  • Turner [2013] Turner, D.Z.. A non-local model for fluid-structure interaction with applications in hydraulic fracturing. International Journal for Computational Methods in Engineering Science and Mechanics 2013;14(5):391–400.
  • Lai et al. [2015] Lai, X., Ren, B., Fan, H., Li, S., Wu, C., Regueiro, R.A., Liu, L.. Peridynamics simulations of geomaterial fragmentation by impulse loads. International Journal for Numerical and Analytical Methods in Geomechanics 2015;39(12):1304–1330.
  • Song and Khalili [2018] Song, X., Khalili, N.. A peridynamics model for strain localization analysis of geomaterials. International Journal for Numerical and Analytical Methods in Geomechanics 2018;.
  • Menon and Song [2019] Menon, S., Song, X.. Coupled analysis of desiccation cracking in unsaturated soils through a non-local mathematical formulation. Geosciences 2019;9(10):428.
  • Katiyar et al. [2014] Katiyar, A., Foster, J.T., Ouchi, H., Sharma, M.M.. A peridynamic formulation of pressure driven convective fluid transport in porous media. Journal of Computational Physics 2014;261:209–229.
  • Song and Silling [2020] Song, X., Silling, S.. On the peridynamic effective force state and multiphase constitutive correspondence principle. Journal of the Mechanics and Physics of Solids 2020;:104161.
  • Gropp et al. [1996] Gropp, W., Lusk, E., Doss, N., Skjellum, A.. A high-performance, portable implementation of the mpi message passing interface standard. Parallel computing 1996;22(6):789–828.
  • Gabriel et al. [2004] Gabriel, E., Fagg, G.E., Bosilca, G., Angskun, T., Dongarra, J.J., Squyres, J.M., Sahay, V., Kambadur, P., Barrett, B., Lumsdaine, A., et al. Open mpi: Goals, concept, and design of a next generation mpi implementation. In: European Parallel Virtual Machine/Message Passing Interface Users’ Group Meeting. Springer; 2004:97–104.
  • Terzaghi [1951] Terzaghi, K.. Theoretical soil mechanics. Chapman And Hall, Limited.; London; 1951.
  • Liakopoulos [1964] Liakopoulos, A.C.. Transient flow through unsaturated porous media. Ph.D. thesis; University of California, Berkeley; 1964.
  • Mandel [1953] Mandel, J.. Consolidation des sols (étude mathématique). Geotechnique 1953;3(7):287–299.
  • Nuth and Laloui [2008] Nuth, M., Laloui, L.. Effective stress concept in unsaturated soils: Clarification and validation of a unified framework. International Journal for Numerical and Analytical Methods in Geomechanics 2008;32(7):771–801.
  • Lu et al. [2010] Lu, N., Godt, J.W., Wu, D.T.. A closed-form equation for effective stress in unsaturated soil. Water Resources Research 2010;46(5).
  • Cao et al. [2018] Cao, J., Jung, J., Song, X., Bate, B.. On the soil water characteristic curves of poorly graded granular materials in aqueous polymer solutions. Acta Geotechnica 2018;13(1):103–116.
  • Van Genuchten [1980] Van Genuchten, M.T.. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal 1980;44(5):892–898.
  • Borja [2004] Borja, R.I.. Cam-clay plasticity. part v: A mathematical framework for three-phase deformation and strain localization analyses of partially saturated porous media. Computer Methods in Applied Mechanics and Engineering 2004;193(48-51):5301–5338.
  • Bear [2013] Bear, J.. Dynamics of fluids in porous media. Courier Corporation; 2013.
  • Coussy [2004] Coussy, O.. Poromechanics. John Wiley & Sons; 2004.
  • Gunzburger and Lehoucq [2010] Gunzburger, M., Lehoucq, R.B.. A nonlocal vector calculus with application to nonlocal boundary value problems. Multiscale Modeling & Simulation 2010;8(5):1581–1598.
  • Silling [2006] Silling, S.. Emu user’s manual. Sandia National Laboratories, Albuquerque, NM 2006;.
  • Silling and Askari [2005] Silling, S.A., Askari, E.. A meshfree method based on the peridynamic model of solid mechanics. Computers & structures 2005;83(17-18):1526–1535.
  • Li and Liu [2002] Li, S., Liu, W.K.. Meshfree and particle methods and their applications. Appl Mech Rev 2002;55(1):1–34.
  • Liu and Gu [2005] Liu, G.R., Gu, Y.T.. An introduction to meshfree methods and their programming. Springer Science & Business Media; 2005.
  • Simo and Hughes [1998] Simo, J.C., Hughes, T.J.. Computational inelasticity; vol. 7. Springer Science & Business Media; 1998.
  • Nova et al. [2003] Nova, R., Castellanza, R., Tamagnini, C.. A constitutive model for bonded geomaterials subject to mechanical and/or chemical degradation. International Journal for Numerical and Analytical Methods in Geomechanics 2003;27(9):705–732.
  • Borja [2013] Borja, R.I.. Plasticity: modeling & computation. Springer Science & Business Media; 2013.
  • Song and Menon [2018] Song, X., Menon, S.. Modeling of chemo-hydromechanical behavior of unsaturated porous media: a nonlocal approach based on integral equations. Acta Geotechnica 2018;:1–21.
  • Biot [1941] Biot, M.A.. General theory of three-dimensional consolidation. Journal of Applied Physics 1941;12(2):155–164.
  • Cryer [1963] Cryer, C.. A comparison of the three-dimensional consolidation theories of biot and terzaghi. The Quarterly Journal of Mechanics and Applied Mathematics 1963;16(4):401–412.
  • Verruijt [2013] Verruijt, A.. Theory and problems of poroelasticity. Delft University of Technology 2013;71.