Elastoplastic Coupled Model of Saturated Soil Consolidation under Effective Stress

Abstract

Soil’s consolidation is a geotechnical problem resulting from a stress-transfer process that initiates when the load is applied to the water contained in the soil, producing a reduction in pore water pressure and rearranging the solid particles, and thus causing a decrease in soil volume. Therefore, consolidation is a coupled flow–mechanical problem. Coupled models have been developed to simulate this phenomenon while considering different theories, providing consistent results. This paper presents an elastoplastic coupled model of consolidation under Terzaghi’s effective stress formulated using the equations of transient flow, balance moment, motion, and the critical state model that considered elastoplastic strains. The coupled model algorithm provided fast and easy results due to its flexibility, as it allowed combinations in loading and boundary conditions. Additionally, it considered the external/internal water flow as an inflow or outflow, which modified the pore water pressure and produced changes in the horizontal and vertical displacements. The numerical results obtained showed an appropriate behavior of the consolidation phenomenon, as well as the evolution of the vertical $U_y$ and horizontal $U_x$ displacements, water pressure $p_w$, volumetric ${\epsilon }_v$ and deviatoric ${\epsilon }_q$ strain, mean ${\sigma }_p$ and deviatoric ${\sigma }_q$ stress, volumetric variation $∆{\epsilon }_v$, and elastic/plastic behavior of the finite elements while considering the yield surface of the critical state.

1. Introduction

Any infrastructure development such as buildings and roads on saturated compressible cohesive soils (solid particles and water) should be avoided [1] because implications are generated for infrastructure performance and stability [2] due to their low bearing capacity [1], consolidation phenomenon [3], and settlement [4]. Consequently, the study of consolidation is one of the main problems of geotechnical engineering [5,6]. Consolidation is a decrease in soil volume that occurs over a period as a result of the water expulsion from the pores caused by the applied load [1,6,7,8,9,10,11,12,13] Volume changes result in differential settlement, leading to pavement cracking, pipe failure and foundation damage [1,14,15,16]. These damages can be more expensive than those associated with natural disasters [17]. The first theory of saturated soil consolidation was proposed by Terzaghi in 1930 [6,7,8,9,12,13,18,19] to predict the change in excess water pore pressure and settlement of the subsoil [10]. Based on this theory, the soil is saturated, Darcy’s law is applicable, the relationship between effective stress and strain is linear, and permeability (saturated hydraulic conductivity) and compressibility remain constant [9].

Moreover, most terrestrial materials exhibit behaviors inconsistent with saturated soils [4,10]; therefore, researchers extended the theory of consolidation of saturated soils to unsaturated soils [19]. These soils consist of three phases: solid, water, and air [7,20,21], and are exposed to climatic changes [22]. Their hydromechanical behavior is sensitive to moisture changes [23]; i.e., if the soil becomes wet it will swell and its strength will decrease. If it dries out, the soil will shrink and crack [23,24]. Changes in water content depend on hydraulic behavior [16], which governs infiltration, water flow (hydraulic conductivity), and soil water absorption (retention curve) [25]. Soil water infiltration helps replenish groundwater and depends on topography, density, texture, and initial water content of the soil, and serves to measure the soil water retention curve (SWRC) [2,26]; while the hydraulic conductivity represents the ease with which water can pass through a soil and depends on the fluid properties (viscosity and chemistry), soil pore size, and temperature [27,28]. In addition, the SWRC relates the degree of saturation to soil suction [29]; this has been used to determine the water retention of an unsaturated soil [30]. These are not unique to a soil; its shape and position are not the same and its wetting and drying paths are different, showing the hysteresis phenomenon, which means that for one value of suction, there are two values of degree of saturation [29,31]. For this reason, the hysteresis of the SWRC should be considered to estimate the behavior of unsaturated soils [32].

On the other hand, the consolidation of saturated soils is governed by the effective stress, as it relates to the soil’s stiffness property [33] The effective stress combines the applied stress and water pore pressure, enabling the conversion of a multiphase porous medium to an equivalent single-phase continuous one [34,35]. Likewise, [36] established that effective stress can also describe the consolidation and SWRC behavior of unsaturated soils. Therefore, the consolidation problem has been formulated by coupling soil strain and pore water pressure [6]. Due to this, two approaches to consolidation are known: (a) the uncoupled theory, which cannot model all aspects of saturated soils consolidation but provides sufficiently accurate solutions; and (b) the coupled theory, which provides a coupling between the magnitude and progress of settlement, resulting in more accurate solutions [4,7,19]. Similarly, there are two methodologies to investigate the effect of hydromechanical coupling for the consolidation of unsaturated soils: (1) uncoupled modeling using independent stress state variables; and (2) coupled modeling of air and water flow and strain [9]. Based on what has been described, the use of effective stress is simple, requires fewer parameters to describe the behavior of unsaturated soil, and allows for natural hydromechanical coupling [37]. Thus, a proper analysis requires the coupling of flow with consolidation and the inclusion of SWRC hysteresis effects [2,38] and a viscous component [39].

This motivated several researchers to analyze the problem of saturated consolidation and how it was affected by different variables. These works were listed in Dong et al. [9], in which it was observed that Fox and Berles [40] and Gibson et al. [41] analyzed large strains; stratified soils were studied by Fox Patrick et al. [42] and accretion layers by Fox [43]; vertical and radial flows were applied in Fox Patrick et al. [44]; a constant strain rate was considered by Pu et al. [45] and Wissa Anwar et al. [46]; secondary compression effects appeared in Brandenberg Scott [47]; solute transport was discussed by Fox Patrick [48]; and variation in hydraulic conductivity and compressibility was taken into account by Indraratna et al. [49]. Self-weight was included in Qi et al. [2], Wang et al. [18], and Wang, et al. [50]. The study of unsaturated soils has progressed from the theory of Blight [51] and Liakopoulos [52], but Fredlund and Hasan [53] considered the nonlinear behavior of air and water flow. This theory was extended in 2D models by Dakshanamurthy and Fredlund [54] and Darkshanamurthy et al. [8] and 3D models by Zhou and Zhao [4] and Ho et al. [11]. In addition, one-dimensional unsaturated consolidation was studied by Zhou and Zhao [4], Shan et al. [12], Zhou et al. [19], and Qin et al. [21], while axisymmetric unsaturated consolidation was studied by Ho et al. [11].

During the last three decades, constitutive models have been developed based on elastoplastic theory, as it allows the modeling of the nonlinear and hysteretic behavior of soil [55,56]. These aimed to carry out hydromechanical coupling, such as in the Basic Barcelona Model (BBM) and the Drucker–Prager Model (DPM), which were defined with respect to mean and deviatoric stress [57]; while Arroyo and Rojas [34] developed a fully coupled hydromechanical model under effective stresses that accurately predicted the transition between saturated and unsaturated states as well as between elastic and elastoplastic states, but with limitations in the prediction of volumetric strains. Due to this, there are no analytical solutions for hydromechanical coupling in unsaturated soils; therefore, it is necessary to employ numerical techniques such as the finite element (FE) method [6,7,13,38] in order to simplify the nonlinear behavior of the unsaturated consolidation theory, assuming that water permeability and air transmittance are constant [12]. As a result, several coupled models have been developed to simulate saturated soils consolidation [58,59,60,61]. These models were based on: (a) the soil being saturated; (b) water phase being incompressible; (c) solid phase being incompressible but with a compressible arrangement; and (d) Darcy’s law governing the behavior of the flow through the soil. In addition, some authors contributed improvements by incorporating certain theories: the theory of small strains was used by Bentler [58]; the principle of virtual work was used by Manzolillo et al. [61]; and the mass conservation law was considered in the models in [59] and in Manzolillo et al. [61]. The last two employed Terzaghi’s effective stress and Galerkin’s solution method of weighted residuals, respectively. Additionally, Krishnamoorthy [60] considered the nonlinear behavior of soil and used the hyperbolic relationship of Duncan and Chang [62]. All these models reported consistent results with those obtained in the domain and in the laboratory while showing advantages and disadvantages, suggesting that they can be improved.

