Next Article in Journal
Gas Permeability Behavior in Frozen Sand Controlled by Formation and Dissociation of Pore Gas Hydrates
Next Article in Special Issue
Parametric Study of Lateral Loaded Blade Pile in Clay
Previous Article in Journal
Recycling and Reuse of Mine Tailings: A Review of Advancements and Their Implications
Previous Article in Special Issue
Evaluating the Effect of Different Stress Path Regimes on Borehole Deformation Using Convergence Measuring Device
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Coupled Large Deformation Finite Element Formulations for the Dynamics of Unsaturated Soil and Their Application

by
Nadarajah Ravichandran
1,* and
Tharshikka Vickneswaran
2
1
Glenn Department of Civil Engineering, Clemson University, 202A Lowry Hall, Clemson, SC 29634, USA
2
Civil and Environmental Engineering, University of Louisville, 218 Eastern Pkwy, Louisville, KY 40208, USA
*
Author to whom correspondence should be addressed.
Geosciences 2022, 12(9), 320; https://doi.org/10.3390/geosciences12090320
Submission received: 21 June 2022 / Revised: 22 August 2022 / Accepted: 23 August 2022 / Published: 27 August 2022
(This article belongs to the Collection New Advances in Geotechnical Engineering)

Abstract

:
Unsaturated soil is a three-phase medium with three interfaces, and the mathematical equations that represent its behavior must be developed in a fully coupled manner for accurately predicting its hydromechanical behavior. In this paper, a set of fully coupled governing equations was developed for the dynamics of unsaturated soil, considering the interaction among the bulk phases and interfaces. In addition to implementing the complete governing equations, a simplified formulation was developed for practical applications. The derivation of the finite element formulation considering all the terms in the partial differential equations resulted in a formulation called complete formulation and was solved for the first time in this paper. Another formulation called reduced formulation was derived by neglecting the relative accelerations and velocities of water and air in the governing equations. In addition, small and large deformation theories were developed and implemented for both formulations. To show the applicability of the proposed models, the dynamic behavior of an unsaturated soil embankment was simulated using both small and large deformation formulations by applying minor and severe earthquakes. The reduced formulation was found to be computationally efficient and numerically stable. The smaller displacements predicted by large deformation theories show that the results are consistent with the expected behavior. Large deformation theories are considered suitable when the geotechnical system undergoes large deformation and may lead to accurate prediction.

1. Introduction

The hydromechanical behavior of geotechnical structures, such as earth dams, levees, and near-surface soil that supports superstructures, is influenced by the degree of saturation (DOS) of the soil. The DOS of soil varies with climatic and hydrological parameters. In general, the soil is in an unsaturated state, and the saturated (DOS = 100%) and dry (DOS = 0) states are special cases of an unsaturated state. During hydrological and/or climatic events, such as rainfall, flood, and drought, the state of the soil changes continuously, resulting in variations in the behavior of geotechnical systems. To better design and construct sustainable and resilient geotechnical systems, the variation in the behavior of the soil when it transitions from one state to the other must be understood. In this study, numerical models that predict the hydromechanical behavior of soil are developed, considering the soil as unsaturated soil, and implemented.
Unsaturated soil consists of three bulk phases (solid, water, and air) and three interfaces (solid–water, water–air, and air–solid). The hydromechanical behavior of the soil is influenced by the bulk phases and interfaces and also by the interaction among them. Of the three interfaces, the water–air interface (contractile skin) has a significant influence on the hydromechanical behavior of unsaturated soil [1], and it must be incorporated at the governing equation level for to accurately predicting the behavior of soil. The behavior of the water–air interface is controlled by the amount of water in the unsaturated soil, particularly at the particle contacts. The contractile skin helps maintain the pressure difference between the water and air phases and results in the water pressure in an unsaturated soil system being always negative. These factors complicate studying unsaturated soil compared to fully saturated and/or dry soils.
The partial differential equations for the dynamics of unsaturated soil systems can be derived based on physical laws, such as the balance of mass, linear momentum, angular momentum, and the first and second laws of thermodynamics. When developing governing equations, the macroscopic quantities representing microscopic quantities, such as the area density of the water–air interface, must be identified and properly used in the derivation. A rigorous volume-averaging technique is used to accomplish this task [2].
The partial differential equations are typically simplified to develop numerically stable formulations without compromising the accuracy of the prediction. In the case of unsaturated soils, the relative accelerations of the water and air phases with respect to the solid phase are typically neglected in the solution procedure [3,4]. However, such an assumption has not been proved to be valid. Schrefler et al. [3] and Wei [5] derived finite element formulations by neglecting the relative accelerations of the water and air phases. They solved the finite element equations, considering solid displacement (u), water pressure (pl), and air pressure (pg) as the primary nodal unknowns ( u p l p g formulation). In the formulation, they used a two-dimensional, four-node quadrilateral element with continuous bilinear displacement (u) and pressure (pl and pg). However, it is shown that such formulation violates the Babuska–Brezzi conditions or the simpler Zienkiewicz–Taylor patch test [6] for solvability and convergence [7].
In general, small deformation theories are commonly used in the finite element simulation of soils, assuming that the geotechnical engineering structures experience small deformations. The true response of the geotechnical structures undergoing large deformations cannot be predicted correctly with small deformation theories. Therefore, large deformation theories must be incorporated while solving the partial differential equations using the finite element method [8,9]. Sanavia et al. [10] and Gawin and Sanavia [11] developed a formulation for partially saturated soils undergoing large deformations. From the example simulations, they showed strain and negative pore pressure localization and their effects on the hydromechanical behavior of partially saturated soils.
In this paper, a set of fully coupled partial differential equations for the dynamics of unsaturated soils was derived, considering the interaction among the bulk phases and interfaces. Next, two different finite element formulations (complete and reduced) of the governing equations were derived for analyzing unsaturated soil mechanics problems undergoing large deformations. In addition, small and large deformation theories were developed and implemented for both formulations. Finally, the dynamic behavior of an unsaturated soil embankment was simulated using both small and large deformation formulations by applying minor and severe earthquake acceleration–time histories.

