Next Article in Journal
Conceptual Design, Development, and Preliminary Safety Evaluation of a PWR Dry Storage Module for Spent Nuclear Fuel
Previous Article in Journal
The Influence of Exercise, Nutritional Status, and Disease on the Functional Ability to Undertake Activities of Daily Living and Instrumental Activities of Daily Living in Old Taiwanese People
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

An Explicit Finite Element Method for Saturated Soil Dynamic Problems and Its Application to Seismic Liquefaction Analysis

1
School of Civil Engineering, North China University of Technology, Beijing 100144, China
2
The Key Laboratory of Urban Security and Disaster Engineering of Ministry of Education, Beijing University of Technology, Beijing 100124, China
*
Author to whom correspondence should be addressed.
Appl. Sci. 2022, 12(9), 4586; https://doi.org/10.3390/app12094586
Submission received: 26 March 2022 / Revised: 16 April 2022 / Accepted: 21 April 2022 / Published: 30 April 2022
(This article belongs to the Topic Advances in Dynamics of Building Structures)

Abstract

:
An explicit finite element method is proposed to solve the u-p-formed dynamic equation of saturated porous media. In this method, a special discretization is implemented to discretize the global computational domain into local node systems of the individual nodes, and the dynamic equation of each local node system corresponding to each node is discretized by the completely explicit integration method in the time domain. By cycling through all the nodes, the dynamic responses of the global computational domain are obtained. In addition, a viscoelastic artificial boundary is added in the method. In numerical examples, the proposed FEM is verified by the good agreements between the results obtained from the proposed method and the analytical and numerical solutions of existing methods, respectively. After being embedded in OpenSees software, the proposed method is implemented into analysis of the seismic responses of liquefiable site.

1. Introduction

Fluid-saturated porous media is a widespread saturated soil, which is a combination of soil skeleton and fluid filled in the pores. The porous media theory and computational schemes have an important influence in geomechanics. Specifically, porous media analyses govern the physical behavior of cohesive soils consolidation and failure [1,2,3,4] and cohesionless soils [5,6] response when subjected to static and dynamic loading. The mechanical properties of the pore fluid and the fluid–solid couple interaction due to the relative movement between the pore fluid and the soil skeleton [7,8] will both directly affect the overall mechanical properties of the saturated porous media; in particular, the saturated sand is prone to liquefaction under dynamic loading [9]. Previous seismic investigations all indicated that ground failure and lateral large deformation due to soil liquefaction are the main causes of damage to the structure [10,11,12,13,14,15,16,17].
The macroscopic mechanical properties of saturated porous media are much more complex than single-phase media, so it is necessary to establish the equations to describe the mechanical behavior of the soil skeleton and pore water, respectively, and to consider the fluid–solid coupling interaction. In 1956, Biot [18] proposed a wave propagation theory of saturated porous media based on a quasi-static consolidation model, considering the inertial effects of the soil skeleton and pore fluid [19]. On that basis, Zienkiewicz et al. [20] rewrote the equations by using different simplifying assumptions, such as the u-U-p, u-U, u-p, and u-w models. Among them, the u-p formulation is a mixed equation of the vector displacement u and the scalar pore pressure p. Compared with the full vector formulation of the dynamic equation (e.g., u-U), the u-p equation has fewer degree of freedom and can be directly solved for the fluid pressure variables without the need of secondary calculations. It is suitable for describing the dynamic problems of saturated porous media under low-frequency loading.
Usually, the dynamic equations of saturated porous media can be solved by the analytical method only under a few special boundary conditions and load forms [21]. Therefore, numerical methods represented by finite element method (FEM) are the main methods for solving the dynamic problems of two-phase media. In the time domain, the integration methods can be categorized into implicit methods and explicit methods, such as the Newmark method [22], the Wilson-θ method [23], the generalized-α method [24,25], and the central difference method. The u-U and u-w formulation are often solved by using the staggered explicit–explicit [26,27], explicit–implicit [28,29] and implicit–implicit [28] methods. However, the implicit–implicit or explicit–implicit scheme is basically used to solve the u-p formulation of saturated porous media. This is because in the u-p formulation, the fluid phase equation is formally a diffusion equation in terms of pore pressure p and is usually solved by an implicit integration method. For example, Park [30] and Zienkiewicz et al. [31] established unconditionally stable staggered implicit–implicit solution procedures for the u-p dynamic equations of two-phase media. Huang and Zienkiewicz [32] improved the Zienkiewicz method by introducing a pore pressure correction method, which obtains both displacement and pore pressure through the estimated-corrected value calculation. Since the u-p formulation is a mixed form of dynamic equations, using the same order of shape functions of displacement and pore pressure would lead to instability of the results. To overcome this problem, Zienkiewicz et al. [33] divided the calculation of acceleration in one time step into two parts before spatial discretization, and proposed an explicit–implicit algorithm for the case of zero permeability assumption to effectively avoid the instability of the mixed formulations calculation. By introducing a parameter δ, Huang et al. [34] gave an implicit–implicit algorithm of the u-p formulation under zero permeability assumption, resulting in computational stability. Pastor et al. [35] proposed a fractional step algorithm based on velocity and pore pressure, allowing the shape functions of the same order for both displacement and pore pressure. To reduce computation and ensure computational stability in a relatively large time step, an iterative stabilized fractional algorithm [36] and an iterative pressure-stabilized fractional step algorithm [37] were proposed successively based on the algorithm of Pastor et al. Park and Tak [38] proposed an FEM-based stable algorithm that combines multiple time steps, remeshing, and subiteration. Soares et al. [39] divided the soil into solid and fluid subdomains and established an uncoupled model for each subdomain and solved each model separately. The dynamic equilibrium equation of the solid phase was calculated using the implicit Newmark method. For the liquid phase, a wave equation of pore pressure in scalar form was established by using its equilibrium equation and the fluid continuity equation, and it was calculated using the explicit Green–Newmark method. The coupling interaction between solid–fluid two phases, which had not yet been considered due to the use of uncoupled models, was achieved by applying interface forces on the two subdomains. This algorithm was finally proved unconditionally stable. By diagonalizing both mass and fluid compressibility matrices, Xu et al. [40] proposed a fully explicit iterative integration method with high efficiency. Considering the low computational efficiency of using constant time step in the temporal computation of the FEM, Tang et al. [41] used the mixed error of soil displacement and pore pressure as a criterion for adjusting the time step, and proposed an adaptive time-stepping method for calculating the fluid–solid coupling interaction or liquefaction problems based on the FEM-FDM method.
In this paper, an explicit FEM for solving the u-p equation of saturated porous media is proposed. In the space domain, a spatial discretization is employed that divides the whole computational domain into local node systems for all nodes, i.e., the dynamic responses of any node are related only to its surrounding nodes. In the time domain, combined with a fully explicit integration method for the u-p equation [40], the displacement and pore pressure responses of the soil in any time step are directly calculated using the results from the previous two time steps. In addition, the proposed method is combined with the viscoelastic local artificial boundary [42] to establish an effective FEM with decoupling characteristics for solving the wave problem of saturated porous media. In the end, two examples are verified the proposed method, and the method is embedded into OpenSees software to analyze the seismic liquefaction of the saturated sand.

