A meshless singular boundary method for elastic wave propagation in 2D partially saturated poroelastic media

Engineering Analysis with Boundary Elements 113 (2020) 82–98

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

A meshless singular boundary method for elastic wave propagation in 2D partially saturated poroelastic media Linlin Sun a,∗, Xing Wei b, Bin Chen c a

Department of Computational Science and Statistics, School of Science, Nantong University, Nantong, Jiangsu, 226019, PR China College of Civil Engineering and Architecture, East China Jiaotong University, Nanchang, Jiangxi, 330013, PR China c College of civil engineering and urban construction, Jiujiang University, Jiujiang, Jiangxi Province, 332005, PR China b

a r t i c l e

i n f o

Keywords: Singular boundary method Meshless Partially saturated media Elastic wave Fundamental solution

a b s t r a c t This paper presents a boundary meshless method, singular boundary method (SBM), in conjunction with the exponential window method (EWM) to simulate the dynamic response of the 2D partially saturated poroelastic media. The problem is firstly solved in the frequency-domain resulted from the time-domain governing equations via the Fourier transformation. Then the SBM approximates the solutions with a linear combination of fundamental solutions with respect to the source points on the physical boundary. To successfully apply the SBM, two issues are addressed. Firstly, the kernel functions, namely the fundamental solutions of frequency-domain governing equations, are derived based on the eigen-analysis. Secondly, the source singularities of the fundamental solutions are desingularized with analytical formulas. After obtaining the SBM formulation, the time-domain solutions can be obtained by the EWM, where the frequency domain windowing technique is introduced to stabilize the oscillations caused by damping parameters. Four numerical examples are carried out to show the validity and accuracy of the proposed method.

1. Introduction Wave propagation in porous media is a significant topic in the geomechanics, seismology, petroleum engineering and civil engineering. Usually, the porous media, such as soil, is supposed as an assemblage of solid skeleton and a pore space. Different from a monophasic continuum in classical elastodynamics, the interaction between the fluid phase in the pore space and solid skeleton is an important portion. Biot [1,2] explored this effect and described the wave propagation in a porous saturated medium, i.e., a solid skeleton fully saturated with a fluid or air, with Biot’s theory. Based on the silhouette drawn in Biot’s theory, the dynamic behavior of the saturated porous media has been widely investigated [3–5]. Nevertheless under natural condition, the pore space likely contains more than one type of fluid. For example, it may contain water and gas or water and oil. This kind of porous material is defined as unsaturated material or partially saturated material. Some assumptions of saturated media are no longer adequate for such media. Therefore, Biot’s theory for saturated porous media fails to predict the mechanical behaviors of unsaturated materials. In the past two decades, several constitutive models are proposed for unsaturated media to take the forces and energies associated with multiphase interactions into account. A comprehensive review on the mathematical models for unsaturated soils can be found in [6]. ∗

To better understand the mechanical property of the unsaturated material, it is of great importance to predict and simulate unsaturated material behavior. Analytical solutions can only be performed for few problems with simple boundary conditions and simple computational domain. Shan, Ling and Ding [7] obtained exact time-domain solutions for one-dimensional transient response of unsaturated single-layer porous media for three types of nonhomogeneous boundary conditions. Li and Schanz [8] derived the semi-analytical solution for the dynamic problem of 1D partially saturated poroelastic column. For general real life problems, numerical methods are essential tools. The finite element method (FEM) is well established in various engineering problems with complex-shaped structures. A number of different finite element formulations [9–11] are proposed for unsaturated media. It has been proved that the FEM can be successfully applied to the wave propagation problems in poroelastic media, especially in bounded domains. However, the wave propagation is mostly observed in semiinfinite or infinite media. A suitable numerical method is required to fulfill the so-called Sommerfeld radiation condition. In the FEM, two ideas are employed to remedy this difficulty. One is the infinite element method [12], and the other is the non-reflecting artificial boundary [13]. The boundary element method (BEM) [14] is another well-known numerical method, which is regarded as an important complementarity to the FEM. The BEM implicitly satisfies the Sommerfeld radiation

Corresponding author. E-mail address: [email protected] (L. Sun).

https://doi.org/10.1016/j.enganabound.2019.12.019 Received 9 August 2019; Received in revised form 22 December 2019; Accepted 23 December 2019 0955-7997/© 2019 Elsevier Ltd. All rights reserved.

L. Sun, X. Wei and B. Chen

Engineering Analysis with Boundary Elements 113 (2020) 82–98

condition [15], which makes the BEM suitable in unbounded domain problems [16]. Another advantage of the BEM is that only boundary meshes are required even for complicated objects. Maghoul, Gatmiri and Duhamel [17] established boundary integral formulations as well as fundamental solutions for the dynamic analysis of unsaturated soils. Gatmiri and Jabbari [18,19] derived the two-dimensional and threedimensional closed form Green’s functions of the governing differential equations for an unsaturated deformable porous medium with linear elastic behavior for a symmetric polar domain in both Laplace-transform and time domains. Li and Schanz [20] established a time domain boundary element formulation and analyzed the dynamic behaviors of 3D partially saturated porous media. Nevertheless, it is time-consuming and tedious for the treatment of singular and near-singular integrals in the BEM when a collocation point is on an integration element or close to the element. The method of fundamental solutions (MFS) [21–24] inherits the boundary-only property of the BEM by using the fundamental solutions as basis functions. To eliminate the troublesome singular and nearsingular integrals, the MFS places the source points on the auxiliary boundary outside the physical domain. The successful application of the MFS to the acoustic wave scattering problem in poroelastic media [25] shows that the MFS is valid for poroelastodynamic problems. However, the MFS suffers the choice of the location of the auxiliary boundary. Much effort [26–28] has been devoted to searching optimal location of the auxiliary boundary. Nevertheless, determining an appropriate auxiliary boundary for the problem in complex-shaped domains is a nontrivial work. The singular boundary method (SBM) [29] can be considered as an improved MFS, which employs the fundamental solutions in terms of source points as kernel functions but places the source points on the physical boundary. To desingularize the source singularity terms, the so-called origin intensity factor (OIF) is introduced in the SBM. The determination of the OIFs affects the accuracy and convergence of the SBM. At first, the SBM calculated the OIFs with an inverse interpolation technique [30], which is cumbersome and time-consuming. Later some analytical and empirical formulas were derived to quickly evaluate the OIFs for various problems, such as potential problems [31–33], heat conduction [34–36], elastic wave propagation [37,38], wave scattering [39–42], water wave [43], dynamic poroelastic problems [44], plate vibration [45], photonic crystals [46,47], and electromagnetic scattering problems [48], etc. In light of the fundamental solutions and OIFs, the SBM inherits the advantages of the MFS, while overcoming the auxiliary boundary issue of the MFS. Besides the above mentioned merits of the SBM, the shortcomings of the SBM also need to be noticed, such as the issues of inhomogeneous partial differential equations, problems without known fundamental solutions, dense coefficient matrices, etc. For inhomogeneous problems, other numerical methods are introduced to obtain the particular solutions, such as the multiple reciprocity method (MRM) [33] or dual reciprocity method (DRM) [49]. For problems without known fundamental solutions, including non-linear problems, the SBM can cooperate with the analog equation method [50]. For large scale-problems, the storage and computational cost of the SBM will be prohibitive because of dense coefficient matrices. To alleviate the drawback, the fast algorithms such as fast multipole method (FMM) [39], precorrected fast Fourier transform (pFFT) [51] or adaptive cross approximation (ACA) [52] can be used to accelerate the SBM. In this study, the SBM is adopted to simulate the dynamic interaction of the solid, air and fluid phases within the partially saturated porous media. Because it is not trivial to obtain the time-domain fundamental solutions, the problem is transformed to the frequency-domain with the Fourier transformation. The frequency-domain fundamental solutions are derived firstly based on the eigen-analysis. Then the SBM is formulated to deal with the frequency-domain problem, where the source singularities of fundamental solutions are studied and desingularized with OIFs. Finally, the exponential window method (EWM) [53] is introduced to perform the inverse discrete Fourier transformation and obtain

the time-domain results. Appendix A gathers the symbols and notations used throughout the paper. 2. Governing equations Consider a triphasic partially saturated poroelastic media consisting of the solid skeleton (solid phase), pore water (wetting fluid phase) and pore air (non-wetting fluid phase). To establish the governing equations for dynamic partially saturated poroelasticity, it is essential to make the following assumptions [20]: 1) There is no phase transformation during the fluids penetrate into the solid skeleton. The evaporation of the wetting fluid and dissolving of non-wetting fluid are ignored; 2) All phases have the same temperature irrespective of the temperature variation. In the following, the mathematical model suggested by Lewis and Schrefler [8,20,54] will be presented to describe the dynamic behavior of partially saturated poroelastic media. 2.1. Basic equations [

Let 𝒖𝑠 = (𝑢𝑠1 , 𝑢𝑠2 )𝑇 and 𝝈 𝑠 =

𝑠 ] 𝜎12 respectively represent the dis𝑠 𝜎22

𝑠 𝜎11 𝑠 𝜎21

placement vector and stress tensor of the solid skeleton. 𝒖𝑓 = (𝑢𝑓1 , 𝑢𝑓2 )𝑇 , pf (f = w,a) respectively denote the relative displacement and pressure of fluids, where the superscripts w and a stand for the pore water and pore air, respectively. The total stresses 𝜎 ij are given by 𝜎𝑖𝑗 = 𝜎𝑖𝑗𝑠 − 𝛿𝑖𝑗 𝛼(𝑆𝑤 𝑝𝑤 + 𝑆𝑎 𝑝𝑎 )

( ) ( ) 𝜕𝑢𝑠𝑗 𝜕𝑢𝑠𝑖 2 𝑠 = 𝐾 − 𝐺 𝛿𝑖𝑗 (∇ ⋅ 𝒖 ) + 𝐺 + 3 𝜕 𝑥𝑗 𝜕 𝑥𝑖 − 𝛿𝑖𝑗 𝛼(𝑆𝑤 𝑝𝑤 + 𝑆𝑎 𝑝𝑎 ), (𝑖, 𝑗 = 1, 2),

(1)

where x1 ,x2 stand for the horizontal and vertical coordinates and ∇ = ( 𝜕 𝑥𝜕 , 𝜕 𝑥𝜕 )𝑇 denotes the gradient operator. Sw and Sa are the degree of 1

2

saturation of water and air, respectively. Sw +Sa = 1. K and G represent drained bulk modulus of the mixture and the shear modulus. 𝛿 ij is the Kronecker delta. 𝛼 = 1 − 𝐾∕𝐾𝑠 with Ks the bulk modulus of the solid skeleton. The balance equation of momentum for each phase can be written as [8,20] solid ∶ ∇ ⋅ 𝝈 + 𝑭 𝑠 = ρ

wet t ing Fluid ∶ 𝑛𝑆𝑤

𝜕 2 𝒖𝑠 𝜕 𝑡2

𝜕 𝒖𝑤 𝜕𝑡

+ 𝑛𝑆𝑤 ρ𝑤

𝜕 2 𝒖𝑤 𝜕 𝑡2

+ 𝑛𝑆𝑎 ρ𝑎

( = −κ𝑤 ∇𝑝𝑤 + ρ𝑤

non − wet t ing Fluid ∶ 𝑛𝑆𝑎

𝜕 𝒖𝑎 𝜕𝑡

( = − κ𝑎

𝜕 2 𝒖𝑎 𝜕 𝑡2

,

𝜕 2 𝒖𝑤 𝜕 2 𝒖𝑠 + ρ𝑤 2 𝜕𝑡 𝜕 𝑡2

(2) ) ,

𝜕 2 𝒖𝑎 𝜕 2 𝒖𝑠 ∇𝑝 + ρ𝑎 + ρ𝑎 𝜕 𝑡2 𝜕 𝑡2 𝑎

(3) ) ,

(4)

where Fs denotes bulk body force, 𝜌s , 𝜌w and 𝜌a stand for the density of solid skeleton, pore water and pore air, respectively. n is the porosity. 𝜌 = (1 − n)𝜌s + nSw 𝜌w + nSa 𝜌a denotes the averaged density of the mix𝐾 𝑘̄ ture. κ𝑓 = 𝑟𝑓 (𝑓 = 𝑤,𝑎) is the phase permeability of fluid where 𝑘̄ de𝜂𝑓

notes the intrinsic fluid permeability. 𝜂 f represents fluid viscosity. Krf stands for relative fluid phase permeability and has the following form for the water and air, respectively: (2+3𝜗)∕𝜗

𝐾𝑟𝑤 = 𝑆𝑒

