A Wavelet Galerkin Finite-Element Method for the Biot Wave Equation in the Fluid-Saturated Porous Medium

Abstract

A wavelet Galerkin finite-element method is proposed by combining the wavelet analysis with traditional finite-element method to analyze wave propagation phenomena in fluid-saturated porous medium. The scaling functions of Daubechies wavelets are considered as the interpolation basis functions to replace the polynomial functions, and then the wavelet element is constructed. In order to overcome the integral difficulty for lacking of the explicit expression for the Daubechies wavelets, a kind of characteristic function is introduced. The recursive expression of calculating the function values of Daubechies wavelets on the fraction nodes is deduced, and the rapid wavelet transform between the wavelet coefficient space and the wave field displacement space is constructed. The results of numerical simulation demonstrate that the method is effective.

1. Introduction

The fluid-saturated porous medium is modeled as a two-phase system consisting of a solid and a fluid phase. It is assumed that the solid phase is homogenous, isotropic, elastic frame and the fluid phase is viscous, compressible, and filled with the pore space of solid frame. Compared with the single-phase medium theory, fluid-saturated porous medium theory can describe the formation underground more precisely and the fluid-saturated porous medium elastic wave equation can bring more lithology information than ever. For these reasons, fluid-saturated porous medium theory can be used widely in geophysics exploration and engineering surveying.

In 1956, a theory was developed for the propagation of stress waves in a porous elastic solid containing compressible viscous fluid by Biot [1, 2]. Biot described the Second-Kind P wave in fluid-saturated porous medium firstly. Since then, many researchers paid their attention to the propagation characters of elastic wave in saturated porous medium and obtained many achievements [3, 4]. Complicated equations given in Biot dynamic theory can be solved by analytical methods with some simple boundary conditions. Most dynamic problems in fluid-saturated porous medium are solved using numerical methods, especially using finite-element method. Ghaboussi and Wilson [5] first proposed a multidimensional finite element numerical scheme to solve the linear coupled governing equations. Prevose [6] proposed an efficient finite element procedure to analyze wave propagation phenomena in fluid saturated porous medium and presented some numerical results which demonstrate the versatility of the proposed procedure. Simon et al. [7, 8] presented an analytical solution for a transient analysis of a one-dimensional column of a fluid saturated porous elastic solid and presented a comparison of this exact closed-form solution with finite-element method for several transient problems in porous media. Yazdchi et al. [9, 10] combined the finite element method with the boundary element method and the infinite element method, constructed the finite-infinite element method and the finite-boundary element method to deal with the two-phase model in lateral extensive field and obtained better result. Zhao et al. [11] proposed an explicit finite element method for Biot dynamic formulation in fluid-saturated porous medium. It does not need to assemble a global stiffness matrix and solve a set of linear equations in each time step by using the decoupling-technique. For the problem of local high gradient, finite element method improves the calculation precision by employing the higher-order polynomial or the denser mesh. However, the increment of polynomial order and mesh knots inevitably needs more computational work. Meanwhile, the condition of numerical dissipation will limit the frequency range that can be obtained. To overcome these disadvantages, wavelet analysis is introduced to the finite-element method in this paper. As a new method, the development of wavelet analysis is recent fairly in many fields. Its desirable advantages are the multiresolution analysis property and various basis functions for structure analysis. According to different requirement, the corresponding scaling functions and wavelet functions can be adopted to improve the numerical calculation precision. Especially, those wavelets with compactly supported property and orthogonality, such as Daubechies wavelets, can play an important role in many problems [12]. Because of the compactly supported property, if the Daubechies wavelets are considered as the interpolation functions of the finite element method, the coefficient matrices obtained are sparse matrices and their condition number can be proved independent of the dimension [13]. Moreover, a new method could be provided because of the existence of various basis functions, which can increase the resolution without changing mesh.

In this paper, the wavelet Galerkin finite element method is applied to the direct simulation of the wave equation in the fluid-saturated porous medium. The scaling functions of Daubechies wavelets are considered as the interpolation basis functions instead of the polynomial functions and the wavelet element is constructed. Because a kind of characteristic function is introduced, the integral difficulty for lacking of the explicit expression for the Daubechies wavelets is solved. Based on the recursive expression of calculating the function values of Daubechies wavelets on the fraction nodes, the rapid wavelet transform between the wavelet coefficient space and the wave field displacement space is constructed and reduces the computational cost. The results of numerical simulation demonstrate the method is effective.

2. Wavelet Galerkin Finite-Element Method

2.1. Wavelet Galerkin Finite-Element Method

For purpose of constructing the wavelet Galerkin finite element method, we consider a typical boundary value problem:

𝐿(𝑢,𝑣)=𝑓,𝐵(𝑢,𝑣)|_𝜕Ω=𝑔, (2.1)(2.2) where𝐿(⋅,⋅) is differential operator,𝐵(⋅,⋅) is boundary operator, 𝑢,𝑣 are the unknown functions in the solving domain Ω, and 𝜕Ω is the boundary.

Supposing 𝑢,𝑣 are the exact solutions of (2.1) and (2.2), then one gets

𝐿(𝑢,𝑣)−𝑓≡0, (2.3) and if 𝐿(𝑢,𝑣)and 𝑓 are continuous, (2.3) is equal to

