Analytical Solution for Thermal Flutter of Laminates in Supersonic Speeds

Lintong Dai, Qiaozhen Sun

Institute of Solid Mechanics, Beihang University (BUAA), Beijing, China.

DOI: 10.4236/jamp.2020.88118                

Abstract

As a basic component of engineering fields such as aeronautics, astronautics and shipbuilding, panel structure has been widely used in engineering and scientific research. It is of great theoretical and practical significance to study the vibration of panels. The panel flutter problem has caused widely concerned by researchers at home and abroad during to the emergence of high-speed aircrafts. With regard to the eigenvalue problem of rectangular panels, it is generally believed that it is difficult to obtain a closed form eigen solution in the case of an adjacent boundaries clamped-supported or a free boundary that cannot be decoupled. Aiming at the problem, this paper studies the two-dimensional symmetric orthogonal laminated plate structure in the hypersonic flow in the thermal environment, and combines the first-order piston aerodynamic theory to study a high-precision separation variable method. Through this method, analytical solution to the closed form of the thermal flutter problem of rectangular panels can be obtained under any homogeneous boundary conditions.

Keywords

Two-Dimensional Panel, Thermal Flutter, Separation Variable Method, Analytical Solution

1. Introduction

Panel flutter is one kind of typical self-excited vibration in aero-elastics that can cause fatigue damage to the structure. This phenomenon was first observed in 1940’s [1], and was clearly observed in experiments in 1950’s [2], Mei [3] gave a summary of which before 1999. Note that piston theory which was developed to approximate gas pressure by Lighthill [4] in 1953, Forsching [5] summarized three available conditions of using piston theory ω*2M_a2 >> 1, ω*2M_a >> 1, M_a2 >> 1, where ω* is the reduced frequency, M_a is the Mach number. Dowell [6] [7] gave a typical investigation on nonlinear panel flutter of 2D and three dimensional (3D) isotropic and simply-supported panels. Impact on 3D panel flutter properties of different boundary conditions (SSSS, SCSC and SFSF) were investigated using the assumed mode method and finite element method in 2014 [8]. Yufeng Xing gave the overall assessment of closed-form solution methods for free vibrations of rectangular thin plate [9]

With the increase of the flight speed of modern aircrafts, the Mach number can reach more than five Mach. Therefore, it is urgent to systematically give the flutter calculation results under all constraint boundaries of the two-dimensional plate model. At the same time, due to the increase of the Mach number, the aerodynamic thermal effect is also not negligible. Based on this, this paper studies a high-precision separation variable method based on two-dimensional symmetric orthogonal laminates, and obtains the exact solution of the two-dimensional panel thermal flutter problem under various homogeneous boundaries (SS, GG, CC, FF, GS, SG, SC, SF, GC, CG, GF and CF). The thermal flutter characteristics of two-dimensional panels are analyzed from the perspective of eigen roots. Finally, the research work on the eigenvalue problem of two-dimensional panel flutter is summarized.

2. Establishment of Basic Equations

Figure 1 is a two-dimensional symmetric orthotropic laminate model with a chord length of a,a thickness of h, and a plate density of ρm. The upper surface of the panel has airflow, and the airflow density, velocity and Mach number are respectively recorded as ρ_a, V, M_a. The temperature is evenly distributed when the panel reaches a steady state after being heated. The laminate is composed of five orthotropic and aeolotropic layers of equal thickness and uniform as shown in Figure 1(b). The two-dimensional panel model is of an infinite length in

Figure 1. Two-dimensional panel subjected to aerodynamic loading over one surface. (a) Plate configuration; (b) Ply stacking sequence.

spanwise, elements per spanwise unit length can be used in the process of analysis The xy plane of cartesian coordinate system is established in the middle of the panel, and the origin point O is built in the corner point of the unit, as shown in Figure 1.

Using the classical laminates theory, which satisfies the Kirchhoff hypothesis, for symmetric cross-ply composite laminates, the motion equation of the panel is

$D_{11}\frac{{\partial }^{4}w}{\partial x^4}-\left(N_x-N_x{}^{\left(T\right)}\right)\frac{{\partial }^{\text{2}}w}{\partial x^2}+{\rho }_{m}h\frac{{\partial }^{2}w}{\partial t^2}+q_a=0$ (1)

