A compact Green's function for general dispersive anisotropic poroelastic media in a full-frequency regime is presented for the first time. First, starting in a frequency domain, the anisotropic dispersion is exactly incorporated into the constitutive relationship, thus avoiding fractional derivatives in a time domain. Then, based on the Radon transform, the original three-dimensional differential equation is effectively reduced to a one-dimensional system in space. Furthermore, inspired by the strategy adopted in the characteristic analysis of hyperbolic equations, the eigenvector diagonalization method is applied to decouple the one-dimensional vector problem into several independent scalar equations. Consequently, the fundamental solutions are easily obtained. A further derivation shows that Green's function can be decomposed into circumferential and spherical integrals, corresponding to static and transient responses, respectively. The procedures shown in this study are also compatible with other pertinent multiphysics coupling problems, such as piezoelectric, magneto-electro-elastic and thermo-elastic materials. Finally, the verifications and validations with existing analytical solutions and numerical solvers corroborate the correctness of the proposed Green's function.

2019 The Author(s) Published by the Royal Society. All rights reserved.

The fundamental solution, i.e. Green's function, is the solution of a physical field (e.g. displacement, electrical fields, etc.) recorded from a receiver due to a unit excitation at a source point. Green's function plays an indispensable role in the realm of integral equations, whose numerical implementation in matrix form is usually called the method of moments in computational electromagnetics [1, p. 71]. The integral equations include (i) surface integral equations, based on the surface equivalence principle [2,3], and (ii) volume integral equations, based on the volume equivalence theorem [1, p. 71]. In computational mechanics, the numerical implementation of the former is widely known as the boundary element method [4,5; 6, p. 3] which is widely used in fracture detection [7] and crack propagation simulations [8][9][10][11]. On the other hand, an important application of volume integral equations is the homogenization technique. For composite problems, the micro-structures are usually on a scale much smaller than a wavelength, requiring computationally intractable dense meshes for conventional numerical solvers; instead, a homogenization technique is preferred to provide the effective material properties [12][13][14]. The homogenization is contingent upon the Eshelby tensor/depolarization dyadic, which is an integral of Green's function [13,15,16; 17, p. 3-3; 18, ch. 4].

The applicability of integral equations depends on (i) the existence of a Green's function, which is usually restricted to a full-space problem or stratified media, and (ii) the calculation complexity, which is usually accelerated by the conjugate-gradient fast Fourier transform method [19] or the fast multipole method [20]. Therefore, a compact and efficient expression of Green's function is always desirable.

However, to the best of the authors' knowledge, a Green's function for anisotropic dispersive poroelastic media is not available yet. There are three difficulties in relation to this: anisotropy, dispersion and poroelasticity. First, anisotropy naturally arises from the aligned microstructures of the crystals [16]. The fundamental solution becomes non-trivial owing to the fact that anisotropy undermines the simplicity of those derivations based on the isotropy assumption [21][22][23][24]. For instance, an explicit expression of Green's function for isotropic media is just a scalar function [1, pp. 10, 198], whereas extensive efforts have to be devoted to deriving a Green's function for anisotropic materials. Typically, the Fourier transform is used pervasively [25, ch. 2.10; 26,27]. However, the underlying problem of the Fourier transform-based methods requires an infinite integral. The elimination of this integral, leading to an explicit expression, can be achieved by means of the Cauchy residue theorem for static problems, such as anisotropic elastic solids [28][29][30] and piezoelectric solids [31]. Fortunately, the Radon transform-based Green's function was developed, where only a finite integral is involved [32][33][34]. Second, dispersive materials exhibit frequency-dependent properties. In time domain, one has to confront either convolutional integrals or fractional derivatives [35][36][37][38]. In the frequency domain, complexvalued material parameters are introduced, thus leading to distinct properties, with respect to the original real-valued system [39]. From an experimental study, it has been found that serpentinite rock requires 21 independent parameters, in the quality factor matrix, to fully depict the anisotropic attenuations in elastic media [40]. From a mathematical perspective, the elasticity matrix becomes non-Hermitian, but is still symmetric, thus the eigenvectors are no longer orthogonal. In addition, the roots of the characteristic polynomial for the Kelvin-Christoffel matrix are no longer in complex conjugate pairs. Therefore, the existing fundamental solutions, with the underlying assumption of a non-dispersive material, may require a prudent modification [33,[41][42][43]. Third, poroelastic materials require many input parameters [44][45][46]. In addition, the pore fluid behaves differently in low-and high-frequency regimes. Especially in the high-frequency regime, the pore fluid is a dispersive medium involving fractional derivatives [47,48].