,

] ( )2 [ (2+𝜗)∕𝜗 , 𝐾𝑟𝑎 = 1 − 𝑆𝑒 1 − 𝑆𝑒 83

(5)

(6)

L. Sun, X. Wei and B. Chen

Engineering Analysis with Boundary Elements 113 (2020) 82–98

where ϑ is the pore size distribution index, Se represents the effective wetting fluid saturation as ⎧0 , ⎪ ⎪ 𝑆𝑤 − 𝑆𝑟𝑤 , 𝑆𝑒 = ⎨ ⎪ 𝑆𝑟𝑎 − 𝑆𝑟𝑤 ⎪1 , ⎩

After applying the Fourier transformation to Eqs. (3) and (4), the relative fluid displacements are recast as ) ( (16) 𝒖̃ 𝑤 = 𝛽 ∇𝑝̃𝑤 − ρ𝑤 𝜔2 𝒖̃ 𝑠 ,

𝑆𝑤 ≤ 𝑆𝑟𝑤 𝑆𝑟𝑤 ≤ 𝑆𝑤 ≤ 𝑆𝑟𝑎 ,

( ) 𝒖̃ 𝑎 = 𝛾 ∇𝑝̃𝑎 − ρ𝑎 𝜔2 𝒖̃ 𝑠 ,

(7)

𝑆𝑤 ≥ 𝑆𝑟𝑎

where 𝛽 =

𝑠

−1∕𝜗

𝑆𝑢 =

are

capillary

𝜗+1 𝜗(𝑆 −𝑆 ) 𝑆 −𝑆 − 𝑟𝑎𝑝𝑑 𝑟𝑤 ( 𝑆 𝑤 −𝑆𝑟𝑤 ) 𝜗 , 𝑟𝑎 𝑟𝑤

pressure

and

𝑝𝑐

𝐴8 ∇ ⋅ 𝒖̃ 𝑠 + 𝐴9 𝑝̃𝑤 + 𝐴10 ∇2 𝑝̃𝑎 + 𝐴11 𝑝̃𝑎 = 𝐼̃𝑎 ,

(20)

3.1. Derivation of fundamental solutions In the SBM, the fundamental solutions play a pivotal role. Up to now, there exist many reports on the derivation of the fundamental solutions of dynamic partially saturated poroelastic problems [17–20], which mainly concern with Laplace-domain and time-domain fundamental solutions. The Laplace-domain fundamental solutions are usually derived by using the method of Hörmander. The main idea of this method is to use the determinant of the adjoint operator to obtain a higher order Helmholtz equation, whose fundamental solutions can be easily derived, then apply the co-factors of the adjoint operator to the fundamental solutions of the higher order Helmholtz equation and result in the fundamental solutions for partially saturated poroelasticity. For the case of time domain problem, the inverse Laplace transformation technique can be employed to transform the Laplace-domain fundamental solutions into time-domain. In this part, we focus on the derivation of the frequency-domain fundamental solutions with the wave decomposition method. Different from the method of Hörmander, this method decouples the governing equations with singular source terms into four scalar equations which can be expressed in the matrix form. With the help of the eigenvalues and eigenvectors of the coefficient matrix, the solutions of the scalar equations can be determined, and thus the fundamental solutions can be obtained. In the following, the detailed derivations of fundamental solutions are provided. The fundamental solutions 𝑢̄ 𝑠𝑖𝑘 (𝒙, 𝒚 ), 𝑡̄𝑠𝑖𝑘 (𝒙, 𝒚 ), 𝑝̄𝑤 (𝒙, 𝒚 ) and 𝑝̄𝑎 (𝒙, 𝒚 ) are the solutions of Eqs. (18) to (20) at point x with a singular load at point y. Three different kinds of singular loads including the source terms of the solid phase equilibrium equation at both x1 and x2 directions and the source terms of the equations for wetting and nonwetting fluid phases will be investigated.

(10)

𝜕 𝑆𝑤 𝜕 𝑝𝑐

(19)

3. Fundamental solutions

where Iw , Ia represent wetting and non-wetting fluid source terms, 𝜕𝑆 𝜕𝑆 𝑛 𝜁 = 𝛼− , 𝑆𝑤𝑤 = 𝑆𝑤 + 𝑝𝑐 𝜕 𝑝𝑤𝑐 , 𝑆𝑎𝑎 = 𝑆𝑎 − 𝑝𝑐 𝜕 𝑝𝑤𝑐 , 𝑝𝑐 = 𝑝𝑎 − 𝑝𝑤 = 𝐾 𝑝𝑑 𝑆𝑒

𝐴4 ∇ ⋅ 𝒖̃ 𝑠 + 𝐴5 ∇2 𝑝̃𝑤 + 𝐴6 𝑝̃𝑤 + 𝐴7 𝑝̃𝑎 = 𝐼̃𝑤 ,

The detailed forms of Ai (i = 1, 2, …, 11) in above formulas are provided in Appendix B.

𝜕∇ ⋅ 𝒖𝑠

𝜕 𝑝𝑤 + (𝜁 𝑆𝑤𝑤 𝑆𝑎 + 𝑆𝑢 ) 𝜕𝑡 𝜕𝑡 ( ) 𝑎 𝜕∇ ⋅ 𝒖𝑎 𝜕𝑝 𝑛 + 𝜁 𝑆𝑎𝑎 𝑆𝑎 + 𝑆𝑎 − 𝑆𝑢 + 𝑛𝑆𝑎 = 𝐼 𝑎, 𝐾𝑎 𝜕𝑡 𝜕𝑡

(17) √ j = −1 is the imaginary

(8)

To ensure the existence of solutions, the continuity equations for fluid phases [20,54] also should be considered: ( ) 𝑤 𝜕∇ ⋅ 𝒖𝑠 𝜕𝑝 𝑛 𝛼𝑆𝑤 + 𝜁 𝑆𝑤𝑤 𝑆𝑤 + 𝑆𝑤 − 𝑆𝑢 𝜕𝑡 𝐾𝑤 𝜕𝑡 𝜕∇ ⋅ 𝒖𝑤 𝜕 𝑝𝑎 + (𝜁 𝑆𝑎𝑎 𝑆𝑤 + 𝑆𝑢 ) + 𝑛𝑆𝑤 = 𝐼 𝑤, (9) 𝜕𝑡 𝜕𝑡

𝛼𝑆𝑎

𝛾=

κ − j𝜔𝑛𝑆 −𝑎κ ρ 𝜔2 , 𝑎 𝑎 𝑎

unit, 𝜔 denotes angular frequency. Then applying Fourier transformation to Eqs. (8), (9) and (10) and substituting Eqs. (16) and (17) into the resulted equations, we obtain: ( ) 𝐺 𝑠 𝐺∇2 𝒖̃ 𝑠 + 𝐾 + (18) ∇∇ ⋅ 𝒖̃ 𝑠 + 𝜔2 𝐴1 𝒖̃ 𝑠 − 𝐴2 ∇𝑝̃𝑤 − 𝐴3 ∇𝑝̃𝑎 = −𝑭̃ , 3

with Srw the residual wetting fluid saturation and Sra the non-wetting fluid entry saturation. Inserting Eq. (1) into Eq. (2) leads to the equation of motion for solid phase as ( ) 1 𝐺∇2 𝒖𝑠 + 𝐾 + 𝐺 ∇∇ ⋅ 𝒖𝑠 − 𝛼(𝑆𝑤 ∇𝑝𝑤 + 𝑆𝑎 ∇𝑝𝑎 ) + 𝑭 𝑠 3 𝜕 2 𝒖𝑤 𝜕 2 𝒖𝑠 𝜕 2 𝒖𝑎 =ρ + 𝑛𝑆𝑤 ρ𝑤 + 𝑛𝑆𝑎 ρ𝑎 , 2 2 𝜕𝑡 𝜕𝑡 𝜕 𝑡2

κ − j𝜔𝑛𝑆 −𝑤κ ρ 𝜔2 , 𝑤 𝑤 𝑤

= −𝜗(𝑆𝑤 − 𝑆𝑟𝑤 ),

respectively.

The boundary conditions are given by 𝑢𝑠𝑖 = 𝑢̂ 𝑠𝑖 , on Γ𝑠𝑢 ,

(11)

𝑡𝑠𝑖 = 𝜎𝑖𝑠1 𝑛1 +𝜎𝑖𝑠2 𝑛2 = 𝑡̂𝑠𝑖 , on Γ𝑠𝑡 ,

(12)

for the solid phase, and 𝑢𝑓𝑖 = 𝑢̂ 𝑓𝑖 , on Γ𝑓𝑢 ,

(13)

𝑝𝑓 = 𝑝̂𝑓 , on Γ𝑓𝑝 ,

(14)

for the fluid phase (f = w, a) with the initial conditions within domain defined as 𝜕 𝒖𝑠 || | | 𝒖𝑠 |𝑡=0 = = 0, 𝒖𝑓 | = 0, and 𝑝𝑓 | = 0 (𝑓 = 𝑤, 𝑎), (15) |𝑡=0 |𝑡=0 𝜕𝑡 ||𝑡=0 where n = (n1 ,n2 ) is the outward normal vector, 𝑢̂ 𝑠𝑖 , 𝑡̂𝑠𝑖 , 𝑢̂ 𝑓𝑖 and 𝑝̂𝑓 are the prescribed solid displacements, tractions, relative fluid displacements and pressure on the boundary, respectively. Based on abovementioned equations, the dynamic response of the partially saturated poroelasticity can be simulated. For more details of parameters, the readers are referred to [8].

3.1.1. Solid loads First, we consider the solid phase. The source terms for this case are 𝑠 𝑭̃ = 𝛿(𝒙 − 𝒚 )𝒆𝑘 (k=1,2), 𝐼̃𝑤 = 𝐼̃𝑎 = 0, where e1 and e2 are the unit vectors along x1 and x2 direction. According to reference [55], the solutions of the equations can be written as ( ( ) ) (21) 𝒖̃ 𝑠 = ∇∇ ⋅ 𝐴𝐿 𝒆𝑘 − ∇ × ∇ × 𝐴𝑇 𝒆𝑘 ,

2.2. Governing equation in frequency-domain It is not a trivial work to directly obtain the time-domain solution for the considered problem. A more general idea is to transform the timedomain governing equation into frequency-domain, and then obtain its frequency-domain solutions with numerical methods, which could be inversely transformed to time-domain solutions.

𝑝̃𝑤 = ∇ ⋅ (𝐴𝑤 𝒆𝑘 ), 𝑝̃𝑎 = ∇ ⋅ (𝐴𝑎 𝒆𝑘 ), where AL , AT , Aw and Aa are the variables to be determined. 84

(22)

L. Sun, X. Wei and B. Chen

Engineering Analysis with Boundary Elements 113 (2020) 82–98

Let 𝜏 be the fundamental solution of the Laplace operator, we have ( ) ( ) ( ) ∇2 𝜏𝒆𝑘 = ∇∇ ⋅ 𝜏𝒆𝑘 − ∇ × ∇ × 𝜏𝒆𝑘 = −𝛿(𝒙 − 𝒚 )𝒆𝑘 ,

Combing Eqs. (35) and (36) with Eqs. (21), (22) and (32), the fundamental solutions can be derived as

(23)

Taking Eqs. (21) and (22) into Eqs. (18), (19) and (20), we get ( ) 4 (24) 𝐾 + 𝐺 ∇2 𝐴𝐿 + 𝜔2 𝐴1 𝐴𝐿 − 𝐴2 𝐴𝑤 − 𝐴3 𝐴𝑎 = 𝜏, 3

𝑢̄ 𝑠𝑖𝑘 = 𝐴𝛿𝑖𝑘 − 𝐵 𝑟,𝑖 𝑟,𝑘 , (𝑖, 𝑘 = 1, 2),

(37)

𝑝̄𝑤 𝑘 = 𝐶 𝑟,𝑘 , (𝑘 = 1),

(38)

(39)

𝐺∇2 𝐴𝑇 + 𝜔2 𝐴1 𝐴𝑇 = 𝜏,

(25)

𝑝̄𝑎𝑘 = 𝐷𝑟,𝑘 , (𝑘 = 1, 2),

𝐴4 ∇2 𝐴𝐿 + 𝐴5 ∇2 𝐴𝑤 + 𝐴6 𝐴𝑤 + 𝐴7 𝐴𝑎 = 0,

(26)

where

(27)

𝐴=

𝐴8 ∇ 𝐴𝐿 + 𝐴9 𝐴𝑤 + 𝐴10 ∇ 𝐴𝑎 + 𝐴11 𝐴𝑎 = 0, 2

2

[ ( )] ∑ 𝐾 (𝑧 ) 𝐾 (𝑧 ) 1 − 𝑔𝑚 1 𝑚 + 𝑔4 𝐾0 (𝑧4 ) + 1 4 , 2𝜋 𝑧𝑚 𝑧4 𝑚=1,2,3 [ ] ∑ 1 𝐵 = − 𝑔𝑚 𝐾2 (𝑧𝑚 ) + 𝑔4 𝐾2 (𝑧4 ) , 2𝜋 𝑚=1,2,3

Reformulate Eqs. (24), (26) and (27) into the matrix form ⎛ 𝐴𝐿 ⎞ ⎛𝐴𝐿 ⎞ ⎛𝜏 ⎞ 𝑀1 ∇2 ⎜𝐴𝑤 ⎟ + 𝑀2 ⎜𝐴𝑤 ⎟ = ⎜0⎟, ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎝ 𝐴𝑎 ⎠ ⎝ 𝐴𝑎 ⎠ ⎝ 0 ⎠

(28)

𝐶 =

j ∑ 𝑟′ 2𝑚 𝑔𝑚 ( ) 𝐾1 𝑧𝑚 , 2𝜋 𝑚=1,2,3 𝑘𝑚

𝐷=

j ∑ 𝑟′ 3𝑚 𝑔𝑚 𝐾1 (𝑧𝑚 ). 2𝜋 𝑚=1,2,3 𝑘𝑚

where ⎡𝐾 + 4∕3𝐺 𝑀1 = ⎢ 𝐴4 ⎢ 𝐴8 ⎣ Multiplying M1

−1

0 ⎤ ⎡𝜔 2 𝐴 1 0 ⎥, 𝑀2 = ⎢ 0 ⎥ ⎢ 𝐴10 ⎦ ⎣ 0

0 𝐴5 0

−𝐴2 𝐴6 𝐴9

−𝐴3 ⎤ 𝐴 7 ⎥. ⎥ 𝐴11 ⎦

(29)

in which zm = jkm r, 𝑟 =

on both sides and we get:

⎛ 𝐴𝐿 ⎞ ⎛ 𝐴𝐿 ⎞ ⎛1⎞ ∇2 ⎜𝐴𝑤 ⎟ + 𝑀 ⎜𝐴𝑤 ⎟ = 𝜏𝑀1−1 ⎜0⎟ = 𝜏𝒈, ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎝ 𝐴𝑎 ⎠ ⎝ 𝐴𝑎 ⎠ ⎝0⎠

(30)

(31)

( [( ) ) ) 𝐵 𝐵 ( 𝑡̄𝑠𝑖𝑘 = 𝜆 𝐴′ − 𝐵 ′ − 𝑟,𝑘 𝑛𝑖 + 𝜇 𝐴′ − 𝑟,𝑛 𝛿𝑖𝑘 + 𝑟,𝑖 𝑛𝑘 𝑟 𝑟 ( ) ] 2𝐵 2𝐵 𝑟,𝑘 𝑛𝑖 + 2 −𝐵 ′ + 𝑟,𝑖 𝑟,𝑘 𝑟,𝑛 , 𝑖, 𝑘 = 1, 2, − 𝑟 𝑟

In light of the eigen-analysis of matrix M, the solutions of Eq. (30) can be expressed as ⎛𝐴 𝐿 ⎞ ⎛ 1 ⎞ ⎛ 1 ⎞ ⎛ 1 ⎞ ⎜𝐴 ⎟ = 𝜙 ⎜𝑟′ ⎟ + 𝜙 ⎜𝑟′ ⎟ + 𝜙 ⎜𝑟′ ⎟, 1 2 3 𝑤 21 22 ⎜ ⎟ ⎜ ′ ⎟ ⎜ ′ ⎟ ⎜ ′ 23 ⎟ ⎝𝐴 𝑎 ⎠ ⎝𝑟 31 ⎠ ⎝𝑟 32 ⎠ ⎝𝑟 33 ⎠

𝑟′ 3𝑖 =

𝑚21 𝑚33 − 𝑚23 𝑚31 − 𝑘2𝑖 𝑚21 𝑘2𝑖 (𝑚22

+ 𝑚33 ) −

𝑘4𝑖

− 𝑚22 𝑚33 + 𝑚23 𝑚32

𝑚22 𝑚31 − 𝑚21 𝑚32 − 𝑘2𝑖 𝑚31 𝑘2𝑖 (𝑚22

+ 𝑚33 ) − 𝑘4𝑖 − 𝑚22 𝑚33 + 𝑚23 𝑚32

(40)

where r, n = r, 1 n1 + r, 2 n2 , 𝜆 = 𝐾 − 23 𝐺, 𝜇 = G, and the sign ()′ in A and B representing the derivative with respect to r. For the fluid, the derivatives of pressure in the direction of the outward normal vector are ( ) 𝜕 𝑝̄𝑤 𝐶 𝐶 𝑞̄𝑘𝑤 = 𝑘 = 𝑛𝑘 + 𝐶 ′ − (41) 𝑟 𝑟 , 𝑘 = 1, 2, 𝜕𝒏 𝑟 𝑟 ,𝑘 ,𝑛

(32)

where (1, r′2i ,r′3i )T (i=1,2,3) is the eigenvector of M and r′2i , r′3i (i=1,2,3) are given by 𝑟′ 2𝑖 =

𝑥 −𝑦

the collocation point x = (x1 ,x2 ) and source point y = (y1 ,y2 ). 𝑟,𝑖 = 𝑖 𝑟 𝑖 , Ki is the modified Bessel function of ith order. Then the traction tensor 𝑡̄𝑠𝑖𝑘 arising from the fundamental solutions can be derived by differentiating the displacement tensor and substituting in the constitutive relationship as follows:

with 𝑀 = 𝑀1−1 𝑀2 and 1 ⎛ ⎞ 1 ⎜ −𝐴 ∕𝐴 ⎟. 𝒈= 4 5 ⎟ 𝐾 + 4∕3𝐺 ⎜ ⎝−𝐴8 ∕𝐴10 ⎠

√ (𝑥1 − 𝑦1 )2 + (𝑥2 − 𝑦2 )2 is the distance between

𝜕 𝑝̄𝑎𝑘

( ) 𝐷 𝐷 𝑟 𝑟 , 𝑘 = 1, 2. 𝑛 + 𝐷′ − 𝑟 𝑘 𝑟 ,𝑘 ,𝑛

,

𝑞̄𝑘𝑎 =

,

3.1.2. Fluid load Two types of fluids, namely wetting and non-wetting fluid phases, are considered in this part. Firstly, the source term of the continuity equation of wetting fluid 𝑠 phase is treated as singular load, that is 𝑭̃ =0, 𝐼̃𝑤 = −𝛿(𝒙 − 𝒚 ),𝐼̃𝑎 = 0. Referring to reference [55], the solution of the solid displacement for this case is expressed as:

where 𝑘2𝑖 (i = 1,2,3) is the eigenvalue of M, and mik denotes the entry in the ith row and the kth column of M. The details of 𝑘𝑖 are provided in Appendix C. Taking Eq. (32) into Eq. (30), we obtain ∇2 𝜙𝑖 + 𝑘2𝑖 𝜙𝑖 = 𝑔𝑖 𝜏,

(33)

𝜕𝒏

=

𝒖̃ 𝑠 = ∇𝜑 − ∇ × 𝝍 .

where gi (i = 1, 2, 3) satisfy ⎛ 1 ⎞ ⎛ 1 ⎞ ⎛ 1 ⎞ 𝒈 = 𝑔1 ⎜𝑟′ 21 ⎟ + 𝑔2 ⎜𝑟′ 22 ⎟ + 𝑔3 ⎜𝑟′ 23 ⎟. ⎜ ′ ⎟ ⎜ ′ ⎟ ⎜ ′ ⎟ ⎝𝑟 31 ⎠ ⎝𝑟 32 ⎠ ⎝𝑟 33 ⎠

(34)

The solution of Eq. (33) is ) 𝑔 ( 𝜙𝑖 (𝑟) = − 𝑖 ln 𝑟 + 𝐾0 (j𝑘𝑖 𝑟) , 2 2𝜋𝑘𝑖

(35)

(42)

(43)

Inserting Eqs. (22) and (43) into governing Eqs. (18), (19) and (20), we can obtain the following equations: ( ) 4 (44) 𝐾 + 𝐺 Δ𝜑 + 𝜔2 𝐴1 𝜑 − 𝐴2 𝑝𝑤 − 𝐴3 𝑝𝑎 = 0, 3 𝐺 ∇2 𝝍 + 𝜔2 𝐴1 𝝍 = 0,

(45)

𝐴4 ∇2 𝜑 + 𝐴5 ∇2 𝑝𝑤 + 𝐴6 𝑝𝑤 + 𝐴7 𝑝𝑎 = −𝛿(𝒙 − 𝒚 ),

(46)

𝐴8 ∇2 𝜑 + 𝐴9 𝑝𝑤 + 𝐴10 ∇2 𝑝𝑎 + 𝐴11 𝑝𝑎 = 0.

(47)

where K0 is the modified Bessel function of zero order. 𝜔2 𝐴

Let 𝑘24 = 𝐺 1 , 𝑔4 = 𝐺1 . Then the solution of 2D shear wave associated with solid phase Eq. (25) is ) 𝑔 ( 𝐴𝑇 = − 4 ln 𝑟 + 𝐾0 (j𝑘4 𝑟) . (36) 2𝜋𝑘24 85

L. Sun, X. Wei and B. Chen

Engineering Analysis with Boundary Elements 113 (2020) 82–98

It can be observed that Eq. (45) is uncoupled with other equations. Thus we assume 𝝍 = 0. Then we recast Eqs. (44), (46) and (47) in matrix form as ⎛𝜑⎞ ⎛𝜑⎞ ∇2 ⎜𝑝𝑤 ⎟ + 𝑀 ⎜𝑝𝑤 ⎟ = −𝛿𝒉𝑤 , ⎜ 𝑎⎟ ⎜ 𝑎⎟ ⎝𝑝 ⎠ ⎝𝑝 ⎠

3.2. Singular behaviors of the fundamental solutions Whent a collocation point approaches a source point, the fundamental solutions become singular. Before isolating the singularities of fundamental solutions and obtaining the OIFs in the SBM, it is essential to investigate the singular behaviors of the fundamental solutions. Based on the asymptotic expansion of the modified Bessel function, when a collocation point x approaches a source point y, we obtain the following equations: [ ] 1 𝑢̄ 𝑠𝑖𝑘 = 𝜒1 − ln 𝑟 𝛿𝑖𝑘 − 𝜒2 𝛿𝑖𝑘 + 𝜒3 𝑟,𝑖 𝑟,𝑘 , 𝑢̄ 𝑠𝑖3 = 0, 𝑢̄ 𝑠𝑖4 = 0, (𝑖, 𝑘 = 1, 2), (64) 2𝜋

(48)

where M is the same as the one in Eq. (30), and 𝒉𝑤 = (0, 1∕𝐴5 , 0)𝑇 . Similar to Eq. (30), we construct the solution as ⎛𝜑⎞ ⎛ 1 ⎞ ⎛ 1 ⎞ ⎛ 1 ⎞ ⎜𝑝𝑤 ⎟ = 𝜙 ⎜𝑟′ ⎟ + 𝜙 ⎜𝑟′ ⎟ + 𝜙 ⎜𝑟′ ⎟, 1 2 3 ⎜ 𝑎⎟ ⎜ ′ 21 ⎟ ⎜ ′ 22 ⎟ ⎜ ′ 23 ⎟ ⎝𝑝 ⎠ ⎝𝑟 31 ⎠ ⎝𝑟 32 ⎠ ⎝𝑟 33 ⎠

(49)

Taking Eq. (49) into Eq. (48) leads to −ℎ𝑤 𝑖 𝛿(𝒙 −

∇ 𝜙𝑖 +

𝑘2𝑖 𝜙𝑖

=

where

ℎ𝑤 𝑖 (𝑖

= 1, 2, 3) satisfies

2

𝑤

𝒉 =

⎛ 1 ⎞

⎜𝑟′ ⎟ + ℎ𝑤 1 ⎜ 21 ⎟ ⎝𝑟′ 31 ⎠

𝑡̄𝑠𝑖𝑘 =

𝒚 ),

⎛ 1 ⎞

⎜𝑟′ ⎟ + ℎ𝑤 2 ⎜ 22 ⎟ ⎝𝑟′ 32 ⎠

(50)

⎛ 1 ⎞

⎜𝑟′ ⎟. ℎ𝑤 3 ⎜ 23 ⎟ ⎝𝑟′ 33 ⎠

(51) ℎ𝑤 𝑖

The solutions of Eq. (50) is 𝜙𝑖 (𝑟) = 2𝜋 𝐾0 (j𝑘𝑖 𝑟). Combing the solutions with Eqs. (22), (43) and (49), we obtain the following fundamental solutions: 𝑢̄ 𝑠𝑖3 = 𝐴𝑤 𝑟,𝑖 , (𝑖 = 1, 2), 𝑝̄𝑤 = 3 𝑝̄𝑎3 =

1 2𝜋

∑ 𝑚=1,2,3

(53)

𝑟′ 3𝑚 ℎ𝑤 𝑚 𝐾0 (𝑧𝑚 ),

𝑞̄3𝑎 = −

j𝑟,𝑛 ∑

(54)

2𝜋

(56)

𝑟′ 3𝑚 ℎ𝑤 𝑚 𝑘𝑚 𝐾1 (𝑧𝑚 ),

(57)

𝑚=1,2,3

j𝑟,𝑛 ∑ 2𝜋

𝑟′ 2𝑚 ℎ𝑤 𝑚 𝑘𝑚 𝐾1 (𝑧𝑚 ),

𝑚=1,2,3

[ ( ) ] [ ( ) ] 1 1 𝑡̄𝑠𝑖3 = 𝜒4 − ln 𝑟 + 𝜒5 𝑛𝑖 , 𝑡̄𝑠𝑖4 = 𝜒6 − ln 𝑟 + 𝜒7 𝑛𝑖 , (𝑖 = 1, 2), 2𝜋 2𝜋

(66)

[ ] [ ] 1 1 𝑤 𝑤 𝑝̄𝑤 𝑘 = 0, 𝑝̄3 = 𝜒8 − 2𝜋 ln 𝑟 + 𝜒9 , 𝑝̄4 = 𝜒10 − 2𝜋 ln 𝑟 + 𝜒11 , (𝑘 = 1, 2),

(67)

[ ] [ ] 1 1 𝑝̄𝑎𝑘 = 0, 𝑝̄𝑎3 = 𝜒14 − ln 𝑟 + 𝜒15 , 𝑝̄𝑎4 = 𝜒16 − ln 𝑟 + 𝜒17 , (𝑘 = 1, 2), 2𝜋 2𝜋 (69)

−j ∑ where 𝐴𝑤 = 2𝜋 𝑚=1,2,3 𝑘𝑚 ℎ𝑤 𝑚 𝐾1 (𝑧𝑚 ). Then the traction tensor and outward normal derivatives of pressure are ( [ ] ) 𝐴𝑤 𝐴𝑤 𝑡̄𝑠𝑖3 = (𝜆 + 2𝜇) + 𝜆(𝐴𝑤 )′ 𝑛𝑖 + 2𝜇 − + (𝐴𝑤 )′ 𝑟𝑖 𝑟𝑛 , (𝑖 = 1, 2), (55) 𝑟 𝑟

𝑞̄3𝑤 = −

( ) ( ) [ ( ) ] 1 1 𝑟,𝑛 1 𝑟,𝑛 𝑞̄𝑘𝑎 = 𝜒18 − , 𝑞̄4𝑎 = 𝜒16 − , ln 𝑟 + 𝜒19 𝑛𝑘 , 𝑞̄3𝑎 = 𝜒14 − 2𝜋 2𝜋 𝑟 2𝜋 𝑟 (70) where 𝑔 + 𝑔2 + 𝑔3 + 𝑔4 𝜒1 = 1 , 2

𝜒2 =

For the case of singular load for continuity equation of non-wetting fluid phase, the fundamental solutions can be obtained in a similar way as the case of the wetting fluid phase. The fundamental solutions are 𝑢̄ 𝑠𝑖4 = 𝐴𝑎 𝑟,𝑖 , (𝑖 = 1, 2), 𝑝̄𝑤 = 4 𝑝̄𝑎4

(58)

1 ∑ ′ 𝑟 ℎ𝑎 𝐾 (𝑧 ), 2𝜋 𝑚=1,2,3 2𝑚 𝑚 0 𝑚

𝑞̄4𝑎

=−

j𝑟,𝑛 ∑ 2𝜋

𝑚=1,2,3

j𝑟,𝑛 ∑ 2𝜋

𝑟

𝑚=1,2,3

where 𝐴𝑎 =

−j 2𝜋

∑

(

) 𝑔1 + 𝑔2 + 𝑔3 − 𝑔4 , 2 ∑

𝑚=1,2,3

) 𝑘2𝑚 ℎ𝑤 𝑚