_Ω(𝐿(𝑢,𝑣)−𝑓)𝜙_𝑘d𝑥d𝑦=0. (2.4)

In fact, because of the derivation of one-dimensional wavelet basis element facilitates a straightforward discussion of multidimensional tensor product wavelet basis element and multiresolution analysis property of wavelet function [12], the functions 𝑢,𝑣 can be assumed to consist of a superposition of scaling functions at 𝑗level and wavelet functions at the same and higher levels:

𝑢(𝑥,𝑦)=𝑈₁(𝑥)𝑈₂(𝑦),𝑣(𝑥,𝑦)=𝑉₁(𝑥)𝑉₂(𝑦), (2.5)(2.6) where

𝑈₁(𝑥)=_𝑘𝑎_𝑗,𝑘𝜙_𝑗𝑘(𝑥)+_{𝑖≥𝑗,𝑘}𝑎_𝑖,𝑘𝜓_𝑖𝑘𝑈(𝑥),₂(𝑦)=_𝑘𝑏_𝑗,𝑘𝜙_𝑗𝑘(𝑦)+_{𝑖≥𝑗,𝑘}𝑏_𝑖,𝑘𝜓_𝑖𝑘𝑉(𝑦),₁(𝑥)=_𝑘𝑐_𝑗,𝑘𝜙_𝑗𝑘(𝑥)+_{𝑖≥𝑗,𝑘}𝑐_𝑖,𝑘𝜓_𝑖𝑘𝑉(𝑥),₂(𝑦)=_𝑘𝑑_𝑗,𝑘𝜙_𝑗𝑘(𝑦)+_{𝑖≥𝑗,𝑘}𝑑_𝑖,𝑘𝜓_𝑖𝑘(𝑦). (2.7)

Upon substituting (2.4) and (2.5) into (2.3), we can obtain an equation system of wavelet coefficients, whose coefficient matrix consists of the following integrals:

_Ω𝜙_𝑗𝑞𝜓_𝑖𝑟,_Ω𝜙_𝑗𝑞𝜙_𝑗𝑟,_Ω𝜓_𝑖𝑞𝜓_𝑙𝑟,_Ω𝜙^(𝑚)_𝑗𝑞𝜓^(𝑛)_𝑖𝑟,_Ω𝜙^(𝑚)_𝑗𝑞𝜙^(𝑛)_𝑗𝑟,_Ω𝜓^(𝑚)_𝑖𝑞𝜓^(𝑛)_𝑙𝑟. (2.8)

In conventional finite element method, these integrals would be calculated by Gauss quadrature formulae. However, it is not feasible for most wavelet functions. In many cases, there is no explicit expression for the function, in this paper, we choose the Daubechies wavelet as the basis function, and they cannot be integrated numerically due to their unusual smoothness characteristics. Moreover, the wavelet function is defined in terms of scaling function, so these integrals can be rewritten in terms of scaling function alone.

Define the connection coefficients [14–16]:

Γ^0,0_𝑝,𝑟=_ΩΓ𝜙(𝑥−𝑝)𝜙(𝑥−𝑟)d𝑥,_{𝑚,𝑛𝑝,𝑟}=_Ωd^(𝑚)𝜙d𝑥^𝑚d(𝑥−𝑝)^(𝑛)𝜙d𝑥^𝑛=(𝑥−𝑟)d𝑥,_Ω𝜙^(𝑚)(𝑥−𝑝)𝜙^(𝑛)(𝑥−𝑟)d𝑥. (2.9)

Once these integrals can be calculated, all the integrals in (2.8) can be obtained and eventually construct the stiffness matrix and load matrix of wavelet Galerkin finite element method.

2.2. The Calculation of Wavelet Connection Coefficients

From what has been discussed earlier, the quality matrix, stiffness matrix, and the load matrix are composed of the integral values of Daubechies wavelets. However, it is well known that Daubechies wavelets have no explicit expression. In order to solve this problem, a kind of characteristic function is introduced:

𝜒_[0,1](𝑥)=10≤𝑥≤1,0otherwise. (2.10) Set 𝜉=2𝑥then 𝜒_[0,1](𝜉2)=10≤𝜉≤2,0otherwise. (2.11)

So the trivial two-scale equation of characteristic function is obtained:

𝜒_[0,1](𝜉2)=𝜒_[0,1](𝜉)+𝜒_[1,2](𝜉)=𝜒_[0,1](𝜉)+𝜒_[0,1](𝜉−1). (2.12)

Set 𝜏^0,0_𝑘,𝑠=_𝑅𝜒_[0,1](𝑥)𝜙(𝑥−𝑘)𝜙(𝑥−𝑠)d𝑥. (2.13)

Substituting 𝜙(𝑥)=∑_𝑘𝑎_𝑘𝜙(2𝑥−𝑘) into (2.13), one obtains

