## 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 [PERSON] [32; 33]; regularization [34], which was also implemented by [PERSON] [35], and generalized modified shift-splitting (GMSS) [36], whose performance was also demonstrated by [PERSON] [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\leq x\leq 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,\) 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\leq x\leq 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 2: Boundary element model of the considered problem. 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,\text{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\leq x\leq 0.05\), then decreased to a minimum value in the range \(0.05\leq x\leq 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(\text{Fourier},\left(\tau_{q}=\tau_{T}=0\right)\right)\), single-phase-lag \(\left(\text{SPL},\left(\tau_{q}=0\text{and}\tau_{T}=25\right)\right)\), and dual-phase-lag \(\left(\text{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\leq x\leq 0.85\) in the isotropic case, then decreased in the range \(0.85\leq x\leq 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 4: Distribution of \(\sigma_{11}\) along the \(x\)-axis for various graded parameter values. Figure 3: Distribution of \(\sigma_{11}\) along the \(x\)-axis for various fractional-order parameter values.