A Computational Model for Nonlinear Biomechanics Problems of FGA Biological Soft Tissues

Abstract

The principal objective of this work was to develop a semi-implicit hybrid boundary element method (HBEM) to describe the nonlinear fractional biomechanical interactions in functionally graded anisotropic (FGA) soft tissues. The local radial basis function collocation method (LRBFCM) and general boundary element method (GBEM) were used to solve the nonlinear fractional dual-phase-lag bioheat governing equation. The boundary element method (BEM) was then used to solve the poroelastic governing equation. To solve equations arising from boundary element discretization, an efficient partitioned semi-implicit coupling algorithm was implemented with the generalized modified shift-splitting (GMSS) preconditioners. The computational findings are presented graphically to display the influences of the graded parameter, fractional parameter, and anisotropic property on the bio-thermal stress. Different bioheat transfer models are presented to show the significant differences between the nonlinear bio-thermal stress distributions in functionally graded anisotropic biological tissues. Numerical findings verified the validity, accuracy, and efficiency of the proposed method.

1. Introduction

Because of advancements in microwave, laser, focused ultrasound, and radiofrequency technologies, many modern thermo-therapeutics are now widely used in clinical treatment. For example, in thermal therapy, an objective lens focuses a laser on a tumor. One of the most difficult problems in thermal therapy is delivering the right amount of heat energy to the diseased tissue while avoiding damage to healthy tissue. As a result, there is a pressing need to comprehend how temperature/stress fields influence the kinetics of thermal therapy. Van and Gybels [1] demonstrated that deformation caused by heating and cooling can cause pain sensation. Accurate predictions of heat and mechanical reactions, as well as thermal damage in biological tissue, are thus essential for treatment planning and the development of novel therapeutic heating systems. Because of the inherent properties of biological tissue, heat transfer analysis in living biological tissue is a complex physiological procedure. Heat transfer analysis in living biological tissue is a difficult physiological process due to biological tissue’s inherent properties, such as blood circulation, sweating, metabolic heat generation, and heat dissipation via hair or fur. Pennes [2] used his bioheat transfer model to simulate the temperature profile in the human forearm to adequately characterize this complex event. Other bioheat transfer models were developed to address the shortcomings of Pennes’ model, which does not consider the blood velocity field or the geometry of blood vessels [3,4,5,6,7]. Although these models are more detailed and precise than Pennes’, their complexities make them difficult to apply in practice. On the other hand, because Pennes’ model is simple and has a small number of material parameters, it has piqued the interest of many researchers [8]. Because of its non-homogeneous structure and composition, soft tissue’s mechanical response in vivo is heterogeneous, anisotropic, nonlinear, and viscoelastic. It is influenced by a variety of factors, including age, gender, location, and hydration. Cattaneo [9] and Vernotee [10] proposed a single-phase lag (SPL) related to heat flow to address a limitation of the linear bioheat transfer model. The thermal wave bioheat transfer (TWBHT) model was obtained by combining the SPL model with the energy equation. Tzou [11,12] extended the SPL model by incorporating another phase lag caused by a temperature gradient, resulting in the dual-phase-lag (DPL) model. Furthermore, standard electrometric constitutive models are insufficient for describing soft tissue’s complex mechanical behavior. A thorough understanding of the thermal [13,14,15,16] and biothermo–mechanical [17,18,19] responses of biological tissues will aid in the design and optimization of heat transfer processes in biological tissues. Several models were proposed to study the mechanical behavior of biological tissues [20,21]. Pioletti and Rakotomanana [22] studied the non-linear viscoelastic laws for soft biological tissues. A generic physics-informed neural-network-based constitutive model for soft biological tissues was established by Liu et al. [23]. Furthermore, Miller and Gasser [24] studied the non-linear and time-dependent properties of soft biological tissues containing collagen.

The primary goal of this study was to introduce the hybrid boundary element method (HBEM) scheme for describing nonlinear fractional bio-thermomechanical interactions in functionally graded anisotropic (FGA) soft tissues, which is based on the local radial basis function collocation method (LRBFCM) and GBEM to solve the bioheat governing equation and on the boundary element method (BEM) to solve the poroelastic governing equation. The calculation results were graphed to show the effects of the anisotropic property and graded and fractional parameters on bio-thermal stress. These findings also supported the proposed hybrid scheme’s validity and efficiency.

2. Formulation of the Problem

Let us consider two-dimensional $\left\{x=\left(x_1,x_2\right)=\left(x,y\right):0\le x_1=x\le 1, 0\le x_2=y\le 1\right\}$ functionally graded tissue with a boundary $\mathsf{\Gamma }$. The governing equations that model the nonlinear fractional bio-thermomechanical interactions in FGA soft tissues are as follows [25,26]:

$${\sigma }_{ij,j}+\rho F_i=\rho {\ddot{u}}_i+\varphi {\rho }_{ℱ}{\ddot{v}}_i$$

(1)

where ${\sigma }_{ij}$, $\rho$, ${\rho }_{ℱ}$, $F_i$, $\varphi$, $u_i$, and $v_i$ are the mechanical stress tensor, bulk density, fluid density, bulk body forces, porosity, solid displacement, and fluid–solid displacement, respectively.

$$\dot{\zeta }+q_{i,i}={\overline{ℂ}}_i$$

(2)

where $\zeta$, $q$, and ${\overline{ℂ}}_i$ are the fluid volume variation, instantaneous flux, and source term, respectively.

Furthermore:

$${\sigma }_{ij}={\left(x+1\right)}^m\left[C_{ijkl} e{\delta }_{ij}-A{\delta }_{ij}p-{\beta }_{ij} T\right]$$

(3)

$$q_i=-\overline{k}{\left(x+1\right)}^m\left(p_{,i}+{\rho }_{ℱ}{\ddot{u}}_i+\frac{{\rho }_0+\varphi {\rho }_{ℱ}}{\varphi }{\ddot{v}}_i\right), \zeta =\left[A{u}_{k,k}+\frac{{\varphi }^2}{R}p \right]$$

(4)

in which ${\epsilon }_{ij}=\frac{1}{2}\left(u_{i,j}+u_{j,i}\right)$, $e={\epsilon }_{ii}$, and $A=\varphi \left(1+\frac{Q}{R}\right)$. $T$, $\lambda$, $C_{ijkl}$, $A$, $p$,${\beta }_{ij}$, $\overline{k}$, $T_0$, $Å$, and $m$ are the temperature, thermal conductivity, constant elastic moduli, Biot’s effective stress coefficient, fluid pressure, stress–temperature coefficients, permeability, reference temperature, unified parameter, and functionally graded parameter, respectively. $Q$ and $R$ are solid–fluid coupling parameters, ${\rho }_0=\eta \varphi {\rho }_{ℱ}$, and $\eta$ is the shape factor.

The fractional order bioheat transfer equation without dual-phase lag can be expressed as

$$D_{\tau }^{a}T\left(x, \tau \right)=\xi \nabla \left[\check{K} \nabla T\left(x, \tau \right)\right]+Q_m\left(x, \tau \right), \xi =\frac{1}{c\rho }$$

(5)

where $a$, $c$, $\rho$, $\tau$, $\check{K}$, and $Q_m$ are the fractional parameter, specific heat, density, time, thermal conductivity, and metabolic heat production, respectively. The schematic flowchart representation of the considered problem is described in Figure 1.

3. Hybrid Technique Implementation

3.1. LRBFCM—GBEM Implementation for the Temperature Field

3.1.1. LRBFCM Implementation for the Time-Fractional-Order Bioheat Equation without Dual-Phase Lag

According to the finite difference scheme of Fahmy [16], the temperature time derivative can be expressed as

$$\dot{T}\left(x, \tau \right)=\frac{T^{f+1}\left(x\right)-T^f\left(x\right)}{\Delta \tau }+O\left(\Delta \tau \right), {\tau }^f=f\Delta \tau , f=0, 1, 2, \dots , F, \dot{T}\left(x, \tau \right)\in \left[{\tau }^f, {\tau }^{f+1}\right]$$

where $D_{\tau }^a$ is the Caputo derivative of order $a$, which can be written as [9]

$$D_{\tau }^{a}T\left(x, \tau \right)=\frac{1}{\mathsf{\Gamma }\left(1-a\right)}\underset{0}{\overset{\tau }{{\displaystyle \int }}}\frac{\partial T\left(r, s\right)}{\partial s}\frac{ds}{{\left(\tau -s\right)}^a}, 0<a<1$$

(6)

Based on the Caputo derivative (6), the following formula can be presented:

$$D_{\tau }^a{T}^{f+1}+D_{\tau }^a{T}^f\approx {\displaystyle \sum }_{j=0}^k{W}_{a,j}\left(T^{f+1-j}\left(x\right)-T^{f-j}\left(x\right)\right),\left(f=1, 2, \dots , F\right)$$

(7)

in which

$$W_{a,0}=\frac{{\left(\Delta \tau \right)}^{-a}}{\mathsf{\Gamma }\left(2-a\right)}$$

(8)

$$W_{a,j}=W_{a,0}\left({\left(j+1\right)}^{1-a}-{\left(j-1\right)}^{1-a}\right), j=1, 2, \dots , F$$

(9)

