Three-Dimensional Boundary Element Strategy for Stress Sensitivity of Fractional-Order Thermo-Elastoplastic Ultrasonic Wave Propagation Problems of Anisotropic Fiber-Reinforced Polymer Composite Material

A new three-dimensional (3D) boundary element method (BEM) strategy was developed to solve fractional-order thermo-elastoplastic ultrasonic wave propagation problems based on the meshless moving least squares (MLS) method. The temperature problem domain was divided into a number of circular sub-domains. Each node was the center of the circular sub-domain surrounding it. The Laplace transform method was used to solve the temperature problem. A unit test function was used in the local weak-form formulation to generate the local boundary integral equations, and the inverse Laplace transformation method was used to find the transient temperature solutions. Then, the three-dimensional elastoplastic problems could be solved using the boundary element method (BEM). Initial stress and strain formulations are adopted, and their distributions are interpolated using boundary integral equations. The effects of the fractional-order parameter and anisotropy are investigated. The proposed method’s validity and performance are demonstrated for a two-dimensional problem with excellent agreement with other experimental and numerical results.


Introduction
All fiber-reinforced polymer (FRP) composite materials, which have significant potential for a wide range of infrastructure applications, contain thermosetting or thermoplastic resins as well as glass and/or carbon fibers. The load-bearing component of the composite is provided by the fiber network, while the resin aids in load transfer and fiber orientation. The resin regulates the manufacturing process and processing variables. Resins also protect the fabrics from environmental factors such as relative humidity-elevated temperatures and chemical attacks.
Significant research has been conducted on the development of FRP composite materials and their novel applications. Many efforts have yielded materials with improved structural properties. Because of their superior corrosion resistance, excellent thermomechanical properties, and high strength-to-weight ratio, FRP composite materials are being promoted as twenty-first-century materials. In terms of their embodied energy, FRP composite materials are also "greener" than traditional materials such as concrete and steel. The use of FRP composite materials in civil and military infrastructure can improve innovation, productivity, and performance while also providing longer service lives, resulting in lower life-cycle costs. These efforts demonstrate that the use of innovative composite materials and designs have significant potential to reduce infrastructure vulnerability.
The BEM with internal collocation nodes has been used to solve thermo-elastoplastic problems [1,2]. However, the BEM's advantage of ease of data preparation is lost in this scenario. Therefore, several BEM strategies have been proposed. Nowak and Neves [3] developed the multiple-reciprocity boundary element method, which cannot be used to analyze thermo-elastoplastic materials. The dual-reciprocity BEM was developed to solve thermoelastoplastic problems with an arbitrary heat source [4]. Eigenvalue analysis can be carried out using the real-part boundary element approach [5,6]. The local boundary element method was used by Sladek and Sladek [7] to solve elastoplastic problems without internal cells. For elastoplastic difficulties, Ochiai and Kobayashi [8] presented the triple-reciprocity BEM, which does not require internal cells. This method allows for a very accurate solution to be produced using only fundamental low-order solutions and reduces the requirements for data preparation. Ochiai [9] applied the triple-reciprocity BEM to solve 2D thermoelastoplastic problems with an arbitrary distributed heat source [10] and three-dimensional elastoplastic problems with initial strain formulas [10]. Recently, Fahmy et al. [11][12][13][14] developed fractional BEM schemes to solve certain thermoelastic problems.
In this paper, a new BEM strategy is developed to solve three-dimensional thermoelastoplastic wave propagation problems with an arbitrary distributed heat source. Boundary elements and arbitrary internal points are used in this strategy. For elastoplastic analysis, the initial strain or stress distribution is interpolated using boundary integral equations. Strong singularities in the calculation of stresses at internal sites become weak using this method. The impacts of anisotropy and the fractional-order parameter are examined. The validity and performance of the suggested method for a two-dimensional problem are demonstrated, showing excellent agreement with existing experimental and numerical results.

