A completely explicit finite element method for solving dynamic u-p equations of fluid-saturated porous media

Soil Dynamics and Earthquake Engineering 97 (2017) 364–376

Contents lists available at ScienceDirect

Soil Dynamics and Earthquake Engineering journal homepage: www.elsevier.com/locate/soildyn

A completely explicit finite element method for solving dynamic u-p equations of fluid-saturated porous media

MARK

⁎

Chengshun Xu, Jia Song, Xiuli Du , Zilan Zhong The key laboratory of Urban Security and Disaster Engineering of Ministry of Education, Beijing University of Technology, Beijing 100124, china

A R T I C L E I N F O

A BS T RAC T

Keywords: Fluid-saturated porous media u-p formed equations Diagonalization method Completely explicit finite element method

A completely explicit finite element method is developed to efficiently solve the dynamic u-p equations of fluidsaturated porous media (FSPM). In the spatial domain, a diagonalization method is applied to both the mass and the fluid compressibility matrices by neglecting the interactions between adjacent nodes in the finite element model of FSPM for the first time. In the time domain, an explicit-explicit algorithm consisting of the central and the unilateral difference methods is applied to solve the dynamic u-p equations. The stability and accuracy of the proposed explicit-explicit algorithm are discussed. Moreover, five numerical examples were presented in this paper. The good agreements between the results obtained from the proposed method and the existing explicit-implicit method indicate that the diagonalization method for the fluid compressibility matrix and the proposed completely explicit method are validity in solving the dynamic u-p equations of FSPM.

1. Introduction The fully saturated soil is generally idealized as fluid-saturated porous media (two-phase media) with deformable porous solid skeleton and fluid saturated in the solid pore structure. Different from single-phase media, the fluid-saturated porous media have much more complicated macro mechanical properties due to the interaction between the solid skeleton and the pore fluid (the fluid-solid coupling), which is especially important in the dynamic analyses. Therefore, based on the quasi static consolidation model [1], Biot [2] originally proposed the theory of wave propagation in saturated porous media, which took account for the inertia of soil skeleton and pore fluid. Zienkiewicz et al. [3,4] developed various approximations of Biot's theory for engineering practice, and also indicated the corresponding limits of those approximations in different engineering problems such as the consolidation problem and the undrained behavior. Several representative formulations are presented by Zienkiewicz et al. [3,4] such as: u-U-p，u-U， u-p，u-w forms, where, u is the displacement of the soil skeleton, U is the absolute displacement of the pore fluid, w is the pore fluid displacement with respect to the soil skeleton, and p is the pore pressure, respectively. Among these different formulations, the u-p formed equations have relatively simpler formulation and higher computational efficiency because of the following advantages. First, the terms in the u-p formed equations involving the pore fluid acceleration with respect to the soil skeleton,ẅ , can be neglected if the motion of fluid, ẇ , is relatively small. Based on this reasonable ⁎

assumption, the u-p formed equations can be rewritten in a simpler formulation. Second, the u-p formed equations, consisting of an unknown vector of displacement, u, and an unknown scalar of pore pressure, p, have fewer degrees of freedom compared to the fullyvectored formulations (i.e. u-U form). Finally, as one of the most important parameters in soil dynamics, the pore pressure p can be obtained directly from the u-p formulation without iterations. Therefore, the u-p formed equations are more appropriate and convenient for describing the dynamic problems of the saturated porous media. For the dynamic equations of the saturated porous media, the completely analytical solutions can only be obtained under a few specific boundary and initial conditions. Therefore, the numerical methods such as the finite element method (FEM) associated with appropriate time integration schemes become the primary tools to solve dynamic equations in most situations. The entire analysis process can be summarized by two steps: first, discretization of the infinite continuous media in spatial domain; second, solving the dynamic equations of the discrete system in discrete time domain. Moreover, the integration of the differential equations of system in time domain can be divided into two types: the implicit and the explicit methods. Compared to implicit methods, explicit methods are more computationally efficient with decoupled equations, especially in solving the complicated nonlinear dynamic problems with large amount of degrees of freedom. The central difference method (CDM) is one of the most commonly used explicit integration methods. However, for dynamic

Corresponding author. E-mail address: [email protected] (X. Du).

http://dx.doi.org/10.1016/j.soildyn.2017.03.016 Received 22 November 2015; Received in revised form 4 February 2017; Accepted 13 March 2017 0267-7261/ © 2017 Published by Elsevier Ltd.

Soil Dynamics and Earthquake Engineering 97 (2017) 364–376

C. Xu et al.

method is implemented into both of the mass and fluid compressibility matrices to create an important precondition of an explicit-explicit temporal algorithm. The assumption and detail procedure of the diagonalization method are also illustrated in Section 3. In Section 4, an explicit-explicit algorithm consisted of the CDM and UDM, is proposed to solve the u-p formed equations in time domain. In Section 5, five numerical examples are given to demonstrate the effectiveness of the diagonalization for the fluid compressibility matrix and the accuracy of the proposed completely explicit FEM. In appendix, the stability conditions and theoretical accuracy of the complete explicit FEM are discussed in detail. And conclusions drawn from this study were presented at the end of this paper.

