Generalized Thermo-poroelasticity Equations and Wave Simulation

Abstract

We establish a generalization of the thermoelasticity wave equation to the porous case, including the Lord–Shulman (LS) and Green–Lindsay (GL) theories that involve a set of relaxation times ($\tau _i, \ i = 1, \ldots , 4$). The dynamical equations predict four propagation modes, namely, a fast P wave, a Biot slow wave, a thermal wave, and a shear wave. The plane-wave analysis shows that the GL theory predicts a higher attenuation of the fast P wave, and consequently a higher velocity dispersion than the LS theory if $\tau _1 = \tau _2 > \tau _3$, whereas both models predict the same anelasticity for $\tau _1 = \tau _2 = \tau _3$. We also propose a generalization of the LS theory by applying two different Maxwell–Vernotte–Cattaneo relaxation times related to the temperature increment ($\tau _3$) and solid/fluid strain components ($\tau _4$), respectively. The generalization predicts positive quality factors when $\tau _4 \ge \tau _3$, and increasing $\tau _4$ further enhances the attenuation. The wavefields are computed with a direct meshing algorithm using the Fourier pseudospectral method to calculate the spatial derivatives and a first-order explicit Crank–Nicolson time-stepping method. The propagation illustrated with snapshots and waveforms at low and high frequencies is in agreement with the dispersion analysis. The study can be useful for a comprehensive understanding of wave propagation in high-temperature high-pressure fields.

Appendices

Appendix 1: Plane-wave analysis

We consider a 1D medium to analyze the phase velocity and attenuation of the different waves modes involved in the propagation, because the medium is isotropic. The S wave is not affected by the temperature effects, and its complex velocity is that of Biot theory (Carcione 2014):

$$\begin{aligned} v_c (\text{ S } \text{ wave}) = v_c = \sqrt{ \frac{\mu }{\rho - \rho _f^2 [ m - {\mathrm{i}}\eta / (\omega \kappa ) ]^{-1}} } . \end{aligned}$$

(17)

In 1D case, the field vector becomes $\mathbf{v} = [ v, q, \sigma , p , T]^\top$. Considering a plane wave of the form $\exp [ {\mathrm{i}}(\omega t - k x)]$, where $\omega$ is the angular frequency and k is the complex wavenumber, Eqs. (1), (3) and (5) in 1D case reduce to

$$\begin{aligned} \begin{array}{l} -k \sigma = \omega \rho v + \omega \rho _f q , \\ k p = \omega \rho _f v + \omega m q - ({\mathrm{i}}\eta / \kappa ) q , \\ \omega \sigma = -k E v - \alpha k M ( \alpha v + q) - \omega \beta {\bar{\tau }}_1 T, \\ \phi \omega p = \phi k M (\alpha v + q) + \omega \beta _f {\bar{\tau }}_2 T , \\ {\mathrm{i}}k T_0 \beta {\bar{\tau }}_4 (v+q) = [\gamma k^2 + {\mathrm{i}}\omega c {\bar{\tau }}_3] T , \\ {\bar{\tau }}_l = 1 + {\mathrm{i}}\omega \tau _l , \ \ l=1, \ldots , 4 , \ \ \ \end{array} \end{aligned}$$

(18)

where $E = \lambda + 2 \mu$.

This is a homogeneous system of linear equations whose solution is not zero if the determinant of matrix $\mathbf{A}$ is zero, whose components are

$$\begin{aligned} \begin{array}{l} a_{11} = \omega \rho , \ \ \ a_{12} = \omega \rho _f , \ \ \ a_{13} = k , \ \ \ a_{14} = 0 , \ \ \ a_{15} = 0, \\ a_{21} = \omega \rho _f, \ \ \ a_{22} = \omega m - {\mathrm{i}}\eta /\kappa , \ \ \ a_{23} = 0, \ \ \ a_{24} = - k, \ \ \ a_{25} = 0, \\ a_{31}= - k (E + \alpha ^2 M), \ \ \ a_{32} = - \alpha k M, \ \ \ a_{33} = - \omega , \ \ \ a_{34} = 0, \ \ \ a_{35} = - \omega \beta {\bar{\tau }}_1 , \\ a_{41} = \phi \alpha k M, \ \ \ a_{42} = \phi k M, \ \ \ a_{43} = 0 , \ \ \ a_{44} = - \omega \phi , \ \ \ a_{45} = \omega \beta _f {\bar{\tau }}_2, \\ a_{51} = - {\mathrm{i}}k T_0 \beta {\bar{\tau }}_4, \ \ \ a_{52} = a_{51}, \ \ \ a_{53} = 0, \ \ \ a_{54} = 0, \ \ \ a_{55} = {\mathrm{i}}\omega c {\bar{\tau }}_3 + \gamma k^2. \end{array} \end{aligned}$$