BEM Implementation for the Temperature Field
The heat conduction equation of a nonhomogeneous anisotropic fiber-reinforced polymer composite in the presence of the distributed heat source W [1]s (q) can be expressed as [15] ρ(x)c(x)D α t θ(x, t) = k ij (x)θ ,j (x, t) ,i + Q(x, t), where the parameters are defined in the Nomenclature Table at the end of this paper.
In the BEM formulation of 3D problems, the distributed heat source function W S 1 (q) is interpolated using the following equations [16]: In 3D problems, the polyharmonic function with the volume distribution T [ f ]A (p, q) is introduced to achieve smooth interpolation and can be described as [17] where r denotes the distance between observation point p and loading point q. On the basis of Caputo's finite difference technique, at ( f + 1)∆τ and f ∆τ, the following formula can be established [18]: where By employing Equation (6), the fractional nonlinear heat conduction Equation (1) is transformed into the following equation [19]: Let Ω be the analyzed domain of the considered problem and the initial condition be The MLS approximates Thus, the following definitions can be deduced: Now, by implementing the Laplace transformation to Equation (1), the following equation is obtained: where Q(x, s) = 1−R x 0 e xa x 0 J(s) , and J(s) = J 0 (s+τ 1 ) 2 , s > τ 1 . The local weak form of Equation (11) can be described as (13) in which θ * (x) and ∂Ω a s are the weight function and local sub-domain boundary, respectively. Applying the Gauss theorem to Equation (13) yields where q(x, s) = k lj (x)θ ,j (x, s)n l (x). (15) and Based on the fundamental solution of (8), the local weak form (14) yields the following boundary integral representation: The MLS is employed to compute the heat flux as On the basis of [20], Equation (17) can be re-expressed as Considering the following representations The inverse Laplace transform [21] has now been implemented to obtain the physical quantities in time domain.

BEM Implementation for the Elastoplastic Field
Now, our purpose is to solve the following boundary integral equation [1,2]: [1] jki (P, q) . ε [1] I jk (q)dΩ where c ij , . ε [1] I jk (q), . u j (Q), and . p j (Q) are the free coefficient, initial strain rate, displacement rate, and surface traction rate, respectively. However, r, Γ, and Ω are the distance between the observation point and loading point, the boundary, and domain, respectively.
ijkl (p, q) is calculated using Equation (51) and the displacement-stress relationship as ijkl (p, q)/∂n and ε [3]A ijkl (p, q) are also obtained as ∂r ∂n The first thermal load is T S , the final thermal load is T 0 , and the number of iterations is N. Then, the incremental load is (T 0 − T S )/N.
The following iterative relationship is used to solve the current thermo-elastoplastic problem: where σ k 0 , σ k+1 0 , H, and dε P e are yield stress at k, yield stress at k + 1, strain hardening, and equivalent plastic strain increment, respectively. Based on the von Mises yield criterion, the stresses rate in Equation (72) yields the deviatoric stress tensor S ij , and the equivalent stress σ e can be computed as where The following Prandtl-Reuss equation is employed to calculate the plastic strain increment dε p ij as dε