2. Wave Equations of Saturated Porous Media in u-p Form

Zienkiewicz et al. [20] proposed the u-p formulation for fluid-saturated porous media, with the following governing equations:
L T D L u α p ρ u ¨ = 0
T u ˙ T k ¯ p + 1 Q b p ˙ = 0
where u is the displacement of solid-phase soil skeleton, p is the pore pressure, α is the compressibility coefficient of soil skeleton, ρ is the density of two-phase media, k ¯ is the dynamic permeability coefficient, and Q b is the compressibility coefficient of pore fluid. In two-dimensional case, T = [ x y ] ; D = [ λ + 2 G λ 0 λ λ + 2 G 0 0 0 G ] , and L T = [ x 0 y 0 y x ] , where λ and G are the Lame constants of the soil skeleton.

3. Explicit Finite Element Method for the u-p Formulation

Taking the two-dimensional (2D) plane strain problem of saturated porous media as an example, the entire computational domain Ω is discretized into a plane domain represented by 2D four-node elements. The shape functions of u, p, fu, and fq at any point in element e are defined as follows:
{ u e = N u u e ;             p e = N p p e f u e = N u f u e ;         f q e = N p f q e
where N u and N p are the shape functions of displacement and pore pressure in element e, as shown in Equation (4), u e and p e are the displacements and pore pressures inside the element e, respectively, f u e and f q e are the external forces and flow rates inside the element e, and u e , p e , f u e , and f q e are the nodal displacements, pore pressures, external forces, and flow rates, respectively, in element e, as shown in Equation (5).
{ N u = [ N u 1 N u 2 N u 3 N u 4 ] ;             N u i = [ N u i x 0 0 N u i y ] ,       ( i = 1 , 2 , 3 , 4 ) N p = [ N p 1 N p 2 N p 3 N p 4 ]
{ u e T = [ u 1 T , u 2 T , u 3 T , u 4 T ] ; u i T = [ u i x , u i y ] , ( i = 1 , 2 , 3 , 4 )   p e T = [ p 1 T , p 2 T , p 3 T , p 4 T ]       f u e T = [ f u 1 T , f u 2 T , f u 3 T , f u 4 T ] ;   f u i T = [ f u x , f u y ] , ( i = 1 , 2 , 3 , 4 )         f q e T = [ f q 1 T , f q 2 T , f q 3 T , f q 4 T ]      
Finite element discretization of Equations (1) and (2) is performed using the Galerkin method [43] to obtain the following Galerkin weak form of the u-p equation:
e Ω e [ N u T L T D L N u u e N u T N p p e + ρ N u T N u u ¨ e ] d Ω = e Γ e N u T f u e d Γ
e Ω e [ ( N p ) T N u u ˙ e +   k ¯ ( N p ) T ( N p ) p e + 1 Q N p T N p p ˙ e ] d Ω = e Γ e N p T f q e d Γ
where Ω e is the computational domain of any element and Γ e is the domain boundary of this element.

3.1. Spatial Discretization

The dynamic responses of node k only are associated with the responses of the neighboring nodes. On that basis, taking the 2D plane strain problem as an example, the whole computational domain is discretized into local node systems for all nodes [44], as shown in Figure 1. For node k1 at the corner of the computational domain, the four nodes of element 1 adjacent to node k are renumbered as nodes 1 to 4; similarly, for the node k2 on the boundary and the node k3 located inside the computational domain, the local node systems are renumbered by the method shown in Figure 1. In the computational domain, for the local node system composed of any node, an individual finite element equation of the local node system is constructed and is computed in time domain. The dynamic responses for the whole computational domain are obtained by cycling all the nodes.
According to Figure 1, it is assumed that there are n adjacent elements for node k, and each element has four nodes. Based on Equations (6) and (7), the finite element equations of the local node system of node k are expressed as follows:
L = 1 n j = 1 4 M k j L u ¨ j L + L = 1 n j = 1 4 K k j L u j L L = 1 n j = 1 4 Q k j L p j L = L = 1 n f k , u L
L = 1 n j = 1 4 S k j L p ˙ j L + L = 1 n j = 1 4 H k j L p j L + L = 1 n j = 1 4 ( Q k j L ) T u ˙ j L = L = 1 n f k , q L
where L represents the influence of the adjacent element L on node k, and j = 1–4 represent the influence of the four nodes of the element L on node k. The matrices in the above equations are defined as follows:
M i j L = 1 1 1 1 ρ N u i N u j | J | d ξ d η K i j L = 1 1 1 1 N u i T L T D L N u j | J | d ξ d η Q i j L = 1 1 1 1 N u i T N p j | J | d ξ d η H i j L = 1 1 1 1 ( N p i ) T   k ¯ ( N p j ) | J | d ξ d η S i j L = 1 1 1 1 N p i T 1 Q b N p j | J | d ξ d η f i , u L = Γ L N u T T s d r f i , q L = Γ L N p T q s d r
where M is the mass matrix of saturated porous media, K is the stiffness matrix of the solid phase, Q is the coupling matrix, S is the compressibility matrix of the fluid phase, H is the permeability matrix of the fluid phase, fu and fq are the external load of and the flow rate, respectively, Ts is the distributed force on the boundary of element L in the x and y directions, and qs is the flow rate perpendicular to the boundary of element L.
In fact, Equations (8) and (9) are an alternative expression of Equations (6) and (7). Equations (6) and (7) are the global equations established by assembling the coefficient matrices of the global computational domain. In comparison, Equations (8) and (9) are the equations established for the local node system. The dynamic responses of the whole domain are obtained by cycling all nodes in the computational domain and solving the equation of local node system of each node. Therefore, for the dynamic problem with large degree of freedom, the dimension of equation of the local node system is very small and computational efficiency is higher.

3.2. Explicit Integration Method in Time Domain

To increase the computational efficiency, the finite element Equations (8) and (9) are computed using the explicit integration method [40] in time domain. The equations of the local node system of node k at the t + 1 time step can be expressed as:
u k t + 1 = 2 u k t u k t 1 + Δ t 2 1 M d k L = 1 n f k , u L , t + Δ t 2 1 M d k L = 1 n j = 1 4 Q k j L p j L , t Δ t 2 1 M d k L = 1 n j = 1 4 K k j L u j L , t
p k t + 1 = p k t + Δ t 1 S d k L = 1 n f k , q L , t 1 S d k L = 1 n j = 1 4 ( Q k j L ) T ( u j t u j t 1 ) Δ t 1 S d k L = 1 n j = 1 4 H k j L p j L , t
where the lumped mass matrix M d k e and the lumped fluid compressibility matrix S d k e of element e are processed by the same diagonalization method [40], as shown in Equations (13) and (14).
( M d k e ) i j = { a ( M k e ) i i = a V e ρ e N i T N i d V e             ( j = i )                                           0                                               ( j i )       ,                   i = 1 n e ( M d k e ) i i = a i = 1 n e ( M k e ) i i = W e I d = ρ e A e I d
( S d k e ) i j = { b ( S k e ) i i = b V e 1 Q b e N i T N i d V e         ( j = i )                                           0                                               ( j i )       ,                   i = 1 n e ( S d k e ) i i = b i = 1 n e ( S k e ) i i = V e I d = 1 Q b A e I d
where M k e and S k e are the consistent mass matrix and the consistent fluid compressibility matrix of element e, a and b are corresponding scaling coefficients, We and Ve are the mass and compressive deformation of element e, respectively, Ae is the area of element e, and Id is unit matrix.

4. Implementation of the Explicit Finite Element Method in OpenSees

The FEM been embedded in the open-source OpenSees software. The calculation does not require solving simultaneous equations but rather solving three uncoupled equations corresponding to the three DOFs (ukx, uky, p) of the node. The dynamic equations AX = B [45,46] can be computed in OpenSees by selecting a suitable algorithm and integrator. The explicit temporal integration method UPExplicitdifference, the algorithm UPExplicitNR, and the equation solver ExplicitLinSOESolver are also added, respectively, and FourNodeQuadUP element for diagonalizing the mass matrix M and fluid compressibility matrix S of elements are added in the software.
The equivalent equations of Equations (11) and (12) are expressed as follows:
K ˜ k   v ˜ k t + 1     = F ˜ k
where
F ˜ k = F k R k M ¯ k v ˜ ¨ t   C ¯ k v ˜ ˙ t  
K ˜ k = [ M d k 0 0 S d k ] , C ¯ k = [ 0 Q k Q k T H k ] , F k = { f u k t f q k t }
v ˜ k t + 1     = { u k t + 1     Δ t 2 p k t + 1 Δ t } , v ˜ ˙ t   = { u t u t 1 Δ t p t } , v ˜ ¨ t   = { u t 1 2 u t Δ t 2 p t Δ t }
where K ˜ in Equation (15) is a diagonal matrix. The dynamic responses of node k are calculated by only the three equations corresponding to the three DOFs (uxk, uyk, pk) in Equation (15):
{ M d k , x   v ˜ k , x t + 1     = F ˜ k , x M d k , y   v ˜ k , y t + 1       = F ˜ k , y S d k     v ˜ k , p t + 1     = F ˜ k , p
Equation (19) is the equivalent form to Equations (11) and (12). The displacements and pore pressure of the node k at time t + 1 are obtained by Equation (19).
{ u k , x t + 1       u k , y t + 1       p k t + 1 } = { Δ t 2 v ˜ k , x t + 1     Δ t 2 v ˜ k , y t + 1       Δ t v ˜ k , p t + 1     }
The dynamic responses of the global domain at any time can be obtained by solving of Equation (19) of all nodes. The computational flowchart of the proposed explicit FEM in OpenSees is shown in Figure 2. In space domain, the proposed algorithm UPExplicitNR is used to construct individual finite element equation of the local node system, and results of the whole system are obtained by cycling all the nodes based on UPExplicitNR. In time domain, the explicit temporal integration method UPExplicitdifference is used to compute the results in every time step.
The proposed algorithm has a much lower computational efficiency than implicit methods for problems with a small DOFs. However, for the problems with a large DOFs or strong nonlinear problems, the proposed method does not require the solving of simultaneous equations, only computing three decoupled equations for each node, and no nonlinear iterative computation is performed within a suitable time step. Therefore, the efficiency of the proposed method is improved significantly, and the computational memory is saved.

5. Viscoelastic Artificial Boundary

For the problems of wave propagation in saturated porous media, a viscoelastic artificial boundary for the u-p equations is adopted. It is expressed as follows:
f u B l = K B l u B l C B l u ˙ B l Q B l p B l
f q B l = H B l p B l + S B l p ˙ B l
where ∞ denotes the exterior domain, Bl is the node l on the boundary, f u B l and f q B l are the stress and flow velocity of the node l on the boundary, respectively, and u B l and p B l are the displacement and pore pressure of the node l on the boundary, respectively. The coefficients in the above equations are expressed as follows:
K B l = W T L l [ K B N 0 0 K B T ] W ;   C B l = W T L l [ C B N 0 0 C B T ] W ;   Q B l = W T L l [ Q B 0 ]
H B l = L l H B S B l = L l S B
where Ll is the length corresponding to node l on the boundary, and W is the coordinate transformation matrix for transforming the local coordinates to the global coordinates. The parameters in the above equations are specifically
W = [ l r x l r y l x l y ]
K B N = λ + 2 G 2 r ( 1 + A ) ;   K B T = G 2 r ( 1 + A )
C B N = B ρ C P 1 r ;   C B T = B ρ C S r
Q B = α ;   H B = k ¯ 2 r ( 1 + A ) ;   S B = k ¯ C P 1 r
In Equation (25), l r x = cos ( r , x ) is the cosine of the angle between the positive direction of the local r axis and the positive direction of the global x-axis. Other parameters in Equation (25) are defined similarly to lxy. A = 0.8 and B = 1.1.
The artificial boundary condition as shown in Equations (21) and (22) actually changes only the diagonal values of the coefficient matrix corresponding to boundary nodes in local node system. Assume that there are m nodes on the boundary in the local node system of node k. Equations (21) and (22) are substituted into Equations (8) and (9) to obtain the finite element equations considering the viscoelastic artificial boundary:
M d k u ¨ k + l = 1 m C B l u ˙ B l + ( L = 1 n j = 1 4 K k j L u j L + l = 1 m K B l u B l ) ( L = 1 n j = 1 4 Q k j L p j L l = 1 m Q B l p B l ) = L = 1 n f k , u L
( S d k p ˙ k + l = 1 m S B l p ˙ B l ) + ( L = 1 n j = 1 4 H k j L p j L + l = 1 m H B l p B l ) + L = 1 n j = 1 4 ( Q k j L ) T u ˙ j L = L = 1 n f k , q L
Compared to Equation (8), Equation (29) has an additional damping term and can be solved by an explicit integration method [47], and the equations of the local node system of node k at the t + 1 time step can be expressed as:
u k t + 1 = u k t + Δ t 2 2 1 M d k L = 1 n f k , u L , t + Δ t 2 2 1 M d k ( L = 1 n j = 1 4 Q k j L p j L l = 1 m Q B l p B l ) Δ t 2 2 1 M d k ( L = 1 n j = 1 4 K k j L u j L + l = 1 m K B l u B l ) + Δ t 2 1 M d k l = 1 m C B l u B l t Δ t 2 1 M d k l = 1 m C B l u B l t 1 Δ t 2 1 M d k l = 1 m C B l u ˙ B l t + Δ t u ˙ B k t
p k t + 1 = p k t + Δ t S d k + S B k L = 1 n f k , q L , t 1 S d + S B k L = 1 n j = 1 4 ( Q k j L ) T ( u j t u j t 1 ) Δ t S d + S B k ( L = 1 n j = 1 4 H k j L p j L + l = 1 m H B l p B l )
where, the subscript Bk of u ˙ B k t and S B k represents the corresponding value when the target node k is also the node on the boundary.

6. Validation of Method and Comparison of Computational Efficiency

Example 1.
The computational model for a one-dimensional (1D) semi-infinite domain of saturated porous media is shown in Figure 3. The left and right sides of the model are impermeable boundaries, and their horizontal displacement is fixed. The bottom of the model is assigned a viscoelastic artificial boundary, and the top surface of the model is a free drainage boundary applied with a constant uniform load of 1.0 Pa. The material parameters of the model listed in Table 1 [21,48] as Material No. 1 are taken from the reference [21,48]. The analytical solutions, as the reference solutions, for the 1D problem of saturated porous media proposed by Simon et al. [21].
In Figure 4, the numerical solutions and analytical solutions are compared at three dimensionless time parameter τ = 20, 40, and 60, respectively. The numerical solutions of the displacement and pore pressure of the 1D soil column are in relatively good agreement with the analytical solutions. In addition, the numerical and analytical solutions of the displacement and pore pressure dynamic time-history responses at nodes 5 m and 45 m from the center of the load are compared in Figure 4, which shows consistency between the two results. The comparison results verify the accuracy of proposed method.
Example 2.
Computational efficiency is compared in this example. A 2D computational model with size of 50 m × 50 m is shown in Figure 5. The top surface of the model is set as a drainage boundary and is subjected to a uniformly distributed sinusoidal load, and each of the other three sides are set as artificial boundaries. Three target nodes (1, 2, and 3) are selected for analysis. The material properties are listed in Table 1 as Material No. 2. Xu’s method and Newmark method are used to compare the responses and computational efficiency with the proposed method, respectively.
Figure 6 shows the displacement and pore pressure time-history curves of the target nodes. The red and blue curves in the figure represent the results obtained from the proposed method and the Xu’s method [40], respectively. The results of the two methods have a good consistence. The consistency further validates the proposed method.
By refining grids of the model in Figure 5, the total computational times of the proposed method and Newmark [31] method are compared in Figure 7. The model is discretized into approximately 1000, 2000, 3000, 4000, and 5000 nodes, respectively. Comparison of the computation time of the two algorithms show that when the model has a relatively small number of nodes, the implicit Newmark method leads to much less computation time than the proposed method, because it is unconditionally stable and allows a relatively large time step. However, as the number of nodes increases to approximately 3000, the total computation time of the two methods is basically similar, and the computational efficiency of the proposed method is significantly improved. When the number of nodes increases to approximately 4000, the computational efficiency of the proposed explicit method is far higher than the implicit Newmark method, and such an advantage becomes more significant as the number of nodes continues to increase.

7. Analysis of Nonlinear Problems

7.1. Finite Element Model

Figure 8 shows the 2D model of a liquefiable site, which is composed of a 15 m dense sand layer, a 15 m liquefiable loose sand layer, and an 8 m-dense sand embankment from bottom to top. The seawater has a depth of 2 m. The top surface of the model is set as a drainage boundary. To simplify the model, the horizontal and vertical displacements at the bottom of the model are fixed, and the horizontal displacements of the left and right side boundaries are also fixed.
The Wolong station seismic record of the Wenchuan earthquake is selected as the input motion. It has a peak acceleration of 0.3 g. The acceleration time history curve is shown in Figure 9. Attention is fixed on the numerical results at 6.2 s, 16.6 s, and 30 s, respectively.
The saturated soils are simulated using the quadrilateral plane strain element element quadUP [49,50]. The parameters of soil elements are listed in Table 2. Each node of the element has three DOFs, corresponding to horizontal displacement, vertical displacement, and pore pressure, respectively. In view of using a 2D plane strain model to simulate a 3D problem, the out of plane thickness of the model is set to be 1.0 m. The constitutive model of saturated soil chooses the PressureDependMutiYield02 (PDMY) command [51,52], with the material parameters [49] listed in Table 3.

7.2. Numerical Results and Analysis

The seismic responses of the saturated porous media are analyzed by the proposed method and the implicit Newmark method, respectively. Figure 10 shows the time history curves of horizontal and vertical displacements at measuring points D1–D4. The horizontal displacement and vertical displacement at each measuring point gradually increase with time and gradually decrease with burial depth, and reach the maximum of 0.1 m and 0.12 m, respectively, at the top of the embankment. In addition, the time history curves of both horizontal and vertical displacement show the stepwise increase in amplitude, respectively, at 6.2 s and 16.6 s.
Figure 11 shows the excess pore pressure time history curves at measuring points W1–W4. These curves exhibit basically the same variation trends, and their peaks occur at basically consistent time points. The amplitude of pore pressure at each measuring point gradually increases with the increase in the burial depth. The excess pore pressure of the dense sand layer reaches a maximum value of 75 kPa.
Moreover, Figure 10 and Figure 11 compare the results of the proposed method and the implicit Newmark method. The comparison shows that the numerical results of each measuring point obtained by the two methods agree well, indicating the correctness of the proposed method.
To observe the overall lateral deformation and settlement of the site, Figure 12 shows the contour plots of horizontal displacement of the site at 6.2 s, 16.6 s, and 30.0 s, as well as the contour plots of vertical displacement at 30 s. The deformation in the figure is magnified 30 times. As can be seen from the figure, the horizontal displacement at the nodes on the left and right sides of the embankment is the largest; the embankment gradually collapses under the earthquake, causing the embankment soil to move toward both sides, with the maximum horizontal and vertical displacements reaching 0.13 m and 0.19 m, respectively, at 30 s.
Figure 13 shows the contour plots of pore pressure ratios of the ground soil. The loose sand layer near the drainage boundary liquefies at 6.2 s and 16.6 s. As the decrease in the input acceleration, the pore water is drained upwards, and the soil particles are redistributed. The soil tends to a stable state.
Figure 12 and Figure 13 also compare the contour plots of displacement and pore pressure ratio computed by the proposed method and the implicit Newmark implicit method, respectively. The soil deformation patterns obtained by the two methods are similar, and the magnitudes of displacement and pore pressure computed by the two methods are basically consistent, further validating the correctness of the proposed method.

8. Conclusions

An explicit FEM for solving the dynamic responses of saturated porous media is proposed based on the u-p formulation. In spatial domain, the global system is divided into the local node systems of all nodes, which means that the responses of each node only are related to the responses of the surrounding nodes. In time domain, the local node system of each node can be discretized by a completed explicit integration method, and further be simplify into three independent equations corresponding to each DOF of the node. The dimension of FE equations of each local node system is small and is solved fast. By cycling through all the nodes, the dynamic responses of the global computational domain can be obtained. In this way, the proposed method has the high efficiency for solving larger DOFs problem and strong nonlinear problem of the fluid-saturated porous media. Moreover, a viscoelastic artificial boundary is added to the method, and changes only the diagonal values of the coefficient matrix corresponding to boundary nodes and does not change the form of each coefficient matrix corresponding to each local node system. Additionally, the proposed method is validated by comparing with the analytical and numerical solutions of existing methods, and it has high computational efficiency. In the end, the application to the seismic responses of liquefiable site indicates that the proposed method is an effective method for solving the dynamic problems of saturated porous media.
In this method, a complete explicit integration method with only first-order accuracy is chosen. However, based on the same spatial discretization, other suitable integration methods with high order accuracy can also be used in the FE method, such as the integration method used in the viscoelastic artificial boundary section. In addition, the computational efficiency of the proposed method only is verified by comparing with implicit Newmark method. The generality of the proposed needs to be confirmed for various types of analyses and comparisons.

Author Contributions

Software, C.F.; Validation, F.W.; Writing—original draft, J.S.; Writing—review & editing, C.X. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by National Natural Science Foundation of China (Grant No. 51808006 and No. 51722801); Natural Science Foundation of Beijing, China (Grant No. 8192012) and Beijing Municipal Education Commission Scientific Research Project (110052972027/146).

Acknowledgments

The authors would like to thank the National Natural Science Foundation of China (Grant No. 51808006 and No. 51722801) and the Natural Science Foundation of Beijing, China (Grant No. 8192012) for funding the work presented in this paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Houmadi, Y.; Benmoussa, M.Y.C.; Cherifi, W.N.E.H.; Rahal, D.D. Probabilistic analysis of consolidation problems using subset simulation. Comput. Geotech. 2020, 124, 103612. [Google Scholar] [CrossRef]
  2. Savvides, A.A.; Papadrakakis, M. A probabilistic assessment for porous consolidation of clays. SN Appl. Sci. 2020, 2, 2115. [Google Scholar] [CrossRef]
  3. Ahmed, A.; Soubra, A.H. Probabilistic analysis of strip footings resting on a spatially random soil using subset simulation approach. Georisk Assess. Manag. Risk Eng. Syst. Geohazards 2012, 6, 188–201. [Google Scholar] [CrossRef] [Green Version]
  4. Savvides, A.A.; Papadrakakis, M. A computational study on the uncertainty quantification of failure of clays with a modified Cam-Clay yield criterion. SN Appl. Sci. 2021, 3, 659. [Google Scholar] [CrossRef]
  5. Mercado, V.; Ochoa-Cornejo, F.; Astroza, R.; El-Sekelly, W.; Abdoun, T.; Pastén, C.; Pastén, F. Uncertainty quantification and propagation in the modeling of liquefiable sands. Soil Dyn. Earthq. Eng. 2019, 123, 217–229. [Google Scholar] [CrossRef]
  6. Najjar, S.S.; Sadek, S.; Alcovero, A. Quantification of model uncertainty in shear strength predictions for fiber-reinforced sand. J. Geotech. Geoenviron. Eng. 2013, 139, 116–133. [Google Scholar] [CrossRef]
  7. Zienkiewicz, O.C.; Chan, A.H.C.; Pastor, M.; Paul, D.K.; Shiomi, T. Static and dynamic behaviour of soils: A rational approach to quantitative solutions. I. Fully saturated problems. Proc. R. Soc. Lond. A Math. Phys. Eng. Sci. R. Soc. 1990, 429, 285–309. [Google Scholar]
  8. Zienkiewicz, O.C.; Taylor, R.L. The Finite Element Method; McGraw-Hill: London, UK, 1977. [Google Scholar]
  9. Marcuson, W.F. Definition of terms related to liquefaction. J. Geotech. Eng. Div. 1978, 104, 1197–1200. [Google Scholar]
  10. Tang, L. Study on p-y Curve Models for Pile-Soil Dynamic Interactions in Liquefied Ground; Harbin Institute of Technology: Harbin, China, 2010. (In Chinese) [Google Scholar]
  11. Mori, S.; Namuta, A.; Miwa, S. Feature of liquefaction damage during the 1993 Hokkaido Nanseioki earthquake. In Proceedings of the 29th Annual Conference of Japanese Society of Soil Mechanics and Foundation Engineering, New Delhi, India, June 1994; pp. 1005–1008. [Google Scholar]
  12. Cao, Z.; Hou, L.; Yuan, X.; Sun, R.; Wang, W.; Chen, L. Characteristics of liquefaction-induced damages during Wenchuan Ms 8.0 earthquake. Rock Soil Mech. 2010, 31, 3549–3555, (In Chinese with English Abstract). [Google Scholar]
  13. Cao, Z.; Youd, T.L.; Yuan, X. Gravelly soils that liquefied during 2008 Wenchuan, China earthquake, Ms = 8.0. Soil Dyn. Earthq. Eng. 2011, 31, 1132–1143. [Google Scholar] [CrossRef]
  14. Huang, Y.; Jiang, X.M. Field-observed phenomena of seismic liquefaction and subsidence during the 2008 Wenchuan earthquake. Nat. Hazards 2010, 54, 839–850. [Google Scholar] [CrossRef]
  15. Ishihara, K.; Yasuda, S.; Yoshida, Y. Liquefaction-induced flow failure of embankment and residual strength of silty sand. Soils Found 1990, 30, 69–80. [Google Scholar] [CrossRef] [Green Version]
  16. Hamada, M.; Sato, H.; Kawakami, T. A consideration of the mechanism for liquefaction-related large ground displacement. In Proceedings of the Fifth U.S.–Japan Workshop on Earthquake Resistant Design of Lifeline Facilities and Countermeasures Against Soil Liquefaction, Salt Lake City, UT, USA, 29 September–1 October 1994; pp. 217–232. [Google Scholar]
  17. Finn, W.D.L.; Fujita, N. Piles in liquefiable soils: Seismic analysis and design issues. Soil Dyn. Earthq. Eng. 2002, 22, 731–742. [Google Scholar] [CrossRef]
  18. Biot, M.A. General solutions of the Equations of elasticity and consolidation for a porous material. Appl. Mech. 1956, 23, 91–96. [Google Scholar] [CrossRef]
  19. Biot, M.A. Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range. Acoust. Soc. Am. 1956, 28, 168–178. [Google Scholar] [CrossRef]
  20. Zienkiewicz, O.C.; Chang, C.T.; Bettess, P. Drained, undrained, consolidating, and dynamic behavior assumptions in soils: Limits of validity. Geotechnique 1980, 30, 385–395. [Google Scholar] [CrossRef]
  21. Simon, B.R.; Zienkiewicz, O.C.; Paul, D.K. An analytical solution for the transient response of saturated porous elastic solids. Numer. Anal. Methods Geomech. 1984, 8, 381–398. [Google Scholar] [CrossRef]
  22. Newmark, N.M. A method of computation for structural dynamics. J. Eng. Mech. Div. 1959, 85, 67–94. [Google Scholar] [CrossRef]
  23. Wilson, E.L.; Farhoomand, I.; Bathe, K.J. Nonlinear dynamic analysis of complex structures. Earthq. Eng. Struct. Dyn. 1972, 1, 241–252. [Google Scholar] [CrossRef]
  24. Chung, J.; Hulbert, G.M. A time integration algorithm for structural dynamics with improved numerical dissipation. Gen.-Method J. Appl. Mech. 1993, 60, 371–375. [Google Scholar] [CrossRef]
  25. Kontoe, S.; Zdravkovic, L.; Potts, D.M. An assessment of time integration schemes for dynamic geotechnical problems. Comput. Geotech. 2008, 35, 253–264. [Google Scholar] [CrossRef] [Green Version]
  26. Zhao, C.; Li, W.; Wang, J. An explicit finite element method for Biot dynamic formulation in fluid-saturated porous media and its application to a rigid foundation. J. Sound Vib. 2005, 282, 1169–1181. [Google Scholar] [CrossRef]
  27. Wang, J.; Zhang, C.; Du, X. An explicit integration scheme for solving dynamic problems of solid and porous media. J. Earthq. Eng. 2008, 12, 293–311. [Google Scholar] [CrossRef]
  28. Simon, B.R.; Wu, J.S.S.; Zienkiewicz, O.C.; Paul, D.K. Evaluation of u-w and u-π finite element methods for the dynamic response of saturated porous media using one-dimensional models. Int. J. Numer. Anal. Methods Geomech. 1986, 10, 461–482. [Google Scholar] [CrossRef]
  29. Prevost, J.H. Wave propagation in fluid-saturated porous media: An efficient finite element procedure. Soil Dyn. Earthq. Eng. 1985, 4, 183–201. [Google Scholar] [CrossRef]
  30. Park, K.C. Stabilization of partitioned solution procedure for pore fluid-soil interaction analysis. Int. J. Numer. Methods Eng. 1983, 19, 1669–1673. [Google Scholar] [CrossRef]
  31. Zienkiewicz, O.C.; Paul, D.K.; Chan, A.H.C. Unconditionally stable staggered solution procedure for soil-pore fluid interaction problems. Int. J. Numer. Methods Eng. 1988, 26, 1039–1055. [Google Scholar] [CrossRef]
  32. Huang, M.S.; Zienkiewicz, O.C. New unconditionally stable staggered solution procedures for coupled soil-pore fluid dynamic problems. Int. J. Numer. Methods Eng. 1998, 45, 1029–1052. [Google Scholar] [CrossRef]
  33. Zienkiewicz, O.C.; Huang, M.S.; Wu, J.; Wu, S.M. A new algorithm for the coupled soil-pore fluid problem. Shock. Vib. 1993, 1, 3–14. [Google Scholar] [CrossRef]
  34. Huang, M.S.; Wu, S.M.; Zienkiewicz, O.C. Incompressible or nearly incompressible soil dynamic behavior—A new staggered algorithm to circumvent restrictions of mixed formulation. Soil Dyn. Earthq. Eng. 2001, 21, 169–179. [Google Scholar] [CrossRef]
  35. Pastor, M.; Li, T.; Liu, X.; Zienkiewicz, O.C.; Quecedo, M. A fractional step algorithm allowing equal order of interpolation for coupled analysis of saturated soil problems. Mech. Cohesive-Frict. Mater. 2000, 5, 511–534. [Google Scholar] [CrossRef]
  36. Li, X.K.; Han, X.H.; Pastor, M. An iterative stabilized fractional step algorithm for finite element analysis soil dynamics. Comput. Methods Appl. Mech. Eng. 2003, 192, 3845–3859. [Google Scholar] [CrossRef]
  37. Li, X.K.; Zhang, X.; Han, X.L.; Sheng, D.C. An iterative pressure-stabilized fractional step algorithm in saturated soil dynamics. Int. J. Numer. Anal. Methods Geomech. 2010, 34, 733–753. [Google Scholar] [CrossRef]
  38. Teahyo, P.; Moonho, T. A new coupled analysis for nearly incompressible and inpermeable saturated porous media on mixed finite element method: I. Proposed method. KSCE J. Civ. Eng. 2010, 14, 7–16. [Google Scholar]
  39. Soares, D.; Rodrigues, G.G.; Gonçalves, K.A. An efficient multi-time-step implicit-explicit method to analyze solid-fluid coupled systems discretized by unconditionally stable time-domain finite element procedures. Comput. Struct. 2010, 88, 387–394. [Google Scholar] [CrossRef]
  40. Xu, C.S.; Song, J.; Du, X.L.; Zhong, Z.L. A completely explicit finite element method for solving dynamic u-p Equations of fluid-saturated porous media. Soil Dyn. Earthq. Eng. 2017, 97, 364–376. [Google Scholar] [CrossRef]
  41. Tang, X.W.; Zhang, X.W.; Uzuoka, R. Novel adaptive time stepping method and its application to soil seismic liquefaction analysis. Soil Dyn. Earthq. Eng. 2015, 71, 100–113. [Google Scholar] [CrossRef]
  42. Xu, C.S.; Song, J.; Du, X.L.; Zhao, M. A local artificial-boundary condition for simulating transient wave radiation in fluid-saturated porous media of infinite domains. Int. J. Numer. Methods Eng. 2017, 112, 529–552. [Google Scholar] [CrossRef]
  43. Zienkiewicz, O.C.; Chan, A.H.C.; Pastor, M.; Schrefler, B.A.; Shiomi, T. Computational Geomechanics; Wiley: Chichester, UK, 1999. [Google Scholar]
  44. Zhao, C.; Li, W.; Wang, J. An explicit finite element method for dynamic analysis in fluid saturated porous medium-elastic single-phase medium-ideal fluid medium coupled systems and its application. J. Sound Vib. 2005, 282, 1155–1168. [Google Scholar] [CrossRef]
  45. Mazzoni, S.; McKenna, F.; Scott, M.H.; Fenves, G.L. OpenSees Command Language Manual. Pacific Earthquake Engineering Research Center; University of California: Berkeley, CA, USA, 2007. [Google Scholar]
  46. McKenna, F.T. Object-Oriented Finite Element Programming: Frameworks for Analysis, Algorithms and Parallel Computing; University of California: Berkeley, CA, USA, 1997. [Google Scholar]
  47. Du, X.L.; Wang, J.T. An explicit difference formulation of dynamic response calculation of elastic structure with damping. Eng. Mech. 2000, 17, 37–43. [Google Scholar]
  48. Akiyoshi, T.; Fuchida, K.; Fang, H.L. Absorbing boundary conditions for dynamic analysis of fluid-saturated porous media. Soil Dyn. Earthq. Eng. 1994, 13, 387–397. [Google Scholar] [CrossRef]
  49. Yang, Z.; Lu, J.; Elgamal, A. OpenSees Soil Models and Solid-Fluid Fully Coupled Elements User’s Manual; University of California: San Diego, CA, USA, 2008. [Google Scholar]
  50. Yang, Z.H.; Elgamal, A. Command Manual and User Reference for OpenSees Soil Models and Fully Coupled Element Developed; University of California: San Diego, CA, USA, 2003. [Google Scholar]
  51. Elgamal, A.; Yang, Z.H.; Parra, E.; Ragheb, A. Modeling of cyclic mobility in saturated cohesionless soils. Int. J. Plast. 2003, 19, 883–905. [Google Scholar] [CrossRef]
  52. Yang, Z.; Elgamal, A. Numerical Modeling of Earthquake Site Response Including Dilation and Liquefaction; Columbia University: New York, NY, USA, 2000. [Google Scholar]
Figure 1. Local node systems.
Figure 1. Local node systems.
Applsci 12 04586 g001
Figure 2. Computational flowchart of the explicit method.
Figure 2. Computational flowchart of the explicit method.
Applsci 12 04586 g002
Figure 3. Numerical model for 1D saturated porous media.
Figure 3. Numerical model for 1D saturated porous media.
Applsci 12 04586 g003
Figure 4. Comparison of numerical and analytical solutions for the 1D model of saturated porous media.
Figure 4. Comparison of numerical and analytical solutions for the 1D model of saturated porous media.
Applsci 12 04586 g004
Figure 5. Numerical model for 2D saturated porous media.
Figure 5. Numerical model for 2D saturated porous media.
Applsci 12 04586 g005
Figure 6. Dynamic responses of the 2D model for saturated porous media.
Figure 6. Dynamic responses of the 2D model for saturated porous media.
Applsci 12 04586 g006
Figure 7. Comparison of computation time.
Figure 7. Comparison of computation time.
Applsci 12 04586 g007
Figure 8. Computational model of 2D liquefiable site.
Figure 8. Computational model of 2D liquefiable site.
Applsci 12 04586 g008
Figure 9. Acceleration time history curve of the Wolong station seismic record of Wenchuan earthquake.
Figure 9. Acceleration time history curve of the Wolong station seismic record of Wenchuan earthquake.
Applsci 12 04586 g009
Figure 10. Time history curves of horizontal and vertical displacements at measuring points D1–D4.
Figure 10. Time history curves of horizontal and vertical displacements at measuring points D1–D4.
Applsci 12 04586 g010
Figure 11. Pore pressure time history curves of the model at measuring points W1–W4.
Figure 11. Pore pressure time history curves of the model at measuring points W1–W4.
Applsci 12 04586 g011
Figure 12. Contour plots of horizontal and vertical displacements of the model.
Figure 12. Contour plots of horizontal and vertical displacements of the model.
Applsci 12 04586 g012
Figure 13. Contour plots of pore pressure ratio of the model.
Figure 13. Contour plots of pore pressure ratio of the model.
Applsci 12 04586 g013
Table 1. Parameters of linear elastic saturated soil.
Table 1. Parameters of linear elastic saturated soil.
Basic Material Parameters
ParameterValue
ES3000 PaYoung’s modulus
ρ0.306 kg/m3Density of two-phase media
ρf0.2977 kg/m3Density of fluid-phase
n0.333Porosity
ν0.2Poisson’s ratio
k ¯ 0.004883 m3s/kgDynamic permeability coefficient
λ833.3 PaLame’s constants of the soil skeleton
G1250 PaShear modulus
Variable material parameters considered in the example
Parameter ParametersNo. 1No. 2
Kf0.3999 × 105 Pa0.6106 × 105 PaBulk modulus of pore fluid
Ks0.5005 × 104 PaBulk modulus of soil skeleton
Qb0.1201 × 106 Pa0.1385 × 105 PaCompressibility coefficient of pore fluid
Wave propagation velocity in saturated soil
ParameterActual wave velocity
No. 1No. 2
Cp635.12 m/s176.15 m/sP wave velocity
Cs63.92 m/s89.69 m/sS wave velocity
Table 2. Soil element parameters.
Table 2. Soil element parameters.
Loose Sand LayerDense Sand Layer
Element thickness1.0 m1.0 m
Vertical gravitational acceleration−9.81 m/s2−9.81 m/s2
Liquid-phase undrained bulk modulus4.4 × 106 kPa7.3 × 106 kPa
Horizontal permeability coefficient1 × 10−4 m/s1 × 10−5 m/s
Vertical permeability coefficient1 × 10−4 m/s1 × 10−5 m/s
Table 3. Material parameters of the model.
Table 3. Material parameters of the model.
ParametersLoose Sand LayerDense Sand Layer
Density   ρ 1.7 ton/m31.9 ton/m3
Reference   shear   modulus   G r 3.57 × 104 kPa2.59 × 104 kPa
Reference   bulk   modulus   B r 8 × 104 kPa6 × 104 kPa
Friction   angle   ϕ 33.531
Octahedral   peak   shear   strain   γ max 0.10.1
Pressure coefficient n0.50.5
Dilatancy   angle   ϕ P T 25.531
Shear   contraction   coefficient   c 1 0.0450.087
Shear   contraction   coefficient   c 3 0.150.18
Dilatancy   coefficient   d 1 0.060.0
Dilatancy   coefficient   d 3 0.150.0
Number of yield surface2020
Shear   contraction   coefficient   c 2 55
Dilatancy   coefficient   d 2 33
Liquefaction   coefficient   l 1 1 kPa1 kPa
Liquefaction   coefficient   l 2 0.00.0
Initial void ratio e0.60.6
Cs 10.90.9
Cs20.020.02
Cs 30.70.7
Standard pressure101101
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Song, J.; Xu, C.; Feng, C.; Wang, F. An Explicit Finite Element Method for Saturated Soil Dynamic Problems and Its Application to Seismic Liquefaction Analysis. Appl. Sci. 2022, 12, 4586. https://doi.org/10.3390/app12094586

AMA Style

Song J, Xu C, Feng C, Wang F. An Explicit Finite Element Method for Saturated Soil Dynamic Problems and Its Application to Seismic Liquefaction Analysis. Applied Sciences. 2022; 12(9):4586. https://doi.org/10.3390/app12094586

Chicago/Turabian Style

Song, Jia, Chengshun Xu, Chaoqun Feng, and Fujie Wang. 2022. "An Explicit Finite Element Method for Saturated Soil Dynamic Problems and Its Application to Seismic Liquefaction Analysis" Applied Sciences 12, no. 9: 4586. https://doi.org/10.3390/app12094586

APA Style

Song, J., Xu, C., Feng, C., & Wang, F. (2022). An Explicit Finite Element Method for Saturated Soil Dynamic Problems and Its Application to Seismic Liquefaction Analysis. Applied Sciences, 12(9), 4586. https://doi.org/10.3390/app12094586

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