The literature shows that consolidation of compressible cohesive soils should be analyzed by means of a fully coupled numerical FE model under the concept of effective stresses able to generate the transition between the saturated and unsaturated states, as well as between the elastic and elastoplastic states, allowing researchers to obtain an adequate consolidation behavior under different boundaries conditions, applied stresses, and water flow into the soil. The above is done to reduce the time and cost of this type of analysis.

The objective of this paper was to formulate and develop the elastoplastic coupled model (flow–mechanical–critical state) for the saturated consolidation of soils under effective stresses while considering the “Moment Balance Laws”, “transient flow” (flow) and “Motion” (mechanical) equations in which Terzaghi’s “effective stress” was introduced. The inner product between vector functions was used while considering the “Principle of Virtual Work” associated with the “Variational Principle of Minimum Potential Energy” in combination with the finite element method (FEM) and Galerkin’s method, which includes a time step that allows the evolution of displacements and water pore pressure as an approximate solution. Additionally, the model was formulated to extend it to unsaturated consolidation, as it permitted the introduction of Bishop’s effective stress, which included the Terzaghi effective stress (σ-u_a) and the matrix stress (χψ) representing the product of the chi parameter and the soil suction this last parameter was obtained from the SWRC. The variation of permeability with respect to soil suction (hydraulic conductivity) could also be obtained from the SWRC.

2. Governing Equations

The applied stress to a saturated soil will cause it to undergo a reduction in soil volume resulting from the expulsion of water from the pores in a time-dependent manner, generating the transient process called consolidation. However, much of the planet is covered by three-phase (unsaturated) soils: solid, liquid, and air. These soils experience variations in their water content produced by climatic conditions that subject them to wetting and drying cycles. Therefore, a numerical procedure is required to obtain the solution to a coupled consolidation problem (saturated–unsaturated) using the finite element method. Thus, the proposed coupled elastoplastic model was formulated and prepared considering the consolidation of saturated soil, but it could be modified for the consolidation of unsaturated soil after coupling the variation in hydraulic conductivity with respect to the degree of saturation obtained from the soil water retention curve (SWRC). Hence, soil mechanical deformation and the flowing water pressure constitute a coupled field problem. Considering the above, the following assumptions were considered: (a) complete saturation and incompressibility of water; (b) incompressible solid phase (soil skeleton); and (c) Darcy’s law governed the flow behavior through the soil. Additionally, the theory of small deformations was used. According to continuum mechanics theory, the governing equations are provided below:

2.1. Equilibrium Equations

According to the formulation of the elastoplastic model shown in Appendix A of this paper and the balance laws of continuum mechanics and considering the Piola–Kirchhoff stress tensor $\mathit{S}$ in the Lagrangian configuration, this tensor provides the measured force per unit area in the reference configuration [63]:

$$div\mathit{S}+{\mathit{b}}_0=\ddot{\mathsf{\chi }}{\rho }_0,$$

(1)

In this study, $div$ is the divergence operator that measures the difference between the outgoing and incoming flux of a vector field on the surface surrounding a control volume. The body forces ${\mathit{b}}_0=\mathbf{0}$ and acceleration field $\ddot{\mathsf{\chi }}=\mathbf{0}$ were considered because the consolidation process could be regarded as a quasistatic process.

2.2. Continuity Equations

$$div\mathit{v}-\dot{\epsilon }=0,$$

(2)

In Equation (2), $\mathit{v}$ is the field velocity and $\dot{\epsilon }$ is the strain derivative of the field. Thus, Equations (1) and (2) are the basic phenomenological equations of the coupled problem (see Appendix A). This coupling procedure depends on the application of the principle of virtual work, which indicates that the work developed by the internal forces within a system equals the work performed by the external forces acting on the system [64]. This principle is associated with the variational principle of the minimum potential energy [65].

2.3. Darcy’s Law

$$v_i=k_{ij}\left({\gamma }_w{h}_j\right)/{\gamma }_w,$$

(3)

where $v_i$ is the flow velocity field, $k_{ij}$ is the soil permeability coefficient, ${\gamma }_w$ is the volumetric mass of water, and $h_i$ is the piezometric level, which can be evaluated according to the reference level z; i.e., $h_i=p_w/{\gamma }_w+z$.

2.4. Boundary Conditions

The boundary conditions are related to the pore water pressure at all boundaries and to the applied stress load at the surface. Two additional conditions were added: surface water flow produced by meteorological conditions and water flow into the soil due to pipe failure. All boundaries $\mathsf{\Omega }$ evolve over time and include Dirichlet ${\mathsf{\Omega }}_D$ and Newman ${\mathsf{\Omega }}_N$ boundaries. Hence, the boundary conditions can be written in terms of prescribed values:

(a)

Pore pressure conditions:

$$\mathit{u}=\widehat{\mathit{u}}\mathrm{for}\partial {\mathsf{\Omega }}_D,$$

(4)

(b)

Applied stress to the surface:

$$\mathit{S}\mathit{n}=\mathit{f}\mathrm{for}\partial {\mathsf{\Omega }}_N,$$

(5)

(c)

Superficial water flow:

$$\mathit{Q}=\widehat{\mathit{Q}}\mathrm{for}\partial {\mathsf{\Omega }}_f,$$

(6)

(d)

Water flow into the soil:

$${\mathit{Q}}_{\mathit{s}}=\widehat{{\mathit{Q}}_{\mathit{s}}}\mathrm{for}\partial {\mathsf{\Omega }}_f,$$

(7)

(e)

Displacements on the soil boundary:

$$\mathit{U}=\widehat{{\mathit{U}}_{\mathit{s}}}\mathrm{for}\partial {\mathsf{\Omega }}_u,$$

(8)

3. Constitutive Equations

3.1. Phenomenological Case of the Solid in the Soil

The inner product was set in Equation (1) while integrating the volume Ω of the body under study. Note that $div\mathit{S}\in V$ (vector space $V$) and such inner product associates a work. For simplicity of writing, ϕ^T is indicated at this point as φ, showing that it is a vector function and the product $\phi div\mathit{S}$ is a scalar function; rewriting in alternate notation gives:

$${{\displaystyle \int }}_{V}div\mathit{S}\phi dV={{\displaystyle \int }}_V{\varphi }^{T}divSdV={{\displaystyle \int }}_V\phi div\mathit{S}dV={{\displaystyle \int }}_V{\phi }_i{\sigma }_{ij,j}dV=0,$$

(9)

Since ${\sigma }_{ij,j}$ is equivalent to $\mathit{S}$, ${\sigma }_{ij,j}=\partial {\sigma }_{ij}/x_j$ is the divergence of the tensor $\mathit{S}$. Moreover, using the derivative in $j$ of the expression ${\phi }_i{\sigma }_{ij,j}$ according to the derivative product rule, Equation (9) becomes (10):

$${{\displaystyle \int }}_V{\left({\phi }_i{\sigma }_{ij}\right)}_{j}dV-{{\displaystyle \int }}_{\mathit{V}}{\phi }_{i,j}{\sigma }_{ij}dV=0,$$

(10)

Applying the “Divergence Theorem” [63] for an arbitrary scalar field φ with Γ is given as:

$${{\displaystyle \int }}_V{\left({\phi }_i{\sigma }_{ij}\right)}_{j}dV={{\displaystyle \int }}_{\mathsf{\Gamma }}{\left({\phi }_i{\sigma }_{ij}\right)}_j\mathit{n}dA={{\displaystyle \int }}_{\mathsf{\Gamma }}{\phi }_i{\sigma }_{ij}{n}_{j}dA\mathbf{,}$$

(11)

Substituting (11) into (10) gives:

$${{\displaystyle \int }}_V{\phi }_{i,j}{\sigma }_{ij}dV={{\displaystyle \int }}_{\mathsf{\Gamma }}{\phi }_i{\sigma }_{ij}{n}_{j}dA,$$

(12)

