Semianalytical Research on Aerothermoelastic Behaviors of Functionally Graded Plates under Arbitrary Temperature Fields in Hypersonic Vehicles

: Hypersonic vehicles are susceptible to considerable aerodynamic heating and noticeable aerothermoelastic effects during flight due to their high speeds. Functionally graded materials (FGMs), which enable continuous changes in material properties by varying the ratio of different materials, provide both thermal protection and load-bearing capabilities. Therefore, they are widely used in thermal protection structures for hypersonic vehicles. In this work, the aerothermoelastic behaviors of functionally graded (FG) plates under arbitrary temperature fields are analyzed via a semianalytical method. This research develops a method considering the influence of thermal loading, specifically the decrease in stiffness due to thermal stresses, as well as the correlation between material properties and temperatures under arbitrary temperature fields, based on Ritz’s method. The classical plate theory, von–Karman’s large defection plate theory and piston theory are employed to formulate the strain energy, kinetic energy and external work functions of the system. This paper presents a novel analysis of static aerothermoelasticity of FG plates, in addition to the linear/nonlinear flutter under arbitrary temperature fields, such as uniform, linear and nonlinear temperature fields. In addition, the effects of the volume fraction index, dynamic pressure, and temperature increase on the aerothermoelastic characteristics of FG plates are analyzed.


Introduction
Hypersonic vehicles, defined as vehicles with a flight Mach number greater than 5, have emerged as a crucial area of interest in aerospace applications due to their rapid flight speed, robust defense capability, and extended flight range.A substantial amount of aerodynamic heat is generated during flight because of the significant air friction.There is a complex coupling effect with the structure, resulting in aerothermoelastic behaviors.Notably, panel flutter occurs in supersonic vehicles and can drastically affect their safety and durability.Several scholars have conducted groundbreaking studies on panel flutter.In earlier researches, linear analysis was used to study panel flutter [1].Later, Dowell [2] conducted a study of two-and three-dimensional plates undergoing limit cycle oscillations using von Karman's large deflection plate theory and quasisteady aerodynamic theory.Additionally, the Galerkin method was used to solve discretely considered in-plane and static pressure differentials.The findings indicated that the limit cycle deflection, stress, and frequency are functions of the dynamic pressure, air/panel mass ratio, static pressure differential, in-plane load and length-to-width ratio.The findings indicated that four to six modes can yield accurate results.There are also numerous methods for the nonlinear analysis of panel flutter, such as the finite element frequency domain method proposed by Xue and Mei [3] and the aeroelastic modes proposed by Guo and Mei [4].McNamara and Friedmann concluded the aeroelastic and aerothermoelastic analyses of hypersonic flow in their review paper [5].Zhou [6] conducted an aerothermoelastic analysis of moderately thick orthotropic plates with general boundary conditions.Quan [7] conducted an analysis of the post-flutter aerothermoelastic characteristics of hypersonic skin panels using a CFDbased approach.
Due to severe aerodynamic heating, it is often crucial to apply a thermal protective layer on the surface to safeguard the structures and devices inside a vehicle.The thermal insulation layer and the structural layer consist of disparate materials with distinct mechanical characteristics, which can cause stress concentration, delamination, peeling, and other interfacial phenomena, impairing structural integrity and safety.In order to avoid this phenomenon, Japanese researchers proposed the concept of functionally graded materials (FGMs), which involves the continuous alteration of the material properties of two or more materials by varying their compositions and structures.This ensures that there is no discernible interface or structural mutation.This continuous variation in material properties promotes the use of FGMs as crucial elements of thermal protective systems for hypersonic vehicles, which has been extensively researched by scholars investigating the phenomenon of panel flutter.It is beneficial to investigate panel flutter for suppression of flutter and optimization.Song [8] moderated the mode of localized airflow in the direction of the flow, for supersonic airflows, by designing axial FGMs to alter the thickness and material properties according to a specific law.The introduction of material property changes led to a significant enhancement in comparison to solely adjusting the thickness, resulting in the passive suppression of flutter.Wei [9] examined the flutter characteristics of axial functionally graded column shells, which could be improved by adjusting the thickness and material properties in the axial direction.Muc and Flis [10] studied the flutter and free vibration characteristics of functionally graded porous plates using analytical and Rayleigh-Ritz methods.Optimization can effectively increase the natural frequency and critical dynamic pressure according to the outcome of a variable parameter analysis.In recent years, novel structures based on the concept of FGMs have been created and analyzed.These structures include functionally graded graphene nanoplatelet-reinforced composite (FG-GPLRC) cylindrical panels [11][12][13], cracked functionally graded plates [14], FG multilayered graphene platelet-reinforced polymer composites (GPL-RPC) [15,16], FG porous doubly curved panels with piezoelectric layers [17], imperfect FG cylindrical shells [18], magnetic-electric-thermoelastic functionally graded plates [19], carbon-based nanohybrid composite panels [20], FG variable thickness blades [21], and cylindrical sandwich panels with FG-saturated porous cores [22], etc.
Thermal loads significantly impact the structural integrity of hypersonic vehicles, so there have been numerous investigations into the panel flutter of FGM structures in thermal environments.Several plate theories are commonly used in FGM modeling, including classical plate theory (CPT), first-order shear deformation theory (FSDT) [23][24][25], higher-order shear deformation theory (HSDT) [26], and other novel shear deformation theories, such as trigonometric shear deformation theory [27].Equations governing sheardeformation theories of plates are typically more complicated and can take shear effects into consideration, which should be not neglected in thick plates [28].However, the shear deformation theories lead to more complicated and less intuitive results than CPT that can characterize physical laws of plates directly.In aircrafts, thin-wall structures are commonly used.For thin plates, shear effects are so tiny that CPT can satisfy the accuracy need through fewer computations.Different methods are used to analyze the panel flutter of FG plates, such as finite element method (FEM) [29,30], Rayleigh-Ritz methods [10], Galerkin method [31,32], meshless method [33] and generalized differential quadrature method (GDQM) [34,35].Praveen and Reddy [36] employed a finite element approach that considered shear strain to assess both the static and dynamic responses of ceramic-metal FG plates under different transient temperature fields.Ibrahim [37][38][39] employed a finite element technique based on first-order shear deformation theory and classical plate theory, respectively, to analyze the nonlinear flutter and thermal buckling issues of FG plates subjected to aerodynamic, thermal, and random acoustic loads.Navazi and Haddadpour [31,32] examined the impact of various parameters on the nonlinear aeroelastic behavior of semi-infinite FG plates using the Galerkin method.It was found that FG plates made of temperature-sensitive materials using the Mori-Tanaka model demonstrated less stability than those utilizing a simple power-law model.Su [40] proposed a unified method for the consideration of vibrations and flutter in elastically constrained FG plates under nonlinear temperature fields.The method involved the use of a modified Fourier series as a displacement basis function to express the displacement field of the reinforcement through the displacement variable of the plate, thereby eliminating the need for additional displacement variables.Lee and Kim [41] investigated the thermal postbuckling and snap-through instabilities of FG plates in hypersonic flows.
Most studies on panel flutter of FG plates are conducted in uniform or nonlinear temperature fields.However, the temperature field inside the structure of hypersonic vehicles can often become complex due to the influence of thermal loads.Therefore, it is necessary to establish an aerothermoelastic analysis method that can accommodate arbitrary temperature fields.Additionally, thermal stress and the relationship between properties and temperature should be taken into account simultaneously in thermal environments.
However, among the various aerothermoelastic behaviors, current studies have mostly concentrated on panel flutter for FG plates.There is a glaring lack of static aerothermoelastic studies examining aeroelastic responses in a certain temperature field in hypersonic flows.To investigate the influence of different variables on aerothermoelastic behaviors of FGM plates and provide a basis for design and optimization of FGM structures, concise solutions are needed.Therefore, aerothermoelastic analysis based on CPT and Ritz method is helpful to acquire a semianalytical solution for the FGM plate response.
In this study, an aerothermoelastic analysis was conducted via a semianalytical method.Static aerothermoelastic, linear panel flutter, and nonlinear panel flutter behaviors were imposed on simply supported FG plates.Equilibrium equations were established using classical plate theory and the Ritz method under arbitrary temperature fields in general forms, such as uniform, linear, and nonlinear temperature fields.Thermal stress and temperature-dependent properties were focused on in analytical expressions.The von Karman's strain-displacement relationship was used to consider geometric nonlinearity, and the stress function is used for simplification in case of jointly solved the multivariable partial differential equations.The effects of the volume fraction index, aerodynamic pressure, and thermal load were investigated.