According to Equation (7), we can write the bioheat transfer as Equation (5):

$$\begin{array}{c}{W}_{a,0}{T}^{f+1}\left(\mathit{x}\right)-\check{K}\left(\mathit{x}\right)T_{,II}^{f+1}\left(\mathit{x}\right)-{\check{K}}_{,I}\left(\mathit{x}\right)T_{,I}^{f+1}\left(\mathit{x}\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} =W_{a,0}{T}^f\left(\mathit{x}\right)-\check{K}\left(\mathit{x}\right)T_{,II}^f\left(\mathit{x}\right)-{\check{K}}_{,I}\left(\mathit{x}\right)T_{,J}^f\left(\mathit{x}\right)-{\displaystyle {\displaystyle \sum }_{j=1}^f}{W}_{a,j}\left(T^{f+1-j}\left(\mathit{x}\right)-T^{f-j}\left(\mathit{x}\right)\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} +Q_m^{f+1}\left(\mathit{x},\tau \right)+Q_m^f\left(\mathit{x},\tau \right),f=0,1,2,\dots ,F\hfill \end{array}$$

(10)

Krahulec et al. [27] used the LRBFCM to solve (10) as follows:

$$\begin{array}{c}\left(W_{a,0}E-\mathsf{\Lambda }{\check{\mathsf{\Phi }}}_{II}{\check{\mathsf{\Phi }}}^{-1}-{\mathsf{\Lambda }}_I{\check{\mathsf{\Phi }}}_I{\check{\mathsf{\Phi }}}^{-1}\right){\check{T}}^{f+1}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} =\left(W_{a,0}E+\mathsf{\Lambda }{\check{\mathsf{\Phi }}}_{II}{\check{\mathsf{\Phi }}}^{-1}+{\mathsf{\Lambda }}_I{\check{\mathsf{\Phi }}}_I{\check{\mathsf{\Phi }}}^{-1}\right){\check{T}}^f-{\displaystyle {\displaystyle \sum }_{j=1}^f}{W}_{a,j}\left(T^{f+1-j}\left(x\right)-T^{f-j}\left(x\right)\right)+Q_m^{f+1}\left(\mathit{x}, \tau \right)+Q_m^f\left(\mathit{x}, \tau \right), f\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} =0, 1, 2, \dots , F\hfill \end{array}$$

3.1.2. GBEM Implementation for the Dual-Phase-Lag Bioheat Equation without a Fractional-Order Derivative

The dual-phase-lag bioheat equation in the absence of a fractional derivative may be written as [28]

$$c\rho \left[\frac{\partial \mathit{T}}{\partial \tau }+{\tau }_q\frac{{\partial }^2\mathit{T}}{\partial {\tau }^2}\right]=\check{K}{\nabla }^2\mathit{T}+\check{K}{\tau }_{\mathit{T}}\frac{\partial }{\partial \tau }\left({\nabla }^2\mathit{T}\right)+W_b{C}_b\left({\breve{\mathit{T}}}_b-\mathit{T}\right)+Q_m-W_b{C}_b{\tau }_q\frac{\partial \mathit{T}}{\partial \tau }$$

(11)

where ${\tau }_q$, ${\tau }_T$, $W_b$, $C_b$, and ${\breve{\mathit{T}}}_b$ are the heat-flux phase lag, temperature-gradient phase lag, perfusion rate, specific heat, and arterial temperature, respectively.

In the proposed model, we considered the following conditions:

$$\mathit{T}\left(\mathit{x},0\right)=T_0, {\frac{\partial \mathit{T}\left(\mathit{x},\tau \right)}{\partial \tau }|}_{\tau =0}=0$$

(12)

$$\mathit{T}\left(\mathit{x},\tau \right)={\mathit{T}}_b\left(\mathit{x},\tau \right) \mathrm{for} \mathit{x}\in {\mathsf{\Gamma }}_1$$

(13)

$${\mathit{q}}_b\left(\mathit{x},\tau \right)+{\tau }_q\frac{\partial {\mathit{q}}_b\left(\mathit{x},\tau \right)}{\partial \tau }=-\check{K}\left[\frac{\partial \mathit{T}\left(\mathit{x},\tau \right)}{\partial \mathit{n}}+{\tau }_{\mathit{T}}\frac{\partial }{\partial \tau }\left(\frac{\partial \mathit{T}\left(\mathit{x},\tau \right)}{\partial \mathit{n}}\right)\right] \mathrm{for} \mathit{x}\in {\mathsf{\Gamma }}_2$$

(14)

where $T^0\left(\mathit{x}\right)=T^1\left(\mathit{x}\right)=T_0$ and ${\tau }^{f-1}\to {\tau }^f \left(f\ge 2\right)$.