The stress tensor $\mathit{S}$ considers Terzaghi’s effective stress equation involving stresses in both the soil skeleton and water contained in soil pores; hence, in index notation it is given as in Equation (13):

$$\mathit{S}=\mathit{S}\prime +p_w\mathit{I},$$

(13)

where $\mathit{S}$ is the soil total stress tensor, $\mathit{S}\prime$ is the soil effective stress, $p_w$ is the pore water pressure (hydrostatic stresses), and $\mathit{I}$ is the identity tensor $={\displaystyle \sum }_i{e}_i\otimes e_i$.

$${\sigma }_{ij}=\sigma {\prime }_{ij}+p_w{\delta }_{ij},$$

(14)

Substituting (14) into (12) results in:

$${{\displaystyle \int }}_V{\phi }_{i,j}\sigma {\prime }_{ij}dV+{{\displaystyle \int }}_V{\phi }_{i,j}{p}_w{\delta }_{ij}dV={{\displaystyle \int }}_{\mathsf{\Gamma }}{\phi }_i{\sigma }_{ij}{n}_{j}dA,$$

(15)

On the other hand, ${\phi }_{i,j}$ associates the derivative of a vector function corresponding in turn to a tensor, whereas any tensor $\mathit{A}=\mathit{E}+\mathit{W}$ can be represented as the sum of a symmetric tensor plus an antisymmetric tensor [63], which in index notation are $\mathit{E}={\phi }_{\left(i,j\right)}$ and $\mathit{W}={\phi }_{\left[i,j\right]}$; therefore, ${\phi }_{i,j}={\phi }_{\left(i,j\right)}+{\phi }_{\left[i,j\right]}$ when substituting into the first integral of (15) and recalling that the inner product $\mathit{W}·\mathit{S}=0$, since $\mathit{W}\in Asim$.

$${{\displaystyle \int }}_V{\phi }_i\sigma {\prime }_{ij}dV={{\displaystyle \int }}_V{\phi }_{\left(i,j\right)}\sigma {\prime }_{ij}dV,$$

(16)

The term ${\phi }_{\left(i,j\right)}$ relates the strain tensor in the index notation ${\epsilon }_{ij}^{\phi }={\phi }_{\left(i,j\right)}$.

The field equations of elastostatics use the constitutive equation $\mathit{S}=\mathit{C}\left[\mathit{E}\right]$, which relates the stress tensor to the infinitesimal strain tensor, which in index notation is $\sigma {\prime }_{ij}=\mathit{C}{\epsilon }_{ij}$ (derived from Hooke’s law). Thus, considering the strain tensor and Hooke’s law in (16) obtains:

$${{\displaystyle \int }}_V{\phi }_i\sigma {\prime }_{ij}dV={{\displaystyle \int }}_V{\phi }_{\left(i,j\right)}\mathit{D}{\prime }_{ijkl}{\epsilon }_{kl}^{u}dV,$$

(17)

By substituting (17) into (15), it is possible to obtain:

$${{\displaystyle \int }}_V{\phi }_{\left(i,j\right)}\mathit{D}{\prime }_{ijkl}{\epsilon }_{kl}^{u}dV+{{\displaystyle \int }}_V{\phi }_{i,j}{p}_w{\delta }_{ij}dV={{\displaystyle \int }}_{\mathsf{\Gamma }}{\phi }_i{\sigma }_{ij}{n}_{j}dA\mathbf{,}$$

(18)

According to the strain tensor, ${\epsilon }_{ij}^{\phi }={\phi }_{\left(i,j\right)}$ and it applies for $p_w{\delta }_{ij}$, as ${\delta }_{ij}=0$ for $i\ne j$. In addition, recalling that ${\sigma }_{ij}{n}_j=\mathit{S}\mathit{n}$ in accordance with the “Cauchy hypothesis” of stress existence [63], $\mathit{S}\mathit{n}=\mathit{s}\left(\mathit{n}\right)=\mathit{t}\equiv$ vector of surface tractions, so Equation (18) takes the form:

$${{\displaystyle \int }}_V{\phi }_{\left(i,j\right)}\mathit{D}{\prime }_{ijkl}{\epsilon }_{kl}^{u}dV+{{\displaystyle \int }}_V{\epsilon }_{ij}^{\phi }{p}_{w}dV={{\displaystyle \int }}_{\mathsf{\Gamma }}{\phi }_i{t}_{j}dA\mathbf{,}$$

(19)

The finite element approximation (Galerkin) requires:

$$\begin{gathered} \phi =\left[N_{\phi }\right]\left\{\varphi \right\},\hspace{1em}\hspace{1em}\hspace{1em}u=\left[N_u\right]\left\{U\right\}, \\ u=\left[N_u\right]\left\{U\right\},\hspace{1em}\hspace{1em}\hspace{1em}{\epsilon }_{kl}^u=\left[B_u\right]\left\{U\right\}, \\ {\epsilon }_{ij}^{\phi }=\left[B_{\phi }\right]\left\{\varphi \right\},\hspace{1em}\hspace{1em}\hspace{1em}{p}_w=\left[N_p\right]\left\{p_w\right\}, \end{gathered}$$

(20)

where $\left[N\right]$ is the matrix of shape functions and $\left[B\right]$ is the matrix of derivatives of shape functions.

$${({\epsilon }_{ij}^{\phi })}^T={\left\{\varphi \right\}}^T{\left[B_{\phi }\right]}^T,\hspace{1em}\left\{{\epsilon }_{ii}^{\phi }\right\}=\left[B_{\phi i}\right]\left\{\varphi \right\},\hspace{1em}{\left\{{\epsilon }_{ii}^{\phi }\right\}}^T={\left\{\varphi \right\}}^T{\left[B_{\phi i}\right]}^T,$$

(21)

Equations (20) and (21) allow their variations to be set in the following form (where ${\delta }_{ij}$ is the “Kronecker Delta”):

$$\begin{gathered} \left\{\delta \phi \right\}=\left[N_{\phi }\right]\left\{\delta \varphi \right\},\hspace{1em}\hspace{1em}\hspace{1em}{\left\{\delta \phi \right\}}^T={\left\{\delta \varphi \right\}}^T{\left[N_{\phi }\right]}^T, \\ \left\{\delta {\epsilon }_{ij}^{\phi }\right\}=\left[B_{\phi }\right]\left\{\delta \varphi \right\},\hspace{1em}\hspace{1em}\hspace{1em}{\left\{\delta {\epsilon }_{ij}^{\phi }\right\}}^T={\left\{\delta \varphi \right\}}^T{\left[B_{\phi }\right]}^T, \\ \left\{\delta {\epsilon }_{ii}^{\phi }\right\}=\left[B_{\phi i}\right]\left\{\delta \varphi \right\},\hspace{1em}\hspace{1em}\hspace{1em}{\left\{\delta {\epsilon }_{ii}^{\phi }\right\}}^T={\left\{\delta \varphi \right\}}^T{\left[B_{\phi i}\right]}^T, \end{gathered}$$

(22)

Applying the “Virtual Work Theorem” [63] and substituting the variations in the fields $\phi$ and $\epsilon$ in Equation (19) gives:

$${{\displaystyle \int }}_V{(\delta {\epsilon }_{ij}^{\phi })}^T\mathit{D}\prime \left[B_u\right]\left\{U\right\}dV+{{\displaystyle \int }}_V{\left(\delta {\epsilon }_{ii}^{\phi }\right)}^T\left[N_p\right]\left\{p_w\right\}dV={{\displaystyle \int }}_{\mathsf{\Gamma }}{\left(\delta \phi \right)}^T{\mathit{t}}_{j}dA\mathbf{,}$$

(23)

Equation (23) establishes the expression of the virtual work in which a small variation in the function $\phi$ was applied to generate such work. When substituting the variations according to Equations (22) and the unitary isotropic tensor $m$ [64], which correlates the degrees of freedom of the adjoint functions ${\left[B_{\phi }\right]}^T$ and $\left[N_p\right]$, it can be written as:

$${{\displaystyle \int }}_V{\left\{\delta \varphi \right\}}^T{\left[B_{\phi }\right]}^T\left[D\right]\left[B_u\right]\left\{U\right\}dV+{{\displaystyle \int }}_V{\left\{\delta \varphi \right\}}^T{\left[B_{\phi }\right]}^T\langle m\rangle \left[N_p\right]\left\{p_w\right\}dV={{\displaystyle \int }}_{\mathsf{\Gamma }}{\left\{\delta \varphi \right\}}^T{\left[N_{\phi }\right]}^T{\mathit{t}}_{j}dA\mathbf{,}$$

(24)

because ${\left\{\delta \varphi \right\}}^T$ does not involve integration variables. Since it is a common term on both sides of the equality, Equation (24) becomes:

$${{\displaystyle \int }}_V{\left[B_{\phi }\right]}^T\left[D\right]\left[B_u\right]dV\left\{U\right\}+{{\displaystyle \int }}_V{\left[B_{\phi }\right]}^T{\begin{array}{ccc}\langle 1& 1& 0\rangle \end{array}}^T\left[N_p\right]dV\left\{p_w\right\}={{\displaystyle \int }}_{\mathsf{\Gamma }}{\left[N_{\phi }\right]}^T{\mathit{t}}_{j}dA\mathbf{,}$$

(25)

The identification of the matrices with the subscripts $\phi$, $u$, and $p$ allowed us to recognize where their approximation originated.

3.2. Phenomenological Case of the Flow in the Soil

By integrating the phenomenological Equation (2) in the volume of the body under study and proceeding analogously to the case of soil solids, applying an inner product with an arbitrary function $\theta$ to generate a work gives:

$${{\displaystyle \int }}_{V}divv\theta dV-{{\displaystyle \int }}_V\dot{\epsilon }\theta dV={{\displaystyle \int }}_V{\theta }^{T}divvdV-{{\displaystyle \int }}_V{\theta }^T\dot{\epsilon }dV=0,$$

(26)

Denoting for the moment, as was done in (9) $\vartheta ={\theta }^T$, and reminding that $divv$ in index notation is $divv=v_{i,i}$, then (26) transforms into:

$${{\displaystyle \int }}_V\vartheta v_{i,i}dV-{{\displaystyle \int }}_V\vartheta {\dot{\epsilon }}_{ii}dV=0,$$

(27)

Considering and applying the derivative product rule in (27) gives:

$${{\displaystyle \int }}_V{\left(\vartheta v_i\right)}_{,i}dV-{{\displaystyle \int }}_V{\vartheta }_{,i}{v}_{i}dV-{{\displaystyle \int }}_V\vartheta {\dot{\epsilon }}_{ii}dV=0,$$

(28)

According to the “Divergence Theorem” [63] in the first term of (28), where $v_i{n}_i=q_i$, and considering Equation (3), which establishes Darcy’s law, results in:

$${{\displaystyle \int }}_V{\vartheta }_{,j}\left[k_{ij}{\left({\gamma }_{w}h\right)}_{,j}/{\gamma }_w\right]dV+{{\displaystyle \int }}_V\vartheta {\dot{\epsilon }}_{ii}dV={{\displaystyle \int }}_{\Gamma }\vartheta {\mathit{q}}_{i}dA,$$

(29)

Here, ${\epsilon }_{ii}=u_{i,i};$ therefore, $\dot{{\epsilon }_{ii}}=\dot{u_{i,i}}$. Furthermore, the piezometric level $h=\left(p_w/{\gamma }_w\right)+z$ can be shown in terms of the water pore pressure $p_w$ and the reference head $z$:

$${{\displaystyle \int }}_V{\vartheta }_{,j}\left(k_{ij}/{\gamma }_w\right)\left\{p_w\right\}dV+{{\displaystyle \int }}_V{\vartheta }_{,j}{k}_{ij}{z}_{,j}dV+{{\displaystyle \int }}_V\vartheta {\dot{u}}_{i,i}dV={{\displaystyle \int }}_{\Gamma }\vartheta {\mathit{q}}_{i}dA,$$

(30)

The finite element approximation (Galerkin) requires:

$$\begin{gathered} \vartheta =\left[N_{\theta }\right]\left\{\theta \right\},\hspace{1em}\hspace{1em}\hspace{1em}{\vartheta }^T={\left\{\theta \right\}}^T{\left[N_{\theta }\right]}^T, \\ {\vartheta }_{,j}=\left[B_{\theta }\right]\left\{\theta \right\},\hspace{1em}\hspace{1em}\hspace{1em}{\vartheta }_{,j}^T={\left\{\theta \right\}}^T{\left[B_{\theta }\right]}^T, \\ \left\{p_w\right\}=\left[N_p\right]\left\{p\right\},\hspace{1em}\hspace{1em}\hspace{1em}{\left\{p_w\right\}}_{,j}=\left[B_p\right]\left\{p\right\}, \end{gathered}$$

(31)

$$u=\left[N\right]\left\{U\right\},\hspace{1em}\hspace{1em}\hspace{1em}{\epsilon }_{ii}=u=\left[B_{ui}\right]\left\{U\right\}$$

(32)

When deriving Equation (32), the result is:

$$\left\{\dot{{\epsilon }_{ii}}\right\}=\left[B_{ui}\right]\{\dot{U}\},$$

(33)

indicating that the associated variations in (34) apply the “Principle of Virtual Work” in Equation (30):

$$\begin{gathered} \left\{\delta \vartheta \right\}=\left[N_{\theta }\right]\left\{\delta \theta \right\},\hspace{1em}\hspace{1em}\hspace{1em}{\left\{\delta \vartheta \right\}}^T={\left\{\delta \theta \right\}}^T{\left[N_{\theta }\right]}^T, \\ \left\{\delta {\vartheta }_{,j}\right\}=\left[B_{\theta }\right]\left\{\delta \theta \right\},\hspace{1em}\hspace{1em}\hspace{1em}{\left\{\delta {\vartheta }_{,j}\right\}}^T={\left\{\delta \theta \right\}}^T{\left[B_{\theta }\right]}^T, \end{gathered}$$

(34)

When substituting (34) into (30), the unit isotropic tensor $m$ gives:

$$\begin{gathered} {{\displaystyle \int }}_V{\left\{\delta \vartheta \right\}}^T{\left[B_{\theta }\right]}^T\left(k_{ij}/{\gamma }_w\right)\left[B_p\right]\left\{p_w\right\}dV+{{\displaystyle \int }}_V{\left\{\delta \theta \right\}}^T{\left[N_{\theta }\right]}^{T}m\left[B_{ui}\right]\{\dot{U}\}dV= \\ {{\displaystyle \int }}_{\Gamma }{\left\{\delta \theta \right\}}^T{\left[N_{\theta }\right]}^T{\mathit{q}}_{i}dA, \end{gathered}$$

(35)

In Equation (35), the term ${{\displaystyle \int }}_V{\vartheta }_{,j}{k}_{ij}{z}_{,j}dV$ was not considered because it depends on the reference level at which the piezometric level is evaluated, which for the problem reviewed was not included. On the other hand, substituting the matrix $\left[B_{ui}\right]$ for the product of the unit isotropic tensor in the coupling matrix and extracting the term ${\left\{\delta \vartheta \right\}}^T$ from the integrals obtains:

$${{\displaystyle \int }}_V{\left[B_{\theta }\right]}^T\left(k_{ij}/{\gamma }_w\right)\left[B_p\right]dV\left\{p_w\right\}+{{\displaystyle \int }}_V{\left[N_{\theta }\right]}^T{\begin{array}{ccc}1& 1& 0\end{array}}^T\left[B_{ui}\right]dV\{\dot{U}\}={{\displaystyle \int }}_{\Gamma }{\left[N_{\theta }\right]}^{T}dA\left\{q_i\right\},$$

(36)

Here, $\left[B_{\theta }\right]$ has a correspondence analogous to $\left[B_p\right]$. In addition, $k_{ij}$ is associated with the permeability matrix $\left[k\right]$.

The final discrete coupled Equations (25) and (36) can thus be written as follows:

$$\left[K\right]\left\{U\right\}+{\left[K_V\right]}^T\left\{p_w\right\}=\left\{F\right\},$$

(37)

$$\left[K_V\right]\left\{\dot{U}\right\}+\left(1/{\gamma }_w\right)\left[K_h\right]\left\{p_w\right\}=\left\{Q\right\},$$

(38)

where $\left[K\right]$ is the soil stiffness matrix, ${\left[K_V\right]}^T$ is the transpose coupling matrix, $\left\{U\right\}$ is the nodal displacement vector, $\left\{p_w\right\}$ is the nodal pore water pressure vector, $\left\{F\right\}$ is the external force vector, $\left[K_V\right]$ is the coupling matrix, $\left[K_h\right]$ is the drainage matrix, $\left\{Q\right\}$ is the external/internal water flow vector, and $\{\dot{U}\}$ is the nodal displacement velocity vector.

The resulting matrices in Equations (37) and (38) suggest the following integral forms (Equation (39)):

$$\begin{gathered} \left[K\right]={{\displaystyle \int }}_V{\left[B_{\phi }\right]}^T\left[D\right]\left[B_u\right]dV, \\ {\left[K_V\right]}^T={{\displaystyle \int }}_V{\left[B_{\phi }\right]}^{T}1 1 0^T\left[N_p\right]dV, \\ \left[K_h\right]={{\displaystyle \int }}_V{\left[B_{\theta }\right]}^T\left[k\right]\left[B_p\right]dV, \end{gathered}$$

(39)

Matrices $\left[B_i\right]$ contain the derivatives of the shape functions; $\left[k\right]$ is the permeability matrix, and a vector is needed to couple the values corresponding to the degrees of freedom of the displacement $\left\{U\right\}$ and pore water pressure $\left\{p_w\right\}$.

The coupled model provides a solution for one-dimensional consolidation of saturated soils, which were developed and expressed as Equations (37) and (38), and can be formulated based on the concept of effective stresses. Therefore, coupling the critical state model (Cam-clay) to the above model also yields the mean effective ${\sigma }_p$ and deviatoric ${\sigma }_q$ stresses in the elastoplastic behavior of saturated soils [66].

The strain in the critical state model can be evaluated while considering the yield surface $f$ as follows:

$$f={\sigma }_q^2-M^2{\sigma }_p\left({\sigma }_0-{\sigma }_p\right),$$

(40)

$$M=\frac{6\sin \phi }{3-\sin \phi },$$

(41)

where $M$ is the slope of the critical state line or stress path in the CD triaxial test method, $\phi$ is the soil internal friction angle, and ${\sigma }_0$ is the preconsolidation stress, which can be considered a material hardening parameter. In the event of plastic strains, the stress can evolve the hardening law in Equation (42):

$${\left(∆{\sigma }_0\right)}_i=\frac{\upsilon {\sigma }_0}{\left(\lambda -\kappa \right)}{\left(∆{\epsilon }_v^P\right)}_i,$$

(42)

where $∆{\epsilon }_v^P$ is the variation in the plastic volumetric strain and $\lambda$ and $\kappa$ are the slopes of the virgin and unloading sections, respectively, of the compressibility curve of the saturated consolidation test method.

According to Castro Barco [67], the yield surface ($f=0$) represents an ellipse in the ${\sigma }_p$-${\sigma }_q$ (${\sigma }_q$ is the deviatoric stress). The yield surface is the limit of the elastic stress states. These stresses within this limit only give rise to elastic strain increments; therefore, in this case the stiffness matrix is the stiffness matrix of the elastic component $\left[K\right]=\left[K_e\right]$. Conversely, while stresses that tend to exceed the limit give rise to increases in plastic strains, in this case, both elastic and plastic strains occur, where the stiffness matrix is a composite of elastic and plastic components $\left[K\right]=\left[K_e\right]+\left[K_p\right]$. Therefore, the terms of Equations (37) and (38) are maintained and only the stiffness matrix $\left[K\right]$ is modified. Then, Equation (37) can be rewritten as Equation (43):

$$\left(\left[K_e\right]+\left[K_p\right]\right)\left\{U\right\}+{\left[K_V\right]}^T\left\{p_w\right\}=\left\{F\right\},$$

(43)

The coupled model established in Equations (37) and (38) was formulated in terms of Cartesian stresses $\left({\sigma }_x,{\sigma }_y,\mathrm{and}{\sigma }_z\right)$ and strains $\left({\epsilon }_x,{\epsilon }_y,\mathrm{and}{\epsilon }_z\right)$, while the critical state provides constitutive equations formulated while considering both volumetric ${\epsilon }_v$ and deviatoric ${\epsilon }_q$ strains. Due to the above, the constitutive Equation (43) must be transformed as a function of both the stresses $\left({\sigma }_x,{\sigma }_y,\mathrm{and}{\tau }_{xy}\right)$ and strains $\left({\epsilon }_x,{\epsilon }_y,\mathrm{and}{\gamma }_{xy}\right)$. Hence, since this case involves plane strain, then ${\epsilon }_z=0$ even though ${\sigma }_z$ exists.

$$\left\{\begin{array}{c}{\epsilon }_v\\ {\epsilon }_q\end{array}\right\}=\left[\begin{array}{cc}\frac{1}{K}& 0\\ 0& \frac{1}{3G}\end{array}\right]\left\{\begin{array}{c}{\sigma }_p\\ {\sigma }_q\end{array}\right\}\to \left\{\begin{array}{c}{\sigma }_p\\ {\sigma }_q\end{array}\right\}={\left[\begin{array}{cc}\frac{1}{K}& 0\\ 0& \frac{1}{3G}\end{array}\right]}^{-1}\left\{\begin{array}{c}{\epsilon }_v\\ {\epsilon }_q\end{array}\right\},$$

(44)

$${\sigma }_p=\frac{{\sigma }_x+{\sigma }_y+{\sigma }_z}{3},\hspace{1em}{\sigma }_z=\frac{\upsilon E}{\left(1+\upsilon \right)\left(1-2\upsilon \right)}\left({\epsilon }_x+{\epsilon }_y\right),\hspace{1em}{\sigma }_q={\sigma }_1-{\sigma }_2,$$

(45)

$$K=\frac{E}{3\left(1-2\upsilon \right)},G=\frac{E}{2\left(1+\upsilon \right)}=\frac{{\tau }_{xy}}{{\gamma }_{xy}}=\frac{{\tau }_{max}}{{\gamma }_{max}},E=\frac{3\left(1+e_0\right)\left(1-2\upsilon \right){\sigma }_{apl}}{\kappa e_f},$$

(46)

where $\upsilon$ and $E$ are Poisson’s ratio and the elastic modulus, respectively, of the soil; ${\sigma }_1,{\sigma }_2,\mathrm{and}{\sigma }_3$ are the principal stresses applied in the CD triaxial test process corresponding to the stresses ${\sigma }_x,{\sigma }_y,\mathrm{and}{\sigma }_z$, respectively; ${\epsilon }_x,{\epsilon }_y,\mathrm{and}{\epsilon }_z$ are the strains along the corresponding directions; and $K$ and $G$ are the volumetric (bulk) and shear moduli, respectively. The last one relates shear stress and angular strain, $e_0$ and $e_f$ are the initial and final void ratios, respectively, and ${\sigma }_{apl}$ is the applied surface stress.

Similarly, the stiffness matrix of the plastic component originates from the determination of the mean effective $\sigma {\prime }_p$ and deviatoric $\sigma {\prime }_q$ stresses as follows:

$$\left\{\begin{array}{c}d{\sigma }_p\\ d{\sigma }_q\end{array}\right\}=-\frac{1}{M^2{\sigma }_p\frac{v{\sigma }_0}{\left(\lambda -\kappa \right)}}{\left[\begin{array}{cc}{M}^2\left(2{\sigma }_p-{\sigma }_0\right)& 2{\sigma }_q\\ 2{\sigma }_q& \frac{4{\sigma }_q^2}{M^2\left(2{\sigma }_p-{\sigma }_0\right)}\end{array}\right]}^{-1}\left\{\begin{array}{c}d{\epsilon }_v^P\\ d{\epsilon }_q^P\end{array}\right\},$$

