Boundary Element Modeling and Simulation Algorithm for Fractional Bio-Thermomechanical Problems of Anisotropic Soft Tissues

Abstract

The main purpose of this chapter is to propose a novel boundary element modeling and simulation algorithm for solving fractional bio-thermomechanical problems in anisotropic soft tissues. The governing equations are studied on the basis of the thermal wave model of bio-heat transfer (TWMBT) and Biot’s theory. These governing equations are solved using the boundary element method (BEM), which is a flexible and effective approach since it deals with more complex shapes of soft tissues and does not need the internal domain to be discretized, also, it has low RAM and CPU usage. The transpose-free quasi-minimal residual (TFQMR) solver are implemented with a dual-threshold incomplete LU factorization technique (ILUT) preconditioner to solve the linear systems arising from BEM. Numerical findings are depicted graphically to illustrate the influence of fractional order parameter on the problem variables and confirm the validity, efficiency and accuracy of the proposed BEM technique.

1. Introduction

Human body is a complex thermal system, Arsene d’Arsonval and Claude Bernard have shown that the temperature difference between arterial blood and venous blood is due to oxygenation of blood [1]. A large number of research papers in bio-heat transfer over the past few decades have focused on an understanding of the impact of blood flow on the temperature distribution within living tissues. Pennes [2] was the first attempt to explain the temperature distribution in human tissue with blood flow effect. The improvement of numerical models for determination of temperature distribution in living tissues has been a topic of interest for numerous researchers. Askarizadeh and Ahmadikia [3] introduced analytical solutions for the transient Fourier and non-Fourier bio-heat transfer equations. Li et al. [4] studied the bio-thermomechanical interactions within the human skin in the context of generalized thermoelasticity.

Analytical solutions for the current problem [5, 6] are very difficult to obtain, so numerical methods have become the main way for solving these problems [7, 8, 9, 10]. The boundary element method (BEM) [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21] is one of the numerical methods used to solve the current general problem [22, 23, 24, 25, 26, 27, 28, 29, 30, 31]. Generally, Laplace-domain fundamental solutions are easier to obtain than time-domain fundamental solutions for engineering and scientific problems [32, 33]. consequently, the BEM is very helpful for problems that did not have time-domain fundamental solutions, because it requires the Laplace-domain fundamental solutions of the problem’s governing equations. So, BEM expands the range of engineering problems that can be solved with the classical time-domain BEM.

The main aim of this chapter is to propose a new boundary element fractional model for describing the bio-thermomechanical properties of anisotropic soft tissues. The dual reciprocity boundary element method has been used to solve the TWMBT for obtaining the temperature distribution, and then the BEM has been used to obtain the displacement and stress at each time step. The linear systems from BEM were solved by the TFQMR solver with the ILUT preconditioner which reduces the number of iterations and the total CPU time.

A brief summary of the chapter is as follows: Section 1 introduces the background and provides the readers with the necessary information to books and articles for a better understanding of bio-thermomechanical problems in anisotropic soft tissues Section 2 describes the BEM modeling of the bio-thermomechanical interactions and introduces the partial differential equations that govern its related problems. Section 3 outlines the dual reciprocity boundary element method (DRBEM) for temperature field. Section 4 discusses the convolution quadrature boundary element method (CQBEM) for poro-elastic field. Section 5 presents the new numerical results that describe the bio-thermomechanical problems in anisotropic soft tissues.

2. Formulation of the problem

Consider an anisotropic soft tissue in the Cartesian coordinate system ${\mathit{Ox}}_1{x}_2{x}_3$ as shown in Figure 1. It occupies the region $\mathrm{\Omega }=\left\{\left(x_1,x_2,x_3\right):0<x_1<\underset{¯}{\alpha },0<x_2<\underset{¯}{\beta },0<x_3<\underset{¯}{\gamma }\right\}$ with boundary $\mathrm{\Gamma }$ that is subdivided into two non-intersective parts ${\mathrm{\Gamma }}_{\mathrm{D}}$ and ${\mathrm{\Gamma }}_N$.

Figure 1.

Geometry of the current problem.

The governing equations which model the fractional bio-thermomechanical problems in anisotropic soft tissues can be written as follows [34, 35].

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

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

$$\sigma =\left(C_{\mathit{ai}\mathit{lg}}\mathrm{tr}\in -\mathit{Ap}\right)I-\mathfrak{B}\theta ,$$E3

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

$$\zeta =A\mathrm{tr}\in +\frac{{\varphi }^2}{R}P$$E5

where the fluid was modeled by the following Darcy’s law [36].

$$q=-K\left(\nabla p+{\rho }_f\ddot{u}+\frac{{\rho }_a+{\mathit{\varphi \rho }}_f}{\varphi }\left({\ddot{u}}_f-\ddot{u}\right)\right)$$E6

