Thermoelastic Analysis of Functionally Graded Cylindrical Panels with Piezoelectric Layers

: We propose a coupled thermoelastic approach based on the Lord-Shulman (L-S) and Maxwell’s formulations to study the wave propagation in functionally graded (FG) cylindrical panels with piezoelectric layers under a thermal shock loading. The material properties of the FG core layer feature a graded distribution throughout the thickness and vary according to a simple power law. A layerwise di ﬀ erential quadrature method (LW-DQM) is combined with a non-uniform rational B-spline (NURBS) multi-step time integration scheme to discretize the governing equations both in the spatial and time domains. The compatibility conditions of the physical quantities are enforced at the interfaces to describe their structural behavior in a closed form. A validation and comparative analysis with the available literature, together with a convergence study, show the e ﬃ ciency and stability of the proposed method to handle thermoelastic problems. Numerical applications are herein performed systematically to check for the sensitivity of the thermoelastic response to the material graded index, piezoelectric layer thickness, external electrical voltage, opening angle, and shock thermal loading, which would be very helpful for practical engineering applications. coupled governing equations of solved using hybrid method built upon a discrete LW-DQM method and a NURBS-based multi-step time integration scheme. The geometric and physical compatibility conditions have been enforced exactly at the interfacial layers. The model is able to accurately simulate the variations of the field variables throughout the thickness, as demonstrated by the fast rate of convergence of the numerical results. A wide study has considered the effects of


