3D Boundary Element Model for Ultrasonic Wave Propagation Fractional Order Boundary Value Problems of Functionally Graded Anisotropic Fiber-Reinforced Plates

: This paper proposes a three–dimensional (3D) local boundary element model based on meshless moving least squares (MLS) method for ultrasonic wave propagation fractional order boundary value problems of functionally graded anisotropic (FGA) ﬁber-reinforced plates. The problem domain is split into several circular sub-domains. The nodal points are randomly distributed across the examined region. Each node is the focal point of a circular sub-domain that encircles it. The Laplace-transform approach is used to solve dynamic issues. In the local weak form of the governing equations for the converted quantities, a unit test function is utilized. The Gauss divergence theorem to the weak-form is used to produce local boundary-domain integral equations. A meshless approximation is achieved using the MLS method. To ﬁnd time-dependent solutions, an inverse Laplace-transform approach is used. The effects of the fractional order parameter, functionally graded material, anisotropy, and the time characteristic of the laser pulse are investigated. The proposed method’s validity and performance are demonstrated for a two-dimensional problem with excellent agreement with the ﬁnite element method.


Introduction
Recently, fractional calculus has gained popularity as a method for studying the theory and applications of arbitrary non-integer order derivatives and integrals. This mathematical branch has recently emerged as a useful and powerful tool for mathematical modeling in a variety of engineering, industrial, and materials-science applications [1]. Fractional-order operators are useful in expressing the memory and heredity properties of many materials and processes due to their nonlocal nature. According to the associated literature published by prominent fractional calculus journals [2], the primary focus of the investigation had shifted from traditional integer-order models to fractional-order models.
Fractional calculus is used in many fields, including hereditary solid mechanics, fluid dynamics, viscoelasticity, heat conduction modeling and identification, biology, food engineering, econophysics, biophysics, biochemistry, robotics and control theory, signal and image processing, electronics, electric circuits, wave propagation, nanotechnology, and many others [3].
Many mathematicians have contributed to the history of fractional calculus, which begins in 1730 with Euler's first mention of interpolating between integral orders of a derivative. In 1812, Laplace used an integral to define a fractional derivative.
In 1819, Lacroix introduced the first fractional order derivative into calculus, expressing the fractional order derivative law as follows where u is any arbitrary positive or negative order. The first use of fractional operations has been established by Abel in 1823, where he used the following formula Liouville introduced the foundations of fractional calculus in 1832, where he supposed that d v dx v (e ax ) = a v e ax for v > 0 to derive the fractional derivative as In 1848, Hargreave introduced the first the generalization of Leibniz's n th derivative of a product, as follows where D (n) , D (v−n) are ordinary differentiation and fractional operation, respectively. In 1876, Riemann introduced his fractional integration theory based on Taylor series as (6) in which ψ(x) is the complementary function. Laurent defined arbitrary order integration in 1884, based on the Cauchy's integral law for complex-valued analytic functions as where c D v x represents differentiation of order v of the function f along the x-axis. Cauchy's fractional order derivative is denoted as [4] f (α) Caputo proposed the following fractional order derivative in 1967 Lancaster and Salkauskas [5] introduce the MLS method for surface construction, and the corresponding error analysis is discussed in [6]. The MLS method provides a best approximation in a weighted least squares sense, and it emphasizes the compacted support of the weight function especially, so it has local characteristics. It does, however, have some limitations, such as complex computation and the absence of the Kronecker delta function property. The classical theory is written in terms of out-of-plane displacement and its derivatives, or, as demonstrated in [7], plate rotations can be represented by an irrotational field. The integral equation has four boundary parameters and satisfying two boundary conditions is required to obtain a single solution. When considering a polygonal plate, it should be noted that, in addition to the two boundary unknowns, a concentrated reaction is placed at each corner as an additional parameter in the boundary value problem. The classical theory's inaccuracy turns out to be of practical interest in the edge zone of a plate and around holes with a diameter no larger than the plate's thickness. To overcome the above-mentioned characteristics of the classical theory's one-displacement dependence. Mindlin [8] and Reissner [9] proposed similar theories based on shear deformation. Unless those theories are designed to deal with thin or thick plates, they are commonly referred to as thick-plate theories in numerical analysis. In the constitutive equations, curvatures are not directly related to out-of-plane displacement derivatives, and three boundary conditions should be satisfied in the boundary value problem rather than the two in classical theory.
The boundary element method (BEM) has emerged as a viable numerical solution for plate problems. Wang and Huang [10] were the first to use BEM to model orthotropic thick plates. Meshless approaches to continuum mechanics problems have received a lot of attention in the last decade [11]. For such plates, meshless approaches with continuous stress approximation are more convenient [12]. The element-free Galerkin approach was used by Krysl and Belytschko [13] to offer the first use of a meshless method to plate problems. Their results showed excellent convergence, however, their formulation is not applicable to shear deformable plate problems. Fahmy [14] applied BEM to threetemperature distributions in carbon nanotube fiber-reinforced plates with inclusions.
Finding an analytical solution to a problem is extremely difficult in general; thus, several engineering papers devoted to numerical methods have studied such problems in various thermoelasticity topics, such as thermoelastic metal and alloy discs [15,16], generalized magneto-thermoelasticity [17], and micropolar magneto-thermoviscoelasticity [18]. However, several papers have used the boundary element method in general, for example, to solve micropolar FGA composites problems [19], Photothermal waves [20], and magneto-thermo-viscoelasticity [21]. Because the trial and test functions can be chosen from different functional spaces, the meshless local Petrov-Galerkin (MLPG) method [22] is a fundamental foundation for the derivation of many meshless formulations. The method has also been applied successfully to plate problems [23][24][25].
In the present paper, the local BEM based on MLS method has been successfully applied to solve dynamic problems of FGA fiber-reinforced plates. The Laplace-transform technique is used to solve the governing equations system of elastodynamic Reissner bending theory. The local boundary-domain integral equations are derived by applying the Gauss divergence theorem to the local weak-form governing equations. For each time instant under consideration, boundary value problems must be treated using a variety of Laplace-transform parameter values. The transformed quantities can be calculated in time domain by using the numerical inverse Laplace transformation method.