problems with non-diagonalized damping matrix, the unilateral difference method (UDM) is usually used to solve the damping term, which limits the accuracy of combined method to be the first order. To improve the accuracy of calculation, Li et al. [5], Du and Wang [6] established two explicit difference methods based on the CMD for the dynamic equations of damped systems, respectively; Wang and Du [7] developed another explicit difference method which approximate the displacement by the quadratic Lagrangian interpolation functions. In consequence, the use of the explicit method for the saturated porous media, consisting more governing dynamic equations and unknown parameters, leads to more high-efficient computation. In last three decades, many researchers proposed a variety of timeintegration methods incorporated in the FEM for solving the dynamic coupled problems of fluid-saturated porous media. The dynamic equations in u-U and u-w form are usually solved by the explicitexplicit [8,9], explicit-implicit [10,11] or implicit-implicit [11] algorithms combined with a staggered procedure [12]. For the u-p form, Park [13] and Zienkiewicz et al. [14] established two different types of unconditionally stable implicit-implicit staggered algorithms respectively. Huang and Zienkiewicz [15] improved Zienkiewicz's method [14] by introducing a pressure correction method. However, because up formed equations are a mixed form consisted of vector, u, and scalar, p, in the case of nearly incompressible pore fluid and small permeability, the numerical results usually have significant fluctuations and even lead to non-convergence if interpolation functions with the same order are selected for u and p during the discretization in spatial domain. Therefore, the Babuska-Brezzi stability criteria [16,17] should be satisfied for interpolation functions of u and p to limit the fluctuations in the results, but it complicates the coding in the FEM [18]. Zienkiewicz et al. [18] introduced a new explicit-implicit algorithm based on operator splitting before spatial discretization which effectively overcome the limitation and improved the stability of the algorithm. Huang et al. [19] further developed an unconditionally stable implicit-implicit algorithm by introducing a splitting factor δ to adjust the time step to guarantee a stable solution. Park and Tak [20] established a stable algorithm consisting of the multiple time steps, the remeshing step and the sub-iteration step to remove the instability in the solutions during the finite element analyses. Soares et al. [21] divided the global two-phase model into solid and fluid sub-domains as uncoupled models and performed finite element analyses on each subdomain independently. The equation for solid sub-domain was solved by implicit Newmark method, while the equation for fluid sub-domain depicted by a new scalar pore pressure wave formulation was solved by the explicit Green-Newmark method. The interactions between the solid and fluid sub-domains were achieved by applying a couple of surface forces in the two sub-domains. Markert et al. [22] compared the implicit monolithic scheme and the semi-explicit-implicit splitting scheme, and provided reliable recommendations for the applications of these two schemes for different dynamic problems of saturated porous media. It should be observed that those time integration methods mentioned above for solving the u-p formed equations are always explicit-implicit or implicit-implicit algorithms, and explicit-explicit algorithm has never proposed by any researchers. Because the fluid dynamic equation is in the form of a diffusion equation, the fluid compressibility matrix in the fluid dynamic equation is never formed into a lumped matrix by a diagonalization method, which always is used in the mass matrix. This limits the choice to the implicit method for the fluid dynamic equation, leading to larger computational costs and lower efficiency. In this paper, a completely explicit FEM, that involves the standard FEM with diagonalization of mass and fluid compressibility matrices in spatial domain and an explicit-explicit algorithm in time domain, is proposed to solve the dynamic fluid-solid coupled problems for the saturated porous media in u-p formulation. In Section 2, a brief summary of the u-p formed equations and the standard Galerkin FEM are introduced in spatial domain. In Section 3, a diagonalization

2. u-p formed Equations and FEM The u-p formed equations for the fluid-saturated porous media are derived in by Zienkiewicz et al. [3,4]. The general governing equations are shown below: ⎧ LT σ + ρb − ρü = 0 ⎪ ⎪ ⎪− k ∇p + k ρf b = ẇ ⎨ σ = σ ′ − αmp = Dε − αmp ⎪ 1 ⎪ T T ⎪− ∇ w = Q p + αm ε ⎩ b

dynamic equilibrium equation of solid −fluid mixture a) dynamic equilibrium equation of fluid b) constitutive relationship

c)

mass balance equation of the flow

d)

(1)

The definitions of the symbols used in Eq. (1) are listed in Table 1. Substituting Eqs. (1c) and (1d) into Eqs. (1a) and (1b), and neglecting the body force acceleration, b, the u-p formed equations can be rewritten as: T ⎧ a) ⎪ L (σ ′ − αmp ) − ρü = 0 ⎨ T 1 ⎪∇ u̇ − ∇T k ∇p + p ̇ = 0 b ) Qb ⎩

(2)

The boundary and the initial conditions are given by Eqs. (3) and (4), respectively. Table 1 Definitions of the symbols. Symbols

Definitions

u w p ρ ρf

Displacement of soil skeleton Displacement of pore fluid with respect to the soil skeleton Pore pressure Density of the solid-fluid mixture Density of the fluid

n b k k ε

Porosity Body force acceleration Permeability coefficient of pore fluid Dynamic permeability coefficient Strain of soil skeleton

σ σ′ α Qb

k = k /ρf g

Geometric equations of solid phase: ε = Lu Total stress of the two-phase media Constitutive relationship: Effective stress of the soil skeleton σ = σ ′ − α mp ; Compressibility α = 1 − KD /KS , 1/Qb = n /Kf + (α − n )/KS , coefficient of solid KD = E /[3(1 − 2ν )] = λ + 2G /3; where, Compressibility KD 、KS 、Kf are the bulk modulus of soil coefficient of skeleton, soil particles and pore fluid fluid respectively.For the material of soil, soil

Matrices D, L, ∇ and m in twodimension:

particles are nearly incompressible relative to the soil skeleton, thus they have the relations of KS ≫ KD , α ≈ 1. ⎡∂ ∂ ⎤ ⎡ λ + 2G λ 0⎤ ⎢ ∂x 0 ∂y ⎥ D=⎢ λ λ + 2G 0 ⎥ ; LT = ⎢ ⎥; ∂ ∂ ⎢ ⎥ ⎢0 ⎥ ⎣ 0 0 G⎦ ∂y ∂x ⎦ ⎣

⎡∂ ∇T = ⎢ ∂x ⎣

∂ ⎤ ; ∂y ⎥ ⎦

mT = [1 1 0 ]; where, λ and G are the Lame constants of soil skeleton.

365

Soil Dynamics and Earthquake Engineering 97 (2017) 364–376

C. Xu et al.

on Γu

Table 2 Mechanical parameters utilized in the study.

on Γt on Γp on Γq

(3)

u = u0 ⎧ ⎪ ⎨ u̇ = u̇0 ⎪ ⎩ p = p0

(4)

The Galerkin method [23] is applied to discretize the numerical model of fluid-saturated porous media in spatial domain. The spatial approximation of displacements of soil skeleton, pore pressure and external loading can be achieved by the following shape functions:

u = Nu u , p = Np p , fu = Nu fue , fq = Np fqe

(5)

(6)

where,

⎧ M = ∫ NT ρN dΩ u ⎪ Ω u ⎪K = T DLN dΩ ( ) LN ∫ u u ⎪ d Ω ⎪Q = T ∇ dΩ N N ∫ p ⎪ Ω u ⎪ 1 ⎨ S = ∫Ω NpT Q Np dΩ b ⎪ ⎪ J = ∫ (∇Np)T k (∇Np) dΩ Ω ⎪ ⎪ fu = ∫ NuT ⋅∼ fu ⋅dΓ Γ ⎪ ∼ T ⎪ fq = ∫ Np ⋅fq ⋅dΓ ⎩

(7)

Γ

As listed above, M is the mass matrix of two-phase media, Kd is the stiffness matrix of soil skeleton, Q is the coupling matrix, S is the compressibility matrix of the pore fluid, J is the permeability matrix of the pore fluid, and fu , fq are respectively the external mechanical loadings.

(a)

Value

Young’s modulus Bulk modulus of pore fluid Bulk modulus of solid Density of solid and pore fluid Density of pore fluid porosity Poisson’s ratio Hydraulic conductivity

ES Kf Ks ρs ρf n ν k

70 MPa 2.24×103 Mpa 1×1017 Mpa 1700 kg/m3 1000 kg/m3 0.4 0.3 1×10−3 m/s