𝜏^0,0_𝑘,𝑠=_𝑅𝜒_[0,1](𝑥)_𝑙𝑎_𝑙𝜙(2(𝑥−𝑘)−𝑙)_𝑚𝑎_𝑚=𝜙(2(𝑥−𝑠)−𝑚)d𝑥_𝑅𝜒_[0,1](𝑥)_𝑙𝑚𝑎_𝑙𝑎_𝑚=1𝜙(2𝑥−2𝑘−𝑙)𝜙(2𝑥−2𝑠−𝑚)d𝑥2_𝑅𝜒_[0,1](𝜉2)_𝑙𝑚𝑎_𝑙𝑎_𝑚=1𝜙(𝜉−2𝑘−𝑙)𝜙(𝜉−2𝑠−𝑚)d𝜉2_𝑙𝑚𝑎_𝑙𝑎_𝑚𝑅(𝜒_[0,1](𝜉)+𝜒_[0,1]=1(𝜉−1))𝜙(𝜉−2𝑘−𝑙)𝜙(𝜉−2𝑠−𝑚)d𝜉2_𝑙𝑚𝑎_𝑙𝑎_𝑚𝑅𝜒_[0,1]+1(𝜉)𝜙(𝜉−2𝑘−𝑙)𝜙(𝜉−2𝑠−𝑚)d𝜉2_𝑙𝑚𝑎_𝑙𝑎_𝑚𝑅𝜒_[0,1]=1(𝜉−1)𝜙(𝜉−2𝑘−𝑙)𝜙(𝜉−2𝑠−𝑚)d𝜉2_𝑙𝑚𝑎_𝑙𝑎_𝑚𝜏^0,0_{2𝑘+𝑙,2𝑠+𝑚}+12_𝑙𝑚𝑎_𝑙𝑎_𝑚𝜏^0,0_{2𝑘+𝑙−1,2𝑠+𝑚−1}=12_𝑝𝑟(𝑎_𝑝−2𝑘𝑎_𝑟−2𝑠+𝑎_{𝑝−2𝑘+1}𝑎_{𝑟−2𝑠+1})𝜏^0,0_𝑝,𝑟. (2.14)

It is not difficult to show that we will require the solution of an eigenvalue problem having the form

𝜏^0,0_𝑘,𝑠=12𝐴^′𝜏^′0,0_𝑝,𝑟,2𝑁−1≤𝑘,𝑠≤0,2𝑁−1≤𝑝,𝑟≤0, (2.15) where 𝐴^′ is a (2𝑁−1)×(2𝑁−1) partitioned matrix, each submatrix is also a (2𝑁−1)×(2𝑁−1) matrix, in which𝑎^′_𝑠𝑟=𝑎_𝑝−2𝑘𝑎_𝑟−2𝑠+𝑎_{𝑝−2𝑘+1}𝑎_{𝑟−2𝑠+1}.

Considering the requirement of numerical simulation set

1≤𝑖=𝑠+8≤8,1≤𝑗=𝑟+8≤8,1≤𝑚=𝑘+8≤8,1≤𝑛=𝑝+8≤8, (2.16) then 𝐴^′_𝑘𝑝 is changed to𝐴^′_𝑚𝑛, in which 𝑎^′_𝑖𝑗=𝑎_{𝑛−2𝑚+8}𝑎_{𝑗−2𝑖+8}+𝑎_{𝑛−2𝑚+9}𝑎_{𝑗−2𝑖+9}.

However, the eigenvalue problem does not uniquely define the solution, it is essential to introduce an additional condition to define the solution uniquely.

It is well known that the Daubechies wavelets satisfy

1=_𝑘𝜙(𝑥−𝑘). (2.17)

By multiplying (2.17) by itself, and subsequently multiplying the product by the characteristic function𝜒_[0,1](𝑥), one obtains

1=_𝑘𝑙𝜒𝜙(𝑥−𝑘)𝜙(𝑥−𝑙),_[0,1](𝑥)=_𝑘𝑙𝜒_[0,1](𝑥)𝜙(𝑥−𝑘)𝜙(𝑥−𝑙). (2.18)

Now, a single integration yields a first normalization condition:

1=_𝑘𝑙𝜏^0,0_𝑘,𝑙. (2.19)

So, the unique solution of the eigenvalue problem is defined.

The same step can be followed to calculate

𝜏^𝑚,𝑛_𝑘,𝑠=_𝑅𝜒_[0,1](𝑥)𝜙^(𝑚)(𝑥−𝑘)𝜙^(𝑛)(𝑥−𝑠)d𝑥. (2.20)

Substituting 𝜙^(𝑚)(𝑥)=2^𝑚∑_𝑙𝑎_𝑙𝜙^(𝑚)(2𝑥−𝑘) and 𝜙^(𝑛)(𝑥)=2^𝑛∑_𝑞𝑎_𝑞𝜙^(𝑛)(2𝑥−𝑘) into (2.20), one gets

