An accurate hyper-singular boundary integral equation method for dynamic poroelasticity in two dimensions
AAn accurate hyper-singular boundary integral equation methodfor dynamic poroelasticity in two dimensions
Lu Zhang ∗ , Liwei Xu † , Tao Yin ‡ August 18, 2020
Abstract
This paper is concerned with the boundary integral equation method for solving the exterior Neu-mann boundary value problem of dynamic poroelasticity in two dimensions. The main contribution ofthis work consists of two aspescts: the proposal of a novel regularized boundary integral equation, andthe presentation of new regularized formulations of the strongly-singular and hyper-singular bound-ary integral operators. Firstly, turning to the spectral properties of the double-layer operator andthe corresponding Calderón relation of the poroelasticity, we propose the novel low-GMRES-iterationintegral equation whose eigenvalues are bounded away from zero and infinity. Secondly, with the helpof the Günter derivatives, we reformulate the strongly-singular and hyper-singular integral operatorsinto combinations of the weakly-singular operators and the tangential derivatives. The accuracy andefficiency of the proposed methodology are demonstrated through several numerical examples.
Keywords:
Poroelasticity, hyper-singular operator, Calderón relation, regularized integral equa-tion
We investigate the application of the boundary integral equation method (BIEM) to solve the dy-namic poroelastic scattering problem [14,24–28] in an unbounded exterior domain, and this problem is ofgreat importance in many fields of applications such as oil and gas exploration, materials science, seismicanalysis, etc. The poroelastic problem can be characterized by the Biot model [8–13, 16], and the Neu-mann boundary condition will be considered in this work. Compared with the volumetric discretizationmethods [19, 20, 23, 29], the BIEM possesses such advantages as requiring discretization of domains oflower dimensionality and taking into account the radiation condition at infinity in a direct manner, andit has been widely studied for the numerical solutions of scattering problems [2, 3, 5–7, 15, 16, 18, 30].Regarding to the time-harmonic elastic wave scattering problems in an unbounded exterior domain,a combination form [2, 7, 17] of single-layer and double-layer potentials is usually used to represent thesolution, and the resulting combined boundary integral equation (CBIE) potentially ensures the validityof unique solvability corresponding to all frequencies. Although the unique solvability of the CBIE forthe poroelastic scattering problem remains unsolved, the CBIE still provides an efficient numerical toolfor the solution of the problem imposed on the unbounded domain. Since the integral operators resultingfrom the action of the traction operator on both the single-layer and the double-layer potentials containstrongly-singular and hyper-singular kernels, they are only well defined in the sense of Cauchy principlevalue and Hadamard finite part [21], respectively. Meanwhile, the appearance of the strongly-singular and ∗ School of Mathematical Sciences, University of Electronic Science and Technology of China, Chengdu, Sichuan 611731,China. Email: [email protected] † School of Mathematical Sciences, University of Electronic Science and Technology of China, Chengdu, Sichuan 611731,China. Email: [email protected] ‡ Department of Computing & Mathematical Sciences, California Institute of Technology, 1200 East California Blvd.,CA 91125, United States. Email: [email protected] a r X i v : . [ phy s i c s . c o m p - ph ] A ug yper-singular integral operators in the CBIE leads to some difficulties related to the spectral characterand the accurate evaluation of these operators. Firstly, it is known that the eigenvalues of the hyper-singular operators accumulate at infinity. Therefore, solving the CBIE by means of Krylov-subspaceiterative solvers, for instance the GMRES, often requires a relatively large number of iterations for theconvergence of numerical solution. Secondly, the evaluation of the associated Cauchy principle value (resp.Hadamard finite part) of the strongly-singular (resp. hyper-singular) integrals has also been remaining asignificant challenge.To reduce the number of GMRES iterations required in the process of solving the CBIE, an efficientmethodology which was originally proposed in [6,15] for acoustic and electromagnetic scattering problemsutilizes the Calderón relation together with a regularized operator with a form similar to a single-layeroperator. The derived regularized boundary integral equations (RBIEs) are of the second-kind Fredholmtype. This approach has been extended to the homogeneous elastic cases [5, 7], and it can be shown thatthe eigenvalues of the RBIE are bounded away from zero and infinity. However, the poroelastic Calderónformulas have not been studied in open literatures, and the main difficulty comes from the fact that theporoelastic double-layer integral operator (which plays important roles in the Calderón relations) is notcompact. On the basis of the spectral property studied in [1] for the static-elastic double-layer operator,it can be proved (see Theorem 3.1) that the poroelastic double-layer integral operator is polynomiallycompact, and then the composition of the single-layer and hyper-singular integral operators can beexpressed as the sum of a multiple of the identity operator and a compact operator. Relying on thetheoretical results, we propose a new RBIE method for solving the dynamic poroelastic problem, andverify numerically that the eigenvalues of the RBIE accumulate at fixed points only depending on theLamé parameters of elastic media.In this work, the classical Nyström method, which has been widely used for the acoustic and elasticproblems [17, 18, 22], is employed for the numerical implementation of the proposed RBIE. As applyingthe method, we encounter the challenge of accurate evaluation of the strongly-singular and hyper-singularintegrals. In light of the novel regularized formulations presented in [3,30] for the elastic and thermoelasticproblems, it can be shown (see Lemmas 3.2-3.6) that the strongly-singular and hyper-singular integralscan be re-expressed as compositions of weakly-singular integrals and tangential-derivative operators bymeans of the Günter derivative and integration-by-parts. As a result, the Nyström method allows usto evaluate the weakly-singular integrals with spectral accuracy for analytic surfaces, and to calculatethe tangential-derivative of a given function via fast Fourier transform (FFT) in GMRES iterations.Numerical tests show that the proposed scheme demonstrate a simpler and more efficient performancethan some alternative numerical treatments [18, 22].The remainder of this paper is organized as follows. Section 2 describes the dynamic poroelastic prob-lem (Section 2.1) and the classical combined field integral equation (Section 2.2). Section 3.1 theoreticallyand numerically studies the spectral properties of the poroelastic integral operators and the correspond-ing Calderón relation, and then a new regularized integral equation is proposed in Section 3.2. Exactreformulations of the strongly-singular and hyper-singular operators are presented in Section 3.3. TheNyström method for numerical evaluation of the integral operators is briefly described in Section 4. Sec-tion 5 presents the numerical examples to demonstrate the high-accuracy and efficiency of the proposedmethod. Finally, we present a conclusion in Section 6. Let Ω ⊂ R be a bounded domain with smooth boundary Γ := ∂ Ω . Assume that the exteriordomain Ω c = R \ Ω ⊂ R is occupied by a linear isotropic poroelastic medium. Following the Biot’stheory [8, 10, 11] to model wave propagation in poroelastic medium, the dynamic poroelastic problem infrequency-domain to be considered in this work is characterized by the governed equations of the solid2isplacements u = ( u , u ) (cid:62) and the pore pressure p that are given by ∆ ∗ u + ( ρ − βρ f ) ω u − ( α − β ) ∇ p = 0∆ p + qp + iωγ ∇ · u = 0 in Ω c , (2.1)or in an operator notation LU = 0 , L = (cid:20) ∆ ∗ + ( ρ − βρ f ) ω I − ( α − β ) ∇ iωγ ∇· ∆ + q (cid:21) , U = ( u , u , p ) (cid:62) , where β = ωφ ρ f κiφ + ωκ ( ρ a + φρ f ) , q = ω φ ρ f βR , γ = − iωρ f ( α − β ) β . Here, ω denotes the frequency, I is the identity operator and ∆ ∗ is the Lamé operator defined by ∆ ∗ := µ ∆ + ( λ + µ ) ∇∇· with ∆ being the Laplacian operator, and ∇ being the gradient operator. The material parameters usedin (2.1) are listed in Table 1. In addition, we consider the Neumann boundary condition on Γ given by (cid:101) T ( ∂, ν ) U := (cid:34) T ( ∂, ν ) − αν − iωβν (cid:62) iβωρ f ∂ ν (cid:35) U = F, (2.2)in which the traction operator T ( ∂, ν ) is defined as T ( ∂, ν ) u := 2 µ∂ ν u + λν ∇ · u + µν ⊥ ( ∂ u − ∂ u ) , ν ⊥ = ( − ν , ν ) (cid:62) , where ν = ( ν , ν ) (cid:62) denotes the outward unit normal to the boundary Γ and ∂ ν := ν · ∇ is the normalderivative. If the scattered field is induced by an incident field U inc , the boundary data is determined as F = − (cid:101) T ( ∂, ν ) U inc . Table 1: The material parameters in poroelasticity.Notation Physical meaning λ, µ ( µ > , λ + µ > Lamé parameters ν p Poisson ratio ν u undrained Poisson ratio B Skempton porepressure coefficient ρ s solid density ρ f fluid density ρ a apparent mass density φ porosity κ permeability coefficient ρ = (1 − φ ) ρ s + φρ f bulk density α = ν u − ν p ) B (1 − ν p )(1+ ν u ) compressibility R = φ µB (1 − ν p )(1+ ν u ) ν u − ν p )(1 − ν u ) constitutive coefficient Let E ( x, y ) be the fundamental solution of the adjoint operator L ∗ of L in R given by E ( x, y ) = (cid:20) E ( x, y ) E ( x, y ) E (cid:62) ( x, y ) E ( x, y ) (cid:21) , x (cid:54) = y, E ( x, y ) = 1 µ γ k s ( x, y ) I + 1( ρ − βρ f ) ω ∇ x ∇ (cid:62) x (cid:34) γ k s ( x, y ) − k p − k k − k γ k ( x, y ) + k p − k k − k γ k ( x, y ) (cid:35) ,E ( x, y ) = iωγ ( λ + 2 µ )( k − k ) ∇ x [ γ k ( x, y ) − γ k ( x, y )] ,E ( x, y ) = − γ ( λ + 2 µ )( k − k ) ∇ x [ γ k ( x, y ) − γ k ( x, y )] ,E ( x, y ) = iρ f ωβ ( k − k ) [( k p − k ) γ k ( x, y ) − ( k p − k ) γ k ( x, y )] , in which γ k t ( x, y ) = i H (1)0 ( k t | x − y | ) , x (cid:54) = y, t = s, p, , , denotes the fundamental solution of the Helmholtz equation in R with wave number k t . Here, k p and k s , referred as the compressional and shear wave numbers, respectively, are given by k p := ω (cid:115) ρ − βρ f λ + 2 µ , k s := ω (cid:115) ρ − βρ f µ . The wave numbers k , k , satisfying Im ( k i ) ≥ , i = 1 , , are the roots of the characteristic system k + k = q (1 + (cid:15) ) + k p , k k = qk p , (cid:15) = iωγ ( α − β ) q ( λ + 2 µ ) , and it follows that k = (cid:114) (cid:110) k p + q (1 + ε ) + (cid:113) [ k p + q (1 + ε )] − qk p (cid:111) ,k = (cid:114) (cid:110) k p + q (1 + ε ) − (cid:113) [ k p + q (1 + ε )] − qk p (cid:111) . From the potential theory, the unknown function U in Ω c can be represented as a combination of thesingle-layer and double-layer potentials U ( x ) = ( D − iηS )( ϕ )( x ) , x ∈ Ω c , Re ( η ) (cid:54) = 0 , (2.3)where S ( ϕ )( x ) := (cid:90) Γ ( E ( x, y )) (cid:62) ϕ ( y ) ds y , (2.4) D ( ϕ )( x ) := (cid:90) Γ ( (cid:101) T ∗ ( ∂ y , ν y ) E ( x, y )) (cid:62) ϕ ( y ) ds y , (2.5)denote the single-layer and double-layer potentials, respectively. Here (cid:101) T ∗ denotes the correspondingNeumann boundary operator of L ∗ given by (cid:101) T ∗ ( ∂, ν ) = (cid:34) T ( ∂, ν ) − iωαν − βν (cid:62) iβωρ f ∂ ν (cid:35) . (2.6)The combination form of solution representation (2.3) has been widely used for the corresponding acousticand elastic scattering problems [7, 17], and the resulting boundary integral equation can ensure unique4olvability for all frequencies. Operating with the boundary operator (cid:101) T on (2.3), taking the limit as x → Γ , the CBIE [ iη ( I − K (cid:48) ) + N ]( ϕ ) = F on Γ (2.7)is obtained. Here I denotes the identity operator and the boundary integral operators K (cid:48) and N aredefined by K (cid:48) ( ϕ )( x ) = (cid:101) T ( ∂ x , ν x ) (cid:90) Γ ( E ( x, y )) (cid:62) ϕ ( y ) ds y , (2.8)and N ( ϕ )( x ) = (cid:101) T ( ∂ x , ν x ) (cid:90) Γ (cid:16) (cid:101) T ∗ ( ∂ y , ν y ) E ( x, y ) (cid:17) (cid:62) ϕ ( y ) ds y , (2.9)in the sense of Cauchy principal value and Hadamard finite part [21], respectively, in view of the stronglysingular and hyper-singular character of the corresponding kernels. Remark 2.1.
It is well-known that the unkown solution U in Ω c can also be expressed simply as asingle-layer potential U ( x ) = S ( ψ )( x ) , x ∈ Ω c , (2.10) or as a double layer potential U ( x ) = D ( ψ )( x ) , x ∈ Ω c . (2.11) Operating with the boundary operator (cid:101) T on (2.10) and (2.11), taking the limit as x → Γ , we can obtainthe following boundary integral equations ( − I K (cid:48) )( ψ ) = F on Γ , (2.12) and N ( ψ ) = F on Γ , (2.13) for the unknown density ψ in (2.10) and (2.11), respectively. Remark 2.2.
Unfortunately, the unique solvability of integral equation (3.3) can not be derived followingthe classical approach to prove the corresponding unique solvability of combined field integral equationsfor acoustic and elastic problems [2, 17]. The main reason is that due to the special Neumann boundaryoperator (cid:101) T ( ∂, ν ) , there is no appropriate Green’s first identity for the considered poroelastic problems andonly the following Green’s second identity holds (cid:90) Ω ( LU · V − U · L ∗ V ) dx = (cid:90) Γ (cid:16) (cid:101) T ( ∂, ν ) U · V − U · (cid:101) T ∗ ( ∂, ν ) V (cid:17) ds. The uniqueness of integral equation (3.3) still remains open, however, as discussed in Section 5, the de-terminant of the stiffness matrix resulting from the discretization of (3.3) does not contain any significantshape trough which generally can indicate the existence of eigenfrequency [30].
Remark 2.3.
Note that the eigenvalues of the hypersingular integral operator N accumulate at infinity,see Figure 1 in which Γ = {| x | = 1 , x ∈ R } is considered. As a result, the solution of the integral equation(2.7) and (2.13) by means of Krylov-subspace iterative solvers such as GMRES generally requires largenumbers of iterations. N for a circular scatterer. To avoid the difficulty arising from the use of hypersingular integral operator in (2.7), we proposea new RBIE for solving the poroelastic problem in Section 2.1. Here, two types of “regularization” areemployed. By means of introducing a regularized operator R and studying the spectrum properties ofporoelastic integral operators (Section 3.1), we derive a new boundary integral equation which correspondsto a linear system with a better convergence property after the discretization compared to that of (2.7),see Section 3.2. Meanwhile, the strongly-singular and hyper-singular integral operators are re-expressedinto compositions of weakly-singular operators and differentiation operators in directions tangential tothe boundary using the Günter derivative and integration by parts (Section 3.3) and we call this as aregularization procedure. The spectra of the integral operators K (cid:48) is concluded in the following theorem. Theorem 3.1.
Let Γ denote a smooth closed surface in two-dimensional space. Then K (cid:48) − (cid:20) C λ,µ I
00 0 (cid:21) : H / (Γ) → H / (Γ) is compact, where C λ,µ is a constant given by C λ,µ = µ λ + 2 µ ) < . Furthermore, the spectrum of K (cid:48) consists of three nonempty sequences of eigenvalues which accumulateat , C λ,µ and − C λ,µ respectively.Proof. Recall the definition of K (cid:48) that K (cid:48) ( U )( x ) = (cid:101) T ( ∂ x , ν x ) (cid:90) Γ ( E ( x, y )) (cid:62) U ( y ) ds y = (cid:20) K (cid:48) K (cid:48) K (cid:48) K (cid:48) (cid:21) (cid:20) up (cid:21) ( x ) , x ∈ Γ , K (cid:48) j , j = 1 , · · · , are denoted as K (cid:48) ( u )( x ) = (cid:90) Γ (cid:0) T ( ∂ x , ν x ) E − αν x E (cid:62) (cid:1) u ( y ) ds y ,K (cid:48) ( p )( x ) = (cid:90) Γ ( T ( ∂ x , ν x ) E − αν x E ) p ( y ) ds y ,K (cid:48) ( u )( x ) = (cid:90) Γ (cid:18) − iωβν (cid:62) x E + iβρ f ω ∂ ν x E (cid:62) (cid:19) u ( y ) ds y ,K (cid:48) ( p )( x ) = (cid:90) Γ (cid:18) − iωβν (cid:62) x E + iβρ f ω ∂ ν x E (cid:19) p ( y ) ds y . Let K (cid:48) denote the static ( ω = 0) boundary integral operator corresponding to K (cid:48) K (cid:48) = (cid:101) T ( ∂ x , ν x ) (cid:90) Γ E ( x, y ) U ( y ) ds y = (cid:20) K (cid:48) , K (cid:48) , K (cid:48) , K (cid:48) , (cid:21) (cid:20) up (cid:21) ( x ) , x ∈ Γ where (cid:101) T ( ∂ x , ν x ) = (cid:20) T ( ∂ x , ν x ) − αν x κ∂ ν x (cid:21) and E ( x, y ) = (cid:20) E , E , E , E , (cid:21) = (cid:34) E e, − α ( x − y )2( λ +2 µ ) ln | x − y | − π ln | x − y | (cid:35) is the fundamental solution of static poroelastic problem with E e, ( x, y ) being the fundamental solutionof Lamé equation which is given by E e, ( x, y ) = λ + 3 µ πµ ( λ + 2 µ ) (cid:26) − ln | x − y | I + λ + µλ + 3 µ | x − y | ( x − y )( x − y ) (cid:62) (cid:27) . Thus we can obtain that K (cid:48) , ( u )( x ) = (cid:90) Γ T ( ∂ x , ν x ) E , u ( y ) ds y ,K (cid:48) , ( p )( x ) = (cid:90) Γ ( T ( ∂ x , ν x ) E , − αν x E , ) p ( y ) ds y ,K (cid:48) , ( u )( x ) = 0 ,K (cid:48) , ( p )( x ) = (cid:90) Γ κ∂ ν x E , p ( y ) ds y . From [1], it is known that K (cid:48) , − C λ,µ I is compact. It can be easily deduced that the kernels of K (cid:48) j, , j = 2 , are weakly-singualar implying that K (cid:48) j, , j = 2 , are compact. Therefore, K (cid:48) − (cid:20) C λ,µ I
00 0 (cid:21) = K (cid:48) ( K (cid:48) − K (cid:48) ) + ( K (cid:48) − K (cid:48) ) K (cid:48) + (cid:20) K (cid:48) , − C λ,µ I K (cid:48) , K (cid:48) , + K (cid:48) , K (cid:48) , K (cid:48) , (cid:21) is compact due to the fact that K (cid:48) − K (cid:48) has a weakly-singular kernel and is a compact operator. Theinequality < C λ,µ < / can be obtained easily from the conditions λ + µ > , µ > . This completesthe proof. 7rom Theorem 3.1, it can be seen that the accumulation points of the eigenvalues of K (cid:48) are independentof the frequency. In addition, we can conclude from the Calderón relation N S = − I K (cid:48) (3.1)that the spectrum of the composite operator N S , which plays an essential role in the regularized integralequations proposed in the following section, consists of two nonempty sequences of eigenvalues whichaccumulate at − and − + C λ,µ .In order to verify the above results numerically, we consider the problem of poroelastic scatteringby a circular scatterer of radius one, and choose the same values of coefficients as in Section 5 whichgives C λ,µ = 0 . . Figure 2 displays the eigenvalue distribution of the integral operators K (cid:48) and N S from which the eigenvalues of K (cid:48) and N S are seen to accumulate at the points predicted by ourtheoretical results. Here, the eigenvalue computation, on a basis of the exact regularized formulationsfor the strongly-singular operator K (cid:48) and the hypersingular operator N given in Section 3.3, has beenimplemented by means of the high-order Nyström methodology (see Section 4) together with FFT forevaluation of tangential derivatives and choosing a sufficiently large number of discretization points.(a) K (cid:48) (b) N S
Figure 2: Eigenvalue distribution of the operators K (cid:48) and N S for a circular scatterer.
Relying on the spectra studies of the poroelastic integral operators presented in Section 3.1, we proposein this section a RBIE by utilizing a regularization operator R in addition to the aforementioned double-layer and hypersingular operators K (cid:48) and N . In light of the regularized integral equation method discussedin [6] for acoustic problems, R is specified as the static ( ω = 0 ) single-layer operator S corresponding to S in this section (see the comparison between R = S and other selections of the regularization operator R in Section 5). Then replacing the solution presentation (2.3) by U ( x ) = ( D R − iηS )( ϕ )( x ) , x ∈ Ω c . (3.2)we obtain the following RBIE (cid:20) iη ( I − K (cid:48) ) + N R (cid:21) ( ϕ ) = F on Γ , (3.3)instead of the classical CBIE (2.7). It follows from the spectra results in Section 3.1 that the spectrum ofthe regularized combined field integral operator on the left hand side of (3.3) consists of three non-emptysequences of eigenvalues which converge to − / iη/ , − / C λ,µ + iη (1 / C λ,µ ) and − / C λ,µ + iη (1 / − C λ,µ ) , respectively, and these values are all bounded away from zero and infinity, see Figure 3.8igure 3: Eigenvalue distribution of the integral operator i ( I − K (cid:48) ) + N S for a circular scatterer.
As aforementioned, the integral operators K (cid:48) and N are strongly-singular and hyper-singular, respec-tively. In this section, following the techniques proposed in [30], we re-expresses these two operators ascombinations of weakly-singular integral operators and tangential derivative. (The main approach to de-rive the regularized formulations will be shown in Appendix.) Applying these new formulations togetherwith the Nyström method to be described in Section 4 and the linear algebra solver GMRES to equation(3.3) would lead then leads to the proposed solvers for the boundary value problem in Section 2.1.From [30], it is known that the traction operator T ( ∂, ν ) can be rewritten as T ( ∂, ν ) u ( x ) = ( λ + 2 µ ) ν ( ∇ · u ) + µ∂ ν u + µM ( ∂, ν ) u, where the operator M ( ∂, ν ) , whose elements are also called Günter derivatives [21], admits M ( ∂, ν ) u ( x ) = A duds , A = (cid:20) −
11 0 (cid:21) . For the strongly-singular integral operator K (cid:48) , the following regularized formulation can be obtainedusing the notations in Theorem 3.1. Lemma 3.2.
The boundary integral operators K (cid:48) j , j = 1 , · · · , can be expressed as K (cid:48) j = K (cid:48) j + dds K (cid:48) j , j = 1 , · · · , , (3.4)9 here K (cid:48) ( u )( x ) = (cid:90) Γ ν x ∇ (cid:62) x (cid:20) − ( γ k s ( x, y ) − γ k ( x, y )) + k − qk − k ( γ k s ( x, y ) − γ k ( x, y )) (cid:21) u ( y ) ds y + (cid:90) Γ (cid:2) ∂ ν x γ k s ( x, y ) I − αν x E (cid:62) ( x, y ) (cid:3) u ( y ) ds y ,K (cid:48) ( u )( x ) = A (cid:90) Γ [2 µE ( x, y ) − γ k s ( x, y ) I ] u ( y ) ds y ,K (cid:48) ( p )( x ) = (cid:90) Γ (cid:20) γk − k ( k γ k ( x, y ) − k γ k ( x, y )) ν x − αν x E ( x, y ) (cid:21) p ( y ) ds y ,K (cid:48) ( p )( x ) = − µγ ( λ + 2 µ )( k − k ) A (cid:90) Γ ∇ x [ γ k ( x, y ) − γ k ( x, y )] p ( y ) ds y ,K (cid:48) ( u )( x ) = (cid:90) Γ (cid:2) − iωβν (cid:62) x E ( x, y ) + C ( k γ k ( x, y ) − k γ k ( x, y )) ν (cid:62) x (cid:3) u ( y ) ds y ,K (cid:48) ( u )( x ) = C (cid:90) Γ ∇ (cid:62) x [ γ k ( x, y ) − γ k ( x, y )] Ap ( y ) ds y ,K (cid:48) ( p )( x ) = (cid:90) Γ ν (cid:62) x (cid:34) iωβE ( x, y ) + ∇ x (cid:32) k p − k k − k γ k ( x, y ) − k p − k k − k γ k ( x, y ) (cid:33)(cid:35) p ( y ) ds y ,K (cid:48) ( p )( x ) = 0 , with the constant C being C = βγρ f ( λ + 2 µ )( k − k ) . Now we consider the hyper-singular operator N . For ψ = ( u (cid:62) , p ) (cid:62) , denote N ( ψ )( x ) = (cid:20) N N N N (cid:21) (cid:20) up (cid:21) ( x ) , x ∈ Γ . The regularized formulations for the operators N j , j = 1 , · · · , are given in the following Lemmas 3.3-3.6. Lemma 3.3.
The hyper-singular operator N can be expressed as N = N + dds N dds + dds N + N dds , (3.5) where N ( u )( x ) = − ( ρ − βρ f ) ω (cid:90) Γ γ k s ( x, y ) (cid:0) ν x ν (cid:62) y − ν (cid:62) x ν y I − J ν x ,ν y (cid:1) u ( y ) ds y , + (cid:90) Γ [ C γ k ( x, y ) − C γ k ( x, y )] ν x ν (cid:62) y u ( y ) ds y N ( u )( x ) = 4 µ (cid:90) Γ [ γ k s ( x, y ) + AE ( x, y ) A ] u ( y ) ds y ,N ( u )( x ) = (cid:90) Γ A ∇ y [ − µ ( γ k s ( x, y ) − γ k ( x, y )) + C ( γ k ( x, y ) − γ k ( x, y ))] ν (cid:62) y u ( y ) ds y ,N ( u )( x ) = (cid:90) Γ ν x ∇ (cid:62) x [ − µ ( γ k s ( x, y ) − γ k ( x, y )) + C ( γ k ( x, y ) − γ k ( x, y ))] Au ( y ) ds y , with the constant C i , i = 1 , , being C = k ( k − q )( λ + 2 µ ) − iωαγk ( k − k ) , C = k ( k − q )( λ + 2 µ ) − iωαγk ( k − k ) = 2 µk − k ( − iωγαλ + 2 µ + k − q ) . Lemma 3.4.
The hyper-singular operator N can be expressed as N = N + dds N dds + dds N , (3.6) where N ( p )( x ) = β (cid:90) Γ ν x ∇ (cid:62) x (cid:20) ( γ k s ( x, y ) − γ k ( x, y )) − k − qk − k ( γ k ( x, y ) − γ k ( x, y )) (cid:21) ν y p ( y ) ds y − β (cid:90) Γ (cid:2) ∂ ν x γ k s ( x, y ) − αν x E (cid:62) ( x, y ) ν y (cid:3) p ( y ) ds y + iβγρ f ω ( k − k ) (cid:90) Γ (cid:2) ∂ ν y ( k γ k ( x, y ) − k γ k ( x, y )) (cid:3) ν x p ( y ) ds y + αk − k (cid:90) Γ ∂ ν y (cid:2) ( k p − k ) k γ k ( x, y ) − ( k p − k ) k γ k ( x, y ) (cid:3) ν x p ( y ) ds y ,N ( p )( x ) = − iµβγρ f ω ( λ + 2 µ )( k − k ) (cid:90) Γ ∇ x [ γ k ( x, y ) − γ k ( x, y )] ds y ,N ( p )( x ) = − β (cid:90) Γ A [2 µE ( x, y ) − γ k s ( x, y ) I )] ν y p ( y ) ds y − iµβγρ f ω ( λ + 2 µ )( k − k ) (cid:90) Γ (cid:2) k γ k ( x, y ) − k γ k ( x, y ) (cid:3) Aν y p ( y ) ds y . Lemma 3.5.
The hyper-singular operator N can be expressed as N = N + dds N dds + N dds , (3.7) where N ( u )( x ) = − iωβ (cid:90) Γ ν (cid:62) x ∇ x (cid:20) γ k s ( x, y ) − γ k ( x, y ) − k − qk − k ( γ k ( x, y ) − γ k ( x, y )) (cid:21) ν (cid:62) y u ( y ) ds y − iωβ (cid:90) Γ (cid:2) ∂ ν y γ k s ( x, y ) ν (cid:62) x − iωαν (cid:62) x E ν (cid:62) y (cid:3) u ( y ) ds y − βγρ f ( k − k ) (cid:90) Γ ∂ ν x (cid:2) k γ k ( x, y ) − k γ k ( x, y ) (cid:3) ν (cid:62) y u ( y ) ds y + iωαk − k (cid:90) Γ ∂ ν x (cid:2) ( k p − k ) γ k ( x, y ) − ( k p − k ) γ k ( x, y ) (cid:3) ν (cid:62) y u ( y ) ds y ,N ( u )( x ) = − µβγρ f ( λ + 2 µ )( k − k ) (cid:90) Γ ∇ (cid:62) x [ γ k ( x, y ) − γ k ( x, y )] u ( y ) ds y N ( u )( x ) = − iωβ (cid:90) Γ ν (cid:62) x [2 µE ( x, y ) − γ k s ( x, y ) I ] Au ( y ) ds y + 2 µβγρ f ( λ + 2 µ )( k − k ) (cid:90) Γ (cid:2) k γ k ( x, y ) − k γ k ( x, y ) (cid:3) ν (cid:62) x Au ( y ) ds y . Lemma 3.6.
The hyper-singular operator N can be expressed as N = N + dds N dds + dds N , (3.8)11 here N ( p )( x ) = iωβ (cid:90) Γ ν (cid:62) x E ( x, y ) p ( y ) ds y + 2 iβγω ( λ + 2 µ ) C (cid:90) Γ (cid:2) k γ k ( x, y ) − k γ k ( x, y ) (cid:3) ν (cid:62) x ν y p ( y ) ds y N ( p )( x ) = iβγω ( λ + 2 µ ) C (cid:90) Γ [ γ k ( x, y ) − γ k ( x, y )] p ( y ) ds y − C (cid:90) Γ (cid:2) ( k p − k ) γ k ( x, y ) − ( k p − k ) γ k ( x, y ) (cid:3) p ( y ) ds y N ( p )( x ) = iβγω ( λ + 2 µ ) C (cid:90) Γ ∇ (cid:62) x [ γ k ( x, y ) − γ k ( x, y )] Aν y p ( y ) ds y − C (cid:90) Γ (cid:2) ( k p − k ) k γ k ( x, y ) − ( k p − k ) k γ k ( x, y ) (cid:3) ν (cid:62) x ν y p ( y ) ds y , with C = iβρ f ω ( k − k ) . Remark 3.7.
It can be verified that the operators K (cid:48) ij , i = 1 , , j = 1 , · · · , , N ij , i = 1 , · · · , , j = 1 , · · · , and N are all weakly-singular for smooth boundary Γ . The derived results of the regularized formulationfor the integral operators K (cid:48) and N can also be extended to the case of Lipshitz boundary in terms of theproperties of Günter derivative given in [4]. According to the regularized formulations of integral operators K (cid:48) and N given in Section 3.3, thenumerical implementation of the RBIE (3.3) can be converted into the evaluation of multiple operatorsof two types, (i) Integral operators of the form H ( ϕ )( x ) = (cid:90) Γ H ( x, y ) ϕ ( y ) ds y (4.1)in which the kernel H ( x, y ) is weakly-singular, and (ii) Tangential derivative d/ds of a given smoothfunction defined on Γ . This section presents algorithms for numerical evaluation of operators of thesetypes by utilizing the well-known Nyström method [17].Assume that the boundary curve Γ is analytic and is given through Γ = { x ( t ) = ( x ( t ) , x ( t )) : 0 ≤ t ≤ π } (4.2)in counterclockwise orientation where x : R → R is analytic and 2 π -periodic with | x (cid:48) ( t ) | > for all t .Then the outward unit normal at x ∈ Γ is given by ν t = ν x ( t ) = ( x (cid:48) ( t ) , − x (cid:48) ( t )) | x (cid:48) ( t ) | and the integral (4.1) canbe transformed into a parametric form H ( ϕ )( x ( t )) = (cid:90) π (cid:101) H ( t, τ ) ψ ( τ ) dτ, (4.3)where (cid:101) H ( t, τ ) = H ( x ( t ) , x ( τ )) | x (cid:48) ( τ ) | and ψ ( t ) := ϕ ( x ( t )) . For the weakly-singular kernel (cid:101) H ( t, τ ) , theNyström method requires splitting it into (cid:101) H ( t, τ ) = (cid:101) H ( t, τ ) ln(4 sin t − τ (cid:101) H ( t, τ ) , (4.4)12here the terms (cid:101) H , (cid:101) H are analytic. Choosing an equidistant mesh t j := jπn , j = 0 , · · · , N − N ∈ N ) ,the quantities H ( ϕ )( x ( t i )) , i = 0 , · · · , N − can be approximated with exponential decaying errors usingthe weighted trigonometric interpolation quadratures and the trapezoidal rule, that is, H ( ϕ )( x ( t i )) ≈ N − (cid:88) j =0 R ( n ) j ( t i ) (cid:101) H ( t i , t j ) ψ ( t j ) + πN N − (cid:88) j =0 H ( t i , t j ) ψ ( t j ) . For the regularized integral equation (3.3), it can be summarized from the regularized formulations derivedin Section 3.3 that the weakly-singular kernel (cid:101) H ( t, τ ) only takes four types of forms that are listed asfollows, in which c e = 0 . · · · denotes the Euler constant, together with the corresponding splittingterms (cid:101) H ( t, τ ) , (cid:101) H ( t, τ ) . • Type 1. (cid:101) H ( t, τ ) = H (1)0 ( k | x ( t ) − x ( τ ) | ) | x (cid:48) ( τ ) | , (cid:101) H ( t, τ ) = (cid:40) iπ J ( k | x ( t ) − x ( τ ) | ) | x (cid:48) ( τ ) | , t (cid:54) = τ, iπ | x (cid:48) ( t ) | , t = τ, (cid:101) H ( t, τ ) = (cid:40) (cid:101) H ( t, τ ) − (cid:101) H ( t, τ ) ln(4sin ( t − τ )) , t (cid:54) = τ, (1 + iπ ( c e + ln k | x ( τ ) (cid:48) | )) | x (cid:48) ( t ) | , t = τ. • Type 2. (cid:101) H ( t, τ ) = ˆ kH (1)1 (ˆ k | x ( t ) − x ( τ ) | ) − ˜ kH (1)1 (˜ k | x ( t ) − x ( τ ) | ) | x ( t ) − x ( τ ) | | x (cid:48) ( τ ) | , (cid:101) H ( t, τ ) = (cid:40) iπ ˆ kJ (ˆ k | x ( t ) − x ( τ ) | ) − ˜ kJ (˜ k | x ( t ) − x ( τ ) | ) | x ( t ) − x ( τ ) | | x (cid:48) ( τ ) | , t (cid:54) = τ, π ( k − k ) | x (cid:48) ( t ) | , t = τ, (cid:101) H ( t, τ ) = (cid:40) (cid:101) H ( t, τ ) − (cid:101) H ( t, τ ) ln(4sin ( t − τ )) , t (cid:54) = τ, (cid:104) k − k (1 + iπ (ln | x (cid:48) ( t ) | + c e ) − iπ ) + iπ (cid:0) k ln k − k ln k (cid:1)(cid:105) | x (cid:48) ( t ) | , t = τ. • Type 3. (cid:101) H ( t, τ ) = ˆ k H (1)2 (ˆ k | x ( t ) − x ( τ ) | ) − ˜ k H (1)2 (˜ k | x ( t ) − x ( τ ) | ) | x ( t ) − x ( τ ) | | x (cid:48) ( τ ) | , (cid:101) H ( t, τ ) = (cid:40) iπ ˆ k J (ˆ k | x ( t ) − x ( τ ) | ) − ˜ k J (˜ k | x ( t ) − x ( τ ) | ) | x ( t ) − x ( τ ) | | x (cid:48) ( τ ) | , t (cid:54) = τ , t = τ (cid:101) H ( t, τ ) = (cid:40) (cid:101) H ( t, τ ) − (cid:101) H ( t, τ ) ln(4sin ( t − τ )) , t (cid:54) = τ, − iπ ( k − k ) | x (cid:48) ( t ) | , t = τ. • Type 4. (cid:101) H ( t, τ ) = ik ν t · [ x ( τ ) − x ( t )] H (1)1 ( k | x ( t ) − x ( τ ) | ) | x ( t ) − x ( τ ) | | x (cid:48) ( τ ) || x (cid:48) ( t ) | , (cid:101) H ( t, τ ) = (cid:40) − k π ν t · [ x ( τ ) − x ( t )] J ( k | x ( t ) − x ( τ ) | ) | x ( t ) − x ( τ ) | | x (cid:48) ( τ ) | | x (cid:48) ( t ) | , t (cid:54) = τ, , t = τ, (cid:101) H ( t, τ ) = (cid:40) (cid:101) H ( t, τ ) − (cid:101) H ( t, τ ) ln(4sin ( t − τ )) , t (cid:54) = τ, π ν t · x (cid:48)(cid:48) ( t ) | x (cid:48) ( t ) | , t = τ. ϕ ( x ( t )) , t ∈ [0 , π ] with values ϕ ( x ( t j )) , j = 0 , · · · , N − ,applying the Fourier series approximation ϕ ( x ( t )) = N − (cid:88) m = − N ϕ m e imt , t ∈ [0 , π ] , which can be obtained using FFT, we can compute the tangential derivative of ϕ at x ( t j ) , j = 0 , · · · , N − as dϕds (cid:12)(cid:12)(cid:12) x = x ( t j ) = 1 | x (cid:48) ( t j ) | N − (cid:88) m = − N imϕ m e imt j which can also be evaluated by FFT. In this section, various numerical examples are presented to demonstrate the accuracy and efficiencyof the proposed RBIE method for solving the poroelastic problems in two dimensions. We always choosethe values of parameters in (2.1) as µ = 1 , ν p = 0 . , ν u = 0 . , B = 0 . , φ = 0 . , ρ f = 1 , ρ s = 3 , ρ a = 0 . , κ = 0 . . The fully complex version of the iterative solver GMRES would be utilized toproduce the solutions of the integral equations. The relative maximum errors presented in this sectionare calculated in accordance with the expression ε ∞ := max x ∈ S {| U num ( x ) − U exa ( x ) |} max x ∈ S {| U exa ( x ) |} , (5.1)where S = {| x | = 6 , x ∈ R } ⊂ Ω c and U exa is the exact solution. All of the numerical tests are obtainedusing MATLAB.Table 2: Numerical errors of solutions for the problem of poroelastic scattering by a kite-shaped orrounded-triangle-shaped obstacle. GMRES tol: − . ω N Boundary (5.2) Boundary (5.3)CBIE RBIE CBIE RBIE20 . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × −
10 150 . × − . × − . × − . × − . × − . × − . × − . × − In our first example, we consider a kite-shaped or rounded-triangle-shaped obstacle with boundary Γ characterized by x ( t ) = (2 + 0 . t )(cos t, sin t ) , t ∈ [0 , π ) , (5.2)and x ( t ) = (cos t + 0 .
65 cos 2 t − . , . t ) , t ∈ [0 , π ) , (5.3)respectively, and the exact solution is given by u ( x ) = E ( x, z ) , p ( x ) = E ( x, z ) , x ∈ Ω c , ( u ) (b) Re ( u ) (c) Re ( p ) (d) Im ( u ) (e) Im ( u ) (f) Im ( p ) Figure 4: Comparison of the exact and numerical solutions with N = 256 for ω = 20 , and the relativeerror is . × − .with z = (0 , . (cid:62) ∈ Ω . Table 2 displays the numerical errors of solutions using the un-regularized integralequation (2.7) or the regularized integral equation (3.3) with respect to N , which demonstrate the highaccuracy and rapid convergence of newly proposed method. The numerical and exact solutions on S forthe cases of kite-shaped obstacle, N = 256 , ω = 20 are presented in Figure 4 and the relative maximumerror is . × − .Table 3: Iterations required for the problem of poroelastic scattering by a kite-shaped or rounded-triangle-shaped obstacle. GMRES tol: − . ω N Boundary (5.2) Boundary (5.3)CBIE RBIE CBIE RBIE1 20 38 22 48 175 100 132 52 242 3910 200 285 89 468 6620 400 590 175 781 12630 600 893 234 1005 171Table 3 presents the iteration numbers required to achieve the GMRES tolerance − in solvingthe boundary integral equations (2.7) and (3.3). It can be seen that the regularized integral equation(3.3) requires much smaller number of iterations than the un-regularized integral equation (2.7). It ismentioned in Remark 2.1 that the single-layer potential (2.4) can be used to represent the solution andthe integral equation (2.12) is obtained for the unknown potential. In addition, according to the Caldrónrelation studied in Section 3.1, the regularized operator R can also be selected as the single-layer operators S , and the resulted regularized boundary integral equation enjoys the same spectral property satisfied by(3.3). Table 4 lists the iteration numbers required by various integral equations to achieve the GMREStolerance − and it shows that the regularized integral equations (3.3) with R = S or R = S requirealmost the same number of iterations which is smaller than that required by equation (2.12). In Figure 5,we plot the determinants of the stiffness matrixes resulting from the discretization of the single-layeroperator S and the regularized integral equation with R = S or R = S . It can be seen that significanttroughs appear at the same frequencies in Figure 5(a,b) which generally shows the existence of interior15able 4: Iterations required by various integral equations for the problem of poroelastic scattering by akite-shaped obstacle. GMRES tol: − . ω N Eqn. (2.12) RBIE with R = S RBIE with R = S R = S can ensures the unique solvability of (3.3)for all frequencies from numerical point of view to a certain context.(a) S (b) Eqn. (3.3) with R = S (c) Eqn. (3.3) with R = S Figure 5: The determinant of the stiffness matrix with respect to frequency.Finally, we consider the scattering of an incident point source U inc = ( u inc (cid:62) , p inc ) (cid:62) in the form u inc ( x ) = E ( x, z ) , p inc ( x ) = E ( x, z ) by a kite-shaped obstacle where z denotes the location of the point source. The numerical solutions in Ω c with ω = 20 are presented in Figures 6 for z = (3 , and z = ( − , , respectively. We propose in this work a novel regularized integral equation method for solving the exterior Neumannboundary value problem of poroelastic wave scattering in two dimensions. On a basis of the spectralproperties of poroelastic double-layer integral operator and the corresponding Caldrón relation, and thenew regularized formulations of strongly-singular and hyper-singular operators, the proposed approachpossesses spectral convergence (using Nyström discretization method) and reduces the number of GMRESiterations consistently across various geometries and frequency regimes. Theoretical analysis of the uniquesolvability of combined integral equation (2.7), the combination of the regularized integral equationmethod with other popular fast solvers such as fast multipole method, the application of the regularizedintegral equation method for three dimensional problems and solving problems with non-smooth surfacesare left for future work.
Acknowledgement
LWX is partially supported by Key Project of the Major Research Plan of NSFC (91630205) and agrant of NSFC (11771068 ). 16a) | u | (b) | p | (d) | u | (e) | p | Figure 6: The total field U for the scattering of an incident point source by a kite-shaped obstacle. (a,b): z = ( − , ; (c,d): z = (3 , . Appendix A. Proofs of Lemmas 3.2-3.6
This appendix presents the main idea for the proofs of Lemmas 3.2-3.8. Recall that the tractionoperator can be rewritten as T ( ∂, ν ) u ( x ) = ( λ + 2 µ ) ν ( ∇ · u ) + µ∂ ν u + µM ( ∂, ν ) u. (A.1)It follows that ∂ ν ∇ − M ( ∂, ν ) ∇ = ν ∆ , (A.2)and T ( ∂, ν ) ∇ = ( λ + 2 µ ) ν ∆ + 2 µM ( ∂, ν ) ∇ . (A.3) Proof of Lemma 3.2.
We first investigate the operator K (cid:48) . It is known from [3] that for x (cid:54) = y , T ( ∂ x , ν x ) E ( x, y ) = − ν x ∇ (cid:62) x [ γ k s ( x, y ) − γ k ( x, y )] + k − qk − k ν x ∇ (cid:62) x [ γ k ( x, y ) − γ k ( x, y )] (A.4) + ∂ ν x γ k s ( x, y ) I + M ( ∂ x , ν x ) [2 µE ( x, y ) − γ k s ( x, y ) I ] , and T ( ∂ y , ν y ) E ( x, y ) = − ν y ∇ (cid:62) y [ γ k s ( x, y ) − γ k ( x, y )] + k − qk − k ν y ∇ (cid:62) y [ γ k ( x, y ) − γ k ( x, y )] (A.5) + ∂ ν y γ k s ( x, y ) I + M ( ∂ y , ν y ) [2 µE ( x, y ) − γ k s ( x, y ) I ] . K (cid:48) ( u )( x ) = (cid:90) Γ ν x ∇ (cid:62) x (cid:20) − ( γ k s ( x, y ) − γ k ( x, y )) + k − qk − k ( γ k s ( x, y ) − γ k ( x, y )) (cid:21) u ( y ) ds y + (cid:90) Γ (cid:20) ∂ ν x γ k s ( x, y ) I − αν x E (cid:62) ( x, y ) + dds x A (2 µE ( x, y ) − γ k s ( x, y ) I ) (cid:21) u ( y ) ds y . It follows from (A.3) that T ( ∂ x , ν x ) E ( x, y ) = − γ ( λ + 2 µ )( k − k ) T ( ∂ x , ν x ) ∇ x [ γ k ( x, y ) − γ k ( x, y )]= − γ ( k − k ) ν x ∆ x [ γ k ( x, y ) − γ k ( x, y )] − µγ ( λ + 2 µ )( k − k ) M ( ∂ x , ν x ) ∇ x [ γ k ( x, y ) − γ k ( x, y )]= γ ( k − k ) (cid:2) k γ k ( x, y ) − k γ k ( x, y ) (cid:3) ν x − µγ ( λ + 2 µ )( k − k ) A dds x ∇ x [ γ k ( x, y ) − γ k ( x, y )] . Therefore, K (cid:48) ( p )( x ) = (cid:90) Γ (cid:20) γk − k ( k γ k ( x, y ) − k γ k ( x, y )) − αν x E ( x, y ) (cid:21) p ( y ) ds y − (cid:90) Γ µγ ( λ + 2 µ )( k − k ) A dds x ∇ x [ γ k ( x, y ) − γ k ( x, y )] p ( y ) ds y . On the other hand, we obtain from (A.2) that ∂ ν x E (cid:62) ( x, y ) = iωγ ( λ + 2 µ )( k − k ) ∂ ν x ∇ (cid:62) x [ γ k ( x, y ) − γ k ( x, y )]= iωγ ( λ + 2 µ )( k − k ) ν x ∆ x [ γ k ( x, y ) − γ k ( x, y )]+ iωγ ( λ + 2 µ )( k − k ) { M ( ∂ x , ν x ) ∇ x [ γ k ( x, y ) − γ k ( x, y )] } (cid:62) = − iωγ ( λ + 2 µ )( k − k ) (cid:2) k γ k ( x, y ) − k γ k ( x, y ) (cid:3) ν x − iωγ ( λ + 2 µ )( k − k ) dds x ∇ (cid:62) x [ γ k ( x, y ) − γ k ( x, y )] A. Then we have K (cid:48) ( u )( x ) = (cid:90) Γ (cid:2) − iωβν (cid:62) x E ( x, y ) + C ( k γ k ( x, y ) − k γ k ( x, y )) ν (cid:62) x (cid:3) u ( y ) ds y + (cid:90) Γ C dds x ∇ (cid:62) x [ γ k ( x, y ) − γ k ( x, y )] Ap ( y ) ds y . The formula for K (cid:48) can be obtained directly from its definition and this completes the proof of Lemma 3.2. Proof of Lemmas 3.3-3.6.
Note that the hyper-singular integral operator of N is given by N ( ψ )( x ) = (cid:20) N N N N (cid:21) (cid:20) up (cid:21) ( x ) , x ∈ Γ , ψ = ( u (cid:62) , p ) (cid:62) , N ( u )( x ) = (cid:90) Γ (cid:2) T ( ∂ x , ν x )( T ( ∂ y , ν y ) E ( x, y )) (cid:62) − iωαT ( ∂ x , ν x ) E ( x, y ) ν (cid:62) y (cid:3) u ( y ) ds y + (cid:90) Γ (cid:2) − αν x ( T ( ∂ y , ν y ) E ( x, y )) (cid:62) + iωα E ( x, y ) ν x ν (cid:62) y (cid:3) u ( y ) ds y N ( p )( x ) = (cid:90) Γ (cid:20) − βT ( ∂ x , ν x ) E ( x, y ) ν y + iβρ f ω ∂ ν y T ( ∂ x , ν x ) E ( x, y ) (cid:21) p ( y ) ds y + (cid:90) Γ (cid:20) αβν x E (cid:62) ν y − iαβρ f ω ∂ ν y E ( x, y ) ν x (cid:21) p ( y ) ds y N ( u )( x ) = − (cid:90) Γ (cid:2) iωβν (cid:62) x ( T ( ∂ y , ν y ) E ( x, y )) (cid:62) + ω αβν (cid:62) x E ν (cid:62) y (cid:3) u ( y ) ds y + (cid:90) Γ (cid:20) iβρ f ω ∂ ν x ( T ( ∂ y , ν y ) E ( x, y )) (cid:62) + αβρ f ∂ ν x E ( x, y ) ν (cid:62) y (cid:21) u ( y ) ds y N ( p )( x ) = (cid:90) Γ (cid:20) iωβ ν (cid:62) x E ( x, y ) ν y + β ρ f ν (cid:62) x ∂ ν y E ( x, y ) (cid:21) p ( y ) ds y − (cid:90) Γ (cid:34) iβ ρ f ω ∂ ν x E (cid:62) ( x, y ) ν y + β ρ f ω ∂ ν x ( ∂ ν y E ( x, y )) (cid:35) p ( y ) ds y . We first consider the following term T ( ∂ x , ν x ) (cid:90) Γ ( T ( ∂ y , ν y ) E ( x, y )) (cid:62) u ( y ) ds y . Set f ( x ) = (cid:90) Γ ∇ y [ γ k s ( x, y ) − γ k ( x, y )] ν (cid:62) y u ( y ) ds y ,f ( x ) = (cid:90) Γ ∇ y [ γ k ( x, y ) − γ k ( x, y )] ν (cid:62) y u ( y ) ds y ,f ( x ) = (cid:90) Γ ∂ ν y γ k s ( x, y ) u ( y ) ds y ,f ( x ) = (cid:90) Γ [2 µE ( x, y ) − γ k s ( x, y ) I ] M ( ∂ y , ν y ) u ( y ) ds y , and g i ( x ) = µ∂ ν x f i ( x ) + ( λ + µ ) ν x ∇ x · f i ( x ) + µM ( ∂ x , ν x ) f i ( x ) . Then we obtain from (A.1) that g ( x ) = ( λ + 2 µ ) (cid:90) Γ (cid:2) k s γ k s ( x, y ) − k γ k ( x, y ) (cid:3) ν x ν (cid:62) y u ( y ) ds y +2 µ dds x (cid:90) Γ A ∇ y [ γ k s ( x, y ) − γ k ( x, y )] ν (cid:62) y u ( y ) ds y , (A.6)and g ( x ) = ( λ + 2 µ ) (cid:90) Γ (cid:2) k γ k ( x, y ) − k γ k ( x, y ) (cid:3) ν x ν (cid:62) y u ( y ) ds y +2 µ dds x (cid:90) Γ A ∇ y [ γ k ( x, y ) − γ k ( x, y )] ν (cid:62) y u ( y ) ds y . (A.7)19t follows from [3] that g ( x ) can be expressed as g ( x ) = µ dds x (cid:90) Γ γ k s ( x, y ) duds y ds y + µk s (cid:90) Γ γ k s ( x, y ) ν (cid:62) x ν y u ( y ) ds y +( λ + µ ) (cid:90) Γ ν x ∇ (cid:62) x ∂ ν y γ k s ( x, y ) u ( y ) ds y + µ dds x (cid:90) Γ ∂ ν y γ k s ( x, y ) u ( y ) ds y . (A.8)For g ( x ) , we know from (A.5) that g ( x ) = µ (cid:90) Γ ν (cid:62) x ∇ x γ k s ( x, y ) A duds y ds y − µ (cid:90) Γ ν x ∇ (cid:62) x [ γ k s ( x, y ) − γ k ( x, y )] A duds y ds y + 2 µ ( k − q ) k − k (cid:90) Γ ν x ∇ Tx [ γ k ( x, y ) − γ k ( x, y )] A duds y ds y +4 µ dds x (cid:90) Γ AE ( x, y ) A duds y ds y + 3 µ dds x (cid:90) Γ γ k s ( x, y ) duds y ds y − ( λ + µ ) (cid:90) Γ ν x ∇ (cid:62) x γ k s ( x, y ) A duds y ds y . (A.9)Therefore, (A.6)-(A.9) yields T ( ∂ x , ν x ) (cid:90) Γ ( T ( ∂ y , ν y ) E ( x, y )) (cid:62) u ( y ) ds y (A.10) = − g ( x ) + k − qk − k g ( x ) + g ( x ) + g ( x )= − ( ρ − βρ f ) ω (cid:90) Γ γ k s ( x, y )( ν x ν (cid:62) y − ν (cid:62) x ν y I − J ν x ,ν y ) u ( y ) ds y + (cid:90) Γ (cid:20) k ( λ + 2 µ )( k − q ) k − k γ k ( x, y ) − k ( λ + 2 µ )( k − q ) k − k γ k ( x, y ) (cid:21) ν x ν (cid:62) y u ( y ) ds y +4 µ dds x (cid:90) Γ γ k s ( x, y ) duds y ds y + 4 µ dds x (cid:90) Γ AE ( x, y ) A duds y ds y − µ (cid:90) Γ ν x ∇ (cid:62) x [ γ k s ( x, y ) − γ k ( x, y )] A duds y ds y − µ dds x (cid:90) Γ A ∇ y [ γ k s ( x, y ) − γ k ( x, y )] ν (cid:62) y u ( y ) ds y + 2 µ ( k − q ) k − k (cid:90) Γ ν x ∇ (cid:62) x [ γ k ( x, y ) − γ k ( x, y )] A duds y ds y + 2 µ ( k − q ) k − k dds x (cid:90) Γ A ∇ y [ γ k ( x, y ) − γ k ( x, y )] ν (cid:62) y u ( y ) ds y . On the other hand, we obtain from (A.3) that (cid:90) Γ T ( ∂ x , ν x ) E ( x, y ) ν (cid:62) y u ( y ) ds y (A.11) = − γk − k (cid:90) Γ (cid:2) k γ k ( x, y ) − k γ k ( x, y ) (cid:3) ν x ν (cid:62) y u ( y ) ds y − µγ ( λ + 2 µ )( k − k ) dds x (cid:90) Γ A ∇ y [ γ k ( x, y ) − γ k ( x, y )] ν (cid:62) y u ( y ) ds y , (cid:90) Γ ν x ( T ( ∂ y , ν y ) E ( x, y )) (cid:62) u ( y ) ds y (A.12) = − γk − k (cid:90) Γ (cid:2) k γ k ( x, y ) − k γ k ( x, y ) (cid:3) ν x ν (cid:62) y u ( y ) ds y − µγ ( λ + 2 µ )( k − k ) dds x (cid:90) Γ ν x ∇ (cid:62) x ( γ k ( x, y ) − γ k ( x, y )) A duds y ds y . Then Lemma 3.3 can be proved by combining (A.11)-(A.13). The proofs of Lemma 3.4-3.6 are analogous,and thus are omitted here.