In Eq. (7), the mass matrix M and the fluid compressibility matrix S assembled by the conventional finite element discretization method are consistent matrices. The inversions of M and S are always full matrices indicating that there are interactions between different nodes. However, those simultaneous equations are generally solved by the implicit method, which takes up lots of memory of the computers, costs more computing times, and has low efficiency especially in solving nonlinear problems with numerous degrees of freedom. Therefore, in order to implement the explicit method, a diagonalization method is used to simplify the matrices M and S by assuming that the total mass and the compressive deformation of the solid-fluid system are concentrated at each node to eliminate the coupling between adjacent nodes. The following assumptions are made for the mass matrix, M: (1) the mass of the whole system is conservative; (2) the differences of the inertia forces among the nodes in the same element can be neglected after the nodal vibration reaching a steady state. Based on these assumptions, the mass of the solid-fluid system is concentrated at the nodes, which indicates the inertia terms of adjacent nodes in Eq. (6a) are independent to each other. Similarly, it is also assumed that the total compressive deformation of the whole system is constant and the differences of the compressive deformations of the nodes in a same element can be neglected after the nodal vibration reaching a steady state. Therefore, the compressive deformations of the system are concentrated at the nodes, and compressive deformations of adjacent nodes are uncoupled. The detailed procedure of the diagonalization methods [24,25] is presented as follows.

a) b)

Symbol

3. Diagonalization of the mass and fluid compressibility matrices

After discretization of the numerical model, the Galerkin weak form of the u-p formed equations for fluid-saturated porous media can be obtained as follows [3,4]:

⎧ ⎪ Mu ̈ + K d u − Qp = f u ⎨ T ⎪ ⎩ Sp ̇ + Jp + Q u̇ = fq

Parameters

(b)

(c)

2

Load [kPa]

∼ ⎧u = u ⎪ ∼ ⎪ t = σ ′n − pn = t ⎨ ∼ ⎪p = p ⎪ (−k ⋅∇p )T n = q∼ ⎩

Step Loading

1

0

0

15

30

Time [s]

Load [kPa]

2

Sine loading

1

0

0

15

30

Time [s] Fig. 1. (a) Two-dimensional computational rectangle model; (b) Two-dimensional computational trapezoid model; and (c) the two types of applied loadings involving the step loading and the sine loading.

366

Soil Dynamics and Earthquake Engineering 97 (2017) 364–376

C. Xu et al.

Fig. 2. Comparison of numerical results using three integration methods. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article).

For the lumped mass matrix of each element after diagonalization, Mle , the entries on its diagonal are proportional to those on the diagonal of the consistent mass matrix of each element, Me , while the rest of the entries in Mle are set to be zeros, as shown in Eq. (8). Moreover, the arithmetic sum of all the entries of Mle remains the same with that of Me as shown by Eq. (9). Same rules are applied to diagonalize the fluid compressibility matrices Se as shown in Eqs. (10) and (11):

⎧ a (Me ) = a ∫ ρNT N dV ( j = i ) ii i i Ve (Mle )ij = ⎨ ( j ≠ i) ⎩ 0

Table 3 Mechanical parameters of saturated porous media. Linear elastic mechanical parameters of the saturated porous media

⎪ ⎪

(8)

p(surface, t) = 0, σ(surface, t) = σ0

Free drainage

L=

=100

Mechanical parameters

Value

Description

ES Kf Ks ρ

3000 Pa 0.6106×105 Pa 0.5005×104 Pa 0.3060 kg/m3

ρf n ν k λ G Qb

0.2977 kg/m3 0.333 0.2 0.004883 m/s 833.3 Pa 1250 Pa 0.1385×105 Pa

Young’s modulus Bulk modulus of pore fluid Bulk modulus of solid Density of solid and pore fluid Density of pore fluid Porosity Poisson’s ratio Permeability Lame constant Shear modulus Fluid compressibility

Wave velocities of the saturated soil Parameters Non-dimensional velocity Cp 1.0 Cs 0.5092

No flow, no expansion laterally constrained

ne

Actual velocity

Description

176.2 m/s 89.69 m/s

Velocity of P wave Velocity of S wave

ne

∑ (Mle)ii = a ∑ (Me)ii = WId = ρVe Id i =1

i =1

⎧ b ( Se ) = b ∫V e Q1b NiT Ni dV ( j = i ) ii (Sle )ij = ⎨ ( j ≠ i) ⎩ 0

No flow, fixed

(9)

⎪

⎪

ne

ux (L, t) = 0, uy (L, t) = 0

ne

∑ (Sle)ii = b ∑ (Se)ii = VId =

Fig. 3. One-dimensional numerical model of saturated porous media.

i =1

367

i =1

1 Ve Id Qb

(10)

(11)

Soil Dynamics and Earthquake Engineering 97 (2017) 364–376

C. Xu et al.

Fig. 4. Displacement and pore pressure time histories of one-dimensional model at different moments.

where, Mle and Me are the lumped and consistent mass matrices in an element, respectively; Sle and Se are the lumped and consistent fluid compressibility matrices in an element, respectively; a, b are corresponding scaling coefficients; W and V are the mass and compressive deformation in an element, respectively; Ve is the area of an element in two-dimension; and Id is unit matrix.

It can be found that the displacements, u, and pore pressure, p, in the equations above are decoupled, and the solutions of each time step can be obtained based on the solutions of the last two steps. Therefore, the proposed completely explicit FEM, consisting of the diagonalization method in spatial domain and the explicit-explicit algorithm in time domain, has significant advantages compared to the traditional methods in performing dynamic analyses on the fluid saturated porous media by saving computers’ memory, reducing computing time and improving computational efficiency especially.

4. Explicit-explicit algorithm in time domain An explicit-explicit algorithm is proposed to solve the u-p formed equations in time domain as shown in Eq. (6). In the algorithm, the CDM is applied to solve Eq. (6a) and the UDM to solve Eq. (6b). The detailed procedures are presented in this section. In the dynamic equilibrium equation of the solid-fluid mixture, Eq. (6a), the acceleration uk̈ of soil skeleton at the kth time step can be expressed by Eq. (12) using the CDM:

uk̈ = (uk +1 − 2uk + uk −1)/ Δt 2

5. Numerical examples Five numerical models with linear elastic porous materials are used in this section as examples to illustrate the effectiveness of the completely explicit FEM. The numerical results obtained from the proposed method are validated by comparing to those from Zienkiewicz's explicit-implicit method and the analytical solutions developed by Simon et al. [26].

(12)

where, Δt is the length of the time step. By substituting Eq. (12) into Eq. (6a), the recursion relation for the two-step method of displacement uk+1 at the (k+1)th time step can be written as follows:

uk +1 = (2I − Ml −1Kd Δt 2 ) uk − uk −1 + Ml −1QΔt 2pk + Ml −1Δt 2fuk

5.1. Example 1 Two planar models of saturated soil, rectangular- and trapezoidalshaped models are illustrated in Fig. 1. For the model, the top surface is a free drainage boundary; the bottom surface is an impermeable boundary and is fixed in both horizontal and vertical directions; and the two lateral boundaries are impermeable boundaries, where the nodes are fixed in the horizontal direction but are allowed to move freely in the vertical direction. The trapezoidal-shaped model has the same boundary conditions as the rectangle model except for the right lateral boundary, where the nodes can move freely in both horizontal and vertical directions. Same soil properties adopted from literature [27,28] are listed in Table 2. Two different surface loadings are considered in this example: a step loading of 1.0 kPa and a sinusoidal loading, σ = 1.0 + 0.5 sin(2.0t ) kPa, as shown in Fig. 1c, are imposed at the top surface of two models with a duration of 30 s. For both loading conditions, the length of time step is selected as Δt=1×10−4 s during the numerical analyses. Three nodes A, B, and C of the two models with corresponding coordinates as high-lightened in Figs. 1a and b are selected as the target nodes. For the case of sinusoidal loading, the displacement and pore pressure time histories of the three nodes in rectangular-shaped model are presented in Figs. 2a and b. The red and blue curves in Figs. 2a and b represent the results obtained from the proposed method using diagonalized and non-diagonalized fluid compressibility matrix, respectively. It can be observed that the red curves of displacement and pore pressure are in good agreement with the blue curves, which validate the application of the diagonalization method to the fluid

(13)

where, I is a unit matrix with the same order as K matrix. Ml is lumped mass matrix. For the dynamic equilibrium equation of fluid, Eq. (6b), the UDM is applied to both the pore pressure term and the velocity term. Specifically, the first-order derivative of the pore pressure pk̇ and the velocity uk̇ at the kth time step can be rewritten by the forward and the backward difference method, respectively, as shown in Eqs. (14) and (15):

pk̇ = [ pk +1 − pk ]/ Δt

(14)

uk̇ = [uk − uk−1]/ Δt

(15)

By substituting Eqs. (14) and (15) into Eq. (6b), the pore pressure pk+1 at the (k+1)th time step can be expressed as:

pk +1 = (I − Sl −1JΔt ) pk + Sl −1QT (−uk + uk −1) + Sl −1Δtfqk

(16)

where Sl is lumped fluid compressibility matrix. Eqs. (13) and (16) can be written together in a monolithic formulation as shown in Eq. (17):

⎧ uk +1 ⎫ ⎡ 2I − Ml −1Kd Δt 2 Ml −1QΔt 2 ⎤ ⎧ uk ⎫ ⎡ − I 0 ⎤ ⎧ uk −1 ⎫ ⎨ ⎬=⎢ ⎬ ⎥⎨ ⎬ + ⎢ ⎥⎨ ⎩ pk +1 ⎭ ⎣ − Sl −1QT I − Sl −1JΔt ⎦ ⎩ pk ⎭ ⎣ Sl −1QT 0 ⎦ ⎩ pk −1 ⎭ ⎧f ⎪ ⎫ ⎡ M −1Δt 2 0 ⎤⎪ uk ⎨ ⎪ ⎬ ⎥⎪ +⎢ l −1 Sl Δt ⎦ ⎩ fqk ⎭ ⎣ 0

(17) 368

Soil Dynamics and Earthquake Engineering 97 (2017) 364–376

C. Xu et al.

Impulsiive Loading P (0, 0) 0

X

50 m

A (0, 2) 2

Z B (0, 25)

Viscouss-Spring Bounndaries 50 m Fig. 5. Numerical model under local impulsive loading with high frequency.

time histories of the three nodes in trapezoidal-shape model, are presented in Figs. 2c and d. In the two figures, the good agreements between the results obtained from both the proposed method and Zienkiewicz's explicit-implicit method indicate that the proposed diagonalization method is applicable for the model with irregular (in the trapezoidal-shaped model) shapes of elements. It should be pointed out that severe fluctuation can be observed at the beginning stage of the displacement and pore pressure time histories obtained from both the proposed method and Zienkiewicz's explicit-implicit method for the trapezoidal-shape model. It is primarily caused by the oscillations of the unrestrained right boundary in its horizontal direction.

compressibility matrix. Furthermore, the curves obtained from the proposed method using the diagonalized fluid compressibility matrix are smoother. Moreover, the displacement and pore pressure time histories obtained from Zienkiewicz's explicit-implicit algorithm [29] (green curves) and the proposed method (red curves) are compared in Figs. 2a and b, respectively. The consistency of both the displacement and pore pressure time histories obtained from the two methods validates the proposed method in solving dynamic u-p equations of fluid-saturated porous media. For the case of step loading, the displacement and pore pressure

Fig. 6. Numerical results of the model under high frequency impulsive loading. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article).

369

Soil Dynamics and Earthquake Engineering 97 (2017) 364–376

C. Xu et al.

η H L

z

Wavee

d

Dynamiic Pressure Pb

x

0

A (100, -10)

Sand d

200 m

Viscous-Spring V g Boundary B (100, -100)

200 m Fig. 7. The seabed model involving the wave loading.

Table 3. The numerical solutions of displacement and pore pressure calculated from the proposed method at three dimensionless moments, τ=20, 40 and 60, are compared with Simon's analytical solutions in Fig. 4. These representative moments are selected so that reflected waves do not present in the region of interest in the finite element model. The good agreements between the numerical solutions and the analytical solutions again verify the accuracy of the proposed method. However, in the area closed to the top surface, ξ = 0~20 , the numerical solution of the displacements and the pore pressures have a certain gap with the analytical solution indicate a relative low accuracy for the proposed method.

Table 4 Wave and soil characteristics used in example 3. Wave characteristics

Value

Soil characteristics

Value

Wave period T Wave height H

8.0 s 2.0 m

Water depth d

20 m

Wave length L

88.88 m

Shear modulus G Poisson’s ratio ν Permeability k Density of soil ρs Density of water ρf Porosity n Bulk modulus of pore fluid Kf

79.6 MPa 0.4 1×10−4 m/s 2650 kg/m3 1000 kg/m3 0.4 2.24×103 Mpa