2. Summary of Governing Equations

The partial differential equations for the dynamics of unsaturated soils were derived based on the laws of physics, such as the momentum balance, mass balance, energy balance, and the second law of thermodynamics. Since the amount of water in the soil system is directly related to the matric suction (air pressure minus the water pressure), a constitutive equation that relates the amount of water to the matric suction was established. In this study, the volume fraction of the water phase is assumed to be a function of the solid skeleton’s matric suction and volumetric strain, as shown in Equation (1).
n l = n l ( S , ε v )
where n l is the volume fraction of the water phase, ε v is the volumetric strain of the solid skeleton, and S is the matric suction defined by S = p g p l , where p g is the pore air pressure and p l is the pore water pressure.

2.1. Mass Balance Equations

When deriving the governing equations, the volume spanned by the solid phase is considered as the representative elementary volume (REV) and its motion is given by φ s ( X , t ) , where X is the material coordinate and t is the time. The water and air phases can move in (influx) and out (outflux) of the REV. This phenomenon is graphically shown in Figure 1. Therefore, there will be net flow across the closed REV. Such influx and outflux were considered in this study when deriving the governing equations.
The balance of mass in an REV for a bulk phase α is given by Equation (2).
m α ( Ω ) = Ω t + Δ t ( n α ρ α ) d Ω
where m α is the mass of the α -phase, n α is the volume fraction of the α -phase, and ρ α is the density of the α -phase.
  • Mass balance for the solid phase:
The following equation (Equation (3)) for the solid phase was derived, assuming that the solid particles are incompressible:
n ˙ + ( 1 n ) div ( v s ) = 0
where n is the porosity of the unsaturated soil system and v s is the velocity vector of the solid phase.
  • Mass balance for the water phase:
The final mass balance equation for the water phase was derived by incorporating the mass balance equation for the solid phase into the mass balance equation for the water phase, as shown in Equation (4).
( n l ε v ) u ˙ i , i + n l U ˙ i , i l + ( n l Γ l n l S ) p ˙ l + ( n l S ) p ˙ g = 0
where u is the displacement of the solid phase, U l is the displacement of the water phase, and Γ l is the bulk modulus of the water phase.
  • Mass balance for the air phase:
Similarly, the final mass balance for the air phase was derived by incorporating the mass balance equation for the solid phase into the mass balance equation for the air phase, as shown in Equation (5).
( 1 n n l ε v ) u ˙ i , i + n g U ˙ i , i g + ( n l S ) p ˙ l + ( n g Γ g n l S ) p ˙ g = 0
where U g is the displacement of the air phase and Γ g is the bulk modulus of the air phase.

2.2. Momentum Balance Equations

Momentum balance equations can fully describe the motion of the unsaturated soil system for the soil (solid, water, and air mixture), water phase, and air phase. The momentum balance equations for the water and air phases are essentially the generalized Darcy’s flow equations in the flow direction. At the micromechanical level, the primary resistance to the flow of these fluids is the drag forces on the solid skeleton and the primary driving force is the fluid pressure gradients. The final momentum balance equations for soil, water, and air were derived, as shown in Equations (6)–(8), respectively.
  • Linear momentum balance for the mixture:
n s ρ s u ¨ j + n l ρ l U ¨ j l + n g ρ g U ¨ j g σ i j , i ρ b j = 0
  • Linear momentum balance for the water:
ρ l U ¨ j l ( k ^ i j l n l ) u ˙ i + ( k ^ i j l n l ) U ˙ i l + ( δ i j p l ) , i ρ l b j = 0
  • Linear momentum balance for the air:
ρ g U ¨ j g ( k ^ i j g n g ) u ˙ i + ( k ^ i j g n g ) U ˙ i g + ( δ i j p g ) , i ρ g b j = 0
where σ i j is the total stress tensor, b j is the gravitational acceis the gravitational acceleration vector, leration vector, k ^ i j l is the inverted permeability tensor of the water phase, k ^ i j g is the inverted permeability tensor of the air phase, and δ i j is the Kronecker delta.

3. Updated Lagrangian Formulation of Governing Equations

3.1. Application of the Principle of Virtual Work

The principle of virtual work requires that the virtual work performed when the soil body undergoes a virtual displacement δ u be equal to the external work performed by the body force and traction, i.e.,
W t + Δ t int = Ω t + Δ t σ t + Δ t i j δ e t + Δ t i j d t + Δ t Ω
W t + Δ t e x t = Ω t + Δ t ρ t + Δ t t + Δ t b i δ u i d t + Δ t Ω Ω t + Δ t n t + Δ t s ρ t + Δ t s u ¨ t + Δ t i δ u i d t + Δ t Ω Ω t + Δ t n t + Δ t l ρ t + Δ t l U ¨ t + Δ t i l δ u i d t + Δ t Ω Ω t + Δ t n t + Δ t g ρ t + Δ t g U ¨ t + Δ t i g δ u i d t + Δ t Ω + S t + Δ t T f t + Δ t i t δ u i d t + Δ t S
where σ t + Δ t i j are the Cartesian components of the Cauchy total stress tensor and e t + Δ t i j are the Cartesian components of the strain tensor. In the case of unsaturated soils, the body force, the inertial force of the solid skeleton, the inertial forces of the pore fluids, and the surface traction contribute to the external work. Similar equations were derived for the motion of the pore fluids. It should be noted that there are two major difficulties that exist when applying these equations for developing large deformation theories, which involve rotation and a change in configuration. First, the configuration at time t + Δ t is unknown, and the integrals over the volume Ω t + Δ t and surface S t + Δ t cannot be evaluated before calculating the equilibrium position at time t + Δ t . The second difficulty is the presence of total stress in the internal work equation. Since the total stress does not directly influence the mechanical behavior of the soil, the principle of net stress, which is expressed in terms of Cauchy stress, was applied.