Numerical Results and Discussion
The proposed BEM method is general because it can be used to deal with a wide range of fractional thermo-elastoplastic problems affecting anisotropic fiber-reinforced polymer composite materials. Additionally, it is simple because only the surface of the domain needs to be discretized.
The pure anisotropic fiber-reinforced behavior satisfies where a ≡ (a 1 , a 2 , a 3 ), a 2 1 + a 2 2 + a 2 3 Additionally, the isotropic behavior satisfies α = β = (µ L − µ T ) = 0. As illustrated in Figure 1, the domain of the considered 3D problem includes 40 boundary nodes and 81 internal nodes. Additionally, we assumed that the wave direction is parallel to the x 1 -axis.
where ≡ (a , a , a ), a + a + a Additionally, the isotropic behavior satisfies = ̅ = ( ̅ − ̅ ) = 0. As illustrated in Figure 1, the domain of the considered 3D problem includes 40 boundary nodes and 81 internal nodes. Additionally, we assumed that the wave direction is parallel to the -axis..                                     Figure 8 explains the distribution of the strain ε 11 sensitivity along the x 1 -axis, which, in isotropic and anisotropic cases, begins with a negative value. It can be seen from this figure that the distribution of the strain ε 11 sensitivity initially increases and then decreases along the x 1 -axis. Additionally, it has α = 0.7 > α = 0.4 > α = 1.0 > α = 0.1 in anisotropic cases but α = 0.4 > α = 0.7 > α = 1.0 > α = 0.1 for isotropic cases. This figure demonstrates that the fractional-order parameter has a significant effect on the strain ε 11 sensitivity in both isotropic and anisotropic cases. The strain ε 11 sensitivity curves at the upper ( α = 1.0) and lower ( α = 0.1) values of the fractional parameter are also close to each other, and we notice that they are closer in the isotropic case than in the anisotropic case. It is demonstrated that the strain ε 11 sensitivity curves at the interface values diverge from each other, as they are further away in the isotropic case than in the anisotropic case.
Polymers 2022, 14, x FOR PEER REVIEW 18 of 25 in anisotropic cases but α = 0.4 > α = 0.7 > α = 1.0 > α = 0.1 for isotropic cases. This figure demonstrates that the fractional-order parameter has a significant effect on the strain sensitivity in both isotropic and anisotropic cases. The strain sensitivity curves at the upper (α = 1.0) and lower (α = 0.1) values of the fractional parameter are also close to each other, and we notice that they are closer in the isotropic case than in the anisotropic case. It is demonstrated that the strain sensitivity curves at the interface values diverge from each other, as they are further away in the isotropic case than in the anisotropic case.    Figure 9 illustrates the distribution of the strain ε 12 sensitivity along the x 1 -axis in the context of isotropic and anisotropic fiber-reinforced polymer composites for various fractional-order values. It can be noticed from this figure that the strain ε 12 sensitivity increases as x 1 increases at small x 1 values. Additionally, it has α = 0.4 > α = 1.0 > α = 0.1 > α = 0.7 in anisotropic cases, but it has α = 0.7 > α = 0.4 > α = 1.0 > α = 0.1 in isotropic cases, which are close to the approximate values as x 1 tends to infinity. This figure demonstrates that the fractional-order parameter has an important effect on the strain ε 12 sensitivity in both isotropic and anisotropic cases. The strain ε 12 sensitivity curves at the upper ( α = 1.0) and lower ( α = 0.1) values of the fractional parameter are congruent in both cases. It is demonstrated that the strain ε 12 sensitivity curves at the interface values diverge from each other, as they are further away in the anisotropic case than in the isotropic case.   Figure 10 explains the distribution of the strain ε 22 sensitivity along the x 1 -axis, which starts near zero at x 1 = 0 in the context of both isotropic and anisotropic cases. It is noticed that distribution of the strain ε 22 sensitivity first decreases then increases as x 1 increases at small x 1 values. Additionally, it has α = 0.7 > α = 0.1 > α = 1.0 > α = 0.4 in isotropic and anisotropic cases.   This figure demonstrates that the fractional-order parameter has a significant effect on the strain sensitivity in both isotropic and anisotropic cases. The strain sensitivity curves at the upper (α = 1.0) and lower (α = 0.1) values of the fractional parameter are also close to each other, and we notice that they are closer in the anisotropic case than in the isotropic case. It is demonstrated that the strain sensitivity curves at the interface values diverge from each other, as they are further away in the anisotropic case than in the isotropic case. Figure 11 depicts the distribution of the strain sensitivity along the -axis, which starts from zero at = 0 in the context of isotropic and anisotropic cases. It noticed that the strain sensitivity is increases first and decreases and then increases again was increases. Additionally, it has α = 0.1 > α = 1.0 > α = 0.4 > α = 0.7 for This figure demonstrates that the fractional-order parameter has a significant effect on the strain ε 22 sensitivity in both isotropic and anisotropic cases. The strain ε 22 sensitivity curves at the upper ( α = 1.0) and lower ( α = 0.1) values of the fractional parameter are also close to each other, and we notice that they are closer in the anisotropic case than in the isotropic case. It is demonstrated that the strain ε 22 sensitivity curves at the interface values diverge from each other, as they are further away in the anisotropic case than in the isotropic case. Figure 11 depicts the distribution of the strain ε 13 sensitivity along the x 1 -axis, which starts from zero at x 1 = 0 in the context of isotropic and anisotropic cases. It noticed that the strain ε 13 sensitivity is increases first and decreases and then increases again was x 1 increases. Additionally, it has α = 0.1 > α = 1.0 > α = 0.4 > α = 0.7 for isotropic cases and α = 0.7 > α = 0.1 > α = 1.0 > α = 0.4 for anisotropic cases. This figure demonstrates that the fractional-order parameter has a significant effect on the strain ε 13 sensitivity in both isotropic and anisotropic cases. The strain ε 13 sensitivity curves at the upper ( α = 1.0) and lower ( α = 0.1) values of the fractional parameter are also close to each other, and we notice that they are closer in the anisotropic case than in the isotropic case. It is demonstrated that the strain ε 13 sensitivity curves at the interface values diverge from each other, as they are further away in the anisotropic case than in the isotropic case. sensitivity in both isotropic and anisotropic cases. The strain sensitivity curves at the upper (α = 1.0) and lower (α = 0.1) values of the fractional parameter are also close to each other, and we notice that they are closer in the anisotropic case than in the isotropic case. It is demonstrated that the strain sensitivity curves at the interface values diverge from each other, as they are further away in the anisotropic case than in the isotropic case. Figure 11. Distribution of the sensitivity along -axis in isotropic and anisotropic FRP composites for various fractional-order values. Figure 12 explains the distribution of the strain sensitivity along the -axis, which starts near zero at = 0 in the context of isotropic and anisotropic fiber-reinforced polymer composites for various fractional-order values. It can be seen from this figure that the distribution of strain sensitivity initially increases and then decreasing along the -axis. Additionally, it has α = 0.7 > α = 0.4 > α = 1.0 > α = 0.1 in isotropic cases but α = 0.4 > α = 1.0 > α = 0.1 > α = 0.7 in anisotropic cases. This figure demonstrates that the fractional-order parameter has a significant effect on the strain sensitivity in both isotropic and anisotropic cases. The strain sensitivity curves at the upper (α = 1.0) and lower (α = 0.1) values of the fractional parameter are also close to each other, and we notice that they are closer in the anisotropic case than in the isotropic case. It is demonstrated that the strain sensitivity curves at the interface values diverge from each other, as they are further away in the anisotropic case than in the isotropic case. Figure 11. Distribution of the ε 13 sensitivity along x 1 -axis in isotropic and anisotropic FRP composites for various fractional-order values. Figure 12 explains the distribution of the strain ε 23 sensitivity along the x 1 -axis, which starts near zero at x 1 = 0 in the context of isotropic and anisotropic fiber-reinforced polymer composites for various fractional-order values. It can be seen from this figure that the distribution of strain ε 23 sensitivity initially increases and then decreasing along the x 1 -axis. Additionally, it has α = 0.7 > α = 0.4 > α = 1.0 > α = 0.1 in isotropic cases but α = 0.4 > α = 1.0 > α = 0.1 > α = 0.7 in anisotropic cases. This figure demonstrates that the fractional-order parameter has a significant effect on the strain ε 23 sensitivity in both isotropic and anisotropic cases. The strain ε 23 sensitivity curves at the upper ( α = 1.0) and lower ( α = 0.1) values of the fractional parameter are also close to each other, and we notice that they are closer in the anisotropic case than in the isotropic case. It is demonstrated that the strain ε 23 sensitivity curves at the interface values diverge from each other, as they are further away in the anisotropic case than in the isotropic case.    Figure 13 depicts the distribution of strain ε 33 , which starts from zero at x 1 = 0 in the context of isotropic and anisotropic cases. It noticed that the distribution decreases and then increases as x 1 increases at small x 1 values. Additionally, it has α = 0.1 > α = 1.0 > α = 0.4 > α = 0.7 in both isotropic and anisotropic cases. This figure demonstrates that the fractional-order parameter has a significant effect on the strain ε 23 sensitivity in both isotropic and anisotropic cases. The strain ε 23 sensitivity curves at the upper ( α = 1.0) and lower ( α = 0.1) values of the fractional parameter are also close to each other, and we notice that they are closer in the anisotropic case than in the isotropic case. It is demonstrated that the strain ε 23 sensitivity curves at the interface values diverge from each other, as they are further away in the anisotropic case than in the isotropic case.  Figure 13 depicts the distribution of strain , which starts from zero at = 0 in the context of isotropic and anisotropic cases. It noticed that the distribution decreases and then increases as increases at small values. Additionally, it has α = 0.1 > α = 1.0 > α = 0.4 > α = 0.7 in both isotropic and anisotropic cases. This figure demonstrates that the fractional-order parameter has a significant effect on the strain sensitivity in both isotropic and anisotropic cases. The strain sensitivity curves at the upper (α = 1.0) and lower (α = 0.1) values of the fractional parameter are also close to each other, and we notice that they are closer in the anisotropic case than in the isotropic case. It is demonstrated that the strain sensitivity curves at the interface values diverge from each other, as they are further away in the anisotropic case than in the isotropic case. There are no published results that demonstrate the validity and accuracy of the current BEM method strategy. On the other hand, some studies can be thought of as special cases in the context of this current general study. The special case distributions , , and for the considered BEM combined the finite element method/normal mode method (FEM-NMM) of An et al. [23] and the experimental technique (Exp.) of Solodov et al. [24] and are shown in Figures 14-16 for fractional-order (α = 0.4) anisotropic fiber- There are no published results that demonstrate the validity and accuracy of the current BEM method strategy. On the other hand, some studies can be thought of as special cases in the context of this current general study. The special case distributions σ 11 , σ 12 , and σ 22 for the considered BEM combined the finite element method/normal mode method (FEM-NMM) of An et al. [23] and the experimental technique (Exp.) of Solodov et al. [24] and are shown in Figures 14-16 for fractional-order ( α = 0.4) anisotropic fiber-reinforced polymer composites. These results show that the BEM findings are in excellent agreement with those of FEM-NMM [23] and Exp. [24]. As a result, the validity of the proposed technique was confirmed.
Polymers 2022, 14, x FOR PEER REVIEW 22 of 25 reinforced polymer composites. These results show that the BEM findings are in excellent agreement with those of FEM-NMM [23] and Exp. [24]. As a result, the validity of the proposed technique was confirmed.