Methods
In this section, the relevant theories including structure modeling of FGM, thermal effects, and aerodynamic modeling are briefly presented to serve for theoretical derivations.The semianalytical expressions of aerothermoelasic behaviors of FG plates are also presented.

FGM Properties
Unlike those of conventional homogeneous materials, the properties of FGMs vary because they consist of a blend of two components.For metal-ceramic FGMs, which are commonly used in hypersonic vehicles, modifications in properties are attained based on a continuous component gradient.
The Voigt model is used here to determine the effective material properties.The FGM properties can be represented by the following equation: where P(z), P c , and P m express the properties of the FGM, ceramic and metal, respectively.The subscripts c, m and cm denote the properties of the ceramic and metal and the difference between the ceramic and metal, that is, The FG plate coordinates are presented in Figure 1.The airflow is assumed to be in the x-direction, and the z-coordinate is measured from the mid-plane of the plates, where a, b, and h denote the length, the width, and the thickness of the FG plate, respectively.The relationship between the ceramic volume ratio, V c and the z-coordinate can be expressed as follows: where k is the volume fraction index of the FGM.where P(z), Pc, and Pm express the properties of the FGM, ceramic and metal, respectively.The subscripts c, m and cm denote the properties of the ceramic and metal and the difference between the ceramic and metal, that is, cm c m

P P P = −
(2) The FG plate coordinates are presented in Figure 1.The airflow is assumed to be in the -direction, and the -coordinate is measured from the mid-plane of the plates, where a, b, and h denote the length, the width, and the thickness of the FG plate, respectively.The relationship between the ceramic volume ratio, Vc and the z-coordinate can be expressed as follows: where k is the volume fraction index of the FGM.Therefore, for FG plates, the material properties vary according to the location in the thickness direction as follows: ( )

Constitutive Relations
This study focuses on the modeling and characterization of FG thin plates, which are predominantly utilized in aircrafts.Thin-walled structures with small shear effects are commonly used in aircraft, making the classical plate theory (CPT) an appropriate choice for analysis due to its speed and simplicity of application.Because of the asymmetry of the material in the thickness direction of FG plates, the strain on the geometric middle surface is not 0. Instead, the physical neutral surface has a strain of zero.When subjected to a bending load, only the axial force is applied on this plane, with no moment.The distance between the physical neutral surface and the geometric middle surface, z0 can be expressed as [42] ( ) Therefore, for FG plates, the material properties vary according to the location in the thickness direction as follows:

Constitutive Relations
This study focuses on the modeling and characterization of FG thin plates, which are predominantly utilized in aircrafts.Thin-walled structures with small shear effects are commonly used in aircraft, making the classical plate theory (CPT) an appropriate choice for analysis due to its speed and simplicity of application.Because of the asymmetry of the material in the thickness direction of FG plates, the strain on the geometric middle surface is not 0. Instead, the physical neutral surface has a strain of zero.When subjected to a bending load, only the axial force is applied on this plane, with no moment.The distance between the physical neutral surface and the geometric middle surface, z 0 can be expressed as [42] where E(z), E c , and E m are the elastic moduli of the FGM, the ceramic, and the metal, respectively.The strains on the FG plates can be expressed through the deflection function, w, as According to the plane stress assumption, the stress on the FG plates can be expressed as where υ is Poisson's ratio.

Thermal Effects on FGM
High-temperature environments are frequently encountered when using FGMs.Such environments subject structures to intense conditions, revealing structural flaws that might not surface at lower temperatures.The impacts of thermal environments on structures are typically apparent in two forms: (1) Thermal stresses related to temperature changes lower the overall structural stiffness.
(2) Material properties are subject to changes in response to temperature, often resulting in a decrease in overall material performance.
In the following section, the dual impacts of thermal environments on FGMs are outlined.
The thermal environment in hypersonic vehicles is intricate, necessitating the characterization of the temperature field through a general form.This would facilitate the analysis of thermal loads and aerothermoelastic behaviors for FG plates.The polynomial expansion form can fulfil the requirement of modeling FG plates with volume fraction as a power function, and also represent various types of temperature fields with satisfactory accuracy.Consequently, in this paper, the temperature field is characterized by means of a polynomial expansion.
The temperature increase in the FG plates can be expressed as where z is the nondimensional z-coordinate, which is defined as z = 1/2 + z/h, and c i is the coefficient of the temperature increase with respect to the nondimensional z-coordinate.
The thermal stresses can be expressed as where σ T , α, υ, and ∆T denote the thermal stress, coefficient of thermal expansion, Poisson's ratio and temperature increase caused by thermal loads, respectively.The effect of temperature loading on the structure can be considered by the in-plane compression force as where N T is the thermal force.
From an energy perspective, the work performed by thermal stress can be defined as where W T is the thermal work.
In addition, temperature influences FGMs by altering the material properties and can significantly affect the structure.This phenomenon can be described as the correlation between material properties and an increase in environment temperature, where P 0 , P −1 , P 1 , P 2 , and P 3 denote the coefficients of properties and the subscript i denotes the constituent material, such as 'm' for metal and 'c' for ceramic: To ensure clarity and conciseness, only the primary term impact of temperature on the material properties is explored: The material properties of Ti-6Al-4V and ZrO 2 are shown in Table 1.Changes in material properties primarily impact the elastic modulus and thermal conductivity, which significantly affects the structure.The correlation between material properties and temperature alters the strain energy and the work performed by thermal stresses.The general forms of the thermal forces and the position of the physical neutral surface at arbitrary temperature fields are listed in the Appendix A.
In this paper, the distribution of temperature fields within the FG plate can be designated as uniform, linear or nonlinear.By determining the temperature field within the plate, the thermal force can be obtained.
In the uniform temperature field, the temperature increase in the FG plate is uniform, where the corresponding coefficients in Equation ( 8) are When there is no internal heat source and the temperature increase (∆T) remains constant, known as steady-state heat transfer, Fourier's law dictates that where λ is the thermal conductivity.
When the thermal conductivity is assumed to be constant, the temperature distribution is linear through the thickness.The corresponding coefficients of the linear temperature field in Equation (8) become Due to the nonuniformity of the material properties in the thickness direction, the thermal conductivity is related to the spatial location in the FGM as follows: For FG plates, the differential equation Equation ( 15) is no longer linear, and the temperature field solved is also nonlinear.By using Taylor expansion, the nonlinear temperature field in the FG plates can be acquired.The corresponding coefficients in Equation ( 8) are The thermal forces and the position of the physical neutral surface under uniform, linear and nonlinear temperature fields can be acquired by substituting the coefficients above into the general equations in the Appendix A. Temperature and property distributions of FG plates under uniform, linear and nonlinear temperature fields can be found in the previous study [43].

Aerodynamic Modeling
In this paper, the aerodynamic forces are calculated using piston theory, which is analogous to the propagation of a perturbation generated by a wing in the normal direction at supersonic speeds.When Ma ≫ 1, the first-order piston theory can be expressed as where q is the applied aerodynamic pressure, q d is the dynamic pressure, Ma is the Mach number, α 0 is the initial angle of attack, c ∞ is the velocity of sound and t is the time.
The aerodynamic forces can be decomposed into three parts: first, a spatially relevant term concerning the normal disturbance of the incoming flow caused by the variations in the slope of the plate surface; second, a constant term related to the downwash produced by the initial angle of attack; and third, a time-varying term originating from the oscillations of the plate surface.
The displacement function can be selected as a double sine function under a simply supported boundary condition.To examine the work carried out by the theoretical aerodynamic force, we apply the first-order function in the vertical inflow direction, and two separate ordered functions as the m 1 -and m 2 -order functions in the downstream flow direction: The influence of the applied aerodynamic forces can be quantified by the work of the external force as

