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

1. Introduction

In recent years, many engineering studies have drawn researchers' attention to investigate the mechanical behaviour of viscoelastic porous 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 magnetic field [7], initial stress [8], rotation and gravity [9] on the generalized thermo-viscoelastic diffusion medium. Also, the effects of initial stress and temperature-dependent on the thermo-microstretch elastic solid have been investigated with dual-phase-lag model [10]. Xu et al. [11] studied the effects of viscoelastic dampers with high energy dissipation based on an acrylate rubber matrix. Because 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 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 get 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 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, where the bioheat equation is solved first using the GBEM based on LRBFCM to obtain the temperature distribution, and then the mechanical equation was solved using the CQBEM to obtain the displacements and stresses. The resulting linear systems have been solved by 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 $\mathsf{\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+\mathbf{F}=\rho \ddot{\mathbf{u}}+\varphi {\rho }_f\left({\ddot{\mathbf{u}}}_f-\ddot{\mathbf{u}}\right)$$

(1)

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

(2)

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

(3)

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

(4)

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

(5)

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

(6)

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

$$\left.\begin{array}{ll}{\widehat{B}}_{\tilde{x}}{\widehat{\mathbf{u}}}^g\left(\tilde{\mathbf{x}}\right)=0& \mathrm{f}\mathrm{o}\mathrm{r}\tilde{x}\in \mathsf{\Omega }\\ {\widehat{\mathbf{u}}}^g\left(\mathbf{x}\right)={\widehat{g}}_D& \mathrm{f}\mathrm{o}\mathrm{r}x\in {\Gamma }_D\\ {\widehat{\mathbf{t}}}^g\left(\mathbf{x}\right)={\widehat{g}}_N& \mathrm{f}\mathrm{o}\mathrm{r}x\in {\Gamma }_N\end{array}\right\}$$

(7)

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

$${\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}{∆}_{\tilde{\mathrm{x}}}+\frac{{s\varphi }^2}{R}& 0\end{array}\right],$$

$${\widehat{\mathbf{t}}}^g\left(x\right)=\left[\begin{array}{ccc}{T}_{\mathrm{x}}^e& -\alpha {\mathbf{n}}_{\mathbf{x}}& 0\\ s\beta {\mathbf{n}}_{\mathbf{x}}^T& \frac{\beta }{s{\rho }_f}{\mathbf{n}}_{\mathbf{x}}^T{\nabla }_{\mathrm{x}}& 0\end{array}\right]\left[\begin{array}{c}\widehat{\mathbf{u}}\left(\mathbf{x}\right)\\ \widehat{p}\left(\mathbf{x}\right)\\ \theta \left(\mathbf{x}\right)\end{array}\right],\beta =\frac{{\varphi }^{2}sK{\rho }_f}{{\varphi }^2+sK({\rho }_e+\varphi {\rho }_f)}$$

(8)

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

The three-temperature radiative diffusion equations are as follows