5.2. Example 2 Fig. 3 illustrates a 1-dimensional (1-D) model of saturated soil developed by Simon et al. [26]. Same mechanical properties as specified by Simon et al. [26] in his model are listed in Table 3. The saturated soil is assumed to surround by impermeable, frictionless and rigid boundaries except for its top surface, which is assumed to be a perfect drainage boundary. A constant distribution loading with a dimensionless amplitude of 1.0 is applied at the top surface. The numerical solution obtained from the proposed completely explicit FEM (solid curves) are compared to the analytical solutions derived by Simon et al. [26] (dash curves). In this example, the dimensionless time t index, τ , is calculated by τ = ρk , where t is the time; ρ is the density of the solid-fluid mixture; and k is the permeability coefficient of pore x fluid. And dimensionless depth index, ξ , is calculated by ξ = ρkV , where c x is the depth of the soil; Vc is the velocity of P1 wave (Cp), as shown in

5.3. Example 3 The simplified assumption makes the u-p formulation not suitable for solving the dynamic problems under high frequency loadings. But in this part, a numerical example is employed only to illustrate the ability of explicit algorithm to solve the problems under impulsive loading with high frequency component, rather than to contradict with the limitation of the u-p formulation. The saturated soil model in two-dimensional (2-D) half infinite space is shown in Fig. 5. The lateral and bottom boundaries are set as viscousspring boundaries to absorb the spurious waves reflected from the truncated boundary. The material properties are same with those in Example 1. Moreover, a point impulsive loading with high frequency components is imposed on the top surface of the model as shown in Fig. 5.

Fig. 8. Numerical results of the seabed under wave loading. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article).

370

Soil Dynamics and Earthquake Engineering 97 (2017) 364–376

C. Xu et al.

Fig. 9. The nonlinear layer slope model.

Fig. 6 shows the displacement and pore pressure curves of the A and B points in the 2-D saturated soil model. The red and blue curves represent the results of the proposed method and Zienkiewicz's method, respectively. From Fig. 6, we can see that the proposed method has stable results, and the solutions from both methods are in a good agreement with each other. This shows the ability of the explicit algorithm to solve the problems under high frequency loadings. 5.4. Example 4 In order to analyze the dynamic behavior of saturated seabed under wave loadings, the proposed method and Zienkiewicz's method [29] are employed to compute the model as shown in Fig. 7. The bottom of the seabed is fixed and impermeable boundary, and both of the lateral boundaries are viscous-spring boundaries. The free draining surface of the seabed is imposed on the dynamic pressure based on the secondorder Stokes wave theory [30,31]. Referring to the literatures [31,32], mechanical parameters are listed in Table 4. Fig. 8 shows the displacement and pore pressure time curves of the two target points (A and B). The blue and red curves represent the results obtained from the Zienkiewicz's method and the proposed method, respectively. The results of the two methods have a good consistence. The consistency validates that the proposed explicit method is good way to deal with the wave-induced seabed dynamic problems.

Fig. 10. Time curves of displacement of slope nodes.

5.5. Example 5 It is well known that the explicit method has more significant advantages in the aspects of saving computing time and improving the computational efficiency than the implicit method for the dynamic nonlinear problems with numerous degrees of freedom. Thus, a nonlinear slope of seabed is applied to illustrate the computational efficiency of the proposed explicit method, as shown in Fig. 9. The slope is underlain by a rigid foundation, where the nodes are fixed against horizontal and vertical translation. The freedom columns

Fig. 11. Time curves of pore pressure of slope nodes and total computational time.

Table 5 Mainly material parameters.

Density/(ton/m3) Shear modulus Gr/kPa Bulk modulus Br/kPa Horizontal permeability/(m/s) Vertical permeability/(m/s)

Silt 1

Silt 2

Silt 3

Sand 1

Sand 2

Sand 3

Gravel 4

1.68 14046.9 42140.7 10−5 10−5

1.68 39020.0 117060.0 10−5 10−5

1.68 87793.4 263380.0 10−5 10−5

1.84 42735.4 128206.2 10−4 10−4

2.24 42735.4 128206.2 10−3 10−3

2.24 42735.4 128206.2 10−3 10−3

2.24 42735.4 128206.2 10−3 10−3

371

Soil Dynamics and Earthquake Engineering 97 (2017) 364–376

C. Xu et al.

the proposed method has the advantage of high efficiency.

in Fig. 9 are used to describe the free-field responses at the horizontal boundaries [33]. They are located far enough from the slope to ensure that the critical field is not affected by the horizontal boundaries. Most part of the seabed is submerged below the sea level. The nodes under the water, shown in Fig. 9, are imposed on the added mass. And corresponding nodes upper the sea level are fixed the third degree of freedom, namely the pore pressure is equal to zero, to keep a drainage path for the water. The layer seabed is consisted of cohesive and cohesionless soil, such as the sand, silt and gravel soils. The Pressure Depend Multi Yield 02 and Material object and the Pressure Independ Multi Yield and Material object [34,35] are used to describe the constitutive behaviors for the cohesionless and the cohesive soils, repectively. The mechanical parameters are listed in Table 5. After meshing, the model has 2781 degrees of freedom. The Kobe earthquake recorded at the Kobe University on January 17, 1995 is regarded as the uniform seismic excitation in horizontal direction. Horizontal component of the earthquake is depicted in Fig. 10. The proposed explicit method has been embed in OpenSees platform. The nonlinear seismic problem of the model is solved by the proposed explicit method and the Newmark implicit method, respectively. The displacement and pore pressure curves of key nodes of the slope are shown in the Figs. 10 and 11. In the two figures, the good agreements between the results obtained from both methods further validate the proposed method. The total computational time of both methods are listed in the Fig. 11. In OpenSees, both of the mass and fluid compressive matrices are diagonalized to avoid solving the simultaneous equations effectively. However, the computational mode in OpenSees is relatively regular, which require the stiffness matrix to be integrated in every computational step. In fact, due to this limitation it do not completely realize the explicit computing up to now. Even in the situation, it can be obtained from the Fig. 11 that the total time of the proposed method is the half of the Newmark method, although the time step of the proposed method is only one-tenth of the Newmark method. This illustrates that

6. Conclusions To efficiently solve the dynamic u-p equations of fluid-saturated porous media, the explicit-explicit algorithm composed of the central and unilateral difference methods is developed based on the diagonalization of both solid mass and fluid compressibility matrices of finite element. The stability and accuracy of the proposed method is discussed in addition. The proposed method is validated using five different numerical examples. Two main conclusions from this study are listed as follows: (1) The diagonalization for the fluid compressibility matrix is valid and likely leads to more stable results with relatively small fluctuations comparing to the solutions obtained from the proposed method with non-diagonal fluid compressibility matrix. (2) The completely explicit method is validated by comparing the numerical results with those from Zienkiewicz's explicit-implicit method and the Simon's analytical solutions. In order to illustrate the effect of the diagonalization for the fluid compressibility matrix, the simply explicit-explicit algorithm with only first-order accuracy is employed, which restricts severely accuracy of the complete method. Therefore, other explicit-explicit algorithms with higher accuracy will be studied in future. Acknowledgements This research work was financially supported by the Foundation for Innovative Research Groups of the National Natural Science Foundation of China (No. 51421005), and the National Natural Science Foundation of China (No. 51578026). The supports from the above organizations are gratefully acknowledged.