𝜏^𝑚,𝑛_𝑘,𝑠=_𝑅𝜒_[0,1](𝑥)2^𝑚+𝑛_𝑙𝑎_𝑙𝜙^(𝑚)(2(𝑥−𝑘)−𝑙)_𝑞𝑎_𝑞𝜙^(𝑛)(2(𝑥−𝑠)−𝑞)d𝑥=2^𝑚+𝑛_𝑅𝜒_[0,1](𝑥)_𝑙𝑞𝑎_𝑙𝑎_𝑞𝜙^(𝑚)(2𝑥−2𝑘−𝑙)𝜙^(𝑛)(2𝑥−2𝑠−𝑞)d𝑥=2^{𝑚+𝑛−1}_𝑅𝜒_[0,1](𝜉2)_𝑙𝑞𝑎_𝑙𝑎_𝑞𝜙^(𝑚)(𝜉−2𝑘−𝑙)𝜙^(𝑛)(𝜉−2𝑠−𝑞)d𝜉=2^{𝑚+𝑛−1}_𝑙𝑞𝑎_𝑙𝑎_𝑞𝑅(𝜒_[0,1](𝜉)+𝜒_[0,1](𝜉−1))𝜙^(𝑚)(𝜉−2𝑘−𝑙)𝜙^(𝑛)(𝜉−2𝑠−𝑞)d𝜉=2^{𝑚+𝑛−1}_𝑙𝑞𝑎_𝑙𝑎_𝑞𝑅𝜒_[0,1](𝜉)𝜙^(𝑚)(𝜉−2𝑘−𝑙)𝜙^(𝑛)(𝜉−2𝑠−𝑞)d𝜉+2^{𝑚+𝑛−1}_𝑙𝑞𝑎_𝑙𝑎_𝑞𝑅𝜒_[0,1](𝜉−1)𝜙^(𝑚)(𝜉−2𝑘−𝑙)𝜙^(𝑛)(𝜉−2𝑠−𝑞)d𝜉=2^{𝑚+𝑛−1}_𝑙𝑞𝑎_𝑙𝑎_𝑞𝜏^𝑚,𝑛_{2𝑘+𝑙,2𝑠+𝑞}+2^{𝑚+𝑛−1}_𝑙𝑞𝑎_𝑙𝑎_𝑞𝜏^𝑚,𝑛_{2𝑘+𝑙−1,2𝑠+𝑞−1}=2^{𝑚+𝑛−1}_𝑝𝑟(𝑎_𝑝−2𝑘𝑎_𝑟−2𝑠+𝑎_{𝑝−2𝑘+1}𝑎_{𝑟−2𝑠+1})𝜏^𝑚,𝑛_𝑝,𝑟, (2.21) namely,

𝜏^𝑚,𝑛_𝑘,𝑠=2^{𝑚+𝑛−1}𝐴^′𝜏^𝑚,𝑛_𝑝,𝑟. (2.22)

The polynomial reproducing property is employed to construct the additional condition:

𝑥^𝑚=_𝑘𝑝_𝑘𝑥𝜙(𝑥−𝑘),_𝑛=_𝑙𝑝_𝑙𝜙(𝑥−𝑙), (2.23)(2.24)

Explicit form for calculating the coefficients 𝑝_𝑘𝑝_𝑙 can be found in [17].

By differentiating (2.23) 𝑚 times, one obtains

𝑚!=_𝑘𝑝_𝑘𝜙^(𝑚)(𝑥−𝑘). (2.25)

By differentiating (2.24) 𝑛 times, one gets

𝑛!=_𝑙𝑝_𝑙𝜙^(𝑛)(𝑥−𝑙). (2.26)

However (2.25) can be multiplied by (2.26), and subsequently multiplying the product by the characteristic function𝜒_[0,1](𝑥),

𝑚!𝑛!𝜒_[0,1](𝑥)=_𝑘𝑙𝜒_[0,1](𝑥)𝑝_𝑘𝑝_𝑙𝜙^(𝑚)(𝑥−𝑘)𝜙^(𝑛)(𝑥−𝑙). (2.27)

By integrating (2.27), one obtains the additional condition.

𝑚!𝑛!=_𝑘𝑙𝑝_𝑘𝑝_𝑙𝜏^𝑚,𝑛_𝑘,𝑙. (2.28)

Then, the unique solution of the eigenvalue problem is defined.

3. Wavelet Galerkin Finite-Element Solution of 1D Elastic Wave Equation in Fluid-Saturated Porous Medium

From the Biot theory, the 1D differential equation governing wave propagation in the fluid-saturated porous medium, without fluid viscosity, can be expressed as

𝜕𝜕𝑥((𝜆+2𝜇+𝛼𝑀)𝜕𝑢𝜕𝜕𝑥)+𝜕𝑥(𝛼𝑀𝜕𝜔̈𝜕𝑥)=𝜌𝑢+𝜌_𝑓𝜔̈−𝑓₁,𝜕𝜕𝑥(𝛼𝑀𝜕𝑢𝜕𝜕𝑥)+𝜕𝑥(𝑀𝜕𝜔𝜕𝑥)=𝜌_𝑓𝑢̈+𝑚𝜔̈−𝑓₂, (3.1) where 𝑢 is the solid displacement and𝜔is the relative fluid to solid displacement. 𝛽 is the porosity, 𝜌=(1−𝛽)𝜌_𝑠+𝛽𝜌_𝑓 is the bulk density of solid-fluid mixture, and 𝜌_𝑠and 𝜌_𝑓 are the densities of solid and fluid, respectively. Also 𝑡 is time and 𝜆_𝑏,𝜇 are the Lame coefficients,𝜆=𝜆_𝑏+𝛼²𝑀, where 𝛼 is the effective stress parameter and 𝑀 is the compressibility of pore fluid. 𝛼=1−𝐾_𝑏/𝐾_𝑠, 𝑀=𝐾_𝑠/[𝛼+𝛽(𝐾_𝑠/𝐾_𝑓−1)] where 𝐾_𝑠,𝐾_𝑓,𝐾_𝑏 are the bulk change modulus of the solid, fluid, and skeleton, respectively. Moreover𝐾_𝑏=𝜆_𝑏+2𝜇/3,𝑚=𝜌_𝑓/𝛽, Finally 𝑓is seismic focus, and 𝑓₁=(1−𝛽)𝑓, 𝑓₂=𝛽(2𝛽−1)𝑓.