Introduction
Functionally graded (FG) cylindrical panels are largely used as structural members in many industrial applications such as aerospace vehicles, nuclear equipment, petro-chemical structures, among others. FG structures usually operate under high thermal conditions, which can yield a vibratory motion, especially when subjected to a sudden variation of the thermal conditions. On the other hand, the mechanical behavior and the vibration control of FG cylindrical panels can be improved by introducing piezoelectric layers at their inner and/or outer surfaces. This requires an accurate evaluation of the thermoelastic properties of FG cylindrical panels, with attached piezoelectric layers under a thermal shock loading, for design and manufacturing purposes.
In this context, it is well known that a conventional Fourier heat conduction law is unable to predict the temperature distribution and heat wave propagation within a continuum at the initial phase. On the other hand, an uncoupled classical thermoelastic approach cannot accurately predict the thermoelastic behavior of a body under a thermal shock. Countermeasures against these limitations have encouraged the development of different coupled thermoelasticity theories [1][2][3]. Among them, the Lord and Shulman (L-S) generalized thermoelasticity theory is one of the simplest formulations applied in the literature [4][5][6][7][8][9][10][11][12][13][14], which includes a relaxation time and a coupling between thermal and mechanical energies [1]. In this framework, some interesting thermoelastic studies based on an uncoupled classical thermoelasticity theory of FG cylindrical panels and shells can be found in [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31]. These studies considered the possible presence of piezoelectric layers, and estimated the temperature distribution of a body by means of a Fourier heat conduction law. On the other hand, many works, in literature, have applied a non-Fourier heat transfer law or coupled theories to study the thermoelastic response of FG structures [32][33][34][35][36][37]. More specifically, Entezari et al. [32] explored the capabilities of a refined finite element approach based on the L-S generalized theory of thermoelasticity and Carrera Unified Formulation (CUF) for the 3D thermoelastic analysis of FG disks exposed to a thermal shock loading. Zhang et al. [33] employed a generalized fractional Cattaneo-type heat conduction model to analyze the thermal shock problem of elastic FG strips with internal cracks, while using a Fourier and Laplace transforms for an analytical computation of the thermoelastic field unknowns. Sator et al. [34] presented a unified formulation for the bending response of FG plates under a transient thermal loading in the framework of the L-S generalized thermoelasticity and two-dimensional plate theories. In their work, the authors applied a meshless method and a Wilson time integration scheme to solve numerically the governing equations of the problem. In line with the previous works, Pourasghar and Chen [35] applied the non-Fourier heat conduction equations to determine the temperature distribution and the induced thermal stresses and deformations, in carbon nanotube-reinforced composite cylindrical panels subjected to a heat pulse. They assumed temperature-dependent material properties and applied a differential quadrature method (DQM) to discretize the governing equations and solve the problem numerically. Heydarpour and Malekzadeh [36] studied the propagation of thermoelastic waves in laminated spherical shells with temperature-dependent material properties and FG layers under thermal boundary conditions according to the L-S generalized thermoelasticity. Heydarpour et al. [37] investigated the thermo-mechanical wave propagation in composite spherical shells reinforced with multilayer functionally graded graphene platelets under a thermo-mechanical shock loading. They formulated the problem in the context of L-S thermoelasticity together with a layerwise (LW) DQM, where the governing equations were solved in the time domain both analytically and numerically. The non-linear behavior of the shell structures under different loading conditions, makes the wave propagation problem particularly sensitive to the selected mechanical and geometrical parameters. This represents a key aspect of higher-order numerical approaches for shell structures [38][39][40][41][42][43].
To date, however, there has been a general lack of works applying coupled thermoelastic approaches to FG cylindrical panels with surface piezoelectric layers under a thermal shock. Due to its practical importance and usefulness from a research and practical engineering standpoint, this issue is investigated in the present paper. In order to describe the high transient nature of a thermal loading, especially at the initial phase of its application, the problem is formulated according to the L-S thermoelasticity theory, with a special emphasis on the heat wave propagation. The material properties of FG panels are assumed to vary throughout the thickness, and a layerwise (LW) approach is adopted to simulate the variation of field variables throughout the thickness, due to its high accuracy and computational efficiency compared to a conventional finite element method (FEM) or finite difference method (FDM), as verified for different engineering problems. The LW theory, indeed, has been largely proved to be an efficient tool for the study of the response and transverse the shear and normal stress distributions of laminated composites. It represents a refined approach that considers the effects of thickness at a reduced computational cost. Based on a LW theory, the shear strains are discontinuous, and hence the shear stresses can be continuous at the interface of two adjacent layers. This represents a clear advantage of the LW approach with respect to higher order shear deformation theories. At the same time, the successful employment of the DQM for different structural and fluid problems has largely been proved in literature [11,20,36,[43][44][45][46][47][48][49][50][51][52].
Accordingly, a LW-DQM is proposed to discretize the governing equations in the spatial domain, which are solved using a multi-step Non-Uniform-Rational-B-Spline (NURBS)-based technique. The correctness and accuracy of the present approach is endowed with a convergence study and comparative analysis of results with respect to the available predictions from the literature in the limit cases. Afterward, the effect of the material gradient index, relaxation time, opening angle, piezoelectric layer thickness, and external electrical voltage on the thermoelastic behavior of the system is investigated parametrically.
The outline of the paper is as follows: in Section 2 we describe the theoretical basics and the governing equations of the problem. In Section 3 we present the numerical procedure of a NURBS-based multi-step method and DQM, whose numerical results are discussed and compared in Section 4 for different mechanical, geometrical, and thermal parameters. The main conclusions of our work are, finally, summarized in Section 5.

Governing Equations
Let us describe the mathematical problem and its governing equations for a long FG cylindrical panel with an inner and outer piezoelectric layer. The laminated panel features an opening angle θ 0 , an inner radius R i , and an outer radius R o , whereby h c is the FG layer thickness, h p is the piezoelectric layer thickness, and for a total thickness h = h c + 2h p (see Figure 1). A cylindrical coordinate system (r, θ, z) is chosen to identify the material points of the laminated panel in the undeformed configuration, as visible in Figure 1. The FG core layer is assumed to change gradually throughout the radial direction, namely, the mechanical properties vary continuously and smoothly from a ceramic-rich material (at the inner surface) to a metal-rich material (at the outer surface) according to a power law distribution. Thus, the effective material property P of the core is defined, accordingly, as where n is a positive real number associated with the power law index (or the material property gradient index).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 26 piezoelectric layer thickness, and external electrical voltage on the thermoelastic behavior of the system is investigated parametrically. The outline of the paper is as follows: in Section 2 we describe the theoretical basics and the governing equations of the problem. In Section 3 we present the numerical procedure of a NURBSbased multi-step method and DQM, whose numerical results are discussed and compared in Section 4 for different mechanical, geometrical, and thermal parameters. The main conclusions of our work are, finally, summarized in Section 5.

Governing Equations
Let us describe the mathematical problem and its governing equations for a long FG cylindrical panel with an inner and outer piezoelectric layer. The laminated panel features an opening angle 0 θ , an inner radius Ri, and an outer radius Ro, whereby hc is the FG layer thickness, hp is the piezoelectric layer thickness, and for a total thickness h = hc + 2hp (see Figure 1). A cylindrical coordinate system ,θ is chosen to identify the material points of the laminated panel in the undeformed configuration, as visible in Figure 1. The FG core layer is assumed to change gradually throughout the radial direction, namely, the mechanical properties vary continuously and smoothly from a ceramic-rich material (at the inner surface) to a metal-rich material (at the outer surface) according to a power law distribution. Thus, the effective material property P of the core is defined, accordingly, as where n is a positive real number associated with the power law index (or the material property gradient index).  Long panels follow the plane strain assumptions, as employed thereinafter. Accordingly, the linear constitutive relations for the core (i.e., the FG layer) and the piezoelectric layers take the following form, respectively [30]: Long panels follow the plane strain assumptions, as employed thereinafter. Accordingly, the linear constitutive relations for the core (i.e., the FG layer) and the piezoelectric layers take the following form, respectively [30]: where (i, j = r, θ) is the stress tensor components in the FG core (e = c) and piezoelectric layers (e = p); u and v are the displacement components at an arbitrary material point of the panel along the r and θ-directions, respectively; α denotes the thermal expansion coefficient; (∆T) e = T e − T 0 is the temperature variation; T 0 is the shell stress free temperature; C e ij (i, j = 1, 2, 5) is the material elastic coefficients of the core and piezoelectric layers; D i and E i (i = r, θ) are, the electric displacement and electric field in the ith direction, respectively; and e i (i = 2, 3, 4), K ii (i = 2, 3) and P 3 refer to the piezoelectric, dielectric, and pyroelectric constants, respectively.
The following relations are introduced between the electric field E i and the electric potential ψ in piezoelectric layers: The FG cylindrical panel is subjected to an initial temperature T 0 before the variation of the thermal boundary conditions at the two external surfaces occurs. Under the plane strain assumptions, the governing equations are independent of the axial coordinate variable z. Thus, for the FG core layer and piezoelectric layers, the linearized coupled energy balance equation according to the L-S thermoelasticity theory together with the equation of motion can be defined as [1].
Energy balance equation: where α e = α e (3λ e + 2µ e ), λ e and µ e are the Lamé constants of the layers (e = p, c); τ e 0 (e = p, c) is the relaxation time of the L-S theory, and t is the time.
Moreover, the piezoelectric layers are governed by the following electroelastic equation: Among different possibilities of assuming the thermal boundary conditions of the panel, in the present work we consider a drastic variation for the transient temperature (at the inner surface) and for the convection heat transfer (at the outer surface). In detail, at where T i is the temperature of the inner surface over a long period of time (i.e., t >> t 0 ), t 0 is a time constant, T ∞ is the ambient temperature, andĥ c denotes the convective heat transfer coefficient. By neglecting the heat transfer from the straight edges of the panel (i.e., θ = 0 and θ = θ 0 ), the corresponding thermal condition becomes The mechanical boundary conditions at the inner and outer surfaces of the panels (r = R i , R o ) are defined as follows: Moreover, the following thermal and mechanical compatibility conditions are enforced at the interface between the FG core and piezoelectric layers of the panel (i.e., In this study, clamped boundary conditions are assumed at the edges θ = 0 and θ = θ 0 , namely, whereby the initial thermal and mechanical conditions are assumed, respectively, as follows:

Solution Procedure
The analytical solution of the abovementioned governing equations of the problem with variable coefficients is clearly difficult and cumbersome, whereas a semi-analytical or numerical method should be more convenient and suitable for use. Due to the computational efficiency and accuracy of the DQM demonstrated in the literature for different heat transfer and structural problems [11,20,36,[43][44][45][46][47][48][49][50][51][52], we combined this method with the NURBS-based multi-step method in the current work to discretize the governing equations in the spatial and time domains, respectively. Based on the DQM, the FG core and piezoelectric layers are discretized into a set of N e ξ (e = c, p) and N θ grid points along the rand θ-directions, respectively. In this regard, the Gauss-Lobatto-Chebyshev grid generation rule is employed [11,20,36,[43][44][45][46][47][48][49][50][51][52]. Thus, the spatial derivatives involved in the differential equations are discretized together with the corresponding thermo-mechanical boundaries and compatibility conditions at the interface between the core and piezoelectric layers. For the sake of brevity, we describe only the discretized form for the of energy balance from Equation (4), as also repeated for all the other equations of the problem.
where e = c, p; represent the first-and second-order weighting coefficients for the β-direction (β = r, θ), respectively, which is computed here through the Lagrange polynomials.
Once the DQM discretization procedure is completed, we get the following system of ordinary differential equations (ODEs) in the time domain: where D is the vector of degrees of freedom (i.e., temperature and displacement components at grid points); M, C, and K represent the mass, thermoelastic damping, and stiffness matrices, respectively; and f(t) is the load vector. A NURBS-based multi-step method [37] is here adopted to solve the system of ODE (17) according to the initial conditions. The higher performance of this technique with respect to the conventional finite difference-based schemes, such as the Newmark's time integration scheme, has been recently demonstrated by Heydarpour et al. [37]. To apply this method, the second-order system of Equation (17) is transformed into a first-order one, as follows: where Ψ 1 = D and Ψ 2 = dD dt . Based on the present approach, a different multi-step scheme can be generated by changing the order and/or the weighting coefficients (w i ) of the NURBS curves. In this study, a four-step scheme [37] with the weight coefficients w 1 = 0.001, w 2 = 0.001, w 3 = 2 and w 4 = 3 is considered to solve Equation (18a,b), the discretized form of which can be expressed as follows: where ∆t is the time step size, andΨ The solution process starts with the initial conditions to determine the unknown field variables at the first four points. Subsequently, by employing the multi-step scheme, Equation (18a,b) is converted into a system of linear algebraic equations whose solution is related to the field variables at the end of each time step. By repeating the process, we are able to determine the displacement components and temperature at the grid points of FG cylindrical panel.

Numerical Results
We start this section with a validation of the proposed approach, which is followed by a parametric study that investigates the sensitivity of the thermoelastic response of the panel to the mechanical, thermal, or geometrical input parameters. Unless otherwise mentioned, we consider Ti-6Al-4V and ZrO 2 in the numerical examples to be the metal and ceramic phases, respectively (Table 1), and we introduce the following dimensionless parameters whereα m (= k m /ρ m c m ) is the thermal diffusivity of the metal, and the following geometry and physical parameters are assumed: R m = 1 (m) , h = 0.35 m , T 0 = T ∞ = 300K, T i = 800K,ĥ c = 10 W/m 2 K, t 0 = 2s. The piezoelectric layers are made of lead zirconate titanate (PZT-4) with the mechanical properties given in Table 2, in accordance with [30]. As a first example, we compare the CPU time necessary for a thermoelastic study of a FG cylindrical panel, as provided by a NURBS-based multi-step time integration scheme and a Newmark's method. Table 3 gathers the convergence results for the dimensionless displacement, hoop stress, and temperature at point ξ = 0.5, θ = 0.5θ 0 versus the number of time steps N t . As visible in Table 3, a NURBS-based multi-step approach requires a lower CPU time compared to a classical Newmark's method. Accordingly, Figures 2 and 3 depict the convergence results from a transient thermoelastic study of FG cylindrical panels with a piezoelectric layer under a thermal loading condition. In these figures, the variations along the thickness direction of the dimensionless temperature, radial displacement, and hoop stress components are reported versus the DQ number of grid points throughout the thickness N e ξ = N ξ , e = c, p and circumferential direction (N θ ). It is worth observing that nine grid points in the radial direction per layer (i.e., N ξ = 9) and eleven grid points in the circumferential direction yielded acceptable results. Hence, the values of N ξ = 9, N θ = 11, and N t = 300 are assumed hereafter within the numerical investigation. Table 1. Material properties of ceramic (ZrO 2 ) and metal (Ti-6Al-4V) [11].

Material
Ti-6Al-4V    Table 3. Comparison between a NURBS-based multi-step technique and the Newmark's scheme n = 1, h c /h p = 5, θ 0 = 75 • , ψ = 0.03, ξ = 0.5, θ = 0.5θ 0 , F 0 = 0.8, τ = 0.1 .    For validation purposes, we analyzed the FG hollow cylindrical shells under a thermal loading, as a limit case of cylindrical panels. Therefore, by setting θ 0 = 2π and enforcing the geometric and physical compatibility conditions at θ = 0 and θ = θ 0 = 2π, a FG cylindrical panel reverted to a cylindrical shell. This case study is here selected in agreement with [53], where a semi-analytical finite element method was applied for the analysis of FG cylindrical shells under the following thermal boundary and initial conditions:

NURBS-Based Multi−Step Method Newmark's Method (Galerkin Scheme)
As far as the material properties are concerned, they vary throughout the shell thickness from a ceramic-rich material, at the inner surface, to a metal-rich material, at the outer surface, in accordance with Equation (1). The metal and ceramic mechanical properties are listed in Table 4. The dimensionless temperature and displacement distributions across the shell thickness are depicted in Figure 4 for a varying material graded index n, where the accuracy of the proposed method against the uncoupled theory proposed in [53] is confirmed by the very good agreement between results.
As far as the material properties are concerned, they vary throughout the shell thickness from a ceramic-rich material, at the inner surface, to a metal-rich material, at the outer surface, in accordance with Equation (1). The metal and ceramic mechanical properties are listed in Table 4. The dimensionless temperature and displacement distributions across the shell thickness are depicted in Figure 4 for a varying material graded index n, where the accuracy of the proposed method against the uncoupled theory proposed in [53] is confirmed by the very good agreement between results.    The effect of the material gradient index n on the through-the-thickness variation and time histories of the thermoelastic field variables is plotted in Figures 5 and 6 for FG cylindrical panels with piezoelectric layers under a thermal shock at the inner surface. The heat wave front can easily be seen in Figures 5a and 6a, where panels with a pure metal core (i.e., n = 0) feature a higher value of heat wave speed compared to the ones with nonhomogeneous cores (i.e., n 0). In addition, the heat wave speed increases by increasing the metal content, or equivalently with a decreasing value of n.    (a-c): Effect of the material gradient index n on the time histories of the thermoelastic field variables for the FG cylindrical panels with piezoelectric layers h c /h p = 5, θ 0 = 75 • , ψ = 0.03, ξ = 0.5, θ = 0.5θ 0 , τ = 0.2 . Figures 7 and 8 also plot the effect of the dimensionless relaxation time τ on the thermoelastic response of FG cylindrical panels with piezoelectric layers, revealing the through-the-thickness variations and time histories of the thermoelastic field variables of a material point (ξ = 0.5, θ = 0.5θ 0 ), respectively. Based on Figures 7 and 8, it is clear that the use of either a classical thermoelasticity theory (i.e., τ = 0) or the L-S theory (i.e., the case with τ 0) leads to significant differences in the numerical results, evaluated here in terms of temperature, radial displacement, and hoop stress. Moreover, the oscillating behavior of the temperature time histories for τ 0, (i.e., for a non-Fourier heat conduction law) plotted in Figure 8a, indicates the necessity of using the L-S theory for thermal shock loading problems.    In Figure 10, the effect of the core thickness-to-piezoelectric thickness ratio h c /h p on the thermoelastic behavior of FG cylindrical panels is described, where an increased value of h c /h p yields a general increase of all the field variables. A further systematic study considers the effect of the panel opening angle θ 0 on the through-the-thickness variation of the dimensionless thermoelastic field variables, as depicted in Figure 11a-c. As expected, by increasing θ 0 , the radial displacement increases, and the hoop stress component decreases, although the maximum values of the field variables involve the same locations of the panel independently of θ 0 . Furthermore, in Figure 12a-c, we evaluate the effect of the opening angle θ 0 on the time histories for all thermoelastic field variables while considering the material points ξ = 0.5 and θ = 0.5θ 0 . Also visible in Figures 11a and 12a, the distribution and time histories of the temperature seem to be unaffected by a variation in θ 0 .      In the last parametric study, we checked the sensitivity of the thermoelastic response of FG cylindrical panels with piezoelectric layers to varying external electrical voltages (Figures 13 and 14). More specifically, we evaluated the through-the-thickness variations and time histories of the dimensionless thermoelastic field variables for three different values of external electrical voltage 0; 0.03; 0.06 ψ = . As clearly visible in Figures 13a and 14a, an increase in the external electrical voltage did not affect the temperature variations, but it led to an overall decrease in radial displacement and an increase of the hoop stress component. In the last parametric study, we checked the sensitivity of the thermoelastic response of FG cylindrical panels with piezoelectric layers to varying external electrical voltages (Figures 13  and 14). More specifically, we evaluated the through-the-thickness variations and time histories of the dimensionless thermoelastic field variables for three different values of external electrical voltage ψ = 0; 0.03; 0.06. As clearly visible in Figures 13a and 14a, an increase in the external electrical voltage did not affect the temperature variations, but it led to an overall decrease in radial displacement and an increase of the hoop stress component.

Conclusions
In this paper, we have proposed a coupled thermoelasticity theory of Lord-Shulman and Maxwell's equation, to examine thermoelastic wave propagation for FG cylindrical panels with inner and outer surfaces bonded by piezoelectric layers, and subjected the panels to a thermal shock loading. The coupled governing equations of the problem have been solved using a computationally efficient and accurate numerical hybrid method built upon a discrete LW-DQM method and a NURBS-based multi-step time integration scheme. The geometric and physical compatibility conditions have been enforced exactly at the interfacial layers. The model is able to accurately simulate the variations of the field variables throughout the thickness, as demonstrated by the fast rate of convergence of the numerical results. A wide parametric study has considered the effects of the material graded index, piezoelectric layer thickness, external electrical voltage, and opening angle

Conclusions
In this paper, we have proposed a coupled thermoelasticity theory of Lord-Shulman and Maxwell's equation, to examine thermoelastic wave propagation for FG cylindrical panels with inner and outer surfaces bonded by piezoelectric layers, and subjected the panels to a thermal shock loading. The coupled governing equations of the problem have been solved using a computationally efficient and accurate numerical hybrid method built upon a discrete LW-DQM method and a NURBS-based multi-step time integration scheme. The geometric and physical compatibility conditions have been enforced exactly at the interfacial layers. The model is able to accurately simulate the variations of the field variables throughout the thickness, as demonstrated by the fast rate of convergence of the numerical results. A wide parametric study has considered the effects of the material graded index, piezoelectric layer thickness, external electrical voltage, and opening angle on the thermoelastic behavior of FG cylindrical panels, while providing the following promising conclusions for practical applications:

•
The model can simulate thermoelastic wave propagation in multilayer panels.

•
An increased material gradient index n decreases the heat wave speed.

•
The classical thermoelasticity theory based on a null value of relaxation time is unable to accurately predict the thermoelastic behavior of laminated panels under a thermal shock loading.

•
An increased core thickness-to-piezoelectric thickness ratio h c /h p yields an increased value of all the thermoelastic field variables. Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.