Therefore, the main goal of this study is to first succinctly propose a unified formulation of Green's function for anisotropically dispersive fluid-saturated porous media. Even though the fundamental solutions for anisotropic/viscoelastic/poroelastic media have been around for a while [43,46,49], they are more from a collection of derived formulations, thus lacking a comprehensive procedure for complete anisotropically dispersive poroelastic materials. It is therefore desirable to give a compact formulation of the corresponding Green's function in this paper, equipped with the following three steps. First, the dispersion is exactly incorporated into the frequency-domain constitutive relationship, thus the temporal fractional derivatives are avoided. As a result, we propose a full anisotropic-Q matrix to quantitatively depict the anisotropic dispersion in porous media. Second, we take full advantage of the Radon transform, by transforming the original three-dimensional space differential equation into a one-dimensional system in local coordinates [50]. Third, inspired by previous work on the characteristic analysis of hyperbolic equations [21,35,45,[51][52][53], the eigenvector diagonalization strategy is adopted to decouple the components of a vector into separate scalar equations. Consequently, a closedform of Green's function is succinctly expressed as a circumferential integral and a spherical integral, corresponding to static and transient responses, respectively. Finally, the results are verified and validated with analytical and numerical solutions. It is worth mentioning that the above procedure is also compatible with other pertinent multi-physics coupling problems, such as piezoelectric, magneto-electro-elastic and thermo-elastic materials [41,54]. This will also smooth the way for new numerical solver verifications.

To summarize, this research delivers the following novelties:

(i) Green's function for general anisotropically dispersive poroelastic media is presented for the first time. (ii) The Radon transform is applied to effectively reduce the three-dimensional spatial derivatives into one dimension. (iii) The eigenvalue diagonalization is further applied to simplify the non-Hermitian anisotropic problem. (iv) The results are verified and validated with existing analytical solutions and numerical solvers.

This paper is organized as follows. Section 2 lists the notations used in the paper. The Fourier transform and Radon transform are introduced in §3, which are preparatory for the subsequent sections. In §4, the governing equations for poroelastic waves are provided in a fullfrequency range. Then, § §5-7 elaborate the full anisotropic dispersion, the Radon transform-based fundamental problem and the eigenvector diagonalization technique, respectively. Extensions to relevant problems are illustrated in §8. Section 9 elucidates the numerical implementation, verification and validation. Concluding remarks are given in §10.

For convenience, we list the various notations used in this paper in alphabetical order below.

-The summation convention is not applied unless specified.

-A is a matrix, either symmetric or non-symmetric, either Hermitian or non-Hermitian, either complex or real. -B i is magnetic flux density.

-C is the elasticity matrix with size 6 × 6.

-D is the elasticity matrix for poroelastic media.

-D i is the electric flux density.

-D u IJ is the stiffness tensor of the undrained (i.e. unjacketed) material, subject to the Voigt notation [18, p. 24; 55], where I, J = 1, 2, 3, 4, 5, 6.

-F is the force matrix in a diagonalized system. -F is the Fourier transform operator.

-G is the dyadic Green's function.

-Ĝ is the Radon transform dyadic Green's function.

-G R is the transient part of the dyadic Green's function. -G S is the static part of the dyadic Green's function.

-I is the identity matrix.

-K is a diagonal matrix, the generalized wavenumber.

-L is the first spatial derivative operator matrix.

-M is the fluid-solid coupling modulus.

-M is a modulus, any entry of D.

-M r are the reference moduli measured at ω r .

-P i 1/2 is the Pride number.

-Q is the quality factor -R is the Radon transform operator.

-R is the right eigenvector matrix, assembled in a column manner.

-T i is the tortuosity vector of the solid matrix.

-U is the left eigenvector matrix, assembled in a row manner.

-W is the dependent variable matrix in a diagonalized system. -α I is the generalized Biot's effective-stress coefficient. Specifically, it is the ratio of the pore-fluid pressure, which causes the same amount of strain as the the total stress [55, p. 225]. -b is the component of n along e.

δ(•) is the Dirac delta function.

δ IP is the Kronecker delta.

-(e, t 1 , t 2 ) are the local coordinate bases.

ε is the electrical permittivity matrix with size 3 × 3.

ij is the skeleton strain

where i, j = 1, 2, 3. -f c is the central frequency for the source time function.

-f i is the body force density imposed in the elastic system. -f (1) i is the body force density imposed in the matrix frame. -f (2) i is the body force density imposed in the pore fluid.

-f e is the static electrical charge density source.

-f m is the static magnetic charge density source.

-Γ (•) is the Kelvin-Christoffel matrix.

-j is the imaginary number.

κ i is the permeability of the solid matrix.

λ I is the diagonal entry of Λ.

-Λ is the right/left eigenvalue matrix.

-Λ i /2 is interpreted, approximately, as the pore-volume to pore-surface ratio.

μ is the magnetic permeability matrix with size 3 × 3.

n is the Radon transform projection direction, a unit vector.

ν is the fluid viscosity.

ω is the angular frequency.

ω c i is the critical frequency,

ω r is the reference angular frequency.

-p is the fluid pressure for the pore.

p is the third-rank piezoelectric tensor, written in matrix form with size 3 × 6.

φ is the matrix porosity.

-Ψ i is the time-domain viscodynamic operator in the x i direction.

-Ψi is the frequency-domain viscodynamic operator in the x i direction.

q is the third-rank piezomagnetic tensor, written in matrix form with size 3 × 6. ρ is the density of the poroelastic medium: ρ = (1 -φ)ρ s + φρ f , where φ is the porosity of the poroelastic medium. -ρ is the density matrix.

ρ w i is fluid inertia along axis i,

ρ f is the density of the fluid.

ρ s is the density of the solid frame.

-s is the Radon transform space.

-t is the time axis.

t is the orthogonal unit vector with respect to e.

τ ij is the total stress of the poroelastic medium.

θ is the azimuthal angle in the local coordinate.

-u i is the dependent variable, particle displacement of the solid frame.

-u f i is the dependent variable, filtration displacement, i.e. the relative particle displacement, with respect to the solid frame, for the pore fluid.

v i is the particle velocity of the solid frame.

v f i is the filtration velocity, i.e. the relative particle velocity, with respect to the solid frame, for the pore fluid,

i is the particle velocity of the pore fluid. -ϕ is the electric scalar potential.

ς is the variation (increment) of the fluid content, a measurement of the fluid amount flowing in/out of the porous media,

ϑ is the magnetic scalar potential.

x := (x, y, z) .

-x i is the spatial Cartesian coordinate with i = 1, 2, 3 corresponding to the x, y, z directions.

ξ is the magneto-electric matrix with size 3 × 3.

-Ξ is a diagonal matrix.

ζ is the electro-magnetic matrix with size 3 × 3.

For a function f with a real independent variable t, the Fourier transform is defined as (R → C)

for any real number ω, and the inverse Fourier transform is (C → R)

Based on these two definitions, we have the following immediate property:

For a function f with a real independent variable x in three-dimensional space, the Radon transform and its inverse are defined as

respectively, which correspond to a planar integral of the function f (x) over n • x = s and a spherical surface integral of the function f (s, n) over |n| = 1, respectively. Based on these two definitions, we have the following two immediate properties:

and

which transfer the multiple direction derivatives into the one-dimensional s-space derivative.

In figure 1, for a fixed x, usually a vector pointing from the source to the receiver, we can establish a new local coordinate, with the base (e, t 1 , t 2 ) [33, appendix A]. We have the x-aligned unit vector e defined as

The other two unit orthogonal vectors can be freely chosen. Equivalently, on the orthogonal plane, we can introduce the local azimuthal angle θ with unit vector t as

Then, the unit vector n is expressed as

Therefore, the spherical integral in §3b can be calculated as

Based on the Euler-Lagrange equation, we can derive the following wave equations for the poroelastic equation in time domain [55, p. 254]:

and

where '*' denotes a time convolution and (x, y, z, t) are the independent variables. Note that the summation convention is applied over subscript 'j'. According to §3a, we apply the Fourier transform over the t-axis, leading to the frequency-domain equations

and

with independent variables (x, y, z, ω) and dependent variables (u i , u f i ). Written in matrix form, the equations read [55, p. 282]

with

)