3.2. Stress–Strain Relationship for the Solid Skeleton undergoing Large Deformation

For the large deformation formulation, the rotation of the element must be separated from its deformation to accurately calculate the strain. During the elastoplastic deformation of the body from the reference configuration to the current configuration, the material undergoes elastic reversible deformation and plastic irreversible deformation. The motion of the body from Ω t to Ω t + Δ t (a reference configuration Ω t , a virtual intermediate configuration Ω i , and a current configuration Ω t + Δ t ) is considered in two steps, the motion of the body from Ω t to Ω i and then from Ω i to Ω t + Δ t , as shown in Figure 2. The motion from Ω t to Ω i is purely plastic and irreversible. Therefore, the configuration Ω i can be considered an unstressed configuration. The motion from Ω i to Ω t + Δ t is purely elastic and reversible. The deformation gradient for the motion from Ω t to Ω i is denoted by F p , and the deformation gradient for the motion from Ω i to Ω t + Δ t is denoted by F e in matrix form. When the motion from Ω t to Ω t + Δ t is continuous, the deformation gradient has the following non-cumulative representation in its plastic and elastic parts [12,13]:
F = x x ¯ x ¯ X = F e F p
where x is the spatial coordinate, X is the material coordinate, and x ¯ is the intermediate coordinate. Since the deformation measures are not linearly expressed in terms of displacements, the elastic and plastic components are not summable. Choosing a representation of the elastic part of deformation independent of rigid body motion (rotation), the deformation rate tensor D is given by the symmetric part of the velocity gradient L [13,14].
L = F ˙ e ( F e ) 1 + F e F ˙ p ( F p ) 1 ( F e ) 1 and D = D e + ( F e D p F e 1 ) s + ( F e W p F e 1 ) s
where L = D + W
where W is the skew-symmetric part of the velocity gradient. The subscript s indicates the symmetric part. If the elastic components of the total strain are assumed to be minor, which is true for most soils, F e 1 and the last term in Equation (12) is also small. Then, the rate of deformation reduces to the additive decomposition, as shown in Equation (13).
D = D e + D p
where D e and D p are the elastic and plastic parts of the total strain rate, respectively.
The corotational form of the stress–strain relationship for the elastoplastic material was expressed in the following form (Equation (14)):
σ ˙ = C e p : D
where C e p is the tangential elastoplastic stiffness tensor, which may be a function of the current state of the net stress, matric suction, strains, and some other internal variables.
Among the many forms of objective stress rates, the Green–Naghdi stress rate was used in this study, as shown in Equation (15).
σ = D D t σ Ω σ σ Ω T
where the angular velocity is given by Ω = R ˙ R T
The rotation tensor R was calculated using the polar decomposition theorem, which states that any deformation gradient tensor F can be multiplicatively decomposed into a product of an orthogonal matrix R and a symmetric tensor U , called the right stretch tensor. The deformation gradient in the index notation is given by F i j = x i X j = R i k U k j .
By rearranging the objective rates and applying the net stress principle to the total objective rate, the net stress is expressed as:
D σ i j D t = σ i j n δ i j p ˙ g + Ω i k σ k j + σ i k Ω k j T   and   D σ i j D t = C i j k l e p D k l δ i j p ˙ g + Ω i k σ k j + σ i k Ω k j T

3.3. Stress-State Variables for Unsaturated Soils

The net stress and matric suction are widely accepted as the stress-state variables to define the mechanical behavior of unsaturated soils [1] and were also used in this study. For large deformation analysis, an objective stress measure must be used to consider the effect of rotation. Therefore, the objectivity of the net stress must be derived before using appropriate objective stress measures to the net stress. The net stress is given by
σ = σ n p g I
where σ is the total stress tensor, σ n is the net stress tensor, and p g is the pore air pressure. The conventional solid mechanics sign convention is used in the above equation. Since the net stress is defined in terms of the Cauchy stress tensor, which is not an objective measure of stress, a suitable rate form was established for the net stress equation. Equation (18) was derived by taking the time derivative of Equation (17).
D D t σ = D D t σ n D D t p g I
where D D t indicates the material time derivative.
The net stress equation was rewritten in corotational form, as shown in Equation (19).
σ = σ n p ˙ g I
where σ is the objective total stress rate, σ n is the objective net stress rate, and p ˙ g is the pore air pressure. The objective form of the net stress equation was incorporated into the virtual work equation.

3.4. Boundary Conditions

Solid, water, and air displacements; solid traction; and pore pressure boundary conditions were considered in the derivation. These boundary conditions are specified on different portions of the boundary surface S t + Δ t of the soil body at a generic time t + Δ t and are defined as follows:
  • Solid displacement boundary condition:
u t + Δ t i = u ¯ t + Δ t i on S t + Δ t u
where u ¯ t + t i is the specified value of solid displacement on the boundary surface S t + Δ t u at time t + Δ t .
  • Water displacement boundary condition:
U t + Δ t i l = U ¯ t + Δ t i l on S t + Δ t U l
where U ¯ t + t i is the specified value of water displacement on the boundary surface S t + Δ t U at time t + Δ t .
  • Air displacement boundary condition:
U t + Δ t i = U ¯ t + Δ t i on S t + Δ t U g
where U ¯ t + t i g is the specified value of air displacement on the boundary surface S t + Δ t U g at time t + Δ t .
  • Traction boundary condition:
σ t + Δ t i j n j = f t + Δ t i t on S t + Δ t T
where f t + t i t is the specified value of traction on the boundary surface S t + Δ t T at time t + Δ t , n j is the unit normal, and σ t + Δ t i j is the total Cauchy stress tensor acting on the neighborhood of S t + Δ t T .
  • Pore water pressure boundary condition:
p t + t l = p ¯ t + t l on S t + t p l
where p ¯ t + t l is the specified value of the pore water pressure at time t + Δ t .
  • Pore air pressure boundary condition:
p t + t g = p ¯ t + t g on S t + t p g
where p ¯ t + t g is the specified value of the pore air pressure at time t + Δ t .

3.5. Incremental Equations and Newton’s Method

By knowing the equilibrium state of the soil body at time t , the state Π was defined by the known stresses, tractions, deformations, and history of the soil body. Let the right and left sides of the virtual work equation in the reference configuration be I ( Π ) and E ( Π ) , respectively. At time t + Δ t , a new equilibrium state was established for the body. Let Δ Π be the change in state, which is the solution of [ I ( Π + Δ Π ) E ( Π + Δ Π ) ] δ v = 0 . Denoting Π ¯ as a guess for the new equilibrium state, the above equation can be expanded around the new guessed state, as follows:
[ I ( Π + Δ Π ) E ( Π + Δ Π ) ] δ v = [ I ( Π ¯ ) E ( Π ¯ ) ] δ v + [ I ( Π ¯ ) Π ¯ Π ¯ E ( Π ¯ ) Π ¯ Π ¯ ] δ v +
where Π ¯ is the increment between the correct equilibrium state Π + Δ Π and the guessed equilibrium state Π ¯ and δ v is the virtual velocity. By taking first-order approximation, the above equation reduces to Equation (21), as follows:
[ I ( Π ¯ ) Π ¯ Π ¯ E ( Π ¯ ) Π ¯ Π ¯ ] δ v = [ I ( Π ¯ ) E ( Π ¯ ) ] δ v
The successive solution of the above equation for various trial states Π ¯ can be found until the right side of the equation becomes zero, i.e., Π ¯ equal Π + Δ Π and the equilibrium is satisfied.
By substituting the stress measures and replacing the rate form with the incremental form, we obtain the following equations:
S ˙ i j = J F i k 1 ( C k l m n e p D m n + Ω k m σ m l + σ k m Ω m l T ) F l j T   and   S i j = J F i k 1 ( C k l m n e p D m n + Ω k m σ m l + σ k m Ω m l T ) F l j T
Equation (22) was substituted into Equation (10), and the incremental internal virtual work equation was derived, as shown in Equation (23).
W t + Δ t int = Ω t ( J F i k 1 ( C k l m n e p D m n + Ω k m σ m l + σ k m Ω m l T ) F l j T ) δ ε t + Δ t i j d t Ω Ω t h t t + Δ t i j δ ε t t + Δ t i j d t Ω
Equation (23) was simplified by choosing the current configuration to coincide with the reference configuration. This led to the deformation gradient becoming the identity tensor, while all the stress measures remained the same. This choice of the reference configuration is called the updated Lagrangian method [13,14,15]. This method is straightforward to use in computer codes because it requires only the coordinates of the body to be updated after each iteration so that the current configuration is also the reference configuration.

4. Finite Element Forms of Complete and Reduced Governing Equations

Closed-form solutions to the partial differential equations for many real-world problems do not exist. Therefore, a numerical technique, such as the finite element method, must be used to find approximate solutions for the system of equations. The finite element method is a powerful and widely used numerical method in geotechnical engineering to solve complex partial differential equations. To investigate the influence of relative accelerations and velocities of pore fluids on the overall dynamic behavior of unsaturated soils, two different formulations were derived in this study. The system of equations that considers the relative accelerations and velocities is called the complete formulation. In contrast, the system of equations simplified by neglecting the relative accelerations and velocities of the pore fluids is called the reduced formulation. The predictions from these two equations are compared in later sections.

4.1. Complete Formulation