Nonlinear Characterization
Dowell [2] derived a form of geometrical nonlinearity in differential equations by employing a stress function.The impact of von Karman's nonlinear strain-displacement relation on the structure manifests in the integrated stress value throughout the plate thickness, influencing the neutral surface.Because of the variation in the elastic modulus through the thickness of FGMs, the internal force calculation is slightly different from that for homogenous plates.
As in Dowell's calculation in [2], the effect of the geometric nonlinear term can be expressed in energy form as shown in Equation (23), where A, B, C, D, E, F, and G are integration coefficients.The specific expression is shown in Equation (A5) in Appendix A, where E FGM is the equivalent elastic modulus of the FGM:

Static Aerothermoelastic Analysis
In this paper, the aerothermoelastic behaviors of the FG plate are determined with the Ritz method and classical plate theory, in which only one unknown variable, the deflection function, is necessary to express other physical quantities.The double sine function can be used as the basis function based on simply supported boundary conditions.The deflection function can be expressed as follows: where c mn is an unknown coefficient.The partial differential equations can be transformed into algebraic equations or ordinary differential equations via the Ritz method.By substituting the deflection function, the strain energy, thermal strain energy, and kinetic energy can be obtained.These quantities are determined under uniform temperature conditions and when the material properties are temperature independent as follows: where In aeroelastic analyses, it is essential to consider the impact of the work by the applied aerodynamic forces.In the static analysis, the time-dependent term in Equation ( 20) is not considered.When two separate ordered functions are selected, the deflection function can then be expressed as where c w1 and c w2 are unknown coefficients.
The general applied aerodynamic forces can be expressed as In this case, the constant component of the odd-order term is effective, whereas the even-order term is not.The aerodynamic forces of the same order of the spatially relevant term do not influence the deformations of the same order, and there is a coupling effect between the different orders.
From Hamilton's principle, the governing equation of static aerothermoelastic behaviors can be derived as Using two orders as an example, we can derive the algebraic equation as where when the temperature-independent FG plates are under a uniform temperature field, the corresponding coefficients can be expressed as shown in Equation (35).
The solution of the unknown coefficients is given by

Linear Analysis of the Thermal Panel Flutter
In linear analyses of thermal panel flutter, with the assumption of harmonic oscillations, the deflection function can be expressed as Then the kinetic energy can be expressed as Linear thermal panel flutter analysis can be formulated as a generalized eigenvalue problem to acquire the flutter modes and frequencies.In this analysis, only the homogeneous terms are retained, and the constant term is eliminated.Consequently, the generalized aerodynamic force can be expressed as From Hamilton's principle, the governing equation can be derived, as A generalized eigenvalue problem can be formulated as Ac = λBc with the corresponding eigenvalues λ m = ω 2 m , in which, ω m denotes the flutter frequency for the corresponding mode shape number m in the x direction.For temperature-independent FG plates under a uniform temperature field, this problem can be solved by expressing the i-th diagonal term of the coefficient matrix A as The nondiagonal term primarily comprises coupling terms of different orders of the generalized aerodynamic forces.The i-th row and j-th column can be expressed as Matrix B, which represents the coefficient of mass, consists solely of diagonal entries.The i-th element can be expressed as

Nonlinear Analysis for the Thermal Panel Flutter
In nonlinear analysis, the effect of geometric nonlinearity is considered.The aerodynamic forces consist of constant, time-dependent and spatially dependent terms and are explicitly defined as The deflection function can be defined as Then the generalized aerodynamic forces are expressed as From Hamilton's principle, the governing equation can be derived as When the temperature-independent FG plates are under a uniform temperature field, a set of M differential equations for deflection coefficients can be derived.The m-th differential equation is shown in Equation (48).where D FGM is the equivalent flexural rigidity:

Numerical Results of Static Aerothermoelasticity
In this section, vibration analysis and static aerothermoelastic analysis of simple supported FG thin plates are conducted under a uniform/linear/nonlinear temperature field.The impact of the thermal environment is particularly explored.
Case 1: A ceramic/metal FG plate composed of ZrO 2 /Ti-6Al-4V with a length of a = 0.5 m, a width of b = 0.5 m and a thickness of h = 0.05 m.The material properties of ZrO 2 and Ti-6Al-4V are shown in Table 1.
The following numerical results are based on the model in Case 1.The results are displayed as dimensionless variables, including the dimensionless time, dynamic pressure, and air-to-plate mass ratio, as shown in Equation (49).

