Three-Temperature Boundary Element Modeling of Ultrasound Wave Propagation in Anisotropic Viscoelastic Porous Media

Abstract

The main goal of this work is to develop a novel boundary element method (BEM) model for analyzing ultrasonic wave propagation in three-temperature anisotropic viscoelastic porous media. Due to the problems of the strong nonlinearity of ultrasonic wave propagation in three-temperature porous media, the analytical or numerical solutions to the problems under consideration are always challenging, necessitating the development of new computational techniques. As a result, we use a new BEM model to solve such problems. A time-stepping procedure based on the linear multistep method is obtained after solving the discretized boundary integral equation with the quadrature rule. The calculation of a double integral is required to obtain fundamental solutions, but this increases the total BEM computation time. Our proposed BEM technique is used to solve the current problem and improve the formulation efficiency. The numerical results are graphed to demonstrate the effects of viscosity and anisotropy on the nonlinear ultrasonic stress waves in three-temperature porous media. The validity, accuracy, and efficiency of the proposed methodology are demonstrated by comparing the obtained results to a corresponding solution obtained from the finite difference method (FDM).

1. Introduction

In the numerical study of laser-driven implosion in inertial confinement fusion experiments, two-dimensional three-temperature radiation diffusion equations coupled with photon, electron, and ion temperatures are extensively utilized. The equations describe the evolution of electron, ion, and photon radiation in diverse materials. As a result, they exhibit a high degree of nonlinearity and are based on complex computational domains. Meshes are used for discretization in traditional numerical methods, such as the finite difference method and the finite element method.

In recent years, many engineering studies have drawn the attention of researchers to investigate the mechanical behavior of thermoelastic materials [1,2,3,4,5,6] due to the positive results obtained in applied science, engineering, and technological applications such as biology, biophysics, biomechanics, geotechnical engineering, reservoir geomechanics, mining and petroleum engineering, geothermal engineering, thermal insulation, and lightweight structural design. Some researchers have investigated the effects of the magnetic field [7], initial stress [8], rotation, and gravity [9] on the generalized thermo-viscoelastic diffusion medium. In addition, the effects of initial stress and temperature dependence on the thermo-microstretch elastic solid have been investigated with a dual-phase-lag model [10]. Xu et al. [11] studied the effects of viscoelastic dampers with high levels of energy dissipation based on an acrylate rubber matrix. As analytical solutions to the current situation are extremely difficult to achieve, numerical methods have emerged as the primary tool for resolving these problems, such as in Pei et al. [12], Ooi et al. [13], Zhou et al. [14], Ng et al. [15], Majchrzak and Turchan [16], Bottauscio et al. [17], Deng and Liu [18], and Partridge and Wrobel [19]. The concept of treating biological tissue as a porous medium has been deemed more acceptable for incorporating blood flow through arteries implanted in the tissue. One of the computational strategies used to solve the bioheat transfer problems of biological tissues [20,21,22] is the boundary element method (BEM) [23,24,25,26,27,28,29,30,31]. For physical and technological problems, Laplace-domain fundamental solutions are generally easier to achieve than time-domain fundamental solutions [32,33]. As a result, because it requires the Laplace-domain fundamental solutions of the problem’s governing equations, the CQBEM is highly useful for problems that did not have time-domain fundamental solutions. As a result, CQBEM broadens the spectrum of engineering problems that can be tackled using the traditional time-domain BEM.

The primary goal of this paper is to present a new boundary element model for explaining thermomechanical interactions in three-temperature anisotropic viscoelastic porous media. The uncoupled governing equations are solved independently: the bioheat equation is solved first using the GBEM based on LRBFCM to obtain the temperature distribution, and then the mechanical equation is solved using the CQBEM to obtain the displacements and stresses. The resulting linear systems have been solved by a communication-avoiding Arnoldi (CA-Arnoldi) preconditioner, which reduces the number of iterations and the total CPU time. The numerical results demonstrate the validity, efficiency, and accuracy of the proposed model.

2. Formulation of the Problem

In the Cartesian system $\left(x, y, z\right)$, we consider a region $\mathsf{\Omega }=\left\{0<x<\underset{¯}{\alpha ,}, 0<y<\underset{¯}{\beta ,}0<z<\underset{¯}{\gamma ,}\right\}$ with a boundary $\Gamma$ occupied by an anisotropic viscoelastic porous media, as shown in Figure 1.

According to Biot’s model [34] and Darcy’s law [35], the thermo-poroelastic governing equations can be expressed as:

$${\left({\nabla }^T\sigma \right)}^T+\mathrm{F}=\rho \ddot{\mathrm{u}}+\varphi {\rho }_f\left({\ddot{\mathrm{u}}}_f-\ddot{\mathrm{u}}\right)$$

(1)

$$\dot{\zeta }+{\nabla }^T\mathrm{q}=0$$

(2)

$$\sigma =\left(C_{ajlg}\left(\tau \right) \chi \mathrm{tr} ϵ-A\overline{p}\right)\mathrm{I}-\mathfrak{B} \theta$$

(3)

$$ϵ=\frac{1}{2}(\nabla {\mathrm{u}}^T+{(\nabla {\mathrm{u}}^T)}^T)$$

(4)

$$\zeta =A \mathrm{tr} ϵ+\frac{{\varphi }^2}{R}\overline{p}$$

(5)

$$\mathrm{q}=-K\left(\nabla \overline{p}+{\rho }_f\ddot{\mathrm{u}}+\frac{{\rho }_e+\varphi {\rho }_f}{\varphi }\left({\ddot{\mathrm{u}}}_f-\ddot{\mathrm{u}}\right)\right)$$

(6)

On the basis of Bonnet [36], the governing equations can be written as follows [37]:

$$\left.\begin{array}{c}{\widehat{B}}_{\tilde{x}}{\hat{\mathrm{u}}}^g\left(\tilde{\mathrm{x}}\right)=0 \mathrm{for} \tilde{x}\in \mathsf{\Omega } \\ {\hat{\mathrm{u}}}^g\left(\mathrm{x}\right)={\widehat{g}}_D \mathrm{for} x\in {\Gamma }_D\\ {\hat{\mathrm{t}}}^g\left(\mathrm{x}\right)={\widehat{g}}_N \mathrm{for} x\in {\Gamma }_N\end{array}\right\}$$