(47)

Finally, the constitutive equations of the elastoplastic coupled model of saturated soil consolidation under effective stresses can be obtained as:

$$\left(\left[K_e\right]+\left[K_p\right]\right)\left\{U\right\}+{\left[K_V\right]}^T\left\{p_w\right\}=\left\{F\right\},$$

(48)

$$\left[K_V\right]\left\{\dot{U}\right\}+\left(1/{\gamma }_w\right)\left[K_h\right]\left\{p_w\right\}=\left\{Q\right\},$$

(49)

The resulting matrices in Equations (48) and (49) thus suggest the following integral forms (Equation (50)):

$$\begin{gathered} \left[K_e\right]={{\displaystyle \int }}_V{\left[B_{\phi }\right]}^T\left[D_e\right]\left[B_{\phi }\right]dV \\ \left[K_p\right]={{\displaystyle \int }}_V{\left[B_{\phi }\right]}^T\left[D_p\right]\left[B_{\phi }\right]dV \\ {\left[K_V\right]}^T={{\displaystyle \int }}_V{\left[B_{\phi }\right]}^T{\langle 1 0\rangle }^T\left[N_p\right]dV \\ \left[K_h\right]={{\displaystyle \int }}_V{\left[B_{\theta }\right]}^T\left[k\right]\left[B_{\theta }\right]dV, \end{gathered}$$

(50)

To solve Equations (48) and (49), a solution procedure can be applied. The variables include the soil displacement $U$ and the pore water pressure $p_w$ at that time. One of the procedures used in this paper entailed the application of the finite difference method in time combined with Galerkin’s method to obtain theta values. The approximate values refer to the external forces $F$ acting on the soil surface, external/internal water flow in the soil $Q$, pore water pressure $p_w$, displacement of the soil mass, and rate of movement in terms of the velocity.

$$\begin{gathered} \left\{F\right\}=\left\{F_a\right\}+\theta \left\{\left\{F_b\right\}-\left\{F_a\right\}\right\} \\ \left\{Q\right\}=\left\{Q_a\right\}+\theta \left\{\left\{Q_b\right\}-\left\{Q_a\right\}\right\} \\ \left\{p_w\right\}=\left(1-\theta \right)\left\{p_{wb}\right\}+\theta \left\{p_{wa}\right\} \\ \left\{U\right\}=\left(1-\theta \right)\left\{U_b\right\}+\theta \left\{U_a\right\} \\ \left\{\dot{U}\right\}=\left(\left\{U_b\right\}-\left\{U_a\right\}\right)/∆t, \end{gathered}$$

(51)

where $\left\{F_a\right\}$, $\left\{Q_a\right\}$, $\left\{p_{wa}\right\}$, and $\left\{U_a\right\}$ are the initial vectors; $\left\{F_b\right\}$, $\left\{Q_b\right\}$, $\left\{p_{wb}\right\}$, and $\left\{U_b\right\}$ are the vectors after time interval $∆t$; $\theta$ is the value of Galerkin’s parameter in the numerical time evolution solution; and $∆t$ is the time step to avoid numerical oscillations according to the following relationship (for triangular finite elements) [68]:

$$∆t\le \frac{\alpha }{1-\theta }, {\alpha }_{min}=\frac{2A}{9C_V},$$

(52)

where $A$ is the area of the linear triangular finite element and $C_V$ is the soil consolidation coefficient.

Finally, the elastoplastic coupled model algorithm was developed in Fortran code to visualize the saturated consolidation phenomenon and its temporal evolution.

4. Numerical Examples and Results

To validate the elastoplastic coupled model, two numerical examples were examined: (1) a foundation slab covering the entire soil surface, in which vertical flow was allowed only at the bottom boundary, representing the simulation of typical Terzaghi’s one-dimensional consolidation; and (2) an isolated footing partially covering the surface, in which vertical flow was allowed at the bottom and top boundaries. In both cases, lateral displacement was restricted.

4.1. Case I: Vertical Flow (Bottom Boundary) under Lateral Restrictions

The typical case of one-dimensional consolidation due to a rectangular foundation slab of length $L=5.0\mathrm{m}$ and infinite width (plane strain state) transferring a uniform load $F=-80\mathrm{kPa}$ $(F=-8.0\mathrm{Ton}/{\mathrm{m}}^2$) onto a clayey soil stratum with thickness $H=3.0\mathrm{m}$ was analyzed. The $U_x$ and $U_y$ displacements were restricted both at the lateral boundaries and at the bottom boundary of the soil. In addition, a sand layer occurred below the clay layer, permitting the bottom boundary of the clay layer to drain and indicating a pore pressure of zero.

Similarly, the lateral boundaries did not allow water flow. Therefore, the vertical flow (bottom boundary) was studied. The mesh spacing was set to 0.25 m to use the calibration curves obtained. Figure 1 shows the soil domain under study and

Table 1. Properties of the saturated soil for mechanical analysis of the foundation.

Soil Property	Symbol	Magnitude
Elastic modulus	$E$	20,000 kPa  2000 Ton/m2
Poisson’s ratio	$\upsilon$	0.35
Coefficient of permeability (x direction)	$k_x$	1.184 × 10−4 m/day  1.37 × 10−9 m/s
Coefficient of permeability (y direction)	$k_y$	1.184 × 10−4 m/day  1.37 × 10−9 m/s
Preconsolidation stress	${\sigma }_0$	−120.0 kPa  −12.0 Ton/m2
Internal friction angle	$\phi$	30°
Lambda	$\lambda$	0.20
Kappa	$\kappa$	0.02
Time for analysis	$t$	1379.07 days
Time step	$∆t$	13.93 days

lists the soil properties needed for the mechanical analysis. In this analysis, the time step was $∆t=13.93$ days.

Figure 2 shows the pore water pressure distribution at time steps of 13.93, 682.57, and 1379.07 days. The pore water pressure distribution gradually decreased with increasing depth because the change in the applied stress distribution also decreased with increasing depth. Thus, at 13.93 days, the pore pressure value reached −79.8 kPa (−7.98 Ton/m²) below the foundation and was very close to the applied stress value of −80 kPa (−8 Ton/m²); however, when approaching the lower boundary, the pore pressure reached zero because flow was permitted. Therefore, when considering only one impermeable boundary, the water contained in the soil domain needed a longer drainage time and covered a greater distance while the filtration rate was very low, resulting in a required analysis time of 1379.07 days.

Additionally, when comparing Figure 2a–c, we observed that the pore pressure decreased over time until a value close to zero was reached. This behavior was appropriate and correct when analyzing the soil consolidation problems, which evolved over time.

In addition, the Mandel–Cryer phenomenon as described by Conte [7] was observed due to the sensitivity of the model, which caused the excess pressure to increase to infinity. For this reason, it was required to calibrate the model and obtain calibration curves to reduce this effect by determining the imminent point at which the pore water pressure began to dissipate.

Figure 3 shows the symmetrical behavior of horizontal displacement $U_x$. The foundation produced displacement to the right (positive) and left (negative) from the centerline of the domain. The largest displacements occurred in the top third (near the top boundary) while displacement was restricted at the lateral and bottom boundaries; hence, the displacements were very small or negligible. As shown in Figure 3a–c, the maximum value of −0.0007 m (−0.700 mm) was observed at 13.93 days, but over time, this value increased to −0.000713 m (−0.713 mm). This indicated that the largest displacements occurred at the beginning of the consolidation process; i.e., at the time of foundation emplacement. Subsequently, the displacements decreased due to soil rearrangement, which was attributed to pore pressure dissipation throughout the domain.