Validation
To validate the FGM modeling method in this study, a vibration analysis is conducted for comparison with the FEM.The finite element model is partitioned into 50 sections in the length and width directions, and 10 segments in the thickness direction, resulting in 250,000 solid elements with the simple-supported boundary conditions, as shown in Figure 2. The first-order vibration frequencies, which correspond to the bending modes, are compared between the semianalytical method and the FEM results solved by sol 103 in MSC.Nastran, as shown in Figure 3. Firstly, the stiffness matrix under thermal loading is solved by sol 153.Based on the updated stiffness matrix as a static pre-stress field, the vibration characteristics are solved by sol 103. (

Numerical Results of Static Aerothermoelasticity
In this section, vibration analysis and static aerothermoelastic analysis of simple supported FG thin plates are conducted under a uniform/linear/nonlinear temperature field.The impact of the thermal environment is particularly explored.
Case 1: A ceramic/metal FG plate composed of ZrO2/Ti-6Al-4V with a length of a = 0.5 m, a width of b = 0.5 m and a thickness of h = 0.05 m.The material properties of ZrO2 and Ti-6Al-4V are shown in Table 1.
The following numerical results are based on the model in Case 1.The results are displayed as dimensionless variables, including the dimensionless time, dynamic pressure, and air-to-plate mass ratio, as shown in Equation (49).

Validation
To validate the FGM modeling method in this study, a vibration analysis is conducted for comparison with the FEM.The finite element model is partitioned into 50 sections in the length and width directions, and 10 segments in the thickness direction, resulting in 250,000 solid elements with the simple-supported boundary conditions, as shown in Figure 2. The first-order vibration frequencies, which correspond to the bending modes, are compared between the semianalytical method and the FEM results solved by sol 103 in MSC.Nastran, as shown in Figure 3. Firstly, the stiffness matrix under thermal loading is solved by sol 153.Based on the updated stiffness matrix as a static pre-stress field, the vibration characteristics are solved by sol 103.

Numerical Results of Static Aerothermoelasticity
In this section, vibration analysis and static aerothermoelastic analysis of simple supported FG thin plates are conducted under a uniform/linear/nonlinear temperature field.The impact of the thermal environment is particularly explored.
Case 1: A ceramic/metal FG plate composed of ZrO2/Ti-6Al-4V with a length of a = 0.5 m, a width of b = 0.5 m and a thickness of h = 0.05 m.The material properties of ZrO2 and Ti-6Al-4V are shown in Table 1.
The following numerical results are based on the model in Case 1.The results are displayed as dimensionless variables, including the dimensionless time, dynamic pressure, and air-to-plate mass ratio, as shown in Equation (49).

Validation
To validate the FGM modeling method in this study, a vibration analysis is conducted for comparison with the FEM.The finite element model is partitioned into 50 sections in the length and width directions, and 10 segments in the thickness direction, resulting in 250,000 solid elements with the simple-supported boundary conditions, as shown in Figure 2. The first-order vibration frequencies, which correspond to the bending modes, are compared between the semianalytical method and the FEM results solved by sol 103 in MSC.Nastran, as shown in Figure 3. Firstly, the stiffness matrix under thermal loading is solved by sol 153.Based on the updated stiffness matrix as a static pre-stress field, the vibration characteristics are solved by sol 103.Table 2 shows the dimensionless frequencies (ω 1 ) of the FG plates, defined in Equation (50) under different temperature fields obtained with the Ritz-based method and FEM at k = 1 and a temperature increase of 30 K, which means that ∆T = 30 K under a uniform temperature field and T m = 30 K and T cm = 30 K under a linear and nonlinear temperature fields.In addition, the results under room temperature without the influence of thermal gradient are listed as a reference.The relative errors between the semianalytical method and the finite element method are trivial, thus confirming the accuracy of the proposed approach.Errors arise from the use of the semianalytical method based on the CPT, which neglects the shear effect, while the FEM utilizes three-dimensional solid element modeling to account for this effect.In addition, the semianalytical method involves a continuous theoretical approach in the thickness direction whereby the material properties and individual response quantities vary continuously.Conversely, the FEM provides discrete results and hence does not accurately represent the material properties at that specific location. (50)

Static Aerothermoelastic Analysis
First, the static aerothermoelastic behaviors are analyzed under a uniform temperature field.The first six orders are considered in the Ritz function.The conditions are assumed to be as follows: the flight Mach number is Ma = 5, the dimensionless mass-to-Mach number ratio is µ/Ma = 0.01, and the dimensionless dynamic pressure is λ d = 200, the temperature increase is ∆T = 20 K, and the initial angle of attack is α 0 = 0.01.
When k = 1, Figure 4 reveal the deformation of the FG plate and the corresponding applied aerodynamic pressures, respectively.The deformation is symmetric perpendicularly to the flow direction and exhibits a half-period of the sinusoidal function, as the first-order deflection function is the only one applied.In the direction of the incoming flow, the first-order coefficient dominates and generates a slight forward tilt of the deformation due to aerodynamic forces, resembling a sinusoidal function.In terms of the aerodynamic distribution on the FG plates, the distribution in the x-direction follows a cosine function law since the spatial correlation term is proportional to the first-order derivative of the de-flection function.The front half of the FG plate is designated the positive pressure zone, while the rear half is the negative pressure zone, and the maximum positive pressure zone occurs at the front end.
Research has revealed that the FG plate response is temperature sensitive, particularly in terms of the variation in the response versus the volume fraction index.Figure 5 shows the variance patterns of the maximum deformation with respect to the volume fraction index at ∆T = 0, 2 K, 4 K, 6 K, 8 K and 20 K.As the temperature increases, the correlation between the maximum deformation and the volume fraction index shifts from a monotonically increasing trend to a monotonically decreasing trend, revealing a highly intricate pattern of change.This is because structural stiffness is associated with mechanical stiffness, which leads to a positive contribution, and thermal stiffness related to the temperature increase, which leads to a reduction in stiffness.As the volume fraction index increases, metals with a lower elastic modulus and coefficient of thermal expansion than those of ceramic become more dominant.Simultaneously, the absolute values of both the mechanical stiffness and the thermal stiffness decrease.As the temperature increases, the reduction influence of the thermal stiffness increases, while the influence of the mechanical stiffness remains constant.Therefore, at lower temperatures, the mechanical stiffness dominates, and the maximum deformation grows with increasing volume fraction index.On the other hand, at higher temperatures, the thermal stiffness takes over, resulting in a decrease in the maximum deformation.Due to the contrasting trends of mechanical stiffness and thermal stiffness with respect to the volume fraction index, which have varying rates of change, the trend with respect to the volume fraction index becomes more intricate when the temperature increases moderately.Research has revealed that the FG plate response is temperature sensitive, particularly in terms of the variation in the response versus the volume fraction index.Figure 5 shows the variance patterns of the maximum deformation with respect to the volume fraction index at ΔT = 0, 2 K, 4 K, 6 K, 8 K and 20 K.As the temperature increases, the correlation between the maximum deformation and the volume fraction index shifts from a monotonically increasing trend to a monotonically decreasing trend, revealing a highly intricate pattern of change.This is because structural stiffness is associated with mechanical stiffness, which leads to a positive contribution, and thermal stiffness related to the temperature increase, which leads to a reduction in stiffness.As the volume fraction index increases, metals with a lower elastic modulus and coefficient of thermal expansion than those of ceramic become more dominant.Simultaneously, the absolute values of both the mechanical stiffness and the thermal stiffness decrease.As the temperature increases, the reduction influence of the thermal stiffness increases, while the influence of the mechanical stiffness remains constant.Therefore, at lower temperatures, the mechanical stiffness dominates, and the maximum deformation grows with increasing volume fraction index.On the other hand, at higher temperatures, the thermal stiffness takes over, resulting in a decrease in the maximum deformation.Due to the contrasting trends of mechanical stiffness and thermal stiffness with respect to the volume fraction index, which have varying rates of change, the trend with respect to the volume fraction index becomes more intricate when the temperature increases moderately.Research has revealed that the FG plate response is temperature sensitive, particularly in terms of the variation in the response versus the volume fraction index.Figure 5 shows the variance patterns of the maximum deformation with respect to the volume fraction index at ΔT = 0, 2 K, 4 K, 6 K, 8 K and 20 K.As the temperature increases, the correlation between the maximum deformation and the volume fraction index shifts from a monotonically increasing trend to a monotonically decreasing trend, revealing a highly intricate pattern of change.This is because structural stiffness is associated with mechanical stiffness, which leads to a positive contribution, and thermal stiffness related to the temperature increase, which leads to a reduction in stiffness.As the volume fraction index increases, metals with a lower elastic modulus and coefficient of thermal expansion than those of ceramic become more dominant.Simultaneously, the absolute values of both the mechanical stiffness and the thermal stiffness decrease.As the temperature increases, the reduction influence of the thermal stiffness increases, while the influence of the mechanical stiffness remains constant.Therefore, at lower temperatures, the mechanical stiffness dominates, and the maximum deformation grows with increasing volume fraction index.On the other hand, at higher temperatures, the thermal stiffness takes over, resulting in a decrease in the maximum deformation.Due to the contrasting trends of mechanical stiffness and thermal stiffness with respect to the volume fraction index, which have varying rates of change, the trend with respect to the volume fraction index becomes more intricate when the temperature increases moderately.When the temperature increase is 20 K or 2 K, the deformation of the FG plate and the applied aerodynamic pressure in relation to the volume fraction index are illustrated in Figures 6 and 7.The volume fraction index predominantly affects the overall stiffness of the structure.However, the influence does not directly alter the distribution pattern of the aerodynamic force.Additionally, the deformation in the flow direction remains consistent with the prescribed rules.The deformation and aerodynamic force of the FG plate fall between those of the metal and ceramic plates.At a temperature of 20 K, as the volume fraction index increases, the decrease in thermal stiffness is greater than the increase in mechanical stiffness, resulting in an overall stiffer FG plate.This leads to reduced deformation, decreased aerodynamic force due to aeroelastic deformation, and, consequently, a decreased aerodynamic force.When the temperature is 2 K, the overall stiffness increases since the reduction in thermal stiffness is less than the increase in mechanical stiffness, in contrast to the trend observed at ∆T = 20 K.The deformation and aerodynamic force of the FG plate have direct relationships with the dynamic pressure.Figures 8 and 9 depict the relationships between the deformation, aerodynamic force and deflection coefficients at each order and the dynamic pressure.When the dynamic pressure increases from λ = 100 to λ = 1000, the maximum dimensionless deformation shifts from 0.321 to 0.638 as the position of appearance moves forward, from 0.420 to 0.215.Additionally, the peak of the aerodynamic force rapidly increases, appearing at the leading edge.The first-order coefficient initially increases rapidly, reaching a maximum value at a dimensionless pressure of 0.478 before gradually decreasing, while the remaining coefficients increase as the dynamic pressure increases.As the dynamic pressure increases, the energy from the first order slowly spreads to subsequent orders, causing the other order coefficients to gradually increase while the first order becomes less dominant.This results in a gradual change in the deformation and the subjected aerodynamic force pattern, with the first order initially dominating and then spreading to the second and third orders.
fraction index increases, the decrease in thermal stiffness is greater than the increase mechanical stiffness, resulting in an overall stiffer FG plate.This leads to reduced def mation, decreased aerodynamic force due to aeroelastic deformation, and, consequen a decreased aerodynamic force.When the temperature is 2 K, the overall stiffness creases since the reduction in thermal stiffness is less than the increase in mechanical st ness, in contrast to the trend observed at ΔT = 20 K.The deformation and aerodynam force of the FG plate have direct relationships with the dynamic pressure.Figures 8 an depict the relationships between the deformation, aerodynamic force and deflection co ficients at each order and the dynamic pressure.When the dynamic pressure increa from λ = 100 to λ = 1000, the maximum dimensionless deformation shifts from 0.321 0.638 as the position of appearance moves forward, from 0.420 to 0.215.Additionally, peak of the aerodynamic force rapidly increases, appearing at the leading edge.The fi order coefficient initially increases rapidly, reaching a maximum value at a dimensionl pressure of 0.478 before gradually decreasing, while the remaining coefficients increase the dynamic pressure increases.As the dynamic pressure increases, the energy from first order slowly spreads to subsequent orders, causing the other order coefficients gradually increase while the first order becomes less dominant.This results in a grad change in the deformation and the subjected aerodynamic force pattern, with the first der initially dominating and then spreading to the second and third orders.fall between those of the metal and ceramic plates.At a temperature of 20 K, as the volu fraction index increases, the decrease in thermal stiffness is greater than the increase mechanical stiffness, resulting in an overall stiffer FG plate.This leads to reduced def mation, decreased aerodynamic force due to aeroelastic deformation, and, consequen a decreased aerodynamic force.When the temperature is 2 K, the overall stiffness creases since the reduction in thermal stiffness is less than the increase in mechanical st ness, in contrast to the trend observed at ΔT = 20 K.The deformation and aerodynam force of the FG plate have direct relationships with the dynamic pressure.Figures 8 an depict the relationships between the deformation, aerodynamic force and deflection co ficients at each order and the dynamic pressure.When the dynamic pressure increa from λ = 100 to λ = 1000, the maximum dimensionless deformation shifts from 0.321 0.638 as the position of appearance moves forward, from 0.420 to 0.215.Additionally, peak of the aerodynamic force rapidly increases, appearing at the leading edge.The fi order coefficient initially increases rapidly, reaching a maximum value at a dimensionl pressure of 0.478 before gradually decreasing, while the remaining coefficients increase the dynamic pressure increases.As the dynamic pressure increases, the energy from first order slowly spreads to subsequent orders, causing the other order coefficients gradually increase while the first order becomes less dominant.This results in a grad change in the deformation and the subjected aerodynamic force pattern, with the first der initially dominating and then spreading to the second and third orders.The maximum deformation and maximum stress of the FG plate under varying temperature fields, with a volume fraction index of k = 1 and temperature increase of 20 K and 2 K, are outlined in Table 3.The results show that the maximum deformation and maximum stress are the largest under the uniform temperature field, but are the smallest under the nonlinear temperature field.This is because at the same temperature increase, the FG plate experiences the highest internal temperature, the highest thermal stress, the lowest total stiffness, and the largest deformation and stress within the uniform temperature field.Conversely, this effect is associated with the lowest internal temperature and the smallest deformation and stress within the nonlinear temperature field.Considering the temperature-dependent properties, when subjected to a larger temperature increase of 20 K, the coefficient of thermal expansion of the ceramic component decreases more, resulting in less thermal stiffness but higher total stiffness.Conversely, a smaller temperature increase of 2 K causes the coefficient of thermal expansion and elastic modulus of the metal component to increase further, leading to decreased thermal stiffness and lower total stiffness.The maximum deformation and maximum stress of the FG plate under varying temperature fields, with a volume fraction index of k = 1 and temperature increase of 20 K and 2 K, are outlined in Table 3.The results show that the maximum deformation and maximum stress are the largest under the uniform temperature field, but are the smallest under the nonlinear temperature field.This is because at the same temperature increase, the FG plate experiences the highest internal temperature, the highest thermal stress, the lowest total stiffness, and the largest deformation and stress within the uniform temperature field.Conversely, this effect is associated with the lowest internal temperature and the smallest deformation and stress within the nonlinear temperature field.Considering the temperature-dependent properties, when subjected to a larger temperature increase of 20 K, the coefficient of thermal expansion of the ceramic component decreases more, resulting in less thermal stiffness but higher total stiffness.Conversely, a smaller temperature increase of 2 K causes the coefficient of thermal expansion and elastic modulus of the metal component to increase further, leading to decreased thermal stiffness and lower total stiffness.

Research on Aeroelastic Effects
Furthermore, this study examines the impact of aeroelastic effects.The investigation is conducted in the absence of aeroelastic effects, i.e., when the aerodynamic force contains only the constant term originating from the angle of attack and not the component caused by deformation.Figures 10-12 display the deformation of the structure and its corresponding variations with the volume fraction index and dynamic pressure.The deformation approximates a sinusoidal function and is symmetrically distributed when the aeroelastic effect is not considered due to the constant aerodynamic force on the FG plate.As the dynamic pressure increases, the stringwise deformation increases, while the deformation distribution remains symmetric.With a temperature increase of ∆T = 20 K, the maximum deformation decreases with increasing volume fraction index, but still exceeds that considering aeroelastic effects.Therefore, the aeroelastic effect primarily reduces structural deformation and alters its form to direct the deformation toward the leading edge.In addition, deformations display variations under distinct dynamic pressures, with increased pressures moving the deformation increasingly forward.mation distribution remains symmetric.With a temperature increase of ΔT = 20 K, the maximum deformation decreases with increasing volume fraction index, but still exceeds that considering aeroelastic effects.Therefore, the aeroelastic effect primarily reduces structural deformation and alters its form to direct the deformation toward the leading edge.In addition, deformations display variations under distinct dynamic pressures, with increased pressures moving the deformation increasingly forward.As the dynamic pressure increases, the stringwise deformation increases, while the deformation distribution remains symmetric.With a temperature increase of ΔT = 20 K, the maximum deformation decreases with increasing volume fraction index, but still exceeds that considering aeroelastic effects.Therefore, the aeroelastic effect primarily reduces structural deformation and alters its form to direct the deformation toward the leading edge.In addition, deformations display variations under distinct dynamic pressures, with increased pressures moving the deformation increasingly forward.

Numerical Results of Panel Flutter in Thermal Environments
In this section, linear and nonlinear analyses of the panel flutter of FG plates in thermal environments are conducted under a uniform/linear/nonlinear temperature field based on Case 1.

Linear Analysis
First, the results of the Ritz-based method in this study are compared with those of Song's method [44] for a metal plate.The relationship between the natural frequencies and dynamic pressure of the plates when k = 0 in Case 1 is depicted in Figure 13.It can be seen that the curves are very close to each other, which verifies the correctness of the method used in this paper.

Numerical Results of Panel Flutter in Thermal Environments
In this section, linear and nonlinear analyses of the panel flutter of FG plates in thermal environments are conducted under a uniform/linear/nonlinear temperature field based on Case 1.

Linear Analysis
First, the results of the Ritz-based method in this study are compared with those of Song's method [44] for a metal plate.The relationship between the natural frequencies and dynamic pressure of the plates when k = 0 in Case 1 is depicted in Figure 13.It can be seen that the curves are very close to each other, which verifies the correctness of the method used in this paper.

Linear Analysis
First, the results of the Ritz-based method in this study are compared with those of Song's method [44] for a metal plate.The relationship between the natural frequencies and dynamic pressure of the plates when k = 0 in Case 1 is depicted in Figure 13.It can be seen that the curves are very close to each other, which verifies the correctness of the method used in this paper.In this section, linear analysis is conducted for the FG plate in Case 1 using the first six Ritz orders.Figure 14a illustrates the relationship between the dynamic pressure and flutter frequency for each order at a Mach number of Ma = 5, a temperature increase of ΔT = 20 K of each order experiences different changes.Low-order frequency changes are more pronounced, while the first-order frequency gradually increases, and the high-order frequency slightly decreases.At a dimensionless dynamic pressure of 424, the first two orders of frequency coincide, and the eigenvalues of the first two orders of the system transform from two real roots into a pair of conjugate complex roots, thus inducing panel flutter.This is due to the mutual coupling between different orders of aerodynamic forces in the spatial correlation term in piston theory.This coupling results in the higher-order term providing additional aerodynamic stiffness for the lower-order term, improving the lower-order natural frequencies.Similarly, the lower-order term weakens the aerodynamic stiffness of the higher-order term.As the dynamic pressure increases, the frequency of the first coupled orders becomes the lowest-order term and increases with increasing dynamic pressure, while the higher-order dynamic pressure decreases with increasing In this section, linear analysis is conducted for the FG plate in Case 1 using the first six Ritz orders.Figure 14a illustrates the relationship between the dynamic pressure and flutter frequency for each order at a Mach number of Ma = 5, a temperature increase of ∆T = 20 K of each order experiences different changes.Low-order frequency changes are more pronounced, while the first-order frequency gradually increases, and the high-order frequency slightly decreases.At a dimensionless dynamic pressure of 424, the first two orders of frequency coincide, and the eigenvalues of the first two orders of the system transform from two real roots into a pair of conjugate complex roots, thus inducing panel flutter.This is due to the mutual coupling between different orders of aerodynamic forces in the spatial correlation term in piston theory.This coupling results in the higher-order term providing additional aerodynamic stiffness for the lower-order term, improving the lower-order natural frequencies.Similarly, the lower-order term weakens the aerodynamic stiffness of the higher-order term.As the dynamic pressure increases, the frequency of the first coupled orders becomes the lowest-order term and increases with increasing dynamic pressure, while the higher-order dynamic pressure decreases with increasing dynamic pressure.The nondimensional frequency can be defined as the ratio of the frequency to the first-order natural frequency of ceramic plates, ω 0 , as Aerospace 2024, 11, x FOR PEER REVIEW 19 of 29 dynamic pressure.The nondimensional frequency can be defined as the ratio of the frequency to the first-order natural frequency of ceramic plates, ω0, as ( ) To examine the flutter characteristics at various temperatures, Figure 14b shows the fluctuations in the frequencies in the first two orders, corresponding to dynamic pressure increases at different temperatures.Furthermore, higher temperatures are shown to have a more notable effect in this regard.The results indicate that as the temperature increases, the overall stiffness of the FG plate decreases, causing a reduction in the frequency and flutter dynamic pressure.Due to the presence of aerodynamic forces, the first-order stiffness increases, allowing the FG plate to regain a stable state within a certain range even after buckling under no aerodynamic loading.In order to investigate the influence of the Mach number, Figure 15 show the variation of flutter frequencies with the nondimensional dynamic pressure and the actual dynamic pressure under various Mach numbers, respectively.It can be observed that the flutter characteristics are independent of the dimensionless dynamic pressure, as the Mach number is already accounted for in the nondimensionalization, which is also shown in Equation (49).In contrast, the flutter characteristics are dependent of the actual dynamic pressure.The flutter dynamic pressure is higher at the higher Mach number when the dynamic pressure is dimensional.To examine the flutter characteristics at various temperatures, Figure 14b shows the fluctuations in the frequencies in the first two orders, corresponding to dynamic pressure increases at different temperatures.Furthermore, higher temperatures are shown to have a more notable effect in this regard.The results indicate that as the temperature increases, the overall stiffness of the FG plate decreases, causing a reduction in the frequency and flutter dynamic pressure.Due to the presence of aerodynamic forces, the first-order stiffness increases, allowing the FG plate to regain a stable state within a certain range even after buckling under no aerodynamic loading.
In order to investigate the influence of the Mach number, Figure 15 show the variation of flutter frequencies with the nondimensional dynamic pressure and the actual dynamic pressure under various Mach numbers, respectively.It can be observed that the flutter characteristics are independent of the dimensionless dynamic pressure, as the Mach number is already accounted for in the nondimensionalization, which is also shown in Equation (49).In contrast, the flutter characteristics are dependent of the actual dynamic pressure.The flutter dynamic pressure is higher at the higher Mach number when the dynamic pressure is dimensional.In order to investigate the influence of the Mach number, Figure 15 show the variation of flutter frequencies with the nondimensional dynamic pressure and the actual dynamic pressure under various Mach numbers, respectively.It can be observed that the flutter characteristics are independent of the dimensionless dynamic pressure, as the Mach number is already accounted for in the nondimensionalization, which is also shown in Equation (49).In contrast, the flutter characteristics are dependent of the actual dynamic pressure.The flutter dynamic pressure is higher at the higher Mach number when the dynamic pressure is dimensional.Table 4 shows the flutter dynamic pressure under various temperature fields.With the same temperature increase, the internal temperature of the uniform temperature field plate is the highest, while the internal temperature of the nonlinear temperature field plate is the lowest and the opposite is true for the flutter dynamic pressure.Additionally, considering the temperature-dependent properties of FG plates, the total stiffness is lower, leading to larger deformations and a lower critical flutter dynamic pressure at lower temperatures.Figure 16 displays the modes of each order for the FG plate at dynamic pressures of λ = 200 (before the flutter boundary) and λ = 424 (near the flutter boundary).The aerodynamic forces lead to distinct modes at varying dynamic pressures.When the dynamic pressure does not reach the flutter boundary, each mode is directly influenced by its respective sinusoidal order, and the other sinusoidal functions have a lesser impact.However, when the dynamic pressure approaches the flutter boundary, the modes of each order exhibit changes compared to when the plate is not near the boundary.Specifically, the first few orders experience more dramatic changes, while the first two modes display almost identical shapes, indicating that they are coupled.This coupling results in the occurrence of flutter.
spective sinusoidal order, and the other sinusoidal functions have a lesser impact.How ever, when the dynamic pressure approaches the flutter boundary, the modes of each o der exhibit changes compared to when the plate is not near the boundary.Specifically, th first few orders experience more dramatic changes, while the first two modes display a most identical shapes, indicating that they are coupled.This coupling results in the occu rence of flutter.

Nonlinear Analysis
First, the results of the Ritz-based method in this study are compared with those of Dowell's method [2] at a dimensionless mass-to-Mach number ratio of µ/Ma = 0.01 and Mach number of Ma = 5.The relationship between the limit-cycle amplitudes and dynamic pressure of the ceramic plates when k = 0 in Case 1 is depicted in Figure 17.The results of this study resemble those of Dowell, indicating the accuracy of the method used in this study.
pressure of the ceramic plates when k = 0 in Case 1 is depicted in Figure 17.The results of this study resemble those of Dowell, indicating the accuracy of the method used in this study.
A nonlinear analysis of panel flutter in thermal environments is conducted for the FG plate in Case 1 at a dimensionless mass-to-Mach number ratio of µ/Ma = 0.01, Mach number of Ma = 5 and angle of attack of α0 = 0.01.With the initial condition that all order deflection coefficients have displacement and velocity values of 0, the differential equations described in Equation ( 48) can be numerically integrated in the time domain to determine various types of response quantities.A nonlinear analysis of panel flutter in thermal environments is conducted for the FG plate in Case 1 at a dimensionless mass-to-Mach number ratio of µ/Ma = 0.01, Mach number of Ma = 5 and angle of attack of α 0 = 0.01.With the initial condition that all order deflection coefficients have displacement and velocity values of 0, the differential equations described in Equation ( 48) can be numerically integrated in the time domain to determine various types of response quantities.
The dynamic responses and phase charts of the FG plate under different dynamic pressures at representative coordinates chosen as x/a = 0.75 and y/b = 0.5 are shown in Figure 18.The time-domain response characteristics vary significantly under various dynamic pressures.When the dynamic pressure is low, the FG plate remains dynamically stable.Initially, the perturbations decrease gradually while the system dissipates the perturbation energy, yielding a finite value.The phase diagram shows the convergence of the displacement and velocity of the representative coordinates toward a single point.As the dynamic pressure increases, the initial perturbation becomes increasingly larger.However, the system is unable to dissipate the initial perturbation energy.Consequently, the perturbation gradually accumulates and increases, resulting in the displacement and velocity of the representative coordinates converging to those of the limit cycle at λ d = 426.5.The system response then moves slowly from the center of the phase diagram to the outside and ultimately reaches the limit cycle.As the pressure increases, the amplitude of the limit cycle also increases because of the geometrically nonlinear effect of the structure.As the deformation becomes larger, the stiffness increases, which restricts further amplitude expansion.The system begins with a slight amplitude vibration near the equilibrium position.The system response swiftly expands from the center of the phase diagram to the limit ring and ultimately converges.
The limit-cycle amplitude for the dynamic pressure is illustrated in Figure 19a.When the critical dynamic pressure is not reached, the maximum deformation is minor, resulting solely from static forces.However, when the critical dynamic pressure is achieved, nonlinear panel flutter occurs due to the mutual influence of aerodynamic forces and the structure, and the limit-cycle amplitude grows rapidly.The relationship between the stress and dynamic pressure is illustrated in Figure 19b.The stresses can be divided into three parts: bending stresses, which are symmetrically related to the neutral surface; membrane stresses, which act on the neutral surface due to geometrical nonlinearities; and thermal stresses due to temperature loading.In bending stresses, the maximum tensile stresses are reached at the top while the minimum compressive stresses are reached at the bottom.On the other hand, in a uniform temperature field, the tensile membrane stresses and compressive thermal stresses are uniform in the FG plates.Both bending and membrane stresses are influenced by deformation, and the geometrically nonlinear effect becomes more prominent in large deformation cases.Prior to reaching the critical dynamic pressure, the stresses remain low due to minimal deformation, which is primarily controlled by thermal stresses.The maximum stresses at the top and bottom of the plate are similar and are both compressive.Once the flutter dynamic pressure is reached, rapid deformation and membrane stress expansion occur.Furthermore, due to the upward direction of the aerodynamic force in most cases, the maximum upward deformation is greater than the maximum downward deformation.Consequently, the absolute value of the maximum tensile stress exceeds that of the compressive stress, with a ratio of approximately 3:1.Therefore, the location of the maximum stress is at x/a = 3/4, y/b = 1/2, and z/h = 1/2.The stress is expressed in a dimensionless form, where The dynamic responses and phase charts of the FG plate under different dynamic pressures at representative coordinates chosen as x/a = 0.75 and y/b = 0.5 are shown in Figure 18.The time-domain response characteristics vary significantly under various dynamic pressures.When the dynamic pressure is low, the FG plate remains dynamically stable.Initially, the perturbations decrease gradually while the system dissipates the perturbation energy, yielding a finite value.The phase diagram shows the convergence of the displacement and velocity of the representative coordinates toward a single point.As the dynamic pressure increases, the initial perturbation becomes increasingly larger.However, the system is unable to dissipate the initial perturbation energy.Consequently, the perturbation gradually accumulates and increases, resulting in the displacement and velocity of the representative coordinates converging to those of the limit cycle at λd = 426.5.The system response then moves slowly from the center of the phase diagram to the outside and ultimately reaches the limit cycle.As the pressure increases, the amplitude of the limit cycle also increases because of the geometrically nonlinear effect of the structure.As the deformation becomes larger, the stiffness increases, which restricts further amplitude expansion.The system begins with a slight amplitude vibration near the To investigate the relationships between the limit-cycle amplitude and the dynamic pressure, volume fraction index, and temperature increase, the characteristics of nonlinear panel flutter are analyzed based on the basic state at λ = 427, k = 1 and ∆T = 20 K, which is near the critical conditions, as shown in Figure 20.Before reaching the critical dynamic pressure, the FG plate remains dynamically stable with minimal deformation due to static aerodynamic forces.However, after the critical dynamic pressure is reached, the limit-cycle amplitude expands rapidly.The amplitude is minimal near the leading edge, where it is downward, and it increases steadily toward the trailing edge.This is due to the increase in the limit-cycle amplitude, which causes the deformation to move toward the center.As the mechanical properties of the metal and ceramic components are quite similar in this example, the volume fraction index has less influence.However, as the volume fraction index increases and the proportion of the metal component increases, the total stiffness of the FG plate also increases, leading to an increase in the flutter dynamic pressure.The response of the FG plate is not always situated in the center between the metal plates and ceramic plates, such that the maximum deformation of the FG plate at k = 0.2 is slightly greater than that of the ceramic plates.This is because the responses of the plates are impacted by both the mechanical and thermal stiffness terms, which are related to the elastic modulus and thermal strains, respectively.Furthermore, at different volume fraction indexes, the various proportions of ceramic and metal components have different properties following a certain temperature increase, resulting in complex stiffness variations.At low temperatures and given dynamic pressures, the critical dynamic pressure is not reached.However, with an increase in temperature, the stiffness of the FG plate decreases, leading to deteriorated stability.Gradually, the plate becomes dynamically unstable, and the limit-cycle amplitude increases progressively.
linear panel flutter occurs due to the mutual influence of aerodynamic forces and the structure, and the limit-cycle amplitude grows rapidly.The relationship between the stress and dynamic pressure is illustrated in Figure 19b.The stresses can be divided into three parts: bending stresses, which are symmetrically related to the neutral surface; membrane stresses, which act on the neutral surface due to geometrical nonlinearities; and thermal stresses due to temperature loading.In bending stresses, the maximum tensile stresses are reached at the top while the minimum compressive stresses are reached at the bottom.On the other hand, in a uniform temperature field, the tensile membrane stresses and compressive thermal stresses are uniform in the FG plates.Both bending and membrane stresses are influenced by deformation, and the geometrically nonlinear effect becomes more prominent in large deformation cases.Prior to reaching the critical dynamic pressure, the stresses remain low due to minimal deformation, which is primarily controlled by thermal stresses.The maximum stresses at the top and bottom of the plate are similar and are both compressive.Once the flutter dynamic pressure is reached, rapid deformation and membrane stress expansion occur.Furthermore, due to the upward direction of the aerodynamic force in most cases, the maximum upward deformation is greater than the maximum downward deformation.Consequently, the absolute value of the maximum tensile stress exceeds that of the compressive stress, with a ratio of approximately 3:1.Therefore, the location of the maximum stress is at x/a = 3/4, y/b = 1/2, and z/h = 1/2.The stress is expressed in a dimensionless form, where To investigate the relationships between the limit-cycle amplitude and the dynamic pressure, volume fraction index, and temperature increase, the characteristics of nonlinear panel flutter are analyzed based on the basic state at λ = 427, k = 1 and ΔT = 20 K, which is near the critical conditions, as shown in Figure 20.Before reaching the critical dynamic pressure, the FG plate remains dynamically stable with minimal deformation due to static aerodynamic forces.However, after the critical dynamic pressure is reached, the limitcycle amplitude expands rapidly.The amplitude is minimal near the leading edge, where it is downward, and it increases steadily toward the trailing edge.This is due to the increase in the limit-cycle amplitude, which causes the deformation to move toward the center.As the mechanical properties of the metal and ceramic components are quite similar in this example, the volume fraction index has less influence.However, as the volume fraction index increases and the proportion of the metal component increases, the total stiffness of the FG plate also increases, leading to an increase in the flutter dynamic pressure.The response of the FG plate is not always situated in the center between the metal plates and ceramic plates, such that the maximum deformation of the FG plate at k = 0.2 is slightly greater than that of the ceramic plates.This is because the responses of the plates are impacted by both the mechanical and thermal stiffness terms, which are related to the elastic modulus and thermal strains, respectively.Furthermore, at different volume fraction indexes, the various proportions of ceramic and metal components have different properties following a certain temperature increase, resulting in complex stiffness variations.At low temperatures and given dynamic pressures, the critical dynamic pressure is not reached.However, with an increase in temperature, the stiffness of the FG plate decreases, leading to deteriorated stability.Gradually, the plate becomes dynamically unstable, and the limit-cycle amplitude increases progressively.

Discussion
In this study, a semianalytical aerothermoelastic analysis based on the energy method is conducted, encompassing static aerothermoelastic, linear and nonlinear analyses of the panel flutter of functionally graded (FG) plates in thermal environments.The influence of the thermal environment is focused on the thermal stiffness and thermal degradation under different temperature fields.In addition, the impacts of the volume fraction index, temperature increase, and dynamic pressure on the aerothermoelastic behaviors of the FG plates are also investigated.The method can theoretically aid in investigating the influences of different variables and in the efficient analysis and optimization of FG structures.

Discussion
In this study, a semianalytical aerothermoelastic analysis based on the energy method is conducted, encompassing static aerothermoelastic, linear and nonlinear analyses of the panel flutter of functionally graded (FG) plates in thermal environments.The influence of the thermal environment is focused on the thermal stiffness and thermal degradation under different temperature fields.In addition, the impacts of the volume fraction index, temperature increase, and dynamic pressure on the aerothermoelastic behaviors of the FG plates are also investigated.The method can theoretically aid in investigating the influences of different variables and in the efficient analysis and optimization of FG structures.
The results obtained using the semianalytical method in this paper are compared to those obtained using the finite element method and with the results reported in the literature.The comparison reveals a high degree of agreement, which provides evidence of the effectiveness of the method presented in this paper.However, due to the limitations of the conditions under which the work was conducted, the results have not been verified experimentally.Consequently, further verification of the method can be achieved through the ground vibration test of a FG plate and other methods.

Conclusions
(1) The paper presents a semianalytical analysis method for solving the unknown deflection coefficients to acquire the explicit solutions for static areothermoelastic responses and to solve the eigenvalues for the coefficient matrix and a system of ordinary differential equations in linear and nonlinear analysis for panel flutter, respectively.(2) According to the static aerothermoelastic analysis, the aeroelastic effect reduces the maximum deformation and the peak of the applied aerodynamic force.The coupling effect between the orders of aerodynamic forces shifts the deformation of the structure toward the leading edge and reduces the maximum deformation.(3) FG plates exhibit sensitivity to thermal loading, particularly concerning the correlation between stiffness and the volume fraction index.It is important to note that the ratio of mechanical stiffness to thermal stiffness was not consistent at different temperatures, as both of these stiffnesses progressively decrease with an increase in the volume fraction index despite contributing to the total stiffness in opposite directions.At lower temperatures, mechanical stiffness prevails, resulting in a decrease in total stiffness with an increase in the volume fraction index.(4) The general form of arbitrary temperature fields proposed in this paper can deal with complex thermal loads and provides a foundation for practical aerothermoelastic engineering problems.With the same temperature increase, the internal temperature of the uniform temperature field plate is the highest, while the internal temperature of the nonlinear temperature field plate is the lowest.
When the properties of FGM are temperature-dependent, the thermal forces become more complex in form, as shown in Equation (A2) where T ij (i = 1, 2, 3, j = 0, 1, 2) denotes the intermediate coefficients, as shown in Equation (A3).Moreover, the position of the physical neutral surface differs from the one not concerning the temperature-dependent elastic modulus of FGM, which can be expressed as

Figure 1 .
Figure 1.Coordinate system and geometry of the FG plates.

Figure 1 .
Figure 1.Coordinate system and geometry of the FG plates.

Figure 2 .
Figure 2. Finite element model of the FG plates in Case 1.

Figure 3 .
Figure 3.The first-order mode of the FG plates in Case 1 solved by FEM.

Figure 2 .
Figure 2. Finite element model of the FG plates in Case 1.

Figure 2 .
Figure 2. Finite element model of the FG plates in Case 1.

Figure 3 .
Figure 3.The first-order mode of the FG plates in Case 1 solved by FEM.Figure 3. The first-order mode of the FG plates in Case 1 solved by FEM.

Figure 3 .
Figure 3.The first-order mode of the FG plates in Case 1 solved by FEM.Figure 3. The first-order mode of the FG plates in Case 1 solved by FEM.

Figure 4 .
Figure 4.The static aerothermoelstic responses of FG plates under a uniform field: (a) static aerothermoelastic deformation; (b) aerodynamic pressure.

Figure 5 .
Figure 5. Maximum deformation of the FG plates vs. the volume fraction index at different tem-

Figure 4 .
Figure 4.The static aerothermoelstic responses of FG plates under a uniform field: (a) static aerothermoelastic deformation; (b) aerodynamic pressure.

Figure 4 .
Figure 4.The static aerothermoelstic responses of FG plates under a uniform field: (a) static aerothermoelastic deformation; (b) aerodynamic pressure.

Figure 5 .
Figure 5. Maximum deformation of the FG plates vs. the volume fraction index at different temperatures.

Figure 5 .
Figure 5. Maximum deformation of the FG plates vs. the volume fraction index at different temperatures.

Figure 8 .
Figure 8. Static aerothermoelastic responses of the FG plate under different dynamic pressures: (a) dimensionless deformation; (b) applied aerodynamic pressure.

Figure 8 .
Figure 8. Static aerothermoelastic responses of the FG plate under different dynamic pressures: (a) dimensionless deformation; (b) applied aerodynamic pressure.

Figure 10 .
Figure 10.Deformation of an FG plate without aeroelastic effects.

Figure 11 .
Figure 11.Dimensionless deformation of an FG plate under different dynamic pressures without aeroelastic effects.

Figure 10 .
Figure 10.Deformation of an FG plate without aeroelastic effects.

Figure 10 .
Figure 10.Deformation of an FG plate without aeroelastic effects.

Figure 11 .
Figure 11.Dimensionless deformation of an FG plate under different dynamic pressures without aeroelastic effects.

Figure 11 .
Figure 11.Dimensionless deformation of an FG plate under different dynamic pressures without aeroelastic effects.Aerospace 2024, 11, x FOR PEER REVIEW 18 of 29

Figure 12 .
Figure 12.Comparison of maximum deformation vs. the volume fraction index with and without aeroelastic effects.

Figure 12 .
Figure 12.Comparison of maximum deformation vs. the volume fraction index with and without aeroelastic effects.

Figure 13 .
Figure 13.Comparison of natural frequencies between the Ritz-based method in this study and Song's method.

Figure 13 .
Figure 13.Comparison of natural frequencies between the Ritz-based method in this study and Song's method.

Figure 14 .
Figure 14.Flutter frequencies vs. nondimensional dynamic pressure: (a) the first six orders at ΔT = 20 K; (b) the first two orders at different temperature increases.

Figure 14 .
Figure 14.Flutter frequencies vs. nondimensional dynamic pressure: (a) the first six orders at ∆T = 20 K; (b) the first two orders at different temperature increases.

Figure 14 .
Figure 14.Flutter frequencies vs. nondimensional dynamic pressure: (a) the first six orders at ΔT = 20 K; (b) the first two orders at different temperature increases.

Figure 15 .
Figure 15.The variation of flutter frequency under different Mach number: (a) with the nondimensional dynamic pressure; (b) with the dynamic pressure.Figure 15.The variation of flutter frequency under different Mach number: (a) with the nondimensional dynamic pressure; (b) with the dynamic pressure.

Figure 15 .
Figure 15.The variation of flutter frequency under different Mach number: (a) with the nondimensional dynamic pressure; (b) with the dynamic pressure.Figure 15.The variation of flutter frequency under different Mach number: (a) with the nondimensional dynamic pressure; (b) with the dynamic pressure.

Figure 17 .
Figure 17.Comparison of limit-cycle amplitudes between the Ritz-based method in this study and Dowell's method.Figure 17.Comparison of limit-cycle amplitudes between the Ritz-based method in this study and Dowell's method.

Figure 17 .
Figure 17.Comparison of limit-cycle amplitudes between the Ritz-based method in this study and Dowell's method.Figure 17.Comparison of limit-cycle amplitudes between the Ritz-based method in this study and Dowell's method.

Figure 18 .
Figure 18.Dynamic responses and phase diagrams of the FG plate at representative coordinates.

Figure 18 .
Figure 18.Dynamic responses and phase diagrams of the FG plate at representative coordinates.

Figure 19 .
Figure 19.Different variables of the FG plate at representative coordinates with nondimensional dynamic pressure: (a) nondimensional limit-cycle amplitude; (b) nondimensional stresses.

Figure 19 .
Figure 19.Different variables of the FG plate at representative coordinates with nondimensional dynamic pressure: (a) nondimensional limit-cycle amplitude; (b) nondimensional stresses.

Figure 20 .
Figure 20.Deformations of the FG plate at representative coordinates under different cases: (a) different nondimensional dynamic pressure; (b) different volume fraction indexes; (c) different temperature increases.

Figure 20 .
Figure 20.Deformations of the FG plate at representative coordinates under different cases: (a) different nondimensional dynamic pressure; (b) different volume fraction indexes; (c) different temperature increases.

Table 1 .
Material properties of FGM.

Table 2 .
Comparison of dimensionless first-order frequencies (ω 1 ) under different temperature fields determined by the semianalytical method and FEM.

Table 3 .
Static aerothermoelastic responses under various temperature fields.

Table 3 .
Static aerothermoelastic responses under various temperature fields.

Table 4 .
Flutter dynamic pressure under various temperature fields.

Table 4 .
Flutter dynamic pressure under various temperature fields.