(19)

Based on it, we obtain the dispersion relation for P waves:

$$\begin{aligned} a_6 v_c^6 + a_4 v_c^4 + a_2 v_c^2 + a_0 = 0 , \end{aligned}$$

(20)

where

$$\begin{aligned} a_0&= {} {\mathrm{i}}\omega ^2 \phi \gamma M E , \nonumber \\ a_2&= {} - \omega \left\{ \phi \left[ b \gamma E_G + {\mathrm{i}}\omega \gamma ( m E_G + M (\rho -2 \alpha \rho _f)) + c M E {\bar{\tau }}_3 \right] \right. \nonumber \\&\quad \left. + \beta T_0 {\bar{\tau }}_4 \left[ \beta _f E {\bar{\tau }}_2 + (1-\alpha )M ( \phi \beta {\bar{\tau }}_1 - \alpha \beta _f {\bar{\tau }}_2) \right] \right\} , \nonumber \\ a_4&= {} b \phi \left[ \omega \rho \gamma - {\mathrm{i}}c E_G {\bar{\tau }}_3 - {\mathrm{i}}\beta ^2 T_0 {\bar{\tau }}_1 {\bar{\tau }}_4 \right] + \omega \left[ {\mathrm{i}}\omega \gamma \phi (m \rho -\rho _f^2) + c \phi {\bar{\tau }}_3 ( m E_G \right. \nonumber \\&\quad \left. + M (\rho -2 \alpha \rho _f)) + \beta T_0 {\bar{\tau }}_4 [ \beta \phi (m-\rho _f) {\bar{\tau }}_1 + \beta _f (\rho -\rho _f) {\bar{\tau }}_2 ] \right] , \nonumber \\ a_6&= {} c \phi {\bar{\tau }}_3 \left[ {\mathrm{i}}b \rho + \omega (\rho _f^2 - m \rho ) \right] , \nonumber \\ E_G&= {} E + \alpha ^2 M , \ \ \ b = \displaystyle \frac{\eta }{\kappa } . \end{aligned}$$

(21)

When $\beta = \beta _f$ = 0, we obtain a quadratic equation about $v_c$, corresponding to Biot velocities for the fast and slow P waves:

$$\begin{aligned} (-{\mathrm{i}}b \rho + \omega m \rho - \omega \rho _f^2) v_c^4 + ({\mathrm{i}}b E_G - \omega m E_G - \omega M \rho + 2 \omega \alpha M \rho _f) v_c^2 + \omega M E = 0 , \end{aligned}$$

(22)

and an additional root:

$$\begin{aligned} v_c = \sqrt{\frac{ {\mathrm{i}}\omega a^2}{{\bar{\tau }}_3}} , \ \ \ a = \sqrt{\frac{\gamma }{c}} , \end{aligned}$$

(23)

for the decoupled thermal wave, where a is the thermal diffusivity (Carcione et al. 2019a). It is evident that, at low frequencies, this velocity is zero.

The plane-wave analysis performed here is similar to the 1D time-periodic solutions obtained by Ignaczak and Ostoja-Starzewski (2010), whose complex wavenumbers are given by their Eq. (11.1.32), corresponding to the thermal and elastic waves. The phase velocity and attenuation factor can be obtained from the complex velocity as

$$\begin{aligned} v_p = \left[ {\mathrm{Re}} \left( v_c^{-1} \right) \right] ^{-1} \ \ \ {\mathrm{and}} \ \ \ Q = \dfrac{{\mathrm{Re}} (v_c^2)}{{\mathrm{Im}} (v_c^2)}, \end{aligned}$$

(24)

respectively (e.g., Carcione 2014).

Appendix 2: Crank–Nicolson explicit scheme

The Crank–Nicolson explicit scheme has been implemented by Carcione and Quiroga-Goode (1995) to solve the equations of poroelasticity and by Carcione et al. (2019b) to solve the thermoelasticity equations. The scheme, applied to the thermo-poroelasticity equations, is