The system of equations describing the dynamics of unsaturated soils, as listed below, consists of five equations and five unknowns.
n s ρ s u ¨ j + n l ρ l U ¨ j l + n g ρ g U ¨ j g σ i j , i ρ b j = 0
ρ l U ¨ j l ( k ^ i j l n l ) u ˙ i + ( k ^ i j l n l ) U ˙ i l + ( δ i j p l ) , i ρ l b j = 0
ρ g U ¨ j g ( k ^ i j g n g ) u ˙ i + ( k ^ i j g n g ) U ˙ i g + ( δ i j p g ) , i ρ g b j = 0
( n l ε v ) u ˙ i , i + n l U ˙ i , i l + ( n l Γ l n l S ) p ˙ l + ( n l S ) p ˙ g = 0
( 1 n n l ε v ) u ˙ i , i + n g U ˙ i , i g + ( n l S ) p ˙ l + ( n g Γ g n l S ) p ˙ g = 0
By eliminating the water and air pressures in the momentum balance equations using the mass balance equations, the final set of equations (Equations (25a)–(25c)) was derived in terms of the displacement, velocity, and acceleration of solid, water, and air phases. In the finite element solution, the displacement fields of all three phases were considered the primary nodal unknowns, as shown in Figure 3a.
n s ρ s u ¨ j + n l ρ l U ¨ j l + n g ρ g U ¨ j g σ i j , i ρ b j = 0
ρ l U ¨ j l ( k ^ i j l n l ) u ˙ i + ( k ^ i j l n l ) U ˙ i l + μ 11 u i , i j + μ 12 U i , i j l + μ 13 U i , i j g ρ l b j = 0
ρ g U ¨ j g ( k ^ i j g n g ) u ˙ i + ( k ^ i j g n g ) U ˙ i g + μ 21 u i , i j + μ 22 U i , i j l + μ 23 U i , i j g ρ g b j = 0
where
μ 11 = ( a 12 b 21 ( a 22 a 11 a 12 a 21 ) a 22 b 11 ( a 22 a 11 a 12 a 21 ) )
μ 12 = ( a 22 b 12 ( a 22 a 11 a 12 a 21 ) )
μ 13 = ( a 12 b 22 ( a 22 a 11 a 12 a 21 ) )
μ 21 = ( a 21 b 11 ( a 22 a 11 a 12 a 21 ) a 11 b 21 ( a 22 a 11 a 12 a 21 ) )
μ 22 = ( a 21 b 12 ( a 22 a 11 a 12 a 21 ) )
μ 23 = ( a 11 b 22 ( a 22 a 11 a 12 a 21 ) )
a 11 = ( n l Γ l n l S ) , a 12 = ( n l S ) , b 11 = ( n l ε v ) , b 12 = n l
a 21 = ( n l S ) , a 22 = ( n g Γ g n l S ) , b 21 = ( 1 n n l ε v ) and b 22 = ( n g )
The water and air pressures can be calculated using
p l = μ 11 u k , k + μ 12 U k , k l + + μ 13 U k , k g   and   p g = μ 21 u k , k + μ 22 U k , k l + μ 23 U k , k g
Figure 3. Nodal variables for the (a) u-U-U and (b) u formulations in 2 dimensions.
Figure 3. Nodal variables for the (a) u-U-U and (b) u formulations in 2 dimensions.
Geosciences 12 00320 g003

4.2. Simplified Formulation

The reduced formulation was derived by neglecting the pore fluids’ relative accelerations and velocities. Such simplification results in three equations (the momentum balance equation for the mixture and the mass balance equations for the water and air phases). In this case, the momentum balance equation was solved considering the solid displacements as the primary nodal unknowns, as shown in Figure 3b, and the water pressure and air pressure were calculated using the mass balance equations. In this formulation, the changes in water and air pressures are directly related to the volumetric deformation of the solid skeleton since the flow of fluids does not occur. This formulation simulates the undrained behavior of unsaturated soils.
The final set of equations is summarized below (Equations (27)–(29)).
ρ u ¨ j σ i j , i ρ g j = 0   in   Ω t + Δ t
( n l + n l ε v ) u ˙ i , i + ( n l Γ l n l p c ) p ˙ l + ( n l p c ) p ˙ g = 0   in   Ω t + Δ t
( 1 n l n l ε v ) u ˙ i , i + ( n l p c ) p ˙ l + ( n g Γ g n l p c ) p ˙ g = 0   in   Ω t + Δ t

4.3. Finite Element Formulation

The finite element formulations of the two formulations were derived using Galerkin’s weighted residual method, considering isoparametric four-node quadrilateral elements. The time integration was performed using the Hilber–Hughes–Taylor-α method, which is a standard procedure in structural and geotechnical earthquake engineering. It allows for energy dissipation and second-order accuracy, which is not possible with the regular Newmark’s method. The readers may refer to Wang and Chester [16], Wang et al. [17], and Ravichandran [8] for details on the derivation of weak form, shape functions, element-level residuals, and tangents.

5. Example Simulations

5.1. Finite Element Model

The finite element mesh for the unsaturated soil embankment used in this study is shown in Figure 4. The base of the embankment was assumed to be impermeable and fixed in both horizontal and vertical directions during the dynamic analyses. For the other three sides of the embankment, no displacement boundary conditions were applied so that they could move in any direction during shaking, and closed boundaries (no flow) were applied for the flow of water and air.

5.2. Constitutive models and model parameters

Stress–Strain Relationship of the Soil Skeleton

The stress–strain behavior of the solid skeleton was represented by an elastoplastic constitutive model based on the bounding surface concept. The schematic of the bounding surface model on stress-invariant space is shown in Figure 5. The original three-surface model for saturated cohesive soils [18] was modified for unsaturated soils by Muraleetharan and Nedunuri [19]. Additional parameters related to matric suction have been incorporated into the original model. This model uses the net stress ( σ i j p g δ i j ) and matric suction ( S ) as the stress-state variables. The modifications to the base model to incorporate the suction effects were based on the concepts proposed by Alonso et al. [20], Wheeler and Sivakumar [21], and Wheeler [22] for unsaturated soils. The simulations shown in this paper were carried out for an Oklahoma soil called Minco Silt (liquid limit = 28.0, plastic limit = 20.0, and USCS classification = CL). The model parameters for Minco Silt were obtained using the laboratory tests performed by Ananthanathan [23] and Vinayagam [24]. A complete list of the model parameters for Minco Silt and their values are listed in Table 1.

5.3. Soil Water Characteristic Curve (SWCC) and Model Parameters