x , f

(1)

x , f

(2)

)

Note that (4.3) will be closed by relating the stress vector to the strain vector based on the constitutive relationship in §5. The pore fluid behaves differently in low-/high-frequency regimes. Specifically, for a lowfrequency regime (ω < ω c i ), the fluid has Poiseuille-type behaviour; therefore, we have

Note that ω 2 δ(ω) ≡ 0, so we can write

For a high-frequency regime (ω > ω c i ), the fluid is governed by the Johnson-Koplik-Dashen model [47, p. 27]; therefore, we have

with The detailed expressions are given below:

Note that the above equation is valid in both the time and frequency domains, given that the porous medium is purely elastic.

(b) Quality factor matrix and dispersive poroelastic media

In the case of dispersive media, Kjartansson [57] introduced the quality factor definition, based on the Caputo fractional derivative.

In time domain, the constitutive relationship (5.1) has to be amended: the modulus M , any entry of D in (5.1), involves a Caputo fractional derivative operator C 0 D 2γ t := ∂ 2γ /∂t 2γ as

In frequency domain, the constitutive relationship (5.1) is still valid, but with a complex-valued elasticity matrix: M becomes a much simpler complex value with a fractional exponent

Note the strain definition = L u. (6.1)

Then we have the governing equation in the frequency domain as