Appendix A A.1. Stability analysis This section discusses the stability of the completely explicit finite element method by examining the eigenvalues of the transfer matrix. Eqs. (13) and (16) can be combined together as follows:

⎧f ⎫ ⎧ vk +1 ⎫ ⎧v ⎫ ⎪ uk ⎪ ⎨ ⎬ = T ⎨ k ⎬ + L⎨ 0 ⎬ p p ⎩ k +1 ⎭ ⎩ k⎭ ⎪f ⎪ ⎩ qk ⎭

(A1)

where,

⎧ vk +1 = { uk +1 uk }T , vk = { uk uk −1}T a) ⎪ ⎡ 2I − Ml −1Kd Δt 2 −I ⎪ Ml −1QΔt 2 ⎤ ⎥ ⎪T = ⎢ I 0 0 ⎢ ⎥ b) ⎪ ⎢⎣ Sl−1 QT 1 − Sl −1JΔt ⎥⎦ −Sl −1QT ⎨ ⎪ ⎡ Ml −1Δt 2 0 0 ⎤ ⎪ ⎥ ⎪L = ⎢ 0 c) 0 0 ⎥ ⎢ ⎪ ⎢⎣ 0 0 Sl −1Δt ⎥⎦ ⎩

(A2)

In Eq. (A2b), T is the transfer matrix. The stability conditions of a numerical method in time domain depend on the spectral radius of the transfer matrix, denoted as ρ (T ) = ( λ1 , λ2 , λ3 , ⋯)max , where λi are the eigenvalues of the transfer matrix T. The numerical method is stable only if ρ (T ) < 1 [36,37]. In order to make the derivation straightforward and compact, the matrices Ml, Kd, Q, J, and Sl are regarded as scalars m, k, q, j and s. The eigenvalues of the scalar formed transfer matrix can be obtained as follows:

372

Soil Dynamics and Earthquake Engineering 97 (2017) 364–376

C. Xu et al.

λ−2+ λE − T =

k Δt 2 m

−1

1

−

q Δt 2 m

λ 0 j q − s λ − 1 + s Δt

q s

=0 (A3)

where, E is a three-order unit matrix. By setting

λ3

+

Aλ2

Ω2

=

k Δt 2 ，ϕ1 m

=

j Δt ，αΩ 2 s

=

q2 Δt 2 , ms

where α =

q2 , ks

Eq. (A3) can be rewritten as: (A4)

+ Bλ + C = 0

where,

⎧ A = Ω2 + ϕ − 3 1 ⎪ ⎨ B = ϕ1 Ω 2 − Ω 2 − 2ϕ1 + αΩ 2 + 3 ⎪ ⎩C = ϕ1 − αΩ 2 − 1

λ=

(A5)

Moreover, the stability requirement of a numerical method, λ ≤ 1, is transformed into an equivalent negative plane domain (z ≤ 0 ) by setting 1+z and we arrive at: 1−z

a 0 z 3 + a1 z 2 + a2 z + a3 = 0

(A6)

where,

⎧ a 0 = 8 − 4ϕ − 2Ω 2 + 2αΩ 2 + ϕ Ω 2 1 1 ⎪ ⎪ a1 = 4ϕ1 − 4αΩ 2 − ϕ1 Ω 2 ⎨ ⎪ a2 = 2Ω 2 + 2αΩ 2 − ϕ1 Ω 2 ⎪ a = ϕ Ω2 ⎩ 3 1

(A7)

Based on the Routh-Hurwitz [38] criterion, Eq. (A6) has negative real solutions if and only if the coefficients ai in Eq. (A7) meet the following requirements:

⎧ a0 > 0 ⎪ ⎨ ai ≥ 0(i = 1, 2, 3) ⎪Δ = a a − a a ≥ 0 ⎩ 2 1 2 0 3

(A8)

By substituting Eq. (A7) into Eq. (A8), the necessary and sufficient conditions for the stability of the proposed method are obtained as follows:

⎧8 − 4ϕ − 2Ω 2 + 2αΩ 2 + ϕ Ω 2 > 0 1 1 ⎪ ⎪ ϕ1 Ω 2 ≥ 0 ⎪ ⎨ 2Ω 2 + 2αΩ 2 − ϕ1 Ω 2 ≥ 0 ⎪ 2 2 ⎪ 4ϕ1 − 4αΩ − ϕ1 Ω ≥ 0 ⎪ ϕ − Ω 2 − αΩ 2 ≥ 0 ⎩ 1

(A9)

Moreover, because the CDM and the UDM are adopted to solve Eqs. (13) and (16), respectively, the following stability requirements of the CDM and the UDM should also be satisfied by the proposed method.

⎧ Ω2 ≤ 4 ⎨ ⎩ ϕ1 ≤ 2

(A10)

Eq. (A9) can be furthermore rewritten in a simplified form as shown below:

⎧ (4 − Ω 2 )(2 − ϕ ) + 2αΩ 2 > 0 1 ⎪ ⎪ ϕ1 Ω 2 ≥ 0 ⎪ ⎨ 2(1 + α ) ≥ ϕ1 ⎪ 2 2 ⎪ 4ϕ1 − 4αΩ ≥ ϕ1 Ω ⎪ 4ϕ − 4αΩ 2 ≥ 4Ω 2 ⎩ 1

a) b) c) d) e)

(A11)

Given the requirements shown in Eq. (A10), the first three stability restrictions, Eq. (A11a, b and c) are automatically satisfied. Moreover, it is found that Eq. (A11e) has stricter restriction than Eq. (A11d) in the given range of ϕ1 ≤ 2 as shown in Eq. (A10). Therefore, the stability requirements of the proposed algorithm can be further simplified as follows:

⎧Ω ≤ 2 ⎪ ⎨ ϕ1 ≤ 2 ⎪ (1 + α ) Ω 2 ≤ ϕ ⎩ 1

a) b) c)

(A12)

Based on the stability restrictions (A12), the stable zone of the proposed algorithm is also illustrated in Fig. A1. The space enclosed by the red surface, the blue plane and the vertical plane of Ω =0 in Fig. A1 indicates the stable zone of the proposed method. It can be observed from Fig. A1 that the third restriction (i.e. Eq. (A12c)) is the critical one as shown by the red surface. k

j

Finally, substitutions of Ω 2 = m Δt 2 ，ϕ1 = s Δt ，αΩ 2 = needs to satisfy the following restriction:

q2 Δt 2 ms

into stability restriction (A12c), the time step length, Δt , about the stable condition

373

Soil Dynamics and Earthquake Engineering 97 (2017) 364–376

