Aerothermoelastic Analysis of Functionally Graded Plates Using Generalized Differential Quadrature Method

In the present paper, the aerothermoelastic behavior of Functionally Graded (FG) plates under supersonic airflow is investigated using Generalized Differential Quadrature Method (GDQM). The structural model is considered based on the classical plate theory and the von Karman strain-displacement relations are utilized to involve the nonlinear behavior of the plate. To consider the supersonic aerodynamic loads on the plate, the first order piston theory is applied. The material properties of the FG panel are assumed to be temperature independent and alter in the thickness direction according to a power law distribution. The temperature distribution on the surface of the plate is assumed to be constant and in the thickness direction is obtained by one-dimensional steady conductive heat transfer equation. The discretized governing equations via GDQM are solved by the fourth order Runge-Kutta method. Comparison of the obtained results with those available in literature confirms the accuracy and ability of the GDQM to perform the aerothermoelastic analysis of FG plates. Also, the effect of some important parameters such as Mach number, inplane thermal load, plate aspect ratio and volume fraction index on the plate aerothermoelastic behavior is examined.


Latin American Journal of Solids and
The first studies on the flutter behavior of a panel can be traced back to the works of Houbolt (1958), Bolotin (1963) and Dowell (1966Dowell ( , 1975)).Schaeffer and Heard (1965) studied aeroelastic response of a flat panel exposed to a nonlinear temperature distribution with simply support boundary conditions.Xue and Mei (1993a) investigated the nonlinear flutter response of isotropic panels under thermal effects using FEM.To account for structural nonlinearity, von Karman's large detion plate theory was considered.They studied the effects of nonuniform temperature distribution, panel aspect ratio, and boundary conditions on the flutter behaviors of rectangular and triangular panels.Also, Xue and Mei (1993b) considered the fatigue life of isotropic panels in frequency domain using FEM and investigated the effect of dynamic pressure and temperature on the fatigue life of a panel.
In the modern structures, a new material with desired mechanical and strength properties is generated by combining of different material layers.Thus, Isotropic materials have been replaced by composite materials.In the recent decades, Functionally Graded Materials (FGMs) have attracted considerable attention as materials for various advanced purposes.Functionally graded materials are a new type of inhomogeneous materials and advanced composites.A functionally graded material is usually a combination of two materials or phases that have gradual transition of properties from one side of sample to another side.This gradual transition allows the thermo-mechanically interface problems in a composite structure such as sharp local stress concentration, delamination and weak thermal resistance, can be impressively decreased (Miyamoto et al., 1999).Across the consecutive development of FGMs, there have been many research works as follows.Praveen and Reddy (1998) investigated the static and dynamic response of functionally graded plates using a plate finite element that accounted for the effects of the transverse shear strains, rotary inertia and the moderately large rotations in the von Karman sense.They showed that the response of the plates with material properties between those of the ceramic and metal is not intermediate to the responses of the ceramic and metal plates.Javaheri and Eslami (2002) derived stability equations of a rectangular functionally graded plate using with the classical plate theory.Sohn and Kim (2008) analyzed structural stability of functionally graded panels under simultaneous thermal and aerodynamic loads.In this work, thermal post-buckling behaviors and stability boundaries for clamped and simply support FG panels with uniform temperature gradients were studied.Sohn and Kim (2009) studied the aerothermoelastic instability of FGM panels in supersonic flow and showed the effects of the volume fraction distributions, temperature changes, aerodynamic pressures and the boundary conditions on the panel flutter.Navazi and Haddadpour (2011) analyzed the nonlinear flutter behavior and stability boundaries of isotropic and FGM plates in supersonic airflow, and showed that under real flight conditions and using coupled model, the aerodynamic heating is very severe and the type of instability is divergence.Janghorban and Zare (2011) examined the vibration analysis of a functionally graded plate with cutouts and skew boundary.In this work, the role of different parameters such as cutout size, type of loading and different boundary conditions on the vibration of the plate was reported.Taj and Chakrabarti (2013) investigated a finite element formulation based on Reddy's higher order theory to investigate the dynamic response of a FG skew shell.
So far, different analysis techniques have been used for panel flutter studies.However, the methods having less computational complexity (and effort) and greater accuracy are of most interest to researchers.One of them is DQM which first presented by Bellman and Casti (1971).The Latin American Journal of Solids and Structures 13 (2016) 796-818 main idea behind the DQM is that the derivative of a function with respect to a space variable at a given point is approximated as a weighted linear sum of the function values at all discrete points along the domain of that variable.
The DQM has been applied to solve various structural elements such as beams, plates and shells.Bert et al. (1988) applied the DQM to investigate static and dynamic response of structures for the first time, and afterwards it was improved by Bert and Malik (1996).Also, Bert et al. (1989) used DQM for composite plates for the first time and analyzed nonlinear bending of orthotropic rectangular plates.Then, Shu and Richards (1992) presented the GDQM to simplify the computation of the weighting coefficients.Shu and Wang (1999) applied GDQM for vibration analysis of a rectangular plate with combined and non-uniform boundary conditions.Fazelzadeh et al. (2007) investigated the vibration of a rotating thin walled-blade made of functionally graded materials operating under high temperature supersonic gas flow with DQM.Talebitooti et al. (2013) presented the effects of boundary conditions and axial loading on the frequency characteristics of rotating laminated conical shells with orthogonal stiffeners using GDQM.
This paper extends the application of the DQM to aerothermoelastic analysis of a flat plate in supersonic flow.In this regard the governing differential aeroelastic equations of a FG plate as first discretized using GDQM and then aerothermoelastic response of the plate is studied by the fourth order Runge-Kutta method.To demonstrate the accuracy and fidelity of the GDQM, the dynamic stability boundaries of the plate are validated with available results presented by other researchers.Also, the effects of some important parameters such as Mach number, in-plane thermal load, plate aspect ratio and volume fraction index on the plate aerothermoelastic behavior are investigated.