Multiplying both sides of the fluid-saturated porous medium wave equation by the Daubechies wavelets basis function𝜙_𝑗𝑘(𝑥)=2^𝑗/2𝜙(2^𝑗𝑥−𝑘), and integrating them at[0,𝐿], we can get

∫^𝐿₀(𝜕𝜕𝑥((𝜆+2𝜇+𝛼𝑀)𝜕𝑢𝜕𝜕𝑥)+𝜕𝑥(𝛼𝑀𝜕𝜔𝜕𝑥))𝜙_𝑗𝑘(𝑥)d𝑥=∫^𝐿₀(𝜌𝑢̈+𝜌_𝑓𝜔̈−𝑓₁)𝜙_𝑗𝑘(𝑥)d𝑥_𝐿0(𝜕𝜕𝑥(𝛼𝑀𝜕𝑢𝜕𝜕𝑥)+𝜕𝑥(𝑀𝜕𝜔𝜕𝑥))𝜙_𝑗𝑘(𝑥)d𝑥=^𝐿₀(𝜌_𝑓𝑢̈+𝑚𝜔̈−𝑓₂)𝜙_𝑗𝑘(𝑥)d𝑥. (3.2)

By using integration by part

(𝜆+2𝜇+𝛼𝑀)𝜕𝑢𝜙𝜕𝑥^𝑗𝑘(𝑥)|^𝐿₀−^𝐿₀(𝜆+2𝜇+𝛼𝑀)𝜕𝑢𝜕𝑥𝜕𝜙_𝑗𝑘(𝑥)𝜕𝑥d𝑥+𝛼𝑀𝜕𝜔𝜙𝜕𝑥^𝑗𝑘(𝑥)|^𝐿₀−^𝐿₀𝛼𝑀𝜕𝜔𝜕𝑥𝜕𝜙_𝑗𝑘(𝑥)𝜕𝑥d𝑥=^𝐿₀(𝜌𝑢̈+𝜌_𝑓𝜔̈−𝑓₁)𝜙_𝑗𝑘(𝑥)d𝑥,𝛼𝑀𝜕𝑢𝜙𝜕𝑥^𝑗𝑘(𝑥)|^𝐿₀−^𝐿₀𝛼𝑀𝜕𝑢𝜕𝑥𝜕𝜙_𝑗𝑘(𝑥)𝜕𝑥d𝑥+𝑀𝜕𝜔𝜙𝜕𝑥^𝑗𝑘(𝑥)|^𝐿₀−^𝐿₀𝑀𝜕𝜔𝜕𝑥𝜕𝜙_𝑗𝑘(𝑥)=𝜕𝑥d𝑥_𝐿0(𝜌_𝑓𝑢̈+𝑚𝜔̈−𝑓₂)𝜙_𝑗𝑘(𝑥)d𝑥. (3.3)(3.4)

Set

𝑢(𝑥,𝑡)=⁰_{𝑙=2−2𝑁−2𝑗𝐿}𝑎_𝑙(𝑡)𝜙_𝑗𝑙(𝑥),𝜔(𝑥,𝑡)=⁰_{𝑙=2−2𝑁−2𝑗𝐿}𝑏_𝑙(𝑡)𝜙_𝑗𝑙(𝑥). (3.5)

Upon substituting (3.5) into (3.3) and (3.4), one gets

(𝜆+2𝜇+𝛼𝑀)𝜕𝑢𝜙𝜕𝑥^𝑗𝑘(𝑥)|^𝐿₀+𝛼𝑀𝜕𝜔𝜙𝜕𝑥^𝑗𝑘(𝑥)|^𝐿₀−^𝐿₀((𝜆+2𝜇+𝛼𝑀)⁰_{𝑙=2−2𝑁−2𝑗𝐿}𝑎_𝑙(𝑡)𝜕𝜙_𝑗𝑙(𝑥)𝜕𝑥𝜕𝜙_𝑗𝑘(𝑥)𝜕𝑥+𝛼𝑀⁰_{𝑙=2−2𝑁−2𝑗𝐿}𝑏_𝑙(𝑡)𝜕𝜙_𝑗𝑙(𝑥)𝜕𝑥𝜕𝜙_𝑗𝑘(𝑥)=𝜕𝑥)d𝑥_𝐿0(𝜌⁰_{𝑙=2−2𝑁−2𝑗𝐿}𝑎_𝑙^′′(𝑡)𝜙_𝑗𝑙(𝑥)+𝜌_𝑓⁰_{𝑙=2−2𝑁−2𝑗𝐿}𝑏_𝑙^′′(𝑡)𝜙_𝑗𝑙(𝑥)−𝑓₁)𝜙_𝑗𝑘(𝑥)d𝑥,𝛼𝑀𝜕𝑢𝜙𝜕𝑥^𝑗𝑘(𝑥)|^𝐿₀+𝑀𝜕𝜔𝜙𝜕𝑥^𝑗𝑘(𝑥)|^𝐿₀−^𝐿₀(𝛼𝑀⁰_{𝑙=2−2𝑁−2𝑗𝐿}𝑎_𝑙(𝑡)𝜕𝜙_𝑗𝑙(𝑥)𝜕𝑥𝜕𝜙_𝑗𝑘(𝑥)𝜕𝑥+𝑀⁰_{𝑙=2−2𝑁−2𝑗𝐿}𝑏_𝑙(𝑡)𝜕𝜙_𝑗𝑙(𝑥)𝜕𝑥𝜕𝜙_𝑗𝑘(𝑥)𝜕𝑥)d𝑥=^𝐿₀(𝜌_𝑓⁰_{𝑙=2−2𝑁−2𝑗𝐿}𝑎_𝑙^′′(𝑡)𝜙_𝑗𝑙(𝑥)+𝑚⁰_{𝑙=2−2𝑁−2𝑗𝐿}𝑏_𝑙^′′(𝑡)𝜙_𝑗𝑙(𝑥)−𝑓₂)𝜙_𝑗𝑘(𝑥)d𝑥. (3.6)