,

[ ( ] ( ( ))) ∑ 𝑘2𝑚 ℎ𝑤 ∑ 𝑘2𝑚 ℎ𝑤 j𝑘𝑚 1 𝑚 𝑚 𝜒5 = 2(𝜆 + 𝜇) 𝜏̄ + ln −𝜇 , 2𝜋 2 2 2 𝑚=1,2,3 𝑚=1,2,3

(61)

( 𝜒6 = −(𝜆 + 𝜇)

(62)

𝑎 3𝑚 ℎ𝑚 𝑘𝑚 𝐾1 (𝑧𝑚 ),

𝑎 𝑚=1,2,3 𝑘𝑚 ℎ𝑚 𝐾1 (𝑧𝑚 ),

1 2𝜋

(

(60)

𝑟′ 2𝑚 ℎ𝑎𝑚 𝑘𝑚 𝐾1 (𝑧𝑚 ), ′

𝜒3 = −

𝜒4 = −(𝜆 + 𝜇)

( [ ] ) 𝐴𝑎 𝐴𝑎 𝑡̄𝑠𝑖4 = (𝜆 + 2𝜇) + 𝜆(𝐴𝑎 )′ 𝑛𝑖 + 2𝜇 − + (𝐴𝑎 )′ 𝑟,𝑖 𝑟,𝑛 , (𝑖 = 1, 2), 𝑟 𝑟 𝑞̄4𝑤 = −

( ) ( ) ( ) [ 𝑔 j𝑘3 𝑔 j𝑘1 𝑔 j𝑘2 1 𝑔1 + 𝑔2 + 𝑔3 + 𝑔4 𝜏̄ + 1 ln + 2 ln + 3 ln 2𝜋 2 2 2 2 2 2 2 ( ) ] 𝑔 + 𝑔2 + 𝑔3 − 𝑔4 𝑔 j𝑘4 + 4 ln − 1 , 2 2 4

(59)

1 ∑ ′ = 𝑟 ℎ𝑎 𝐾 (𝑧 ), 2𝜋 𝑚=1,2,3 3𝑚 𝑚 0 𝑚

(65)

( ) ( ) [ ( ) ] 1 1 𝑟,𝑛 1 𝑟,𝑛 𝑞̄𝑘𝑤 = 𝜒12 − , 𝑞̄4𝑤 = 𝜒10 − , ln 𝑟 + 𝜒13 𝑛𝑘 , 𝑞̄3𝑤 = 𝜒8 − 2𝜋 2𝜋 𝑟 2𝜋 𝑟 (68)

(52)

1 ∑ ′ 𝑟 ℎ𝑤 𝐾 (𝑧 ), 2𝜋 𝑚=1,2,3 2𝑚 𝑚 0 𝑚

{ [ ] 1 (1 − 2𝑣) 𝑟,𝑘 𝑛𝑖 − 𝑟,𝑖 𝑛𝑘 4𝜋(1 − 𝑣)𝑟 [ ] } − (1 − 2𝑣)𝛿𝑖𝑘 + 2𝑟,𝑖 𝑟,𝑘 𝑟,𝑛 , (𝑖, 𝑘 = 1, 2),

∑

𝑚=1,2,3

𝜒7 =

(63)

𝜒8 =

and ℎ𝑎𝑖 (𝑖 = 1, 2, 3) satisfies

𝜒9 = 86

,

[ ( ] ( ( ))) ∑ 𝑘2𝑚 ℎ𝑎𝑚 ∑ 𝑘2𝑚 ℎ𝑎𝑚 j𝑘𝑚 1 2(𝜆 + 𝜇) 𝜏̄ + ln −𝜇 , 2𝜋 2 2 2 𝑚=1,2,3 𝑚=1,2,3 ∑ 𝑚=1,2,3

⎛ 0 ⎞ ⎛ 1 ⎞ ⎛ 1 ⎞ ⎛ 1 ⎞ ⎜ 0 ⎟ = ℎ𝑎 ⎜𝑟′ ⎟ + ℎ𝑎 ⎜𝑟′ ⎟ + ℎ𝑎 ⎜𝑟′ ⎟. 21 22 1 2 3 ⎜ ⎟ ⎜ ′ ⎟ ⎜ ′ ⎟ ⎜ ′ 23 ⎟ ⎝1∕𝐴10 ⎠ ⎝𝑟 31 ⎠ ⎝𝑟 32 ⎠ ⎝𝑟 33 ⎠

) 𝑘2𝑚 ℎ𝑎𝑚

𝑟′ 2𝑚 ℎ𝑤 𝑚,

[ ( ) ] ∑ ∑ j𝑘𝑚 1 ′ 𝑤 − 𝑟′ 2𝑚 ℎ𝑤 𝜏 ̄ − 𝑟 ℎ ln , 2𝑚 𝑚 𝑚 2𝜋 2 𝑚=1,2,3 𝑚=1,2,3

L. Sun, X. Wei and B. Chen

Engineering Analysis with Boundary Elements 113 (2020) 82–98

+

𝜒14 =

] 𝑟′ 21 𝑔1 j𝑘1 𝑟′ 22 𝑔2 j𝑘2 𝑟′ 23 𝑔3 j𝑘3 ln + ln + ln , 2 2 2 2 2 2

∑ 𝑚=1,2,3

𝜒15 =

𝜒16 =

) ] [ ( ∑ ∑ j𝑘𝑚 1 ′ 𝑤 − 𝑟′ 3𝑚 ℎ𝑤 𝜏 ̄ − 𝑟 ℎ ln , 3 𝑚 𝑚 𝑚 2𝜋 2 𝑚=1,2,3 𝑚=1,2,3 ∑ 𝑚=1,2,3

𝜒17 =

Fig. 1. Problem description.

𝜒10 =

∑ 𝑚=1,2,3

𝑟′ 2𝑚 ℎ𝑎𝑚 ,

[ ( 𝜒11 =

𝜒12 =

1 − 2𝜋

∑

𝑚=1,2,3

𝜒18 = ) 𝑟′ 2𝑚 ℎ𝑎𝑚 𝜏̄ −

∑ 𝑚=1,2,3

𝑟′ 3𝑚 ℎ𝑎𝑚 ,

[ ( ) ] ∑ ∑ j𝑘 1 − 𝑟′ 3𝑚 ℎ𝑎𝑚 𝜏̄ − 𝑟′ 3𝑚 ℎ𝑎𝑚 ln 𝑚 . 2𝜋 2 𝑚=1,2,3 𝑚=1,2,3 𝑟′31 𝑔1 + 𝑟′32 𝑔2 + 𝑟′33 𝑔3 2

]

𝑟′ 2𝑚 ℎ𝑎𝑚 ln

j𝑘𝑚 , 2

[ 𝜒19

𝑟′ 21 𝑔1 + 𝑟′ 22 𝑔2 + 𝑟′ 23 𝑔3 , 2

1 =− 2𝜋 +

𝜒13 = −

𝑟′ 3𝑚 ℎ𝑤 𝑚,

[ ′ ( ) ′ ′ 1 𝑟 21 𝑔1 + 𝑟 22 𝑔2 + 𝑟 23 𝑔3 1 𝜏̄ − 2𝜋 2 2

𝑟′31 𝑔1 + 𝑟′32 𝑔2 + 𝑟′33 𝑔3 (

𝑟′33 𝑔3 2

,

In

𝑗 𝑘3 2

]

2

𝜏−

) 𝑟′ 𝑔1 𝑗 𝑘 𝑟′ 𝑔2 𝑗 𝑘 1 + 31 In 1 + 32 In 2 2 2 2 2 2

with 𝜏̄ = 0.57721566490153286. Fig. 2. Time history of (a) displacement 𝑢𝑠2 at 𝑠 (0.5, 1), (b) stress 𝜎22 at (0.5, 0), (c) pressure pa at (0.5, 0.5), and (d) pressure pw at (0.5, 1.6) calculated by the SBM.

87

L. Sun, X. Wei and B. Chen

Engineering Analysis with Boundary Elements 113 (2020) 82–98 s Fig. 3. Time history of 𝜎22 (0.5, 1) generated by the exact solutions and the SBM with different sampling number N𝜔

Fig. 4. Time history of the absolute error of 𝑢s2 (0.5, 1.5) calculated by the SBM with different 𝜅.

4. Singular boundary method 𝑠 In this study, the homogenous problems with 𝑭̃ = 0,𝐼̃𝑤 = 0,𝐼̃𝑎 = 0 are considered. The SBM approximates the solutions at a collocation point xc in the domain Ω by a linear combination of the fundamental solutions as 𝑁 𝑁 ( ) ∑ ( ) ∑ ( ) 𝑢̃ 𝑠𝑖 𝒙𝑐 = 𝛼̄ 1𝑑 𝑢̄ 𝑠𝑖1 𝒙𝑐 , 𝒚 𝑑 + 𝛼̄ 2𝑑 𝑢̄ 𝑠𝑖2 𝒙𝑐 , 𝒚 𝑑 𝑑=1

+

𝑁 ∑ 𝑑=1

𝑑=1

𝑁 ( ) ∑ ( ) 𝛼̄ 3𝑑 𝑢̄ 𝑠𝑖3 𝒙𝑐 , 𝒚 𝑑 + 𝛼̄ 4𝑑 𝑢̄ 𝑠𝑖4 𝒙𝑐 , 𝒚 𝑑 ,

(71)

𝑑=1

𝑁 𝑁 ( ) ∑ ( ) ∑ ( ) 𝑝̃𝑤 𝒙𝑐 = 𝛼̄ 1𝑑 𝑝̄𝑤 𝛼̄ 2𝑑 𝑝̄𝑤 𝒙𝑐 , 𝒚 𝑑 + 𝒙𝑐 , 𝒚 𝑑 1 2

Fig. 5. The mean absolute error of 𝑢s2 (0.5, 1.5) versus damping coefficient calculated by the SBM with different sampling number N𝜔 .

𝑑=1

+

𝑑=1

𝑁 ∑ 𝑑=1

( 𝑤

)

𝛼̄ 3𝑑 𝑝̄3 𝒙𝑐 , 𝒚 𝑑 +

𝑁 ∑ 𝑑=1

( ) 𝒙𝑐 , 𝒚 𝑑 , 𝛼̄ 4𝑑 𝑝̄𝑤 4

(72)

𝑁 𝑁 ( ) ∑ ( ) ∑ ( ) 𝑝̃𝑎 𝒙𝑐 = 𝛼̄ 1𝑑 𝑝̄𝑎1 𝒙𝑐 , 𝒚 𝑑 + 𝛼̄ 2𝑑 𝑝̄𝑎2 𝒙𝑐 , 𝒚 𝑑 𝑑=1

+

𝑁 ∑ 𝑑=1

𝑑=1

𝑁 ( ) ∑ ( ) 𝛼̄ 3𝑑 𝑝̄𝑎3 𝒙𝑐 , 𝒚 𝑑 + 𝛼̄ 4𝑑 𝑝̄𝑎4 𝒙𝑐 , 𝒚 𝑑 ,

(73)

𝑑=1

where yd is the dth source points on boundary Γ, 𝛼̄ 𝑖𝑑 (𝑖 = 1, 2, 3, 4) are the unknowns to be determined. Furthermore, the relative fluid displacements in Eqs. (16)and (17) can be obtained from Eqs. (71)–(73). To determine the unknown coefficients, boundary conditions are exerted through Eqs. (71)–(73) as Fig. 6. The mean absolute error of 𝑢s2 (0.5, 1.5) versus damping coefficient calculated by the SBM with different number of boundary points N.

𝑢̃ 𝑠𝑖 (𝒚 𝑐 ) = 88

4 ∑ 𝑁 ∑ 𝑘=1 𝑑≠𝑐

𝛼̄ 𝑘𝑑 𝑢̄ 𝑠𝑖𝑘 (𝒚 𝑐 , 𝒚 𝑑 ) +

4 ∑ 𝑘=1

𝑠 𝛼̄ 𝑘𝑐 𝑈̄ 𝑖𝑘 (𝒚 𝑐 , 𝒚 𝑐 ),

(74)

L. Sun, X. Wei and B. Chen

Engineering Analysis with Boundary Elements 113 (2020) 82–98

Fig. 7. Time history of (a) pressures pw at (0.9, 0), and (b) stress 𝜎𝑟𝑠 at (0, 0.8) calculated by the SBM.

Fig. 8. The mean absolute error of 𝑢𝑠𝑟 (0.8, 0) versus damping coefficient calculated by the SBM with different sampling number N𝜔 .

𝑞̃𝑎 (𝒚 𝑐 ) =

4 ∑ 𝑁 ∑

𝑘=1 𝑑≠𝑐 4 ∑ 𝑁 ∑

( ) 𝑝̃𝑤 𝒚 𝑐 =

𝑞̃𝑤 (𝒚 𝑐 ) =

𝑝̃𝑎 (𝒚 𝑐 ) =

𝛼̄ 𝑘𝑑 𝑡̄𝑠𝑖𝑘 (𝒚 𝑐 , 𝒚 𝑑 ) +

𝑘=1 𝑑≠𝑐

4 ∑ 𝑘=1

𝛼̄ 𝑘𝑐 𝑇̄𝑖𝑘𝑠 (𝒚 𝑐 , 𝒚 𝑐 ),

4 ( ) ∑ ( ) 𝛼̄ 𝑘𝑑 𝑝̄𝑤 𝛼̄ 𝑘𝑐 𝑃̄𝑘𝑤 𝒚 𝑐 , 𝒚 𝑐 , 𝑘 𝒚𝑐 , 𝒚𝑑 +

4 ∑ 𝑁 ∑ 𝑘=1 𝑑≠𝑐

𝛼̄ 𝑘𝑑 𝑝̄𝑎𝑘 (𝒚 𝑐 , 𝒚 𝑑 ) +

𝑟

and 𝑄𝐿 (𝒙, 𝒚 ) = − 21𝜋 𝑟,𝑛 , and traction fundamental solutions in elastostatics. Fortunately, the OIFs for abovementioned fundamental solutions are available, which can be applied to the remove the source singularities for Eqs. (74)–(79). For the OIFs of 2D Laplace problems, the analytical OIF for Dirichlet boundary condition is derived by Wei and Chen [32,34] as ( ) 𝑙𝑐 1 𝐺̂ 𝐿 (𝒚 𝑐 , 𝒚 𝑐 ) = − , (80) ln 2𝜋 2𝜋

(75) (76)

𝑘=1

4 𝑁 4 ∑ 𝜕 𝑝̃𝑤 (𝒚 𝑐 ) ∑ ∑ = 𝛼̄ 𝑘𝑑 𝑞̄𝑘𝑤 (𝒚 𝑐 , 𝒚 𝑑 ) + 𝛼̄ 𝑘𝑐 𝑄̄ 𝑤 𝑘 (𝒚 𝑐 , 𝒚 𝑐 ), 𝜕 𝒏𝑦 𝑐 𝑘=1 𝑑≠𝑐 𝑘=1 4 ∑ 𝑘=1

𝛼̄ 𝑘𝑐 𝑃̄𝑘𝑎 (𝒚 𝑐 , 𝒚 𝑐 ),

(79)

According to the nature of fundamental solutions, the basis functions in above formulas are singular when a collocation point overlaps a source point. Such singular terms requires special treatment. The SBM introduces the so-called origin intensity factors (OIFs), namely, 𝑠 (𝒚 , 𝒚 ), 𝑇̄ 𝑠 (𝒚 , 𝒚 ), 𝑃̄ 𝑤 (𝒚 , 𝒚 ), 𝑄 ̄ 𝑤 (𝒚 𝑐 , 𝒚 𝑐 ), 𝑃̄ 𝑎 (𝒚 𝑐 , 𝒚 𝑐 ) and 𝑄̄ 𝑎 (𝒚 𝑐 , 𝒚 𝑐 ) 𝑈̄ 𝑖𝑘 𝑐 𝑐 𝑐 𝑐 𝑖𝑘 𝑐 𝑐 𝑘 𝑘 𝑘 𝑘 in the Eqs. (74)–(79), to replace the singular terms in straightforward calculations. When the collocation point approaches the source point, the asymptotic behaviors in Eqs. (64)–(70) can be utilized to derive the OIFs. It can be observed that, as x → y, the fundamental solutions for partially saturated poroelastic problem are mainly associated with fundamental solutions of Laplace operator and its derivatives, i.e. 𝐺𝐿 (𝒙, 𝒚 ) = − 21𝜋 ln 𝑟

Fig. 9. The mean absolute error of pa (0,0.8) versus damping coefficient calculated by the SBM with different sampling number N𝜔 .

𝑡̃𝑠𝑖 (𝒚 𝑐 ) =

4 𝑁 4 ∑ 𝜕 𝑝̃𝑎 (𝒚 𝑐 ) ∑ ∑ = 𝛼̄ 𝑘𝑑 𝑞̄𝑘𝑎 (𝒚 𝑐 , 𝒚 𝑑 ) + 𝛼̄ 𝑘𝑐 𝑄̄ 𝑎𝑘 (𝒚 𝑐 , 𝒚 𝑐 ). 𝜕 𝒏𝑦 𝑐 𝑘=1 𝑑≠𝑐 𝑘=1

(77)

where lc is a half of the arc length between yc − 1 and yc+1 on the boundary. The OIFs for Neumann boundary conditions are derived by subtracting and adding-back technique [56] as:

(78)

𝑄̂ 𝐿 (𝒚 𝑐 , 𝒚 𝑐 ) = − 89

𝑁 ∑ 𝑙𝑑 𝜕 𝐺𝐿 (𝒚 𝑐 , 𝒚 𝑑 ) . 𝑙 𝜕 𝒏𝒚 𝑑 𝑑≠𝑐 𝑐

(81)

L. Sun, X. Wei and B. Chen

Engineering Analysis with Boundary Elements 113 (2020) 82–98

Fig. 12. Annular-sector domain problem. Fig. 10. Ring-shaped domain problem.

[ ( ) ] ( ) 𝑤 ( ) ̂ ̄𝑤 ̂ ̄ ̂ 𝑄̄ 𝑤 𝑘 = 𝜒12 𝐺𝐿 𝒚 𝑐 , 𝒚 𝑐 + 𝜒13 𝑛𝑘 , 𝑄3 = 𝜒8 𝑄𝐿 𝒚 𝑐 , 𝒚 𝑐 , 𝑄4 = 𝜒10 𝑄𝐿 𝒚 𝑐 , 𝒚 𝑐 ,

For elastostatics, the traction fundamental solutions are given by 𝑡̂𝑖𝑘 =

1 4𝜋(1 − 𝑣)𝑟

{ [ ] [ ] } (1 − 2𝑣) 𝑟,𝑘 𝑛𝑖 − 𝑟,𝑖 𝑛𝑘 − (1 − 2𝑣)𝛿𝑖𝑘 + 2𝑟,𝑖 𝑟,𝑘 𝑟,𝑛 .

(90)

(82)

The OIFs [57] are 𝑁 ∑ ( ) ) 𝑙𝑑 ( 𝑇̂𝑖𝑘 𝒚 𝑐 , 𝒚 𝑐 = − 𝑡̂𝑖𝑘 𝒚 𝑑 , 𝒚 𝑐 . 𝑙 𝑐 𝑑=1

𝑃̄𝑘𝑎 = 0, 𝑃̄3𝑎 = 𝜒14 𝐺̂ 𝐿 (𝒚 𝑐 , 𝒚 𝑐 ) + 𝜒15 , 𝑃̄4𝑎 = 𝜒16 𝐺̂ 𝐿 (𝒚 𝑐 , 𝒚 𝑐 ) + 𝜒17 , (𝑘 = 1, 2), (91)

(83)

Besides the above singular terms, the non-singular term r, i r, k in Eq. (64) also needs to be determined. As described in [57], this term at yc is related to the tangential vector (𝜏̃1 , 𝜏̃2 ) = (𝑛2 , −𝑛1 ), and can be given as 𝜏̂𝑖𝑘 = 𝑟,𝑖 𝑟,𝑘 =

𝜏̃𝑖 𝜏̃𝑘 𝜏̃12 + 𝜏̃22

.

[ ( ) ] ( ) ( ) 𝑄̄ 𝑎𝑘 = 𝜒18 𝐺̂ 𝐿 𝒚 𝑐 , 𝒚 𝑐 + 𝜒19 𝑛𝑘 , 𝑄̄ 𝑎3 = 𝜒14 𝑄̂ 𝐿 𝒚 𝑐 , 𝒚 𝑐 , 𝑄̄ 𝑎4 = 𝜒16 𝑄̂ 𝐿 𝒚 𝑐 , 𝒚 𝑐 . (92) With the computed OIFs, the unknowns in Eqs. (74) to (79) can be solved by the collocation method, and subsequently the variables within the domain can be determined.

(84)

Applying Eqs. (80)–(84) to Eqs. (64)–(70), the OIFs for partially saturated poroelastic media problems are [ ] 𝑠 𝑈̄ 𝑖𝑘 = 𝜒1 𝐺̂ 𝐿 (𝒚 𝑐 , 𝒚 𝑐 ) 𝛿𝑖𝑘 − 𝜒2 𝛿𝑖𝑘 + 𝜒3 𝜏̂𝑖𝑘 , (𝑖, 𝑘 = 1, 2), (85) 𝑈̄ 𝑖𝑠3 = 0, 𝑈̄ 𝑖𝑠4 = 0, (𝑖 = 1, 2),

(86)

( ) 𝑇̄𝑖𝑘𝑠 = 𝑇̂𝑖𝑘 𝒚 𝑐 , 𝒚 𝑐 , (𝑖, 𝑘 = 1, 2),

(87)

[ ] [ ] 𝑇̄𝑖𝑠3 = 𝜒4 𝐺̂ 𝐿 (𝒚 𝑐 , 𝒚 𝑐 ) + 𝜒5 𝑛𝑖 , 𝑇̄𝑖𝑠4 = 𝜒6 𝐺̂ 𝐿 (𝒚 𝑐 , 𝒚 𝑐 ) + 𝜒7 𝑛𝑖 , (𝑖 = 1, 2),

(88)

5. Exponential window method For transient problems, the inverse fast Fourier transformation (FFT) [58,59] is attractive in transforming the frequency-domain solutions at a series of discrete frequencies into time-domain solutions. Nevertheless, it fails to deal with non-dissipative system, and is computationally expensive for lightly damped systems. To remedy this problem, the exponential window method (EWM) [60] is introduced to solve undamped (or lightly damped) systems. In this technique, the frequency is shifted to create an artificial damping to damp out the free vibrations within the selected period. The procedure for the EWM can be summarized as: 𝑠 , 𝑢𝑤 , 𝑝𝑤 , 𝑢𝑎 and 𝑝𝑎 ) with (I) Scale the unknown variables 𝜉(𝒙, 𝑡)(𝑢𝑠𝑖 , 𝜎𝑖𝑘 𝑖 𝑖 − 𝜀 t an exponential window function e to obtain a damped variable − 𝜀 t as 𝜉 ew (x,t) = 𝜉(x, t)e , where 𝜀 is the shifting constant. Take 𝜉(x,

( ) ( ) 𝑃̄𝑘𝑤 = 0, 𝑃̄3𝑤 = 𝜒8 𝐺̂ 𝐿 𝒚 𝑐 , 𝒚 𝑐 + 𝜒9 , 𝑃̄4𝑤 = 𝜒10 𝐺̂ 𝐿 𝒚 𝑐 , 𝒚 𝑐 + 𝜒11 , (𝑘 = 1, 2), (89)

Fig. 11. Time history of (a) displacement 𝑢𝑠𝑟 at (1, 0), and (b) fluid pressure pw at (0.4, 0.45) calculated by the SBM.

90

L. Sun, X. Wei and B. Chen

Engineering Analysis with Boundary Elements 113 (2020) 82–98

Fig. 13. Time history of (a) radial displacement 𝑢𝑠𝑟 at (0.4, 0.4), and (b) gas pressure pa at (0.55, 0.7) calculated by the SBM. Fig. 14. The radial stress 𝜎𝑟𝑠 along the radius calculated by the SBM at t = 1.25ms, 2.5ms and 3.75ms for the annular-sector domain problem.

Fig. 15. Problem description.

t) = 𝜉 ew (x,t)e𝜀t into the governing Eqs. (8) to (10), and obtain a frequency-domain boundary value problem (18) to (20) with 𝜔̄ = 𝜔 − j𝜀 as angular frequency and 𝜉̃𝑒𝑤 (𝒙, 𝜔), the Fourier transform of 𝜉 ew (x,t), as the variables to be determined. (II) Let T be the total time duration for the analysis with the given number of sampling frequencies N𝜔 , then the time and circular frequency resolutions are Δ𝑡 = 𝑇 ∕𝑁𝜔 and Δ𝜔 = 2𝜋∕𝑇 .

(III) Choose a suitable shifting constant based on the rule-of-thumb as 𝜀=

𝜅 ln 10 , 𝑇

(93)

where 𝜅 is damping coefficient, and 2 ≤ 𝜅 ≤ 3 is recommended. (IV) Scale the boundary conditions P(x, t) with exponential window function Pew (x,t) = e−𝜀t P(x,t) and then transform the scaled boundary conditions into frequency domain via discretized 91

L. Sun, X. Wei and B. Chen

Engineering Analysis with Boundary Elements 113 (2020) 82–98

Fig. 16. Time history of (a) displacement 𝑢𝑠2 at (0.5, 1), (b) pressure pw at (0.5, 0.8), (c) pres𝑠 sure pa at (0.5, 0.2), and (d) stress 𝜎22 at (0.5, 0) calculated by the SBM for the first kind of boundary conditions.

Fourier transform (DFT) as

Appendix D. For the problems with regular geometries and simple boundary conditions, the exact solutions are derived to compare with the SBM results. To evaluate the accuracy of the present method, the mean absolute error of variable 𝜉 versus time at point x is defined as

𝑁𝜔 −1

( ) 1 ∑ 𝑃̃𝑒𝑤 𝒙, 𝜔̄ 𝑘 = 𝑃 (𝒙, 𝑚Δ𝑡)𝑒−2𝜋j𝑚𝑘∕𝑁𝜔 𝑁𝜔 𝑚=0 𝑒𝑤 =

𝑁𝜔 −1 1 ∑ −𝜀𝑚Δ𝑡 𝑒 𝑃 (𝒙, 𝑚Δ𝑡)𝑒−2𝜋j𝑚𝑘∕𝑁𝜔 , 𝑁𝜔 𝑚=0

(94)

𝑀𝐴𝐸(𝜉) =

where 𝜔̄ 𝑘 = 𝑘Δ𝜔 − j𝜀(𝑘 = 0, 1, … , 𝑁𝜔 − 1). (V) Solve the frequency-domain problems at the frequencies 𝜔̄ 𝑘 (𝑘 = 0, 1, … , 𝑁𝜔 ∕2) to obtain the frequency-domain solutions 𝑅∗𝑒𝑤 (𝜔̄ 𝑘 ). The results for the remaining frequencies could be determined by the conjugate symmetric property of the DFT as ( ) ( ( )) 𝑅∗𝑒𝑤 𝜔̄ 𝑁𝜔 −𝑘 = conj 𝑅∗𝑒𝑤 𝜔̄ 𝑘 , 𝑘 = 𝑁𝜔 ∕2 + 1, 𝑁𝜔 ∕2 + 2, … , 𝑁𝜔 − 1.

6.1. Wave propagation in one dimensional partially saturated poroelastic column

(VI) Utilize the inverse DFT to transfer the frequency-domain solutions into time-domain solutions as: 𝑁𝜔 −1

∑

𝑘=0

𝑊𝑘 𝑅∗𝑒𝑤 (𝜔̄ 𝑘 )𝑒2𝜋j𝑚𝑘∕𝑁𝜔 ,

In the first case, a plane domain with dimensions of 1 m × 2 m (Fig. 1) is considered. A vertical uniform impact force is exerted on the upper boundary and the rest boundaries are sliding: {𝑠 𝑡𝑛 = −𝐻(𝑡)N∕m2 , 𝑡𝜏𝑠̃ = 0N∕m2 , 𝑝𝑤 = 𝑝𝑎 = 0Pa, on the top boudary , 𝑎 𝑢𝑠𝑛 = 0m, 𝑡𝜏𝑠̃ = 0N∕m2 , 𝑢𝑤 𝑛 = 𝑢𝑛 = 0m, on the other boudaries

(96)

where 𝑊𝑘 = 0.5[1 + cos(2𝜋𝑘∕𝑁𝜔 )] is Hanning window functions to alleviate the Gibbs oscillations. (VII) Remove the influence of the damping exponential window function in the first step, and the true time-domain solutions are given by 𝑅(𝑚Δ𝑡) = 𝑒𝜀𝑚Δ𝑡 𝑅𝑒𝑤 (𝑚Δ𝑡).

(98)

where 𝜉 h and 𝜉 are numerical and exact solutions, respectively. All calculations are performed on a single Intel Core(TM) processor running at 2.4GHz with 6MB L3 cache of a quad core CPU workstation having a total of 32GB DDR3 memory. The SBM are implemented in Matlab R2014a.

(95)

𝑅𝑒𝑤 (𝑚Δ𝑡) =

( ( ) )| 𝑁𝜔 −1 1 ∑ || ℎ 𝑘𝑇 𝑘𝑇 | ,𝒙 − 𝜉 , 𝒙 |, |𝜉 | 𝑁𝜔 𝑘=0 || 𝑁𝜔 𝑁𝜔 |

where

𝑡𝑠𝑛 = 𝑡𝑠1 𝑛1 + 𝑡𝑠2 𝑛2 , 𝑢𝑓𝑛

𝑡𝜏𝑠̃ = 𝑡𝑠1 𝜏̃1 + 𝑡𝑠2 𝜏̃2 ,

𝛕̃ = (𝜏̃1 , 𝜏̃2 ) = (𝑛2 , −𝑛1 ),

𝑢𝑠𝑛 = 𝑢𝑠1 𝑛1 + 𝑢𝑠2 𝑛2 , = + 𝑢𝑓2 𝑛2 (f = w, a) and H(t) is Heaviside step function. The exact solution for this case is available in [8]. The total time duration T for the analysis is 15ms. In this case, the boundary is discretized into 450 boundary nodes. The parameters N𝜔 and 𝜅 for the EWM are selected as 128 and 3. 𝑠 at (0.5, Fig. 2 displays the solid displacement 𝑢𝑠2 at (0.5, 1), stress 𝜎22 0), gas pressure pa at (0.5, 0.5) and fluid pressure pw at (0.5, 1.6). The CPU times for the results in Fig. 2a obtained by the numerical method

(97)

6. Numerical examples In this section, four examples are carried out to verify the validity of the SBM. The material properties for all cases are provided in 92

𝑢𝑓1 𝑛1

L. Sun, X. Wei and B. Chen

Engineering Analysis with Boundary Elements 113 (2020) 82–98 𝑠 Fig. 17. Distribution of stress 𝜎12 [N/m2 ] within the domain at different times for the first kind of boundary conditions.

this part, we choose 𝜅 = 1.5, 2, 2.5, 3 and 3.5 to investigate the effect of the parameter 𝜅. Fig. 4 presents the time response of the absolute error of displacement 𝑢𝑠2 at (0.5, 1.5) with respect to different parameters 𝜅 with 450 boundary points and N𝜔 = 128. In the first 10ms, the solutions for 𝜅 = 2.5, 3 and 3.5 are better than 𝜅 = 1.5 and 2. Nevertheless, the solutions for 𝜅 = 3 and 3.5 exhibit obvious increased errors in the last 5ms. This phenomenon is caused by the fact that when the parameter 𝜅 and the time variable t are large, the errors will be amplified by the function e𝜀t in Eq. (97). The mean absolute errors of 𝑢s2 (0.5, 1.5) in terms of damping coefficients 𝜅 with different number of sampling frequencies and boundary points are shown in Fig. 5 and Fig. 6. The figures reveal that the best damping coefficient lies in intervals [2.3, 2.9] for different cases. Furthermore, it can be observed that although the optimal 𝜅 is different for different cases, 𝜅= 2.5 is a suitable choice because the difference between the error for 𝜅= 2.5 and the minimum error is small.

without and with Hanning window functions are 182s and 183s, respectively. Correspondingly, the CPU time taken by the SBM for each discrete frequency 𝜔̄ 𝑘 = 𝑘Δ𝜔 − j𝜀 is about 2.8s. Furthermore, the CPU times for the results in Fig. 2b–d are the same as Fig. 2a, because the solutions of Fig. 2 are obtained by using the same program. For the results without Hanning window functions in Eq. (96), there exist oscillations in later times. Thus, the frequency windowing technique, Hanning window, is introduced to remove the oscillations. All results from the present method with Hanning window well agree with the exact solutions. The Hanning window will be applied to all of the following cases. Furthermore, the parameters N𝜔 and 𝜅 in the EWM also have a great effect on the accuracy of the numerical solutions. N𝜔 is the number of sampling frequencies in the DFT. According to Fourier sampling theory, it can be known that more sampling frequencies will produce more accurate results, and also bring more time cons at (0.5, 1) with different sumption. In Fig. 3, the numerical results 𝜎22 N𝜔 (64, 128, 256) are studied with 450 boundary points and 𝜅 = 2.5. From Fig. 3, we can see that the results converge to the analytical solution as the number of samplings increases. The reason lies in that more sampling frequencies enhance the resolution of the time-history curves via reducing Δt. To make a trade-off between the time consumption and accuracy, N𝜔 = 128 is an appropriate choice. 𝜅 is the damping parameter, which determines the shifting constant 𝜀. The shifting constant 𝜀 should not be too large or too small. A small shifting constant will produces a lightly damped system. A large shifting constant will lead either to severe loss of numerical precision, or to underflows and overflows in the evaluation of the exponential windows. In

6.2. Circular partially saturated poroelastic domain In the second example, a circle domain with the radius of 1m is considered. On the whole boundary, the following boundary conditions are exerted: 𝑡𝑠𝑛 = −𝐻(𝑡)N∕m2 , 𝑡𝜏𝑠̃ = 0N∕m2 , 𝑝𝑤 = 𝑝𝑎 = 0Pa. To compare and validate our numerical methods, we derive the analytic frequency-domain solutions for this case in Appendix E. The timedomain solutions can be obtained by using the EWM. 93

L. Sun, X. Wei and B. Chen

Engineering Analysis with Boundary Elements 113 (2020) 82–98

Fig. 18. Time history of (a) displacement 𝑢𝑠2 at (0.5, 1), (b) pressure pw at (0.5, 0.8), (c) pres𝑠 sure pa at (0.5, 0.2), and (d) stress 𝜎22 at (0.5, 0) calculated by the SBM for the second kind of boundary conditions.

The boundary is discretized into 400 boundary points. T= 5ms. The parameters in the EWM are selected as N𝜔 = 128, 𝜅 = 2.5. Fig. 7 displays the time history of the fluid pressure pw at (0.9,0) and the stress 𝜎𝑟𝑠 at (0, 0.8). The CPU time taken by the present method for each figure is 130s. Good agreements can be observed between the SBM and analytical solutions. Similar to the first example, we investigate the optimal choice for the damping coefficient 𝜅. Figs. 8 and 9 respectively present the mean absolute error of radial displacement ur at (0.8, 0) and gas pressure pa at (0,0.8) with respect to different 𝜅. From figures, the optimal 𝜅 lies in interval [2.1, 2.7] with different numbers of sampling frequencies. Based on the results obtained in last and this examples, it can be found that 𝜅 = 2.5 is a suitable selection for different cases. Furthermore, we consider a multi-connected ring-shaped domain, a disk of radius 1m with a circle cutout of radius 0.5m, as shown in Fig. 10. The boundary conditions are:

𝑡𝑠𝑛 = −𝐻(𝑡)N∕m2 ,𝑡𝜏𝑠̃ = 0N∕m2 , 𝑝𝑤 = 𝑝𝑎 = 0Pa on Γ2 , 𝑡𝑠𝑛 = 𝑡𝜏𝑠̃ = 0N∕m2 , 𝑝𝑤 = 𝑝𝑎 = 0Pa on Γ4 . With the exact solution in frequency-domain derived in Appendix E, the exact solution in time-domain can be obtained by using the EWM. The SBM employs 269 boundary points. T=5ms, N𝜔 = 256 and 𝜅 = 2.5. Time history of radial displacement 𝑢𝑠𝑟 at (0.4, 0.4) and gas pressure pa at (0.55, 0.7) are depicted in Fig. 13. The CPU time for each figure is 142s. It can be observed that the numerical solutions match the exact solutions well. Furthermore, the radial stresses 𝜎𝑟𝑠 along the radius at t = 1.25 ms, 2.5 ms and 3.75 ms are given in Fig. 14. From the figures, it can be found that the SBM can accurately predict the dynamic response of partially saturated media within complex domains.

𝑡𝑠𝑛 = −𝐻(𝑡)N∕m2 ,𝑡𝜏𝑠̃ = 0N∕m2 , 𝑝𝑤 = 𝑝𝑎 = 0Pa on outer boundary, 𝑡𝑠𝑛 = 𝑡𝜏𝑠̃ = 0N∕m2 , 𝑝𝑤 = 𝑝𝑎 = 0Pa on inner boundary.

6.4. Two dimensional partially saturated poroelastic domain exerted by discontinuous load

The exact frequency-domain solution can be found in Appendix E. 600 boundary points are employed in the SBM. T= 5ms. N𝜔 = 256, 𝜅 = 2.5 are used in the EWM. Fig. 11 plots the time history of the displacement 𝑢𝑠𝑟 at (1,0) and pressure pw at (0.4,0.45) . It takes 540s for the results in each figure. The agreements between the exact solutions and SBM results indicate that the proposed method is applicable to multiconnected domain problem.

A unit square domain ([0, 1] × [0, 1]) subjected to non-smooth boundary conditions, as shown in Fig. 15, is studied. Two kinds of boundary conditions are investigated: 𝑠 2 𝑠 2 𝑤 𝑎 ⎧𝑡𝑛 = −𝑓 (𝑥1 )𝑔 (𝑡)N∕m , 𝑡𝜏̃ = 0N∕m , 𝑝 = 𝑝 = 0Pa, on the top boudary ⎪ (I) ⎨𝑓 (𝑥1 ) = 1 − |1 − 2𝑥1 |, 𝑔 (𝑡) = 𝐻(𝑡) − 𝐻(𝑡 − 0.001) , ⎪ 𝑠 𝑠 2 𝑤 𝑎 ⎩𝑢𝑛 = 0m, 𝑡𝜏̃ = 0N∕m , 𝑢𝑛 = 𝑢𝑛 = 0m, on the other boudaries { ⎧ −𝐻(𝑡)N∕m2 , 𝑥1 < 0.25 ⎪𝑡𝑠 = , 𝑡𝜏𝑠̃ = 0N∕m2 , 𝑝𝑤 = 𝑝𝑎 𝑛 ⎪ 0N∕m2 , 𝑥1 ≥ 0.25 (II) ⎨ , = 0Pa, on the top boudary ⎪ ⎪𝑢𝑠 = 0m, 𝑡𝑠 = 0N∕m2 , 𝑢𝑤 = 𝑢𝑎 = 0m, on the other boudaries 𝑛 𝑛 𝜏̃ ⎩ 𝑛

6.3. Annular-sector partially saturated poroelastic domain In this example, we consider the problem with annular-sector domain as shown in Fig. 12. The boundary conditions are given by: 𝑎 𝑢𝑠𝑛 = 0m,𝑡𝜏𝑠̃ = 0N∕m2 , 𝑢𝑤 𝑛 = 𝑢𝑛 = 0m on Γ1 , Γ3 ,

94

L. Sun, X. Wei and B. Chen

Engineering Analysis with Boundary Elements 113 (2020) 82–98

Fig. 19. Distribution of displacement 𝑢𝑠2 [10−10 m] within the domain at different times for the second kind of boundary conditions.

consistent with each other. The displacement 𝑢𝑠2 and deformation at different time are shown in Fig. 19. At first, the deformation on the left of the top boundary is more obvious. Then the deformation propagates to the right side of the domain. The reason lies in that the load is exerted on the left of the top boundary. With respect to the boundary conditions, the numerical solutions obtained from the SBM in combination with the EWM for two cases seem reasonable.

In this case, there are no available analytical solutions owing to the complexity of the boundary conditions. In order to illustrate the accuracy of numerical results, the SBM is carried out with different numbers of sampling frequencies (N𝜔 = 128, 256) in the FFT and different numbers of boundary nodes (N = 340, 660) in the SBM. T = 8 ms. The damping coefficient 𝜅 is set as 2.5 for two types of boundary conditions. For the boundary condition (I), a symmetric triangular load, discontinuous with respect to time, is exerted on the top boundary. The time responses of displacement 𝑢𝑠2 at (0.5, 1), pressure pw at (0.5, 0.8), pres𝑠 at (0.5, 0) computed by the SBM with sure pa at (0.5, 0.2) and stress 𝜎22 different numbers of boundary nodes, and the EWM with different numbers of frequency sampling points are plotted in Fig. 16. The CPU times for the results in each figure are 736s, 382s and 104s for the cases with N𝜔 = 256, N = 660, N𝜔 = 128, N = 660 and N𝜔 = 128, N = 340, respectively. From the figures, the SBM in conjunction with the EWM with different parameters and boundary nodes achieves nearly the same re𝑠 and deformation of the domain sults. The contour plots of the stress 𝜎12 at different time are displayed in Fig. 17. In the figures, the deformation of the domain is symmetric and the distribution of the shear stresses is anti-symmetric, which agrees well with the boundary conditions. For boundary condition (II), the boundary condition on the top boundary is discontinuous at x = 0.25. Fig. 18 displays the time history of displacement 𝑢𝑠2 at (0.5, 1), pressure pw at (0.5, 0.8), pressure pa 𝑠 at (0.5, 0). The CPU times are 734s, 379s and at (0.5, 0.2) and stress 𝜎22 104s for N𝜔 = 256, N = 660, N𝜔 = 128, N = 660, and N𝜔 = 128, N = 340, respectively. The results of the SBM-EWM with different parameters are

7. Conclusion In this study, the SBM in combination with the EWM (SBM-EWM) is formulated to deal with the 2D dynamic analysis of the partially saturated poroelastic media. The SBM solves the problem in the frequencydomain. In the SBM formulation, two difficulties are solved: 1) The fundamental solutions are derived, which are the basis functions for the SBM. 2) The analytical OIFs are derived to isolate the source singularities of fundamental solutions, which is the key issue in the SBM. Then the SBM solutions are transformed to the time-domain by the EWM, where the artificial damping and Hanning window are introduced. The damping enhances the performance of the fast Fourier transformation in the undamped or lightly damped systems. The Hanning window removes the Gibbs oscillations in the numerical solutions. 95

L. Sun, X. Wei and B. Chen

Engineering Analysis with Boundary Elements 113 (2020) 82–98

The SBM-EWM is tested to four numerical examples with different geometries and boundary conditions, meanwhile, the parameters in the EWM are investigated. Numerical results verified the accuracy and feasibility of the SBM-EWM with different parameters.

Table A.2 The symbols and notations used in fundamental solutions. Description 𝑢̄ 𝑠𝑖𝑘 𝑡̄𝑠𝑖𝑘 𝑝̄𝑤 𝑝̄𝑎 ek AL AT Aw Aa 𝜏 Mi (i = 1, 2) M g gi (i = 1, 2, 3, 4) 𝜙i (i = 1, 2, 3) r′2i (i = 1, 2, 3) r′3i (i = 1, 2, 3) 𝑘𝑖 (𝑖 = 1, 2, 3, 4) Ki (i = 0,1,2,…) A,B,C,D 𝜆 𝜇 𝜑 𝝍 hw ℎ𝑤 𝑖 (𝑖 = 1, 2, 3) Aw Aa ℎ𝑎𝑖 (𝑖 = 1, 2, 3) 𝜒 i (i = 1, 2, …, 17)

Acknowledgments The work described in this paper was supported by the National Natural Science Foundation of China (Grant Nos. 11602114, 11662003, 11762008, 11771225, 11802151). Appendix A. Symbols and notations

Table A.1 The main symbols and notations used in governing equations. Description us 𝝈s uw ua pw pa 𝝈 ∇ Sw Sa K G 𝛿 ij 𝛼 Ks 𝜌s 𝜌w 𝜌a 𝜌 Fs n 𝜅w 𝜅a 𝑘̄ Krw Kra 𝜂w 𝜂a ϑ Se Srw Sra 𝜁 Saa Sww Su pc Iw Ia n 𝒖̃ 𝑠 𝒖̃ 𝑤 𝒖̃ 𝑎 𝑝̃𝑤 𝑝̃𝑎 𝜔 𝛽 𝛾

Displacement vector of the solid skeleton Stress tensor of the solid skeleton Relative displacement vector of pore water Relative displacement vector of pore air Pressure of pore water Pressure of pore air Total stress tensor Gradient operator Degree of saturation of water Degree of saturation of air Drained bulk modulus of the mixture Shear modulus Kronecker delta The coefficient used to describe the compressibility of the solid grains Bulk modulus of the solid skeleton Density of solid skeleton Density of pore water Density of pore air Averaged density of the mixture Bulk body force Porosity Phase permeability of water Phase permeability of air Intrinsic fluid permeability Relative water phase permeability Relative air phase permeability Viscosity of the water Viscosity of the air Pore size distribution index Effective wetting fluid saturation Residual wetting fluid saturation Non-wetting fluid entry saturation Abbreviation used in Eqs. (9) and (10) Abbreviation used in Eqs. (9) and (10) Abbreviation used in Eqs. (9) and (10) Abbreviation used in Eqs. (9) and (10) Capillary pressure Wetting fluid source term Non- wetting fluid source term Outward normal vector Displacement vector of the solid skeleton in frequency domain Relative displacement vector of pore water in frequency domain Relative displacement vector of pore air in frequency domain Pressure of pore water in frequency domain Pressure of pore air in frequency domain Angular frequency Abbreviation used in Eq. (16) Abbreviation used in Eq. (17)

Displacement fundamental solutions of solid phase Traction fundamental solutions of solid phase Pressure fundamental solutions of water Pressure fundamental solutions of air Unit vector along xk direction Variable in Eq. (21) Variable in Eq. (21) Variable in Eq. (22) Variable in Eq. (22) Fundamental solution of the Laplace operator Matrix in Eq. (29) Matrix in Eq. (30) Abbreviation used in Eq. (30) Coefficient in Eq. (33) and (36) Variable in Eq. (32) Coefficient in Eq. (32) Coefficient in Eq. (32) Wave number Modified Bessel function of ith order Variables in Eqs. (37)–(39) First Lamé coeffificient Second Lamé coeffificient Variable in Eq. (43) Variable in Eq. (43) Abbreviation used in Eq. (48) Coefficient in Eq. (50) Variable in Eq. (52) Variable in Eq. (58) Coefficient in Eq. (59) Coefficient in Eqs. (64)–(70)

Table A.3 The symbols and notations used in the singular boundary method. Description 𝛼̄ 𝑖𝑑 (𝑖 = 1, 2, 3, 4) 𝑈̄ 𝑖𝑘𝑠 𝑇̄𝑖𝑘𝑠 𝑃̄𝑘𝑤 𝑄̄ 𝑤 𝑘 𝑃̄𝑘𝑎 ̄ 𝑄𝑎 𝑘

lc 𝐺̂ 𝐿 𝑄̂ 𝐿 𝑡̂𝑖𝑘 𝑇̂𝑖𝑘 𝝉̃ = (𝜏̃1 , 𝜏̃2 ) 𝜏̂𝑖𝑘

Coefficient in Eqs. (71)–(73) Origin intensity factor in Eq. (74) Origin intensity factor in Eq. (75) Origin intensity factor in Eq. (76) Origin intensity factor in Eq. (77) Origin intensity factor in Eq. (78) Origin intensity factor in Eq. (79) Half of the arc length between yc − 1 and yc+1 Origin intensity factor for Dirichlet boundary condition of Laplace problem Origin intensity factor for Neumann boundary conditions of Laplace problem Traction fundamental solution of elastostatic problem Origin intensity factor for traction boundary conditions of elastostatic problem Tangential vector Abbreviation used in Eq. (84)

Table A.4 The symbols and notations used in the exponential window method. Description 𝜀 N𝜔 T Δt Δ𝜔 𝜅 𝜔̄ 𝑘 Wk

96

Shifting constant Number of sampling frequencies Total time duration Time resolution Circular frequency resolution Damping coefficient Discrete complex angular frequency Hanning window function

L. Sun, X. Wei and B. Chen

Engineering Analysis with Boundary Elements 113 (2020) 82–98

Appendix B. The parameters of Section 2.2

Appendix E. Exact solutions for the problems in frequency domain

The following abbreviations are used in Section 2.2: I) For the circle domain problem with radius of R and the boundary conditions: 𝑡𝑠𝑛 = −1N∕m2 , 𝑡𝜏𝑠̃ = 0N∕m2 , pw = pa = 0Pa, that is, 𝜎𝑟𝑠 = −1N∕m2 , pw = pa = 0Pa, the exact solutions in polar coordinate system can be expressed as:

𝐴1 = ρ − 𝑛𝑆𝑤 ρ2𝑤 𝛽𝜔2 − 𝑛𝑆𝑎 ρ2𝑎 𝛾𝜔2 , 𝐴2 = 𝛼𝑆𝑤 − 𝜔2 𝑛𝑆𝑤 ρ𝑤 𝛽, 𝐴3 = 𝛼𝑆𝑎 − 𝜔2 𝑛𝑆𝑎 ρ𝑎 𝛾, ( ) 𝐴4 = j𝜔 𝛼𝑆𝑤 − 𝑛𝑆𝑤 𝛽ρ𝑤 𝜔2 , 𝐴5 = j𝜔𝑛𝑆𝑤 𝛽, ( ) ( ) 𝑛 𝐴6 = j𝜔 𝜁 𝑆𝑤𝑤 𝑆𝑤 + 𝑆𝑤 − 𝑆𝑢 , 𝐴7 = j𝜔 𝜁 𝑆𝑎𝑎 𝑆𝑤 + 𝑆𝑢 , 𝐾𝑤 ( ) ( ) 𝐴8 = j𝜔 𝛼𝑆𝑎 − 𝑛𝑆𝑎 𝛾𝜌𝑎 𝜔2 , 𝐴9 = j𝜔 𝜁 𝑆𝑤𝑤 𝑆𝑎 + 𝑆𝑢 , ( ) 𝑛 𝐴10 = j𝜔𝑛𝑆𝑎 𝛾, 𝐴11 = j𝜔 𝜁 𝑆𝑎𝑎 𝑆𝑎 + 𝑆 − 𝑆𝑢 . 𝐾𝑎 𝑎

𝐴1 𝐴5 𝐴10 𝜔2 + 𝐴2 𝐴4 𝐴10 + 𝐴3 𝐴5 𝐴8 + (𝐾 + 4∕3𝐺)(𝐴5 𝐴11 + 𝐴6 𝐴10 ) , (𝐾 + 4∕3𝐺)𝐴5 𝐴10 ( ) 𝜔2 𝐴1 𝐴5 𝐴11 + 𝐴1 𝐴6 𝐴10 + 𝐴2 𝐴4 𝐴11 − 𝐴3 𝐴4 𝐴9 − 𝐴2 𝐴7 𝐴8 + 𝐴3 𝐴6 𝐴8

𝐶1 = −

(𝐾 + 4∕3𝐺)𝐴5 𝐴10 ( ) (𝐾 + 4∕3𝐺) 𝐴6 𝐴11 − 𝐴7 𝐴9 + , (𝐾 + 4∕3𝐺)𝐴5 𝐴10 ( ) 𝜔2 𝐴1 𝐴7 𝐴9 − 𝐴1 𝐴6 𝐴11 𝐶3 = , (𝐾 + 4∕3𝐺)𝐴5 𝐴10 √ 3 𝐶1 2 , 𝑁2 = , 3 𝑁3 √ 3

−2𝐶13 + 9𝐶1 𝐶2 − 27𝐶3 +

√ 𝑘1 =

𝑁1 + √

𝑘2 =

𝑁1 + √

𝑘3 =

𝑁1 +

𝑁2 𝐶12 3

6 3𝐶2 − 𝐶12 6

√ 2 3 4(−𝐶12 + 3𝐶2 ) + (−2𝐶13 + 9𝐶1 𝐶2 − 27𝐶3 ) ,

√ ) 1−𝑖 3 −

( √ ) 1 1+𝑖 3 , 6𝑁2

( √ ) 𝑁2 1 + 𝑖 3 −

( √ ) 1 1−𝑖 3 . 6𝑁2

𝑁2

(

(E.2)

𝑝𝑎 = 𝑎1 𝑟31 𝐽0 (𝑘1 𝑟) + 𝑎2 𝑟32 𝐽0 (𝑘2 𝑟) + 𝑎3 𝑟33 𝐽0 (𝑘3 𝑟),

(E.3)

⎡ Θ11 ⎢ ⎢𝑟21 𝐽0 (𝑘1 𝑅) ⎢ ⎣𝑟31 𝐽0 (𝑘1 𝑅)

1 − 𝑁2 𝐶2 + , 3𝑁2

3𝐶2 − 𝐶12

𝑝𝑤 = 𝑎1 𝑟21 𝐽0 (𝑘1 𝑟) + 𝑎2 𝑟22 𝐽0 (𝑘2 𝑟) + 𝑎3 𝑟23 𝐽0 (𝑘3 𝑟),

To obtain the specific forms of the exact solutions, the unknown coefficients can be determined through enforcing Eqs. (E.2)–(E.4) to satisfy the boundary conditions, that is, solving the following linear system of equations:

𝑁1 = − 𝑁3 =

(E.1)

where Jk denotes the Bessel function of the first kind of order k, a1 , a2 and a3 are unknown coefficients, ki , r2i and r3i (i=1, 2, 3) are parameters used in the fundamental solutions. With the help of the relationship between the displacements and stress, we can get the following stress [ ] 𝜎𝑟𝑠 = (𝜆 + 2𝜇) −𝑎1 𝑘21 𝐽0 (𝑘1 𝑟) − 𝑎2 𝑘22 𝐽0 (𝑘2 𝑟) − 𝑎3 𝑘23 𝐽0 (𝑘3 𝑟) [ ] 𝐽 (𝑘 𝑟) 𝐽 (𝑘 𝑟) 𝐽 (𝑘 𝑟) + 2𝜇 𝑎1 𝑘1 1 1 + 𝑎2 𝑘2 1 2 + 𝑎3 𝑘3 1 3 𝑟 𝑟 𝑟 ( ) 𝐽1 (𝑘1 𝑟) = 𝑎1 2𝜇𝑘1 − (𝜆 + 2𝜇)𝑘21 𝐽0 (𝑘1 𝑟) 𝑟 ( ) 𝐽 (𝑘 𝑟) + 𝑎2 2𝜇𝑘2 1 2 − (𝜆 + 2𝜇)𝑘22 𝐽0 (𝑘2 𝑟) 𝑟 ( ) 𝐽1 (𝑘3 𝑟) + 𝑎3 2𝜇𝑘3 − (𝜆 + 2𝜇)𝑘23 𝐽0 (𝑘3 𝑟) . (E.4) 𝑟

Appendix C. Wave numbers

𝐶2 =

𝑢𝑟 = −𝑎 1 𝑘1 𝐽1 (𝑘1 𝑟) − 𝑎 2 𝑘2 𝐽1 (𝑘2 𝑟) − 𝑎 3 𝑘3 𝐽1 (𝑘3 𝑟),

Θ12 𝑟22 𝐽0 (𝑘2 𝑅) 𝑟32 𝐽0 (𝑘2 𝑅)

⎤⎡𝑎1 ⎤ ⎡−1⎤ Θ13 ⎥⎢ ⎥ ⎢ ⎥ 𝑟23 𝐽0 (𝑘3 𝑅)⎥⎢𝑎2 ⎥ = ⎢ 0 ⎥, ⎥⎢ ⎥ ⎢ ⎥ 𝑟33 𝐽0 (𝑘3 𝑅)⎦⎣𝑎3 ⎦ ⎣ 0 ⎦

𝐽 (𝑘 𝑅)

(E.5)

with Θ1𝑖 = 2𝜇𝑘𝑖 1 𝑅𝑖 − (𝜆 + 2𝜇)𝑘2𝑖 𝐽0 (𝑘𝑖 𝑅). II) For the ring-shaped problem with the boundary conditions: 𝜎𝑟𝑠 = 0N∕m2 , pw = pa = 0Pa on the inner boundary, and 𝜎𝑟𝑠 = −1N∕m2 pw = pa = 0Pa on outer boundary, the exact solutions can be given by 𝑢𝑟 = − 𝑎 1 𝑘1 𝐽1 (𝑘1 𝑟) − 𝑎 2 𝑘2 𝐽1 (𝑘2 𝑟) − 𝑎 3 𝑘3 𝐽1 (𝑘3 𝑟) − 𝑏1 𝑘1 𝑌1 (𝑘1 𝑟)

Appendix D. Geomechanical characteristics of Massilon sandstone (from [8,20])

− 𝑏2 𝑘2 𝑌1 (𝑘2 𝑟) − 𝑏3 𝑘3 𝑌1 (𝑘3 𝑟),

The parameters of the considered material are:

𝑝𝑤 = 𝑎1 𝑟21 𝐽0 (𝑘1 𝑟) + 𝑎2 𝑟22 𝐽0 (𝑘2 𝑟) + 𝑎3 𝑟23 𝐽0 (𝑘3 𝑟) + 𝑏1 𝑟21 𝑌0 (𝑘1 𝑟) + 𝑏2 𝑟22 𝑌0 (𝑘2 𝑟) + 𝑏3 𝑟23 𝑌0 (𝑘3 𝑟),

Pore size distribution index: ϑ =1.5, Saturation of water: Sw = 0.9, Residual wetting fluid saturation: Srw = 0, Non-wetting fluid entry saturation: Sra = 1, Porosity: n = 0.23, Density of the solid skeleton: 𝜌s = 2650 kg/m3 , Density of the water: 𝜌w = 997 kg/m3 , Density of the air: 𝜌a = 1.10 kg/m3 , Drained bulk modulus of the mixture: K = 1.02 × 109 N/m2 , Shear modulus: G = 1.44 × 109 N/m2 , Bulk modulus of the solid grains: Ks = 3.5 × 1010 N/m2 , Bulk modulus of the water: Kw = 2.25 × 109 N/m2 , Bulk modulus of the air: Ka = 1.10 × 105 N/m2 , Intrinsic permeability:𝑘̄ = 2.5 × 10−12 m2 , Viscosity of the water: 𝜂 w = 1.0 × 10−3 N • s/m2 , Viscosity of the air:𝜂 a = 1.8 × 10−5 N • s/m2 , Gas entry pressure:pd = 5 × 104 N/m2 .

𝑝𝑎 = 𝑎1 𝑟31 𝐽0 (𝑘1 𝑟) + 𝑎2 𝑟32 𝐽0 (𝑘2 𝑟) + 𝑎3 𝑟33 𝐽0 (𝑘3 𝑟) + 𝑏1 𝑟31 𝑌0 (𝑘1 𝑟) + 𝑏2 𝑟32 𝑌0 (𝑘2 𝑟) + 𝑏3 𝑟33 𝑌0 (𝑘3 𝑟), ( ) 𝐽 (𝑘 𝑟) 𝜎𝑟 = 𝑎1 2𝜇𝑘1 1 1 − (𝜆 + 2𝜇)𝑘21 𝐽0 (𝑘1 𝑟) 𝑟 ( ) 𝐽1 (𝑘2 𝑟) + 𝑎2 2𝜇𝑘2 − (𝜆 + 2𝜇)𝑘22 𝐽0 (𝑘2 𝑟) 𝑟 ( ) 𝐽1 (𝑘3 𝑟) + 𝑎3 2𝜇𝑘3 − (𝜆 + 2𝜇)𝑘23 𝐽0 (𝑘3 𝑟) 𝑟 ( ) 𝑌 (𝑘 𝑟) + 𝑏1 2𝜇𝑘1 1 1 − (𝜆 + 2𝜇)𝑘21 𝑌0 (𝑘1 𝑟) 𝑟 ( ) 𝑌1 (𝑘2 𝑟) + 𝑏2 2𝜇𝑘2 − (𝜆 + 2𝜇)𝑘22 𝑌0 (𝑘2 𝑟) 𝑟 97

L. Sun, X. Wei and B. Chen

Engineering Analysis with Boundary Elements 113 (2020) 82–98

( ) 𝑌 (𝑘 𝑟) + 𝑏3 2𝜇𝑘3 1 3 − (𝜆 + 2𝜇)𝑘23 𝑌0 (𝑘3 𝑟) , 𝑟

[25] Nennig B, Perrey-Debain E, Chazot JD. The method of fundamental solutions for acoustic wave scattering by a single and a periodic array of poroelastic scatterers. Eng Anal Boundary Elem 2011;35(8):1019–28. [26] Liu C-S. An equilibrated method of fundamental solutions to choose the best source points for the Laplace equation. Eng Anal Boundary Elem 2012;36(8):1235–45. [27] Chen CS, Karageorghis A, Li Y. On choosing the location of the sources in the MFS. Numer Algor 2016;72(1):107–30. [28] Wang F, Liu C-S, Qu W. Optimal sources in the MFS by minimizing a new merit function: energy gap functional. Appl Math Lett 2018;86:229–35. [29] Li J, Chen W. A modified singular boundary method for three-dimensional high frequency acoustic wave problems. Appl Math Modell 2018;54:189–201. [30] Chen W, Fu Z, Wei X. Potential problems by singular boundary method satisfying moment condition. Comput Model Eng Sci (CMES) 2009;54(1):65–85. [31] Li W, Chen W, Fu Z. Precorrected-FFT Accelerated Singular Boundary Method for Large-Scale Three-Dimensional Potential Problems. Commun Comput Phys 2017;22(2):460–72. [32] Wei X, Chen W, Sun L, Chen B. A simple accurate formula evaluating origin intensity factor in singular boundary method for two-dimensional potential problems with Dirichlet boundary. Eng Anal Boundary Elem 2015;58:151–65. [33] Wei X, Sun L, Yin S, Chen B. A boundary-only treatment by singular boundary method for two-dimensional inhomogeneous problems. Appl Math Modell 2018;62:338–51. [34] Wei X, Chen W, Chen B, Sun L. Singular boundary method for heat conduction problems with certain spatially varying conductivity. Comput Math Appl 2015;69(3):206–22. [35] Gu Y, He X, Chen W, Zhang C. Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method. Comput Math Appl 2018;75(1):33–44. [36] Wei X, Huang A, Sun L. Singular boundary method for 2D and 3D heat source reconstruction. Appl Math Lett 2020;102:106103. [37] Lin J, Zhang C, Sun L, Lu J. Simulation of seismic wave scattering by embedded cavities in an elastic half-plane using the novel singular boundary method. Adv Appl Math Mech 2018;10(2):322–42. [38] Tang ZC, Fu ZJ, J ZD, Huang JD. Singular boundary method to simulate scattering of SH wave by the canyon topography. Adv Appl Math Mech 2018;10(4):912–24. [39] Qu W, Chen W, Gu Y. Fast multipole accelerated singular boundary method for the 3D Helmholtz equation in low frequency regime. Comput Math Appl 2015;70(4):679–90. [40] Fu Z, Chen W, Wen P, Zhang C. Singular boundary method for wave propagation analysis in periodic structures. J Sound Vib 2018;425:170–88. [41] Li J, Qin Q, Fu Z. A dual-level method of fundamental solutions for three-dimensional exterior high frequency acoustic problems. Appl Math Modell 2018;63:558–76. [42] Li W. A fast singular boundary method for 3D Helmholtz equation. Comput Math Appl 2018. [43] Li J, Fu Z, Chen W. Numerical investigation on the obliquely incident water wave passing through the submerged breakwater by singular boundary method. Comput Math Appl 2016;71(1):381–90. [44] Sun L, Chen W, Cheng AHD. Singular boundary method for 2D dynamic poroelastic problems. Wave Motion 2016;61:40–62. [45] Sun L, Wei X. A frequency domain formulation of the singular boundary method for dynamic analysis of thin elastic plate. Eng Anal Boundary Elem 2019;98:77–87. [46] Li W, Chen W. Simulation of the band structure for scalar waves in 2D phononic crystals by the singular boundary method. Eng Anal Boundary Elem 2019;101:17–26. [47] Li W, Chen W. Band gap calculations of photonic crystals by singular boundary method. J Comput Appl Math 2017;315:273–86. [48] Wei X, Sun L. Singular boundary method for 3D time-harmonic electromagnetic scattering problems. Appl Math Modell 2019;76:617–31. [49] Wei X, Chen W, Fu Z. Solving inhomogeneous problems by singular boundary method. J Mar Sci Technol 2013;21(1):8–14. [50] Li X, Zhu J. The method of fundamental solutions for nonlinear elliptic problems. Eng Anal Boundary Elem 2009;33(3):322–9. [51] Li W. A fast singular boundary method for 3D Helmholtz equation. Comput Math Appl 2019;77(2):525–35. [52] Wei X, Chen B, Chen S, Yin S. An ACA-SBM for some 2D steady-state heat conduction problems. Eng Anal Boundary Elem 2016;71:101–11. [53] Kausel E, Roësset JM. Frequency Domain Analysis of Undamped Systems. J Eng Mech 1992;118(4):721–34. [54] Maier G. The finite element method in the static and dynamic deformation and consolidation of porous media - R.W. Lewis and B.A. Schrefler. Meccanica 1999;34(3):231–2. [55] Tanneau O, Lamary P, Chevalier Y. A boundary element method for porous media. J Acoust Soc Am 2006;120(3):1239–51. [56] Gu Y, Chen W, He X-Q. Singular boundary method for steady-state heat conduction in three dimensional general anisotropic media. Int J Heat Mass Transfer 2012;55(17–18):4837–48. [57] Sun L, Chen W, Cheng AHD. Method of fundamental solutions without fictitious boundary for plane time harmonic linear elastic and viscoelastic wave problems. Comput Struct 2016;162:80–90. [58] Phan AV, Gray LJ, Salvadori A. Transient analysis of the dynamic stress intensity factors using SGBEM for frequency-domain elastodynamics. Comput Meth Appl Mech Eng 2010;199(45):3039–50. [59] Ariza MP, Dominguez J. General BE approach for three-dimensional dynamic fracture analysis. Eng Anal Boundary Elem 2002;26(8):639–51. [60] Xiao J, Ye W, Wen L. Efficiency improvement of the frequency-domain BEM for rapid transient elastodynamic analysis. Comput Mech 2013;52(4):903–12.

where Yk denotes the Bessel function of the second kind of order k, a1 , a2 , a3 , b1 , b2 and b3 are unknown coefficients which can be determined through enforcing the above formulas to satisfy the boundary conditions similar to the above circle case. III) For the annular-sector domain problem with the following boundary conditions: 𝑎 𝑢𝑠𝑛 = 0m,𝑡𝜏𝑠̃ = 0N∕m2 , 𝑢𝑤 𝑛 = 𝑢𝑛 = 0m on Γ1 , Γ3 ,

𝑡𝑠𝑛 = −1N∕m2 ,𝑡𝜏𝑠̃ = 0N∕m2 , 𝑝𝑤 = 𝑝𝑎 = 0Pa on Γ2 , 𝑡𝑠𝑛 = 𝑡𝜏𝑠̃ = 0N∕m2 , 𝑝𝑤 = 𝑝𝑎 = 0Pa on Γ4 , where Γ1 ,Γ2 ,Γ3 and Γ4 are shown in Fig. 10, the forms of exact solution are the same as those of the ring-shaped problems. References [1] Biot MA. Thermoelasticity and Irreversible Thermodynamics. J Appl Phys 1956;27(3):240–53. [2] Schanz M. Poroelastodynamics: Linear Models, Analytical Solutions, and Numerical Methods. Appl Mech Rev 2009;62(3) 030803-03-15. [3] Deckers E, Hörlin N-E, Vandepitte D, Desmet W. A Wave Based Method for the efficient solution of the 2D poroelastic Biot equations. Comput Meth Appl Mech Eng 2012;201-204:245–62. [4] Chen J. Time domain fundamental solution to Biot’s complete equations of dynamic poroelasticity. Part I: Two-dimensional solution. Int J Solids Struct 1994;31(10):1447–90. [5] HÖRlin NE, NordstrÖM M, GÖRansson P. A 3-D hierarchical Fe formulation of biot’s equations for elasto-acoustic modelling of porous media. J Sound Vib 2001;245(4):633–52. [6] Sheng D. Review of fundamental principles in modelling unsaturated soil behaviour. Comput Geotech 2011;38(6):757–76. [7] Shan ZD, Ling DS, Ding HJ. Exact solutions for the one-dimensional transient response of unsaturated single-layer porous media. Comput Geotech 2013;48:127–33. [8] Li P, Schanz M. Wave propagation in a 1-D partially saturated poroelastic column. Geophys J Int 2011;184(3):1341–53. [9] Tan T-S, Phoon K-K, Chong P-C. Numerical Study of Finite Element Method Based Solutions for Propagation of Wetting Fronts in Unsaturated Soil. J Geotech Geoenviron Eng 2004;130(3):254–63. [10] Sheng D, Sloan SW, Gens A, Smith DW. Finite element formulation and algorithms for unsaturated soils. Part I: Theory. Int J Numer Anal Methods Geomech 2003;27(9):745–65. [11] Callari C, Abati A. Finite element methods for unsaturated porous solids and their application to dam engineering problems. Comput Struct 2009;87(7):485–501. [12] Dreyer D, Petersen S, von Estorff O. Effectiveness and robustness of improved infinite elements for exterior acoustics. Comput Meth Appl Mech Eng 2006;195(29):3591–607. [13] Liu GR, Quek Jerry SS. A non-reflecting boundary for analyzing wave propagation using the finite element method. Finite Elem Anal Des 2003;39(5):403–17. [14] Sun Y. Indirect Boundary Integral Equation Method for the Cauchy Problem of the Laplace Equation. J Sci Comput 2017;71(2):469–98. [15] Messner M. Fast boundary element methods in acoustics. Graz University of Technology; 2011. [16] Gaul L, Schanz M. A comparative study of three boundary element approaches to calculate the transient response of viscoelastic solids with unbounded domains. Comput Meth Appl Mech Eng 1999;179(1):111–23. [17] Maghoul P, Gatmiri B, Duhamel D. Boundary integral formulation and two-dimensional fundamental solutions for dynamic behavior analysis of unsaturated soils. Soil Dyn Earthquake Eng 2011;31(11):1480–95. [18] Gatmiri B, Jabbari E. Time-domain Green’s functions for unsaturated soils. Part I: two-dimensional solution. Int J Solids Struct 2005;42(23):5971–90. [19] Gatmiri B, Jabbari E. Time-domain Green’s functions for unsaturated soils. Part II: three-dimensional solution. Int J Solids Struct 2005;42(23):5991–6002. [20] Li P, Schanz M. Time domain boundary element formulation for partially saturated poroelasticity. Eng Anal Boundary Elem 2013;37(11):1483–98. [21] Karageorghis A, Fairweather G. The method of fundamental solutions for axisymmetric elasticity problems. Comput Mech 2000;25(6):524–32. [22] Fan CM, Chen CS, Monroe J. The method of fundamental solutions for solving convection-diffusion equations with variable coefficients. Adv Appl Math Mech 2009;1(2):215–30. [23] Wang F, Fan C-M, Hua Q, Gu Y. Localized MFS for the inverse Cauchy problems of two-dimensional Laplace and biharmonic equations. Appl Math Comput 2020;364:124658. [24] Qu W, Fan C-M, Gu Y, Wang F. Analysis of three-dimensional interior acoustic fields by using the localized method of fundamental solutions. Appl Math Modell 2019;76:122–32. 98