The fractional order equation which describes the TWMBT can be expressed as [37].

$$\begin{array}{l}\nabla \mathcal{K}\nabla \mathcal{T}\left(r,\tau \right)+W_b{C}_b\left({\mathcal{T}}_b-\mathcal{T}\right)+Q_{\mathit{met}}+Q_{\mathit{ext}}+\frac{\overline{\tau }}{\alpha !}\left(-W_b{C}_b{D}_{\tau }^{\alpha }\mathcal{T}+D_{\tau }^{\alpha }{Q}_{\mathit{met}}+D_{\tau }^{\alpha }{Q}_{\mathit{ext}}\right)\\ =\mathit{\rho C}\left[\frac{\overline{\tau }}{\alpha !}{D}_{\tau }^{\alpha +1}\mathcal{T}+\frac{\partial \mathcal{T}}{\partial \tau }\right],0<\alpha \le 1\end{array}$$E7

where $\sigma$, $\in$, $C_{\mathit{\text{ajlg}}}$, $\rho ={\rho }_s\left(1-\varphi \right)+{\mathit{\varphi \rho }}_f$, ${\rho }_s$, ${\rho }_f$, $u$, $u_f$, $\mathit{FF}$ and $\mathit{qq}$ are total stress tensor, linear strain tensor, constant elastic moduli, bulk density, solid density, fluid density, solid displacement, fluid displacement, bulk body forces and specific flux of the fluid, respectively, $\mathfrak{B}$ are stress-temperature coefficients, $\mathrm{tr}$ denotes the trace, $A=\varphi \left(1+Q/R\right)$ is Biot’s coefficient, $Q$ and $R$ are the solid–fluid coupling parameters, $p$ is the fluid pressure in the vasculature, $\zeta$ is the fluid volume variation measured in unit reference volume, $\varphi =\frac{V^f}{V}$ is the porosity, $V^f$ is the fluid volume, $V=V^f+V^s$ is the bulk volume, $V^s$ is the solid volume, $\tau$ is the time, $K$ is the permeability, ${\rho }_a=\mathbb{C}{\mathit{\varphi \rho }}_f$ where $\mathbb{C}=0.66$ at low frequency [38], $\mathit{K}$ is the soft tissue thermal conductivity, $W_b$ is the blood perfusion rate, $C_b$ is the blood specific heat, $T_b$ is the arterial blood temperature, $T$ is the soft tissue temperature, $\overline{\tau }$ is the thermal relaxation time $\rho$ is the soft tissue density, $C$ is the soft tissue specific heat, $Q_{\mathit{met}}$ is the metabolic heat generation and $Q_{\mathit{ext}}$ is the external heat generation.

According to Bonnet [39], our problem can be expressed as a matrix system as [40].

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

where

$${\widehat{B}}_{\tilde{x}}=\left[\begin{array}{ccc}{B}_{\tilde{x}}^e+s^2\left(\rho -{\mathit{\beta \rho }}_f\right)I& \left(\alpha -\beta \right){\nabla }_{\tilde{x}}& -\mathfrak{B}{\nabla }_{\tilde{x}}\\ s\left(\alpha -\beta \right){\nabla }_{\tilde{x}}^T& -\frac{\beta }{s{\rho }_f}{\mathrm{\Delta }}_{\tilde{x}}+\frac{s{\varphi }^2}{R}& 0\end{array}\right]$$E9

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

3. Boundary element implementation for bioheat transfer field

Through this chapter, we supposed that $Q_{\mathit{met}}$ and ${\mathcal{T}}_b$ are constants and $\theta \left(r,\tau \right)=\mathcal{T}\left(r,\tau \right)-\mathcal{T}\left(r,0\right)$. Thus, Eq. (7) can be written as

$$\mathit{\rho C}\frac{\overline{\tau }}{\alpha !}{D}_{\tau }^{\alpha +1}\theta +\mathit{\rho C}\frac{\partial \theta }{\partial \tau }+\frac{\overline{\tau }}{\alpha !}{W}_b{C}_b{D}_{\tau }^{\alpha }\theta +W_b{C}_b\theta =\mathcal{K}\frac{{\partial }^2\theta }{\partial x^2}+q,0<\alpha \le 1$$E11

According to finite difference scheme of Caputo [22] and using the fundamental solution of difference equation resulting from fractional bio-heat Eq. (11) [41], we can write the following dual reciprocity boundary integral equation