Formulation of the Problem
The plane-stress constitutive equation for nonhomogeneous anisotropic plates can be written as [26]  The heat conduction equation of nonhomogeneous anisotropic plates can be expressed as [20]: The heat conduction equation of nonhomogeneous anisotropic plates can be expressed as [20]: in which a(Q, To investigate the pure anisotropic fiber-reinforced effect, we assumed that where the reinforcing parameters t, x ̅ and (r̅ w − r̅ s ) introduce significantly anisotropic behavior in the considered structure. Moreover, the isotropic behavior can be achieved under the following condition t = x ̅ = (r̅ w − r̅ s ) = 0 (see Table A1).

BEM Implementation for the Temperature Field
Based on Caputo's finite difference scheme at ( + 1)∆E and ∆E, we can write [7] where By using Equation (15), the fractional nonlinear heat conduction Equation (13) becomes [23] Let the analyzed domain of the studied plate be denoted by Ω with ‡ C and ‡ for the top and bottom surfaces, respectively. It is assumed that the initial condition is We considered the sub-domain Ω $ , the MLS approximates ‰ (Q) of , as ‰ (Q) = Š s (Q)‹(Q) ∀Q ∈ Ω $ , where Š s (Q) = Ž (Q), (Q), … , (Q)• is a vector of complete monomial basis of order , and ‹(Q) is a vector of coefficients * Y (Q), ({ = 1, 2 , … , ), Q = Ž , , L • s . In the considered 3D analysis based on the MLPG method, we have By applying the Laplace transform to Equation (13), we obtain The heat conduction equation of nonhomogeneous anisotropic plates can be expressed as [20]: in which a(Q, To investigate the pure anisotropic fiber-reinforced effect, we assumed that where the reinforcing parameters t, x ̅ and (r̅ w − r̅ s ) introduce significantly anisotropic behavior in the considered structure. Moreover, the isotropic behavior can be achieved under the following condition t = x ̅ = (r̅ w − r̅ s ) = 0 (see Table A1).

BEM Implementation for the Temperature Field
Based on Caputo's finite difference scheme at ( + 1)∆E and ∆E, we can write [7] where By using Equation (15), the fractional nonlinear heat conduction Equation (13) becomes [23] Let the analyzed domain of the studied plate be denoted by Ω with ‡ C and ‡ for the top and bottom surfaces, respectively. It is assumed that the initial condition is We considered the sub-domain Ω $ , the MLS approximates ‰ (Q) of , as ‰ (Q) = Š s (Q)‹(Q) ∀Q ∈ Ω $ , where Š s (Q) = Ž (Q), (Q), … , (Q)• is a vector of complete monomial basis of order , and ‹(Q) is a vector of coefficients * Y (Q), ({ = 1, 2 , … , ), Q = Ž , , L • s . In the considered 3D analysis based on the MLPG method, we have By applying the Laplace transform to Equation (13), we obtain where the reinforcing parameters α, β and (µ L − µ T ) introduce significantly anisotropic behavior in the considered structure. Moreover, the isotropic behavior can be achieved under the following condition α = β = (µ L − µ T ) = 0 (see Table A1).

BEM Implementation for the Temperature Field
Based on Caputo's finite difference scheme at ( f + 1)∆τ and f ∆τ, we can write [7] where By using Equation (15), the fractional nonlinear heat conduction Equation (13) becomes [23] Let the analyzed domain of the studied plate be denoted by Ω with S + and S − for the top and bottom surfaces, respectively. It is assumed that the initial condition is We considered the sub-domain Ω x , the MLS approximates is a vector of complete monomial basis of order m, and a(x) is a vector of coefficients a j (x), (j = 1, 2 , . . . , m), In the considered 3D analysis based on the MLPG method, we have By applying the Laplace transform to Equation (13), we obtain in which where Q(x, s) The local weak form of the governing Equation (20) for x a ∈ Ω a s can be written as in which θ * (x) is a weight function.
Applying the Gauss divergence theorem to Equation (22) we obtain where ∂Ω a s denotes the local sub-domain boundary and If a Heaviside unit step function is used as the test function, then we can write θ * (x) in each subdomain as follows Now, by using the fundamental solution of (17), the following local boundary integral equation is derived from the local weak form (23) The boundary integral formulation will be obtained if zero body heat sources are assumed. The MLS is used to calculate the heat flux q(x, s) as follows According to [26], we can write (26) as follows Now, we consider the following notations By implementing numerical inverse Laplace transformation method [27], the transformed quantities can be calculated in time domain.

BEM Implementation for the Displacement Field
Suppose that the material parameters are graded throughout the functionally graded fiber-reinforced plate thickness as where P, P t , P b and n are the generic property, top face property, bottom face property and functionally graded parameter, respectively. The bending moments M αβ as well as the shear forces Q α are identified as where κ = 5/6 in the Reissner plate theory. Substituting Equation (10) into (31), we can write where In which the material parameters D αβ and C αβ can be written as For a general variation in material properties as a function of plate thickness According to Reissner [28], the equations of motion can be expressed as Now, the Laplace-transform can be defined as By applying the Laplace-transform (38) to (37), we obtain where s denotes the Laplace transform parameter and where w α (x) and w 3 (x) are the displacement initial values and . w α (x) and . w 3 (x) are the displacement initial velocities.
MLPG techniques generate the weak-form circular local sub-domains such as Ω s , which is a small area assigned to each node within the global domain Ω. For x i ∈ Ω i s , the local weak form of the governing Equations (39) and (40) are as follows: where w * αβ (x) and w * (x) are test functions. According to [26], the application of Gauss divergence theorem to Equations (43) and (44) yields where ∂Ω i s is the boundary of the local sub-domain and where the following unit step functions are chosen as the test functions in each sub-domain Fractal Fract. 2022, 6, 247 8 of 23 The local weak forms (45) and (46) are then transformed into the local boundary integral equations based on the unit step functions w * αβ (x) and w * (x) of each sub-domain as follows where the MLS approximations w α (x, s) for rotations and w 3 (x, s) for deflections [26]. Thus, the generalized displacements can be written as Substituting from (51) into (31), we obtain where and N α (x) which are connected to n(x) on ∂Ω s can be expressed as Moreover, B a α can be expressed in terms of the shape functions as (54) Now, we can write [26] Q(x, s) = C(x) where Then, insertion of the MLS-discretized force fields (52) and (55) into the local boundary integral Equations (49) and (50) yields the discretized local integral equations (LIEs) Fractal Fract. 2022, 6, 247 in which According to [26], and using (51) with Equations (57) and (58), we can write where x i is the source point located on the global boundary Γ, ∂Ω i s is the sub-domain boundary, Γ i sM is the global boundary part with prescribed bending moment, Γ i sw is the global boundary part with prescribed rotations or displacements, and By implementing numerical inverse Laplace transformation method [27], we obtain In our numerical analyses, we considered N = 10 and s = i ln 2/t (i = 1, 2 , . . . , N).

Numerical Results and Discussion
We consider a special case of this study for comparison purposes, to demonstrate the accuracy, feasibility, effectiveness, and convergence of the present MLS method, so we define the root mean square error Rew in terms of exact solution w(x i , t), numerical solution w h (x i , t) and M nodes surrounding x as follows [29] In which ξ j (x) can be expressed as the linear combination of the MLS shape functions where v is a constant that can have various values in [0, 1]. Now, we give the results of Rew with different number of collocation points (40, 80, 160) in Table 1. As shown in Table 1, as the collocation increases, the value of Rew decreases, and the results and error analysis agree well. We find that the approximation effect of the present MLS method suggested in this study is the best in all cases when compared to moving least squares and local radial basis functions (MLS-LRBF) [29] and Modified Moving Least Squares (MMLS) [30]. In our numerical calculations, we considered an anisotropic FGM clamped plate under a uniform impact load with a side-length a = 0.254 m, plate thicknesses h/a = 0.05, and a Heaviside time dependence.
The material properties and geometry parameters are as follows: mass density ρ = 7.166 × 10 3 kg m −3 . A quadratic variation of the volume fraction V is considered for the considered plate. For the approximation of rotations and deflections in our numerical calculations, 441 nodes with a regular distribution were used [26]. If s is the distance between two nodes, the radius of the circular subdomain is chosen as r loc = 0.4s and the radius of the support domain for node a is r a = 4r loc .
The following material parameters for anisotropic FGA fiber-reinforced plate are used in numerical analysis: where reinforcement parameters α, β and µ L − µ T introduce anisotropic behavior in the considered functionally graded fiber-reinforced plates.
The domain boundary of the considered problem has been discretized into 42 boundary elements and 68 internal points, as illustrated in Figure 1.  The following material parameters for orthotropic FGM plate are used in numerical analysis: Young's moduli E 2t =0.6895 · 10 10 N m 2 , E 1t = 2E 2t , Poisson's ratios ν 21 = 0.15, ν 12 = 0.3, and shear modulus are A quadratic variation of the volume fraction V is considered with Young's moduli on the bottom side are: The following material parameters for isotropic FGM plate are used in numerical analysis: Young's modulus E 1 =E 2 = 0.6895 · 10 10 N/m 2 , Poisson's ratios ν 21 = ν 12 = 0.3, the thermal expansion coefficients α 11 = α 22 = 1 · 10 −5 deg −1 , and shear modulus are . We considered numerical results of three-dimensional problem in the computational domain that consists of 40 boundary nodes and 81 internal nodes as shown in Figure 1.   As shown in Figure 2, the value of the thermal stress wave always increas a positive value. It grows until it reaches its maximum value in the 0 < < 1 Since then, it has been on a downward trend. Finally, it tends to zero as it move direction of wave propagation. In both FG and H anisotropic fiber-reinforced pla thermal stress wave behaves similarly. As shown in Figure 2, the H reduces t imum value of the thermal stress wave . When the fractional order parame varied, the distributions of thermal stress wave are similar. The thermal stress wave , as shown in Figure 3, has a negative value at and then has a downward trend. It moves in the form of wave propagation after a of rising. The thermal stress wave behaves similarly in the FG anisotropic fib forced plate as it does in the H anisotropic fiber-reinforced plate. As shown in F the FG fiber-reinforced anisotropic plate increases the minimum value of the  The thermal stress wave exhibits the same behavior in FG and H anisotropic fiber-reinforced plates, as shown in Figure 4. It demonstrates that the value of the thermal stress wave reaches a negative value early on and has a downward trend. Since then, it has risen from the lowest to the highest point. Finally, it moves in the direction of wave propagation. Figure 4 shows that the maximum value of the thermal stress wave in FG anisotropic fiber-reinforced plates is greater than that in H anisotropic fiber-reinforced plates. The distributions of the thermal stress wave are similar when the fractional anisotropic fiber-reinforced plate, the value of the thermal stress wave in the range of 0 ≤ ≤ 2.2 is the same as in the FG anisotropic fiber-reinforced plate of Figure 2 From a positive value to its maximum value, it always increases. Then it tends to zero and moves along with the wave. The FG anisotropic fiber-reinforced plate, as shown in Figure  5, increases the amplitude of the thermal stress wave . The distributions of the therma stress wave are similar when the fractional order parameter is varied.  The thermal stress wave has a negative value at = 0, as shown in Figure 6. The behavior of the thermal stress wave in the FG anisotropic fiber-reinforced plate is identical to that shown in Figure 3. In the absence of gravity, it has a downward trend in the range of 0 ≤ ≤ 1.7. It gradually decreases after a period of rising. As shown in Figure 6, the FG anisotropic fiber-reinforced plate increases the amplitude of the thermal stress wave . When the fractional order parameter is changed, the distributions of the thermal stress wave are similar. As shown in Figure 7, the behavior of the thermal stress wave in the FG anisotropic fiber-reinforced plate is identical to that of the FG anisotropic fiber-reinforced plate in Figure 4. In the context of H anisotropic fiber-reinforced plate, it always decreases from a negative value at the start. It goes down to its minimum value in the range of 0 ≤ ≤ 2.9, then up from the minimum to the maximum. Finally, as the distance x increases, it tends to zero. As shown in Figure 7, the amplitude of the thermal stress wave is greater in the FG anisotropic fiber-reinforced plate than in the H anisotropic fiber-reinforced plate. The distributions of the thermal stress wave are similar when the fractional order parameter α is varied. The thermal stress wave has a negative value at = 0, as shown in Figure 6 The behavior of the thermal stress wave in the FG anisotropic fiber-reinforced plate is identical to that shown in Figure 3. In the absence of gravity, it has a downward trend in the range of 0 ≤ ≤ 1.7. It gradually decreases after a period of rising. As shown in Figure 6, the FG anisotropic fiber-reinforced plate increases the amplitude of the therma stress wave . When the fractional order parameter is changed, the distributions of the thermal stress wave are similar. As shown in Figure 7, the behavior of the thermal stress wave in the FG anisotropic fiber-reinforced plate is identical to that of the FG anisotropic fiber-reinforced plate in Figure 4. In the context of H anisotropic fiber-reinforced plate, it always decreases from a negative value at the start. It goes down to its minimum value in the range of 0 ≤ ≤ 2.9, then up from the minimum to the maximum. Finally, as the distance x increases, it tends to zero. As shown in Figure 7, the amplitude of the thermal stress wave is greater in the FG anisotropic fiber-reinforced plate than in the H anisotropic fiber-reinforced plate. The distributions of the thermal stress wave are similar when the fractional order parameter α is varied.  in the functionally graded (FG) and homogeneous (H) plates for anisotropic, orthotropic, and isotropic materials.
As shown in Figure 8, the thermal stress wave always rises from a positive starting point. It grows until it reaches its maximum value. It has remained on a downward trend since then. Finally, it moves along with the wave propagation and tends to zero The thermal stress wave behaves similarly in FG and H fiber-reinforced plates. As illustrated in Figure 8, the homogeneous case reduces the maximum value of the therma stress wave . The distributions of the thermal stress wave are similar when anisotropic, orthotropic, and isotropic fiber-reinforced plates are considered. As shown in Figure 9, the thermal stress wave decreases from zero at x=0 and coincides with the boundary condition. It decreases slightly and then rapidly in the range of 0 ≤ ≤ 1.5. Then it decreases and moves along with the wave. In both the FG and H fiber-reinforced plates, the distributions of the thermal stress wave are similar. The maximum value of the thermal stress wave in the H fiber-reinforced plate is lower than in the FG fiber-reinforced plate, as shown in Figure 9. The distributions of the therma stress wave are similar when anisotropic, orthotropic, and isotropic fiber-reinforced As shown in Figure 2, the value of the thermal stress wave σ 11 always increases from a positive value. It grows until it reaches its maximum value in the 0 < x 1 < 1 range. Since then, it has been on a downward trend. Finally, it tends to zero as it moves in the direction of wave propagation. In both FG and H anisotropic fiber-reinforced plates, the thermal stress wave σ 11 behaves similarly. As shown in Figure 2, the H reduces the maximum value of the thermal stress wave σ 11 . When the fractional order parameter ∼ α is varied, the distributions of thermal stress wave σ 11 are similar.
The thermal stress wave σ 12 , as shown in Figure 3, has a negative value at x 1 = 0 and then has a downward trend. It moves in the form of wave propagation after a period of rising. The thermal stress wave σ 12 behaves similarly in the FG anisotropic fiber-reinforced plate as it does in the H anisotropic fiber-reinforced plate. As shown in Figure 3, the FG fiber-reinforced anisotropic plate increases the minimum value of the thermal stress wave σ 12 . The distributions of the thermal stress wave σ 12 are similar when the fractional order parameter ∼ α is varied. The thermal stress wave σ 22 exhibits the same behavior in FG and H anisotropic fiber-reinforced plates, as shown in Figure 4. It demonstrates that the value of the thermal stress wave σ 22 reaches a negative value early on and has a downward trend. Since then, it has risen from the lowest to the highest point. Finally, it moves in the direction of wave propagation. Figure 4 shows that the maximum value of the thermal stress wave σ 22 in FG anisotropic fiber-reinforced plates is greater than that in H anisotropic fiber-reinforced plates. The distributions of the thermal stress wave σ 22 are similar when the fractional order parameter ∼ α is varied. In the context of FG anisotropic fiber-reinforced plate, the behavior of the thermal stress wave σ 13 always goes up from a positive value, as shown in Figure 5. In the FG anisotropic fiber-reinforced plate, the value of the thermal stress wave σ 13 in the range of 0 ≤ x 1 ≤ 2.2 is the same as in the FG anisotropic fiber-reinforced plate of Figure 2. From a positive value to its maximum value, it always increases. Then it tends to zero and moves along with the wave. The FG anisotropic fiber-reinforced plate, as shown in Figure 5, increases the amplitude of the thermal stress wave σ 13 . The distributions of the thermal stress wave σ 13 are similar when the fractional order parameter is varied.
The thermal stress wave σ 23 has a negative value at x 1 = 0, as shown in Figure 6. The behavior of the thermal stress wave σ 23 in the FG anisotropic fiber-reinforced plate is identical to that shown in Figure 3. In the absence of gravity, it has a downward trend in the range of 0 ≤ x 1 ≤ 1.7. It gradually decreases after a period of rising. As shown in Figure 6, the FG anisotropic fiber-reinforced plate increases the amplitude of the thermal stress wave σ 23 . When the fractional order parameter is changed, the distributions of the thermal stress wave σ 23 are similar.
As shown in Figure 7, the behavior of the thermal stress wave σ 33 in the FG anisotropic fiber-reinforced plate is identical to that of the FG anisotropic fiber-reinforced plate in Figure 4. In the context of H anisotropic fiber-reinforced plate, it always decreases from a negative value at the start. It goes down to its minimum value in the range of 0 ≤ x 1 ≤ 2.9, then up from the minimum to the maximum. Finally, as the distance x increases, it tends to zero. As shown in Figure 7, the amplitude of the thermal stress wave σ 33 is greater in the FG anisotropic fiber-reinforced plate than in the H anisotropic fiber-reinforced plate. The distributions of the thermal stress wave σ 33 are similar when the fractional order parameter  in the functionally graded (FG) and homogeneous (H) plates for anisotropic, orthotropic, and isotropic materials.
As shown in Figure 8, the thermal stress wave always rises from a positive starting point. It grows until it reaches its maximum value. It has remained on a downward trend since then. Finally, it moves along with the wave propagation and tends to zero The thermal stress wave behaves similarly in FG and H fiber-reinforced plates. As illustrated in Figure 8, the homogeneous case reduces the maximum value of the thermal stress wave . The distributions of the thermal stress wave are similar when anisotropic, orthotropic, and isotropic fiber-reinforced plates are considered. As shown in Figure 9, the thermal stress wave decreases from zero at x=0 and coincides with the boundary condition. It decreases slightly and then rapidly in the range of 0 ≤ ≤ 1.5. Then it decreases and moves along with the wave. In both the FG and H fiber-reinforced plates, the distributions of the thermal stress wave are similar. The maximum value of the thermal stress wave in the H fiber-reinforced plate is lower than in the FG fiber-reinforced plate, as shown in Figure 9. The distributions of the thermal stress wave are similar when anisotropic, orthotropic, and isotropic fiber-reinforced plates are considered.  As shown in Figure 10, the thermal stress wave begins with a positive value. It begins with an upward trend and reaches its maximum value at = 1.1. Finally, it decreases and moves as the wave propagates. It behaves similarly in the H and FG fiberreinforced plates. As illustrated in Figure 10, the H case reduces the maximum value of the thermal stress wave . The distributions of the thermal stress wave are similar when anisotropic, orthotropic, and isotropic fiber-reinforced plates are considered. As shown in Figure 10, the thermal stress wave begins with a positive value. It begins with an upward trend and reaches its maximum value at = 1.1. Finally, it decreases and moves as the wave propagates. It behaves similarly in the H and FG fiberreinforced plates. As illustrated in Figure 10, the H case reduces the maximum value of the thermal stress wave . The distributions of the thermal stress wave are similar when anisotropic, orthotropic, and isotropic fiber-reinforced plates are considered. As shown in Figure 11, the behavior of the thermal stress wave in the FG case is identical to that shown in Figure 8, it always rises from a positive value. It began with an upward trend and then began to decline. Finally, as the distance x increases, it tends to zero. As shown in Figure 11, the FG fiber-reinforced plate increases the amplitude of the thermal stress wave . The distributions of the thermal stress wave are similar when anisotropic, orthotropic, and isotropic fiber-reinforced plates are considered.  The behavior of the thermal stress wave in the FG fiber-reinforced plate is shown in Figure 12, and it is the same as in the FG fiber-reinforced plate shown in Figure 9. It decreases from zero at = 0 to coincide with the boundary condition. It rapidly decreases in the range of 0 ≤ ≤ 2 and then begins to rise. Finally, it decreases and moves as the wave propagates. As shown in Figure 12, the amplitude of the thermal stress wave is greater in the FG fiber-reinforced plate than in the H fiber-reinforced plate case. The distributions of the thermal stress wave are similar when anisotropic, orthotropic, and isotropic fiber-reinforced plates are considered.  The behavior of the thermal stress wave in the FG fiber-reinforced plate is shown in Figure 12, and it is the same as in the FG fiber-reinforced plate shown in Figure 9. It decreases from zero at = 0 to coincide with the boundary condition. It rapidly decreases in the range of 0 ≤ ≤ 2 and then begins to rise. Finally, it decreases and moves as the wave propagates. As shown in Figure 12, the amplitude of the thermal stress wave is greater in the FG fiber-reinforced plate than in the H fiber-reinforced plate case. The distributions of the thermal stress wave are similar when anisotropic, orthotropic, and isotropic fiber-reinforced plates are considered. The behavior of the thermal stress wave in the FG fiber-reinforced plate is the same as in the FG fiber-reinforced plate of Figure 10, as shown in Figure 13. In the H fiberreinforced plate, it falls at first and reaches its lowest point at = 0.7. Following that, it rises and tends to zero as the distance increases. As shown in Figure 13, the FG fiber-   As shown in Figure 14, the behavior of the thermal stress wave in the FG anisotropic fiber-reinforced plate is identical to that shown in Figure 5. It reaches its peak in the range 0 ≤ ≤ 2.2, and the value of the thermal stress wave continues to fall. Following that, it tends to zero and moves in the wave propagation. As shown in Figure 14  As shown in Figure 15, the thermal stress wave in the H anisotropic fiber-reinforced plate has an upward trend in the range of 0 ≤ ≤ 1.9. When it reaches its maximum value, it gradually decreases and tends to zero. As shown in Figure 15  As shown in Figure 8, the thermal stress wave σ 11 always rises from a positive starting point. It grows until it reaches its maximum value. It has remained on a downward trend since then. Finally, it moves along with the wave propagation and tends to zero. The thermal stress wave σ 11 behaves similarly in FG and H fiber-reinforced plates. As illustrated in Figure 8, the homogeneous case reduces the maximum value of the thermal stress wave σ 11 . The distributions of the thermal stress wave σ 11 are similar when anisotropic, orthotropic, and isotropic fiber-reinforced plates are considered.
As shown in Figure 9, the thermal stress wave σ 12 decreases from zero at x = 0 and coincides with the boundary condition. It decreases slightly and then rapidly in the range of 0 ≤ x 1 ≤ 1.5. Then it decreases and moves along with the wave. In both the FG and H fiber-reinforced plates, the distributions of the thermal stress wave σ 12 are similar. The maximum value of the thermal stress wave σ 12 in the H fiber-reinforced plate is lower than in the FG fiber-reinforced plate, as shown in Figure 9. The distributions of the thermal stress wave σ 12 are similar when anisotropic, orthotropic, and isotropic fiber-reinforced plates are considered.
As shown in Figure 10, the thermal stress wave σ 22 begins with a positive value. It begins with an upward trend and reaches its maximum value at x 1 = 1.1. Finally, it decreases and moves as the wave propagates. It behaves similarly in the H and FG fiberreinforced plates. As illustrated in Figure 10, the H case reduces the maximum value of the thermal stress wave σ 22 . The distributions of the thermal stress wave σ 22 are similar when anisotropic, orthotropic, and isotropic fiber-reinforced plates are considered.
As shown in Figure 11, the behavior of the thermal stress wave σ 13 in the FG case is identical to that shown in Figure 8, it always rises from a positive value. It began with an upward trend and then began to decline. Finally, as the distance x increases, it tends to zero. As shown in Figure 11, the FG fiber-reinforced plate increases the amplitude of the thermal stress wave σ 13 . The distributions of the thermal stress wave σ 13 are similar when anisotropic, orthotropic, and isotropic fiber-reinforced plates are considered.
The behavior of the thermal stress wave σ 23 in the FG fiber-reinforced plate is shown in Figure 12, and it is the same as in the FG fiber-reinforced plate shown in Figure 9. It decreases from zero at x 1 = 0 to coincide with the boundary condition. It rapidly decreases in the range of 0 ≤ x 1 ≤ 2 and then begins to rise. Finally, it decreases and moves as the wave propagates. As shown in Figure 12, the amplitude of the thermal stress wave σ 23 is greater in the FG fiber-reinforced plate than in the H fiber-reinforced plate case. The distributions of the thermal stress wave σ 23 are similar when anisotropic, orthotropic, and isotropic fiber-reinforced plates are considered.
The behavior of the thermal stress wave σ 33 in the FG fiber-reinforced plate is the same as in the FG fiber-reinforced plate of Figure 10, as shown in Figure 13. In the H fiber-reinforced plate, it falls at first and reaches its lowest point at x 1 = 0.7. Following that, it rises and tends to zero as the distance x 1 increases. As shown in Figure 13, the FG fiber-reinforced plate increases the amplitude of the thermal stress wave σ 33 . The distributions of the thermal stress wave σ 33 are similar when anisotropic, orthotropic, and isotropic fiber-reinforced plates are considered. Figures 14-19 display the variations of the thermal stress waves σ 11 , σ 12 , σ 22 , σ 13 , σ 23 and σ 33 in the functionally graded (FG) and homogeneous (H) anisotropic fiber-reinforced plates under different values of the time characteristic of the laser pulse τ 1 .  As shown in Figure 14, the behavior of the thermal stress wave in the FG anisotropic fiber-reinforced plate is identical to that shown in Figure 5. It reaches its peak in the range 0 ≤ ≤ 2.2, and the value of the thermal stress wave continues to fall. Following that, it tends to zero and moves in the wave propagation. As shown in Figure 14, the time characteristic of the laser pulse influences the behavior of the thermal stress wave . The distributions of the thermal stress wave are similar when the time characteristic of the laser pulse changes. As shown in Figure 15, the thermal stress wave in the H anisotropic fiber-reinforced plate has an upward trend in the range of 0 ≤ ≤ 1.9. When it reaches its maximum value, it gradually decreases and tends to zero. As shown in Figure 15, the time characteristic of the laser pulse causes the value of the thermal stress wave to increase. The distributions of the thermal stress wave are similar when the time characteristic of the laser pulse changes. As shown in Figure 16, the thermal stress wave decreases slightly and then rapidly increases to its maximum value in the range of 0 ≤ ≤ 1.7. Finally, it approaches zero and moves in the wave propagation. As shown in Figure 16, the time characteristic of the laser pulse causes the maximum value of the FG anisotropic fiber-reinforced plate of thermal stress wave to decrease. The distributions of the thermal stress wave are similar when the time characteristic of the laser pulse changes. As shown in Figure 16, the thermal stress wave decreases slightly and then rapidly increases to its maximum value in the range of 0 ≤ ≤ 1.7. Finally, it approaches zero and moves in the wave propagation. As shown in Figure 16, the time characteristic of the laser pulse causes the maximum value of the FG anisotropic fiber-reinforced plate of thermal stress wave to decrease. The distributions of the thermal stress wave are similar when the time characteristic of the laser pulse changes. The thermal stress wave , as shown in Figure 17, decreases from a positive value in the H anisotropic fiber-reinforced plate. In the H anisotropic fiber-reinforced plate, it decreases to its minimum value in the range of 0 ≤ ≤ 2. Following that, it grows and moves in the wave propagation. As shown in Figure 17, the FG anisotropic fiber-reinforced plate increases the maximum value of the thermal stress wave . The distributions of the thermal stress wave are similar when the time characteristic of the laser pulse changes. As shown in Figure 18, the thermal stress wave decreases from zero to coincide with the boundary condition. It rapidly decreases and then increases. Finally, it moves in the direction of wave propagation. Figure 18 shows that the maximum value of the thermal stress wave in the FG anisotropic fiber-reinforced plate is greater than that in the H fiber-reinforced plate. The distributions of the thermal stress wave are similar when the time characteristic of the laser pulse changes. As shown in Figure 18, the thermal stress wave decreases from zero to coincide with the boundary condition. It rapidly decreases and then increases. Finally, it moves in the direction of wave propagation. Figure 18 shows that the maximum value of the thermal stress wave in the FG anisotropic fiber-reinforced plate is greater than that in the H fiber-reinforced plate. The distributions of the thermal stress wave are similar when the time characteristic of the laser pulse changes. The behavior of the thermal stress wave , as shown in Figure 19, is identical to that shown in Figure 13. In the FG and H anisotropic fiber-reinforced plates, it has a positive value at first. It rises and reaches its highest point at = 0.4. Then it diminishes and   Table 2 shows a comparison of required computer resources for the current BEM re sults, and FEM-NMM results of An et al. [31] for the modeling of ultrasonic wave propa gation fractional order boundary value problems of FGA plates. There were no published results to support the validity of the proposed technique's findings. Some papers, on the other hand, can be regarded as subsets of the larger study under consideration. The variations of the special case thermal stress waves , , and along the axis for BEM and combined finite element method/normal mode method (FEM-NMM) in fractional order ( = 0.6) functionally graded plates are shown in Figures 20-22, respectively. These results show that the BEM findings agree very wel with the FEM-NMM findings of An et al. [31]. As a result, the proposed technique's va lidity was confirmed. We refer the interested readers to the references of fractional deriv ative of the Riemann Zeta function [32][33][34], fractional derivatives in complex planes [35,36] and fractional boundary element method [37,38]. As shown in Figure 14, the behavior of the thermal stress wave σ 11 in the FG anisotropic fiber-reinforced plate is identical to that shown in Figure 5. It reaches its peak in the range 0 ≤ x 1 ≤ 2.2, and the value of the thermal stress wave σ 11 continues to fall. Following that, it tends to zero and moves in the wave propagation. As shown in Figure 14, the time characteristic of the laser pulse influences the behavior of the thermal stress wave σ 11 . The distributions of the thermal stress wave σ 11 are similar when the time characteristic of the laser pulse τ 1 changes.
As shown in Figure 15, the thermal stress wave σ 12 in the H anisotropic fiber-reinforced plate has an upward trend in the range of. 0 ≤ x 1 ≤ 1.9. When it reaches its maximum value, it gradually decreases and tends to zero. As shown in Figure 15, the time characteristic of the laser pulse causes the value of the thermal stress wave σ 12 to increase. The distributions of the thermal stress wave σ 12 are similar when the time characteristic of the laser pulse τ 1 changes.
As shown in Figure 16, the thermal stress wave σ 22 decreases slightly and then rapidly increases to its maximum value in the range of 0 ≤ x 1 ≤ 1.7. Finally, it approaches zero and moves in the wave propagation. As shown in Figure 16, the time characteristic of the laser pulse causes the maximum value of the FG anisotropic fiber-reinforced plate of thermal stress wave σ 22 to decrease. The distributions of the thermal stress wave σ 22 are similar when the time characteristic of the laser pulse τ 1 changes.
The thermal stress wave σ 13 , as shown in Figure 17, decreases from a positive value in the H anisotropic fiber-reinforced plate. In the H anisotropic fiber-reinforced plate, it decreases to its minimum value in the range of 0 ≤ x ≤ 2. Following that, it grows and moves in the wave propagation. As shown in Figure 17, the FG anisotropic fiber-reinforced plate increases the maximum value of the thermal stress wave σ 13 . The distributions of the thermal stress wave σ 13 are similar when the time characteristic of the laser pulse τ 1 changes.
As shown in Figure 18, the thermal stress wave σ 23 decreases from zero to coincide with the boundary condition. It rapidly decreases and then increases. Finally, it moves in the direction of wave propagation. Figure 18 shows that the maximum value of the thermal stress wave σ 23 in the FG anisotropic fiber-reinforced plate is greater than that in the H fiber-reinforced plate. The distributions of the thermal stress wave σ 23 are similar when the time characteristic of the laser pulse τ 1 changes.
The behavior of the thermal stress wave σ 33 , as shown in Figure 19, is identical to that shown in Figure 13. In the FG and H anisotropic fiber-reinforced plates, it has a positive value at first. It rises and reaches its highest point at x = 0.4. Then it diminishes and moves in the wave propagation. As shown in Figure 19, the FG anisotropic fiber-reinforced plate increases the maximum value of the thermal stress wave σ 33 . The distributions of the thermal stress wave σ 33 are similar when the time characteristic of the laser pulse τ 1 takes different values. Table 2 shows a comparison of required computer resources for the current BEM results, and FEM-NMM results of An et al. [31] for the modeling of ultrasonic wave propagation fractional order boundary value problems of FGA plates. There were no published results to support the validity of the proposed technique's findings. Some papers, on the other hand, can be regarded as subsets of the larger study under consideration. The variations of the special case thermal stress waves σ 11 , σ 12 , and σ 22 along the x 1 −axis for BEM and combined finite element method/normal mode method (FEM-NMM) in fractional order (a = 0.6) functionally graded plates are shown in Figures 20-22, respectively. These results show that the BEM findings agree very well with the FEM-NMM findings of An et al. [31]. As a result, the proposed technique's validity was confirmed. We refer the interested readers to the references of fractional derivative of the Riemann Zeta function [32][33][34], fractional derivatives in complex planes [35,36] and fractional boundary element method [37,38].

Conclusions
Some of the conclusions that can be derived from this paper are as follows: 1. A numerical BEM scheme based on MLS is applied to FGA fiber-reinforced plates under thermoelastic loads. The Reissner-Mindlin theory, which considers shear deformation, is used to explain the behavior. 2. The Laplace-transform is used to remove the time variable from the governing equations. 3. The domain under consideration is subdivided into small circular subdomains. The local boundary integral equations are derived using the unit step test function. 4. The MLS scheme has been proposed for treating the domain integrals arising from the inertial term and approximates the physical quantities. 5. Numerical results demonstrate the accuracy, feasibility, effectiveness, and convergence of the present MLS method. 6. The boundary value problems must be solved for a variety of Laplace-transform parameter values chosen for each time instant under consideration. 7. The primary benefit of the current technique is its generality and simplicity. The used test function in the proposed technique is less complicated than the fundamental solution for anisotropic fiber-reinforced plates. 8. Numerical results demonstrate that the fractional order parameter, functionally graded material, anisotropy, and time characteristic of the laser pulse have significant effects on the thermal stresses of FGA fiber-reinforced plates. 9. The numerical results confirm that the proposed technique provides more benefits than other domain discretization methods. 10. The results presented in this paper may provide interesting information for research- Figure 22. Variation of the thermal stress wave σ 22 along x 1 −axis in the special case of the functionally graded anisotropic fiber-reinforced plates for BEM and FEM-NMM.

Conclusions
Some of the conclusions that can be derived from this paper are as follows:

1.
A numerical BEM scheme based on MLS is applied to FGA fiber-reinforced plates under thermoelastic loads. The Reissner-Mindlin theory, which considers shear deformation, is used to explain the behavior.

2.
The Laplace-transform is used to remove the time variable from the governing equations.

3.
The domain under consideration is subdivided into small circular subdomains. The local boundary integral equations are derived using the unit step test function. 4.
The MLS scheme has been proposed for treating the domain integrals arising from the inertial term and approximates the physical quantities.

5.
Numerical results demonstrate the accuracy, feasibility, effectiveness, and convergence of the present MLS method. 6.
The boundary value problems must be solved for a variety of Laplace-transform parameter values chosen for each time instant under consideration. 7.
The primary benefit of the current technique is its generality and simplicity. The used test function in the proposed technique is less complicated than the fundamental solution for anisotropic fiber-reinforced plates.

8.
Numerical results demonstrate that the fractional order parameter, functionally graded material, anisotropy, and time characteristic of the laser pulse have significant effects on the thermal stresses of FGA fiber-reinforced plates. 9.
The numerical results confirm that the proposed technique provides more benefits than other domain discretization methods. 10. The results presented in this paper may provide interesting information for researchers who are working on computer science, material science, mathematical physics, geotechnical engineering, and geothermal engineering as well as for those working on the development of the functionally graded anisotropic fiber-reinforced plates.