The behavior of vertical displacement $U_y$ is shown in Figure 4. It was evident that the vertical displacements at the top boundary were greater and decreased with increasing depth due to the distribution of the applied stress. As shown in Figure 3a, during the first 13.93 days after the beginning of the consolidation process, in the superficial corners, the magnitude of the displacements at −0.01579 m (−1.579 cm) was larger than that of the displacements in the central zones at −0.008755 m (−0.8755 cm), producing a concave profile of superficial displacements due to the large FE size; thus, with finer meshing, this effect could be reduced. This effect could only be observed through the coupled model and could not be observed via classical mathematical functions. In contrast, the $U_y$ displacements near the bottom boundary were already very small or zero. In addition, at 682.57 days (equivalent to 50% of the analysis time), the maximum $U_y$ displacements changed to 3.368 cm, and at the end of the analysis time, the displacements reached 3.373 cm (as shown in Figure 4b). This indicated that significant displacements could occur after 50% of the consolidation time, which was consistent with Terzaghi’s consolidation theory.

Figure 5 shows the evolution of the profile of the vertical displacements of the superficial nodes, revealing that the vertical displacements of the superficial boundary produced by the foundation were negative; i.e., settlement occurred. Figure 5 shows different soil superficial profiles for each time increment (a color line was obtained every 13.93 days); i.e., after 13.93 days, the foundation produced a settlement of −0.80 cm in the center, while in the corner it was −1.50 cm (superior magenta line). At the end of the time analysis, the settlements reached a magnitude of −2.65 cm in the center and −3.30 cm in the corner (bottom magenta line). Based on this finding, we deduced that the soil stratum volume decreased. This behavior was consistent with the type of analysis performed.

Figure 6 shows the finite element mesh of the soil stratum and thus the behavior (elastic and/or plastic) they were under according to the value of the parameter $f$. When considering this meshing, the behavior of the FEs within the soil stratum could be identified at any time step during the application of the foundation-applied stress. As shown in Figure 6a, for Δt = 13.93 days, 90% of the finite elements within the domain exhibited elastic behavior (stress and strain) or blue zones. Over time (as shown in Figure 6b), 75% of the FEs exhibited plastic behavior (red zones), and at the end of the analysis time; i.e., 1379.07 days (as shown in Figure 6c), all FEs exhibited plastic behavior. Because the low permeability impeded the pore pressure dissipation process and due to the lateral restrictions; i.e., at the bottom and top, the soil was subjected to plastic strains (irreversible) at the end of the consolidation process.

4.2. Case II: Vertical Flow (Top and Bottom Boundary) under Lateral Restrictions

Additionally, a rectangular isolated footing of length $L=2.5\mathrm{m}$ and infinite width (plane strain state) was analyzed. The footing was placed at the center below the top boundary of the domain and applied a uniform stress $F=-80\mathrm{kPa}$ $(F=-8.0\mathrm{Ton}/{\mathrm{m}}^2$) on the clayey soil stratum with thickness $H=3.0\mathrm{m}$.

The top boundary permitted water flow in the unloaded zones (to the left and right of the footing). Similarly, the bottom boundary also allowed water flow, indicating that the pore pressure at these boundaries was zero. $U_x$ displacements were restricted at the lateral boundaries, while at the bottom soil boundary, $U_x$ and $U_y$ were also restricted. The meshing was set to 0.25 m in order to use the calibration curves obtained. Figure 7 shows the study domain loaded by the footing, while the soil properties needed for mechanical analysis are listed in

Table 2. Properties of the saturated soil for mechanical analysis of the isolated footing.

Soil Property	Symbol	Magnitude
Elastic modulus	$E$	20,000 kPa  2000 Ton/m2
Poisson’s ratio	$\upsilon$	0.35
Coefficient of permeability (x direction)	$k_x$	1.184 × 10−4 m/day  1.37 × 10−9 m/s
Coefficient of permeability (y direction)	$k_y$	1.184 × 10−4 m/day  1.37 × 10−9 m/s
Preconsolidation stress	${\sigma }_0$	−120.0 kPa  −12.0 Ton/m2
Internal friction angle	$\phi$	30°
Lambda	$\lambda$	0.20
Kappa	$\kappa$	0.02
Time for analysis	$t$	877 days
Time step	$∆t$	4.385 days

. According to the calibration curve (see Figure A3 in Appendix B), to perform the analysis, the time step was set to $∆t=4.385$ days, while the analysis time was set to $t=877$ days.

Through numerical analysis, Figure 8 reveals that the pore water pressure of the soil subjected to the above conditions exhibited a very consistent behavior. After the analysis time and time step of the calibration curve were obtained (see Figure A3 and Figure A4 in Appendix B), it was observed that the pore water pressure did not approach infinity. Conversely, it was evident that the pressure began to decrease over time. Thus, Figure 8a shows that after 4.385 days of consolidation, the pore water pressure was −57.89 kPa (−5.789 Ton/m²), which represented a reduction of 27.64% with respect to the stress applied by the footing (−80 kPa). Similarly, after 434.115 days (as shown in Figure 8b), the pore water pressure decreased 99.65% (−0.276 kPa) compared to the applied stress. At the end of the analysis time; i.e., 872.615 days (as shown in Figure 8c), it was evident that the pore water pressure approached zero, representing a typical consolidation phenomenon.

In addition, Figure 8 shows that the pore water pressure at the top boundary, which was not loaded (to the left and right of the footing), reached zero. Similarly, the bottom boundary exhibited the same behavior because these boundaries allowed water to drain and the pore pressure within the domain to dissipate.

Figure 9 shows the evolution of the horizontal displacements $U_x$ with the consolidation time. Figure 9 shows symmetric behavior of the horizontal displacements $U_x$. On the basis of performed numerical analyses, we obtained, at the center of the top boundary, produced displacements to the left (negative) and right (positive) from the centerline of the stratum. After 4.385 days (as shown in Figure 9a), the maximum value ranged from −0.004102 m (−0.4102 cm) to 0.004205 m (0.4205 cm) just at the strip below the footing. However, with increasing depth, the displacements were very small and almost negligible. Additionally, over time (as shown in Figure 9b,c), this value decreased to 0.00314 m (0.314 cm). This indicated that the largest $U_x$ displacements occurred at the beginning of the consolidation process. Thus, when the consolidation time exceeded 50% of the analysis time, the displacements continued to increase but at a smaller magnitude due to rearrangement of the solid soil particles and pore water pressure dissipation throughout the domain.

Figure 10 shows the distribution and behavior of vertical displacements $U_y$. During the first 4.385 days (as shown in Figure 10a) after the start of the consolidation process, at the top boundary of the domain; i.e., the loaded zone (central strip and interior of the domain), the displacements were negative (typical of loaded zones) and reached a maximum value of −0.01202 m (−1.202 cm). However, in the unloaded zones, the vertical displacements were positive, with a value of +0.006679 m (+0.6679 cm). This indicated that the isolated footing settled across the entire contact surface with the soil, which produced a combination of vertical and horizontal displacements within the domain that generated, in turn, upward (positive) vertical displacements in the unloaded zones at the top boundary. Additionally, after 434.115 (Figure 10b) and 872.615 days (Figure 10c), the footing generated settlement within the central strip of the domain; while at the unloaded surface, lateral boundaries, and bottom boundary, the displacements were zero. Thus, by analyzing the $U_y$ displacements shown in Figure 10, it could be inferred that the domain experienced settlement in the loaded zone and uplift in the unloaded zones, indicating permanent strains.

The above is clearly shown in Figure 11. This figure shows the variation in the soil superficial profile of the vertical displacements of the superficial nodes, providing a clear view of the evolution of the vertical displacements at the top boundary produced by load application to this boundary and pore water pressure dissipation throughout the interior of the domain (soil). Similarly, we observed how the superficial profile was strained, so we could deduce that the loaded area experienced settlement (displacements below the red dotted line), revealing a decrease in volume; while the unloaded area experienced uplifts (displacements above the red line), revealing an increase in volume.