Structural Model
A plate with length a, width b, and thickness h made up of a mixture of ceramic and metal is considered.The airflow is assumed in the x-direction.Volume fraction of the functionally graded material varies continuously through the plate thickness according to a simple power law (Javaheri and Eslami, 2002).Hence: where z is the coordinate in the thickness direction with origin at the plate mid-surface and n is the volume fraction index.Thus, the material properties of the functionally graded plate can be expressed as Based on the Classical Plate Theory (CPT), the displacement field of the plate is: w x y t u x y z t u x y t z x w x y t v x y z t v x y t z y w x y z t w x y t where u 0 and v 0 are the in-plane displacement components, and w 0 is the out-of-plane displacement component measured from the plate's mid-plane.
According to the von Karman nonlinear strain-displacement relations, the nonlinear strains are defined as: where 0  and k are the mid-plane membrane and bending strain vectors, respectively, which can be defined as follows: The thermoelastic constitutive equations of the FG panels are:     Latin American Journal of Solids and Structures 13 (2016) 796-818 where T0, T(z) and ( , ) z T  are reference temperature, temperature distribution in the plate thickness direction and thermal expansion coefficient, respectively.Qij are the elements of material constant matrix and defined by where E(z) is the elastic modulus of a FG panel.
In-plane force resultant and out-of-plate moment resultant are obtained as follows: Here, N T and M T are the thermal in-plane force and moment resultant vectors.Thus while A, B and D are extension, bending-extension coupling, and bending stiffness and are given as follows:

Aerodynamic Model
According to the first order piston theory (Dowell, 1975), the aerodynamic pressure, may be expressed as where q, M and U  represent the dynamic pressure, Mach number, and the free stream velocity, respectively.The linear piston theory is valid for 2 5 M   (Mei et al., 1999).

Temperature Distribution
The temperature distribution on the surface of the plate is assumed to be constant while in the thickness direction it is considered to be variable and may be obtained by solving the onedimensional Fourier equation of the heat conduction, which is Tm and Tc are the temperature of the lower and upper surfaces of the panel, and temperature distribution in the plate thickness direction is obtained by means of polynomial series (Javaheri and Eslami, 2002) and given as follows: where

