Scattering problem for a spherical inclusion in poroelastic media: Application of the asymptotic expansion method

International Journal of Engineering Science 128 (2018) 187–207

Contents lists available at ScienceDirect

International Journal of Engineering Science journal homepage: www.elsevier.com/locate/ijengsci

Scattering problem for a spherical inclusion in poroelastic media: Application of the asymptotic expansion method S. Kanaun, V. Levin, M. Markov∗ Instituto Mexicano del Petróleo, Eje Central Lázaro Cárdenas 152, CDMX, D.F.07730, México

a r t i c l e

i n f o

Article history: Received 8 January 2018 Revised 26 February 2018 Accepted 13 March 2018

a b s t r a c t Scattering of longitudinal monochromatic waves on an isolated inclusion in an infinite poroelastic medium is considered. Wave propagation in both medium and inclusion is described by Biot’s equations of poroelasticity. For the material parameters typical for sedimentary rocks, the corresponding system of equations contains a second order differential operator with a small parameter. As the result, the wave field in the medium consists of a slowly changing part and boundary layer functions concentrated near the inclusion interface. Such a specific form of the wave field should be taken into account by the numerical solution of the problem since application of conventional numerical methods can result in substantial errors. In the present paper, the method of matched asymptotic expansions is applied to the solution of the scattering problem for a spherical inclusion. The systems of equations for a slowly changing part of the wave field and for boundary layer functions are derived. The solutions of these systems for the fields in the vicinity of the inclusion and for far wave fields are obtained and analyzed. Comparisons of the results of the straightforward solution and the asymptotic expansion are presented. © 2018 Elsevier Ltd. All rights reserved.

1. Introduction The theory of wave propagation in fluid-saturated porous media was developed by Biot (1956, 1962) (see also Bourbie, Coussy, & Zinzner, 1987). An important area of application of this theory is wave propagation in geologic structures. For many actual structures, the system of equations of poroelasticity contains a differential operator with a small parameter (Section 2). As a result, in vicinities of the interfaces of scatterers (inclusions) there appear narrow boundary layers where wave fields change fastly. Apparently, the existence of such layers was firstly indicated in Mei and Foda (1981). Such a specific structure of the wave fields requires appropriate methods for numerical solution of the scattering problems, since conventional methods can result in substantial errors. In the current work, the problem of scattering of a longitudinal monochromatic wave on a spherical inclusion is considered. The existence of a differential operator with a small parameter in the equations of poroelasticity is taken into account. The theory of partial differential equations with a small parameter in the highest order differential operators was developed in 70th and 80th of the last century by many authors (see, e.g., Il’yin, 1992 and references therein). The principal strategy of solution of such equations is based on the so-called formal asymptotic expansion with respect to the small parameter and

∗

Corresponding author. E-mail address: [email protected] (M. Markov).

https://doi.org/10.1016/j.ijengsci.2018.03.003 0020-7225/© 2018 Elsevier Ltd. All rights reserved.

188

S. Kanaun et al. / International Journal of Engineering Science 128 (2018) 187–207

then, constructing the boundary layer functions at the vicinity of the region boundary. The sum of the formal asymptotics and boundary layer functions provides the complete solution. For homogeneous spherical inclusions, the scattering problem of poroelasticity can be solved by using expansion of the wave fields in series of spherical harmonics (Berryman, 1985; Ciz, Gurevitch, & Markov, 2006; Krutin, Markov, & Yumatov, 1984; Liu, Greenhalgh, & Zhou, 2009; Morochnic & Barded, 1996; Zimmerman & Stern, 1993). Formally, this method can be applied to the solution of the scattering problem in the whole range of the parameters of poroelastic media. However, if the above mentioned parameter is small, on the one hand, the conditional numbers of the matrices of the system for the coefficients of these series turn out to be huge, and a reliable solution is not secured. On the other hand, the number of terms in the series should be taken very large in order to describe fast changing fields near the inclusion interface with sufficient accuracy. The structure of the paper is as follows. In Section 2, the system of equations of poroelasticity is revised and its coefficients for actual geological structures are assessed. In Section 3, a straightforward serial solution of the scattering problem for a spherical inclusion is presented. In Section 4, the formal asymptotic expansion of the solution with respect to a small parameter and the corresponding boundary layer functions are constructed. The zero-order term of the asymptotic solution is obtained in Section 5. Numerical examples of the wave field in the vicinity of the inclusion are presented in Section 6. In Section 7, the wave fields far from the inclusion are considered. An important characteristic of the energy flux of the scattered field (the inclusion differential scattering cross-section) is obtained and analyzed in Section 8. 2. Propagation of plane monochromatic waves in poroelastic media For an isotropic homogeneous medium, equations of dynamic poroelasticity in terms of displacement vector U of the solid skeleton and fluid pressure P have the form (Biot, 1956; 1962)

 P = 0, (λ + μ )∇ div U + μU − ρt U¨ − α∇

(1)

1 ¨ − P + β P¨ = 0. div U α ρ

(2)

In these equations, λ and μ are effective Lame constants of the solid skeleton filled with fluid,

ρt = φρ f + (1 − φ )ρs −

ρ 2f ρ =α− f, , α ρ ρ

(3)

ρ f and ρ s are the densities of the fluid and solid phases; φ is the medium porosity; α and β are Biot’s parameters

α =1−

K , Ks

β=

α−φ Ks

+

φ Kf

,

(4)

 where Ks and Kf are the bulk moduli of the solid and fluid phases; K is the effective bulk modulus of the solid skeleton; ρ is the operator related to viscoelastic properties of the fluid and permeability of the medium. In Eq. (1), Ü and P¨ are the second time derivatives of the displacement vector and pressure. For monochromatic waves of frequency ω, U(y, t ) = u(y )eiωt and P (y, t ) = p(y )eiωt , where y(y1 , y2 , y3 ) is a point of the medium, and Eqs. (1) and (2) can be rewritten in terms of amplitudes u(y) and p(y) as follows

 p = 0, (λ + μ )∇ div u + μu + ρt ω2 u − α∇  1  div u + α  + β p = 0.  2 ρω

(5) (6)

 is a parameter of the dimension of density (Pride & Haartsen, 1996) In Eq. (6), ρ

ρ = −

iη , ωχ

(7)

where η is the fluid viscosity (Pa · sec); χ is the permeability of the medium (m2 ). Further, the scattering problem of monochromatic longitudinal waves in a poroelastic medium containing an inclusion with characteristic size a is considered. In this case, one can introduce dimensionless coordinates xi = yi /a and functions (displacements and pressure)

 u (x ) =

u (x ) p( x ) ,  p( x ) = , a μ

(8)

and rewrite system (5) and (6) in the form

(λ + μ )   u +  u + qt2 u − α∇ p = 0, ∇ div  μ   2 div  α u + iε  + μβ  p = 0.

(9) (10)

S. Kanaun et al. / International Journal of Engineering Science 128 (2018) 187–207

189

Here, the gradient ∇ , divergence div, and Laplacian  are calculated with respect to x-variables,

qt2 =

iμ ω2 ρt 2 2 μχ a , ε =− = .  2 a2 μ ρω ω η a2

(11)

For typical geological structures (see, e.g., Mavko, Mukerji, & Dvorkin, 2009), the following assessments of the coefficients in Eqs. (9) and (10) hold

λ+μ  = O(1 ), μβ ∼ O(1 ), ∼ O ( 1 ), α μ

(12)

μ ∼ 1010 Pa, ρt ∼ 103 kg/m3 , χ ∼ 10−13 ÷ 10−14 m2 , η ∼ 10−2 ÷ 10−4 Pa · s. As a result, for inclusions with characteristic size   seismic experiments ω ∼ 103 s−1 , we obtain

a ∼ 10÷102 m

(13)

and for typical frequencies of propagating waves used in

ε 2 ∼ 10−3 ÷ 10−4 .

(14)

Thus, the parameter ε 2 is smaller than other coefficients of the system (9) and (10), and it decreases as the wave frequency increases or the medium permeability decreases. Monochromatic longitudinal plane waves with the amplitudes of displacements U0 and the pressure P0 propagating in poroelastic media have the form

U(x, t ) = U0 m exp(i(ωt − qm · x )),

P (x, t ) = P0 exp(i(ωt − qm · x )),

(15)

where q is the wave number; qm is the wave vector; m · x is scalar product of the vectors m and x. The characteristic equation for the dimensionless wave number q of longitudinal waves follows from the system (9) and (10)

ε 2 q4 + i

μ M



2 q2 − iqt2 ε 2 + Mβ + α

iμqt2 (μβ ) = 0, M

(16)

M = λ + 2 μ.

(17)

Two roots of this equation with negative imaginary parts correspond to two types of longitudinal waves (f- and s-waves) that can propagate in the medium. The first terms of the expansion of these roots with respect to the parameter ε have the form



q f = qt

μβ

(1 − i ) + O ( ε 2 ), qs = 2 + Mβ ε α

   2 + Mβ μ α 2M

+ O ( ε ).

(18)

The corresponding velocities vf , vs and the attenuation factors γ f , γ s of these waves are assessed as follows

vf =

ωa = O ( 1 ), Re q f

vs =

ωa = O ( ε ), Re qs

γ f a = Im q f = O(1 ), γs a = Im qs = O(ε −1 ).

(19) (20)

Thus, the velocity vs is much smaller than vf , and attenuation factor γ s is much higher than γ f . Note that to applicate the equation of poroelasticity, the frequency of the propagating waves should not exceed the socalled Biot frequency ωc

ω < ωc ,

ωc =

ηφ , ρf χ

(21)

and for mentioned parameters of poroelastic media, ωc ∼ 107 s−1 . 3. Scattering of plane longitudinal waves on a spherical inclusion in poroelastic media Let us consider a homogenous poroelastic medium with a spherical inclusion of radius a = 1 in dimensionless coordinates xi . The poroelastic parameters of the inclusion are indicated by the index 1

1 , β1 , ρt1 , ε12 , λ1 , μ1 , α

(22)

and the parameters of the medium – by the index 2

2 , β2 , ρt2 , ε22 . λ2 , μ2 , α

(23)

We introduce Cartesian (x1 , x2 , x3 ) and spherical (r, θ , ϕ ) coordinate systems with the origins at the center of the inclusion, r = |x|, the polar axis of the spherical system coincides with the x3 -axis (Fig. 1).

190

S. Kanaun et al. / International Journal of Engineering Science 128 (2018) 187–207

Fig. 1. Spherical inclusion in a poroelastic medium subjected to a longitudinal incident wave field directed along the x3 -axis.

The incident (fast) longitudinal wave of displacements and pressure propagates along the x3 -axis and it is scattered on the inclusion. In the Cartesian system, the incident displacement and pressure waves have the form

u(i ) = U0 e3 exp(−iq f 2 x3 ), p(i ) = −

2 μ2 iq f 2 α U exp(−iq f 2 x3 ), γ2 − μ2 β2 0

γ2 =

ρt2 q2f 2 · 2 . ρ2 qt2

(24)

On the inclusion interface (r = 1 ) the following boundary conditions hold (Derisiewicz & Skalak, 1963)

ur1 = ur2 ,

uθ 1 = uθ 2 ,

(25)

σr r 1 = σr r 2 , σr θ 1 = σr θ 2 ,

(26)

p1 = p2 ,

(27)

wr1 = wr2 .

Here, ur , uθ are the components of the displacement vector in the spherical coordinate system; σ rr , σ rθ are the components of the stress tensor

σ = (λdiv u − α p)I + μ(∇  u + u  ∇ ), Ii j = δi j ,

(28)

and wr is the radial component of the vector w that is the relative displacement of the fluid with respect to the solid skeleton (Biot, 1962)

w=

1

ω2 ρ

∇p−

ρf u. ρ

(29)

Using expansion of the field in the series of spherical harmonics, the radial and tangential components of the displacement field in the inclusion are presented in the form (Morochnic & Barded, 1996; Zimmerman & Stern, 1993)

ur1 = uθ 1 =

∞  1

Bn f f11n (q f 1 r ) + Bns f11n (q2s1 r ) + Bnt f21n (qt1 r ) Pn (cos θ ), r n=0 ∞

1 r



Bn f f31n (q f 1 r ) + Bns f31n (q2s1 r ) + Bnt f41n (qt1 r )

n=0

(30)

 dPn (cos θ ) . dθ

(31)

The pressure p1 (x) in the inclusion takes the form

p1 = −

∞

μ1

r2



Bn f f91n (q f 1 r, s f 1 ) + Bns f91n (qs1 r, ss1 ) Pn (cos θ ), s f 1 =

n=0

qt1 qt1 , ss1 = . qf1 qs1

(32)

The components of the displacement vector in the medium are presented as follows

ur2 = −

∞ U0 (−i )n+1 (2n + 1 ) f11n (q f 2 r )Pn (cos θ ) q f 2r n=0

∞  1

+ Cn f f12n (q f 2 r ) + Cns f12n (qs2 r ) + Cnt f22n (qt2 r ) Pn (cos θ ), r n=0

uθ 2 = −

∞ U0 dP (cos θ ) (−i )n+1 (2n + 1 ) f31n (q f 2 r ) n q f 2r dθ n=0

(33)

S. Kanaun et al. / International Journal of Engineering Science 128 (2018) 187–207

+

∞  dPn (cos θ ) 1

Cn f f32n (q f 2 r ) + Cns f32n (qs2 r ) + Cnt f42n (qt2 r ) . r dθ

191

(34)

n=0

and the pressure in the host medium is

p2 =

μ2 q f 2r −

U 2 0

∞

(−i )n+1 (2n + 1 ) f9n (q f 2 r, s f 2 )Pn (cos θ )

n=0

∞

μ2

r2



Cn f f92n (q f 2 r, s f 2 ) + Cns f92n (qs2 r, ss2 ) Pn (cos θ ),

(35)

n=0

i Here, Pn (t) are Legendre polynomials; the functions fmn are presented in Appendix as well as the series for stresses σ rr , σ rθ and displacement wr . The first sums in Eqs. (33)–(35) are the expansions of the incident fields (24) in the series of spherical harmonics; the second sums are the fields scattered on the inclusion. The constants Cnf , Cns , Cnt and Bnf , Bns , Bnt are to be found from conditions (25)–(27). The corresponding system can be written in the canonical form as follows

M ( n )X (n ) = F ( n ),

⎡

M(n )= M2 (n )  M1 (n ), n = 0, 1, 2, . . .

f12n (q f 2 )

f12n (qs2 )

f22n (kt2 )

⎤

⎢ f 2 (q f 2 ) f32n (qs2 ) f42n (qt2 )⎥ ⎢ 3n ⎥ ⎢ f 2 (q , s ) ⎥ 2 2 f ( q , s ) f ( q ) s 2 s 2 t2 ⎢ ⎥ f2 f2 5n 6n M 2 ( n ) = ⎢ 5n 2 ⎥, f72n (qs2 ) f82n (qt2 )⎥ ⎢ f 7n ( q f 2 ) ⎢ 2 ⎥ ⎣ f 9n ( q f 2 , s f 2 ) f92n (qs2 , ss2 ) 0 ⎦ 2 2 f10 f10 0 n (q f 2 , s f 2 ) n ( qs2 , ss2 ) ⎡ ⎤ − f11n (q f 1 ) − f11n (qs1 ) − f21n (qt1 ) ⎢ − f 1 (q f 1 ) − f31n (qs1 ) − f41n (qt1 ) ⎥ 3n ⎢ ⎥ ⎢−ξ f 1 (q , s ) −ξ f 1 (q , s ) −ξ f 1 (q )⎥ s1 s1 ⎢ f1 f1 f1 ⎥ 5n 5n 6n M1 ( n ) = ⎢ ⎥, −ξ f71n (qs1 ) −ξ f81n (qt1 ) ⎥ ⎢ −ξ f71n (q f 1 ) ⎢ ⎥ ⎣−ξ f91n (q f 1 , s f 1 ) −ξ f91n (qs1 , ss1 ) ⎦ 0 1 1 − f10n (q f 1 , s f 1 ) − f10n (qs1 , ss1 ) 0 where

ξ=

μ1 qti qti ,s = ,s = , i = 1, 2; μ2 f i q f i si qsi ⎡C ⎤ nf

⎡

f11n (q f 2 )

f2

(37)

(38)

⎤

⎢ f 1 (q f 2 ) ⎥ ⎢ Cns ⎥ ⎢ 3n ⎥ ⎢ ⎥ ⎢ f 1 (q , s ) ⎥ ⎢ Cnt ⎥ 2 n + 1 ⎢ ⎥ f 2 f 2 n+1 ⎥ X (n ) = ⎢ ⎢ 5n 1 ⎥. ⎢Bn f ⎥, F (n ) = U0 (−i ) q f ( q ) ⎢ ⎥ f2 f2 7n ⎢ ⎥ ⎢ ⎥ ⎣ Bns ⎦ 1 ⎣ f (q , s ) ⎦ 9n

(36)

(39)

f2

1 f10 n (q f 2 , s f 2 )

Bnt

For the numerical solution, the number N of terms in series (30)–(35) should be taken finite, so that sufficient accuracy can be achieved. An appropriate choice of N was discussed in the book of Bohren and Huffman (1983), where the electromagnetic scattering problem for a spherical inclusion was considered. The authors proposed that number N should be taken √ equal to the integer nearest to q + 4 3 q + 2, where q is a dimensionless wave number. In the case of poroelasticity, it is natural to take q equal to maximal of the absolute values of the wave numbers qf , qt , and qs . Let us consider an example of a real geological structure. The inclusion material consists of a solid skeleton with Young’s modulus Es1 = 10 GPa and Poisson’s ratio νs1 = 0.2. The skeleton porosity φ1 = 0.3; the bulk modulus of the fluid in the pores K f 1 = 1.5 GPa; the viscosity η1 = 10−3 Pa · s; the densities of the solid and fluid phases ρs1 = 2.5 · 103 kg/m3 and ρ f 1 = 1.2 · 103 kg/m3 , respectively; the permeability of the inclusion material χ1 = 10−13 m2 . The host medium has skeleton elastic parameters: Young’s modulus Es2 = 50 GPa; Poisson’s ratio νs2 = 0.2; the porosity φ2 = 0.2. The fluid bulk modulus K f 2 = 2 GPa ; the viscosity η2 = 10−3 Pa · s; the densities of the solid and fluid phases ρs2 = 2.65 · 103 kg/m3 and ρ f 2 = 103 kg/m3 , respectively; the host medium permeability χ2 = 10−14 m2 . These parameters correspond to a typical sedimentary rocks at a depth of 1.5÷2.0 km (Mavko et al., 2009). The effective bulk K and shear μ moduli of the host medium and inclusion materials are calculated by the equations of the effective field method (Kanaun & Levin, 2008)



K = Ks 1 +

 φ (K f − Ks ) , Ks + (1 − φ )S1 (K f − Ks )

(40)

192

S. Kanaun et al. / International Journal of Engineering Science 128 (2018) 187–207

 μ = μs 1 − S1 =

φ

 ,

1 − ( 1 − φ )S2

(41)

3Ks 1 , S2 = ( 3 − S1 ). 3Ks + 4μs 5

For the inclusion of radius a = 1 m, we have the following values of the dimensionless wave numbers, and the number of terms N in series (30)–(35) The imaginary parts of the wave numbers qf2 and qt2 is of the order of 10−7 . A great number of terms is required for accurate descriptions of fast changing wave fields with wave number qs2 . Another problem relates to the solution of linear algebraic Eqs. (39). Conditional numbers d(n) of the be defined as the ratio of the maximal and minimal singular values ek (n) of these matrices (Golub & Ch,

d (n ) =

maxk ek (n ) , k = 1, 2, . . . , 6. mink ek (n )

(42) the parameter ε22 , in series (30)–(35) matrices M(n) can 1989)

(43)

In Table 2, the wave numbers of the transverse waves qt2 = 0.1, 0.5, 1, 2, and the decimal logarithms of the conditional numbers of the matrices M(n) for n = 0, 1, 2, . . . , 6 are presented. Huge magnitudes of the conditional numbers of the matrices M(n) show that sufficient accuracy of the calculation of the coefficients in series (30)–(35) cannot be achieved. In the next Table 3, the values of the constants Cnf , Cns , Cnt and Bnf , Bns , Bnt in Eqs. (30)–(35) for the wave number qt2 = 0.1 and n = 0, 1, 2 are presented. Great values of Cns and small values of Bns are the results of solution of ill-posed system (36) that does not give reliable values of these constants. 4. Asymptotic expansion of the solution and boundary layer functions The mentioned difficulties in numerical solution of the scattering problem are typical for partial differential equations with small parameters in the differential operator of the highest order. A stable algorithm of solution of such problems is described in Il’yin (1992), and it consists of the following steps. First, the solution is presented in series with respect to small parameter ε for the fields inside and outside the inclusion. In the case of poroelasticiy, these series are

u ( x ) = u (0 ) ( x ) + ε 2 u (2 ) ( x ) + . . . ,

(44)

p(x ) = p(0) (x ) + ε 2 p(2) (x ) + . . . .

(45)

Substituting these expansions into Eqs. (9) and (10) and joining the terms of the same order with respect to ε , we obtain the following system of differential equations for functions u(m) (x), p(m) (x)

 p(0) = 0, Lu(0) + qt2 u(0) − α∇

div u(0) + μβ p(0) = 0 α

(46)

 p(2) = 0, Lu(2) + qt2 u(2) − α∇

div u(2) + μβ p(2) = −i p(0) , α

(47)

where L is the differential operator of the second order

L=

(λ + μ ) ∇ div + . μ

(48)

Substituting p(0) from the second equation of the system (46) for u(0) , p(0) into the first one, we yield



 λ+μ ∇ divu(0) + u(0) + qt2 u(0) = 0, μ

p(0) = −

 2 α α div u(0) , λ = λ + . μβ β

(49) (50)

Series (44) and (45) are called the formal asymptotic expansions of the solution of the original problem (9) and (10). If the scattering problem is considered, the “zero”-terms u(0) , p(0) allow one to satisfy only the boundary conditions (25) and (26). In order to satisfy conditions (27), the so-called boundary layer functions should be introduced. These functions are concentrated near the inclusion interface  and rapidly vanish in the normal direction to . Following the scheme or construction of these functions pointed out in Il’yin (1992), we introduce a local Cartesian coordinate system (z1 , z2 , z3 ) with the origin at the point x0 of the inclusion interface and the z3 -axis directed along the normal to  at this point (Fig. 2). Then, we define extended coordinates ςi = zi /ε , and the complete solution of the problem is presented in the form of sums of the formal asymptotic expansions and boundary layers functions

u (x ) =

∞ k=0

ε 2k u ( 2k ) ( x ) +

∞ m=0

ε m V ( m ) ( x0 , ς ),

(51)

S. Kanaun et al. / International Journal of Engineering Science 128 (2018) 187–207

193

Fig. 2. Global and local Cartesian coordinate systems for the formal asymptotic solution and the boundary layer functions presentation.

p( x ) =

∞

ε 2k p(2k) (x ) +

∞

ε m P ( m ) ( x0 , ς ).

(52)

m=0

k=0

Here, the functions V(m ) (x0 , ς ) = V(m ) (x0 , ς1 , ς2 , ς3 ) and P (m ) (x0 , ς ) = P (m ) (x0 , ς1 , ς2 , ς3 ) vanish, when ς 3 → ∞, and change slowly along the ς 1 , ς 2 -axes. Three first terms of these series that determine the solution with precision of ε have the form

u ( x ) = u ( 0 ) ( x ) + V ( 0 ) ( x0 , ς ) + ε V ( 1 ) ( x0 , ς ),

(53)

p(x ) = p(0) (x ) + P (0) (x0 , ς ) + ε P (1) (x0 , ς ). Substituting these functions into the original system (9) and (10) and taking into account that the functions satisfy system (49) and (50), we obtain

(54) u(0) (x),

p(0) (x)

Lζ V0 + ε Lζ V1 − ε α∇ς P (0) + O(ε 2 ) = 0,

(55)

α divς V0 + εα divς V1 + ε (iς + μβ )P (0) + O(ε 2 ) = 0.

(56)

Here, Lζ , divς , and ς are the differential operators (48), of divergence, and Laplace with respect to ς -variables. It is taken into account that

∂ 1 ∂ = . ∂ zi ε ∂ςi

(57)

Joining the terms of the same order with respect to ε , we obtain

Lζ V0 = 0,

α divς V0 = 0,

Lζ V1 − α∇ς P (0) = 0,

(58)

α divς V1 + (iς + μβ )P (0) = 0.

(59)

As the boundary layer functions vanish along the ς 3 -axis faster than along the ς 1 , ς 2 -axes, the following inequality holds

∂ ∂ ∂ ,  . ∂ς1 ∂ς2 ∂ς3

(60)

As a result, Eqs. (58) and (59) take the form

∂ 2V3(0) ∂ V3(0) = 0 , = 0, ∂ς3 ∂ς32

(61)

λ + 2μ ∂ 2V3(1) ∂ P (0 ) −α = 0, 2 μ ∂ς3 ∂ς3

α

  2 ∂ V3(1) ∂ + i + μβ P ( 0 ) = 0. ∂ς3 ∂ς32

(62)

It follows from Eq. (61) and condition at infinity, that the function V0 is equal to zero, and from two Eq. (62)

∂ 2V3(1) αμ ∂ P (0) = , λ + 2μ ∂ς3 ∂ς32

 2  ∂ ∂ α2μ i + μβ + P ( 0 ) = 0. ∂ς3 ∂ς32 λ + 2μ

(63)

194

S. Kanaun et al. / International Journal of Engineering Science 128 (2018) 187–207

The vanishing at infinity solution of the last equation is

P (0) (x0 , ς3 ) = A(x0 ) exp(−iqs ες3 ), qs = M=

(1 − i ) ε



μ 2M

(α 2 + Mβ ),

λ + 2 μ,

(64) (65)

where A(x0 ) is an arbitrary constant. Note that qs in this equation coincides with the principal term of the expansion of the s -wave number of the complete dispersion Eq. (16) with respect to ε (18). Thus, the “zero”-term of the asymptotic solution has the form

u(x ) = u(0) (x ), p(x ) = p(0) (x ) + P (0) (x0 , ς3 ),

(66)

where the functions u(0) (x) and p(0) (x) satisfy the system (49) and (50), and the function P(0) (x0 , ς 3 ) is defined in the local coordinate system at each point x0 of the inclusion interface  by Eq. (64). In the case of the scattering problem, we have to construct similar asymptotic solutions in the host medium and the inclusion

u1 (x ) = u1(0) (x ),

p1 (x ) = p(10) (x ) + A1 (x0 ) exp(−iqs1 |z3 | ),

(67)

u2 (x ) = u2(0) (x ),

p2 (x ) = p(20) (x ) + A2 (x0 ) exp(−iqs2 |z3 | ),

(68)

where

(1 − i ) qs1 = ε1



μ1

(1 − i ) (α + M1 β1 ), qs2 = 2 M1 ε2



2 1

μ2 2 M2

(α12 + M2 β2 ).

(69)

The parameters with the index 1 correspond to the inclusion material and the index 2 corresponds to the host medium. The functions u1(0 ) (x ), p(10 ) (x ) and u2(0 ) (x ), p(20 ) (x ) allow us to satisfy the boundary conditions (25) and (26) on the inclusion interface. Meanwhile, the boundary conditions (27) can be satisfied by the appropriate constants A1 (x0 ) and A2 (x0 ) in Eqs. (66) and (67). Continuation of the pressure p(x) on the inclusion surface yields the equation

p(10) (x0 ) + A1 (x0 ) = p(20) (x ) + A2 (x0 ), x0 ∈ ,

(70)

and the second condition (27) takes the form

1 ω2 ρ1



  (0 )  ρ f 1  0 ρ f 2  0 ∂ p(10) ∂ p2 1 − iqs1 A1 (x0 ) − u1 = 2 + iqs2 A2 (x0 ) − u . 3 ∂ z3 ρ1 ρ2 2 3 ω ρ2 ∂ z3

(71)

Keeping in this equation the principal terms with respect to ε , we obtain

−

qs1 qs2 A (x ) = A ( x ). ρ1 1 0 ρ2 2 0

(72)

Finally, we derive explicit equations for the constants A1 (x0 ) and A2 (x0 ) from (70) and (72)

A1 ( x0 ) =

 1  (0 ) 0 p2 (x ) − p(10) (x0 ) , 1+δ

A2 ( x0 ) = −

δ=

2 qs1 ρ , 1 qs2 ρ

(73)

 δ  (0 ) 0 p2 (x ) − p(10) (x0 ) . 1+δ

(74)

5. Zero-order term of the asymptotic expansion in the case of a spherical inclusion In the case of a spherical inclusion, a solution of the system (49) and (50) can be found in the form of series of spherical harmonics similar to the general case. The components of the displacement field inside the inclusion are presented in the form (0 ) ur1 =

∞  1 0 1 0 0 Bn f f1n (q f 1 r ) + B0nt f21n (qt1 r ) Pn (cos θ ), r

(75)

∞  dPn (cos θ ) 1 0 1 0 0 Bn f f3n (q f 1 r ) + B0nt f41n (qt1 r) . r dθ

(76)

n=0

uθ(01) =

n=0

The pressure in the inclusion takes the form

p(10) = −

∞ μ1

r2

n=0

Bn(0f ) f91n (q0f 1 r, s f 1 )Pn (cos θ ), s f 1 =

0 qt1

q0f 1

.

(77)

S. Kanaun et al. / International Journal of Engineering Science 128 (2018) 187–207

195

Table 1 Wave numbers of transverse qt2 , fast qf2 , and slow qs2 longitudinal waves in the host medium; corresponding values of the parameter ε 2 ; and the number N of terms in the series (30)–(35) of the complete solution of the scattering problem. qt 2

qf2

qs2

ε22

N

0.1 0.5 1 2

0.0597 0.298 0.597 1.194

35.96(1 − i ) 80.41(1 − i ) 113.72(1 − i ) 160.83(1 − i )

5.7 · 10−4 1.1 · 10−4 5.7 · 10−5 2.8 · 10−5

51 100 135 185

The displacements and the pressure in the host medium are (0 ) ur2 =−

∞ U0 (−i )n+1 (2n + 1 ) f11n (q0f 2 r )Pn (cos θ ) q f 2r n=0

∞  a 0 2 0 2 0 + Cn f f1n (q f 2 r ) + Cnt f2n (qt2 r ) Pn (cos θ ), r

(78)

n=0

uθ(02) = −

∞ dP (cos θ ) U0 (−i )n+1 (2n + 1 ) f31n (q0f 2 r ) n q f 2r dθ n=0

∞  dPn (cos θ ) a 0 2 0 0 2 0 + Cn f f3n (q f 2 r ) + Cnt f4n (qt2 r) , r dθ

(79)

n=0

p(20) =

∞ μ2 U0

q0f 2 r 2 a −

∞ μ2

r2

(−i )n+1 (2n + 1 ) f9n (q0f 2 r, s f 2 )Pn (cos θ )

n=0

Cn f f92n (q0f 2 r, s f 2 )Pn (cos θ ).

(80)

n=0

0 are defined by the equations Here, the wave numbers q0f i and qti

q0f i

= aω



ρti λi + 2μi

 ,

qti0

= aω

ρti , i = 1, 2. μi

(81)

0 and B0 , B0 in Eqs. (75)–(80) The boundary conditions (25) and (26) give the following system for the constants Cn0 f , Cnt nt nf

M ( 0 ) ( n )X 0 ( n ) = F 0 ( n ),

(82)

where

⎡ f 2 ( q0 ) 1n f2 ⎢ f32n (q0f 2 ) 0 M (n ) = ⎢ 2 0 ⎣ f 5n ( q f 2 , s f 2 ) f72n (q0f 2 ) ⎡

Cn0 f

⎤

⎢ Cnt0 ⎥ X 0 (n ) = ⎢ 0 ⎥, ⎣B ⎦ nf B0nt

0 f22n (qt2 ) 0 f42n (qt2 )

f62n f82n

( (

0 qt2 0 qt2

) )

− f11n (q0f 1 ) −ξ

− f31n (q0f 1 )

−ξ

(

q0f 1 , s f 1) f71n q0f 1

f51n

(

)

0 − f21n (qt1 )

⎤

⎥ ⎥, )⎦

0 − f41n (qt1 )

)

−ξ −ξ

f61n f81n

⎡ f 1 ( q0 ) ⎤ 1n f2 ⎢ f3n (q0f 2 ) ⎥ 0 n+1 2n + 1 ⎢ ⎥. F (n ) = U0 (−i ) q0f 2 ⎣ f51n (q0f 2 , s f 2 )⎦ f7n (q0f 2 )