Conclusions
The following findings can be drawn from the present paper: 1. Advanced BEM was applied to solve fractional-order thermo-elastoplastic ultrasonic wave propagation problems affecting anisotropic fiber-reinforced polymer composite materials 2. The Laplace transform was used to eliminate the time variable from the governing equations. 3. The unit step test function was used to derive the local boundary integral equations. 4. The MLS scheme was developed to treat the domain integrals and approximate physical quantities. 5. The numerical data demonstrate the current MLS approach's accuracy, feasibility, effectiveness, and convergence. 6. The inverse Laplace transformation method was then used to find the transient temperature solutions. 7. The current technique's main advantage is its generality and simplicity. 8. The initial stress and strain distributions are interpolated using boundary integral equations. 9. Numerical results show that the fractional-order parameter and anisotropy have significant effects on the thermoelastic behavior of fiber-reinforced polymer composites. 10. The numerical results show that the proposed strategy outperforms previous experimental and numerical methods.

Conclusions
The following findings can be drawn from the present paper: 1.
Advanced BEM was applied to solve fractional-order thermo-elastoplastic ultrasonic wave propagation problems affecting anisotropic fiber-reinforced polymer composite materials 2.
The Laplace transform was used to eliminate the time variable from the governing equations.

3.
The unit step test function was used to derive the local boundary integral equations. 4.
The MLS scheme was developed to treat the domain integrals and approximate physical quantities.

5.
The numerical data demonstrate the current MLS approach's accuracy, feasibility, effectiveness, and convergence. 6.
The inverse Laplace transformation method was then used to find the transient temperature solutions. 7.
The current technique's main advantage is its generality and simplicity. 8.
The initial stress and strain distributions are interpolated using boundary integral equations. 9.
Numerical results show that the fractional-order parameter and anisotropy have significant effects on the thermoelastic behavior of fiber-reinforced polymer composites. 10. The numerical results show that the proposed strategy outperforms previous experimental and numerical methods. 11. The findings presented in this paper may be of interest to researchers in material science, mathematical physics, and geothermal engineering as well as those working on the development of anisotropic fiber-reinforced polymer composite materials.