Figure 12 shows the meshing domain of the FEs and indicates whether the behavior of the FEs is elastic or plastic. This determination was based on the yield surface $f$; for $f\le 0$, the FEs exhibited elastic behavior, but for $f>0$, the FEs exhibited plastic behavior. As a result, at a time step of 4.385 days, a significant portion of the finite elements within the domain experienced plastic stress and strain behavior (red zones in Figure 12a). Over time until the analysis time was reached, the finite elements that exhibited plastic behavior (below the unloaded zones) changed to exhibit elastic behavior (change from red to blue in Figure 12b,c).

Based on the above, a direct relationship could be observed between the vertical $U_y$ and horizontal $U_x$ displacements because the area below the foundation (loaded) exhibited plastic behavior. Furthermore, the secondary effects produced by loading also subjected the unloaded areas to plastic behavior. The reason for this phenomenon was that the soil flowed downward in the central loaded area, displacing the soil toward the sides and finally lifting the sides and deforming the surface soil outward in the unloaded areas. When the pore water pressure was very low—practically zero (at 434.115 days)—the finite elements no longer exhibited plastic behavior but changed to exhibit elastic behavior and remained so until the end of the time analysis process. The reason for this phenomenon was that there were no longer any overpressure-producing permanent strains.

Notably, the proposed elastoplastic coupled model also captured the behavior of (a) the volumetric strain ${\epsilon }_v$, (b) the deviatoric strain ${\epsilon }_q$, (c) the mean stress ${\sigma }_p$, (d) the deviatoric stress ${\sigma }_q$, and (e) the volumetric variation $∆{\epsilon }_v$; the latter reflected the final value compared to the total volume change in the domain as either positive (an increase in volume) or negative (a decrease in volume).

Finally, a comparative analysis between the case studies in the present work revealed that under the same soil properties, magnitude of the applied stress, and restrictions at the lateral and bottom boundaries in both cases, but with different permeable boundaries, the analysis time increased from Cases I to II (872.615 and 1393 days, respectively) and the time step also increased from 4.385 to 13.93 days.

Based on

Table 3. Comparative analysis in conditions with and without superficial flow.

With Superficial Flow			Without Superficial Flow
Time Step	Water Pore Pressure		Time Step	Water Pore Pressure
(Days)	(kPa)	(Ton/m2)	(Days)	(kPa)	(Ton/m2)
4.385	−57.89	−5.789090	13.93	−79.87	−7.987030
434.115	−0.132	−0.013244	682.57	−0.219	−0.021974
872.615	−0.000	−0.000027	1379.07	−0.000	−0.000046

, we deduced that when flow was permitted at the top boundary, the time needed for the pore water pressure to decrease or dissipate was less than that in the analysis case that did not allow surficial flow. This result was congruent with the behavior of the consolidation phenomenon under both conditions; i.e., the time needed for the pore water pressure to dissipate was shorter because water flowed toward the two boundaries, which suggested that the travel distance of water through the soil was also shorter; hence, the pore water pressures were lower. In the opposite case, with the top boundary sealed, water traveled a greater distance to be discharged, and as a result, the pore water pressures were higher.

5. Discussion

When analyzing the pore water pressure ($p_w$) behavior, it was possible to observe the Mandel–Cryer phenomenon described by Conte [7], which explains that an excess pressure exists at the beginning of such a process due to the applied stress; this excess pressure increases and then begins to decrease until it is completely dissipated. The proposed model showed this excess pressure but it increased to infinity; for this reason, it was necessary to calibrate the model and obtain the calibration curves shown in Appendix B. Thus, the behavior of the pore water pressure was consistent with those reported in the literature reviewed. However, a direct comparison was not possible because the reviewed models showed the analysis of a specific point of the stratum while the proposed model analyzed the entire soil stratum by means of color maps.

Regarding horizontal displacements ($U_x$), the classical theory of consolidation describes that these displacements are zero; however, the proposed model demonstrated small displacements, but these did exist. Likewise, after analyzing the horizontal displacements, it was observed that in the two cases of study, the settlements were symmetrical when tracing a line center. In addition, it was remarkable to identify higher displacements occurring in the first quarter of the analysis time that were becoming lower and lower due to this. This behavior was correct and consistent with the consolidation phenomenon.

The vertical displacements ($U_y$) in the two case studies showed higher-magnitude displacements occurring in the superficial zone (below the foundation) that decreased as the depth increased. This was consistent with the distribution of the applied stress because it decreased as the depth increased, having less of an effect each time. Likewise, the displacements evolved over time, presenting those of higher magnitude at the beginning of the process, but as time progressed, they decreased considerably. Additionally, due to the Mandel–Cryer phenomenon and the combination of the vertical and horizontal displacements, it was possible to observe the superficial profile of the soil at each time step. The profile revealed that when the soil was subjected to an applied stress, it experienced settlements (negative displacements), but in nonloaded zones (around the foundation), expansions or swellings (positive displacements) were manifested. Such behavior cannot be determined using classical theories, nor did the literature models reviewed report it.

On the other hand, the elastoplastic model evaluated the value of the yield surface f to identify at each time step the behavior (elastic and/or plastic) of each finite element in the mesh. As a result of this analysis, at the beginning of the process in the two case studies, the FEs under the foundation were in elastic behavior (blue zones); then the FEs evolved to plastic behavior (red zones). Therefore, at the end of the analysis process, the loaded zones experienced permanent strains. Likewise, in the periphery of the foundation, the FEs were under plastic behavior, which evidenced that the swellings on the surface were permanent strains. Finally, the calibration curves revealed that when there were two permeable vertical boundaries instead of one, the analysis time and time step were shorter. This indicated that the consolidation process would take less time because the water could flow through two boundaries and permit the water pore pressure to dissipate faster. This behavior is typical in consolidation processes and was in accordance with the classical theory and literature models reviewed.

6. Conclusions

In this paper, an elastoplastic coupled model (flow–mechanical–critical state) was formulated and developed to study the consolidation of saturated soils under effective stresses and their transition from elastic to elastoplastic behavior and/or vice versa. Based on this investigation, the following conclusions could be derived:

The variation in consolidation in saturated soils was determined by the rate at which the fluid pressure in the soil could be reduced. Thus, the process of strain (mechanical behavior) and the distribution of the water pressure (hydraulic behavior) in soils is a coupled flow–mechanical process. Thus, the use of the coupled model represented an integral solution to the coupling problem, yielding approximate solutions that reproduced the hydromechanical behavior of the soil in a complete way.

The proposed elastoplastic coupled model indicated that the resulting coupling equations confirmed that the approach and formulation were adequate and exhibited correct coupling and congruence in the matrix product process.

Simulations using the proposed elastoplastic coupled model algorithm provided fast and easy results due to its flexibility, since it permitted an infinite combination of conditions: (1) boundaries, (2) loads, and (3) flow patterns. The first condition could consider restrictions at the bottom and lateral boundaries or not. The second condition could consider constant or time-varying values. The third condition could allow internal or external flow, which in turn could exhibit constant or evolving values over time.

The numerical results revealed an appropriate behavior of the consolidation phenomenon under the conceptual framework of the effective stress of saturated soil, supporting its feasibility and reliability. The water pore pressure behavior presented an overpressure at the beginning (Mandel–Cryer effect), but as time passed it tended to zero. The vertical and horizontal displacements were greater below the foundation and decreased with depth and decreased with time. The soil superficial profile evolved with time and showed the permanent strains at the surface (settlements and/or swellings). The model showed that the FE in the zones below the foundation were initially in elastic behavior, but at the end of the process, they evolved to the plastic regime and showed permanent strains. Therefore, the model was consistent with Terzaghi’s consolidation theory, and if two boundaries were permeable instead of one, the analysis time was shorter and would permit water to flow more quickly.

7. Future Work

The analysis process does not end here: the different combinations of boundary conditions, loading levels, flow patterns, and soil properties considered in the many simulations we performed to improve the model must be examined. However, due to space limitations, only the most recurrent and significant consolidation problems were included.

The proposed elastoplastic coupled model was formulated and prepared while considering saturated soil consolidation but could be modified for unsaturated soil consolidation after coupling the variation in the hydraulic conductivity with respect to the saturation degree obtained from the soil water retention curve (SWRC). This process should be incorporated into the permeability matrix.