By rearranging, (3.6) and become

(𝜆+2𝜇+𝛼𝑀)𝜕𝑢𝜙𝜕𝑥^𝑗𝑘(𝑥)|^𝐿₀+𝛼𝑀𝜕𝜔𝜙𝜕𝑥^𝑗𝑘(𝑥)|^𝐿₀−((𝜆+2𝜇+𝛼𝑀)⁰_{𝑙=2−2𝑁−2𝑗𝐿}𝑎_𝑙(𝑡)+𝛼𝑀⁰_{𝑙=2−2𝑁−2𝑗𝐿}𝑏_𝑙(𝑡))^𝐿₀𝜕𝜙_𝑗𝑙(𝑥)𝜕𝑥𝜕𝜙_𝑗𝑘(𝑥)𝜕𝑥d𝑥=(𝜌⁰_{𝑙=2−2𝑁−2𝑗𝐿}𝑎_𝑙^′′(𝑡)+𝜌_𝑓⁰_{𝑙=2−2𝑁−2𝑗𝐿}𝑏_𝑙^′′(𝑡))^𝐿₀𝜙_𝑗𝑙(𝑥)𝜙_𝑗𝑘(𝑥)d𝑥−𝑓₁^𝐿₀𝜙_𝑗𝑘(𝑥)d𝑥,𝛼𝑀𝜕𝑢𝜙𝜕𝑥^𝑗𝑘(𝑥)|^𝐿₀+𝑀𝜕𝜔𝜙𝜕𝑥^𝑗𝑘(𝑥)|^𝐿₀−(𝛼𝑀⁰_{𝑙=2−2𝑁−2𝑗𝐿}𝑎_𝑙(𝑡)+𝑀⁰_{𝑙=2−2𝑁−2𝑗𝐿}𝑏_𝑙(𝑡))^𝐿₀𝜕𝜙_𝑗𝑙(𝑥)𝜕𝑥𝜕𝜙_𝑗𝑘(𝑥)𝜕𝑥d𝑥=(𝜌_𝑓⁰_{𝑙=2−2𝑁−2𝑗𝐿}𝑎_𝑙^′′(𝑡)+𝑚⁰_{𝑙=2−2𝑁−2𝑗𝐿}𝑏_𝑙^′′(𝑡))^𝐿₀𝜙_𝑗𝑙(𝑥)𝜙_𝑗𝑘(𝑥)d𝑥−𝑓₂^𝐿₀𝜙_𝑗𝑘(𝑥)d𝑥. (3.7)

If select 𝐿=1,𝑗=0, (3.5) become

𝑢(𝑥,𝑡)=⁰_{𝑙=1−2𝑁}𝑎_𝑙(𝑡)𝜙(𝑥−𝑙)𝜔(𝑥,𝑡)=⁰_{𝑙=1−2𝑁}𝑏_𝑙(𝑡)𝜙(𝑥−𝑙) (3.8)

Set

𝐀=(𝑎_1−2𝑁,𝑎_2−2𝑁⋯𝑎₀)𝐁=(𝑏_1−2𝑁,𝑏_2−2𝑁⋯𝑏₀),𝐑=(𝐀,𝐁)^𝑇=(𝑎_1−2𝑁,𝑎_2−2𝑁⋯𝑎₀,𝑏_1−2𝑁,𝑏_2−2𝑁⋯𝑏₀)^𝑇. (3.9)

Then, (3.7) can be changed into an equation system of coefficient𝐑:

𝐌𝐑̈+𝐏𝐑=𝐅+𝐐, (3.10) where