with the Kelvin-Christoffel matrix as

involving the derivative operators, only with the differential order 2. Therefore, the fundamental solution problem is written as [1, p. 44]

where I, P, K = 1-6. Note the summation convention is applied over the subscript 'P'. Now applying the Radon transform on both sides of (6.4), we get

However, Γ (n, ω) is still a full matrix, so it is not ready to solve the above equation.

Taking insight from solving the Riemann problem in the localized coordinates [21,35,45,[51][52][53], the eigenvalue diagonalization can be used to solve the fundamental solution for (6.5). For a matrix A, either symmetric or non-symmetric, either Hermitian or non-Hermitian, either complex or real, we have AR = RΛ, (

where R and Λ are, respectively, the right-eigenvector and eigenvalue matrices, assembled in a column manner. Similarly, we have the row-assembled left-eigenvector matrix U, subject to

Then, left multiplying (7.1) with U, and right multiplying (7.2) with R, respectively, lead to the second and third identities UAR = URΛ = ΛUR. (7.3) Therefore, UR commutes with a diagonal matrix Λ, that is, UR is also diagonal. This fact implies that the right eigenvector R is always invertible, even though A is a non-Hermitian matrix. Now we write (6.5) as

And we let A(n, ω) = ρ -1 (ω)Γ (n, ω), (7.5) which is a non-symmetric matrix, or more accurately, a non-Hermitian matrix, thus leading to a non-orthogonal eigenvector system. We note that the method proposed in [33] assumes the orthogonality of each individual eigenvector, thus this is not applicable here. We also note that Green's solution, provided in [

and

It is noteworthy that the analytical solution for (7.7) can be easily obtained row by row, inasmuch as Λ is a diagonal matrix. More specifically, we have

According to [45], the dispersive systems have four propagating modes (λ I > 0) and two nonpropagating modes (λ I = 0), respectively.

In the case of λ I > 0, the solution is

In the case of λ I = 0, (7.9) becomes ω 2 W IJ = -δ(s)F IJ , (7.11) and the solution is

Applying the inverse Radon transform leads to

Therefore, the non-propagating modes make no difference to the results.

Considering (7.8), (7.10) and the identity

we have the following expression in matrix form:

with the diagonal matrix

which can be decomposed into two parts as

where the static response is a circumferential integral,

and the response is a spherical integral,

Further applying the calculation method in §3c, we obtain (7.22) which can be easily calculated by applying the Gaussian quadrature rules [51]. Note that for the zero-eigenvalue inverse in Λ, we just let it be zero. To this end, the time-domain displacement expression is obtained,

Note the summation convention is applied over the subscript 'P'.

The above concise procedures invigorate the applications to other multi-physics problems, such as the piezoelectric-magnetic and thermo-elastic systems, as long as the governing equations share the unified form as (6.2). The linear magneto-electro-elasticity system, with dynamic elastic responses, has the following governing equations [54]:

and ρ 5×5 = diag ρI 3×3 , 0, 0 . (8.5)

Note the above can be reduced to the static piezoelectric problem or lossless poroelastic media [41,54]. Furthermore, we can also unifiedly get the Green's function for thermoelastic waves, by assembling the corresponding matrices [41].

Even though no analytical solution is available for the anisotropic dispersive poroelastic waves yet, we can still use the existing isotropic solutions, with just a special case of anisotropic media here, as a reference to verify our results.

From the perspective of the physical meaning, we need four Q values to describe the intrinsic attenuation of an isotropic poroelastic material: Q K and Q μ for the isotropic frame, corresponding to the attenuation of the bulk and shear moduli, respectively; Q K f for the pore fluid, with respect to the bulk modulus; and Q K s implicitly accounting for the solid-fluid coupling. Then the next step is to calculate the complex moduli, by substituting (5.4) into (5.5); the other entries of the matrix D in (5.2) can be obtained by the correspondence principle [55; 59, pp. 102, 291].

The source and receiver are at (0, 0, 0) m and (20, 250, 320) m, respectively. The source has a Ricker wavelet, with a central frequency f c = 60 Hz, a delayed time 1.2/f c and a polarization vector (0, 0, 1, 0, 0, 0) in (4.3). The reference frequency is f r = 10 Hz. Table 1 provides the isotropic poroelastic parameters. The reference solutions of the viscous poroelastic media are based on the extension of the existing isotropic solution for the purely poroelastic situation [60], according to §5b. Figure 2 provides the waveform verifications for different Q value sets in table 2, and the corresponding root mean square (RMS) differences are provided: the agreement is satisfactory, since almost all the RMS differences have a magnitude of 0.01%. The attenuation becomes stronger and the phases of the three wavelets, i.e. fast-P, S, slow-P waves, are left shifted, as the Q value decreases. Moreover, the wavelets are tremendously distorted, when strong attenuation, i.e. a small Q value set such as case 1 in table 2, is added. Table 1. The isotropic poroelastic parameters.

For the poroelastic materials with anisotropic attenuation, the implementation of Green's function is proposed for the first time; nevertheless, we can validate this analytical solution with an existing high-order discontinuous Galerkin (DG) algorithm [45].

The fluid-saturated medium used is an orthorhombic frame with the elastic constants shown in table 3. The supplementary porous material parameters refer to table 4. In table 5, we provide the Q factor values: a Q matrix with size 6 × 6 for anisotropic attenuation and dispersion in the frame, Q K f for the pore fluid and Q K s for the interaction between the solid skeleton and the pore fluid. According to (5.4) and (5.5), with the correspondence principle [55, pp. 102, 291; 59], we can obtain a complex-valued modulus matrix in (5.2) at the reference frequency f c = 10 Hz, Then, figure 3 displays the phase velocity distributions: (a) the real part, corresponding to the wave propagation, and (b) the imaginary part, corresponding to the wave dissipation, for the four types of waves, i.e. fast P, slow P, fast S and slow S waves, respectively [45, eqns. 26, 27]. In this figure, the anisotropy ratio is defined as

It is remarkable that the imaginary parts have distinct distributions versus the real parts in figure 3, which reveal different anisotropies for wave propagation and attenuation in this fluid-saturated material.

A source with a polarization vector (0, 0, 1, 0, 0, 0) in (4.3), having a Ricker wavelet with f c = 36 Hz and strength 1 × 10 15 N m -3 , is located at (0, 0, 0) m. Considering the symmetry feature of the velocity distributions in figure 3, four receivers are evenly deployed along a quarter circle with azimuthal angles, measured from the X-axis, {18 • , 36 , 72 • }, whose centre is at (0, 0, 600) m, and radius is 500 m, on the XY-plane. Figures 456provide the waveform comparisons between a high-order DG solver and the analytical solution proposed in this study, showing that excellent agreement is achieved. Owing to the anisotropy, the waveforms are overwhelmingly different between the different channels in figures 4-6.

In this study, we first propose a compact integral expression of Green's function for general poroelastic materials, incorporating anisotropic dispersion. Two critical techniques are leveraged: localization and diagonalization by applying the Radon transform and eigenvalue decomposition, respectively. Unlike the Fourier transform-based methods, the Radon transform method used in this paper only requires the integral over a finite area. Furthermore, this Radon transform brings the benefit that simplifies the three-dimensional spatial derivative into a onedimensional problem. Inspired by the procedure analysing the characteristics of hyperbolic equations, the originally intractable full anisotropic system is diagonalized, resulting in several simple scalar equations. The verification and validation demonstrate the correctness of the proposed Green's function. It is noteworthy that the analysis shown in this paper is heuristic for solving other pertinent problems, such as the fundamental solutions for piezoelectric-magnetic and thermo-elastic materials.