( (

0 qt1 0 qt1

(83)

)

(84)

The wave numbers of the transverse waves qt2 = 0.1, 0.5, 1, 2, the decimal logarithms of the conditional numbers d0 (n) of the matrix M0 (n) and the number of terms that is necessary to keep in Eqs. (75)–(80) are presented in Table 4. 0 = 0.1 and n = 0, 1, 2 are given in Table 5. The values of the constants Cnf , Cnt , Bnf , Bnt for qt2 Comparison of these Tables with Tables 1–3 shows that the conditional numbers of the matrices M0 (n) are smaller than the ones of the M(n), and the required numbers of terms in the series for displacements and pressure are also much smaller for the asymptotic solution.

196

S. Kanaun et al. / International Journal of Engineering Science 128 (2018) 187–207 Table 2 Decimal logarithms of the conditional numbers d(n) of the matrices M(n) in Eq. (36) for the coefficients in the complete solution of the scattering problem. The incident field frequency corresponds to the transverse wave number in the host medium qt2 = 0.1. qt 2

lgd(0)

lgd(1)

lgd(2)

lgd(3)

lgd(4)

lgd(5)

lgd(6)

0.1 0.5 1 2

22.04 33.36 33.38 33.82

22.24 32.23 33.65 34.88

19.96 30.90 33.54 35.26

17.66 29.74 33.66 34.14

23.87 28.69 33.89 34.89

28.54 26.60 34.41 35.53

32.98 25.81 33.42 36.36

Table 3 Coefficients in the complete solution of the scattering problem (30)–(35). The incident field frequency corresponds to the transverse wave number in the host medium qt2 = 0.1. n

0

1 −8

−4

2 −8

−2.4 · 10 +1.05 · 10 (2.5 − 4.1i)·109

Cnf Cns

8 · 10 − 5.1 · 10 (6.1 − 8.4i ) · 109

Cnt Bnf Bns Bnt

0.67 − 16.8i 8 · 10−4 +5.9i (−1.36 − 5.5i ) · 10−12 −1.03 − 6.09i

i

−4

−7

−8.1 · 10−4 + 2.3 · 10 i (3.6 − 4.1i ) · 1010

i

−7

−7.06 · 10−8 + 2.9i · 10−4 16.59 + 5.9 · 10−4 i (2.3 + 3.9i ) · 10−12 1.76 + 6.03 · 10−6 i

−1.9 · 10−3 + 5.5 · 10 i −3.07 − 22.96i (−1.1 + 0.13i ) · 10−11 0.31 − 0.7i

Table 4 Decimal logarithms of the conditional numbers d0 (n) of the matrices M0 (n) in Eq. (36) for the coefficients in the asymptotic solution of the scattering problem; N is the number of terms in Eqs. (75)-(80) of the asymptotic solution. The incident field frequency corresponds to the transverse wave numbers in the host medium qt2 = 0.1, 0.5, 1, 2. 0 qt2

lgd0 (0)

lgd0 (1)

lgd0 (2)

lgd0 (3)

lgd0 (4)

lgd0 (5)

lgd0 (6)

N

0.1 0.5 1 2

4.08 2.30 1.21 1.18

6.88 3.39 1.97 1.12

10.68 5.79 3.74 2.06

14.75 8.47 5.79 3.26

18.99 11.31 8.02 4.81

23.37 14.29 10.41 6.58

27.88 17.4.45 12.91 8.46

4 5 6 7

Table 5 Coefficients in the asymptotic solution of the scattering problem (75)-(80). The incident field frequency corresponds to the transverse wave number in the host medium qt2 = 0.1. n

0

Cnf Cnt Bnf Bnt

8 · 10−4 − 3.5 · 10 0.672 − 16.77i 2.4 · 10−4 +5.93i −1.03 − 6.08i

−8

i

1

2

−1.9 · 10−8 +1.05 · 10−5 i −5.4 · 10−8 +2.96 · 10−3 i 16.59 + 2.9 · 10−4 i 1.76 + 8.4 · 10−5 i

−8.1 · 10−4 + 1.6 · 10 i −1.9 · 10−3 +3.8 · 10−7 i −3.2 · 10−3 −20.48i −2.6 · 10−4 −0.95i

−7

6. Wave fields in the vicinity of the inclusion In this section, the straightforward solution of the complete system of equations of the scattering problem (30)–(36) and the “zero”-terms of the asymptotic expansion with respect to small parameter ε defined in Eqs. (67), (68), (75)–(80), and (82) are compared. In Figs 3–6, the distributions of the u3 -component of the displacement vector along the x3 -axis are presented for dimensionless wave numbers qt0 = 0.1, 0.5, 1, 2,



qt0 = ωa

ρ2 , μ2

(85)

and μ2 , ρ 2 are the effective shear modulus and density of the host medium. In these figures, solid lines are the displacements calculated from the solution of the complete system (36), and dashed lines are the zero-order-terms of the asymptotic expansion with respect to parameter ε (system (82)). It is seen that these solutions coincide practically inside the inclusion region |x3 | ≤ 1, as well as outside it 2 > |x3 | > 1. For the pressure, the difference in solutions of the complete system and the “zero”-term of the asymptotic expansion can be observed only in the vicinity of the inclusion interface. In Figs 7–10, solid lines correspond to the solution of the complete system, dashed lines are the “zero”-term of asymptotic expansions p(0) (0, 0, x3 )/Kf for pressure without the boundary layer functions; Kf is the bulk modulus of the fluid inside the inclusion.

S. Kanaun et al. / International Journal of Engineering Science 128 (2018) 187–207

197

Fig. 3. Distribution of the u3 -component of the displacement field along the x3 -axis in the vicinity of a spherical inclusion; Re and Im are the real and imaginary parts of u3 (0, 0, x3 )/U0 ; solid lines are the solution of the complete system; dashed lines are the “zero”-term of the asymptotic solution; qt0 = 0.1 is the dimensionless wave number of the transverse field in the host medium.

Fig. 4. The same as in Fig. 3 for qt0 = 0.5.

Comparisons of the solution of the complete system and the sum of p(0) (x) and boundary layer function P(0) (x0 , ς 3 ) in the vicinity of the boundary point x0 = (0, 0, −1 ) are presented in Figs 11–14 for qt0 = 0.1, 0.5, 1, 2; −0.8 > x3 > −1.2. Maximal discrepancy of these two solutions takes place in the vicinity of the inclusion interface x3 = −1. 7. Scattered wave fields far from the inclusion Equations for the scattered fields in the host medium follow from (33)–(35)

ur(s ) =

∞  1

Cn f f12n (q f 2 r ) + Cns f12n (qs2 r ) + Cnt f22n (qt2 r ) Pn (cos θ ), r n=0

(86)

198

S. Kanaun et al. / International Journal of Engineering Science 128 (2018) 187–207

Fig. 5. The same as in Fig. 3 for qt0 = 1.

Fig. 6. The same as in Fig. 3 for qt0 = 2.

uθ(s ) =

∞  dPn (cos θ ) 1

Cn f f32n (q f 2 r ) + Cns f32n (qs2 r ) + Cnt f42n (qt2 r ) , r dθ

p(s ) = −

(87)

n=0

∞

μ2

r2



Cn f f92n (q f 2 r, s f 2 ) + Cns f92n (qs2 r, ss2 ) Pn (cos θ ),

(88)

n=0

Further, we omit the index 2 for the host medium parameters. Using the asymptotic formula for spherical Hankel functions of large argument

hn(2) (z ) ∼ −in−1

exp(−iz ) , z

(89)

S. Kanaun et al. / International Journal of Engineering Science 128 (2018) 187–207

199

Fig. 7. Distribution of the pressure p along the x3 -axis in the vicinity of a spherical inclusion; Re and Im are the real and imaginary parts of p(0, 0, x3 )/Kf ; Kf is the bulk modulus of the fluid in the inclusion; solid lines are the solution of the complete system; dashed lines are the “zero”-term of the asymptotic solution without boundary layer functions; qt0 = 0.1 is the dimensionless wave number of the transverse field in the host medium.

Fig. 8. The same as in Fig. 7 for qt0 = 0.5.

and neglecting the s-waves that have large attenuation factors, we find for r  1

ur(s ) =

∞  1 n

i Cn f e−iq f r + Cnt e−iqt r Pn (cos θ ) + O(r −2 ), r

(90)

∞ e−iqt r n dPn (cos θ ) i Cnt + O(r −2 ), r dθ

(91)

n=0

uθ(s ) =

p(s ) = −

n=0

∞ ρt q2f  q f e−iq f r αμ in+1Cn f Pn (cos θ ) + O(r −2 ), γ = (γ − βμ ) r n=0 ρqt2

(92)

200

S. Kanaun et al. / International Journal of Engineering Science 128 (2018) 187–207

Fig. 9. The same as in Fig. 7 for qt0 = 1.

Fig. 10. The same as in Fig. 7 for qt0 = 2.

As a result, the far scattered field of displacements and pressure can be presented in the form

u (s ) = A

e−iq f r e−iqt r e−iq f r +B , p(s ) = C , r r r

(93)

where

A = nA, B = eθ B, C=−

∞

A= B=

 qf αμ γ − μβ

inCn f Pn (cos θ ),

n=0 ∞

n=0 ∞

inCnt

n = er =

dPn (cos θ ) , dθ

in+1Cn f Pn (cos θ ) = −i

n=0

x , r

(94) (95)

 qf αμ A. γ − μβ

(96)

S. Kanaun et al. / International Journal of Engineering Science 128 (2018) 187–207

201

Fig. 11. Distribution of the pressure p along the x3 -axis in the vicinity of point x0 = (0, 0, −1 ) on the interface of a spherical inclusion; Re and Im are the real and imaginary parts of p(0, 0, x3 )/Kf ; Kf is the bulk modulus of the fluid in the inclusion; solid lines are the solution of the complete system; dashed lines are the “zero”-term of the asymptotic solution with boundary layer functions; qt0 = 0.1 is the dimensionless wave number of the transverse field in the host medium.

Fig. 12. The same as in Fig. 11 for qt0 = 0.5.

Here, er , eθ are the basis vectors in spherical coordinates with the polar axis directed along the x3 -axis. 8. Energy flux of the scattered field Let us consider equations for the densities of kinetic T and potential  energy of a poroelastic medium (Biot, 1962)

T =

=

. . 1 .2 1 .2 ρ u + ρ f ui wi + ρwi , 2 i 2

1 1 1 2 C ei j ekl − α ekk ζ + ζ . 2 i jkl β 2β

(97) (98)

202

S. Kanaun et al. / International Journal of Engineering Science 128 (2018) 187–207

Fig. 13. The same as in Fig. 11 for qt0 = 1.

Fig. 14. The same as in Fig. 11 for qt0 = 2.

Here

ei j = ∂ ( i u j ) ,

ζ = −∂i wi ,

Ci jkl = λδi j δkl + μ(δik δil + δil δik ).

(99)

Using the constitutive equations of poroelasticity

σi j = Ci jkl ekl −

α ζ, β

p=−

α 1 e + ζ, β kk β

(100)

we have from Eqs. (98) and (100)

∂ ∂ = σi j , = p. ∂ ei j ∂ζ

(101)

S. Kanaun et al. / International Journal of Engineering Science 128 (2018) 187–207

203

The total energy of a spherical domain V of the poroelastic medium is the following integral



E=

(T + )dV.

V

(102)

The time derivative of the total energy is

    . .. . .. . .. . .. ∂ . ∂ . ∂T ∂ + dV = e + ζ dV. ρ ui ui + ρ f (wi ui + ui wi ) + ρwi wi + ∂t ∂t ∂ ei j i j ∂ζ V





∂E = ∂t

V

(103)

Using Gauss’s theorem, the volume integral of the time derivative of the potential energy can be transformed as follows

 V

 . .  ∂u ∂w σi j i − p i dV ∂xj ∂ xi V    . . . ∂σ . ∂p ∂ ∂ = ( σi j ui ) − ( pw j ) − ui i j + wi dV ∂xj ∂xj ∂xj ∂ xi V   ˙  . . . ∂σi j . ∂p = (σi j u j − pwi )mi dS − ui + wi dV . ∂xj ∂ xi S V

∂ dV = ∂t



(104)

Here, mi is the normal to surface S of the domain V. Thus, for the time derivative of the total energy, we obtain



dE = dt Ii =

     . .. .. . .. .. ∂σi j ∂p wi + ui ρ ui + ρ f wi − + wi ρ f ui + ρ dV, ∂xj ∂ xi V

 S

Ii mi dS +

.

.

σi j u j − pwi .

(105) (106)

It follows from poroelastic equations of motion (Biot, 1962) that the volume integral in this equation is equal to zero, and as a result, we obtain

dE = dt



S

Ii mi dS.

(107)

Here, Ii is the energy flux vector. For the incident field, the average rate of the energy transfer per unit area normal to the direction of propagation is defined by the equation

I (i ) =



 σi(ji) u˙ (ji) − p(i) w˙ i(i) ni.

(108)

The normalized total scattering cross-section Q(ω) of the inclusion is the ratio of the power flux scattered in all directions, averaged over the oscillation period, to the average intensity of the incident field

Q (ω ) =

1



π I (i )





 S

1

Ii mi dS =   (i )



I



2

S1

Ii mi R dS1 ,

(109)

where · means averaging over the period of oscillations; R is the radius of the domain V; S1 is the surface area of a sphere of unit radius. The dimensionless differential cross-section F for a spherical inclusion is defined by the equation

F=

 R2  Ii mi 2   R =   mi (σi(js) u˙ (js) − p(s) w˙ i(s) , ( i ) ( i ) πI πI

(110)

Since the time dependence of the fields is determined by the factor eiωt , the numerator in this equation is





mi (σi(js ) u˙ i(s ) − p(s ) w˙ i(s ) ) =

where

∗

ω 2



Im mi

 σi(js) u(js)∗ − p(s) wi(s)∗ ,

(111)

means complex conjugation. Substituting ui(s ) from Eq. (93) into Eq. (111) and using constitutive Eq. (99), we obtain

F =− =−

  iq f ω 2 ∗CA∗ +   Im iq f M|A|2 + iμ|B|2 + α C | | ρ∗ ω2 2π I ( i )

   2 μ2 β ω ρt q2f |A|2 + iμ|B|2 , M  = M + |α|   Im iq f M , γ = . ρ qt2 |γ − μβ|2 2π I ( i )

(112)

By deriving this equation, the imaginary parts of the wave numbers qf and qt are neglected. For a longitudinal incident wave, we have

u(i ) = U0 e3 exp(−iq f e3 · x ), p(i ) = P0 exp(−iq f e3 · x ), P0 = −

 iq f αμ U , γ − μβ 0

(113)

204

S. Kanaun et al. / International Journal of Engineering Science 128 (2018) 187–207

Fig. 15. Differential cross-section F(θ ) of a spherical inclusion for an incident longitudinal wave directed along the polar axis of the spherical coordinate system; solid line is the solution of the complete system; dashed line is the “zero”-term of the asymptotic solution; the frequency of the incident field corresponds to the dimensionless transverse wave number in the host medium qt0 = 0.1.

Fig. 16. The same as in Fig. 15 for qt0 = 0.5.

and the time averaged intensity I(i) is

 (i )  I

=

ω 2

=−

 

Im ni

ω 2



σi(ji) u(ji)∗ − p(i) wi(i)∗



∗ P0U0∗ − Im iq f M|U0 |2 + α

 γ∗ ω  2 U0 . |P0 |2 = − Im iq f M iq f μ 2

(114)

Finally, the differential cross-section F(ω, θ ) and the normalized total scattering cross-section Q(ω) take the form|

F=



1





 U2 π Im iq f M 0

Q (ω ) =



4





|A|2 + iqt μ|B|2 , Im iq f M

 U2 Im iq f M 0

Im

∞ n=0

 i  q f M|Cn f |2 + qt μn(n + 1 )|Cnt |2 . 2n + 1

(115)

(116)

S. Kanaun et al. / International Journal of Engineering Science 128 (2018) 187–207

205

Fig. 17. The same as in Fig. 16 for qt0 = 1.

Fig. 18. The same as in Fig. 17 for qt0 = 2.

Graphs of the function F(θ ) for the values of dimensionless wave number qt0 = 0.1, 0.5, 1, 2 are presented in Figs 15– 18. In these figures, solid lines are obtained from the solution of the complete system of equations of poroelasticity, dashed lines correspond to the “zero”-terms of the asymptotic expansion with respect to small parameter ε . Note that the boundary layer functions don’t contribute to the far scattered field characteristics.

206

S. Kanaun et al. / International Journal of Engineering Science 128 (2018) 187–207

9. Conclusion The problem of scattering of plane longitudinal waves on a spherical inclusion in poroelastic media is considered for typical material properties of sedimentary rocks. In this case, the system of equations of poroelasticity contains a differential operator of the highest order with a small parameter. The method of matched asymptotic expansions is used for construction of a reliable solution of the problem. It is shown that the straightforward solution of the original problem and the “zero”term of the asymptotic expansion give close results for the displacements and pressure inside the inclusion, as well as in the host medium. Discrepancy of the two solutions appears for the fluid pressure in a small vicinity of the inclusion interface. This result can be expected, because the small parameter appears in the differential operator that acts only on the pressure. The numerical results presented in the paper were obtained by “Mathematica” package (https://www.wolfram.com/ mathematica/). High precision numerical calculations supported by the package allow one to avoid notable errors in the straightforward solution of the original problem. In spite of large conditional numbers of the matrices for the coefficient of the series, the solution turns out to be close to the zero-order term of the asymptotic expansion. It is not guaranteed that substantial numerical errors cannot appear in the straightforward solution for other medium parameters or by application of other numerical methods. Meanwhile, the method of matched asymptotic expansions provides a stable and more reliable algorithm in a wide range of poroelastic parameters of a medium, inclusions and frequencies of the incident field. Appendix The basic functions used for presentation of the solution of the scattering problem for a spherical inclusion in the form of series of spherical harmonics are defined as follows

f1i n (z ) = nyin (z ) − zyin+1 (z ), f2i n (z ) = n(n + 1 )yin (z ), f3i n (z ) = yin (z ),

f4i n (z ) = (n + 1 )yin (z ) − zyin+1 (z ),





αi i (z, s ) = n − n − yin (z ) + 2zyin+1 (z ) − f (z, s ), 2 2 9n

 i i i f6n (z ) = n(n + 1 ) (n − 1 )yn (z ) − zyn+1 (z ) , f5i n

(zs )2

2

f7i n (z ) = (n − 1 )yin (z ) − zyin+1 (z ),



f8i n (z ) =

n2 − 1 −

i z2 f9i n (z, s ) = α

ρ

 2

z 2

(A.1)

yin (z ) + zyin+1 (z ),

1 − μi βi s2

−1

yin (z ), i  ρ   −1 ρfi i ρti ρti i 2  f10 ( z, s ) = α − s μ β − f ( z ), i i n ρi i ρi ρti 1n y1n (z ) = jn (z ); y2n (z ) = hn(2) (z ). ti

·

Here, jn (z) are the spherical Bessel functions; hn(2 ) (z ) are the Hankel functions of the second kind; i = 1, 2; the index “i = 1” is for the inclusion; the index “i = 2” is for the host medium. The components σ rr and σ rθ of the stress tensor in the host medium have the form

σrr = −

∞

2 μ2 U0 (−i )n+1 (2n + 1 ) f51n (q f 2 r, s f 2 )Pn (cos θ ) 2 q f 2r n=0

∞  2 μ1

+ 2 Bn f f51n (q f 1 r, s f 1 ) + Bns f51n (q2s1 r, ss1 ) + Cnt f61n (qt1 r ) Pn (cos θ ), r n=0

σr θ = −

(A.2)

∞

2 μ2 dP (cos θ ) U0 (−i )n+1 (2n + 1 ) f71n (q f 2 r ) n 2 dθ q f 2r n=0

∞  dPn (cos θ ) 2 μ1

+ 2 Bn f f71n (q f 1 r ) + Bns f71n (q2s1 r ) + Cnt f81n (qt1 r ) , dθ r

(A.3)

n=0

and the component wr of the vector w is ∞

wr =

2 μ2 1 U (−i )n+1 (2n + 1 ) f10 n (q f 2 r, s f 2 )Pn (cos θ ) q f 2 r2 n=0

∞  1

1 1 − Bn f f10 n (q f 1 r, s f 1 ) + Bns f 10n (qs1 r, ss1 ) Pn (cos θ ). r n=0

(A.4)

S. Kanaun et al. / International Journal of Engineering Science 128 (2018) 187–207

207

The stresses inside the inclusion are

σrr(r ) =

∞  2 μ1

Bn f f51n (q f 1 r, s f 1 ) + Bns f51n (q2s1 r, ss1 ) + Cnt f61n (qt1 r ) Pn (cos θ ), 2 r

(A.5)

∞  Pn (cos θ ) 2 μ1

Bn f f71n (q f 1 r ) + Bns f71n (q2s1 r ) + Cnt f81n (qt1 r ) . dθ r2

(A.6)

n=0

σr(θr ) =

n=0

The relative fluid displacement with respect to the solid skeleton in the inclusion is

wr(r ) = −

∞  1

1 1 Bn f f10 n (q f 1 r, s f 1 ) + Bns f 10n (qs1 r, ss1 ) Pn (cos θ ). r

(A.7)

n=0

In these equations

sfi =

qti qti , ssi = , i = 1, 2. qfi qsi

(A.8)

Supplementary material Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.ijengsci.2018.03. 003. References Berryman, J. (1985). Scattering by a spherical inhomogeneity in a fluid-saturated porous medium. Journal of Mathematical Physics, 26, 1408–1419. Biot, M. (1956). Theory of propagation of elastic waves in a fluid-saturated porous solid. i. low-frequency range. Journal Acoustical Society of America, 28, 168–178. Biot, M. (1962). Mechanics of deformation and acoustic propagation in porous media. Journal of Applied Physics, 33/4, 1482–1498. Bohren, C., & Huffman, D. (1983). Absorption and scattering of light by small particles p. 660. NY: Wiley&Sons. Bourbie, T., Coussy, O., & Zinzner, B. (1987). Acoustics of porous media p. 334. Houston: Gulf Publishing Co. Mei, Ch., & Foda, M. (1981). Wave-induced responses in a fluid-filled poro-elastic solid with a free surface - a boundary layer theory. Geophysics Journal of Researches astr. Society, 66, 597–631. Ciz, R., Gurevitch, B., & Markov, M. (2006). Seismic attenuation due to wave-induced fluid flow in a porous rock with spherical heterogeneities. Geophysical Journal International, 165, 957–968. Derisiewicz, H., & Skalak, R. (1963). On uniqueness in dynamic poroelasticity. Bulletin of Seismological Society of America, 53, 783–788. Golub, G., & Ch, V. L. (1989). Matrix computations p. 548. The John Hopkins University Press. Il’yin, A. (1992). Matching of asymptotic expansions of solutions of boundary value problems (p. 281). AMS-Books. Kanaun, S., & Levin, V. (2008). Effective field method for composites, v.1 Static problems (p. 376). Dordrecht: Springer. Krutin, V., Markov, M., & Yumatov, A. (1984). Scattering of longitudinal wave by a spherical cavity with a fluid in an porous saturated medium. Journal of Applied Mathematics and Mechanics (PMM), 48, 238–241. Liu, X., Greenhalgh, S., & Zhou, B. (2009). Scattering of plane transverse waves by spherical inclusions in a poroelastic medium. Geophysical Journal International, 176, 938–950. Mavko, G., Mukerji, T., & Dvorkin, J. (2009). The rock physics handbook. tools for seismic analysis in porous media (2nd ed., p. 400). Cambridge University Press. Morochnic, V., & Barded, J. (1996). Viscoelastic approximation of poroelastic media for wave scattering problems. Soil Dynamics and Earthquake Engineering, 15, 337–346. Pride, S., & Haartsen, M. (1996). Electroseismic wave properties. Journal of Acoustical Society of America, 100/3, 1301–1315. Zimmerman, C., & Stern, M. (1993). Scattering of plane compressional waves by spherical inclusions in poroelastic medium. Journal of Acoustical Society of America, 94/1, 527–536.