Thus, Equation (11) may be approximated as follows [28]:

$$\begin{array}{c}c\rho \left(\frac{{\mathit{T}}^f\left(\mathit{x}\right)-{\mathit{T}}^{f-1}\left(\mathit{x}\right)}{\Delta \tau }+{\tau }_q\frac{{\mathit{T}}^f\left(\mathit{x}\right)-2{\mathit{T}}^{f-1}\left(\mathit{x}\right)+{\mathit{T}}^{f-2}\left(\mathit{x}\right)}{{\left(\Delta \tau \right)}^2}\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\check{K}{\nabla }^2{\mathit{T}}^f\left(\mathit{x}\right)+\frac{\check{K}{\tau }_{\mathit{T}}}{\Delta \tau }\left[{\nabla }^2{\mathit{T}}^f\left(\mathit{x}\right)-{\nabla }^2{\mathit{T}}^{f-1}\left(\mathit{x}\right)\right]+W_b{C}_b\left[{\breve{\mathit{T}}}_b-{\mathit{T}}^f\left(\mathit{x}\right)\right]+Q_m\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-W_b{C}_b{\tau }_q\frac{{\mathit{T}}^f\left(\mathit{x}\right)-{\mathit{T}}^{f-1}\left(\mathit{x}\right)}{\Delta \tau }\hfill \end{array}$$

(15)

Equation (15) may be expressed as

$${\nabla }^2{\mathit{T}}^f\left(\mathit{x}\right)-B{\mathit{T}}^f\left(\mathit{x}\right)+C{\nabla }^2{\mathit{T}}^{f-1}\left(\mathit{x}\right)+D{\mathit{T}}^{f-1}\left(\mathit{x}\right)+E{\mathit{T}}^{f-2}\left(\mathit{x}\right)+F=0$$

(16)

in which

$$B=\frac{\left(c\rho +W_b{C}_b\Delta \tau \right)\left(\Delta \tau +{\tau }_q\right)}{\check{K}\Delta \tau \left(\Delta \tau +{\tau }_T\right)}, C=\frac{{\tau }_T}{\Delta \tau +{\tau }_T}, D=\frac{c\rho \left(\Delta \tau +2{\tau }_q\right)+W_b{C}_b{\tau }_q\Delta \tau }{\check{K}\Delta \tau \left(\Delta \tau +{\tau }_T\right)},$$

$$E=-\frac{c\rho {\tau }_q}{\check{K}\Delta \tau \left(\Delta \tau +{\tau }_T\right)}, F=\frac{\Delta \tau \left(W_b{C}_b{\breve{T}}_b+Q_m\right)}{\check{K}\left(\Delta \tau +{\tau }_T\right)}$$

Equations (13) and (14) may be expressed as

$${\mathit{T}}^f\left(\mathit{x}\right)={\mathit{T}}_b^f\left(\mathit{x}\right) \mathrm{for} \mathit{x}\in {\mathsf{\Gamma }}_1$$

(17)

$${\mathit{Z}}^f\left(\mathit{x}\right)=w_b^f\left(\mathit{x}\right)=-\check{K}\frac{\partial {\mathit{T}}^f\left(\mathit{x}\right)}{\partial \mathit{n}} \mathrm{for} \mathit{x}\in {\mathsf{\Gamma }}_2$$

(18)

where

$$w_b^f\left(\mathit{x}\right)=\frac{\Delta \tau }{\Delta \tau +{\tau }_T}\left(q_b^f\left(\mathit{x}\right)+{\tau }_q{\frac{\partial q_b\left(\mathit{x},\tau \right)}{\partial \tau }|}_{t=t^f}\right)-\check{K}\frac{{\tau }_T}{\Delta \tau +{\tau }_T}\frac{\partial {\mathit{T}}^{f-1}\left(\mathit{x}\right)}{\partial n}$$

By using the same technique used by Fahmy [16], the following boundary integral equation can be obtained:

$$B\left(\xi \right){\mathit{U}}^{\left(1\right)}\left(\xi \right)+\frac{1}{\check{K}}{{\displaystyle \int }}_{\mathsf{\Gamma }}^{}{\mathit{T}}^{*}\left(\xi ,\mathit{x}\right){\mathit{w}}^{\left(1\right)}\left(\mathit{x}\right)d\mathsf{\Gamma }=\frac{1}{\check{K}}{{\displaystyle \int }}_{\mathsf{\Gamma }}^{}{\mathit{q}}^{*}\left(\xi ,\mathit{x}\right){\mathit{U}}^{\left(1\right)}\left(\mathit{x}\right)d\mathsf{\Gamma }+{{\displaystyle \int }}_{\mathsf{\Omega }}^{}R\left[{\mathit{T}}_{k-1}^f\left(\mathit{x}\right)\right]{\mathit{T}}^{*}\left(\xi ,\mathit{x}\right)d\mathsf{\Omega }$$

(19)

The fundamental solutions are

$${\mathit{T}}^{*}\left(\xi ,\mathit{x}\right)=\frac{1}{4\pi r}\exp \left(-r\sqrt{B}\right)$$

(20)

$${\mathit{q}}^{*}\left(\xi ,\mathit{x}\right)=\frac{\check{K}d}{4\pi r^2}\exp \left(-r\sqrt{B}\right)\left(\frac{1}{r}+\sqrt{B}\right), d={\displaystyle \sum }_{e=1}^3\left({\mathit{x}}_e-{\xi }_e\right)\cos {\alpha }_e$$

(21)

After the discretization process, Equation (19) is approximated as follows [16]:

$${\displaystyle \sum }_{j=1}^N{G}_{ij}{W}^{\left[1\right]}\left({\mathit{x}}_j\right)={\displaystyle \sum }_{j=1}^N{H}_{ij}{\mathit{U}}^{\left[1\right]}\left({\mathit{x}}_j\right)+{\displaystyle \sum }_{l=1}^L{P}_{il}R\left[{\mathit{T}}_{k-1}^f\left({\mathit{x}}_l\right)\right]$$

(22)

in which

$$G_{ij}=\frac{1}{\check{K}}{{\displaystyle \int }}_{{\mathsf{\Gamma }}_j}^{}{\mathit{T}}^{*}\left({\xi }_i,\mathit{x}\right)d{\mathsf{\Gamma }}_j$$

$$H_{ij}=\{\begin{array}{c}{{\displaystyle \int }}_{{\mathsf{\Gamma }}_j}^{}{\mathit{q}}^{*}\left({\xi }_i,\mathit{x}\right)d{\mathsf{\Gamma }}_j, i\ne j\\ -0.5, i=j\end{array}$$

$$P_{il}={{\displaystyle \int }}_{{\mathsf{\Omega }}_j}^{}{\mathit{T}}^{*}\left({\xi }_i,\mathit{x}\right)d{\mathsf{\Omega }}_j$$

From (22), we obtain the boundary unknowns ${\mathit{W}}^{\left(1\right)}$ and ${\mathit{U}}^{\left(1\right)}$. Then, we can calculate ${\mathit{U}}^{\left(1\right)}\left({\xi }_i\right)$ using

$${\mathit{U}}^{\left(1\right)}\left({\xi }_i\right)={\displaystyle \sum }_{j=1}^N{H}_{ij}{\mathit{U}}^{\left(1\right)}\left({\mathit{x}}_j\right)-{\displaystyle \sum }_{j=1}^N{G}_{ij}{\mathit{W}}^{\left(1\right)}\left({\mathit{x}}_j\right)+{\displaystyle \sum }_{j=1}^N{P}_{il}R\left[{\mathit{T}}_{k-1}^f\left({\mathit{x}}_l\right)\right]$$

(23)

3.2. BEM Implementation for a Poroelastic Displacement Field

Based on the weighted residual methodology, Equations (1) and (2) can be written as

$${{\displaystyle \int }}_R\left({\sigma }_{ij,j}+U_i\right) u_i^{*} dR=0$$

(24)

$${{\displaystyle \int }}_R\left(q_{i,i}+{\dot{\zeta }}_i-{ℂ}_i\right) p_i^{*} dR=0$$

(25)

where

$${\sigma }_{ij,j}={\left(x+1\right)}^m\left[C_{ijkl} u_{k,lj}-A{\delta }_{ij}{p}_{,j}-{\beta }_{ij} \left(T_{,j}+{\tau }_1{\dot{T}}_{,j}\right)\right]++\frac{m}{x+1}{\sigma }_{ij}$$

$$q_{i,i}=-\overline{k}{\left(x+1\right)}^m\left(p_{,ii}+{\rho }_{ℱ}{\ddot{u}}_{i,i}+\frac{{\rho }_0+\varphi {\rho }_{ℱ}}{\varphi }{\ddot{v}}_{i,i}\right)+\frac{m}{x+1}{q}_i$$

in which $U_i=\rho F_i-\rho {\ddot{u}}_i-\varphi {\rho }_{ℱ}{\ddot{v}}_i$, and $u_i^{*}$ and $p_i^{*}$ are weighting functions.

For the first terms of (24) and (25), integration by parts yields

$$-{{\displaystyle \int }}_R{\sigma }_{ij} u_{i,j}^{*} dR+{{\displaystyle \int }}_R{U}_i u_i^{*} dR=-{{\displaystyle \int }}_{S_2}{\lambda }_i u_i^{*} dS$$

(26)

$$-{{\displaystyle \int }}_{R}q p_{i,j}^{*} dR+{{\displaystyle \int }}_R{\dot{\zeta }}_i p_i^{*} dR-{{\displaystyle \int }}_R{ℂ}_i p_i^{*} dR=-{{\displaystyle \int }}_{S_4}{L}_i p_i^{*} dS$$

(27)

According to Fahmy [29], we can write

$${{\displaystyle \int }}_R{\sigma }_{ij,j}^{*} u_i dR=-{{\displaystyle \int }}_S{u}_i^{*} {\lambda }_i dS-{{\displaystyle \int }}_S{p}_i^{*} L_i dS+{{\displaystyle \int }}_S{\lambda }_i^{*} u_i dS+{{\displaystyle \int }}_S{L}_i^{*} p_i dS$$

(28)

which may be written as

$$C^n{𝕢}^n=-{{\displaystyle \int }}_S{𝕡}^{*} 𝕢dS+{{\displaystyle \int }}_S{𝕢}^{*} 𝕡dS+{{\displaystyle \int }}_S{𝕒}^{*} p dS+{{\displaystyle \int }}_S{𝕓}^{*} \frac{\partial p}{\partial n} dS$$

(29)

in which

$$C^n=\left[\begin{array}{ccc}{C}_{11}& C_{12}& \\ C_{21}& C_{22}& \\ & & \end{array}\right], {𝕢}^{*}=\left[\begin{array}{ccc}{u}_{11}^{*}& u_{12}^{*}& {\omega }_{13}^{*}\\ u_{21}^{*}& u_{22}^{*}& {\omega }_{23}^{*}\\ u_{31}^{**}& u_{32}^{**}& {\omega }_{33}^{**}\end{array}\right], {𝕡}^{*}=\left[\begin{array}{ccc}{\lambda }_{11}^{*}& {\lambda }_{12}^{*}& {\mu }_{13}^{*}\\ {\lambda }_{21}^{*}& {\lambda }_{22}^{*}& {\mu }_{23}^{*}\\ {\lambda }_{31}^{**}& {\lambda }_{32}^{**}& {\mu }_{33}^{**}\end{array}\right]$$

$$𝕢=\left[\begin{array}{c}{u}_1\\ u_2\\ {\omega }_3\end{array}\right], 𝕡=\left[\begin{array}{c}{\lambda }_1\\ {\lambda }_2\\ {\mu }_3\end{array}\right], {𝕒}^{*}=\left[\begin{array}{c}{𝕒}_1^{*}\\ {𝕒}_2^{*}\\ 0\end{array}\right], {𝕓}^{*}=\left[\begin{array}{c}{𝕓}_1^{*}\\ {𝕓}_2^{*}\\ 0\end{array}\right]$$

Now, we introduce the following definitions:

$$𝕢=\psi {𝕢}^j, 𝕡=\psi {𝕡}^j, p={\psi }_0 p^j, \frac{\partial p}{\partial n}={\psi }_0{\left(\frac{\partial p}{\partial n}\right)}^j$$

(30)

Substituting (30) into (29) yields

$$C^n{𝕢}^n={\displaystyle \sum }_{j=1}^{N_e}\left[-{{\displaystyle \int }}_{{\Gamma }_j}{𝕡}^{*}\psi d\Gamma \right]{𝕢}^j+{\displaystyle \sum }_{j=1}^{N_e}\left[{{\displaystyle \int }}_{{\Gamma }_j}{𝕢}^{*}\psi d\Gamma \right]{𝕡}^j+{\displaystyle \sum }_{j=1}^{N_e}\left[{{\displaystyle \int }}_{{\Gamma }_j}{𝕒}^{*}{\psi }_0 d\Gamma \right]p^j+{\displaystyle \sum }_{j=1}^{N_e}\left[{{\displaystyle \int }}_{{\Gamma }_j}{𝕓}^{*}{\psi }_0 d\Gamma \right]{\left(\frac{\partial p}{\partial n}\right)}^j$$

(31)

which may be expressed as

$$C^i{𝕢}^i=-{\displaystyle \sum }_{j=1}^{N_e}{\widehat{ℍ}}^{ij}{𝕢}^j+{\displaystyle \sum }_{j=1}^{N_e}{\widehat{\mathbb{G}}}^{ij}{𝕡}^j+{\displaystyle \sum }_{j=1}^{N_e}{\widehat{𝕒}}^{ij}{p}^j+{\displaystyle \sum }_{j=1}^{N_e}{\widehat{𝕓}}^{ij}{\left(\frac{\partial p}{\partial n}\right)}^j$$

(32)

in which

$${ℍ}^{ij}=\{\begin{array}{c}{\widehat{ℍ}}^{ij} if i\ne j\\ {\widehat{ℍ}}^{ij}+C^i if i=j\end{array}$$

(33)

Now, Equation (32) can be written as

$${\displaystyle \sum }_{j=1}^{N_e}{ℍ}^{ij}{𝕢}^j={\displaystyle \sum }_{j=1}^{N_e}{\widehat{\mathbb{G}}}^{ij}{𝕡}^j+{\displaystyle \sum }_{j=1}^{N_e}{\widehat{𝕒}}^{ij}{p}^j+{\displaystyle \sum }_{j=1}^{N_e}{\widehat{𝕓}}^{ij}{\left(\frac{\partial p}{\partial n}\right)}^j$$

(34)

which may be written as

$$ℍ\mathbb{Q}=\mathbb{G}\mathbb{P}+𝕒𝕚+𝕓𝕛$$

(35)

Substituting the boundary conditions into (35), the  following system can be established:

$$\mathbb{A} \mathbb{X}=\mathbb{B}$$

(36)

in which $\mathbb{X}$ and $\mathbb{B}$ are known matrices, and $\mathbb{A}$ is an unknown matrix.

Breuer et al. [30] implemented a robust and efficient partitioned semi-implicit predictor–corrector coupling algorithm with generalized modified shift-splitting (GMSS) to solve the resulting linear Equation (36) that arises from boundary element discretization [10,11], where poro-thermo-elastic coupling is considered rather than fluid–structure interaction coupling.

4. Numerical Results and Discussion

The principal purpose of this study was to propose an efficient hybrid BEM (HBEM) model to describe the two-dimensional nonlinear fractional biomechanical interactions in FGA biological tissues. The proposed HBEM technique, which is based on the coupling algorithm [13], should be applied to a wide variety of nonlinear fractional bio-thermomechanical problems. The three efficient iterative methods that were applied to solve the linear systems resulting from the proposed hybrid technique were the communication-avoiding Arnoldi (CA-Arnoldi) procedure [31], which was also employed by Fahmy [32,33]; regularization [34], which was also implemented by Fahmy [35], and generalized modified shift-splitting (GMSS) [36], whose performance was also demonstrated by Fahmy [37]. In the current study, the properties of isotropic, transversely isotropic, and anisotropic soft tissues as given in references [38,39] were considered with the limitations 0 ≤ x ≤ 7. Furthermore, in the current problem, 84 boundary elements and 404 internal points were used in the BEM discretization, as shown in Figure 2.

Figure 3 shows the distribution of bio-thermal stress along ${\sigma }_{11}$ the $x$-axis ($y=0.5$) in the functionally graded isotropic and anisotropic biological tissues for various fractional-order parameter $a$ values ($a=0.2,0.6,\mathrm{and}1.0$). It was noticed from this figure that the bio-thermal stress ${\sigma }_{11}$ increased to a maximum value in the range $0\le x\le 0.25$, then decreased, and then increased. This figure also shows that the bio-thermal stress ${\sigma }_{11}$ increased with a fractional parameter increase in anisotropic FGA soft tissues and decreased with a fractional parameter increase in isotropic FGA soft tissues.

Figure 4 shows the distribution of the bio-thermal stress ${\sigma }_{11}$ along the $x$-axis ($y=0.5$) in the functionally graded isotropic and anisotropic biological tissues for various graded parameter $m$ values ($m=0.2,0.6,\mathrm{and}1.0$). It can be seen from this figure that the bio-thermal stress ${\sigma }_{11}$ increased to a maximum value in the range $0\le x\le 0.05$, then decreased to a minimum value in the range $0.05\le x\le 0.9$, and moved as a wave propagation for the isotropic and anisotropic cases. This figure also shows that the bio-thermal stress ${\sigma }_{11}$ increased with a functionally graded parameter $m$ increase in isotropic and anisotropic functionally graded soft tissues.

Figure 5 displays the bio-thermal stress ${\sigma }_{11}$ variation along the $x$-axis ($y=0.5$) for the classical Fourier $\left(\mathrm{Fourier},\left({\tau }_q={\tau }_T=0\right)\right)$, single-phase-lag $\left(\mathrm{SPL},\left({\tau }_q=0\mathrm{and}{\tau }_T=25\right)\right)$, and dual-phase-lag $\left(\mathrm{DPL},\left({\tau }_q={\tau }_T=25\right)\right)$ models. It was noticed from this figure that the bio-thermal stress ${\sigma }_{11}$ began at positive values in the anisotropic case, but it began at negative values in the isotropic case. Furthermore, it increased to a maximum value in the range $0\le x\le 0.85$ in the isotropic case, then decreased in the range $0.85\le x\le 7$. This figure also shows that the maximum value of the bio-thermal stress ${\sigma }_{11}$ occurred in the classical Fourier model in the isotropic and anisotropic functionally graded soft tissues, and the minimum value of the bio-thermal stress ${\sigma }_{11}$ occurred in the dual-phase-lag model in the isotropic and anisotropic functionally graded soft tissues.

Figure 6 displays the bio-thermal stress ${\sigma }_{11}$ variation along the $x$-axis (y = 0.5) for the proposed hybrid technique (Present), finite difference method (FDM) [35], and finite element method (FEM) [36]. It can be seen from this figure that the present hybrid technique results were in excellent agreement with the FDM and FEM results.

Table 1. Processing times and numbers of iterations for the CA-Arnoldi, regularized, and GMSS methods.

Discretization Level	Preconditioning Level	CA-Arnoldi [31,32,33]		Regularized [34,35]		GMSS [36,37]
		Process Time	Number of Iterations	Process Time	Number of Iterations	Process Time	Number of Iterations
1 (32)	0	0.07	6	0.07	6	0.07	6
2 (56)	0	0.2	10	0.2	10	0.2	10
	1	0.16	8	0.16	8	0.16	8
3 (104)	0	0.5	13	0.64	14	0.4	12
	1	0.46	10	0.6	11	0.32	7
	2	0.4	7	0.56	9	0.28	5
4 (200)	0	2.48	14	2.84	18	1.64	14
	1	2.04	12	2.62	16	1.48	9
	2	1.68	8	1.96	11	1.24	7
	3	1.46	7	51.82	4	0.98	3
5 (392)	0	12.2	22	14.26	28	6.85	16
	1	10.06	20	12.36	26	5.79	14
	2	9.38	18	10.4	22	4.98	12
	3	8.28	14	9.48	16	4.04	10
	4	7.62	10	8.1	14	3.64	4
6 (776)	0	40.8	20	46.4	26	36.5	15
	1	38.5	18	42.2	24	32.4	13
	2	36.6	16	40.3	22	30.8	11
	3	32.5	14	36.8	16	26.2	9
	4	30.4	12	34.3	14	20.3	5
	5	28.2	8	32.9	12	18.2	3

shows the processing times and numbers of iterations for the CA-Arnoldi, regularized, and GMSS iterative methods at each discretization level, where the number of equations is written inside the parentheses. It can be shown from this table that the GMSS was more efficient than the CA-Arnoldi and regularized iterative methods.

Table 2. Comparison of computer resource requirements for the considered computations.

	FDM [40]	FEM [41]	Present
Number of nodes	62,000	60,000	64
Number of elements	14,800	12,400	32
CPU time (min)	250	240	4
Memory (MByte)	220	230	1
Disc space (MByte)	310	330	0
Accuracy of results (%)	2.4	2.2	1.2

shows a comparison of the computer resource requirements for modeling the nonlinear fractional bio-thermomechanical interactions in FGA soft tissues for the finite difference method (FDM), finite element method (FEM), and HBEM (Present). This table demonstrates the efficiency of our proposed hybrid technique.

The findings of this paper contribute to the development of mathematical models that can be applied in bio-thermal engineering applications, such as blood perfusion, temperature measurement during cryosurgery, analysis of microvascular heat transfer, thermal injury, cryotherapy, hyperthermia, and the influence of large blood vessels.

5. Conclusions

Some of the inferences that can be drawn from the current study are as follows:

• A new HBEM model was used to describe the nonlinear fractional biomechanical interactions in FGA biological tissues.
• The bioheat governing equation was solved by implementing the LRBFCM and GBEM for obtaining the temperature, and then the poroelastic governing equation was solved using the BEM to calculate the displacement at each time step.
• An efficient partitioned semi-implicit coupling algorithm was implemented with the GMSS to solve equations arising from the boundary element discretization.
• The numerical findings were depicted graphically to display the influences of the graded parameter, fractional parameter, and anisotropic property on the bio-thermal stress.
• The numerical findings also show the differences between the Fourier, single-phase-lag, and dual-phase-lag bioheat models, and verified the validity, accuracy, and effectiveness of the developed HBEM.
• The main advantages of the current HBEM model are its generality and simplicity.
• The numerical findings supported the claim that the proposed method offers more advantages than other domain discretization techniques.

The findings of this paper contribute to the development of mathematical models that can be applied in bio-thermal engineering applications, such as blood perfusion, temperature measurement during cryosurgery, analysis of microvascular heat transfer, thermal injury, cryotherapy, hyperthermia, and the influence of large blood vessels.