$\begin{array}{l}{D}_{11}=\frac{1}{3}{\displaystyle \underset{k=1}{\overset{N}{\sum }}{\left({\bar{Q}}_{11}\right)}_k}\left(z_k^3-z_{k-1}^3\right)\\ z_k=-\left[\frac{h}{N}\left(\frac{N-1}{2}-k\right)+\frac{h}{2N}\right]\end{array}$ (2)

where D₁₁ is the plate bending rigidity, N is the number of stacking layers, and ${({\bar{Q}}_{ij})}_k$ is the transformed reduced stiffness coefficient for the kth layer, and

${\left({\bar{Q}}_{11}\right)}_k=Q_{11}{\cos }^4\theta +2\left(Q_{12}+2Q_{66}\right){\sin }^2\theta {\cos }^2\theta +Q_{22}{\sin }^4\theta$ (3)

For a orthotropic plate,

${\left({\bar{Q}}_{11}\right)}_{\theta =0^{\circ }}=\frac{E_1^2}{E_1-{\nu }_{21}^2{E}_2},{\left({\bar{Q}}_{11}\right)}_{\theta ={90}^{\circ }}=\frac{E_1{E}_2}{E_1-{\nu }_{21}^2{E}_2}$ (4)

when N = 5, the bending rigidity is

$D_{11}=\frac{1}{3}\left(\frac{E_1^2}{E_1-{\nu }_{21}^2{E}_2}\frac{49h^3}{250}+\frac{E_1{E}_2}{E_1-{\nu }_{21}^2{E}_2}\frac{13h^3}{250}+\frac{E_1^2}{E_1-{\nu }_{21}^2{E}_2}\frac{h^3}{500}\right)$ (5)

N_x is the mid-plane compressive force per unit length, $N_x{}^{\left(T\right)}$ is the compressive force caused by the temperature changes $\Delta T$

$N_x{}^{\left(T\right)}={\displaystyle {\int }_{-h/2}^{h/2}{E}_E{\alpha }_E\Delta Tdz=E_{E}h{\alpha }_E\Delta T}$ (6)

$E_E$ is the equivalent elastic modulus and ${\alpha }_E$ is the equivalent thermal expansion coefficient of the laminates.

The supersonic unsteady aerodynamic force is calculated by the piston theory. The aerodynamic load can be expressed by the classical first-order piston theory, and when the M_a is large, it can be approximated as

$q_a=p-p_{\infty }=-\frac{2q}{M_a}\left(\frac{\partial w}{\partial x}+\frac{1}{V}\frac{\partial w}{\partial t}\right)$ (7)

Then Equation (1) can be written as

$D_{11}\frac{{\partial }^{4}w}{\partial x^4}-\left(N_x-N_x{}^{\left(T\right)}\right)\frac{{\partial }^{\text{2}}w}{\partial x^2}+{\rho }_{m}h\frac{{\partial }^{2}w}{\partial t^2}+\frac{2q}{M_a}\left(\frac{\partial w}{\partial x}+\frac{1}{V}\frac{\partial w}{\partial t}\right)=0$ (8)

In this paper, it is the exact eigensolution of Equation (7) that we need to obtain. To solve the control partial differential equation, we need to meet the corresponding boundary conditions. The various classical boundary conditions are shown in Table 1.

Table 1. Classical boundary conditions.

3. Exact Frequency and Mode Functions

In this section, the exact eigensolutions of 2D panel flutter are derived for the cases of SS, GG, CC, FF, SG, SC, SF, GC, GF and CF, in which SS and CC are the most used in previous analysis about 2D panel flutter, CF is rarely used, and the remaining are discussed for the first time.

Let the deflection w be in the form of a separate variable as follows

$w=\varphi (x)\tau (t)=\varphi (x){\text{e}}^{\Omega }{}^t,\Omega =\beta +\text{i}\omega$ (9)

The real part β of Ω represents amplitude variation, while imaginary part ω represents the frequency of principle vibration. If $\beta \ge 0$ the panel flutters. Substituting Equation (8) into Equation (7) yields homogenous characteristic equation.

${\rho }_{m}h\varphi {\Omega }^2+\frac{2q}{M_{a}V}\varphi \Omega +\left[D_{11}\frac{{\text{d}}^4\varphi }{\text{d}{x}^4}-\left(N_x-N_x{}^{\left(T\right)}\right)\frac{{\text{d}}^2\varphi }{\text{d}{x}^2}+\frac{2q}{M_a}\frac{\text{d}\varphi }{\text{d}x}\right]=0$ (10)