C. Xu et al.

Fig. A1. The stable area of the proposed method.

Δt ≤

mj ks + q 2

(A13)

Two aspects should be noticed from the Eq. (A13). One is that time step is proportional to the value of permeability, and the other is the equation cannot reflect how the time step varies with element size. Based on the two aspects, the affections of the permeability and the element size to the stable condition of the proposed method are illustrated. The real minimum time step of three simple models (one element in size of 1 m×1 m, 5 m×5 m, and the model of 25 elements in the element size of 1 m×1 m) and the theory time steps (the Equation (A13))are compared in Table A1, respectively. The sine loading imposed on the top of the three models and the material parameters are shown in Table A1. The scalars of matrices Ml, Kd, Q, J and Sl are chosen approximately according to the material parameters.

k=10−3 m/s

k=10−4 m/s

k=10−5 m/s

k=10−6 m/s

k=10−7 m/s

k=10−8 m/s

k=10−9 m/s

3.4×10−4 s

3.4×10−5 s

3.4×10−6 s

3.4×10−7 s

3.4×10−8 s

3.4×10−9 s

3.4×10−10 s

10−5 s

10−6 s

10−6 s

10−7 s

10−7 s

10−7 s

10−7 s

5m

10−4 s

10−5 s

10−5 s

10−5 s

10−5 s

10−5 s

10−5 s

5m

Table A1 Minimum time steps used in models with different permeability coefficients and element sizes.

5×10−5 s

10−5 s

10−6 s

10−6 s

10−6 s

10−6 s

10−6 s

Theoretical time steps Δt ≤

mj ks + q 2

1m

Real time steps

1m

Uniform Excitation

fx =2sin(4 πt) 5m

1m

5m

Descriptions: 1. Three models are imposed same sine loading fx=2sin(4πt); 2. Both of the first and second model only have one element in size of 1 m×1 m and 5 m×5 m, respectively. The third model has 25 elements and each element is in the size of 1 m×1 m; 3. Each element has four nodes and each node has three degrees of freedom: the horizontal and vertical displacement, and the pore pressure. 4. The bottom of each model is fixed, and surface of each model is free drainaged boundary. 5. Density of the mixture=2 kg/m3; density of fluid=1 kg/m3; shear modulus=6×104 kPa; Fluid-solid combined bulk modulus =2.2×106 kPa; Permeability=10−3-10−9 m/s. Based on these parameters, the scalar values of the coefficients used in the Eq. (A13) respectively are m=0.222222, k=320000, q=0.166667, s=5.05051×10−8, j=6.79579×10−5.

374

Soil Dynamics and Earthquake Engineering 97 (2017) 364–376

C. Xu et al.

From the Table A1, the real computational time steps are closed to the theory time steps in the conditions of higher permeability, such as 10−3 m/s to 10−5 m/s. However, as the permeability is less than 10−6 m/s, the real time steps do not change with the decrease of the permeability. The models with bigger elements sizes have larger time steps. The time steps of model in 5 m×5 m elements size are 10 times larger than the model in 1 m×1 m elements size. This illustrate that the element sizes affect significantly stabilization of the method. In addition, the model of 25 elements has the smaller time step to compare with the model of one element. It needs to be pointed out that the matrices Ml, Kd, Q, J and Sl are all treated as scalars during the stability analyses of the proposed algorithm in this section. However, an analytical expression for the time step length is generally difficult to be obtained for multi-degrees of freedom systems. Therefore, the real time step length is usually determined by empirical judgment or trial-and-error method instead. A.2. Theoretical accuracy The theoretical accuracy of the completely explicit finite element method is discussed in this section. The matrices of M, Kd, Q, J and S are again regarded as scalars as shown in the previous section. In the proposed method, three different explicit methods, the central, the backward and the forward difference methods are used, as shown in Eqs. (13), (15) and (16), respectively, and the scalar forms of Eqs. (13), (15) and (16) are listed below:

⎞ ⎛ q k 1 uk +1 = ⎜2 − Δt 2⎟ uk − uk −1 + pk Δt 2 + fuk Δt 2 ⎝ m ⎠ m m

(A14)

uk̇ = (uk − uk−1)/ Δt

(A15)

⎛ j ⎞ q 1 pk +1 = ⎜1 − Δt ⎟ pk − uk̇ Δt + fqk Δt ⎝ s ⎠ s s

(A16)

a) For Eq. (A14), the local truncation error is shown as follows

⎡⎛ ⎤ ⎞ q k 1 T(u) k +1 = uk +1 − ⎢ ⎜2 − Δt 2⎟ uk − uk −1 + pk Δt 2 + fuk Δt 2⎥ m ⎠ m m ⎣⎝ ⎦

(A17)

Expanding uk+1 and uk−1 in Eq. (A17) in Taylor's series, the local truncation error of Eq. (A14) can be rewritten as follows:

T(u) k +1 =

Δt 2 [muk̈ + kuk − qpk − fuk ] + Ο (Δt 3) m

(A18)

b) For Eq. (A15) after rewritten as Eq. (A19), the local truncation error is shown in Eq. (A20)

uk +1 = uk + uk̇ +1 Δt

(A19)

T(u) k +1 = uk +1 − (uk + uk̇ +1 Δt )

(A20)

Moreover, the displacement uk can be expanded in Taylor's series as follows:

uk = uk +1 − uk̇ +1 Δt +

1 uk̈ +1 Δt 2 + O (Δt 3) 2

(A21)

Considering the continuity condition uk̈ +1 = uk̈ + O (Δt ), the local truncation error of Eq. (A15) is:

1 1 T(u) k +1 = − uk̈ +1 Δt 2 − O (Δt 3) = − uk̈ Δt 2 − O (Δt 3) 2 2

(A22)

c) For Eq. (A16), the local truncation errors is

⎤ ⎡⎛ j ⎞ q 1 T( p ) k +1 = pk +1 − ⎢ ⎜1 − Δt ⎟ pk − uk̇ Δt + fqk Δt ⎥ s ⎠ s s ⎦ ⎣⎝

(A23)

The pore pressure pk+1 in Eq. (A23) is expanded by the Taylor series, and the local truncation error of Eq. (A16) is:

T( p ) k +1 =

Δt 1 [spk̇ + jpk + quk̇ − fqk ] + pk̈ Δt 2 + Ο (Δt 3) s 2

(A24) 1 Ο (Δt 3),− 2 uk̈ Δt 2

1 p ̈ Δt 2 , 2 k

From Eqs. (A18), (A22) and (A24), it can be found that the dominant terms of the local truncation errors are and respectively. Therefore, Eq. (A14) has second-order accuracy, and Eqs. (A15) and (A16) have first-order accuracy, respectively, which results in first-order accuracy for the proposed method.