$$\begin{aligned} D^{1/2} v_x&= {} \beta _{11} (\sigma _{xx,x} + \sigma _{xz,z} - f_x)^n - \beta _{12} p_{,x}^n - \displaystyle \frac{\eta }{\kappa } \beta _{12} A^{1/2}q_x = \varPi _x^n , \nonumber \\ D^{1/2} v_z&= {} \beta _{11} (\sigma _{xz,x} + \sigma _{zz,z} - f_z)^n - \beta _{12} p_{,z}^n - \displaystyle \frac{\eta }{\kappa } \beta _{12} A^{1/2}q_z = \varPi _z^n , \nonumber \\ D^{1/2} q_x&= {} \beta _{21} (\sigma _{xx,x} + \sigma _{xz,z} - f_x) ^n - \beta _{22} p_{,x}^n - \displaystyle \frac{\eta }{\kappa } \beta _{22} A^{1/2}q_x= \varOmega _x^n, \nonumber \\ D^{1/2} q_z&= {} \beta _{21} (\sigma _{xz,x} + \sigma _{zz,z} - f_z)^n - \beta _{22} p_{,z}^n - \displaystyle \frac{\eta }{\kappa } \beta _{22} A^{1/2}q_z = \varOmega _z^n, \nonumber \\ \epsilon _m&= {} (A^{1/2} v_x)_{,x} + ( A^{1/2} v_z)_{,z} , \nonumber \\ \epsilon _f&= {} ( A^{1/2} q_x)_{,x} + (A^{1/2} q_z)_{,z}, \nonumber \\ \epsilon&= {} \alpha \epsilon _m + \epsilon _f , \nonumber \\ {\dot{\epsilon }}_m&= {} (\varPi _x^n)_{,x} + (\varPi _z^n)_{,z}, \nonumber \\ {\dot{\epsilon }}_f&= {} (\varOmega _x^n)_{,x} + (\varOmega _z^n)_{,z}, \nonumber \\ \varDelta _\gamma T^n&= {} c ( A^{1/2} \psi + \tau _3 D^{1/2} \psi ) + \beta T_0 [ (\epsilon _m + \tau _4 {\dot{\epsilon }}_m) + (\epsilon _f + \tau _4 {\dot{\epsilon }}_f) ] + q^n , \nonumber \\ T^{n+1}&= {} T^{n} + dt \ \psi ^{n+1/2} , \nonumber \\ \varPi ^n&= {} (c \tau _3)^{-1} [\varDelta _\gamma T^n - q^n - \beta T_0 (\epsilon _m + \tau _4 ( \varPi _{x,x}^n + \varPi _{z,z}^n ) + \epsilon _f + \tau _4 ( \varOmega _{x,x}^n + \varOmega _{z,z}^n )) ] - \displaystyle \frac{1}{\tau _3} A^{1/2} \psi , \nonumber \\ D^1 \sigma _{xx}&= {} 2 \mu (A^{1/2} v_{x}),x + \lambda \epsilon _m + \alpha M \epsilon - \beta (A^{1/2} \psi + \tau _1 \varPi ^n) + f_{xx} , \nonumber \\ D^1 \sigma _{zz}&= {} 2 \mu (A^{1/2} v_{z}),z + \lambda \epsilon _m + \alpha M \epsilon - \beta (A^{1/2} \psi + \tau _1 \varPi ^n ) + f_{zz} , \nonumber \\ D^1 \sigma _{xz}&= {} \mu [( A^{1/2} v_{x}),z + (A^{1/2} v_{z}),x ) ] + f_{xz}, \nonumber \\ \phi D^1 p&= {} - \phi M \epsilon + \beta _f (A^{1/2} \psi + \tau _2 \varPi ^n) - f_f ,\nonumber \\ \end{aligned}$$

(25)

where

$$\begin{aligned} D^j \phi = \frac{ \phi ^{n+j} - \phi ^{n-j} }{2 j \mathrm{d}t} , \ \ \ \ \text{ and } \ \ \ \ A^j \phi = \frac{ \phi ^{n+j} + \phi ^{n-j} }{2} , \end{aligned}$$

(26)

are the central differences and mean value operators. In this three-level scheme, the particle velocities and $\psi$ at time $(n+1/2)\mathrm{d}t$ and stresses and temperature at time $(n+1)\mathrm{d}t$ are computed explicitly from particle velocities and $\psi$ at time $(n-1/2)\mathrm{d}t$, and stresses and temperature at time $(n-1) \mathrm{d}t$ and ndt, respectively.

For example, by expanding the third equation in (25), we obtain