$$C_i{\theta }_i+{\int }_{\mathrm{\Gamma }}{q}^{\ast }\mathit{\theta d}\mathrm{\Gamma }-{\int }_{\mathrm{\Gamma }}{\theta }^{\ast }\mathit{qd}\mathrm{\Gamma }=\sum _{j=1}^{N+L}{\alpha }_j\left(C_j{\widehat{\theta }}_{\mathit{ij}}+{\int }_{\mathrm{\Gamma }}{q}^{\ast }{\widehat{\theta }}_{j}d\mathrm{\Gamma }-{\int }_{\mathrm{\Gamma }}{\theta }^{\ast }{\widehat{q}}_{j}d\mathrm{\Gamma }\right)$$E12

in which

$$C_i=\frac{\gamma }{2\pi },q=\frac{\partial \theta }{\partial n},q^{\ast }=\frac{\partial {\theta }^{\ast }}{\partial n},{\theta }^{\ast }=ln\left(\frac{1}{r}\right)$$E13

where $n$ is the outward unit normal vector to boundary $\mathrm{\Gamma }$, $r$ is the distance between source point $\mathit{i}$ and considered point $j$, $N$ is the number of boundary nodes and $L$ is the number of internal nodes.

where

$$\alpha ={\mathbb{R}}^{-1}\tilde{f}={\mathbb{R}}^{-1}\left(\tilde{a}\frac{{\partial }^2\theta }{\partial {\tau }^2}+\tilde{b}\frac{\partial \theta }{\partial \tau }+\tilde{c}\theta +\tilde{d}\right)$$E14

The discretization process for Eq. (12) leads to

$$\begin{array}{ll}& C_i{\theta }_i+\sum _{k=1}^N{\int }_{{\mathrm{\Gamma }}_k}{q}^{\ast }\mathit{\theta d}\mathrm{\Gamma }-\sum _{k=1}^N{\int }_{{\mathrm{\Gamma }}_k}{\theta }^{\ast }\mathit{qd}\mathrm{\Gamma }\hfill \\ & =\sum _{j=1}^{N+L}{\alpha }_j\left(C_i{\widehat{\theta }}_{\mathit{ij}}+\sum _{k=1}^N{\int }_{{\mathrm{\Gamma }}_k}{Z}_{\mathit{ik}}{\widehat{\theta }}_{\mathit{kj}}d\mathrm{\Gamma }-\sum _{k=1}^N{\int }_{{\mathrm{\Gamma }}_k}{G}_{\mathit{ik}}{\widehat{q}}_{\mathit{kj}}d\mathrm{\Gamma }\right)\end{array}$$E15

After interpolation and integration processes over boundary elements, Eq. (15) can be expressed as

$$C_i{\theta }_i+\sum _{k=1}^N{Z}_{\mathit{ik}}{\theta }_k-\sum _{k=1}^N{G}_{\mathit{ik}}{q}_k=\sum _{j=1}^{N+L}{\alpha }_j\left(C_i{\widehat{\theta }}_{\mathit{ij}}+\sum _{k=1}^N{Z}_{\mathit{ik}}{\widehat{\theta }}_{\mathit{kj}}-\sum _{k=1}^N{G}_{\mathit{ik}}{\widehat{q}}_{\mathit{kj}}\right)$$E16

The matrix form of Eq. (16) can be written using (14) as

$$\mathit{Z\theta }-\mathit{Gq}=\left(Z\widehat{\mathrm{\Theta }}-G\widehat{Q}\right){\mathbb{R}}^{-1}\left(\tilde{a}\frac{{\partial }^2\theta }{\partial {\tau }^2}+\tilde{b}\frac{\partial \theta }{\partial \tau }+\tilde{c}\theta +\tilde{d}\right)$$E17

which also can be written

$$X\left(\tilde{a}\frac{{\partial }^2\theta }{\partial {\tau }^2}+\tilde{b}\frac{\partial \theta }{\partial \tau }+\tilde{c}\theta +\tilde{d}\right)+\mathit{Z\theta }=\mathit{Gq}$$E18

where

$$X=\left(Z\widehat{\mathrm{\Theta }}-G\widehat{Q}\right){\mathbb{R}}^{-1}$$

The boundary and initial conditions

$$\theta \left(x,y;\tau \right)=0$$E19

$$\frac{\partial \theta \left(x,y;0\right)}{\partial \tau }=\vartheta \left(x,y;0\right)=0$$E20