from which we can solve flutter mode function and flutter frequency for different boundary conditions. The flutter mode function or eigenfunction has the form as

$\varphi \left(\xi \right)=A{\text{e}}^{\lambda \xi }$ (11)

where λ is the eigenvalue with respect to $\eta =x/a$ . Substituting of Equation (10) into Equation (9) yields algebraic eigenvalue equation

${\lambda }^4-\frac{\left(N_x-N_x{}^{\left(T\right)}\right)a^2}{D_{11}}{\lambda }^2+\frac{2q{a}^3}{D_{11}{M}_a}\lambda +\left(\frac{{\rho }_{m}h{\Omega }^2{a}^4}{D_{11}}+\frac{2q{a}^4}{D_{11}{M}_{a}V}\Omega \right)=0$ (12)

That is

${\lambda }^4+R{\lambda }^2+p_l\lambda +k=0$ (13)

Then Equation (12) is the two-dimensional panel flutter Eigen algebraic equations, where R, p_l and k are all nondimensional parameters, and p_l is aerodynamic coefficient, k is frequency prameter.

$\{\begin{array}{l}R=\frac{-\left(N_x-N_x{}^{\left(T\right)}\right)}{a^2}\frac{a^4}{D_{11}}\\ p_l=\frac{2q}{M_{a}a}\frac{a^4}{D_{11}}\\ k=\left({\rho }_{m}h{\Omega }^2+\frac{2q}{M_{a}V}\Omega \right)\frac{a^4}{D_{11}}\end{array}$ (14)

The parameter k in Equation (13) can also be written as

$k={\pi }^4{\left(\frac{\Omega }{{\omega }_{rs}}\right)}^2+{\pi }^4{g}_a\left(\frac{\Omega }{{\omega }_{rs}}\right)$ (15)

where ${\omega }_{\text{rs}}={\left(\pi /a\right)}^2\sqrt{D_{11}/{\rho }_{m}h}$ is the first-order natural frequency of SS panel without aerodynamic force, and aerodynamic damping coefficient $g_a$ is

$g_a=\frac{{\rho }_a{a}_c}{{\rho }_{m}h{\omega }_{rs}}$ (16)

where ρ_a is the mass density of fluid, a_c is the local velocity of sound.

According to Ferrari’s method, the four characteristic roots of Equation (12) can be

$\{\begin{array}{l}{\lambda }_{1,2}=\vartheta \pm \text{i}{\alpha }_1\\ {\lambda }_{3,4}=-\vartheta \pm {\beta }_1\end{array}$ (17)

And the general solution of the eigenfunction can be expressed as

$\varphi (\xi )=A_1{\text{e}}^{{\lambda }_1\xi }+A_2{\text{e}}^{{\lambda }_2\xi }+A_3{\text{e}}^{{\lambda }_3\xi }+A_4{\text{e}}^{{\lambda }_4\xi }$ (18)

Substitute Equation (8) and Equation (18) into the boundary conditions shown in Table 1 to determine frequency equations and the coefficients of $\varphi (\xi )$ , and the method to solve eigensolutions $(\Omega ,\varphi )$ is the same for different combination of boundary conditions, thus the case SS is taken as an example to show the solution procedure. The boundary conditions in terms of $\varphi (\xi )$ are

$\varphi (0)=\varphi (1)=0,{\varphi }^{″}(0)={\varphi }^{″}(1)=0$ (19)

Substitution Equation (18) into Equation (19) results in four homogeneous algebraic equations for undetermined coefficient A₁, A₂, A₃ and A₄ as

$\left[\begin{array}{cccc}1& 1& 1& 1\\ {\lambda }_1^2& {\lambda }_2^2& {\lambda }_3^2& {\lambda }_4^2\\ {\text{e}}^{{\lambda }_1}& {\text{e}}^{{\lambda }_2}& {\text{e}}^{{\lambda }_3}& {\text{e}}^{{\lambda }_4}\\ {\lambda }_1^2{\text{e}}^{{\lambda }_1}& {\lambda }_2^2{\text{e}}^{{\lambda }_2}& {\lambda }_3^2{\text{e}}^{{\lambda }_3}& {\lambda }_4^2{\text{e}}^{{\lambda }_4}\end{array}\right]\left[\begin{array}{c}{A}_1\\ A_2\\ A_3\\ A_4\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right]$ (20)

After substituting the eigenvalue expression Equation (17) into the above equation, the frequency equation and mode function coefficients of two-dimensional simply supported plates can be obtained:

$\begin{array}{l}16{\vartheta }^2{\alpha }_1{\beta }_1\cosh 2\vartheta +\text{i}{\left({\beta }_1+\text{i}{\alpha }_1\right)}^2\left[4{\vartheta }^2-{\left({\beta }_1-\text{i}{\alpha }_1\right)}^2\right]\cosh \left({\beta }_1+\text{i}{\alpha }_1\right)\\ -\text{i}{\left({\beta }_1-\text{i}{\alpha }_1\right)}^2\left[4{\vartheta }^2-{\left({\beta }_1+\text{i}{\alpha }_1\right)}^2\right]\cosh \left({\beta }_1-\text{i}{\alpha }_1\right)=0\end{array}$ (21)

$\begin{array}{l}{A}_1=-A_2-A_3-A_4\\ A_2=-\frac{\left({\lambda }_4^2-{\lambda }_1^2\right)\left({\text{e}}^{{\lambda }_3}-{\text{e}}^{{\lambda }_4}\right)}{\left({\lambda }_2^2-{\lambda }_1^2\right)\left({\text{e}}^{{\lambda }_3}-{\text{e}}^{{\lambda }_2}\right)}{A}_4\\ A_3=-\frac{\left({\lambda }_4^2-{\lambda }_1^2\right)\left({\text{e}}^{{\lambda }_4}-{\text{e}}^{{\lambda }_2}\right)}{\left({\lambda }_3^2-{\lambda }_1^2\right)\left({\text{e}}^{{\lambda }_3}-{\text{e}}^{{\lambda }_2}\right)}{A}_4\end{array}$ (22)

Table 2 eigenvalue properties of different flutter types of all boundaries.

Table 2. Eigenvalue properties of different flutter types.

4. Numerical Analysis

The equivalent elastic modulus of laminated plates is

$E_E=E_1(1-c_f)+E_2{c}_f$ (23)

The equivalent thermal expansion coefficient of laminated plates is

${\alpha }_E=\frac{{\alpha }_1{E}_1(1-c_f)+{\alpha }_2{E}_2{c}_f}{E_1(1-c_f)+E_2{c}_f}$ (24)

And Table 3 shows the Parameters of the panel and supersonic flow.

Table 3. Parameters of the panel and supersonic flow.

4.1. Flutter Frequency and Flutter Type

The relationship among β, ω and M_a can be obtained from solving Equation (21) for case SS as shown in Figure 2 and Figure 3. Figure 2 shows that β is a negative constant when M_a < M_acr = 6.9521, implying that the vibration of panel before flutter is damping vibration.

It can be concluded from Figure 3 that ω_1st and ω_2nd get closer as M_a increases, and comes to an equal at M_acr. Then two frequencies keep equivalence, β begins to rise till β = 0 when M_a = M_af = 6.9896. Due to the existence of aerodynamic damping, M_af > M_acr.

Above qualitative conclusions are for case SS, but they are also correct for the cases of CC, FF, GG and SC etc. and all the flutter types of these panels are coupled-mode as shown in Table 2.

For panel GC, Figure 4 and Figure 5 show its flutter characteristic. The first two order frequencies never coincide as the Mach number increases. When M_a = M_acr, ω_1st = 0, then β = 0, when M_a = M_af, panel flutters, after this moment, panel flutter diverges, and this type of flutter is called zero-frequency flutter or static divergence. GC and CG have different flutter type due to the asymmetry of system

Figure 2. Relation between β and M_a for SS.

Figure 3. Relation between ω and M_a for SS.

Figure 4. Relation between β and M_a for GC.

Figure 5. First two order frequencies for GC.

stiffness, the former is couple-mode flutter while the latter is zero-frequency flutter. Besides, flutter can hardly happen for case CG while it is easy for case GC to have a zero-frequency flutter. The case GS and SG have the different flutter characteristics, M_acr (GS) = 2.0025 while M_acr (SG) = 51.1085.

Figure 6 shows the relationship between temperature and the flutter boundary under several typical boundary conditions of frequency coincidence flutter. When the temperature is lower than the critical buckling temperature (ΔT_cr (CC) = 57.87˚C, ΔT_cr (SS & GG) = 14.47˚C, ΔT_cr (CF) = 3.6˚C), the flutter boundary can be obtained from the linearized model analysis proposed in this

article. The stiffness of the system is reduced due to the temperature rises, so the flutter aerodynamic coefficient decreases.

Temperature can affect the critical Mach number of the panel and flutter boundary. When it reaches the critical thermal buckling temperature, buckling

Figure 6. Flutter boundary versus temperature.

Table 4. Contrast of M_acr (SS) in different ΔT.

Table 5. Comparison between Galerkin’s method and exact solution.

will occur. And only the effective stiffness of the system is changed so as to affect the flutter boundary while no new flutter phenomenon occurs.

4.2. Results Contrast

We know that in the field of flutter theory analysis, the Galerkin’s method is widely used, so we use it to verify the method proposed in this paper furtherly.

Table 4 shows the contrast of M_acr in △T = 0 and ΔT_cr in SS boundry.

Table 5 lists the results of Galerkin’s method and exact solution in this paper. Generally, the third-order panel flutter calculation is more reasonable when using Galerkin’s method. When the third-order mode is selected, the calculation results of the Galerkin’s method coincide the results of this paper well. The relative error of the flutter frequency is 0.77%, and β, ω_1st and ω_2nd converge faster than M_af and ω_f, the relative error of M_af is 2.69%.

5. Summary

In this paper, all possible exact eigen solutions of two-dimensional panel flutter under any homogeneous boundaries are obtained by a unified method. When the critical point of the flutter is reached, the eigen root will become a complex number, so that the system vibrate diverges exponentially with time. And after the critical temperature, the linearized model is no longer suitable. This is because when structural buckling happens, the geometric nonlinear effect becomes the main factors that affect the inherent characteristics of structural vibration, linearization technique can no longer describe the buckling and post-buckling behavior of the structure. The buckling and post-buckling behavior of the structure is outside the scope of this paper.

Although the research in this paper is based on two-dimensional panel and linear theory, the calculation results are still of great use for evaluating numerical methods. The solution steps can provide a reference for other similar stability problems such as three-dimensional panel flutter and wing flutter.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Garrick, I.E. and Reed III, W.H. (1981) Historical Development of Aircraft Flutter. Journal of Aircraft, 18, 897-912. https://doi.org/10.2514/3.57579
[2]	Bohon, H.L. and Dixon, S.C. (1964) Some Recent Developments in Flutter of Flat Panels. Journal of Aircraft, 1, 280-288. https://doi.org/10.2514/3.43594
[3]	Mei, C., Abdel-Motagaly, K. and Chen, R. (1999) Review of Nonlinear Panel Flutter at Supersonic and Hypersonic Speeds. Applied Mechanics Reviews, 52, 321-332. https://doi.org/10.1115/1.3098919
[4]	Lightill, M.J. (1953) Oscillating Airfoils at High Mach Number. Journal of the Aeronautical Sciences, 20, 402-406. https://doi.org/10.2514/8.2657
[5]	Försching, H.W. (1974) Grundlagen der Aeroelastik. Springer-Verlag. https://doi.org/10.1007/978-3-642-48285-4
[6]	Dowell, E.H. (1966) Nonlinear Oscillations of a Fluttering Plate. AIAA Journal, 4, 1267-1275. https://doi.org/10.2514/3.3658
[7]	Dowell, E.H. (1967) Nonlinear Oscillations of a Fluttering Plate. II. AIAA Journal, 5, 1856-1862. https://doi.org/10.2514/3.4316
[8]	Song, Z. and Li, F. (2014) Investigations on the Flutter Properties of Supersonic Panels with Different Boundary Conditions. International Journal of Dynamics and Control, 2, 346-353. https://doi.org/10.1007/s40435-013-0038-5
[9]	Xing, Y., Sun, Q., Liu, B., et al. (2018) The Overall Assessment of Closed-Form Solution Methods for Free Vibrations of Rectangular Thin Plates. International Journal of Mechanical Sciences, S0020740318302923. https://doi.org/10.1016/j.ijmecsci.2018.03.013