𝐌=𝜌𝐄𝜌_𝑓𝐄𝜌_𝑓𝐅𝐄𝑚𝐄,𝐏=−(𝜆+2𝜇+𝛼𝑀)𝐆−𝛼𝑀𝐆−𝛼𝑀𝐆−𝑀𝐆,𝐅=^𝟏𝐅_𝟐,𝐄=¹⁰𝔄𝔄d𝑥¹₀𝔅𝔄d𝑥⋯¹₀𝜙(𝑥)𝔄d𝑥₁₀𝔄𝔅d𝑥¹₀𝔅𝔅d𝑥⋯¹₀𝜙(𝑥)𝔅d𝑥⋮⋮⋮⋮₁₀𝔄𝜙(𝑥)𝑑𝑥¹₀𝔅𝜙(𝑥)𝑑𝑥⋯¹₀,𝜙(𝑥)𝜙(𝑥)d𝑥𝐆=¹⁰ℭℭd𝑥¹₀𝔇ℭd𝑥⋯¹₀𝜙^′(𝑥)ℭd𝑥₁₀ℭ𝔇d𝑥¹₀𝔇𝔇d𝑥⋯¹₀𝜙^′(𝑥)𝔇d𝑥⋮⋮⋮⋮₁₀ℭ𝜙^′(𝑥)d𝑥¹₀𝔇𝜙^′(𝑥)d𝑥⋯¹₀𝜙^′(𝑥)𝜙^′,𝐅(𝑥)d𝑥_𝟏=𝑓₁¹₀𝔄d𝑥,¹₀𝔅d𝑥,…,¹₀𝜙(𝑥)d𝑥^𝑇,𝐅_𝟐=𝑓₂¹₀𝔄d𝑥,¹₀𝔅d𝑥,…,¹₀𝜙(𝑥)d𝑥^𝑇,𝐐=(𝐐_𝟏,𝐐_𝟐)^𝑇,𝐐₁=((𝜆+2𝜇+𝛼𝑀)𝜕𝑢𝜕𝑥+𝛼𝑀𝜕𝜔𝜕𝑥)(𝔄,𝔅⋯𝜙(𝑥))^𝑇|¹₀,𝐐₂=(𝛼𝑀𝜕𝑢𝜕𝑥+𝑀𝜕𝜔𝜕𝑥)(𝔄,𝔅⋯𝜙(𝑥))^𝑇|¹₀, (3) where 𝔄 denote 𝜙(𝑥−1+2𝑁), 𝔅 denote𝜙(𝑥−2+2𝑁), ℭ denote 𝜙^′(𝑥−1+2𝑁) and 𝔇 denote 𝜙^′(𝑥−2+2𝑁).

Using the second-order center difference to approximate the two derivatives in (3.10), we can obtain

𝐌𝐑^𝑛+1−2𝐑^𝑛+𝐑^𝑛−1(Δ𝑡)²+𝐏𝐑^𝑛=𝐅+𝐐. (3.12)

Arranging (3.12), we have

𝐌𝐑^𝑛+1=(2𝐌−(Δ𝑡)²𝐏)𝐑^𝑛−𝐌𝐑^𝑛−1+(Δ𝑡)²𝐅+(Δ𝑡)²𝐐, (3.13) given the initial conditions:

𝑎_𝑘(0)=𝑏_𝑘(0)=0,𝑎_𝑘(1)=𝑏_𝑘(1)=0. (3.14)

So, we can obtain the wavelet coefficients at each time level by solving (3.13) and (3.14) with some boundary conditions, and then substitute the wavelet coefficients into (3.8), the wave field displacements can be obtained.

4. Rapid Wavelet Transform

In order to obtain the wave field displacements conveniently and quickly, the fast wavelet transform between the wavelet coefficients space and the wave field displacements space is constructed as follows:

𝑈=Φ𝑃, (4.1)𝑈 is the wave field displacement vector,𝑃 is the wavelet coefficient vector, Φ is the wavelet transform matrix.

For the sake of simplicity, take the DB2 wavelet as the example. There are 7 nodes in solution field:

1𝑈=(𝑢(81),𝑢(43),𝑢(81),𝑢(25),𝑢(83),𝑢(47),𝑢(8))^𝑇,𝑃=(𝑝₋₂,𝑝₋₁,𝑝₀)^𝑇,1Φ=𝜙(81+2)𝜙(81+1)𝜙(8)1𝜙(41+2)𝜙(41+1)𝜙(4)3𝜙(83+2)𝜙(83+1)𝜙(8)1𝜙(21+2)𝜙(21+1)𝜙(2)5𝜙(85+2)𝜙(85+1)𝜙(8)3𝜙(43+2)𝜙(43+1)𝜙(4)7𝜙(87+2)𝜙(87+1)𝜙(8). (4.2)

It is important for constructing the fast wavelet transform to solve the function values of the Daubechies wavelets on the fraction nodes. So, the recursive expression of calculating the function values of Daubechies wavelets on the fraction nodes is deduced to save the computational cost.

2Φ(^𝑛𝑖+𝑝2^𝑛2)=𝐴Φ(^𝑛−1𝑖+𝑞2^𝑛−1𝑝)if2⋅2^𝑛2≤1𝐵Φ(^𝑛−1𝑖+𝑞2^𝑛−1𝑝)if2⋅2^𝑛>1 (4.3) in which 𝑖=0,1,…2𝑁−2, 𝑝=1∶2∶2^𝑛−1, 𝑞=𝑝mod2^𝑛−1, 𝑛 controls the mesh partition.

5. Numerical Simulation

To verify the correctness and accuracy of the wavelet Galerkin finite element method, two examples are given to compare the results obtained by this method with an analytical solution. An one-dimensional column of length 𝑙as sketched in Figure 1 is considered. It is assumed that the side walls and the bottom are rigid, frictionless, and impermeable. At top, the stress 𝜎_𝑦and the pressure 𝑝 are prescribed. The boundary conditions are

Figure 1 

Model of fluid saturated porous medium.

𝑢|_𝑦=0=𝜔|_𝑦=0=0,𝜎|_𝑦=𝑙=−𝑃₀𝑓(𝑡),𝑝|_𝑦=𝑙=0. (5.1)