$$\theta \left(x,y;0\right)=\left\{\begin{array}{ll}{1}^{\circ }\mathrm{C}& -0.02\le x,y\le 0.02\\ 0& \text{other}x,y\end{array}\right.$$E21

The time discretization has been performed as follows

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

$$\theta =\left(1-{\theta }_u\right){\theta }^m+{\theta }_u{\theta }^{m+1}$$E23

$$\frac{\partial \theta }{\partial \tau }=\frac{1}{\mathrm{\Delta }\tau }\left({\theta }^{m+1}-{\theta }^m\right)$$E24

$$\frac{{\partial }^2\theta }{\partial {\tau }^2}=\frac{1}{\mathrm{\Delta }{\tau }^2}\left({\theta }^{m+1}+{\theta }^{m-1}-2{\theta }^m\right)$$E25

Substituting from Eqs. (22)–(25) into (20), we obtain

$$\begin{array}{c}\left(\frac{X\tilde{a}}{\mathrm{\Delta }{\tau }^2}+\frac{X\tilde{b}}{\mathrm{\Delta }\tau }+X\tilde{c}{\theta }_u+{\theta }_{u}Z\right){\theta }^{m+1}-{\theta }_q{\mathit{Gq}}^{m+1}+X\tilde{d}\\ =\left(\frac{2X\tilde{a}}{\mathrm{\Delta }{\tau }^2}+\frac{X\tilde{b}}{\mathrm{\Delta }\tau }-X\tilde{c}\left(1-{\theta }_u\right)-Z\left(1-{\theta }_u\right)\right){\theta }^m-\frac{X\tilde{a}}{\mathrm{\Delta }{\tau }^2}{\theta }^{m-1}+\left(1-{\theta }_q\right){\mathit{Gq}}^m\end{array}$$E26

Thus, with the temperature $\theta$ determined, the remaining task is to solve the problem (8).

4. Boundary element implementation for the poro-elastic fields

The representation formula of (8) that describes the unknown field ${\widehat{u}}^g$ can be written as

$${\widehat{u}}^g\left(\tilde{x}\right)={\left({\widehat{V}\widehat{t}}^g\right)}_{\Omega }\left(\tilde{x}\right)-{\left({\widehat{K}\widehat{u}}^g\right)}_{\Omega }\left(\tilde{x}\right)\text{for}\tilde{x}\in \mathit{\Omega }$$E27

where

$${\left({\widehat{V}\widehat{t}}^g\right)}_{\Omega }\left(\tilde{x}\right)=\underset{\Gamma }{\overset{\cdot }{\int }}{\widehat{U}}^T\left(y-\tilde{x}\right){\widehat{t}}^g\left(y\right){\mathit{ds}}_y$$E28

$${\left({\widehat{K}\widehat{u}}^g\right)}_{\Omega }\left(\tilde{x}\right)=\underset{\Gamma }{\overset{\cdot }{\int }}{\left({\widehat{T}}_y\widehat{U}\right)}^T\left(y-\tilde{x}\right){\widehat{u}}^g\left(y\right){\mathit{ds}}_y$$E29

For anisotropic case, the Laplace domain fundamental solution $\widehat{U}\left(r\right)$ and the corresponding traction ${\widehat{T}}_v$ can be expressed as [40].

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

where the solid displacement fundamental solution ${\widehat{U}}^s\left(r\right)$ may be expressed as

$${\widehat{U}}^s\left(r\right)=\frac{1}{4\mathit{\pi r}\left(\rho -{\mathit{\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}r}-{\mathbb{R}}_2\frac{\left(k_4^2-k_1^2\right)}{\left(k_1^2-k_2^2\right)}{e}^{-k_{2}r}+\left({\mathit{Ik}}_3^2-{\mathbb{R}}_3\right)e^{-k_{3}r}\right]$$E31

with

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

which can be expressed as [36].

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

The fundamental solution of solid displacement ${\widehat{U}}^s\left(r\right)$ can be dismantled into singular ${\widehat{U}}_s^s\left(r\right)$ and regular ${\widehat{\mathrm{U}}}_{\mathrm{r}}^{\mathrm{s}}\left(\mathrm{r}\right)$ parts as

$$\begin{array}{cc}{\widehat{U}}^s\left(r\right)& ={\widehat{U}}_s^s\left(r\right)+{\widehat{U}}_r^s\left(r\right)\text{with}r≔\left|y-x\right|\\ & =\frac{1}{\mu }\left[I{\mathrm{\Delta }}_y-\frac{\lambda +\mu }{\lambda +2\mu }{\nabla }_y{\nabla }_y^T\right]{\mathrm{\Delta }}_y\widehat{x}\left(r\right)\\ & -\frac{1}{\mu }\left[\left(\left(k_1^2+k_2^2\right){\mathrm{\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(r\right)\end{array}$$E34

in which

$$\begin{array}{cc}\widehat{x}\left(r\right)& =\frac{1}{4\mathit{\pi r}}\left[\frac{e^{-k_{1}r}}{\left(k_2^2-k_1^2\right)\left(k_3^2-k_1^2\right)}+\frac{e^{-k_{2}r}}{\left(k_2^2-k_1^2\right)\left(k_2^2-k_3^2\right)}+\frac{e^{-k_{3}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(r^2\right)\end{array}$$E35

The remaining parts of $\widehat{U}\left(r\right)$ as in (30) can be described as [36].

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

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

$${\widehat{P}}^f\left(r\right)=\frac{s{\rho }^f}{4\mathit{\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\mathit{\pi r\beta }}+O\left(r^0\right)$$E38

On the basis of limiting process $\tilde{x}\in \Omega \to x\in \Gamma$ on (28) we get

$$\underset{\tilde{x}\in \Omega \to x\in \Gamma }{lim}{\left({\widehat{V}\widehat{t}}^g\right)}_{\Omega }\left(\tilde{x}\right)=\left({\widehat{V}\widehat{x}}^g\right)\left(x\right)≔\underset{\Gamma }{\overset{\cdot }{\int }}{\widehat{U}}^T\left(y-x\right){\widehat{t}}^g\left(y\right){\mathit{ds}}_y$$E39

According to limiting process $\tilde{x}\in \Omega \to x\in \Gamma$ on (28) we obtain [42].

$$\underset{\tilde{x}\in \Omega \to x\in \Gamma }{lim}{\left({\widehat{K}\widehat{u}}^g\right)}_{\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)$$E40

where

$$C\left(x\right)=\underset{\epsilon \to 0}{\mathit{lim}}\underset{y\in \Omega :\left|y-x\right|=\epsilon }{\overset{\cdot }{\int }}{\left({\widehat{T}}_y\widehat{U}\right)}^T\left(y-x\right){\mathit{ds}}_y$$E41

and

$$\left({\widehat{K}\widehat{u}}^g\right)\left(x\right)=\underset{\epsilon \to 0}{\mathit{lim}}\underset{\left|y-x\right|\ge \epsilon }{\overset{\cdot }{\int }}{\left({\widehat{T}}_y\widehat{U}\right)}^T\left(y-x\right){\widehat{u}}^g\left(y\right){\mathit{ds}}_y$$E42

By using (39)-(42), we can write

$$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)$$E43

By applying the inverse Laplace transform, we obtain

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

where $\ast$ is the time convolution.

According to [40], the fundamental solution is

$${\left({\widehat{T}}_y\widehat{U}\right)}^T={\left[\left[\begin{array}{cc}{\widehat{T}}_y^e& {\mathit{san}}_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]$$E45

On the basis of Stokes theorem, we obtain

$$\underset{\Gamma }{\overset{\cdot }{\int }}\left({\nabla }_y\times \mathrm{a},n_y\right){\mathit{ds}}_y=-\underset{\partial \Gamma }{\overset{,}{\int }}\left(\mathrm{a},v\right)d{\gamma }_y=-\underset{\varphi }{\overset{,}{\int }}\left(\mathrm{a},v\right)d{\gamma }_y=0$$E46

which can be expressed as

$$\underset{\Gamma }{\overset{\cdot }{\int }}\left(n_y\times {\nabla }_y,\mathrm{a}\right){\mathit{ds}}_y=0$$E47

On the basis of [40], we get

$$\underset{\Gamma }{\overset{\cdot }{\int }}\left(M_y\mathrm{a}\right){\mathit{ds}}_y=0$$E48

in which $M_y={\left({\nabla }_y{\nabla }_y^T\right)}^T-{\nabla }_y{\nabla }_y^T$.

By applying (48) to a formula $\mathrm{a}=\mathit{vu}$ we obtain [43].

$$\underset{\Gamma }{\overset{\cdot }{\int }}\left(M_{y}v\right){\mathit{uds}}_y=-\underset{\Gamma }{\overset{,}{\int }}v\left(M_{y}u\right){\mathit{ds}}_y$$E49

$$\underset{\Gamma }{\overset{\cdot }{\int }}{\left(M_{y}v\right)}^T{\mathit{uds}}_y=-\underset{\Gamma }{\overset{\cdot }{\int }}{v}^T\left(M_{y}u\right){\mathit{ds}}_y$$E50

Making use of (34) and (45), we can express ${\widehat{T}}^s$ as

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

On the basis of [40], we obtain

$${\left({\widehat{T}}^s\right)}^T=\left(\lambda +2\mu \right)n_y{\nabla }_y^T{\widehat{U}}_{\mathit{sin}g}^s-\mu \left(n_y\times \left({\nabla }_y\times {\widehat{U}}_{\mathit{sin}g}^s\right)\right)+2\mu M_y{\widehat{u}}_{\mathit{sin}g}^s+O\left(r^0\right)$$E52

which may be expressed using (34) as

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

By applying (29) (53), we obtain

$${\left(\widehat{k}\widehat{u}\right)}_{\mathrm{\Omega }}^s\left(\tilde{x}\right)={\int }_{\mathrm{\Gamma }}^{\cdot }\left[\left(M_y{\mathrm{\Delta }}_y^2\widehat{X}\right)\widehat{u}+\left(I\left({\mathrm{n}}^T{\nabla }_y\right){\mathrm{\Delta }}_y^2\widehat{X}\right)\widehat{u}+2\mu {\left(M_y{\widehat{U}}_{\mathit{sin}g}^s\right)}^T\widehat{u}+o\left(r^0\right)\widehat{u}\right]{\mathit{ds}}_y$$E54

Based on [42], we have

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

In the in right-side of (55), we can write second term as follows

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

in which

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

According to [40], we can write

$$\underset{\Omega \in \tilde{x}\to x\in \Gamma }{lim}{\left(\widehat{K}\widehat{u}\right)}_{\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)$$E58

By augmenting ${\widehat{U}}_s^s$ to ${\widehat{U}}^s$, we obtain

$${\left(\widehat{k}\widehat{u}\right)}_{\mathrm{\Omega }}^s\left(\tilde{x}\right)=\underset{\Gamma }{\overset{\cdot }{\int }}-{\mathrm{\Delta }}_y^2\widehat{x}\left(M_y\widehat{u}\right)+\left(I\left(n^T{\nabla }_y\right){\mathrm{\Delta }}_y^2\widehat{x}\right)\widehat{u}+2\mu {\widehat{U}}^s\left(M_y\widehat{u}\right)+O\left(r^0\right)\widehat{u}{\mathit{ds}}_y$$E59

According to [41], we get

$$\left(f\ast g\right)\left(t\right)=\underset{0}{\overset{t}{\int }}f\left(t-\tau \right)g\left(\tau \right)\mathit{d\tau }\text{for}t\in \left[0,T\right]$$E60

where

$$\left(f\ast g\right)\left(t_n\right)\approx \sum _{k=0}^n{\omega }_{n-k}^{\mathrm{\Delta }t}\left(\widehat{f}\right)g\left(t_k\right)$$E61

On the basis of Lubich [44, 45], the integration weights ${\omega }_n$ are calculated using Cauchy’s integral formula as

$${\omega }_{n-k}^{\mathrm{\Delta }t}\left(\widehat{f}\right)≔\frac{1}{2\mathit{\pi i}}\underset{\left|z\right|=R}{\int }\widehat{f}\left(\frac{\gamma \left(z\right)}{\mathrm{\Delta }t}\right)z^{-\left(n+1\right)}\mathit{dz}$$E62

Polar coordinate transformation $z={\mathit{Re}}^{-\mathit{i\phi }}$ is used with the trapezoidal rule to approximate the integral (62) as

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

Substitution of Eq. (63) into Eq. (61), we get

$$\begin{array}{cc}\left(f\ast g\right)\left(t_n\right)& \approx \sum _{k=0}^N\frac{R^{-\left(n-k\right)}}{N+1}\sum _{\mathrm{\ell }=0}^N\widehat{f}\left(s_{\mathrm{\ell }}\right){\zeta }^{\mathrm{\ell }\left(n-k\right)}g\left(t_k\right)\\ & \approx \frac{R^{-n}}{N+1}\sum _{\mathrm{\ell }=0}^N\widehat{f}\left(s_{\mathrm{\ell }}\right)\widehat{g}\left(s_{\mathrm{\ell }}\right){\zeta }^{\mathrm{\ell }n}\end{array}$$E64

with

$$\widehat{g}\left(s_{\mathrm{\ell }}\right)=\sum _{k=0}^N{R}^{k}g\left(t_k\right){\zeta }^{-\mathrm{\ell }k}.$$E65

According to the procedure [43], the convolution operator (44) can be expressed as

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

which may be written as

$$C\left(x\right){\widehat{u}}^g\left(x,s_{\mathrm{\ell }}\right)=\left({\widehat{v}\widehat{t}}^g\right)\left(x,s_{\mathrm{\ell }}\right)-\left({\widehat{k}\widehat{u}}^g\right)\left(x,s_{\mathrm{\ell }}\right),\mathrm{\ell }=0\dots \dots ..N$$E67

Let the boundary $\Gamma$ is discretized into $N_e$ boundary elements ${\tau }_e$ as follows

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

Now, we assume that we have

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

$$S_h\left[k\right]\left({\Gamma }_{D,h}\right)≔\mathit{\text{span}}{\left\{{\psi }_i^{\beta }\left[k\right]\right\}}_{j=1}^{\mathbb{j}},\beta \ge 0$$E70

where

$${\widehat{u}}^g\left[k\right]\left(x\right)\approx {\widehat{u}}_h^g\left[k\right]\left(x\right)=\sum _{i=1}^{\mathbb{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)$$E71

$${\widehat{t}}^g\left[k\right]\left(x\right)\approx {\widehat{t}}_h^g\left[k\right]\left(x\right)=\sum _{j=1}^{\mathbb{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)$$E72

where $k=1,2,3,4$ are the poro-elastic degrees of freedom, ${\phi }_i^{\alpha }\left[k\right]$ are $\mathbb{i}$ continuous polynomial shape functions and ${\psi }_i^{\beta }\left[k\right]$ are $\mathbb{j}$ piecewise discontinuous polynomial shape functions.

Thus, from (67), we can write the following $N+1$ algebraic systems of equations

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

5. Numerical results and discussion

In the current study, a Krylov subspace iterative method is used for solving the resulting linear systems. In order to reduce the number of iterations, a dual threshold incomplete LU factorization technique (ILUT) which is one of the well-known preconditioning techniques is implemented as a robust preconditioner for TFQMR (Transpose-free quasi minimal residual) [46] to accelerate the convergence of the solver TFQMR.

To illustrate the numerical calculations computed by the proposed technique, the physical parameters for transversely isotropic soft tissue are given as follows [47]:

The elasticity tensor

$$C_{\mathit{ab}\mathit{lg}}=\left[\begin{array}{llllll}{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]$$E74

in which

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

where

$$v=0.196,v_0=0.163,{\mu }_0=20.98\mathrm{GPa},E=68.34\mathrm{GPa},E_0=51.35\mathrm{GPa}$$E76

and therefore

$$k_1=108.39\mathrm{GPa},k_2=-21.70\mathrm{GPa}$$E77

where $E$ and $E_0$ are the respectively, $v$ and $v_0$ are Poisson’s ratio in the isotropy plane and in the fiber direction respectively, and ${\mu }_0$ is the shear moduli in any direction within a plane perpendicular to isotropy plane.

Since for strongly anisotropic soft tissue both bulk moduli are positive, we used the following physical parameters for anisotropic soft tissue [48].

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

and therefore

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

and other constants considered in the calculations are as follows.

${\rho }_s=1600\mathit{kg}/m^3$,${\rho }_{\mathcal{F}}=1113\mathit{kg}/m^3$, $p=25\mathit{\text{MPap}}=25\mathit{MPa}$,$\varphi =0.15$ and $Q/R=0.65$.(80)

The domain boundary of the current problem has been discretized into 21 boundary elements and 42 internal points as depicted in Figure 2. The computation was done using Matlab R2018a on a MacBook Pro with 2.9GHz quad-core Intel Core i7 processor and 16GB RAM.

Figure 2.

Boundary element model of the current problem.

Figure 3 shows the variation of the temperature $\mathrm{T}$ along $\mathrm{x}‐\text{axis}$ for different values of fractional order parameter. It can be seen from this figure that the fractional order parameter has a significant influence on the temperature.

Figure 3.

Variation of the temperature $\mathit{T}$ along $\mathit{x}‐\mathit{\text{axis}}$.

Figure 4 illustrates the variation of the displacement ${\mathrm{u}}_1$ along $\mathrm{x}‐\text{axis}$ for different values of fractional order parameter. It can be seen from this figure that the fractional order parameter has a significant influence on the displacement ${\mathrm{u}}_1$.

Figure 4.

Variation of the displacement ${\mathit{u}}_{\mathbf{1}}$ along $\mathit{x}‐\mathit{\text{axis}}$.

Figure 5 shows the variation of the displacement ${\mathrm{u}}_2$ along $\mathrm{x}‐\text{axis}$ for different values of fractional order parameter. It can be seen from this figure that the fractional order parameter has a significant influence on the displacement ${\mathrm{u}}_2$.

Figure 5.

Variation of the displacement ${\mathit{u}}_{\mathbf{2}}$ along $\mathit{x}‐\mathit{\text{axis}}$.

Figure 6 shows the variation of the fluid pressure $\mathrm{p}$ along $\mathrm{x}‐\text{axis}$ for different values of fractional order parameter. It can be seen from this figure that the fractional order parameter has a significant influence on the fluid pressure $\mathrm{p}$.

Figure 6.

Variation of the fluid pressure $\mathit{p}$ along $\mathit{x}‐\mathit{\text{axis}}$.

Figure 7 shows the variation of the bio-thermal stress ${\mathrm{\sigma }}_{11}$ along $x‐\text{axis}$ for different values of fractional order parameter. It can be seen from this figure that the fractional order parameter has an important influence on the bio-thermal stress ${\mathrm{\sigma }}_{11}$.

Figure 7.

Variation of the bio-thermal stress ${\mathit{\sigma }}_{\mathbf{11}}$ along $\mathit{x}‐\mathit{\text{axis}}$.

Since there are no findings available for the problem under consideration. Therefore, some literatures may be regarded as special cases from our general problem. In the special case under consideration, the results of the bio-thermal stress caused by coupling between the temperature and displacement fields are plotted in Figure 8 to illustrate the variation of the bio-thermal stress ${\mathrm{\sigma }}_{11}$ along $x‐\text{axis}$ for BEM, FDM and FEM, where the boundary of the special case problem has been discretized into 21, 42 and 84 boundary elements (bes). The validity, accuracy and efficiency of our proposed technique have been confirmed by a graphical comparison of the three different boundary elements (21, 42 and 84) with those obtained using the FDM results of Shen and Zhang [49] and FEM results of Torvi and Dale [50] for the special case under consideration, the increase of BEM boundary elements leads to improve the accuracy and efficiency of the BEM, also, it can be noted that the BEM findings are in excellent agreement with the FDM and FEM results, we refer the interested reader to recent work [51, 52, 53, 54, 55] for understanding the BEM methodology.

Figure 8.

Variation of the bio-thermal stress ${\mathit{\sigma }}_{\mathbf{11}}$ along $\mathit{x}‐\mathit{\text{axis}}$ for BEM, FDM and FEM.

6. Conclusion

1. A novel boundary element model based on the TWMBT and Biot’s theory was established for describing the bio-thermomechanical interactions in anisotropic soft tissues.

2. The bio-heat transfer equation has been solved using the dual reciprocity boundary element method (DRBEM) to obtain the temperature distribution.

3. The mechanical equation has been solved using the convolution quadrature boundary element method (CQBEM) to obtain the displacement and fluid pressure for different temperature distributions at each time step.

4. Due to the advantages of DRBEM and CQBEM such as dealing with more complex shapes of soft tissues and not needing the discretization of the internal domain, also, they have low RAM and CPU usage. Therefore, they are a versatile and powerful methods for modeling of fractional bio-thermomechanical problems in anisotropic soft tissues.

5. The linear systems resulting from BEM have been solved by TFQMR solver with the ILUT preconditioner which reduces the number of iterations and the total CPU time.

6. Numerical findings are presented graphically to show the effect of fractional order parameter on the problem variables temperature, displacements and fluid pressure.

7. Numerical findings confirm the validity, efficiency and accuracy of the proposed BEM technique.

8. The proposed technique can be applied to a wide variety of fractional bio-thermomechanical problems in anisotropic soft tissues.

9. For open boundary problems of soft tissues, such as the considered problem, the BEM users need only to deal with real geometry boundaries. But for these problems, FDM and FEM use artificial boundaries, which are far away from the real soft tissues. Also, these artificial boundaries are also becoming a big challenge for FDM users and FEM users. So, BEM becomes the best method for the considered problem.

10. The presence of fractional order parameter in the current study plays a significant role in all the physical quantities during modeling and simulation in medicine and healthcare.

11. From the research that has been performed, it is possible to conclude that the proposed BEM is an easier, effective, predictable, and stable technique in the treatment of the bio-thermomechanical soft tissue models.

12. It can be concluded from this chapter that Biot’s equations for the dynamic response of poroelastic media can be combined with the bio-heat transfer models to describe the fractional bio-thermomechanical interactions of anisotropic soft tissues.

13. Current numerical results for our complex and general problem may provide interesting information for researchers and scientists in bioengineering, heat transfer, mechanics, neurophysiology, biology and clinicians.

Nomenclature

A=ϕ1+Q/R	Biot’s coefficient
Bx˜e	linear elastostatics operator
Γ	considered boundary
ΓD	Dirichlet boundary
ΓN	Neumann boundary
C	specific heat of soft tissue
C	shape factor
Cb	specific heat of the blood
Cpjkl	specific heat of the blood
F	bulk body forces
ĝD	Dirichlet datum
ĝN	Neumann datum
K	dynamic permeability
K	thermal conductivity of soft tissue
m	iterative parameter
p	pore pressure
P0τ	heating power
Q,R	solid–fluid coupling parameters
Qmet	metabolic heat source
Qext	external heat source
S0	scattering coefficient
T	soft tissue temperature
T⌣b	arterial blood temperature
Txe	traction derivative
u	solid displacement
uf	fluid displacement
V=Vf+Vs	bulk volume
Vf	fluid volume
Vs	solid volume
Wb	blood perfusion rate
B	stress-temperature coefficients
∈	linear strain tensor
ζ	fluid volume variation
ρ=ρs1−ϕ+ϕρf	bulk density
ρs=Cϕρf	mass density of soft tissue
ρf	blood density
σ	total stress tensor
τ	time
τq	phase lag for heat flux
τT	phase lag for temperature gradient
φiαk	continuous polynomial shape functions
ϕ=VfV	porosity
ψjβk	discontinuous polynomial shape functions
Ω	considered region