(7)

in which the operator ${\widehat{B}}_{\tilde{\mathrm{x}}}$ and the tractions ${\widehat{t}}^g$ are defined as:

$$\begin{gathered} {\widehat{B}}_{\tilde{\mathrm{x}}}=\left[\begin{array}{ccc}{B}_{\tilde{\mathrm{x}}}^e+s^2\left(\rho -\beta {\rho }_f\right)I& \left(\alpha -\beta \right){\nabla }_{\tilde{\mathrm{x}}}& -\mathfrak{B}{\nabla }_{\tilde{\mathrm{x}}}\\ s\left(\alpha -\beta \right){\nabla }_{\tilde{\mathrm{x}}}^T& -\frac{\beta }{s{\rho }_f}{\Delta }_{\tilde{\mathrm{x}}}+\frac{s{\varphi }^2}{R}& 0\end{array}\right], \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\hat{\mathrm{t}}}^g\left(x\right)=\left[\begin{array}{ccc}{T}_{\mathrm{x}}^e& -\alpha {\mathrm{n}}_{\mathrm{x}}& 0\\ s\beta {\mathrm{n}}_{\mathrm{x}}^T& \frac{\beta }{s{\rho }_f}{\mathrm{n}}_{\mathrm{x}}^T{\nabla }_{\mathrm{x}}& 0\end{array}\right]\left[\begin{array}{c}\hat{\mathrm{u}}\left(\mathrm{x}\right)\\ \widehat{p}\left(\mathrm{x}\right)\\ \theta \left(\mathrm{x}\right)\end{array}\right], \beta =\frac{{\varphi }^{2}sK{\rho }_f}{{\varphi }^2+sK\left({\rho }_e+\varphi {\rho }_f\right)} \end{gathered}$$

(8)

According to Fahmy [30], the thermomechanical interactions can be found by treating the soft tissue as a thermo-poroelastic medium and implementing the BEM for solving the governing Equations (1) and (7).

The three-temperature radiative diffusion equations are as follows:

$$C_{ve}\frac{\partial {\theta }_e\left(\mathrm{r}, \tau \right)}{\partial \tau }-\frac{1}{\rho }\nabla \left[{\mathbb{K}}_e \nabla {\theta }_e\left(\mathrm{r}, \tau \right)\right]=-{\mathbb{W}}_{ei} \left({\theta }_e-{\theta }_i\right)-{\mathbb{W}}_{er} \left({\theta }_e-{\theta }_r\right)$$

(9a)

$$C_{vi}\frac{\partial {\theta }_i\left(\mathrm{r}, \tau \right)}{\partial \tau }-\frac{1}{\rho }\nabla \left[{\mathbb{K}}_i \nabla {\theta }_i\left(\mathrm{r}, \tau \right)\right]={\mathbb{W}}_{ei} \left({\theta }_e-{\theta }_i\right)$$

(9b)

$$C_{vr}\frac{\partial {\theta }_r\left(\mathrm{r}, \tau \right)}{\partial \tau }-\frac{1}{\rho }\nabla \left[{\mathbb{K}}_r \nabla {\theta }_r\left(\mathrm{r}, \tau \right)\right]={\mathbb{W}}_{er} \left({\theta }_e-{\theta }_r\right)$$

(9c)

in which