For this model, if the permeability tends to infinity, that is,𝜅→∞, the analytical solutions in time domain are [18]

𝑢_𝑦=𝑃₀𝐸(𝑑₁𝜆₂−𝑑₂𝜆₁)^∞_𝑛=0(−1)^−𝑛{𝑑₂[(𝑡−𝜆₁(𝑙(2𝑛+1)−𝑦))𝐻(𝑡−𝜆₁(𝑙(2𝑛+1)−𝑦))−(𝑡−𝜆₁(𝑙(2𝑛+1)+𝑦))𝐻(𝑡−𝜆₁(𝑙(2𝑛+1)+𝑦))]−𝑑₁[(𝑡−𝜆₂(𝑙(2𝑛+1)−𝑦))𝐻(𝑡−𝜆₂(𝑙(2𝑛+1)−𝑦))−(𝑡−𝜆₂(𝑙(2𝑛+1)+𝑦))𝐻(𝑡−𝜆₂𝑃(𝑙(2𝑛+1)+𝑦))]},𝑝=⁰𝑑₁𝑑₂𝐸(𝑑₁𝜆₂−𝑑₂𝜆₁)^∞_𝑛=0(−1)^−𝑛[𝐻(𝑡−𝜆₁(𝑙(2𝑛+1)−𝑦))+𝐻(𝑡−𝜆₁(𝑙(2𝑛+1)+𝑦))−𝐻(𝑡−𝜆₂(𝑙(2𝑛+1)−𝑦))+𝐻(𝑡−𝜆₂(𝑙(2𝑛+1)+𝑦))], (5.2)(5.3) where𝐸is Young modulus, assuming a Heaviside step function as temporal behavior, that is,𝑓(𝑡)=𝐻(𝑡), and together with vanishing initial conditions:

𝑑_𝑖=𝐸𝜆²_𝑖−(𝜌−𝜌_𝑓)(𝛼−𝑄)𝜆_𝑖𝛽(𝑖=1,2),𝑄=²𝜌_𝑓𝜌_𝛼+𝛽𝜌_𝑓𝜌(𝜅→∞),_𝛼=0.66𝛽𝜌_𝑓. (5.4)

However𝜆_𝑖are the characteristic roots of following characteristic equation

𝐸𝑄𝜌_𝑓𝜆⁴𝛽−(𝐸−²𝑀+(𝜌−𝑄𝜌_𝑓)𝑄𝜌_𝑓+(𝛼−𝑄)²)𝜆²+𝛽²(𝜌−𝑄𝜌_𝑓)𝑀=0. (5.5)

Supposing 𝑄𝐴=𝐸𝜌_𝑓𝛽,𝐵=𝐸−²𝑀+(𝜌−𝑄𝜌_𝑓)𝑄𝜌_𝑓+(𝛼−𝑄)²,𝛽𝐶=²(𝜌−𝑄𝜌_𝑓)𝑀, (5.6) one gets

𝜆₁=−𝜆₃=𝐵+𝐵²−4𝐴𝐶,𝜆2𝐴₂=−𝜆₄=𝐵−𝐵²−4𝐴𝐶.2𝐴 (5.7)

In the first example, the length of column is chosen as𝑙=1000m, and three very different materials, a rock (Berea sandstone), a soil (coarse sand), and a sediment (mud) are chosen. The material data are given in Table 1. In Figures 2, 3, 4, we record the pressure𝑝(𝑡,𝑦=995m), five meters behind the excitation (𝑦=𝑙=1000m). The numerical results (plotted with dot) are compared with the analytical solution (5.3), shown as solid lines in Figures 2, 3, 4. In the second example, the length of column is chosen as𝑙=10m. We choose a material-soil, Figures 5, 6 demonstrate the numerical results—the displacements 𝑢_𝑦(𝑡,𝑦=5m)and the pressure𝑝(𝑡,𝑦=5m). All the figures show that the numerical solutions are perfectly close to the analytical solutions, so the method developed in this paper has a very high degree of calculating accuracy.

Table 1 

The parameters of fluid saturated porous medium.

Figure 2 

The pressure of rock (𝑙=1000m,𝑦=995m).

Figure 3 

The pressure of soil (𝑙=1000m,𝑦=995m).

Figure 4 

The pressure of sediment (𝑙=1000m,𝑦=995m).

Figure 5 

The displacement of soil (𝑙=1000m,𝑦=5m).

Figure 6 

The pressure of soil (𝑙=1000m,𝑦=5m).

6. Conclusion

In this article, the wavelet Galerkin finite element method is constructed by combining the finite element method with wavelet analysis, and is applied to the numerical simulation of the fluid-saturated porous medium elastic wave equation. For the beautiful and deep mathematic properties of Daubechies wavelets, such as the compactly supported property and vanishing moment property, the wavelet Galerkin finite element method has the feature of quick iterative rate and high numerical precision. Moreover, contrasts to ℎ- or 𝑝-based FEM, a new refine algorithm can be presented because of the multi-resolution property of the wavelet analysis. The algorithm can increase the numerical precision by adopting various wavelet basis functions or various wavelet spaces, without refining the mesh.

Acknowledgment

This work was supported by the China Postdoctoral Science Foundation, under Grant no. 20080430930 and by the Natural Science Foundation of Guangdong Province, China, under Grant no. 07300059.