$$C_{ve}\frac{\partial {\theta }_e(\mathrm{r},\tau )}{\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(\mathrm{r},\tau )}{\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(\mathrm{r},\tau )}{\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}{ll}{C}_{v\alpha }=C_e& \alpha =e\\ C_{v\alpha }=C_i& \alpha =i\\ {C_{v\alpha }=C}_r{T}_r^3& \alpha =r\end{array}\right.\mathrm{a}\mathrm{n}\mathrm{d}{\mathbb{K}}_{\alpha }=\left\{\begin{array}{ll}{\mathbb{A}}_e{T}_e^{5/2}& \alpha =e\\ {\mathbb{A}}_i{T}_i^{5/2}& \alpha =i\\ {\mathbb{A}}_r{T}_r^{3+\mathbb{B}}& \alpha =r\end{array}\right.$ where $e,i\mathrm{a}\mathrm{n}\mathrm{d}r$ denote electron, ion, and phonon, respectively.

3. Boundary Element Implementation for the Temperature Field

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

$$\nabla \left[{\mathbb{K}}_{\alpha }\nabla {\theta }_{\alpha }(\mathrm{r},\tau )\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}{{\mathbf{x}}_0}{e}^{\left(-\frac{{\mathbf{r}}_a}{{\mathbf{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}}(r,\tau )=\left\{\begin{array}{ll}-\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|}_{{\mathsf{\Gamma }}_1}=0,\alpha =e,i,{\left.{\theta }_p\right|}_{{\mathsf{\Gamma }}_1}=g_2\left(x,\tau \right)$$

(13)

$${\left.{\mathbb{K}}_{\alpha }\frac{\partial {\theta }_{\alpha }}{\partial n}\right|}_{{\mathsf{\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

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

(16)

which can be expressed as follows

$$C{\theta }_{\alpha }={\int }_S^{}\left[{\theta }_{\alpha }{q}^{*}-{\theta }_{\alpha }^{*}q\right]dS-{\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 \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 ={\int }_S^{}\left[{\theta }_{\alpha }{q}^{*}-{\theta }_{\alpha }^{*}q\right]dS+\sum _{j=1}^N{a}^j\left(\tau \right)D^{-1}\left(C{\widehat{\theta }}_{\alpha }^j-{\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)=\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 get

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

(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}{∆{\tau }^m}\le 1$, $∆{\tau }^m={\tau }^{m+1}-{\tau }^m$.

By differentiating (28), we get

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

(30)

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

$$\left(\frac{C}{∆{\tau }^m}+\mathsf{\theta }H\right){\theta }_{\alpha }^{m+1}-\mathsf{\theta }G{Q}^{m+1}=\left(\frac{C}{∆{\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

$$\mathbb{a}\mathsf{{\rm X}}=\mathbb{b}$$

(32)

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

4. Boundary Element Implementation for the Poroelastic Fields

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

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

(33)

where the integral operators are

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

(34)

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

(35)

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

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

(36)

The fundamental solution can be expressed as [35]

$${\widehat{\mathbf{U}}}^s\left(\mathrm{r}\right)=\frac{1}{4\pi \mathrm{r}\left(\rho -\beta {\rho }^f\right)}\left[{\mathbb{R}}_1\frac{\left(k_4^2-k_2^2\right)}{\left(k_1^2-k_2^2\right)}{e}^{-k_1\mathrm{r}}-{\mathbb{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-{\mathbb{R}}_3\right)e^{-k_3\mathrm{r}}\right]$$

(37)

where

$${\mathbb{R}}_j=\frac{3{\nabla }_{\mathbf{y}}\mathrm{r}{\nabla }_{\mathbf{y}}^T\mathrm{r}-I}{{\mathrm{r}}^2}+k_j\frac{3{\nabla }_{\mathbf{y}}\mathrm{r}{\nabla }_{\mathbf{y}}^T\mathrm{r}-I}{\mathrm{r}}+k_j^2{\nabla }_{\mathbf{y}}\mathrm{r}{\nabla }_{\mathbf{y}}^T\mathrm{r}$$

(38)

Equation (37) can be expressed as

$${\widehat{\mathbf{U}}}^s\left(\mathrm{r}\right)=\frac{1}{4\pi \mu \mathrm{r}\left(\lambda +2\mu \right)}\left[\left(\lambda +\mu \right){\nabla }_{\mathbf{y}}r{\nabla }_{\mathbf{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{\mathbf{U}}}^s\left(\mathrm{r}\right)={\widehat{\mathbf{U}}}_s^s\left(\mathrm{r}\right)+{\widehat{\mathbf{U}}}_r^s\left(\mathrm{r}\right) \\ =\frac{1}{\mu }\left[I{∆}_{\mathbf{y}}-\frac{\lambda +\mu }{\lambda +2\mu }{\nabla }_{\mathbf{y}}{\nabla }_{\mathbf{y}}^T\right]{∆}_{\mathbf{y}}\widehat{\mathbf{x}}\left(\mathrm{r}\right) \\ -\frac{1}{\mu }\left[\left(\left(k_1^2+k_2^2\right){∆}_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 }_{\mathbf{y}}{\nabla }_{\mathbf{y}}^T\right]\widehat{\mathbf{x}}\left(\mathrm{r}\right) \end{gathered}$$

(40)

in which

$$\begin{gathered} \widehat{\mathbf{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{\mathbf{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{\mathbf{P}}}^s\left(\mathrm{r}\right)=\frac{{\widehat{\mathbf{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{\mathbf{x}}\in \Omega \to \mathbf{x}\in \Gamma$ to (34) to obtain

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

(45)

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

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

(46)

in which

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

(47)

and

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

(48)

By using equations (45) - (48), the boundary integral equation in Laplace domain can be expressed as

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

(49)

The poroelastodynamic boundary integral equation can be expressed using the inverse Laplace transformation as

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

(50)

The fundamental solution is as follows [37]

$${\left({\widehat{T}}_{\mathbf{y}}\widehat{\mathbf{U}}\right)}^T={\left[\left[\begin{array}{cc}{\widehat{T}}_{\mathbf{y}}^e& s\alpha {\mathbf{n}}_{\mathbf{y}}\\ -\beta {\mathbf{n}}_{\mathbf{y}}^T& \frac{\beta }{s{\rho }_0^f}{\mathbf{n}}_{\mathbf{y}}^T{\nabla }_{\mathbf{y}}\end{array}\right]\left[\begin{array}{cc}{\widehat{\mathbf{U}}}^s& {\widehat{\mathbf{U}}}^f\\ {\left({\widehat{\mathbf{P}}}^s\right)}^T& {\widehat{P}}^f\end{array}\right]\right]}^T={\left[\begin{array}{cc}{\widehat{\mathbf{T}}}^s& {\widehat{\mathbf{T}}}^f\\ {\left({\widehat{\mathbf{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{.}{\int }}\left({\nabla }_{\mathbf{y}}\times \mathrm{a},n_{\mathbf{y}}\right)d{s}_{\mathbf{y}}=-\underset{\partial \Gamma }{\overset{,}{\int }}\left(\mathrm{a},v\right)d{\gamma }_{\mathbf{y}}=-\underset{\varphi }{\overset{,}{\int }}\left(\mathrm{a},v\right)d{\gamma }_{\mathbf{y}}=0$$

(52)

where

$$\underset{\Gamma }{\overset{.}{\int }}\left(n_{\mathbf{y}}\times {\nabla }_{\mathbf{y}},\mathrm{a}\right)d{s}_{\mathbf{y}}=0$$

(53)

We can use (53) to obtain the following formula

$$\underset{\Gamma }{\overset{.}{\int }}\left(M_{\mathbf{y}}\mathrm{a}\right)d{s}_{\mathbf{y}}=0,M_{\mathbf{y}}={\left({\nabla }_{\mathbf{y}}{\nabla }_{\mathbf{y}}^T\right)}^T-{\nabla }_{\mathbf{y}}{\nabla }_{\mathbf{y}}^T,\mathrm{a}=vu$$

(54)

According to [40], we obtain

$$\underset{\Gamma }{\overset{.}{\int }}\left(M_{\mathbf{y}}v\right)\mathbf{u}d{s}_{\mathbf{y}}=-\underset{\Gamma }{\overset{,}{\int }}v\left(M_{\mathbf{y}}\mathbf{u}\right)d{s}_{\mathbf{y}}$$

(55)

$$\underset{\Gamma }{\overset{.}{\int }}{\left(M_{\mathbf{y}}v\right)}^T\mathbf{u}d{s}_{\mathbf{y}}=\underset{\Gamma }{\overset{.}{\int }}{v}^T\left(M_{\mathbf{y}}\mathbf{u}\right)d{s}_{\mathbf{y}}$$

(56)

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

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

(57)

According to [37], we get

$${\left({\widehat{\mathbf{T}}}^s\right)}^T=\left(\lambda +2\mu \right)n_{\mathbf{y}}{\nabla }_{\mathbf{y}}^T{\widehat{\mathbf{U}}}_{sing}^s-\mu \left({\mathit{n}}_{\mathbf{y}}\times \left({\nabla }_{\mathbf{y}}\times {\widehat{\mathbf{U}}}_{sing}^s\right)\right)+2\mu M_{\mathbf{y}}{\widehat{\mathbf{U}}}_{sing}^s+o\left({\mathrm{r}}^0\right)$$

(58)

which can be expressed using (40) as

$${\left({\widehat{\mathbf{T}}}^s\right)}^T=M_{\mathbf{y}}{∆}_{\mathbf{y}}^2\widehat{X}+\mathbf{I}\left({\mathbf{n}}^T{\nabla }_{\mathbf{y}}\right){∆}_{\mathbf{y}}^2\widehat{X}+2\mu {\left(M_{\mathbf{y}}{\widehat{\mathbf{U}}}_{sing}^s\right)}^T+o\left({\mathrm{r}}^0\right)$$

(59)

Using (35) and (59), we obtain

$${\left(\widehat{k}\widehat{\mathbf{u}}\right)}_{\Omega }^s\left(\tilde{\mathbf{x}}\right)={\int }_{┌}^{.}\left[\left(M_{\mathbf{y}}{∆}_{\mathbf{y}}^2\widehat{X}\right)\widehat{\mathbf{u}}+\left(\mathbf{I}\left({\mathbf{n}}^T{\nabla }_{\mathbf{y}}\right){∆}_{\mathbf{y}}^2\widehat{X}\right)\widehat{\mathbf{u}}+2\mu {\left(M_{\mathbf{y}}{\widehat{\mathbf{U}}}_{sing}^s\right)}^T\widehat{\mathbf{u}}+O\left({\mathrm{r}}^0\right)\widehat{\mathbf{u}}\right]{ds}_{\mathbf{y}}$$

(60)

Based on [39], we get

$${\left(\widehat{K}\widehat{\mathbf{u}}\right)}_{\Omega }^s\left(\tilde{\mathbf{x}}\right)=\underset{┌}{\overset{.}{\int }}\left[-{∆}_{\mathbf{y}}^2\widehat{X}\left(M_{\mathbf{y}}\widehat{\mathbf{u}}\right)+\left(\mathbf{I}\left({\mathbf{n}}^T{\nabla }_{\mathbf{y}}\right){∆}_{\mathbf{y}}^2\widehat{x}\right)\widehat{\mathbf{u}}+2\mu {\widehat{U}}_s^s\left(M_{\mathbf{y}}\widehat{\mathbf{u}}\right)+O\left({\mathrm{r}}^0\right)\widehat{\mathbf{u}}\right]{ds}_{\mathbf{y}}$$

(61)

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

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

(62)

where

$$C^s\left(\mathbf{x}\right)=I\left(\mathbf{x}\right)c\left(\mathbf{x}\right)\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}c\left(\mathbf{x}\right)=\frac{\varphi \left(\mathbf{x}\right)}{4\pi }$$

(63)

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

$$\underset{\Omega \ni \tilde{x}\to x\in \Gamma }{\lim }{\left(\widehat{K}\widehat{\mathbf{u}}\right)}_{\Omega }^s\left(\tilde{\mathbf{x}}\right)=-I\left(\mathbf{x}\right)\left[-1+c\left(\mathbf{x}\right)\right]\widehat{\mathbf{u}}\left(\mathbf{x}\right)+{\left(\widehat{K}\widehat{\mathbf{u}}\right)}^s\left(\mathbf{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{\mathbf{u}}\right)}_{\Omega }^s\left(\tilde{\mathbf{x}}\right)=\underset{\Gamma }{\overset{.}{\int }}-{∆}_{\mathbf{y}}^2\widehat{\mathbf{x}}\left(M_{\mathbf{y}}\widehat{\mathbf{u}}\right)+\left(\mathbf{I}\left({\mathbf{n}}^T{\nabla }_{\mathbf{y}}\right){∆}_{\mathbf{y}}^2\widehat{\mathbf{x}}\right)\widehat{\mathbf{u}}+2\mu {\widehat{U}}^s\left(M_{\mathbf{y}}\widehat{\mathbf{u}}\right)+O\left({\mathrm{r}}^0\right)\widehat{\mathbf{u}}{ds}_{\mathbf{y}}$$

(65)

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

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

(66)

In which

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

(67)

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

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

(68)

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

$${\omega }_n^{∆\tau }\left(\widehat{f}\right)\approx \frac{R^{-1}}{L+1}\sum _{\mathcal{l}=o}^L\widehat{f}\left(s_{\mathcal{l}}\right){\zeta }^{\mathcal{l}n}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\zeta =e^{\frac{2\pi i}{L+1}}and{s}_{\mathcal{l}}=\frac{\gamma \left({R\zeta }^{-\mathcal{l}}\right)}{∆\tau }$$

(69)

By plugging Equation (69) into Equation (67), we get

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

(70)

with

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

(71)

Based on [39], we get

$$C\left(\mathbf{x}\right){\mathbf{u}}^g\left(\mathbf{x},\tau \right)=\left({v*\mathbf{t}}^g\right)\left(\mathbf{x},\tau \right)-\left({k*\mathbf{u}}^g\right)\left(\mathbf{x},\tau \right)$$

(72)

that can be expressed in Laplace domain as follows

$$C\left(\mathbf{x}\right){\widehat{\mathbf{u}}}^g\left({\mathbf{x},s}_{\mathcal{l}}\right)=\left(\widehat{v}{\widehat{\mathbf{t}}}^g\right)\left({\mathbf{x},s}_{\mathcal{l}}\right)-\left({\widehat{k}\widehat{\mathbf{u}}}^g\right)\left({\mathbf{x},s}_{\mathcal{l}}\right),$$

(73)

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

$$\Gamma \approx {\Gamma }_h=\bigcup _{e=1}^{{\mathit{N}}_{\mathit{e}}}{\overline{\tau }}_e$$

(74)

Now, we use $\mathbb{i}$ continuous functions ${\phi }_i^{\alpha }\left[k\right]$ and $\mathbb{j}$ 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}^{\mathbb{i}},\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}^{\mathbb{j}},\beta \ge 0$$

(76)

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

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

(77)

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

(78)

Thus, we obtain

$${\left[\genfrac{}{}{0pt}{}{{\widehat{V}}_{DD}}{{\widehat{V}}_{ND}}\genfrac{}{}{0pt}{}{-{\widehat{K}}_{DN}}{-\left(C+{\widehat{K}}_{NN}\right)}\right]}_{\mathcal{l}}{\left[\genfrac{}{}{0pt}{}{{\widehat{\mathbf{t}}}_{D,h}^g}{{\widehat{\mathbf{u}}}_{N,h}^g}\right]}_{\mathcal{l}}={\left[\genfrac{}{}{0pt}{}{-{\widehat{V}}_{DN}}{{-\widehat{V}}_{NN}}\genfrac{}{}{0pt}{}{\left(C+{\widehat{K}}_{DD}\right)}{{\widehat{K}}_{ND}}\right]}_{\mathcal{l}}{\left[\genfrac{}{}{0pt}{}{{\widehat{g}}_{N,h}^g}{{\widehat{g}}_{D,h}^g}\right]}_{\mathcal{l}}\mathcal{l}=0...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 of 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

$$C_{ablg}=\left[\begin{array}{ccc}\begin{array}{cc}\begin{array}{c}\begin{array}{c}{C}_{11}\\ C_{12}\end{array}\\ \begin{array}{c}{C}_{13}\\ 0\end{array}\\ \begin{array}{c}0\\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}{C}_{12}\\ C_{11}\end{array}\\ \begin{array}{c}{C}_{13}\\ 0\end{array}\\ \begin{array}{c}0\\ 0\end{array}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}{C}_{13}\\ C_{13}\end{array}\\ \begin{array}{c}{C}_{33}\\ 0\end{array}\\ \begin{array}{c}0\\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}0\\ C_{44}\end{array}\\ \begin{array}{c}0\\ 0\end{array}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}{C}_{44}\\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}0\\ C_{66}\end{array}\end{array}\end{array}\end{array}\right]$$

(81)

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

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

$$v=0.95,v_0=0.49,{\mu }_0=20.98\mathrm{G}\mathrm{P}\mathrm{a},E=22\mathrm{k}\mathrm{P}\mathrm{a},E_0=447\mathrm{k}\mathrm{P}\mathrm{a}$$

and therefore

$$k_1=1243\mathrm{k}\mathrm{P}\mathrm{a},k_2=442\mathrm{k}\mathrm{P}\mathrm{a}$$

and

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

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

Figure 6, Figure 7 and Figure 8 show the distributions of the nonlinear thermal stress ${\sigma }_{11},{\sigma }_{12},\mathrm{a}\mathrm{n}\mathrm{d}{\sigma }_{22}$ waves along $x_1$-axis for electron, ion, phonon and total 3T for isotropic and anisotropic viscoelastic porous structures. Figure 6 shows the distribution of the nonlinear thermal stress ${\mathsf{\sigma }}_{11}$ waves for 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 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 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 literary works 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{a}\mathrm{n}\mathrm{d}{\sigma }_{22}$ waves along $x_1$-axis for electron, ion, phonon and total 3T using the finite difference method (FDM) [44] and current BEM. Figure 9 shows the distribution of the nonlinear thermal stress ${\mathsf{\sigma }}_{11}$ waves for electron ($\theta ={\theta }_e$), ion ($\theta ={\theta }_i$), phonon ($\theta ={\theta }_r$), and total ($\theta ={\theta }_e+{\theta }_i+{\theta }_r$) for FDM and BEM. Figure 10 shows the distribution of the nonlinear thermal stress ${\mathsf{\sigma }}_{12}$ waves for electron ($\theta ={\theta }_e$), ion ($\theta ={\theta }_i$), phonon ($\theta ={\theta }_r$), and total ($\theta ={\theta }_e+{\theta }_i+{\theta }_r$) for FDM and BEM. Figure 11 shows the distribution of the nonlinear thermal stress ${\mathsf{\sigma }}_{22}$ waves for electron ($\theta ={\theta }_e$), ion ($\theta ={\theta }_i$), phonon ($\theta ={\theta }_r$), and 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 shows a comparison of required computer resources for the current BEM results and FDM results of Hu et al. [44] of modeling of ultrasonic wave propagation problems in three-temperature anisotropic viscoelastic porous media.

6. Conclusion

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 in three-temperature porous media problems, 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 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 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