Equations of Motion
By means of the extended Hamilton's principle, the nonlinear governing equations of motion can be obtained.In the absence of surface shearing forces, body moments and inertial forces in the x and y directions, the aeroelastic equations of a FG plate are (Bert et al., 1989) where By differentiating from Eq. ( 23) and Eq. ( 24) with respect to the variable of x and y, respectively, and summing the resultant equations, we will have By multiplying 11 B to Eq. ( 25) and 11 A to Eq. ( 26) and substituting Eq. ( 25) into Eq.( 26) and after some mathematical manipulations, the following equation is obtained.
The above equation can be written as  In this paper, the plate is considered to be simply supported along all edges and therefore the in-plane displacement components , , , are equal to zero and based on this assumption the in-plane force resultants can be modified as follows (Dowell 1975, Miller et al., 2011): By substituting Eqs. ( 16) and (31) into Eq.( 28), the aerothermoelastic equation is obtained.
In a supersonic regime, the coefficient of the aerodynamic damping in Eq. ( 31) can be approximately expressed as (Liaw and Yang, 1993) In order to express Eq. ( 31) in a non-dimensional form, the following set of dimensionless parameters is defined as follows: Latin American Journal of Solids and Structures 13 (2016) 796-818 By introducing the above parameters into Eq.( 31), the non-dimensional form of aeroelastic equations can be obtained.

DISCRETIZED FORM OF THE GOVERNING EQUATION
In this section, the governing aeroelastic equation is discretized by using the GDQM.This method implies that the rth-order derivative of a function W, at a point s = si , with N discrete points can be estimated by The coefficient

( ) r s ij
A represents the weighting coefficeints.The method for constructing these co- efficients can be found in Chang Shu (2000).
The DQ method may also be used for linear combinations of partial derivatives and integrals (Chang Shu, 2000) as follows: Latin American Journal of Solids and Structures 13 (2016) 796-818 where , l k c c are the weighting coefficients of the one-dimensional integral in the x and y directions respectively, and given by (Chang Shu, 2000) where In the above equations ( ) k r x and ( ) l r y are the Lagrange interpolated polynomials.By using the DQM, the discretized form of the governing equations (Eq.34) can be written as: Latin American Journal of Solids and Structures 13 (2016) 796-818 , , , , , 2 1, 2,..., Also, the assumed boundary conditions may be defined as follows: (0, , ) (0, , ) 0, 0 ( , , ) ( , , ) 0, 0 ( ,0, ) ( ,0, ) 0, 0 ( , , ) ( , , ) 0, 0 w y t w y t x w a y t w a y t x w x t w x t y The discretized form of these boundary conditions will be: Thus, the boundary conditions can be written as: Latin American Journal of Solids and Structures 13 (2016) 796-818 By incorporating the boundary conditions into Eq.( 41) and doing some manipulations, the final aerothermoelastic equation of a FG plate can be drawn as: , , , , , 2 3, 4,..., 2 where C1 to C16 are defined in Appendix A