Among the many SWCCs available in the literature, the model (Equation (30)) proposed by Brooks–Corey [25] was used in this study.
{ Θ = 1 S < a Θ = ( S a ) n S > a
where a and n are the fitting parameters and Θ is the dimensionless water content given by Θ = θ - θ r θ s θ r . The subscripts s and r indicate the saturated and residual values of the volumetric water content, θ , respectively. The SWCC parameters were determined by adjusting the model parameters until the model matched the experimental curve [23]. The model parameters are listed in Table 2. The relative permeability was calculated using the SWCC and the saturated permeability. The saturated permeability of 1 . 02 × 10 - 8 m/s and the initial DOS of 0.43 were used in this analysis. The initial pore water and pore air pressures were −30.0 and 0.0 kPa, respectively. The relative permeability of the water phase calculated using the van Genuchten method [25] was in the order of 10 - 3 . This implies that the unsaturated soil permeability is 1 . 02 × 10 - 11 m/s for the given initial conditions.

5.4. Static Analysis: Determination of the Initial Condition for Dynamic Analysis

The initial stresses for dynamic analysis using an elastoplastic constitutive model were determined through static (gravity) analysis. The static analysis to compute the initial stresses was performed by choosing the time integration parameters as α = 0 , β = 1.0 , and γ = 1.5 [26,27]. The gravity load was increased, as shown in Figure 6, for the static analysis. Only the stresses were transferred to the dynamic analysis, and all the nodal and element fields were reset to zero at the end of the static analysis.

5.5. Dynamic Analysis

For the dynamic analysis, the time integration parameters were changed to α = 0.3 , β = 0.4225 , and γ = 0.8 [28] and the model embankment was shaken with the horizontal and vertical motions, as shown in Figure 7.

5.6. Comparison between Small and Large Deformation Analyses

The response of the unsaturated soil embankment discussed in the previous section was analyzed using both small and large deformation models using only the reduced formulation. The horizontal and vertical displacement–time histories at nodes N073 and N053 are shown in Figure 8 and Figure 9, respectively. The time histories of the pore water pressure, pore air pressure, matric suction, and DOS in element E022 are shown in Figure 10. From these figures, it is clear that the model did not undergo significant deformation during the shaking. To increase the magnitude of the deformation, severe shaking was applied to the model by multiplying the amplitudes of the motions, as shown in Figure 7, by 6, and then the results from the large and small deformation codes were compared. The horizontal and vertical displacements at nodes N073 and N053 are shown in Figure 11 and Figure 12, respectively. The time histories of the pore water pressure, pore air pressure, matric suction, and DOS in element E022 are shown in Figure 13. We found that the large deformation analysis provides smaller vertical and horizontal displacements compared to the small deformation analysis. This is due to the initial stiffness contribution to the large deformation analysis at every time step coming from the updated Lagrangian formulation [15]. We also found that the large deformation analysis shows a minor change in pore water and air pressures. This response is consistent with the deformation response of the embankment.

5.7. Comparison of Computational Efficiency of Complete and Reduced Formulations

In the case of a four-node quadrilateral element for two-dimensional analysis, the complete formulation has 24-element degrees of freedom (DOF) and the reduced formulation has 8-element DOF. In the case of an eight-node brick element for three-dimensional analysis, these numbers are 72 and 24 for the full formulation and the reduced formulation, respectively. These numbers provide an idea of the computational effort needed to solve these formulations. The embankment problem described earlier was run on an Intel-Xeon processor with a 3.0 GHz clock speed processor. The complete formulation took 1160 min, and the reduced formulation took 60 min to analyze 30 s of shaking using the same number of steps in both analyses. The average number of global iterations for the complete formulation was 5 and for the reduced formulation was 3. That is, in this example, the complete formulation required 20 times more computational time than that required by the reduced formulation.

6. Concluding Remarks

The following conclusions were made from this study:
  • The large and small deformation theories for unsaturated soils were implemented within a finite element framework, and a large deformation analysis of the unsaturated soil was performed for the first time.
  • The complete formulation for unsaturated soil was successfully solved for the first time. The effects of relative fluid accelerations and velocities on the overall behavior of unsaturated soils were estimated and compared. The reduced formulation is computationally efficient and numerically stable. It captures the overall behavior reasonably well for the soil considered in this study (Minco Silt) and can be used for the evaluation of earthquake effects on similar soils. However, caution should be exercised in using the reduced formulation for soils with larger permeability values.
  • The small deformation analysis was observed to predict larger displacements compared to the large deformation analysis. Therefore, a large deformation analysis will likely produce an economical design of a geotechnical engineering structure.
  • All predictions made by the computer code TeraDysac should be validated against experimental results. These validations are currently in progress.

Author Contributions

Writing—original draft, N.R.; writing—review and editing, T.V. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Data sharing not applicable. No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Fredlund, D.G.; Rahardjo, H. Soil Mechanics for Unsaturated Soils; A Wiley-Interscience Publication: New York, NY, USA; John Wiley and Sons, Inc.: New York, NY, USA, 1993; p. 517. ISBN 0-471-85008-X. [Google Scholar]
  2. Hassanizadeh, S.M.; Gray, W.G. General Conservation Equation for Multi-Phase System: 1. Averaging Procedure. Adv. Water Resour. 1979, 2, 131–144. [Google Scholar] [CrossRef]
  3. Schrefler, B.A.; Simoni, L.; Xikui, L.; Zienkiwicz, O.C. Mechanics of Partially Saturated Porous Media. In Numerical Methods and Constitutive Modelling in Geomechanics; Desai, C.S., Gioda, G., Eds.; CISM Lecture Notes; Springer: Wein, Vienna, 1990; pp. 169–209. [Google Scholar]
  4. Muraleetharan, K.K.; Wei, C.F. U_DYSAC2: Unsaturated Dynamic Soil Analysis Code for 2-dimensional Problems. In Computer Code; University of Oklahoma: Norman, OK, USA, 1999. [Google Scholar]
  5. Wei, C.F. Static and Dynamic Behavior of Multi-Phase Porous Media: Governing Equations and Finite Element Implementation. Ph.D. Dissertation, University of Oklahoma, Norman, OK, USA, August 2001. [Google Scholar]
  6. Zienkiewicz, O.C.; Taylor, R.L. The Finite Element Method, 5th ed.; Solid Mechanics, Butterworth-Heinemann Publication, Linacre House, Jordan Hill: Oxford, UK, 2000; Volume 2, p. 459. ISBN 0-7506-5055-9. [Google Scholar]
  7. Brezzi, F.; Fortin, M. Mixed and Hybrid Finite Element Methods; Springer: Berlin/Heidelberg, Germany; New York, NY, USA, 2000. [Google Scholar]
  8. Ravichandran, N. A Framework-Based Finite Element Approach to Solving Large Deformation Problems in Multiphase Porous media. Ph.D. Dissertation, School of Civil Engineering and Environmental Science, University of Oklahoma, Norman, Oklahoma, December 2005. [Google Scholar]
  9. Gadala, M.S.; Oravas, G.A.E.; Dokainish, M.A. A consistent eulerian formulation of large deformation problems in statics and dynamics. Int. J. Non-Linear Mech. 1983, 18, 21–35. [Google Scholar] [CrossRef]
  10. Sanavia, L.; Schrefler, B.A.; Steinmann, P. A formulation for unsaturated porous medium undergoing large inelastic strains. Comput. Mech. 2002, 28, 137–151. [Google Scholar] [CrossRef]
  11. Gawin, D.; Sanavia, L. A Unified Approach to Numerical Modeling of Fully and Partially Saturated Porous Materials by Considering Air Dissolved in Water. Comput. Model. Eng. Sci. 2009, 53, 255–302. [Google Scholar] [CrossRef]
  12. Bannmann, D.J.; Johnson, G.C. On the kinematics of finite deformation plasticity. Acta Mech. 1987, 70, 1–13. [Google Scholar] [CrossRef]
  13. Lee, E.H. Some comments on elastic-plastic analysis. Int. J. Solids Struct. 1981, 17, 859–872. [Google Scholar] [CrossRef]
  14. Lubarda, V.A.; Lee, E.H. A correct definition of elastic and plastic deformation and its computational significance. J. Appl. Mech. 1981, 48, 35–40. [Google Scholar] [CrossRef]
  15. Oñate, E.; Carbonell, J.M. Updated lagrangian mixed finite element formulation for quasi and fully incompressible fluids. Comput. Mech. 2014, 54, 1583–1596. [Google Scholar] [CrossRef]
  16. Wang, S.; Chester, S.A. Multi-physics modeling and finite element formulation of corneal uv cross-linking. Biomech. Modeling Mechanobiol. 2021, 20, 1561–1578. [Google Scholar] [CrossRef] [PubMed]
  17. Wang, S.; Saito, K.; Kawasaki, H.; Holland, M.A. Orchestrated neuronal migration and cortical folding: A computational and experimental study. PLoS Comput. Biol. 2022, 18, e1010190. [Google Scholar] [CrossRef]
  18. Dafalias, Y.F.; Herrmann, L.R. Bounding Surface Plasticity II: Application to Isotropic Cohesive Soils. J. Eng. Mech. 1986, 112, 1260–1291. [Google Scholar] [CrossRef]
  19. Muraleetharan, K.K.; Nedunuri, P.R. A Bounding Surface Elastoplastic Constitutive Model for Monotonic and Cyclic Behavior of Unsaturated Soils. In Proceedings of the 12th Engineering Mechanics Conference, La Jolla, CA, USA, 17–20 May 1998; ASCE: Reston, VA, USA; pp. 1331–1334. [Google Scholar]
  20. Alonso, E.E.; Gens, A.; Josa, A.A. Constitutive Model for Partially Saturated Soils. Geotechnique 1990, 40, 405–430. [Google Scholar] [CrossRef]
  21. Wheeler, S.J.; Sivakumar, V. An Elasto-plastic Critical State Framework for Unsaturated Soil. Geotechnique 1995, 45, 35–53. [Google Scholar] [CrossRef]
  22. Wheeler, S.J. Inclusion of Specific Water Volume Within an Elasto-Plastic Model for Unsaturated Soil. Can. Geotech. J. 1996, 33, 42–57. [Google Scholar] [CrossRef]
  23. Ananthanathan, P. Laboratory Testing of Unsaturated Minco Silt. Master’s Thesis, University of Oklahoma, Norman, OK, USA, 2002. [Google Scholar]
  24. Vinayagam, T. Understanding the Stress-Strain Behavior of Unsaturated Minco Silt Using Laboratory Testing and Constitutive Modeling. Master’s Thesis, University of Oklahoma, Norman, OK, USA, 2004. [Google Scholar]
  25. Brooks, R.; Corey, A. Hydraulic Properties of Porous Media. In Hydrology Paper; Colorado State University: Fort Collins, Co, USA, 1964. [Google Scholar]
  26. Pepescu, R.; Prevost, J.H. Numerical Class A Prediction for Models Nos. 1, 2, 3, 4a, 4b, 6, 7, 11, &12. In Verification of Numerical Procedures for the Analysis of Soil Liquefaction Problems; Balkema: Rotterdam, The Netherlands, 1993; pp. 1105–1209. [Google Scholar]
  27. Lacy, S.J.; Prevost, J.H. Nonlinear Seismic Response Analysis of Earth Dams. J. Soil Dyn. Earthq. Eng. 1987, 6, 48–63. [Google Scholar] [CrossRef]
  28. Muraleetharan, K.K.; Mish, K.D.; Arulanandan, K. A Fully Coupled Nonlinear Dynamic Analysis Procedure and its Verification Using Centrifuge Test Results. Int. J. Numer. Anal. Methods Geomech. 1994, 18, 305–325. [Google Scholar] [CrossRef]
Figure 1. Motion of an REV and fluids in an unsaturated soil system.
Figure 1. Motion of an REV and fluids in an unsaturated soil system.
Geosciences 12 00320 g001
Figure 2. Multiplicative decomposition of elastic–plastic deformations.
Figure 2. Multiplicative decomposition of elastic–plastic deformations.
Geosciences 12 00320 g002
Figure 4. Finite element mesh and location of nodes and elements for records (all dimensions are in meters).
Figure 4. Finite element mesh and location of nodes and elements for records (all dimensions are in meters).
Geosciences 12 00320 g004
Figure 5. Schematic of a bounding-surface constitutive model in stress-invariant space.
Figure 5. Schematic of a bounding-surface constitutive model in stress-invariant space.
Geosciences 12 00320 g005
Figure 6. Gravity load–time history for static analysis.
Figure 6. Gravity load–time history for static analysis.
Geosciences 12 00320 g006
Figure 7. Input horizontal and vertical acceleration–time histories at the base of the model.
Figure 7. Input horizontal and vertical acceleration–time histories at the base of the model.
Geosciences 12 00320 g007
Figure 8. Horizontal and vertical displacement–time histories at node N053.
Figure 8. Horizontal and vertical displacement–time histories at node N053.
Geosciences 12 00320 g008
Figure 9. Horizontal and vertical displacement–time histories at node N073.
Figure 9. Horizontal and vertical displacement–time histories at node N073.
Geosciences 12 00320 g009
Figure 10. Fluid pressures, suction, and the degree of saturation–time histories at element E022.
Figure 10. Fluid pressures, suction, and the degree of saturation–time histories at element E022.
Geosciences 12 00320 g010aGeosciences 12 00320 g010b
Figure 11. Horizontal and vertical displacement–time histories at node N053—severe shaking.
Figure 11. Horizontal and vertical displacement–time histories at node N053—severe shaking.
Geosciences 12 00320 g011
Figure 12. Horizontal and vertical displacement–time histories at node N073—severe shaking.
Figure 12. Horizontal and vertical displacement–time histories at node N073—severe shaking.
Geosciences 12 00320 g012
Figure 13. Fluid pressures, suction, and the degree of saturation–time histories in element E022–severe shaking.
Figure 13. Fluid pressures, suction, and the degree of saturation–time histories in element E022–severe shaking.
Geosciences 12 00320 g013aGeosciences 12 00320 g013b
Table 1. Calibrated bounding surface model parameters.
Table 1. Calibrated bounding surface model parameters.
ParameterValue
Slope of virgin consolidation line in e-ln(p) plot, λ0.0954
Slope of swelling line in e-log(p) plot, κ0.0103
Slope of critical state line in compression, Mc1.2678
Ratio of extension to compression for slope of critical state line, Me/Mc1.00
Elastic shear modulus, G ( × 10 6 k P a ) 12.50
Poisson ratio, ν0.20
Mean net stress at which consolidation curve changes from linear in e-ln(p) space to linear in e-p space, PL33.80
Shape parameter in compression corresponding to ellipse 1, Rc4.20
Ratio of extension to compression of shape parameter corresponding to ellipse 1, Re/Rc1.00
Shape parameter in compression corresponding to hyperbola, Ac0.05
Ratio of extension to compression of shape parameter corresponding to hyperbola, Ae/Ac1.00
Shape parameter corresponding to ellipse 2, T0.01
Projection center, C0.00
Elastic zone parameter, Sp1.00
Positive model parameter, m0.02
Degree of hardening in triaxial compression, Hc0.80
Ratio of extension to compression for hardening parameter, He/Hc1.00
Hardening parameter for states in immediate vicinity of I-axis, Ho1.000
Suction-dependent parameter 1, m4.703
Suction-dependent parameter 2, N1.780
Suction-dependent parameter 3, A0.420
Suction-dependent parameter 4, B0.089
Table 2. Calibrated Brooks–Corey model parameters.
Table 2. Calibrated Brooks–Corey model parameters.
ParameterValue
Dry density (kN/m3)14.14
Parameter a5.5
Parameter n0.5
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Ravichandran, N.; Vickneswaran, T. Coupled Large Deformation Finite Element Formulations for the Dynamics of Unsaturated Soil and Their Application. Geosciences 2022, 12, 320. https://doi.org/10.3390/geosciences12090320

AMA Style

Ravichandran N, Vickneswaran T. Coupled Large Deformation Finite Element Formulations for the Dynamics of Unsaturated Soil and Their Application. Geosciences. 2022; 12(9):320. https://doi.org/10.3390/geosciences12090320

Chicago/Turabian Style

Ravichandran, Nadarajah, and Tharshikka Vickneswaran. 2022. "Coupled Large Deformation Finite Element Formulations for the Dynamics of Unsaturated Soil and Their Application" Geosciences 12, no. 9: 320. https://doi.org/10.3390/geosciences12090320

APA Style

Ravichandran, N., & Vickneswaran, T. (2022). Coupled Large Deformation Finite Element Formulations for the Dynamics of Unsaturated Soil and Their Application. Geosciences, 12(9), 320. https://doi.org/10.3390/geosciences12090320

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop