{'fraction_non_alphanumeric': 0.04645622394282311, 'fraction_numerical': 0.017718880285884454, 'mean_word_length': 4.539175257731959, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the hyper-singular boundary integral equation methods for dynamic poroelasticity: three dimensional case

In our previous work [SIAM J. Sci. Comput. 43(3) (2021) B784-B810], an accurate hyper-singular boundary integral equation method for dynamic poroelasticity in two dimensions has been developed. This work is devoted to studying the more complex and difficult three-dimensional problems with Neumann boundary condition and both the direct and indirect methods are adopted to construct combined boundary integral equations. The strongly-singular and hyper-singular integral operators are reformulated into compositions of weakly-singular integral operators and tangential-derivative operators, which allow us to prove the jump relations associated with the poroelastic layer potentials and boundary integral operators in a simple manner. Relying on both the investigated spectral properties of the strongly-singular operators, which indicate that the corresponding eigenvalues accumulate at three points whose values are only dependent on two Lam\'e constants, and the spectral properties of the Calder\'on relations of the poroelasticity, we propose low-GMRES-iteration regularized integral equations. Numerical examples are presented to demonstrate the accuracy and efficiency of the proposed methodology by means of a Chebyshev-based rectangular-polar solver.


Introduction
The dynamic poroelastic problems describing the physical behavior of the wave propagation in the elastic solid and the interstitial fluid can be found in many fields of applications such as petroleum industry, materials science, soil mechanics and biomechanics, etc. In accordance to Biot's theory [6][7][8][9][10]19,22], the dynamic poroelastic problems can be modeled by the coupled equations of the pore pressure and the solid displacement field, and the targeted degrees of freedom can be changed [39]. For the numerical solutions of such kind of wave scattering problems, it is known that the boundary integral equation (BIE) methods [21,27] take advantages over the volumetric discretization methods [23,25,31,42] in the sense of dimensions reduction, discretization of boundary and natural satisfactory of radiation condition, while in particular, the volumetric methods requires introducing appropriate artificial boundary conditions, such as absorbing boundary conditions or perfectly matched layers for the treatment of problems on unbounded domains. As a continuation of our previous work [44] for the two-dimensional poroelastic scattering problems, this work is devoted to proposing efficient BIE methods for solving the three-dimensional problems [18, 33-35, 37, 38, 41] with Neumann boundary condition and it requires more complex technical investigations of the poroelastic boundary integral operators (BIOs).
In the classical BIE theory, both the direct methods based on Green's formula and the indirect methods based on potential theory have been extensively discussed. In practice, the combined boundary integral equations (CBIEs) resulting from a combination of single-layer and double-layer BIOs (for Dirichlet case) or a combination of double-layer and hyper-singular BIOs (for Neumann case) are generally employed to avoid the influence of possible eigenfrequencies. In this work, we employ both the direct method and the indirect method to construct two types of CBIEs for solving the three-dimensional poroelastic problems with Neumann boundary condition. As mentioned in [44], it still remains open to prove the unique solvability of the CBIEs for the poroelastic scattering problem, but these BIEs still can provide efficient numerical tools for the solutions of the problems imposed on unbounded domains. Then analogous to the two-dimensional case [44], the following three issues should be addressed: (i). The jump relations between the layer-potentials and the BIOs in dynamic poroelasticity case are not easy to be observed.
(ii). The double-layer operators K, K (strongly-singular) and hyper-singular operator N are well defined in the sense of Cauchy principle value and Hadamard finite part [27], respectively. Then it requires appropriate solvers for the accurate evaluation of these operators.
(iii). It is known that the eigenvalues of the hyper-singular operator accumulate at infinity and as a result, solving the CBIEs by means of Krylov-subspace iterative solvers, such as GMRES, generally requires a relatively large number of iterations for the convergence of numerical solution. Then low-GMRES-iteration integral formulations are highly desirable.
To resolve the first and second issues, it is necessary to take a comprehensive study on the single kernels of the poroelastic BIOs. To reduce/transform the singularities, some methodologies, for instance, addingand-subtracting appropriate terms and regularization using integration-by-parts, have been discussed in open literatures. Inspired by the idea of reformulating the acoustic/Laplace hyper-singular integral operator into a combination of weakly-singular integral operators and tangential derivatives [27,36], a novel regularization technique using Günter derivative and Stokes formulas has been developed for the elastic and thermoelastic problems [3,4,32,43]. In two-dimensions, the Günter derivative can be simplified as the classical tangential derivative multiplied by a constant matrix, and the regularized formulations for two-dimensional poroelastic BIOs have been investigated in [44]. But the three-dimensional Günter derivative is more complex. Although the thermoelastic problem [30] takes a similar Biot's model as the poroelastic problem, the results presented in [4] can not be extended to the three-dimensional poroelastic case trivially regarding to the more complicated coupled boundary operator and its adjoint, see Section 3. It is proved in Theorems 4.1-4.6 that the three-dimensional strongly-singular and hyper-singular poroelastic integral operators can be re-expressed in terms of multiple weakly-singular integral operators and tangential-derivative operators. Compared with the formulations given in [34,35], the derived regularized expressions in this work take simpler forms, and as a consequence, the jump relations between the layerpotentials and the BIOs in dynamic poroelasticity case can be proved in an extremely simple manner, see Theorem 3.1 and its proof in Section 4.3. In addition, owing to the new regularized expressions, the numerical evaluation of the poroelastic BIOs amounts to the evaluation of weakly-singular type integrals, for which the so-called Chebyshev-based rectangular-polar method proposed in [12] is applicable, and the evaluation of three-dimensional tangential derivatives can be implemented via FFT [15].
The third issue is related to the spectral regularization or preconditioning. In addition to algebraic preconditioning approaches, such as sparse approximate inverse [5,17] and multigrid methods [28], the analytical preconditioning approach based on the Calderón relation, which in fact utilizes the compositions N S, SN of single-layer operator S and hyper-singular operator N , has been discussed for solving wave scattering problems by closed surfaces [11,15,20] or open surfaces [2,13,15]. Due to the weakly-singularity of double-layer operator, it follows easily that the acoustic Calderón relation can be viewed as a compact perturbation of an identity operator for the smooth closed-surface case. But this does not hold trivially in elastic case, and also in poroelastic case (which can be understood naturally since the elastic single kernels are involved in the poroelastic BIOs), on account of the fact that the classical elastic doublelayer operators are not compact. It has been proved in [14] for two-dimensional case and in [15] for three-dimensional case that the elastic double-layer operators K, K are polynomially compact and as a result, the values of the finite accumulation points of the eigenvalues of N S, SN only depend on the Lamé parameters of the elastic medium. The two-dimensional poroelastic case has been discussed in [44], whereas the result does not directly fit for the three-dimensional context. It is shown in this work (see Theorem 3.3) that the three-dimensional poroelastic double-layer operators are compact in the sense of a third-order polynomial and interestingly, the corresponding accumulation points of the eigenvalues are independent of the poroelastic parameters (see Table 1) except the two Lamé parameters. On a basis of the spectral properties of the poroelastic BIOs and analogous to the two-dimensional approach [44], we propose two regularized CBIEs for which the eigenvalues of the combined integral operators are bounded away from zero and infinity and then it leads to significant reductions in the number of GMRES iterations required for convergence to a given residual tolerance over the original CBIEs.
This paper is organized as follows. The dynamic poroelastic scattering problem is introduced in Section 2 and then in Sections 3.1-3.2, we present both the direct and indirect methods to derive the classical CBIEs, respectively. Section 3.3 is arranged to give a theoretical investigation of the spectral properties of the poroelastic integral operators and the corresponding Calderón relation. In Section 4, regularized expressions of the strongly-singular and hyper-singular operators are presented and then the jump relations between the layer-potentials and the BIOs are proved. Section 5.1 proposes two new RBIEs based on the Calderón relation and Section 5.2 briefly describes the numerical discretization method for poroelastic BIOs. Some numerical examples are presented in Section 6 to demonstrate the accuracy and efficiency of the proposed method.

Poroelastic problem
Let Ω be a bounded domain in R 3 with smooth boundary Γ := ∂Ω, and its exterior complement is denoted by Ω c = R 3 \Ω. This work is devoted to studying the numerical solutions of the three-dimensional time-harmonic problems of wave propagation in the domain Ω c which is occupied by a linear isotropic poroelastic medium characterized by the physical parameters listed in Table 1.
Following the Biot's theory [6,8,9], the solid displacements u = (u 1 , u 2 , u 3 ) ∈ H 1 loc (Ω c ) 3 and the pore pressure p ∈ H 1 loc (Ω c ) characterizing the wave propagation in poroelastic medium can be modeled by the following coupled partial differential equations 1) or equivalently, in a matrix form are abbreviations defined to simplify the representation of the problem, ω denotes the frequency, I is the identity operator and ∆ * denotes the Lamé operator given by with σ(u) = λ(∇ · u)I + 2µ ε(u) and ε(u) = 1 2 (∇u + (∇u) ).
Here, we denote by k p and k s the compressional and shear wave numbers, respectively, and they are given by The wave numbers k 1 , k 2 in (2.3), which represent the wave numbers of the fast compressional wave and slow compressional wave in poroelastic medium, respectively, are determined through . In particular, To complete the statement of the poroelastic problem, we assume that the solution U satisfies the following Kupradze radiation conditions as r = |x| → ∞ for l = 1, 2, 3 and j = 1, 2, Remark 2.1. For the poroelastic problem, the degrees of freedom can be determined in different ways [39]. Compared with the formulation in terms of the solid displacement and the fluid displacement, and the formulation in terms of the solid displacement and the seepage displacement, the above model enjoys the lowest number of unknowns. For the uniqueness analysis of the dynamic poroelastic problem, we refer to [24].

Boundary integral equations
In this section, we introduce the hyper-singular BIEs for solving the poroelastic problem together with some theoretical study of the properties of BIOs. Based on the Green's identities and potential theory, direct and indirect boundary integral formulations are derived, respectively. We begin with the first and second Green's identities for the poroelastic problems in Ω (analogous to the problem in Ω c ). For U = (u , p) and V = (v , θ) , the first Green's identity reads while the second Green's identity admits Here, A Ω denotes a bilinear form defined by where L * denotes the adjoint operator of L defined by and T * (∂, ν) is the boundary operator given by It is known [18,38] that the fundamental solution of the operator L * in R 3 is given by is the fundamental solution of the Helmholtz equation in R 3 with wave number k t .

Direct method
It follows from the Green's formulas (3.2) that the solution of (2.1) can be represented in the form where S and D are the single-layer and double-layer potentials given by respectively. Introduce the BIOs for the poroelasticity in the sense of principle value or Hadamard finite part as follows where S, K, K and N are called, respectively, the single-layer, double-layer, transpose of double-layer, and hyper-singular BIOs. Then we conclude the jump relation results associated with the poroelastic layer potentials and BIOs in the following theorem.
Theorem 3.1. For x ∈ Γ, the following jump relations hold: Remark 3.2. The proof of this theorem relies on the study of the regularized expressions of the integral operators that will be presented in Section 4 and thus, will be reported after that.
Now applying the jump conditions, we are led to the BIEs on Γ and Combining the BIEs (3.11)-(3.12) results into the so-called Burton-Miller formulation [16] on Γ where η = 0 is a combination coefficient. Using the boundary condition (2.2), we obtain the direct combined boundary integral equation (DCBIE) (3.14)

Indirect method
The indirect boundary integral formulations can also be used for solving the poroelastic problems, which also allow for a suitable tool to test each operator separately. From the potential theory, the unknown function U of (2.1) can be represented by a combination of the single-layer and double-layer potentials Operating with the boundary operator T (∂, ν) on (3.15), taking the limit as in Theorem 3.1 and applying the boundary condition (2.2), we can obtain the indirect combined boundary integral equation (ICBIE)

Operator properties
Assuming that the boundary Γ is sufficiently smooth, the BIOs are continuous mappings between the following spaces [27,43] S : and the following Calderón relations hold: As analytical preconditioning techniques, the Calderón relations have been investigated and utilized in regularized BIE methods [11,15,44], which require the spectral study of the BIOs, to construct BIE systems possessing highly favorable spectral properties. The main reason is that the eigenvalues of the hyper-singular integral operator N accumulate at infinity. As a result, obtaining the solutions of some integral equations, for example (3.14) and (3.16) in this work, by means of Krylov-subspace iterative solvers such as GMRES generally requires large numbers of iterations. To overcome this difficulty, the spectral properties of the integral operators K, K and the associated Calderón relations N S for two-dimensional poroelasticity are investigated in [44] and then a regularized BIE method is proposed. However, as proved in the following theorem, the three-dimensional integral operators K and K enjoy spectral properties different from those in the two-dimension case. Theorem 3.3. Let I λ,µ be a matrixed operator given by Furthermore, the spectrum of K consists of three non-empty sequences of eigenvalues which converge to 0, C λ,µ and −C λ,µ respectively.
Proof. Analogous to the proof of [44, Theorem 3.1], it is sufficient to consider the static (ω = 0) BIO corresponding to K which can be formulated as being the fundamental solution of Lamé equation. It can be verified that K 3,0 = 0 and the kernels of K j,0 , j = 2, 4 admit weak singularity implying that K j,0 , j = 2, 4 are compact. From [1], we know that The compactness of K −K 0 indicates the compactness of K c and then further implies that K (K 2 − I λ,µ ) is compact. The proof is completed.
The spectral property of the operator K is similar to K . Relying on the Calderón relations [27] we can conclude that the spectrum of bith the composite operators N S and SN consist of two non-empty sequences of eigenvalues which converge to − 1 4 and − 1 4 + C 2 λ,µ , respectively. To verify the above results numerically, we consider the problem of poroelastic scattering by a unit ball, and we choose the same values of parameters as in Section 6. Consequently, the constant C λ,µ = 0.1875. Utilizing the following re-expressions of the BIOs together with the discretization method presented in Sections 4-5, the eigenvalue distributions of the integral operators K , K, N S and SN are displayed in Figure 1 showing an agreement with our theoretical results. Based on these results, the corresponding regularized BIE method is proposed in Section 5 for solving the three-dimensional poroelastic problem.

Strong-singularity and hyper-singularity regularization
As aforementioned, the integral operators K (as well as K ) and N are strongly-singular and hypersingular, respectively. In this section, we will present new (regularized) expressions for these operators. More precisely, with the help of the tangential Günter derivative [4], these operators will be expressed in terms of compositions of weakly-singular integral operators (K i j , K i j , N i j appearing in the following Theorems (4.1)-(4.6)) and tangential-derivative operators (The formulations to derive these regularized expressions will be shown in the supplemented material.). Employing these formulations together with the rectangular-polar quadrature method proposed in [12,15] and the iteration solver GMRES then leads to our boundary integral solver for the poroelastic problem.
We begin with the Günter derivative operator M (∂, ν) defined as From [15], we know that , in which ∇ S is the surface gradient defined as ∇ S = ∇u − ν∂ ν u. Then the traction operator T (∂, ν) can be rewritten as

Strong-singularity regularization
We first consider the operators K and K in the forms of Then the following regularized formulations can be obtained.

Hyper-singularity regularization
In this subsection, we investigate the hyper-singular operator N with The regularized formulations for the operators N j , j = 1, · · · , 4 are given in the following theorems.
Theorem 4.3. The hyper-singular operator N 1 can be expressed as y))] ν y u(y)ds y , with the constants C i , i = 1, 2, 3, 4 being given by and where, for a scalar field p and a vector field v, the operators τ 1 and τ 2 are defined by Theorem 4.4. The hyper-singular operator N 2 can be expressed as Theorem 4.5. The hyper-singular operator N 3 can be expressed as Theorem 4.6. The hyper-singular operator N 4 can be expressed as

Proof of Theorem 3.1
Thanks to the derived regularized formulations of the integral operators, now we can prove the jump conditions stated in Theorem 3.1.
The jump condition for the single-layer potential follows trivially. Relating to the proof of Theorem 4.1, it can be concluded that D(ϕ)(z), x ∈ Γ takes a form similar to (4.1) for z = x ± hν x / ∈ Γ, i.e., Letting h → 0 + (i.e. z → x ∈ Γ), and applying the classical jump relations for acoustic and elastic problems [27], it is only necessary to study the jumps of D 1 1 , D 3 2 , D 2 3 and D 1 4 in the forms associated with K 1 1 , K 3 2 , K 2 3 and K 1 4 , respectively. Note that for x ∈ Γ, It follows that lim h→0 + ,z=x±hνx Therefore, we have Analogously, it can be derived that T (∂ z , ν x )S(ϕ)(z), x ∈ Γ takes a form similar to (4.2) for z = x ± hν x / ∈ Γ. Due to the fact that for x ∈ Γ, lim h→0 + ,z=x±hνx Γ ν x · ∇ z γ ks (z, y)u(y)ds y = ∓ 1 2 u(x) + Γ ∂ νx γ ks (x, y)u(y)ds y , we arrive at the jump conditions It remains to prove the jump conditions for T (∂ z , ν x )D(ϕ)(z) as h → 0 + . We can write it as and it can be proved that N j , j = 1, · · · , 4 take forms similar to N j , j = 1, · · · , 4, respectively. We first study N 1 in the form and thus, In addition, note that y))] ν y u(y)ds y +µ Γ ∂ νy γ ks (z, y)u(y)ds y , and lim h→0 + ,z=x±hνx Therefore, which gives that The proof of the other three jump conditions is analogous and hence is omitted. This completes the proof of Theorem 3.1.

Numerical discretization
In this section, we briefly introduce the application of the Chebyshev-based rectangular-polar solver discussed in [12,15] to the numerical discretization the poroelastic BIOs. Based on a partition of the boundary using non-overlapping parametric curvilinear patches, this approach interpolates the unknowns on a Chebyshev grid on each patch in terms of Chebyshev polynomials. For the corresponding acceleration of this method, we refer to [29].
Let Γ be partitioned into a set of M non-overlapping parametrized (logically-rectangular) patches Γ q , q = 1, ..., M as Introducing the tangential covariant basis vectors and surface normal on Γ q : we can obtain the metric tensor as where g q ij = a q i · a q j and thus, the surface element Jacobian is given by Here, |G q | is the determinant of G q . As a result, the surface gradient of a given density ϕ = ϕ(r q (u, v)) can be expressed as where g ij , i, j = 1, 2 denote the components of the inverse of the matrix G q . Given a density ϕ, it can be approximated on Γ q by the Chebyshev polynomials as and the coefficients ϕ q ij = ϕ(x q ij ) denote the values of the continuous density ϕ at the discretization points Therefore, we can obtain the approximation of the surface gradient as We now discuss the discretizations of the prorelastic BIOs. On a basis of the regularized formulations of integral operators given in Sections 4, the numerical implementations can be converted into evaluating multiple operators of two types, (i)Integral operators with weakly-singular kernels H(x, y) Hϕ(x) = Γ H(x, y)ϕ(y)ds y , and (ii) surface-differentiation operators M (∂, ν), τ 1 , τ 2 appearing in Section 4, which can be extracted from the approximation of surface gradient ∇ S x . Clearly, the integrals Hϕ(x) over Γ can be split into the sum of integrals over each of the M patches, In the "non-adjacent" integration case, in which the target point xq ij is far from the integration patch, the integral H q (xq ij ) is non-singular. Then the classical Fejér's first quadrature rule can be utilized to obtain the approximation where Aq ,q ij,nm = H(xq ij , r q (u n , v m ))J q (u n , v m )w n w m , with the quadrature weights In the "adjacent" integration case, in which the point xq ij either lies within the integration patch or is located very close to it, the integral H q (xq ij ) becomes weakly-singular and nearly-singular, respectively. It is suggested in [12] to constructed a new graded mesh, relying on a smoothing change of variables ξ u (s) (for more details, see [12,15]), around the point which lies closest to xq ij with parameters Then the single integral can be approximated by H(xq ij , r q (ξũq (s), ξṽq (t)))J q (ξũq (s), ξṽq (t)) ×a nm (ξũq (s), ξṽq (t))ξ ũ q (t l1 )ξ ṽ q (t l2 )w l1wl2 .
Here, the quadrature nodes and weights for an order N β are given analogous to the "non-adjacent" case.

Numerical experiments
In this section, several numerical examples, involving three bounded obstacles depicted in Fig. 3, are presented to demonstrate the accuracy and efficiency of the proposed methods for solving threedimensional poroelastic problems. Utilizing the dimensionless technique discussed in [18,40], we set µ = 2, ν p = 0.2, ν u = 0.33, B = 0.62, C = 0.66, φ = 0.333, ρ s = 1, ρ f = 0.5, κ = 1. We additionally choose η = 1. Here, we use the fully complex version of the iterative solver GMRES to produce the solutions of the integral equations and the maximum errors defined by We first test the accuracy of the proposed methods. The exact solution U exa = (u exa , p exa ) is given by with z = (0, 0.5, 0.3) for the obstacle Fig. 3(a) and z = (0, 0, 0) for the obstacles Fig. 3(b,c), which gives the boundary data F = T (∂, ν)U exa on Γ. We first consider ω = π and the Chebyshev grid with M = 6, N = 16 and N β = 100. Fig. 4 shows the errors |U num (x) − U exa (x)| between the numerical solution U num resulting from solving DRBIE (5.2) and the exact solution U exa for the obstacle Fig. 3(a) and x ∈ {x ∈ R 3 : |x| = 2} with a maximum value 1.814 × 10 −8 . Next, we consider the poroelastic problem of scattering by obstacle Fig. 3 Fig. 5. The relative error for this case is ∞ = 6.67 × 10 −6 . In Fig. 6, the numerical errors ∞ with respect to N using DRBIE and IRBIE for the obstacles Fig. 3 are presented while choosing ω = 2π and Chebyshev grid with M = 6. Higher accuracy can be achieved by increasing the parameter N β and treating the evaluation of weakly-singular kernels for small |x − y| with cares. Next, we verify the efficiency of the regularized integral equation methods. Choosing ω = 2π and the Chebyshev grid with M = 6, N = 32 and N β = 200, Fig. 7 displays the history of GMRES residuals as functions of the number of iterations for the method of using DCBIE (3.14), ICBIE (3.16), DRBIE (5.2) and IRBIE (5.4), respectively. The rapid convergence results of regularized methods demonstrate that use of the regularized integral equations is highly beneficial compared to the un-regularized ones.
With ω = 20 and M = 24, Table 6 lists the precomputation time, time per iteration and number of iterations required by the regularized integral equation methods. For the DRBIE method, an accuracy of 7.6 × 10 −3 (resp.2.5 × 10 −6 ) can be achieved by setting N = 8 (resp.N = 16), while an accuracy of 2.1 × 10 −2 (resp.1.7 × 10 −6 ) can be obtained for the IRBIE method. Figure 4: The numerical errors on a sphere for the poroelastic problem of sacttering by the obstacle Fig. 3(a). Figure 5: Comparison of the exact and numerical solutions for the poroelastic problem of scattering by obstacle Fig. 3(b).
Finally, we consider the scattering of an incident point source U inc in the form by the obstacle Fig. 3 Appendix. Regularized expressions of the strongly-singular and hyper-singular operators and proofs.
This appendix presents the main approach for deriving the regularized formulations of the stronglysingular and hyper-singular operators. Analogous technique has been shown in [4] for the threedimensional elastic and thermoelastic problems. But the derivations for three-dimensional poroelastic problem are more complex and thus, we present the full proof to make this appendix individually readable.
Proof. The operator K can be written as Using (A.7) and (A.9), we have : Real parts of the total field U around the obstacle for the scattering of an incident point source by the obstacle Fig. 3(b).
The formula for K 4 (p)(x) can be obtained directly from its definition and this completes the proof of Theorem 4.1.
Following the result (A.17), we have