Truncated transparent boundary conditions for second order hyperbolic systems

Applied Numerical Mathematics 120 (2017) 115–124

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

Truncated transparent boundary conditions for second order hyperbolic systems Ivan Sofronov a,b,∗ a b

Schlumberger, Moscow, Pudovkina 13, Russian Federation MIPT, Moscow region, Dolgoprudny, Institutskii per. 1, Russian Federation

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 29 September 2016 Received in revised form 4 May 2017 Accepted 5 May 2017 Available online 11 May 2017 Keywords: Transparent boundary conditions Local boundary conditions Hyperbolic systems Cylindrical anisotropy Orthotropic elasticity Biot poroelasticity

In [22] we announced equations for yielding differential operators of transparent boundary conditions (TBCs) for a certain class of second order hyperbolic systems. Here we present the full derivation of these equations and consider ways of their solving. The solutions represent local parts of TBCs, and they can be used as approximate nonreflecting boundary conditions. We give examples of computing such conditions called ‘truncated TBCs’ for 3D elasticity and Biot poroelasticity © 2017 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction Transparent boundary conditions (TBCs) for the wave equation were proposed in [18,10] and were further developed, for example, in [20,11,1,25]. Also different aspects of deriving and exploring TBCs are considered in [4,6,12,14,17,26] (note that Th. 2.1 in [26] repeats formulas of kernels derived in [20]). As it is known, on a circle (d = 2) or a sphere (d = 3) TBCs can be written as

∂u ∂u d − 1 c + + u − Q −1 { B k ∗} Q u = 0 ∂t ∂r 2 r

(1)

where Q and Q −1 are the Fourier and inverse Fourier transform operators over the trigonometric (d = 2) or spherical (d = 3) basis, and { B k ∗} denotes the time convolution operator of the k-th Fourier coefficient with kernels B k (t ) derived analytically [20]. In order to provide efficient calculations with reasonable computational recourses the kernels B k (t ) are approximated by sums of exponentials,

B k (t ) ≈

Lk 

al,k ebl,k t ,

Re bl,k ≤ 0.

l =1

*

Correspondence to: Schlumberger, Moscow, Pudovkina 13, Russian Federation. E-mail address: [email protected].

http://dx.doi.org/10.1016/j.apnum.2017.05.002 0168-9274/© 2017 IMACS. Published by Elsevier B.V. All rights reserved.

116

I. Sofronov / Applied Numerical Mathematics 120 (2017) 115–124

Fig. 1. External boundary Γ , local coordinates at the point O .

Note that the operator of (1) is naturally decomposed into a local (differential) part, the first three terms, which is nothing more than the well-known characteristic boundary conditions for the wave equation, and a nonlocal part in space and time; the latter can be treated as a correction to ensure the absolute accuracy of TBCs. The explicit analytical expressions for other examples of TBCs operators can be obtained only for simple boundaries (plane, sphere, etc.) and for specific narrow classes of equations, see the discussion in [13]. Extension to the equations with variable coefficients is possible by using the so-called quasi-analytic TBCs proposed in [21,23]. Operator coefficients of these conditions are found by solving numerically some auxiliary initial-boundary value problems; therefore, it is a costly computational procedure. A compromise between the generality of the governing equations that we consider and the accuracy of the boundary conditions is possible if the goal is to derive only the local (differential) part of TBCs and to use it as an operator of local approximate non-reflecting boundary conditions. We call such conditions truncated TBCs (TTBCs), as they do not contain the nonlocal part of TBCs. A natural and important generalization of the wave equation is the general second-order hyperbolic systems. In [22], we have introduced the equations for evaluating the local part of the TBC operator for such systems. Moreover, we have derived TTBCs for 2D VTI cylindrical anisotropy and for 3D Navier wave equation; a practical application of TTBC for the latter case has been reported in [24]. In the current paper, we provide a detailed derivation of the equations first introduced in the short note [22]. Besides, we present a solution procedure for these equations that yields the operators of the differential part of the TBCs. This is the subject of Section 2. In Section 3, we consider two additional examples of generating TTBCs: for 3D orthotropic elasticity equations in the cylindrical coordinates, and for 3D Biot poroelasticity equations in the Cartesian coordinates. Note that for a general hyperbolic system of first-order pseudodifferential operators with variable coefficients a way of computing the approximate transparent boundary conditions is considered in [2]. Taking advantage of the recursive formulation of Taylor’s method, the authors derive a recursive formula expressing the mth homogeneous symbol of the nonlocal pseudodifferential operator of transparent conditions. In particular, applying to Maxwell’s equations they derive a local condition that they call “a second-order complete radiation condition”. In [3], the authors consider the time dependent 2D Klein–Gordon equation and derive a set of absorbing boundary conditions using microlocal analysis and pseudodifferential operators. It is an interesting question of whether the techniques developed in [2] and [3] can be applied to the secondorder hyperbolic systems we are considering hereafter. 2. Derivation of equations for operators of TTBCs 2.1. Problem formulation Consider a domain Ω ⊂ R 3 with a sufficiently smooth open boundary Γ at the point O . Denote by n the outward normal and by τ = (τ1 , τ2 ) the vector of orthogonal coordinates in a tangent plane, see Fig. 1. Let a vector function u = (u 1 , u 2 , . . . , u N )T satisfies a second-order hyperbolic system of equations [9] written in local coordinates (n, τ ) at the point O in the form

−

∂ 2u ∂u ∂ 2u ∂u + A 2 + ( B ∇τ ) +C + a(∇τ , ∇τ )u + (c ∇τ )u + du = 0. 2 ∂n ∂n ∂t ∂n

(2)

Here we use the following notation:

( B ∇τ )u ≡

2  i =1

Bi

∂u , ∂ τi

a(∇τ , ∇τ )u ≡

2  2  i =1 j =1

ai j

∂ 2u . ∂ τi ∂ τ j

The N × N matrix coefficients A (τ ), B i (τ ) and C (n, τ ), ai j (n, τ ), c i (n, τ ), and d(n, τ ) may depend on τ and on n, τ , respectively. We suppose that the coefficients and solution of (2) are sufficiently smooth in the neighborhood of the point O . Our goal is to obtain a relationship between the solution u and its first derivatives, which is approximately valid at the point O . This relationship will have the form

I. Sofronov / Applied Numerical Mathematics 120 (2017) 115–124

P1

117

∂u ∂u − + P 0u = 0 ∂t ∂n

with some operators P 1 , P 0 that are obtained from the coefficients of (2). Initial condition of the rest at infinity is posed in order to consider the wave solutions u leaving the domain Ω only. 2.2. Derivation of equations for P 1 , P 0 Applying the Laplace transform with respect to time to (2) we have

∂ uˆ ∂ 2 uˆ ∂ uˆ + ( B ∇τ ) +C + a(∇τ , ∇τ )uˆ + (c ∇τ )uˆ + duˆ = 0. 2 ∂n ∂n ∂n After introducing the vanishing parameter ε = 1/s → +0 and local coordinate transformation η = n/ε we obtain the following equation for the function u˜ (η, τ ) ≡ uˆ (εη, τ ): −s2 uˆ + A

A

∂ u˜ ∂ 2 u˜ ∂ u˜ + ε ( B ∇τ ) + ε C˜ + ε 2 a˜ (∇τ , ∇τ )u˜ + ε 2 (˜c ∇τ )u˜ + ε 2 d˜ u˜ − u˜ = 0 ∂η ∂η ∂ η2

(3)

(the tilde in the coefficients denotes that they are functions of η ). Since ∂∂n C (n, τ ) is bounded in physical coordinates (n, τ ), we evidently have

    ∂ ˜ ∂ C (η, τ ) = ε C (n, τ ) = O (ε ) ∂η ∂n η =0 n =0

as

(4)

ε → 0.

Denote by P 1 + ε P 0 the linear portion with respect to following relationship between u˜ and its normal derivative:

ε of the Poincare–Steklov operator for (3), i.e., consider the

  ∂ ˜u (η, τ , ε ) = ( P 1 + ε P 0 )u˜ (η, τ , ε )|η=0 + o(ε ). ∂η η =0

(5)

To obtain P 1 and P 0 we introduce the linear expansion of u˜ with respect to ε :

˜ (η, τ ) + o(ε ) u˜ (η, τ , ε ) = v˜ (η, τ ) + ε w and substitute it for (3). Combining terms of zero and first degrees of

A

ε we obtain:

˜ ∂ v˜ ∂ 2 v˜ ∂2 w ∂ v˜ ˜ + ε ( B ∇τ ) − v + ε A 2 − εw + εC 0 + o(ε ) = 0. 2 ∂η ∂η ∂η ∂η

(6)

Here we use C 0 (τ ) ≡ C˜ (0, τ ) instead of C (n, τ ) according to (4). Taking the limit as ε → 0 in (6), we come to the following necessary conditions:

A

∂ 2 v˜ − v˜ = 0 ∂ η2

(7)

A

 ∂ v˜  ˜ ∂2 w ˜ + ( B ∇τ ) + C 0 −w = 0. 2 ∂η ∂η

(8)

and

The waves leaving Ω correspond to vanishing branch of solutions to (7) as



η → +∞. These solutions have the form



v˜ (η, τ ) = exp − A −1/2 (τ )η v˜ 0 (τ ) where A −1/2 (τ ) is a positive square root of A −1 (τ ), and v˜ 0 (τ ) is a function. Thus we obtain the following Poincare–Steklov map for v˜ :

 ∂ v˜ (η, τ )  = − A −1/2 (τ ) v˜ (η, τ )|η=0 . ∂ η η=0

(9)

˜ supposing that ∂∂ ηv˜ is known. It is the non-uniform (7) with a specific r.h.s. Consider now (8) as equation with respect to w Denoting for brevity



E (η, τ ) ≡ exp − A −1/2 (τ )η



118

I. Sofronov / Applied Numerical Mathematics 120 (2017) 115–124

we transform the r.h.s. as follows:



( B ∇τ ) + C 0

 ∂ v˜ ∂η

 ∂  = ( B ∇τ ) + C 0 E (η, τ ) v˜ 0 (τ ) ∂η = −( B ∇τ ) A −1/2 (τ ) E (η, τ ) v˜ 0 (τ ) − C 0 A −1/2 (τ ) E (η, τ ) v˜ 0 (τ )       = − B ∇τ A −1/2 E v˜ 0 − B A −1/2 ∇τ E v˜ 0 − B A −1/2 E ∇τ v˜ 0 − C 0 A −1/2 E v˜ 0      = − B ∇τ A −1/2 E + B A −1/2 E ∇τ + C 0 A −1/2 E v˜ 0 + O (η).

Evidently E (0, τ ) = I ; therefore, (8) reads at

A

η = 0:

˜ ∂2 w ˜ − D v˜ 0 = 0 −w ∂ η2

(10)

where









D ≡ B ∇τ A −1/2 + B A −1/2 ∇τ + C 0 A −1/2 . It follows from (7) and (10) that

∂2 ˜ + D v˜ 0 ) u˜ (η, τ , ε )|η=0 = A −1 v˜ + ε A −1 ( w ∂ η2

(11)

= A −1 u˜ 0 + ε A −1 D v˜ 0 = A −1 u˜ 0 + ε A −1 D u˜ 0 + o(ε )

where u˜ 0 ≡ u˜ (η, τ , ε )|η=0 . Now notice that (6) is also valid for the normal derivative ∂∂η u˜ (η, τ , ε ) since the coefficients of (6) don’t depend on Therefore, it follows from (5) that

∂ ∂2 u˜ (η, τ , ε )|η=0 = ( P 1 + ε P 0 ) u˜ (η, τ , ε )|η=0 + o(ε ) 2 ∂η ∂η

η.

(12)

= ( P 1 + ε P 0 )2 u˜ (η, τ , ε )|η=0 + o(ε ) = ( P 1 P 1 + ε P 1 P 0 + ε P 0 P 1 )u˜ (η, τ , ε )|η=0 + o(ε ). Using (11) and (12), we obtain

( P 1 P 1 + ε P 1 P 0 + ε P 0 P 1 )u˜ 0 = A −1 u˜ 0 + ε A −1 D u˜ 0 + o(ε ). Equating terms of zero and first degrees of ε , we derive the simple equation for finding P 1 :

P 1 P 1 = A −1 , and the Sylvester equation for finding P 0 :

P 1 P 0 + P 0 P 1 = A −1 D . Hence

P 1 = − A −1/2 ,

(13)

the negative sign follows from (9), and

A −1/2 P 0 + P 0 A −1/2 = − A −1











B ∇τ A −1/2 + B A −1/2 ∇τ + C 0 A −1/2 .

(14)

Coming back to the local coordinates (n, τ ) in (5) we get as s → ∞:

   ∂ uˆ  = (s P 1 + P 0 )uˆ |n=0 + o s0 ∂ n n=0

(15)

(we use here notation 1 = s0 to keep the variable in the asymptotically negligible expansion term). Neglecting the vanishing term we formulate the relationship on the boundary:

∂ uˆ = s P 1 uˆ + P 0 uˆ . ∂n

(16)

I. Sofronov / Applied Numerical Mathematics 120 (2017) 115–124

119

After the inverse Laplace transform of (16) we obtain the following approximate relationship of solution to (2) and its derivatives at the point O :

P1

∂u ∂u − + P 0u = 0 ∂t ∂n

(17)

with operators P 1 , P 0 , defined by (13) and the Sylvester equation (14). Existence of operators P 1 and P 0 in (17) follows from the hyperbolicity of (2). Indeed, introducing the characteristic determinant



2

Q (λ, ξn , ξ1 , ξ2 ) = det −λ I +

A ξn2

+

2 

B i ξi ξn +

i =1

2 2  



a i j ξi ξ j ,

i =1 j =1

see [9,5], where I is the unity N × N matrix, we have the hyperbolicity condition that the equation Q (λ, ξn , ξ1 , ξ2 ) = 0 should have 2N different real roots λ at the point O for arbitrary real ξn , ξ1 , ξ2 , |ξ | = 1. Suppose ξn = 1, ξ1 = ξ2 = 0. Then existence of 2N different real roots λ in the equation





det A − λ2 I = 0 means that all eigenvalues of matrix A are positive [16]. I.e., it has the reciprocal squared root A −1/2 . From positiveness of eigenvalues of matrix A −1/2 it also follows the resolution of the Sylvester equation (14) with respect to P 0 [15]. 2.3. Truncated TBCs Evidently, a further more accurate treatment of the last term in (15) will generate non-local in time operators after the inverse Laplace transform (a trivial case o(s0 ) ≡ 0 corresponds to local TBCs). Thus, (17) completely describes the differential part of the transparent boundary conditions, while any extra terms will have a non-local dependence on time. Also note that all known operators of TBCs for the second order hyperbolic equations [18–20,10,11,13] are the sum of the operator in (17) and a nonlocal operator (cf. (1)). The quasi-analytic TBCs, see details and notation in [21,23], have a more complex but similar form to analytic TBC,

Q−1 P1 Q

 ∂f ∂f − + Q−1 P0 Q f + Q−1 K˜ (t)∗ Q f = 0. ∂t ∂r

Equation (17) shows that one can use in this condition operators P 1 and P 0 instead of Q−1 P1 Q and Q−1 P0 Q, respectively. Therefore we suggest the equation (17) as an approximate local non-reflecting boundary condition for system (2), and call it Truncated TBCs (TTBC). The analysis made in this Section proves the following theorem Theorem. The differential part of the transparent boundary conditions for system (2) has the form

P1

∂u ∂u − + P 0u ∂t ∂n

where

P 1 = − A −1/2 and P 0 is determined by the Sylvester equation

A −1/2 P 0 + P 0 A −1/2 = − A −1











B ∇τ A −1/2 + B A −1/2 ∇τ + C 0 A −1/2 .

2.4. Evaluation of P 1 , P 0 Let us describe details of computations of matrices of TTBC operator. We are looking for operators P 1 and P 0 at the point O . The considered coefficients A, B, and C in (2) can depend on Therefore, let us introduce the matrices with numerical entries:

A O = A|O ,

B O = B|O ,

τ.

C O = C0|O ≡ C |O .

Then (13) is replaced by −1/2

P1 = −AO

.

(18)

Thus the matrix P 1 is found by usual algorithms for positively defined matrices. A more complex situation can arise for P 0 because it contains differential operators. Consider how to solve (14) for sufficiently smooth non-uniform coefficients.

120

I. Sofronov / Applied Numerical Mathematics 120 (2017) 115–124

We represent P 0 by sum of algebraic and differential terms:

P 0 = p + (q∇τ )

(19)

where matrices p and q = (q1 , q2 ) will be separately calculated. Introduce for brevity the following matrices: −1/2

L = AO

−1/2

1 M = − A− O CO AO

,

−1/2

1 i.e., ( K ∇τ ) = − A − O (B O A O

  1 −1/2  − A− , O B O ∇τ A O



−1/2 

1 K = − A− O BO AO

−1/2

1 = − A− O BO AO

,

∇τ ). Then (14) takes the form

 L P 0 + P 0 A −1/2  O = M + ( K ∇τ ). 

Using (19) we obtain





Lp + pL + L (q∇τ ) + q∇τ A −1/2  O = M + ( K ∇τ ). After collecting terms with and without the gradient operator ∇τ we come to the equations:

L (q∇τ ) + (qL ∇τ ) = ( K ∇τ )





Lp + pL = M − q∇τ A −1/2  O . The first one gives the following algebraic equation for q:

Lq + qL = K , and the second one – for p (with already known q). Coming back to the original notation we obtain two matrix equations for sequential resolving with respect to p and q: −1/2

AO

−1/2

AO

−1/2

q + qAO

−1/2

p + pAO

−1/2

1 = − A− O BO AO

−1

,

−1/2

= −AO CO AO

−1 

− A O B O ∇τ A

−1/2

    − q∇τ A −1/2  . O O

(20) (21)

Formulas (18), (19) with the Silvester equations (20), (21) describe the whole process of generating our TTBC operator for (2). Note that it is often convenient to use TTBC in the form resolved with respect to the time derivative:

∂u ∂u − P 1−1 + P 1−1 P 0 u = 0 ∂t ∂n or, denoting p 1 = P 1−1 , p 0 = P 1−1 P 0 ,

∂u ∂u − p1 + p 0 u = 0. ∂t ∂n

(22)

Then we straightforwardly have from (13), (14) the following equations for direct computation of p 1 , p 0 :

p 1 = − A 1/2 ,

     A −1/2 p 0 + p 0 A −1/2 = A −1/2 B ∇τ A −1/2 + B A −1/2 ∇τ + C 0 A −1/2 .

Finally let us note that implementation of TTBCs with a second-order of accuracy with respect to time and space is done immediately when coefficients of operators P 1 , P 0 are known. A reasonable way of deriving a stable explicit difference scheme of TTBCs equations for updating solution at boundary points consists in using already updated solution in the internal points of a grid. 3. Examples of generating TTBCs Two examples of deriving TTBCs according to the proposed method are given in [22] for 2D VTI cylindrical anisotropy and for 3D Navier wave equation. In [24] we reported some results on accuracy of these TTBCs for elastic equations. We reproduce here Fig. A1 from [24] (Fig. 2) where snapshots of the displacements vector absolute value on the central cross-section of a 3D computational domain at the sequential times t = 0.25, 0.45, 0.65, and 0.90 are shown (generated by a pointwise explosion source). One can clearly see in the bottom snapshots that reflections from the left boundary with TTBCs are about four times less compared with ones from the right boundary with the first order Clayton and Engquist absorbing boundary conditions [8]; see details in [24]. About well-posedness of problems with TTBCs. Leaving this important question outside the scope of the paper, we would like to note the following. Clearly, the TBCs guarantee the well-posedness by definition as they are derived from the original problem on unbounded domain. Our conjecture is that the absence of non-local terms in TTBCs influences their accuracy, but not the well-posedness. At any rate, the implementation and use of any specific TTBCs in numerical simulations should be supplemented by the corresponding numerical analysis of accuracy and stability. In particular, the continuation of computations demonstrated in Fig. 2 up to t = 100 shows no signs of instability. Let us consider another two examples of deriving TTBCs for elasticity and poroelasticity.

I. Sofronov / Applied Numerical Mathematics 120 (2017) 115–124

121

Fig. 2. Amplitude of absolute value of solution at times t = 0.25 and t = 0.45, (top panel), t = 0.65 and t = 0.9 (bottom panel). Much bigger reflections are visible from the right boundary (Clayton and Engquist ABCs) comparing to the left boundary (TTBCs). There are no boundaries at the top and bottom.

3.1. Three-dimensional orthotropic elasticity in cylindrical coordinates Here we derive TTBCs on the side cylinder boundary r = const for constant elastic coefficients (formulas at the boundaries z = const can be similarly derived). Let U = (u , v , w )T be the displacements vector in the cylindrical system of coordinates (r , θ, z). The uniform elasticity equations have the following form:

⎧ ∂ τrz σr − σθ ∂ 2 u ∂ σr 1 ∂ τr θ ⎪ ⎪ ⎪ + ρ = + + ⎪ 2 ⎪ ∂ r r ∂θ ∂z r ∂ t ⎪ ⎪ ⎨ 2 ∂ τθ z 2τr θ ∂ τr θ ∂ v 1 ∂ σθ + ρ = + + ⎪ ∂r r ∂θ ∂z r ⎪ ∂t2 ⎪ ⎪ ⎪ ⎪ ∂2 w ∂ τrz 1 ∂ τθ z ∂ σz τrz ⎪ ⎩ρ + = + + ∂r r ∂θ ∂z r ∂t2 where, according to the Hooke’s law for cylindrically orthotropic materials,

122

I. Sofronov / Applied Numerical Mathematics 120 (2017) 115–124

σr = A 11 εr + A 12 εθ + A 13 εz σθ = A 12 εr + A 22 εθ + A 23 εz σz = A 13 εr + A 23 εθ + A 33 εz τθ z = A 44 γθ z ,

τrz = A 55 γrz ,

τr θ = A 66 γr θ ,

ρ is the density, and A 11 ÷ A 66 are the stiffness coefficients; the strains are defined as usual ∂u ∂w 1 ∂v u , εθ = εz = + , ∂r r ∂θ r ∂z ∂v ∂ w ∂u v 1 ∂u γr θ = − + γrz = + , , ∂r r r ∂θ ∂r ∂z

εr =

γθ z =

∂v 1 ∂w + . ∂z r ∂θ

Following notation used in (2), we collect from above equations the terms containing derivatives with respect to r in the second order system:

−

∂2U ∂2U ∂2U ∂2U ∂U + A 2 + Bθ + Bz +C + · · · = 0, 2 ∂θ∂ r ∂ z∂ r ∂r ∂t ∂r

(23)

i.e., obtain the equations

⎧ ∂2 w ∂ 2u ∂ 2u 1 ∂2v 1 ∂u ⎪ ⎪ ⎪ + A 11 + ··· = 0 + ( A 13 + A 55 ) − ρ 2 + A 11 2 + ( A 12 + A 66 ) ⎪ ⎪ r ∂ r ∂θ ∂ r ∂ z r ∂r ∂ t ∂ r ⎪ ⎪ ⎨ ∂2v ∂2v 1 ∂v 1 ∂ 2u . − ρ + A + A + ( A + A ) + ··· = 0 66 66 12 66 ⎪ r ∂r r ∂ r ∂θ ∂t2 ∂ r2 ⎪ ⎪ ⎪ ⎪ ⎪ ∂ 2u ∂2 w ∂2 w 1 ∂w ⎪ ⎩−ρ + A 55 2 + A 55 + ( A 13 + A 55 ) + ··· = 0 2 r ∂r ∂ r∂ z ∂t ∂r

Thus the 3 × 3 matrices in (23) are as follows

⎛

A=

1

ρ

A 11 ⎝ 0 0

⎛

Bθ =

1

ρr

0 A 66 0

⎞

0 0 ⎠, A 55

A 11 ⎝ 0 C= ρr 0

⎝ A 12 + A 66 0

0 0 ⎠, 0

Bz =

Using (13), (14) we calculate the matrix

⎛

√ ⎜ ρ⎝

P 1 = − A −1/2 = −

0 A 66 0

⎞

A 12 + A 66 0 0

0

⎛

1

−1/2

A 11

0

0

0

A 66

0

0

0

−1/2

and obtain the equation for P 0 :

−1/2

ρ

⎛

0 0

⎝

A 13 + A 55

0 0 0

⎟ ⎠,

A 55

⎛

⎛

0 0

−1 −1 ⎝

A

A 13 + A 55

⎛

A 11 1 − ρ −1 A −1 ⎝ 0 r 0

A 12 + A 66 0 0

0 A 66 0

0 0 0

⎞

0 ∂ 0 ⎠ A −1/2 ∂θ 0

⎞

A 13 + A 55 ⎠ A −1/2 ∂ 0 ∂z 0

⎞

0 0 ⎠ A −1/2 . A 55

Due to the diagonal structure of A −1/2 the solution to (24) is easily evaluated:

⎛ ⎜ ⎜

P0 = − ⎜ ⎜

⎝

1 2r A 12 + A 66 1/2 1/2 A 11 A 66 + A 66

1 ∂ r ∂θ

A 13 + A 55 ∂ 1/2 1/2 A 11 A 55 + A 55 ∂ z

⎞

A 13 + A 55 ⎠. 0 0

⎞

0 1 A −1/2 P 0 + P 0 A −1/2 = − ρ −1 A −1 ⎝ A 12 + A 66 r 0

−ρ

1

⎞

0 0 ⎠, A 55

A 12 + A 66 1 ∂ 1/2 1/2 A 11 A 66 + A 11 r ∂θ

A 13 + A 55 ∂ 1/2 1/2 A 11 A 55 + A 11 ∂ z

1 2r

0

0

1 2r

⎞ ⎟ ⎟ ⎟. ⎟ ⎠

(24)

I. Sofronov / Applied Numerical Mathematics 120 (2017) 115–124

123

Therefore, using form (22), we have the following TTBC at r = const:

⎛

⎜ ∂U 1 ⎜ +√ ⎜ ∂t ρ⎜ ⎝

1 A 11 ( 2r + ∂∂r )

A 12 + A 66 1 ∂ 1/2 1/2 A 11 + A 66 r ∂θ

A 12 + A 66 1 ∂ 1/2 1/2 r ∂θ A +A

1 A 66 ( 2r + ∂∂r )

0

A 13 + A 55 ∂ 1/2 1/2 A 11 + A 55 ∂ z

0

1 A 55 ( 2r + ∂∂r )

1/2

11

66

A 13 + A 55 ∂ 1/2 1/2 A 11 + A 55 ∂ z

1/2

1/2

⎞ ⎟ ⎟ ⎟ U = 0. ⎟ ⎠

(25)

For VTI case we can take into account the constrain A 66 = ( A 11 − A 12 )/2. Note that other three VTI constrains ( A 22 = A 11 , A 23 = A 13 , A 55 = A 44 ) do not influence the condition (25). 3.2. Three-dimensional Biot’s poroelasticity equations in Cartesian coordinates We consider how to derive TTBCs at a plane boundary normal to one of the coordinate axes xi , i = 1, 2, 3. Without viscose terms (i.e. at η = 0) the equations (8.34) in [7] of waves propagation in anisotropic media in Cartesian coordinates (i , j = 1, 2, 3) read

   3 3  3 3 3  ∂2 ∂  ∂   μν A i j e μν + M i j ζ = 2 (ρ u i + ρ f w i ) − M μν e μν + M ζ ∂xj ∂ xi ∂t μ=1 ν =1 μ=1 ν =1 j =1   3  ∂2 = 2 ρ f ui + mi j w j . ∂t

(26)

j =1

Here (u 1 , u 2 , u 3 ) and w = ( w 1 , w 2 , w 3 ) are the vectors of solid matrix displacements and relative flow components of the fluid, respectively; ρ and ρ f are the mass densities of the porous material and fluid, respectively; mi j are the so-called efficient densities; e i j are components of the conventional strain tensor of the solid; ζ = − div w; μν A i j , M i j , M are the coefficients in the stress–strain matrix relations. μν

Suppose that all parameters do not depend on time, and the coefficients A i j , M i j , and M are continuously differentiable. Let n be a one of the coordinates (x1 , x2 , x3 ), and τ = (τ1 , τ2 ) be the vector of two remained ones. After some obvious algebra, the system (26) can be written in the form

−J

∂U ∂2U ∂2U  ∂U + C nn 2 + C nτ ∇τ + Cn + ··· = 0 2 ∂n ∂n ∂t ∂n

(27)

where U = (u 1 , u 2 , u 3 , w 1 , w 2 , w 3 )T is the unknown vector;

 ρI ρf I 

is the 6 × 6 matrix composed with the help of 3 × 3 matrices M f = {mi j } and I , the unity matrix; J= ρ I M f f μν nn nτ C , C , and C n are the 6 × 6 matrices with entries obtained from parameters A i j , M i j , M and their derivatives by 2

collecting equation coefficients at terms ∂∂n2 , (·∇τ ) ∂∂n , and ∂∂n , respectively; “. . . ” are the remained terms not used while generating TTBCs. We also suppose that (27) is a hyperbolic system, what is in accordance with the physics of the considered model describing waves’ propagation. Equation (27) differs from (2) by the matrix J at the time derivative. Nevertheless, the theory developed in Section 2 is also applicable to (27). After some algebra we have the following expression for operator P 1 :



P 1 = − J −1 C nn

−1/2

,

and equation for operator P 0 :



J −1 C nn

−1/2



P 0 + P 0 J −1 C nn

−1/2

 −1/2 1/2   −1  nτ  nn −1/2 1/2   nτ  nn −1/2 1/2  + C C . = − C nn C ∇τ C J J ∇τ + C n C nn J

The latter is straightforwardly solved according to formulas (19)–(21). Thus at each boundary point we can evaluate (generally speaking numerically rather than analytically) set of three 6 × 6 matrices for TTBCs in the form (17) or (22).

124

I. Sofronov / Applied Numerical Mathematics 120 (2017) 115–124

4. Conclusion We presented theory of deriving equations for computation of truncated TBCs for second order hyperbolic systems. A way of solving these equations in general case of non-uniform coefficients is developed. The obtained TTBCs are the local (differential) part of the entire TBCs operator. Examples of TTBCs for 3D cylindrically anisotropic elasticity and for inviscid poroelasticity (Biot equations) are considered. Acknowledgements The author is grateful to Schlumberger for permission to publish the work. The study was also supported in part by RSCF grant 15-11-00015. References [1] B. Alpert, L. Greengard, T. Hagstrom, Rapid evaluation of nonreflecting boundary kernels for time-domain wave propagation, SIAM J. Numer. Anal. 37 (2000) 1138–1164. [2] X. Antoine, H. Barucq, Microlocal diagonalization of strictly hyperbolic pseudodifferential systems and application to the design of radiation conditions in electromagnetism, J. Appl. Math. 61 (6) (2001) 1877–1905. [3] X. Antoine, et al., Absorbing boundary conditions for relativistic quantum mechanics equations, J. Comput. Phys. 277 (2014) 268–304. [4] D. Appelö, Absorbing layers and non-reflecting boundary conditions for wave propagation problems, Doctoral Thesis, Stockholm, Sweden, http://www.math.unm.edu/~appelo/preprints/appelothesis.pdf, 2005. [5] F. Bampi, C. Zordan, What are second order hyperbolic systems?, Z. Angew. Math. Phys. 44 (2) (1993) 290–305. [6] A.G. Benedict, S.E. Field, S.R. Lau, Fast evaluation of asymptotic waveforms from gravitational perturbations, Class. Quantum Gravity 30 (5) (2013) 055015. [7] M.A. Biot, Mechanics of deformation and acoustic propagation in porous media, J. Appl. Phys. 33 (4) (1962) 1482–1498. [8] R. Clayton, B. Engquist, Absorbing boundary conditions for acoustic and elastic wave equations, Bull. Seismol. Soc. Am. 6 (1977) 1529–1540. [9] R. Courant, D. Hilbert, Methods of Mathematical Physics, vol. 2: Partial Differential Equations, Wiley-Interscience, New York, 1989, p. 879. [10] M.J. Grote, J.B. Keller, Exact non-reflecting boundary conditions for the time dependent wave equation, SIAM J. Appl. Math. 55 (1995) 280–297. [11] M.J. Grote, J.B. Keller, Exact nonreflecting boundary condition for elastic waves, SIAM J. Appl. Math. 60 (3) (2000) 803–819. [12] M.J. Grote, C. Kirsch, Nonreflecting boundary condition for time-dependent multiple scattering, J. Comput. Phys. 221 (1) (2007) 41–62. [13] T. Hagstrom, Radiation boundary conditions for the numerical simulation of waves, Acta Numer. 8 (1999) 47–106. [14] H. Han, D. Yin, Absorbing boundary conditions for the multidimensional Klein–Gordon equation, Commun. Math. Sci. 5 (3) (2007) 743–764. [15] Kh.D. Ikramov, Nonsymmetric Eigenvalue Problem, Nauka, Moscow, 1991, in Russian. [16] A. Kurosh, Higher Algebra, Mir Publishers, Moscow, 1980, 428 p. [17] S.R. Lau, Rapid evaluation of radiation boundary kernels for time-domain wave propagation on blackholes: theory and numerical methods, J. Comput. Phys. 199 (2004) 376–422. [18] I.L. Sofronov, Conditions for complete transparency on the sphere for the three-dimensional wave equation, Russian Acad. Sci. Dokl. Math. 46 (1993) 397–401. [19] I.L. Sofronov, Non-reflecting inflow and outflow in wind tunnel for transonic time-accurate simulation, J. Math. Anal. Appl. 221 (1998) 92–115. [20] I.L. Sofronov, Artificial boundary conditions of absolute transparency for two- and three-dimensional external time-dependent scattering problems, Eur. J. Appl. Math. 9 (6) (1998) 561–588. [21] I.L. Sofronov, N.A. Zaitsev, Transparent boundary conditions for the elastic waves in anisotropic media, in: Serre Benzoni-Gavage (Ed.), Hyperbolic Problems: Theory, Numerics, Applications, Springer Verlag, 2008, pp. 997–1004. [22] I.L. Sofronov, Differential part of transparent boundary conditions for certain hyperbolic systems of second-order equations, Russian Acad. Sci. Dokl. Math. 79 (3) (2009) 412–414. [23] I.L. Sofronov, N.A. Zaitsev, Numerical generation of transparent boundary conditions on the side surface of a vertical transverse isotropic layer, J. Comput. Appl. Math. 234 (6) (2010) 1732–1738. [24] I. Sofronov, N. Zaitsev, L. Dovgilovich, Multi-block finite-difference method for 3D elastodynamic simulations in anisotropic subhorizontally layered media, Geophys. Prospect. 63 (5) (2015) 1142–1160. [25] I.L. Sofronov, L. Dovgilovich, N. Krasnov, Application of transparent boundary conditions to high-order finite-difference schemes for the wave equation in waveguides, Appl. Numer. Math. 93 (2015) 195–205. [26] L.L. Wang, B. Wang, X. Zhao, Fast and accurate computation of time-domain acoustic scattering problems with exact nonreflecting boundary conditions, SIAM J. Appl. Math. 72 (6) (2012) 1869–1898.