RESULTS AND DISCUSSION
In this section, the critical dynamic pressures for FG plates are calculated to investigate the validity of the GDQM for determining flutter boundaries.In this regard, a simply-supported functionally graded flat plate made of a combination of metal (SUS304) and ceramic (Si3N3) is considered as a test case.The material properties are listed in Table 1 (Navazi and Haddadpour, 2011).
In order to obtain the flutter instability, the governing aerothermoelastic equation (Eq.45) is utilized by considering 11 × 11 sampleing points in the computational domain.However, after some investigations, it was found that the present nonlinear analysis is extremely sensitive to the sampling point distribution.Thus, Chebyshev-Gauss-Lobatto distribution is not a right choice for this type of problems, but the suggested distribution by Tomasiello (1998) can be used as an alternative choice.Then, the resulting ordinary differential equations are solved via the 4th order Runge-Kutta method.Table 2 shows the obtained critical non-dimensional dynamic pressures (λcr) of a FG square plate with a/b=1(aspect ratio), a/h=100 (length/ thickness ratio) along with those reported by Sohn and Kim (2008) under different volume fractions and temperatures.It is clear that the obtained results are in good agreement with Sohn and Kim (2008).Next, to validate the post flutter characteristics of the present computational model, the limit cycle amplitudes of an isotropic square panel for various dynamic pressures, at a specified location ( 0 .7 5   and 0 .5

 
) on the panel are computed and depicted in Fig. 1.As is evident, the DQM results are in satisfactory agreement with those presented by Dowell (1966) and Xue and Mei (1993) which used time integration and finite element method, respectively.Finally, the aerothermoelastic stability boundaries of the panel according to in-plane thermal load are depicted in Fig. 2. As shown in this figure, the developed tool can capture different dynamic behaviors of the plate, including stable, limit cycle oscillation (LCO), Buckled and chaotic, as well as Xue and Mei (1993).the other hand, any increase in volume fraction index leads to reduction in the critical dynamic pressure.However, this reduction is impressive until the value of volume fraction becomes about unity.
The effects of aspect ratio (a/b) on the critical dynamic pressure are presented in Fig. 4. It can be seen that the critical dynamic pressure increases with the plate aspect ratio.Moreover, the variation of critical dynamic pressure with respect to µ/M is more dominant for plates with higher aspect ratio.The limit cycle amplitude of the FG plate at the aforementioned point versus the critical dynamic pressure for various volume fraction indexes and different values of Tc is shown in Figure 5.It should be noted that the value of Tc is considerd to vary from 300 K (the gary lines) to 320 K (the black lines) and Tm is taken to be 300 K.It is found that by increasing the surface temperature of the FG plate, the critical dynamic pressure is decreased.The aerothermoelastic stability margins versus in-plane thermal load for various volume fraction indexes are shown in Fig. 6.It can be seen that as the volume fraction index decreases the critical thermal buckling load, the critical dynamic pressure and the in-plane thermal load, where the chaotic motion begins, all increase.As a result, the bifurcation point shifts to the top right.As shown in Fig. 6, the plate remains stable in zone A. The typical time history of this zone is shown in Fig. 7.In zone B, the panel is dynamically stable, but statically unstable and buckled.The dimensionless deflection time histories in zone B for two specific conditions are shown in Figs. 8 and 9.They reveal that the steady response of the panel increases with increasing of n.They also show that the amplitude of the panel motion decreases with increasing of dynamic pressure.As mentioned before, the limit cycle oscillation occurs in zone D. The deflection time history and the related phase diagram in this zone for two specified conditions are shown in Figs. 12 to 15.However, for moderate Rx as dynamic pressure increases, the periodic limit cycle oscillation changes to the simple harmonic limit cycle oscillation.

CONCLUSIONS
In this study, the GDQM has been used to provide the ordinary differential form of the aerothermoelastic equation of a FG flat plate.To investigate the aerothermoelastic behavior of the plate, the 4 th order Runge-Kutta numerical method has been utilized.The evaluation of the obtained results in comparison with those available in literature shows the accuracy and ability of the GDQM to study the aerothermoelastic behavior of a FG panel in supersonic regime.However, it should be noted that the application of this method is much simpler than other well-known computational methods such as (Galerkin, and FEM).Also, it was found that the nonlinear panel flutter analysis with GDQM is extremely sensitive to the grid point distribution.Therefore, the well-known Chebyshev-Gauss-Lobatto distribution is not suitable for this type of problems, but Tomasiello's distribution (Tomasiello, 1998) can be suggested as an efficient and suitable choice.

Figure 1 :
Figure 1: LCO amplitudes of the isotropic square panel

Figure 2 :
Figure 2: Stability boundaries of the isotropic square panel

Fig. 3
Fig.3shows the critical dynamic pressure of the aforementioned FG panel versus volume fraction index (n).It can be seen that as µ/M decreases the critical dynamic pressure decreases too.On

Figure 3 :Figure 4 :
Figure 3: Effect of µ/M on the critical dynamic pressure of the FG panel (a/h=100)

Figure 5 :
Figure 5: Deflection LCO dimensionless amplitude of the FG panel.

Figure 6 :
Figure 6: Stability boundaries of the FG panel at various n

Table 1 :
Material properties at reference temperature (T= 300 K)

Table 2 :
Critical non-dimensional dynamic pressure