[3] Zienkiewicz OC, Chang CT, Bettess P. Drained, undrained, consolidating, and dynamic behavior assumptions in soils: limits of validity. Geotechnique 1980;30(4):385–95. [4] Zienkiewicz OC. Basic formulation of static and dynamic behaviours of soil and other porous media. J Appl Math Mech 1982;3(4):457–68. [5] Li XJ, Liao ZP, Du XL. An explicit finite difference method for viscoelastic dynamic problem. Earthq Eng Eng Vib 1992;12(4):74–9.

References [1] Biot MA. General solutions of the equations of elasticity and consolidation for a porous material. J Appl Mech 1956;23(1):91–6. [2] Biot MA. Theory of propagation of elastic waves in a fluid‐saturated porous solid. I. Low‐frequency range. J Acoust Soc Am 1956;28(2):168–78.

375

Soil Dynamics and Earthquake Engineering 97 (2017) 364–376

C. Xu et al.

[22] Markert B, Heider Y, Ehlers W. Comparison of monolithic and splitting solution schemes for dynamic porous media problems. Int J Numer Methods Eng 2010;82(11):1341–83. [23] Zienkiewicz OC, Chan AHC, Pastor M, et al. Computational geomechanics. Chichester: Wiley; 1999. [24] Hinton E, Rock T, Zienkiewicz OC. A note on mass lumping and related processes in the finite element method. Earthq Eng Struct Dyn 1976;4(3):245–9. [25] Surana KS. Lumped mass matrices with non− zero inertia for general shell and axisymmetric shell elements. Int J Numer Methods Eng 1978;12(11):1635–50. [26] Simon BR, Zienkiewicz OC, Paul DK. An analytical solution for the transient response of saturated porous elastic solids. Int J Numer Anal Methods Geomech 1984;8(4):381–98. [27] Zienkiewicz OC, Shiomi T. Dynamic behavior of saturated porous media; the generalized Biot formulation and its numerical solution. Int J Numer Anal Methods Geomech 1984;8(1):71–96. [28] Wen F, Jeng DS, Wang JH, Zhou XL. Numerical modeling of response of a saturated porous seabed around an offshore pipeline considering non-linear wave and current interaction. Appl Ocean Res 2012;35:25–37. [29] Zienkiewicz OC, Leung KH, Hinton E, et al. Liquefaction and permanent deformation under dynamic conditions-numerical solution and constitutive relations. In: Pande GN, Zienkiewicz OC, editors. Soil mechanics, transient and cyclic loads. London: Wiley; 1982. p. 71–103. [30] Ye JH, Jeng DS. Response of porous seabed to nature loadings-waves and currents. J Eng Mech 2012;138(6):601–13. [31] Ye JH, Jeng DS. Effects of bottom shear stresses on the wave-induced dynamic response in a porous seabed: poro-wssi (shear) model. Acta Mech Sin 2011;27(6):898–910. [32] Liu GL, Song EX. Visco-elastic transmitting boundary for numerical analysis of infinite saturated soil foundation. Chin J Geotech Eng 2006;28(12):2128–33. [33] Forcellini D, Tarantino AM. Countermeasures assessment of liquefaction-induced lateral deformation in a slope ground system. J Eng 2013;2013:1–9. [34] Yang Z, Elgamal A. Influence of permeability on liquefaction-induced shear deformation. J Eng Mech 2002;128(7):720–9. [35] Yang Z. Numerical modeling of earthquake site response including dilation and liquefaction [Ph.D. dissertation]. New York: Dept. of Civil Engineering and Engineering Mechanics, Columbia University; 2000. p. 47–59. [36] Zhu JQ. On numerical stability in the dynamic analyses of structures. Chin J Theor Appl Mech 1983;4:388–95. [37] Li XJ, Liao ZP. The analysis of the properties of an explicit difference method for solving the nonlinear structural dynamic equations. Eng Mech 1993;10(3):141–7. [38] Pena JM. Characterizations and stable tests for the Routh–Hurwitz conditions and for total positivity. Line Algebra Appl 2004;393:319–32.

[6] Du XL, Wang JT. An explicit difference formulation of dynamic response calculation of elastic structure with damping. Eng Mech 2000;17(5):37–43. [7] Wang JT, Du XL. An explicit difference method for dynamic analysis of a structure system with damping. Eng Mech 2002;19(3):109–12. [8] Zhao CG, Li WH, Wang JT. 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(3):1169–81. [9] Wang JT, Zhang CH, Du XL. An explicit integration scheme for solving dynamic problems of solid and porous media. Int J Earthq Eng 2008;12(2):293–311. [10] Simon BR, Wu JSS, Zienkiewicz OC, Paul DK. Evaluation of u-w and u-π finite element methods for the dynamic response of saturated porous media using onedimensional models. Int J Numer Anal Methods Geomech 1986;10(5):461–82. [11] Prevost JH. Wave propagation in fluid-saturated porous media: an efficient finite element procedure. Soil Dyn Earthq Eng 1985;4(4):183–201. [12] Felippa CA, Park KC. Staggered transient analysis procedures for coupled mechanical systems: formulation. Comput Methods Appl Mech Eng 1980;24(1):61–111. [13] Park KC. Stabilization of partitioned solution procedure for pore fluid-soil interaction analysis. Int J Numer Methods Eng 1983;19(11):1669–73. [14] Zienkiewicz OC, Paul DK, Chan AHC. Unconditionally stable staggered solution procedure for soil-pore fluid interaction problems. Int J Numer Methods Eng 1988;26(5):1039–55. [15] Huang MS, Zienkiewicz OC. New unconditionally stable staggered solution procedures for coupled soil-pore fluid dynamic problems. Int J Numer Methods Eng 1998;43(6):1029–52. [16] Babuska I. Error-bounds for finite element method. Numer Math 1971;16(4):322–33. [17] Brezzi F. On the existence, uniqueness and approximation of saddle-point problems arising from Lagrangian multipliers. Revue française d′automatique, informatique, recherche opérationnelle. Anal Numér 1974;8(2):129–51. [18] Zienkiewicz OC, Huang MS, Wu J, Wu SM. A new algorithm for the coupled soilpore fluid problem. Shock Vib 1993;1(1):3–14. [19] Huang MS, Wu SM, Zienkiewicz OC. Incompressible or nearly incompressible soil dynamic behavior-a new staggered algorithm to circumvent restrictions of mixed formulation. Soil Dyn Earthq Eng 2001;21(2):169–79. [20] Park T, Tak M. A new coupled analysis for nearly incompressible and impermeable saturated porous media on mixed finite element method: I. Proposed method. KSCE J Civ Eng 2010;14(1):7–16. [21] Soares D, Rodrigues GG, Gonçalves KA. 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(5):387–94.

376