$$\begin{aligned} q_x^{n+1/2}\left( 1+\dfrac{b}{2}\beta _{22}dt \right) =q_x^{n-1/2} \left( 1-\dfrac{b}{2}\beta _{22}dt \right) +dt[\beta _{21}(\sigma _{xx,x}+\sigma _{xz,z})^n-\beta _{22}p_{,x}^n]. \end{aligned}$$

(27)

Then,

$$\begin{aligned} v_x^{n+1/2}=v_x^{n-1/2}+dt\big [\beta _{11}(\sigma _{xx,x}+\sigma _{xz,z})^{n} -\beta _{12}p_{,x}^n-\beta _{12}\dfrac{\eta }{2k}(q_x^{n+1/2}+q_x^{n-1/2})\big ]. \end{aligned}$$

(28)

Similarly,

$$\begin{aligned} \sigma _{zz}^{n+1}&= {} \sigma _{zz}^{n-1}+2\mathrm{d}t\mu (v_z^{n+1/2}+v_z^{n-1/2})_{,z}+ 2\mathrm{d}t(\lambda +\alpha ^2 M)\epsilon _m+2\mathrm{d}t\alpha M\epsilon _f \nonumber \\&\quad -\,\mathrm{d}t\beta (\psi ^{n+1/2}+\psi ^{n-1/2})-2dt\beta \tau _1\varPi ^n.\nonumber \\ \end{aligned}$$

(29)

The equations for the other components can be similarly derived, based on equation (25). The numerical algorithms can be implemented on a staggered grid. For example, by defining $v_x$ and $q_x$ at coordinate $(x+d_x/2,z)$, $v_z$ and $q_z$ at coordinate $(x,z+d_z/2)$, $\sigma _{xx},\sigma _{zz}, \psi , T$ and p at $(x+d_x,z+d_z)$, and $\sigma _{xz}$ at $(x+d_x/2,z+d_z/2)$, the spatial derivatives in equation (25) can be obtained with the pseudospectral method as,

$$\begin{aligned} \begin{array}{l} M_{,x}={\mathrm{IFFT}} \big [{\mathrm{FFT}}(M){\mathrm{i}}k_x{\mathrm{exp}}({\mathrm{i}}k_x d_x/2)\big ],\ \ \ \ M = \sigma _{xx}, \ p, \ v_z,\\ M_{,z}= {\mathrm{IFFT}} \big [{\mathrm{FFT}}(M){\mathrm{i}}k_z{\mathrm{exp}}({\mathrm{i}}k_z d_z/2)\big ],\ \ \ \ M = \sigma _{zz}, \ p, \ v_x,\\ N_{,x}= {\mathrm{IFFT}}\big [{\mathrm{FFT}}(N){\mathrm{i}}k_x{\mathrm{exp}}(-{\mathrm{i}}k_x d_x/2)\big ],\ \ \ \ N = \sigma _{xz}, \ v_x, \ q_x,\\ N_{,z}= {\mathrm{IFFT}} \big [{\mathrm{FFT}}(N){\mathrm{i}}k_z{\mathrm{exp}}(-{\mathrm{i}}k_z d_z/2)\big ],\ \ \ \ N = \sigma _{xz}, \ v_z, \ q_z,\\ T_{,zz}={\mathrm{IFFT}}\big [- k_z^2 {\mathrm{FFT}}(T)\big ],\\ T_{,xx}={\mathrm{IFFT}}\big [- k_x^2{\mathrm{FFT}}(T)\big ],\\ \end{array} \end{aligned}$$

(30)

where FFT and IFFT are the forward and backward fast Fourier transforms, respectively, and ${k_x}$ and ${k_z}$ are the wavenumbers along x and z directions, respectively. A list of the algorithm in pseudo-code form is given in Table 2.

Table 2 Crank–Nicolson scheme for the thermo-poroelasticity equations

The stability analysis for similar equations has been studied in Carcione and Quiroga-Goode (1995), with a Von Neumann stability analysis based on the eigenvalues of the amplification matrix. The algorithm has first-order accuracy but possesses the stability properties of implicit algorithms and the solution can be obtained explicitly. Alternatively, implicit methods can also be used when the differential equations are stiff, but are more cumbersome to implement than explicit methods. The instability is mainly due to the presence of the quasi-static mode (the Biot slow wave). While the eigenvalues of the fast waves have a small real part, the eigenvalue of the Biot wave (in the quasi-static regime) has a large real part. Indeed, the simple explicit first-order accurate Crank–Nicolson scheme proposed here provides an efficient scheme to deal with the stiffness, as shown in Carcione and Quiroga-Goode (1995).