$$C_{v\alpha }=\left\{\begin{array}{cc}{C}_{v\alpha }=C_e\hfill & \alpha =e\\ C_{v\alpha }=C_i\hfill & \alpha =i\\ C_{v\alpha }=C_r{T}_r^3\hfill & \alpha =r\end{array}\right.\mathrm{and} {\mathbb{K}}_{\alpha }=\left\{\begin{array}{cc}{\mathbb{A}}_e T_e^{5/2}\hfill & \alpha =e\\ {\mathbb{A}}_i T_i^{5/2}\hfill & \alpha =i\\ {\mathbb{A}}_r T_r^{3+\mathbb{B}}\hfill & \alpha =r\end{array}\right.$$

where $e, i, \mathrm{and} r$ denote the electron, ion, and phonon, respectively.

3. Boundary Element Implementation for the Temperature Field

The two-dimensional (2D) three-temperature (3T) radiation diffusion Equations (9a)–(9c) can be written as:

$$\nabla \left[{\mathbb{K}}_{\alpha } \nabla {\theta }_{\alpha }\left(\mathrm{r}, \tau \right)\right]+\overline{\mathbb{W}}\left(\mathrm{r}, \tau \right)=c_{\alpha }\rho {\delta }_1\frac{\partial {\theta }_{\alpha }\left(\mathrm{r}, \tau \right)}{\partial \tau }+Q\left(\mathrm{r}, \tau \right)$$

(10)

in which

$$Q\left(\mathrm{r}, \tau \right)=\frac{1-R}{x_0}{e}^{\left(-\frac{r_a}{r_0}\right)J\left(\tau \right)}, J\left(\tau \right)=\frac{J_0 \tau }{{\tau }_1^2}{e}^{-\frac{\tau }{{\tau }_1}}, a=1, 2, 3.$$

where

$$\overline{\mathbb{W}}\left(r, \tau \right)=\left\{\begin{array}{c}-\rho {\mathbb{W}}_{ei} \left({\theta }_e-{\theta }_i\right) –\rho {\mathbb{W}}_{er} \left({\theta }_e-{\theta }_r\right), \alpha =e, {\delta }_1=1 \\ \rho {\mathbb{W}}_{ei} \left({\theta }_e-{\theta }_i\right), \alpha =i, {\delta }_1=1 \\ \rho {\mathbb{W}}_{er} \left({\theta }_e-{\theta }_r\right), \alpha =r, {\delta }_1=T_p^3 \end{array}\right.$$

The total energy per unit mass is as follows:

$$P=P_e+P_i+P_r, P_e=c_e{\theta }_e, P_i=c_i{\theta }_i,P_r=\frac{1}{4}{c}_r{\theta }_r^4$$

(11)

The conditions under consideration can be summarized as follows:

$${\theta }_{\alpha }\left(x, y, 0\right)={\theta }_{\alpha }^0\left(x, y\right)=g_1\left(x, \tau \right)$$

(12)

$${\left.{\mathbb{K}}_{\alpha }\frac{\partial {\theta }_{\alpha }}{\partial n}\right|}_{{\Gamma }_1}=0, \alpha =e, i, {\left.{\theta }_p\right|}_{{\Gamma }_1}=g_2\left(x, \tau \right)$$

(13)

$${\left.{\mathbb{K}}_{\alpha }\frac{\partial {\theta }_{\alpha }}{\partial n}\right|}_{{\Gamma }_2}=0, \alpha =e, i, r$$

(14)

Using the fundamental solution to the following differential equation:

$$D{\nabla }^2{\theta }_{\alpha }+\frac{\partial {\theta }_{\alpha }^{*}}{\partial n}=-\delta \left(\mathrm{r}-p_i\right)\delta \left(\tau -\mathrm{r}\right), D=\frac{{\mathbb{K}}_{\alpha }}{\rho c}$$

(15)

The dual-reciprocity boundary integral equation corresponding to (10) can be expressed, as in Fahmy [17], as follows:

$${\theta }_{\alpha }=\frac{D}{{\mathbb{K}}_{\alpha }}{{\displaystyle \int }}_O^{\tau }{{\displaystyle \int }}_S^{}\left[{\theta }_{\alpha }{q}^{* }-{\theta }_{\alpha }^{*}q\right]dS d\tau +\frac{D}{{\mathbb{K}}_{\alpha }}{{\displaystyle \int }}_O^{\tau }{{\displaystyle \int }}_R^{}b{\theta }_{\alpha }^{*}dR d\tau +{{\displaystyle \int }}_R^{}{\left.{\theta }_{\alpha }^i{\theta }_{\alpha }^{*}\right|}_{\tau =0}dR$$

(16)

which can be expressed as follows:

$$C{\theta }_{\alpha }={{\displaystyle \int }}_S^{}\left[{\theta }_{\alpha }{q}^{*}-{\theta }_{\alpha }^{*}q\right] dS-{{\displaystyle \int }}_R^{}\frac{{\mathbb{K}}_{\alpha }}{D} \frac{\partial {\theta }_{\alpha }^{*}}{\partial \tau }{\theta }_{\alpha } dR$$

(17)

We assume that the temperature derivative in (17) is approximated as:

$$\frac{\partial {\theta }_{\alpha }}{\partial \tau }\cong {\displaystyle \sum }_{j=1}^N{f}^j{\left(r\right)}^j{a}^j\left(\tau \right)$$

(18)

Now, we consider:

$${\nabla }^2{\widehat{\theta }}_{\alpha }^j=f^j$$

(19)

Thus, from Equation (17), we obtain:

$$C \theta ={{\displaystyle \int }}_S^{}\left[{\theta }_{\alpha }{q}^{*}-{\theta }_{\alpha }^{*}q\right] dS+{\displaystyle \sum }_{j=1}^N{a}^j\left(\tau \right)D^{-1}\left(C{\widehat{\theta }}_{\alpha }^j-{{\displaystyle \int }}_S^{}\left[{\theta }_{\alpha }^j{q}^{*}-{\widehat{q}}^j{\theta }_{\alpha }^{*}\right] dS\right)$$

(20)

where

$${\widehat{q}}^j=-{\mathbb{K}}_{\alpha }\frac{\partial {\widehat{\theta }}_{\alpha }^j}{\partial n}$$

(21)

and

$$a^j\left(\tau \right)={\displaystyle \sum }_{i=1}^N{f}_{ji}^{-1}\frac{\partial \theta \left(r_i, \tau \right)}{\partial \tau }$$

(22)

where $f_{ji}^{-1}$ are defined as:

$${\left\{F\right\}}_{ji}=f^j\left(r_i\right)$$

(23)

By using Equations (20) and (22), we obtain:

$$C {\dot{\theta }}_{\alpha }+H {\theta }_{\alpha }=G Q$$

(24)

where

$$C=-\left[H {\widehat{\theta }}_{\alpha }-G \widehat{Q}\right]F^{-1}{D}^{-1}$$

(25)

with

$${\left\{\widehat{\theta }\right\}}_{ij}={\widehat{\theta }}^j\left(x_i\right)$$

(26)

$${\left\{\widehat{Q}\right\}}_{ij}={\widehat{q}}^j\left(x_i\right)$$

(27)

Now, we introduce the following functions:

$${\theta }_{\alpha }=\left(1-\mathsf{\theta }\right) {\theta }_{\alpha }^m+\mathsf{\theta } {\theta }_{\alpha }^{m+1}$$

(28)

$$q=\left(1-\mathsf{\theta }\right)q^m+\mathsf{\theta } q^{m+1}$$

(29)

where $0\le \mathsf{\theta }=\frac{\tau -{\tau }^m}{\Delta {\tau }^m} \le 1$, $\Delta {\tau }^m={\tau }^{m+1}-{\tau }^m$.

By differentiating (28), we obtain:

$${\dot{\theta }}_{\alpha }=\frac{d{\theta }_{\alpha }}{d\mathsf{\theta }}\frac{d\mathsf{\theta }}{d\tau }=\frac{{\theta }_{\alpha }^{m+1}-{\theta }_{\alpha }^m}{\Delta {\tau }^m}$$

(30)

Substitution of Equations (31)–(33) into Equation (27), yields:

$$\left(\frac{C}{\Delta {\tau }^m}+\mathsf{\theta }H\right){\theta }_{\alpha }^{m+1}-\mathsf{\theta }G{Q}^{m+1}=\left(\frac{C}{\Delta {\tau }^m}-\left(1-\mathsf{\theta }\right)H\right){\theta }_{\alpha }^m+\left(1-\mathsf{\theta }\right)G{Q}^m$$

(31)

which can be written as:

$$\mathrm{𝕒}\mathsf{{\rm X}}=\mathrm{𝕓}$$

(32)

To solve the resulting linear algebraic systems, the symmetric successive over-relaxation (SSOR) method without matrix inversion [38] is efficiently implemented.

4. Boundary Element Implementation for the Poroelastic Fields

The representation formula for problem (7) is as follows:

$${\widehat{u}}^g\left(\tilde{x}\right)={(\widehat{V}{\widehat{t}}^g)}_{\mathsf{\Omega }}\left(\tilde{x}\right)-{(\widehat{K}{\widehat{u}}^g)}_{\mathsf{\Omega }}\left(\tilde{x}\right) \mathrm{for} \tilde{x}\in \mathbf{\Omega }$$

(33)

where the integral operators are:

$${(\widehat{V}{\widehat{t}}^g)}_{\mathsf{\Omega }}\left(\tilde{x}\right)=\underset{\Gamma }{\overset{.}{{\displaystyle \int }}}{\widehat{U}}^T\left(y-\tilde{x}\right){\widehat{t}}^g\left(y\right)d{s}_y$$

(34)

$${(\widehat{K}{\widehat{u}}^g)}_{\mathsf{\Omega }}\left(\tilde{x}\right)=\underset{\Gamma }{\overset{.}{{\displaystyle \int }}}{({\widehat{T}}_y\widehat{U})}^T\left(y-\tilde{x}\right){\widehat{u}}^g\left(y\right) d{s}_y$$

(35)

In the Laplace domain, the fundamental solution and associated traction are denoted as [9]:

$$\widehat{U}\left(\mathrm{r}\right)=\left[\begin{array}{ccc}{\widehat{U}}^s\left(\mathrm{r}\right)& {\widehat{U}}^f\left(\mathrm{r}\right)& 0\\ {({\widehat{P}}^s)}^T\left(\mathrm{r}\right)& {\widehat{P}}^f\left(\mathrm{r}\right)& 0\end{array}\right], {\widehat{T}}_y=\left[\begin{array}{ccc}{T}_y^e& s\alpha n_y& 0\\ -\beta n_y^T& \frac{\beta }{s{\rho }^f}{n}_y^T\nabla & 0\end{array}\right] \mathrm{with} \mathrm{r}≔\left|y-x\right|$$

(36)

The fundamental solution can be expressed as [35]:

$${\widehat{U}}^s\left(\mathrm{r}\right)=\frac{1}{4\pi \mathrm{r}\left(\rho -\beta {\rho }^f\right)}\left[{ℝ}_1\frac{\left(k_4^2-k_2^2\right)}{\left(k_1^2-k_2^2\right)}{e}^{-k_1\mathrm{r}}-{ℝ}_2\frac{\left(k_4^2-k_1^2\right)}{\left(k_1^2-k_2^2\right)}{e}^{-k_2\mathrm{r}}+\left(I{k}_3^2-{ℝ}_3\right)e^{-k_3\mathrm{r}}\right]$$

(37)

where

$${ℝ}_j=\frac{3{\nabla }_y\mathrm{r}{\nabla }_y^T\mathrm{r}-I}{{\mathrm{r}}^2}+k_j\frac{3{\nabla }_y\mathrm{r}{\nabla }_y^T\mathrm{r}-I}{\mathrm{r}}+k_j^2{\nabla }_y\mathrm{r}{\nabla }_y^T\mathrm{r}$$

(38)

Equation (37) can be expressed as:

$${\widehat{U}}^s\left(\mathrm{r}\right)=\frac{1}{4\pi \mu \mathrm{r}\left(\lambda +2\mu \right)}\left[\left(\lambda +\mu \right){\nabla }_{y}r{\nabla }_y^T\mathrm{r}+I\left(\lambda +3\mu \right)\right]+O\left({\mathrm{r}}^0\right)$$

(39)

The fundamental solution can be expressed as:

$$\begin{gathered} {\widehat{U}}^s\left(\mathrm{r}\right)={\widehat{U}}_s^s\left(\mathrm{r}\right)+{\widehat{U}}_r^s\left(\mathrm{r}\right) \\ =\frac{1}{\mu }\left[I{\Delta }_y-\frac{\lambda +\mu }{\lambda +2\mu }{\nabla }_y{\nabla }_y^T\right]{\Delta }_y\widehat{x}\left(\mathrm{r}\right) \\ -\frac{1}{\mu }\left[\left(\left(k_1^2+k_2^2\right){\Delta }_y-k_1^2{k}_2^2\right)I-\left(k_1^2+k_2^2-k_4^2-\frac{k_1^2{k}_2^2}{k_3^2}\right){\nabla }_y{\nabla }_y^T\right]\widehat{x}\left(\mathrm{r}\right) \end{gathered}$$

(40)

in which

$$\begin{gathered} \widehat{x}\left(\mathrm{r}\right)=\frac{1}{4\pi r}\left[\frac{e^{-k_1\mathrm{r}}}{\left(k_2^2-k_1^2\right)\left(k_3^2-k_1^2\right)}+\frac{e^{-k_2\mathrm{r}}}{\left(k_2^2-k_1^2\right)\left(k_2^2-k_3^2\right)}+\frac{e^{-k_3\mathrm{r}}}{\left(k_1^2-k_3^2\right)\left(k_2^2-k_3^2\right)}\right] \\ =-\frac{1}{\left(k_1^2-k_2^2\right)\left(k_1^2-k_3^2\right)\left(k_3^2-k_2^2\right)}+O\left({\mathrm{r}}^2\right) \end{gathered}$$

(41)

Furthermore, the remaining components of the fundamental solution might be stated as:

$${\widehat{U}}^f\left(\mathrm{r}\right)=\frac{{\rho }^f\left(\alpha -\beta \right){\nabla }_y\mathrm{r}}{4\pi r\beta \left(\lambda +2\mu \right)\left(k_1^2-k_2^2\right)}\left[\left(k_1+\frac{1}{\mathrm{r}}\right)e^{-k_1\mathrm{r}}-\left(k_2+\frac{1}{\mathrm{r}}\right)e^{-k_2\mathrm{r}}\right]=O\left({\mathrm{r}}^0\right)$$

(42)

$${\widehat{P}}^s\left(\mathrm{r}\right)=\frac{{\widehat{U}}^f\left(\mathrm{r}\right)}{s}=O\left({\mathrm{r}}^0\right)$$

(43)

$${\widehat{P}}^f\left(\mathrm{r}\right)=\frac{s{\rho }^f}{4\pi r\beta \left(k_1^2-k_2^2\right)}\left[\left(k_1^2-k_4^2\right)e^{-k_{1}r}-\left(k_2^2-k_4^2\right)e^{-k_{2}r}\right]=\frac{s{\rho }^f}{4\pi r\beta }+O\left({\mathrm{r}}^0\right)$$

(44)

Now, we apply the following limiting $\tilde{x}\in \mathsf{\Omega }\to x\in \Gamma$ to (34) to obtain:

$$\underset{\tilde{x}\in \mathsf{\Omega }\to x\in \Gamma }{\lim }{\left(\widehat{V}{\widehat{t}}^g\right)}_{\mathsf{\Omega }}\left(\tilde{x}\right)=\left(\widehat{V}{\widehat{x}}^g\right)\left(x\right)≔\underset{\Gamma }{\overset{.}{{\displaystyle \int }}}{\widehat{U}}^T\left(y-x\right){\widehat{t}}^g\left(y\right)d{s}_y$$

(45)

In addition, we apply the following limiting method to (35) to obtain [39]:

$$\underset{\tilde{x}\in \mathsf{\Omega }\to x\in \Gamma }{\lim }{\left(\widehat{K}{\widehat{u}}^g\right)}_{\mathsf{\Omega }}\left(\tilde{x}\right)=\left[-I\left(x\right)+C\left(x\right)\right]{\widehat{u}}^g\left(x\right)+\left(\widehat{K}{\widehat{u}}^g\right)\left(x\right)$$

(46)

in which

$$C\left(x\right)=\underset{\epsilon \to 0}{\lim }\underset{y\in \mathsf{\Omega }:\left|y-x\right|=\epsilon }{\overset{.}{{\displaystyle \int }}}{\left({\widehat{T}}_y\widehat{U}\right)}^T\left(y-x\right)d{s}_y$$

(47)

and

$$\left(\widehat{K}{\widehat{u}}^g\right)\left(x\right)=\underset{\epsilon \to 0}{\lim }\underset{\left|y-x\right|\ge \epsilon }{\overset{.}{{\displaystyle \int }}}{\left({\widehat{T}}_y\widehat{U}\right)}^T\left(y-x\right){\widehat{u}}^g\left(y\right)d{s}_y$$

(48)

By using Equations (45)–(48), the boundary integral equation in Laplace domain can be expressed as:

$$C\left(x\right) {\widehat{u}}^g \left(x\right)=\left(\widehat{V}{\widehat{t}}^g\right)\left(x\right)-\left(\widehat{K}{\widehat{u}}^g\right)\left(x\right)$$

(49)

The poro-elastodynamic boundary integral equation can be expressed using the inverse Laplace transformation as:

$$C\left(x\right)u^g\left(x,t\right)=\left(V\ast t^g\right)\left(x,t\right)-\left(K{u}^g\right)\left(x,t\right)$$

(50)

The fundamental solution is as follows [37]:

$${\left({\widehat{T}}_y\widehat{U}\right)}^T={\left[\left[\begin{array}{cc}{\widehat{T}}_y^e& s\alpha n_y\\ -\beta n_y^T& \frac{\beta }{s{\rho }_0^f}{n}_y^T{\nabla }_y\end{array}\right]\left[\begin{array}{cc}{\widehat{U}}^s& {\widehat{U}}^f\\ {\left({\widehat{P}}^s\right)}^T& {\widehat{P}}^f\end{array}\right]\right]}^T={\left[\begin{array}{cc}{\widehat{T}}^s& {\widehat{T}}^f\\ {\left({\widehat{Q}}^s\right)}^T& {\widehat{Q}}^f\end{array}\right]}^T$$

(51)

The Stokes theorem states that the differentiable vector field $\mathrm{a}\left(y\right)$, $\left(y\in \Gamma \right)$ can be represented as:

$$\underset{\Gamma }{\overset{.}{{\displaystyle \int }}}\left({\nabla }_y\times \mathrm{a}, n_y \right)d{s}_y=-\underset{\partial \Gamma }{\overset{,}{{\displaystyle \int }}}\left(\mathrm{a}, v\right)d{\gamma }_y=-\underset{\varphi }{\overset{,}{{\displaystyle \int }}}\left(\mathrm{a}, v\right)d{\gamma }_y=0$$

(52)

where

$$\underset{\Gamma }{\overset{.}{{\displaystyle \int }}}\left(n_y\times {\nabla }_y, \mathrm{a} \right)d{s}_y=0$$

(53)

We can use (53) to obtain the following formula:

$$\underset{\Gamma }{\overset{.}{{\displaystyle \int }}}\left(M_y \mathrm{a}\right)d{s}_y=0,\hspace{1em}{M}_y={\left({\nabla }_y{\nabla }_y^T\right)}^T-{\nabla }_y{\nabla }_y^T,\hspace{1em}\mathrm{a}=vu$$

(54)

According to [40], we obtain

$$\underset{\Gamma }{\overset{.}{{\displaystyle \int }}}\left(M_y v\right)ud{s}_y=-\underset{\Gamma }{\overset{,}{{\displaystyle \int }}}v\left(M_y u\right)d{s}_y$$

(55)

$$\underset{\Gamma }{\overset{.}{{\displaystyle \int }}}{\left(M_y v\right)}^{T}ud{s}_y=\underset{\Gamma }{\overset{.}{{\displaystyle \int }}}{v}^T\left(M_y u\right)d{s}_y$$

(56)

By using (40) and (51), we obtain:

$${\left({\widehat{T}}^s\right)}^T={\left(T_y^e\left({\widehat{U}}_{sing}^s+{\widehat{U}}_{reg}^s\right)\right)}^T+s\alpha {\widehat{P}}^s{n}_y^T={\left(T_y^e {\widehat{U}}_{sing}^s\right)}^T+O\left({\mathrm{r}}^0\right)$$

(57)

According to [37], we obtain:

$${\left({\widehat{T}}^s\right)}^T=\left(\lambda +2\mu \right)n_y{\nabla }_y^T{\widehat{U}}_{sing}^s-\mu \left(n_y\times \left({\nabla }_y\times {\widehat{U}}_{sing}^s\right)\right)+2\mu M_y{\widehat{U}}_{sing}^s+o\left({\mathrm{r}}^0\right)$$

(58)

which can be expressed using (40) as:

$${\left({\widehat{T}}^s\right)}^T=M_y{\Delta }_y^2\widehat{X}+I\left(n^T{\nabla }_y\right){\Delta }_y^2\widehat{X}+2\mu {\left(M_y{\widehat{U}}_{sing}^s\right)}^T+o\left({\mathrm{r}}^0\right)$$

(59)

Using (35) and (59), we obtain:

$${\left(\widehat{k}\widehat{u}\right)}_{\Omega }^s\left(\tilde{x}\right)={{\displaystyle \int }}_{\Gamma }^{.}\left[\left(M_y{\Delta }_y^2\widehat{X}\right)\widehat{u}+\left(I\left(n^T{\nabla }_y\right){\Delta }_y^2\widehat{X}\right)\widehat{u}+2\mu {\left(M_y{\widehat{U}}_{sing}^s\right)}^T\widehat{u}+O\left({\mathrm{r}}^0\right)\widehat{u}\right]d{s}_y$$

(60)

Based on [39], we obtain:

$${\left(\widehat{K}\widehat{u}\right)}_{\Omega }^s\left(\tilde{x}\right)=\underset{\Gamma }{\overset{.}{{\displaystyle \int }}}\left[-{\Delta }_y^2\widehat{X}\left(M_y\widehat{u}\right)+\left(I\left(n^T{\nabla }_y\right){\Delta }_y^2\widehat{x}\right)\widehat{u}+2\mu {\widehat{U}}_s^s\left(M_y\widehat{u}\right)+O\left({\mathrm{r}}^0\right)\widehat{u}\right]d{s}_y$$

(61)

The second term of the integral (61) can be expressed as:

$$\left(n^T{\nabla }_y\right){\Delta }_y^2\widehat{x}\left(\mathrm{r}\right)=\frac{n^T{\nabla }_{y^r}}{4\pi r^2}+O\left({\mathrm{r}}^0\right)$$

(62)

where

$$C^s\left(x\right)=I\left(x\right) c \left(x\right) \mathrm{with} c\left(x\right)=\frac{\varphi \left(x\right)}{4\pi }$$

(63)

Based on [37], the following limit may be rewritten as:

$$\underset{\mathsf{\Omega }\ni \tilde{x}\to x\in \Gamma }{\lim }{\left(\widehat{K}\widehat{u}\right)}_{\mathsf{\Omega }}^s\left(\tilde{x}\right)=-I\left(x\right)\left[-1+c\left(x\right)\right]\widehat{u}\left(x\right)+{\left(\widehat{K}\widehat{u}\right)}^s\left(x\right)$$

(64)

By augmenting ${\widehat{U}}_s^s$ to ${\widehat{U}}^s$ and employing (56) we can write (61) as:

$${\left(\widehat{K}\widehat{u}\right)}_{\mathsf{\Omega }}^s\left(\tilde{x}\right)=\underset{\Gamma }{\overset{.}{{\displaystyle \int }}}-{\Delta }_y^2\widehat{x}\left(M_y\widehat{u}\right)+\left(I\left(n^T{\nabla }_y\right){\Delta }_y^2\widehat{x}\right)\widehat{u}+2\mu {\widehat{U}}^s\left(M_y\widehat{u}\right)+O\left({\mathrm{r}}^0\right)\widehat{u}d{s}_y$$

(65)

By dividing the time interval $\left[0, T\right]$, we obtain the following integral:

$$\left(f\ast g\right)\left(\tau \right)=\underset{0}{\overset{t}{{\displaystyle \int }}}f\left(\tau -t\right)g\left(t\right)dt \mathrm{for} \tau \in \left[0,T\right]$$

(66)

in which

$$\left(f\ast g\right)\left({\tau }_n\right)\approx {\displaystyle \sum }_{k=0}^n{\omega }_{n-k}^{\Delta \tau }\left(\widehat{f}\right)g\left({\tau }_k\right)$$

(67)

According to the Lubich formula [41,42], the integration weights ${\omega }_n$ can be determined as:

$${\omega }_n^{\Delta \tau }\left(\widehat{f}\right)\Delta \frac{1}{2\pi i}\underset{\left|z\right|=R}{\overset{}{{\displaystyle \int }}}\widehat{f}\left(\frac{\gamma \left(z\right)}{\Delta \tau }\right)z^{-\left(n+1\right)}dz$$

(68)

Using $z=R{e}^{-i\phi }$, the integral in Equation (68) may be approximated as:

$${\omega }_n^{\Delta \tau }\left(\widehat{f}\right)\approx \frac{R^{-1}}{L+1}{\displaystyle \sum _{\ell =o}^L}\widehat{f}\left(s_{\ell }\right)\zeta { }^{\ell n} \mathrm{with} \zeta =e^{\frac{2\pi i}{L+1}} and s_{\ell }=\frac{\gamma \left(R{\zeta }^{-\ell }\right)}{\Delta \tau }$$

(69)

By plugging Equation (69) into Equation (67), we obtain:

$$\left(f\ast g\right)\left({\tau }_n\right)\approx {\displaystyle \sum }_{k=0}^N\frac{R^{-\left(n-k\right)}}{N+1}{\displaystyle \sum }_{\ell =0}^N\widehat{f}\left(s_{\ell }\right){\zeta }^{\ell \left(n-k\right)}g\left({\tau }_k\right) \approx \frac{R^{-n}}{N+1}{\displaystyle \sum }_{\ell =0}^N\widehat{f}\left(s_{\ell }\right)\widehat{g}\left(s_{\ell }\right){\zeta }^{\ell n}$$

(70)

with

$$\widehat{g}\left(s_{\ell }\right)={\displaystyle \sum }_{k=0}^N{R}^{k}g\left({\tau }_k\right){\zeta }^{-\ell k}.$$

(71)

Based on [39], we obtain:

$$C\left(x\right)u^g\left(x,\tau \right)=\left(v\ast t^g\right)\left(x,\tau \right)-\left(k\ast u^g\right)\left(x,\tau \right)$$

(72)

which can be expressed in Laplace domain as follows:

$$C\left(x\right){\widehat{u}}^g\left(x, s_{\ell }\right)=\left(\widehat{v}{\widehat{t}}^g\right)\left(x, s_{\ell }\right)-\left(\widehat{k}{\widehat{u}}^g\right)\left(x, s_{\ell }\right),$$

(73)

The discretization of the boundary $\Gamma =\partial \mathsf{\Omega }$ into $N_e$ boundary elements ${\overline{\tau }}_e$ leads to:

$$\Gamma \approx {\Gamma }_h=\underset{e=1}{\overset{N_e}{{\displaystyle \cup }}}{\overline{\tau }}_e$$

(74)

Now, we use $\mathrm{𝕚}$ continuous functions ${\phi }_i^{\alpha }\left[k\right]$ and $\mathrm{𝕛}$ discontinuous functions ${\psi }_j^{\beta }\left[k\right]$ to define the following subspaces:

$$S_h\left[k\right]\left({\Gamma }_{N, h}\right)≔span{\left\{{\phi }_i^{\alpha }\left[k\right]\right\}}_{i=1}^{\mathrm{𝕚}}, \alpha \ge 1$$

(75)

$$S_h\left[k\right]\left({\Gamma }_{D, h}\right)≔span{\left\{{\psi }_j^{\beta }\left[k\right]\right\}}_{j=1}^{\mathrm{𝕛}}, \beta \ge 0$$

(76)

By using (75) and (76), the unknown datum can be approximated as follows:

$${\widehat{u}}^g\left[k\right]\left(x\right)\approx {\widehat{u}}_h^g\left[k\right]\left(x\right)={\displaystyle \sum }_{i=1}^I{\widehat{u}}_{h,i}^g\left[k\right]{\phi }_i^{\alpha }\left[k\right]\left(x\right)\in S_h\left[k\right]\left({\Gamma }_{N,h}\right),$$

(77)

$${\widehat{t}}^g\left[k\right]\left(x\right)\approx {\widehat{t}}_h^g\left[k\right]\left(x\right)={\displaystyle \sum }_{j=1}^J{\widehat{t}}_{h,j}^g\left[k\right]{\psi }_j^{\beta }\left[k\right]\left(x\right)\in S_h\left[k\right]\left({\Gamma }_{D,h}\right),$$

(78)

Thus, we obtain:

$${\left[\begin{array}{c}{\widehat{V}}_{DD} \\ {\widehat{V}}_{ND} \end{array}\begin{array}{c}-{\widehat{K}}_{DN} \\ -\left(C+{\widehat{K}}_{NN}\right)\end{array}\right]}_{\ell }{\left[\begin{array}{c}{\widehat{t}}_{D,h}^g\\ {\widehat{u}}_{N,h}^g\end{array}\right]}_{\ell }={\left[\begin{array}{c}-{\widehat{V}}_{DN}\\ -{\widehat{V}}_{NN}\end{array} \begin{array}{c}\left(C+{\widehat{K}}_{DD}\right)\\ {\widehat{K}}_{ND} \end{array}\right]}_{\ell }{\left[\begin{array}{c}{\widehat{g}}_{N,h}^g\\ {\widehat{g}}_{D,h}^g\end{array}\right]}_{\ell } \ell =0\dots N$$

(79)

where

$${\widehat{S}}_{NN}≔{\widehat{V}}_{ND}{\widehat{V}}_{DD}^{-1}{\widehat{K}}_{DN}-\left(C+{\widehat{K}}_{NN}\right)$$

(80)

5. Numerical Results and Discussion

In the context of analyzing the BEM model results for solving ultrasonic wave propagation problems in three-temperature anisotropic viscoelastic porous media, as shown in Figure 2, the BEM discretization was carried out with 42 boundary elements and 68 internal points.

To demonstrate the numerical results obtained by the proposed technique, the following physical parameters were used:

The elasticity tensor:

$$\begin{gathered} \begin{array}{c}\hfill C_{ablg}=\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& 0& 0& 0\\ C_{12}& C_{11}& C_{13}& 0& 0& 0\\ C_{13}& C_{13}& C_{33}& 0& 0& 0\\ 0& 0& 0& C_{44}& 0& 0\\ 0& 0& 0& 0& C_{44}& 0\\ 0& 0& 0& 0& 0& C_{66}\end{array}\right]\end{array} \\ {C}_{11}=\frac{E^2{v}_0^2-E{E}_0}{\left(1+v\right)\left(2E{v}_0^2+E_0\left(v-1\right)\right)}, C_{12}=-\frac{E^2{v}_0^2+E{E}_{0}v}{\left(1+v\right)\left(2E{v}_0^2+E_0\left(v-1\right)\right)} \\ {C}_{13}=-\frac{E{E}_{0}v}{2E{v}_0^2+E_0\left(v-1\right)}, C_{33}=-\frac{E_0^2\left(v-1\right)}{2E{v}_0^2+E_0\left(v-1\right)} \\ {C}_{44}={\mu }_0, C_{66}=\frac{1}{2}\left(C_{11}-C_{12}\right) \end{gathered}$$

(81)

For the anisotropic viscoelastic porous media, we considered the following physical parameters [43,44]

$$v=0.95, v_0=0.49, {\mu }_0=20.98 \mathrm{GPa}, E=22 \mathrm{kPa}, E_0=447 \mathrm{kPa}$$

(82)

and therefore

$$k_1=1243 \mathrm{kPa}, k_2=442 \mathrm{kPa}$$

(83)

and

$${\rho }_s=1600 \mathrm{kg}/{\mathrm{m}}^3, {\rho }_f=1113 \mathrm{kg}/{\mathrm{m}}^3, p=25 \mathrm{MPa}, \varphi =0.15 \mathrm{and} \frac{Q}{R}=0.65$$

(84)

Figure 3, Figure 4 and Figure 5 show the distributions of the nonlinear thermal stress ${\sigma }_{11}, {\sigma }_{12}, \mathrm{and} {\sigma }_{22}$ waves along the $x_1$-axis for the electron, ion, and phonon and the total 3T with and without viscosity effect. Figure 3 shows the distribution of the nonlinear thermal stress ${\mathsf{\sigma }}_{11}$ waves for the electron ($\theta ={\theta }_e$), ion ($\theta ={\theta }_i$), phonon ($\theta ={\theta }_r$), and total 3T ($\theta ={\theta }_e+{\theta }_i+{\theta }_r$) with and without the effect of viscosity. Figure 4 shows the distribution of the nonlinear thermal stress ${\mathsf{\sigma }}_{12}$ waves for the electron ($\theta ={\theta }_e$), ion ($\theta ={\theta }_i$), phonon ($\theta ={\theta }_r$), and total ($\theta ={\theta }_e+{\theta }_i+{\theta }_r$) with and without the effect of viscosity. Figure 5 shows the distribution of the nonlinear thermal stress ${\mathsf{\sigma }}_{22}$ waves for the electron ($\theta ={\theta }_e$), ion ($\theta ={\theta }_i$), phonon ($\theta ={\theta }_r$), and the total ($\theta ={\theta }_e+{\theta }_i+{\theta }_r$) with and without the effect of viscosity.

Figure 6, Figure 7 and Figure 8 show the distributions of the nonlinear thermal stress ${\sigma }_{11}, {\sigma }_{12}, \mathrm{and} {\sigma }_{22}$ waves along the $x_1$-axis for the electron, ion, phonon, and the total 3T for isotropic and anisotropic viscoelastic porous structures. Figure 6 shows the distribution of the nonlinear thermal stress ${\mathsf{\sigma }}_{11}$ waves for the electron ($\theta ={\theta }_e$), ion ($\theta ={\theta }_i$), phonon ($\theta ={\theta }_r$), and total ($\theta ={\theta }_e+{\theta }_i+{\theta }_r$) for isotropic and anisotropic viscoelastic porous structures. Figure 7 shows the distribution of the nonlinear thermal stress ${\mathsf{\sigma }}_{12}$ waves for the electron ($\theta ={\theta }_e$), ion ($\theta ={\theta }_i$), phonon ($\theta ={\theta }_r$), and total ($\theta ={\theta }_e+{\theta }_i+{\theta }_r$) for isotropic and anisotropic viscoelastic porous structures. Figure 8 shows the distribution of the nonlinear thermal stress ${\mathsf{\sigma }}_{22}$ waves for the electron ($\theta ={\theta }_e$), ion ($\theta ={\theta }_i$), phonon ($\theta ={\theta }_r$), and total ($\theta ={\theta }_e+{\theta }_i+{\theta }_r$) for isotropic and anisotropic viscoelastic porous structures.

The validity of the outcomes of the suggested technique was not supported by any published works. On the other hand, some works in the literature can be seen as special cases of the considered general work.

Figure 9, Figure 10 and Figure 11 show the distributions of the nonlinear thermal stress ${\sigma }_{11}, {\sigma }_{12}, \mathrm{and} {\sigma }_{22}$ waves along the $x_1$-axis for the electron, ion, phonon, and the total 3T using the finite difference method (FDM) [44] and the current BEM. Figure 9 shows the distribution of the nonlinear thermal stress ${\mathsf{\sigma }}_{11}$ waves for the electron ($\theta ={\theta }_e$), ion ($\theta ={\theta }_i$), phonon ($\theta ={\theta }_r$), and the total ($\theta ={\theta }_e+{\theta }_i+{\theta }_r$) for the FDM and the BEM. Figure 10 shows the distribution of the nonlinear thermal stress ${\mathsf{\sigma }}_{12}$ waves for the electron ($\theta ={\theta }_e$), ion ($\theta ={\theta }_i$), phonon ($\theta ={\theta }_r$), and the total ($\theta ={\theta }_e+{\theta }_i+{\theta }_r$) for the FDM and BEM. Figure 11 shows the distribution of the nonlinear thermal stress ${\mathsf{\sigma }}_{22}$ waves for the electron ($\theta ={\theta }_e$), ion ($\theta ={\theta }_i$), phonon ($\theta ={\theta }_r$), and the total ($\theta ={\theta }_e+{\theta }_i+{\theta }_r$) for FDM and BEM. These figures clearly show that the BEM and FDM are in excellent agreement, supporting the validity and precision of our proposed BEM approach.

Table 1. A comparison of the required computer resources for the modeling of ultrasonic wave propagation problems in three-temperature anisotropic viscoelastic porous media.

	BEM	FDM
Number of nodes	66	40,000
Number of elements	36	16,000
CPU time (min)	2	160
Memory (MByte)	1	140
Disk space (MByte)	0	200
Accuracy of results (%)	1	2.0

shows a comparison of required computer resources for the current BEM results and the FDM results of Hu et al. [45] for the modeling of ultrasonic wave propagation problems in three-temperature anisotropic viscoelastic porous media.

6. Conclusions

The main goal of this article is to develop a novel boundary element model for describing ultrasonic thermomechanical interactions in three-temperature anisotropic viscoelastic porous media. Analytical or numerical solutions are always difficult due to the strong nonlinearity of ultrasonic wave propagation problems in three-temperature porous media, necessitating the development of new computational techniques. As a result, we employ a new BEM model to address such problems. The considered BEM model has low RAM and CPU usage due to its advantages, such as dealing with more complex shapes of porous media and not requiring the discretization of the internal domain. As a result, the considered BEM is a powerful and adaptable tool for modeling ultrasonic wave propagation in three-temperature anisotropic viscoelastic porous media. To obtain fundamental solutions, a double integral must be calculated, but this increases the total BEM computation time. To solve the current problem and improve the formulation efficiency, we propose a BEM technique. The numerical results are graphed to show the effects of viscosity and anisotropy on nonlinear ultrasonic stress waves in porous media at three temperatures. The proposed methodology’s validity, accuracy, and efficiency were demonstrated by comparing the obtained results to the corresponding solution obtained using the finite difference method (FDM). The findings of this paper contribute to the development of mathematical models that can be applied in biology, bioengineering and medicine.