Liuqing YE,Zhengyin YE,Xiaochen WANG

School of Aeronautics,Northwestern Polytechnical University,Xi’an 710072,China

KEYWORDS

Abstract The critical conditions for aeroelastic stability and the stability boundaries of a flexible two-dimensional heated panel subjected to an impinging oblique shock are considered using theoretical analysis and numerical computations,respectively.The von-Karman large deflection theory of isotropic flat plates is used to account for the geometrical nonlinearity of the heated panel,and local first-order piston theory is employed in the region before and after shock waves to estimate the aerodynamic pressure.The coupled partial differential governing equations,according to the Hamilton principle,are established with thermal effect based on quasi-steady thermal stress theory.The Galerkin discrete method is employed to truncate the partial differential equations into a set of ordinary differential equations,which are then solved by the fourth-order Runge-Kutta numerical integration method.Lyapunov indirect method is applied to evaluate the stability of the heated panel.The results show that a new aeroelastic instability(distinct from regular panel flutter)arises from the complex interaction of the incident and reflected wave system with the panel flexural modes and thermal loads.What’s more,stability of the panel is reduced in the presence of the oblique shock.In other words,the heated panel becomes aeroelastically unstable at relatively small flight aerodynamic pressure.

1.Introduction

Panel flutter is a phenomenon of self-excited oscillation which may occur from the interaction of the inertial force,elastic force and the aerodynamic loads induced by the supersonic flow.This self-excited oscillation may cause fatigue of the panel or supporting structure,functional failure of equipment attached to the structure,or excessive noise levels in space vehicle compartments near the fluttering panel.1This phenomenon was first observed by Jordan2who suggested that a lot of early V-2 rocket failures might have resulted from panel flutter.However,formal studies have been carried out based on different structural and aerodynamic theories since the 1950 s.Excellent reviews have been provided by Dowell3and Mei et al.4The early studies on the panel flutter are based on linear theory,and the critical flutter velocity is provided.

Nomenclature

When the velocity of flow is above the value,the motion amplitude of the panel increases exponentially with time.Hedgepeth5and Dugundji6studied the problem of panel flutter based on the linear theory.However,the motion of the panel is generally restrained to a bounded limit cycle because of structural nonlinearities.The linear theory can only determine the critical dynamic pressure,frequency of vibration and mode shape during the instability,and give no information about the panel deflection and stress.7Hence,the nonlinear structure theory should be employed in the analysis of the panel flutter8and von Karman plate theory9has been widely adopted to introduce geometrical nonlinearity.

Potential flow theories aimed at providing predictions of the unsteady aerodynamics have been given by Dowell.10,11Vedeneev12has considered a two-dimensional panel flutter problem using potential gas flow theory.However,the exact potential flow theory is not widely used for the unsteady aerodynamics modeling because of its greater complexity compared to piston theory.Many researchers have employed the piston theory in the analysis of panel flutter since the theory was developed by Lighthill13and Ashley et al.14The ‘‘piston theory”is a method for calculating the aerodynamic loads on aircraft in which the load pressure generated by the body’s motion is related to the local normal component of fluid velocity in the same way that these quantities are related at the face of a piston moving in a one-dimensional channel.In 1966,Dowell adopted the von-Karman plate theory and the first order piston theory to study panel flutter of two-and threedimensional plates and the case has become a classical one in panel flutter studies.Lee et al.15used the first-order piston theory to model aerodynamic loads,and the Newton-Raphson iteration method and complex eigenvalue solver with the Linearized Updated Mode/Nonlinear Time Function(LUM/NTF)approximation method approximation method to obtain the post buckled deflection and flutter information,respectively.Koo and Hwang16applied the finite element method to study the effects of distributed structural damping on the flutter boundaries of composite plates with the linear piston theory used to compute the unsteady aerodynamic load in a supersonic flow.Zhao and Cao17employed the third-order piston theory to estimate the nonlinear aerodynamic pressure induced by the supersonic air flow in numerical analysis of the flutter of a stiffened laminate composite panel.Some other principal publications based on piston theory are mentioned by Bolotin18,Dowell19and Algazin and Kiiko.20Euler equations21,22and Navier-Stokes equations23–26are also used in analysis of panel flutter for calculating the unsteady aerodynamics.Recent experimental27and numerical28investigations for single-degree-of-freedom flutter in the transonic flow regime have been performed by Vedeneev et al.

Aerodynamic heating should be considered in the flutter design when the aircraft flies at high Mach number.Temperature increase due to aerodynamic heating or restrained thermal expansion may result in thermal buckling.In general,two simplifying assumptions are used because of the complexity of aerothermoelasticity coupling equations29,30:(A)structural deformation has no effect on the temperature field;(B)response time of the panel flutter is much smaller than that of temperature change.In 1958,Houbolt31first studied the flutter boundaries and buckling instability characteristics of a two-dimensional plate based on the simplifying assumption of uniform temperature field.Yang and Han32employed thefinite element method to study the thermal buckling flutter problem of a two-dimensional plate using the same assumption of a uniform temperature field.Xue and Mei33,34performed finite element analysis of nonlinear flutter response of isotropic two-and three-dimensional plates with arbitrary shapes.

Most of the existing researches focus on classical panel flutter problem when only one surface of the panel is exposed to supersonic flow and no shock is considered.However,in some engineering applications,such as in the scramjet,both surfaces of the thin panel structure of intakes or nozzles are exposed to air flow and also subjected to an impinging oblique shock or even shock train.Zhou et al.35studied the heated panel flutter problem when both surfaces of the panel are exposed to air flow with different aerodynamic pressures but no shock is considered.The critical conditions for aeroelastic stability and the stability boundaries are obtained using theoretical analysis and numerical computations,respectively.The results show that the panel is more prone to becoming unstable when two surfaces are subject to aerodynamic loading.Visbal36considered the dynamics of a flexible two-dimensional panel subjected to an impinging oblique shock.The flow field above the dynamically deformed panel is computed by means of the full compressible Navier-Stokes equations.The computation results showed that the critical dynamic pressure diminishes with increasing shock strength and can be much lower than that corresponding to standard panel flutter in the absence of a shock.In 2014,Visbal37further investigated the complex self-sustained oscillations arising from the interaction of an oblique shock with a flexible panel in both the inviscid and viscous regimes.The flow fields were obtained by solving either the Euler or the full compressible Navier-Stokes equations employing an extensively validated implicit solver.However,Visbal considered the flutter problem with only one surface subjected to the supersonic flow and the numerical analysis requires a large amount of computational time.

In this paper,we conduct theoretical analysis of aeroelastic stability of two-dimensional nonlinear panel subjected to an impinging oblique shock with only one surface being exposed to supersonic flow and also with both surfaces being exposed to supersonic flow,as depicted in Fig.1 where xidenotes the shock impingement location on the panel.In this paper,only xi/l=0.5 is considered.The von-Karman large deflection theory of isotropic flat plates is used to account for the geometrical nonlinear effects of the heated panel,and the local first order piston theory is employed to estimate the aerodynamic pressure induced by the supersonic air flow.The dynamic pressures before and after the shock wave on one surface are specified as different values.Thus,the effect of the oblique shock is introduced into the governing equations.The temperature change is assumed to be steady state.The Galerkin discrete method is employed to truncate the partial differential equations into a set of ordinary differential equations,which are then solved by the fourth-order Runge-Kutta numerical integration method.Lyapunov indirect method is applied to evaluate the stability of the heated panel.The critical conditions for aeroelastic stability and the stability boundaries are studied using theoretical analysis and numerical computations,and the effects of various parameters,such as dynamic pressure below the panel,the incident shock strength and Mach number,on the stability boundaries are studied.

Fig.1 Schematic of flow configuration for oblique shock impinging on a flexible panel.

2.Governing equations

An isotropic heated plate with simply supported conditions is considered.To simplify the analysis,a two-dimensional model of this panel is utilized.

According to the Hamilton principle,the coupled partial differential governing equation of motion for the panel is established.

where D represents the flexural stiffness of the panel.

2.1.Axial loading

Because of the special configuration of panel,the bending deflection of panel is strongly confined by its geometrical boundaries,and panel flutter usually exhibits nonlinear vibration with limited amplitude which has the same order of panel thickness.Generally,panel flutter will cause fatigue damage rather than rapid destructive failure of the structure.Therefore,geometrical nonlinearity due to large deflection should be considered in panel flutter study.

Von Karman’s theory,with large deformation in the strain equations’results,leads to this force that is a main source on nonlinearity in the governing equation.

2.2.Thermal effect

For aircraft in supersonic and hypersonic air flow,thermal stress induced by aerodynamic heating will decrease the bending stiffness of panel,thus leading to thermal buckling,thermal flutter and more complex dynamic behavior.

To simplify the analysis,the temperature in the panel is assumed to be uniform after it is heated,T0is the initial temperature and α is the coefficient of thermal expansion.Thus,the expression of the thermal stress induced by temperature rise is defined as

Here the temperature rise is defined as ΔT=T-T0.

2.3.Aerodynamic loading

In this paper,the flexible two-dimensional panel is subjected to not only aerodynamic loading on both surfaces,but also an impinging oblique shock on the upper surface.Because of the presence of the incident shock and the reflected shock,the flow field on the upper surface is dominated by the shock wave.In shock-dominated flows,accurate and efficient prediction of aerodynamic pressure is particularly challenging,where the flow is discontinuous and highly nonlinear and leads to severe loading.38

2.3.1.Local piston theory

Brouwer et al.38has proved the feasibility of using local piston theory to predict the aerodynamic loads in shock-dominated flows.In this paper,the Local Piston Theory(LPT)is employed,where relevant freestream flow parameters are replaced with spatially local quantities.The corresponding original pressure relation is given as39,40

where γ represents the air specific heat ratio,vnis the resulting transverse flow velocity adjacent to the panel surface,and plocand alocare the local static pressure and speed of sound on the panel surface respectively.These are different ahead of and behind the shock and are assumed to be available from the simplified formula or Computational Fluid Dynamics(CFD)analysis.

2.3.2.Aerodynamic loading in shock-dominated flows

In this paper,the first-order piston theory is utilized.We use the following formula for the aerodynamic loading on the upper surface ahead of the shock:

The aerodynamic loading on the upper surface behind the shock is as follows:

The aerodynamic loading on the bottom surface with no shock is as follows:

In order to make the complicated problem as simple as possible,here we assume that.In this paper,we focus on the effect of unsteady dynamic pressure on stability boundaries of the heated flexible panel and study it in detail.

The pressure differences ahead of and behind the shock between the upper surface and the bottom surface are as follows:

where βu，l, βu，r, βdare the Prandtl-Glauert factors.

2.4.Non-dimensional equations

The governing Eq.(1)can be put in the following nondimensional form:

where the non-dimensional parameters are defined in Appendix A.

We will consider plates simply supported at both edges,and therefore

2.5.Galerkin approach

Galerkin approach can be applied to the continuous system to reduce it to a multi-degree-of-freedom system.

Considering a simply supported boundary condition,the lateral displacement W（ξ，τ）can be defined as

where qi（τ）are generalized coordinates.By substituting this expression into Eq.(11),multiplying by sin（jπξ）,and integrating from 0 to 1,discrete dynamical equations will be obtained.Because of the presence of the incident oblique shock wave and reflected shock,the aerodynamic force can be defined as follows during the integration procedure when the Galerkin Approach is employed.

Thus the discrete dynamical equations in the presence of an impinging(and reflecting)oblique shock are derived as follows:

where

Assume that the panel deflection^qiconsists of two components:the static deflection qiand small dynamic disturbance εi,and then

The linearized differential equation is obtained:

where A is the Jacobi matrix of Eq.(16)in the equilibrium q1=q2=0.

Assume εi= εi0eΩτ,and the characteristic equation of adjoint system is as follows:

Obviously,qi=0(i=1,2)is the equilibrium of the two order aeroelastic system of the heated panel,and at the equilibrium point,the Jacobi matrix is

The characteristic equation of the adjoint system is

3.Stability of flat form of panel

Stability of the flat form of the panel can be studied by using both the linearized partial differential equation and the linearized ordinary differential equations corresponding to the set of equations.Eq.(15)is a high-dimensional nonlinear autonomous system.To analyse the stability of the panel from theoretical analysis,a two-mode(N=2)expansion of the solution set Eq.(14)is considered.

with coefficients

The static equilibrium equation can be achieved by neglecting all the time-varying terms,and the trivial solution(q1=q2=0)corresponds to the equilibrium of initially flat panel.Lyapunov indirect method is employed to evaluate the stability of this equilibrium.

The Routh-Hurwitz criterion is employed to establish conditions that all the roots of Eq.(21)have negative real parts.Put the classical criterion into Eq.(21),and the following inequalities are obtained:

The boundaries from conditions Eq.(23)are defined as divergence instability boundaries.The violation of the last of conditions Eq.(23)means the transfer of one of the roots through the origin of the λ-plane,i.e.,the divergence(buckling)instability of the initial equilibrium.The boundaries from condition Eq.(24)are usually defined as flutter instability boundary.The violation of condition Eq.(24)represents that a couple of complex roots leave the left side half-plane,and it means flutter instability.

In general cases,the air flow density is very small compared with the density of the panel material,and thus,the values of mass parameters ｛RMu，l，RMu，r，RMd｝are very small.From Ref.35,one can conclude that the mass parameters have little effect on the stability boundaries of the panel,and thus neglecting the effect of the mass parameters is a good approximation when the inequalities about the stability boundary are evaluated.

It is easy to see from Eq.(23)that the first condition is satisfied if the dynamic pressureλand the parameters｛RMu，l，RMu，r，RMd｝ are positive.

Here the following notation n is introduced for the ratio of the dynamic pressure before the incident oblique shock wave and the dynamic pressure after the reflected shock wave,and the notation m is the ratio of the dynamic pressure under the panel and the dynamic pressure before the incident oblique shock wave:

Defining the critical buckling temperature elevation of the panel as ΔTcrwhen there is no air flow,we can replace the non-dimensional thermal load parameterby the following temperature ratio:

The second condition of Eq.(23),when solved for the nondimensional dynamic pressure and the temperature elevation ratio,comes to the form as

The last condition of Eq.(23)solved takes the form

The latter inequality can be simplified to the following form if the right side of the inequality is positive

Finally,consider condition Eq.(24)by substituting Eq.(22)into it:

Solved with respect to the same parameter,this gives the following condition for stability of the system:

For n=1,this reduces to a known result.

4.Buckled equilibrium modes

Because of the presence of the incident oblique shock and the reflected shock,the buckled equilibrium modes are also different from the traditional buckled equilibrium modes.In order to analyze the post-critical behavior of heated panel,the equilibrium modes are evaluated from the static aeroelastic equations.Assuming that the two generalized coordinates are constants,one can obtain the static aeroelastic equations of the heated panel:

where

In addition to the zero solution for q1，q2,non-zero solutions occur when the following equality is satisfied

Thus,one can obtain the parameter boundary from the following expression at which buckling modes abruptly transfer into flutter modes.

The stability boundary of a flat heated panel with aerodynamic loading on both surfaces as previously found in Ref.35is as follows:

where λu, λdare respectively non-dimensional dynamic pressure on the upper surface and on the lower surface of the panel.The condition that λu，l(before the oblique shock)equals λu，r(after the reflected shock),i.e.n=1,corresponds to the case of standard panel flutter in which no incident shock is present.By substituting n=1 into Eqs.(27)-(35),the final stability boundaries correspond to those from Eq.(36).

5.Results and discussion

The non-dimensional parameters in this communication are the dynamic pressure λu，lbefore the incident oblique shock wave above the panel,the temperature elevation ratio ΔT/ΔTcr,the ratio n and the ratio m.For the case of traditional heated panel flutter in which only one surface of the panel is subjected to the supersonic air flow,the nondimensional control parameters for the stability boundaries are only two(i.e.λ,ΔT/ΔTcr).For the case of heated panel flutter in which both surfaces of the panel are subjected to the supersonic air flow without oblique shock,the nondimensional control parameters for the stability boundaries are three(i.e.λu，λd，ΔT/ΔTcr).

Different from these previous studies,it is very clear from Eqs.(27),(28),(31),(35)that for the case of heated panel flutter in which an incident shock is present with one surface subjected to air flow,the non-dimensional control parameters for the stability boundaries become three(i.e. λu，l，n，ΔT/ΔTcr).For the case of heated panel flutter in which an oblique shock is present with two surfaces subjected to air flow,the control parameters become four(i.e., λu，l，n，m，ΔT/ΔTcr).In other words,for the case of heated panel flutter in the presence of an impinging shock,the panel will be stable when the nondimensional dynamic pressure λu，l，λd,the temperature elevation ratio ΔT/ΔTcrand the ratio of the dynamic pressure n due to the presence of the incident oblique shock wave are all satisfied on the stability boundary.

5.1.One surface subjected to supersonic air flow

Let us begin with the study of the stability of the heated panel in the presence of an impinging shock.First,one can considerthe situation in which only one surface of the panel is subjected to the supersonic air flow(i.e.λd=0,m=0).In this paper,the dynamic pressure ratio n is obtained by calculating the oblique shock relations.In general,n is larger than 1 in the presence of an impinging shock,and the value of n increases with increasing shock strength,as shown in Table 1.All results for this part of discussion are for dynamic pressure ratios n= λu，r/λu，l=1，2，3，4，5，6，7，8.The value of n= λu，r/λu，l=1 corresponds to the case of standard panel flutter in which no incident shock is present.The stability boundaries have different characteristics when different dynamic pressure ratios are employed.Thus,to simplify the analysis,the dynamic pressure ratios discussed in this part are divided into two parts.The dynamic pressure ratios of n= λu，r/λu，l=1，2，3，4 are in the first part,and n= λu，r/λu，l=5，6，7，8 are in the other part.

Table 1 Dynamic pressure ratio λu，r/λu，l and pressure ratio p u，r/p u，l with different shock angles σ.

The lines AB,DJ and curve CDB represent the dynamic instability boundaries of the traditional flutter(i.e.no shock),shown in Fig.2(a).In the presence of an impinging oblique shock,the dynamic instability boundaries of the case n=4 have also been shown in Fig.2(b).Considering that Eq.(31)is dependent on dynamic pressure ratios,the line segment A1B1denotes the dynamic instability boundary.The static dynamic instability boundaries determined by Eq.(35),dependent on dynamic pressure ratios,are denoted by the line segment D1J1.The curve C1D1B1F1G1represents the boundaries of Eq.(28),and the region outside the curve C1D1B1F1G1satisfies the Eq.(28).However,because of the positive dynamic pressure and the limit of the boundary A1B1,the top part of the curve C1D1B1satisfies stability requirements.Therefore,non-dimensional values of λu，l，n，ΔT/ΔTcrfrom the bounded region A1E1D1H1C1enclose stable results and values outside this region would lead to flutter or divergence.

There is a characteristic similar to the standard stability analysis(no shock)that the region is also divided into four parts shown in Fig.2(b)determined by Eqs.(27),(28),(31),(35).The four regions are:(A) flat and stable,the heated panel will be stable in the case of the flat panel no matter what the initial condition is;(B)basic flutter,the dynamic instability of Limit Cycle Oscillation(LCO)or Chaos occur and there is no static equilibrium that the panel can reach finally;(C)buckling(statically unstable),the heated panel will converge to different static equilibrium positions depending on initial conditions and it is also defined as post-bucking;(D)transition region(D1B1E1),the panel will be stable in the case of the flat panel or in the case of post-buckled panel.

However,it can be observed that the introduction of the oblique shock changes the instability boundaries obviously,as shown in Fig.2(c).The presence of the oblique shock reduces the stable region.In other words,the stable heated panel is more prone to becoming unstable( flutter or buckling)in the presence of the oblique shock.In general,with an increase of the non-dimensional dynamic pressure ratio,the stable parameter region is reduced.The results shown in Fig.2(a)and(c)indicate that the dynamic pressure for which flutter onset occurs decreases when the temperature elevation equals zero.For the case n ≠ 1 shown in Fig.2(a)–(c),the buckling temperature elevation is till 1ΔTcrwhen no dynamic pressure is considered.For the case n=1(no shock),the panel becomes buckled only if the temperature elevation is larger than 1Tcr.However,it should be mentioned that the panel will become buckled when the temperature elevation is less than 1Tcrin the presence of the oblique shock.It means that when the temperature is not very high ΔT/ΔTcr＜ 1,with an decrease of dynamic pressure,the behavior of the heated panel will not become more stable in the presence of the oblique shock.There is still a region in which stable modes abruptly transfer into buckling modes.Increasing the dynamic pressure ratio,the flutter boundaries are still tangent to the divergence boundary.It is also very clear that the boundary lines between the panel flutter and buckling are nearly parallel when the temperature elevation is relatively large.Thus,in this case,the increment of the critical dynamic pressure of aeroelastic instability keeps the same value with the increase of temperature elevation.

Fig.2 Effect of non-dimensional dynamic pressure ratio on instability boundaries.

The dynamic instability boundaries for the case n=5 can be seen from Fig.2(d).The curve C2H2D2B2is tangent to λu，l-axis.The dynamic instability boundaries for the case n=8 can be found from Fig.2(e).The divergence instability boundaries become two parts:H3D3B3and G3C3(Fig.2(e)).The results shown in Fig.2(d)–(g)indicate that when the dynamic pressure ratios become larger,some interesting phenomena emerge.The space is divided into five parts but still four regions.Different from the traditional flat and stable region,the flat and stable region in the presence of the oblique shock becomes two parts:A3E3D3H3and G3C3O.It should be mentioned that the panel will become buckled when the temperature elevation is zero for the higher dynamic pressure ratios.It means when the panel is subjected to an impinging oblique shock,the panel will become buckled when the panel is only subjected to aerodynamic force without any temperature elevation.With an increase of the dynamic pressure ratio,the flat and stable region similar to triangle region(such as G3C3O)becomes smaller than before.What’s more,the buckling region expands to the longitudinal axis.It means that with an increase of the dynamic pressure ratio,the range of dynamic pressure in which the panel will become buckled becomes larger when the temperature elevation is zero.

For the higher dynamic pressure ratios considered shown in Fig.2(d)–(g),there are some other characteristics similar to those shown in Fig.2(a)–(c).With an increase of the nondimensional dynamic pressure ratio,the stable parameter region is reduced,the flutter boundaries are still tangent to the divergence boundary,the buckling temperature elevation is still 1ΔTcrwhen no dynamic pressure is considered,and the boundary lines between the panel flutter and buckling are also nearly parallel when the temperature elevation is relatively large.

5.2.Two surfaces subjected to supersonic air flow

Fig.3 Effects of shock and dynamic pressure below panel on instability boundaries(two surfaces).

The effect of shock wave on the instability boundaries is shown in Fig.3(a)when two surfaces of panel are subjected to supersonic air flow.The results for this part are for m=2(i.e.λd/λu，l=2).It can be observed from Fig.3(a)that the flat and stable region decreases in the presence of an oblique shock.The critical flutter non-dimensional pressure is λu，l=91.321 for the case of ΔT/ΔTcr=0 when the panel is in the absence of the shock.When the panel is subjected to an impinging shock, flutter onset occurs for λu，l＞ 64.462,and the critical flutter non-dimensional pressure is decreased by 26.859.It also can be found that the panel is more prone to becoming buckled even if the temperature elevation is low in the presence of an oblique shock when two surfaces of the panel are subjected to air flow.

Fig.3(b)shows the effect of the dynamic pressure below the panel on the instability boundaries when the heated panel is subjected to an oblique shock wave.All results for this part are for n=5(i.e. λu，r/λu，l=5).It can be seen from Fig.3(b)that the stable region decreases and the flutter region increases with the increase of the dynamic pressure below the panel.λf1is the critical non-dimensional dynamic pressure beyond which buckling modes transfer into flutter modes.λf2is the critical flutter non-dimensional dynamic pressure when the temperature elevation is zero.Fig.3(c)(d)show the effect of the dynamic pressure below the panel on λf1and λf2respectively.The results show that λf1and λf2decrease with the increase of the dynamic pressure below the panel.In general, flutter is more prone to occurring for the heated panel subjected to an impinging shock when two surfaces of the panel are exposed to the supersonic air flow.However,it should be mentioned that the panel is less likely to become buckled with the increase of the dynamic pressure below the panel shown in Fig.3(b).This phenomenon occurs because the pressure below the panel can counteract part of the pressure above the panel.It also can be found that,with an increase of the dynamic pressure below the panel,the critical temperature elevation beyond which the panel can be buckled increases.In other words,the panel is more difficult to become buckled(i.e.higher temperature is needed)when the dynamic pressure below the panel is large.

5.3.Effect of incident shock strength

The incident shock strength is given by the shock angle σ.The non-dimensional dynamic pressure ratios corresponding to different shock angles are listed in Table 1.The table shows that incident shock strength increases with the increase of shock angles in a certain range of shock angles.All results for this communication are for a single Mach number Ma1=3.5.

The effects of different incident shock strength on the stability boundaries of the heated panel with its only one surface being subjected to air flow are shown in Fig.4.When shock angles are not very large,the space is still divided into four parts. Fig. 4(a) shows instability boundaries of σ =20°，22°，24°，27°.When shock angles become large and the incident shock strength becomes large,the space is divided into five parts.Fig.4(b)shows instability boundaries of σ =29°，32°，35°.It can be observed that the stable region decreases with the increase of shock angle shown in Fig.4(a)(b).When temperature elevation is relatively small,the critical flutter dynamic pressure decreases with an increasing temperature elevation.The non-dimensional critical flutter dynamic pressure is defined as λf2when ΔT/ΔTcr=0.It can be seen that in the case of σ =20°,λf2is 217.1.In the case of σ =35°,λf2is 119,i.e.,the non-dimensional critical pressure decreases by 98.1.The computational results provided by Visbal36indicated that the presence of the strong oblique shock has an important effect on the reduction of stability of the panel,and the critical dynamic pressure diminishes with increasing shock strength and can be much lower than corresponding value in the absence of a shock.All of the above characteristics are in qualitative agreement with the previous computations by Visbal(2012)for a Mach number which is equal to 2.0.

Nonlinear ordinary differential Eq.(16)are solved by the fourth-order Runge-Kutta numerical integration method.The numerical results are employed to prove accuracy of theoretical analysis results.The non-dimensional parameters are RMu，l=RMu，r=Rd=0.01 and Δt=0.001 is used.The point observed is xi/l=0.75.

Fig.5 shows simulation results for λu，l=200,ΔT/ΔTcr=1.The theoretical analysis results shown in Fig.4(a)indicate that the panel will be stable for the case of no shock but flutter for the case of σ =20°，22°，24°，27°.The time history and phase plane diagrams for the case of no shock are shown in Fig.5(a)(b)respectively.It can be seen from Fig.5(a)that the transverse vibration finally converges to zero for the case of no shock when λu，l=200, ΔT/ΔTcr=1.On the contrary,the transverse vibration converges at a single harmonic limit cycle for the case σ =20°，22°，24°shown in Fig.5(c)–(h).Fig.5(i)(j)show the time history and phase plane diagrams for the case of σ =27°.It is found that the transverse vibration converges to a multi-harmonic limit cycle.It is also found that the amplitude of the panel increases with the increase of incident shock strength as shown in Fig.5(c)(e)(g).

Fig.4 Stability boundaries with different incident shock strength(on only one surface).

Fig.5 Simulation results for λu，l=200,ΔT/ΔT cr=1.

Fig.6 shows simulation results for λu，l=50,ΔT/ΔTcr=1.The theoretical analysis results shown in Fig.4(a)indicate that the panel will not be buckled for the case of no shock if the temperature elevation is less than 1ΔTcr.When the panel is subjected to an impinging oblique shock,the panel will become buckled even if the temperature elevation is less than 1ΔTcr.The computational results(Fig.6)show that the transverse oscillation converges to a constant static value(but not zero)(i.e. the panel becomes buckled) for the case of σ =20°，22°，24°，27°.

The time history for the case of no shock and σ =20°is shown in Fig.7(a),(b)respectively.Fig.4(a)shows that the non-dimensional critical dynamic pressure beyond which buckling modes transfer into flutter modes decreases with increase of incident shock strength.The computational results show that the transverse vibration finally converges to zero for the case of no shock.On the contrary,the dynamic response of the panel becomes chaotic for the case of σ =20°.

Fig.6 Simulation results for λu，l=50, ΔT/ΔT cr=1.

Fig.7 Simulation results for λu，l=105,ΔT/ΔT cr=4.

The results shown in Fig.4(b)show that when the incident shock strength becomes strong enough,the panel can be buckled when there is no temperature elevation.It can be found from Fig.4(b)that the panel will be stable for the case of σ =29°but buckling for the case of σ =32°，35°when ΔT/ΔTcr=0，λu，l=35.Fig.8 shows the simulation results for the case.It can be seen that the transverse vibration finally converges to a constant(not zero)for the case of σ =32°，35°but zero for the case of σ =29°.

Fig.8 Simulation results for λu，l=35, ΔT/ΔT cr=0.

Fig.9 Stability boundaries with different incident shock strength.

In general,all of the above computational results are in agreement with the theoretical analysis discussed previously.

Fig.9 shows the effect of shock strength on the instability boundaries of the panel when its inner and outer surfaces are subjected to the supersonic flow.It can be observed that the trend of the effect of incident shock strength is similar to that of the heated panel when only one surface is exposed to air flow.However,it should be mentioned that the panel is more difficult to buckle when its inner and outer surfaces are both subjected to supersonic flow,although the shock strength is relatively strong.Comparing the results in Fig.9 with those in Fig.4,it is found that the aerothermoelastic instability boundaries of the panel exposing only one surface to air flow are more sensitive to the change of incident shock strength than those of the panel exposing both surfaces to air flow.

For the traditional panel flutter of a panel in uniform supersonic flow,six or more modes are needed to achieve modal convergence.For such complex phenomena as considered in the present paper,more than two modes may be needed to provide solutions.Firstly,to verify the presented numerical method,the variation of the limit cycle amplitudes with the dynamic pressures is calculated for comparison with its counterpart obtained by Dowell,10as shown in Fig.10.It can be clearly shown that there exists a good agreement between the two results,and thus the numerical model can be validated.

Fig.10 Comparison of limit cycle amplitudes with varying dynamic pressures.

Fig.11 Comparison of stability boundaries using theoretical and numerical analysis.

Fig.11 shows the stability boundaries for the shock strength σ=22°by using the theoretical analysis and numerical simulation with two modes.Because the transition region is influenced by the initial values,this region is not shown by the numerical simulations.It is found that the theoretical results obtained by two modes agree well with the numerical results with two modes.

Fig.12 Stability boundaries for shock strength σ =22° using two,four and six modes.

Fig.12 shows the stability boundaries for the shock strength σ =22°using two,four and six modes.Such a comparison made in Fig.12 shows that the results with four modes are very close to the results with six modes.It is also found that the static buckled instability boundaries and the parameter boundaries at which buckling modes transfer into flutter modes obtained by various modes are very close.There are some differences between the flutter instability boundaries obtained by two,four and six modes,and the critical nondimensional flutter dynamic pressure obtained by four or six modes is obviously larger than that obtained by two modes.However,the quantitative data of results by using two modes may not be very accurate but the qualitative conclusions obtained by two modes are till satisfactory.

The theoretical analysis results using Routh-Hurwitz Criteria with two modes shown in Fig.4(b)indicate that,for a large shock strength even for no thermal stress,the panel may buckle rather than flutter.Here two,four and six modes are used to study this characteristic.It shows that,from the simulation results by using four and six modes,the same conclusions can also be obtained,as shown in Fig.13.For the shock strength σ =22°,the range of dynamic pressure is approximately 20.5–59.0 where buckling can occur when there is no thermal stress for two modes.For four modes,the range of dynamic pressure is approximately 25.0–39.5.For six modes,the range of dynamic pressure becomes 24.6–41.2.

Fig.13 Simulation results for shock strength σ =35° when ΔT/ΔT cr=0.

Fig.14 Lower and upper bound on flutter boundaries for σ =22°.

Considering that there are three standard panel flutter results to which one can compare the results of the model that includes the shock.The following three simple cases of a uniform flow everywhere on the panel are corresponding to(A)conditions below the panel,(B)conditions upstream of the shock above the panel,and(C)conditions downstream of the shock above the panel.These three standard models for comparison will bound the flutter boundary for the case with a shock on the panel.Fig.14 shows the comparison results for the shock strength σ =22°.The instability boundaries of the panel when only one surface is exposed to air flow are shown in Fig.14(a).The instability boundaries of the panel when both surfaces are exposed to air flow are shown in Fig.14(b)and m=3 is selected as the non-dimensional flow parameter below the panel.It is found that one can obtain upper bound on the flutter boundary when using the flow parameters ahead of the shock everywhere on the panel in a standard flutter analysis,and obtain lower bound on the flutter boundary when using the flow parameters behind the shock everywhere on the panel in a standard flutter analysis.

To verify the presented theoretical analysis conclusions obtained by local first-order piston theory,high order piston theory is used to provide corresponding numerical simulation results.Using local first-order piston theory,the theoretical analysis results agree well with the numerical results shown in Fig.11 for the shock strength σ =22°.Because the transition region is influenced by the initial values,this region is not shown by the numerical simulations.The following simulation results with high order piston theory are compared with the numerical results obtained by first-order piston theory.

Using local second-order piston theory to estimate the aerodynamic pressure,the discrete dynamic equations in the presence of an impinging(and reflecting)oblique shock with only one surface exposed to air flow(i.e.λd=0)are derived as follows:

Eq.(37)is then solved by the fourth-order Runge-Kutta numerical integration method.The quadratic terms in Eq.(37)are introduced by local second-order piston theory,and the other parts of Eq.(37)are the same as those of Eq.(16)obtained by local first-order piston theory.A new nondimensional parameter h/l occurs in the discrete dynamic equations obtained by local second-order piston theory.It can be found from Eq.(37)that when the non-dimensional parameter h/l is very small,the discrete dynamic equations obtained by local second-order piston theory are very close to those obtained by local first-order piston theory.The simulation results shown in Fig.15(a)indicate that when the parameter h/l is relatively small （h/l=0.003）,the stability boundaries obtained by local second-order piston theory agree well with those obtained by local first-order piston theory.In Ref.41,the parameter h/l of the panel is about the order of magnitude of 10-3.Thus,the results obtained by local first order piston theory are reliable.In fact,the parameter h/l represents the stiffness of the panel without changing the material properties.When the stiffness of the panel is relatively small,there is little difference in stability boundaries obtained by local first-order piston theory and by local second-order piston theory.When the stiffness of the panel becomes large(the parameter is selected as h/l=0.045）,there will be a relatively large difference in stability boundaries obtained by local first-order piston theory and by local second-order piston theory shown in Fig.15(b).However,the quantitative data of results by using local first-order piston theory may not be very accurate but the qualitative conclusions obtained by local first-order piston theory are till satisfactory.

The theoretical analysis results using Routh-Hurwitz Criteria with local first-order piston theory shown in Fig.4(b)indicate that,for a large shock strength shock even for no thermal stress,the panel may buckle rather than only flutter.Here local second-order piston theory is used to study this characteristic.It shows that,from the simulation results by using local second-order piston theory,the same conclusions can also be obtained,as shown in Fig.16.For the shock strength,the range of dynamic pressure σ =35°is approximately 20.5–59.0 where buckling can occur when there is no thermal stress for local first-order piston theory.For local second-order piston theory,the range of dynamic pressure is approximately 20.3–63.0.

Fig.15 Comparison of stability boundaries by using local first order piston theory and local second-order piston theory for shock strength σ =22°.

Fig.16 Simulation results by using local first-order piston theory and local second-order piston theory for shock strength σ =35° when ΔT/ΔT cr=0,h/l=0.045.

5.4.Effect of Mach number

The Mach number of the supersonic air flow is an important parameter for the aerothermoelastic stability of the heated panel.Here a shock angle σ =35°and Ma=2.5,3,3.5 are considered.

Fig.17 Stability boundaries with different Ma(on only one surface).

Fig.17 shows the aeroelastic stability boundaries of the heated panel subjected to an impinging oblique shock exposing only one surface to supersonic air flow with different Mach numbers.It can be seen that the critical flutter dynamic pressure gradually decreases with increasing Mach number when the temperature elevation keeps constant.In the case of Ma=2.5, ΔT/ΔTcr=0,the basic flutter of the panel will occur when the non-dimensional dynamic pressure λu，l=167.8,whereas in the case of Ma=3.5, ΔT/ΔTcr=0,when the non-dimensional dynamic pressure λu，l＞ 91.4,panel flutter occurs,and the critical flutter dynamic pressure is decreased by 76.4.In the case of Ma=2.5, ΔT/ΔTcr=2,panel flutter occurs when the non-dimensional dynamic pressure is λu，l＞ 100.6, and in the case of Ma=3.5,ΔT/ΔTcr=2,panel flutter occurs when the non-dimensional dynamic pressure is λu，l＞ 54.8,and the critical flutter dynamic pressure is decreased by 45.8.The results indicate that the increment of the critical aerothermoelastic instability dynamic pressure decreases with the increase of temperature elevation when the temperature elevation is relatively small.

Fig.18 Stability boundaries with different Ma(on both surfaces).

Fig.18 shows the aeroelastic stability boundaries of the heated panel subjected to an impinging oblique shock with its inner and outer surfaces exposed to air flow （λd=3λu，l）with different Mach numbers.The trend of the effect of Mach number on the panel with its both surface subjected to air flow is similar to the effect of Mach number on the panel with its only one surface exposed to air flow.However,it should be mentioned that the panel will not be buckled for the case of ΔT/ΔTcr=0 when its inner and outer surfaces are subjected to the supersonic flow.It is also found that the aerothermoelastic instability boundaries of the panel with only one surface exposed to air flow are more sensitive to the change of Mach number than those of the panel with both surfaces subjected to air flow in the presence of an impinging oblique shock.

6.Conclusions

The aerothermoelastic stability of the heated panels subjected to an impinging oblique shock with only one surface exposed to air flow and with both surfaces subjected to air flow was investigated using theoretical analysis and numerical computations.The main conclusions are as follows:

(1)Panel is more prone to becoming unstable when it is subjected to an impinging oblique shock.Compared with the general case in the absence of shock,the special case studied in this paper shows that the heated panel becomes aeroelastically unstable at relatively small aerodynamic pressure.

(2)For the panel subjected to an impinging oblique shock with only one surface exposed to the air flow,only if the non-dimensionaldynamic pressureupstream of theshock impingement location and the non-dimensional dynamic pressure downstream of the shock impingement location of the panel satisfy critical condition of flutter stability,the heated panel will be aeroelastically stable.For the panel subjected to an impinging oblique shock with both surfaces exposed to the air flow,only if the nondimensional dynamic pressure upstream of the shock impingement location,the non-dimensional dynamic pressure downstream of the shock impingement location above the panel and the non-dimensional dynamic pressure below the panel satisfy critical condition of flutter stability,the heated panel will be aeroelastically stable.

(3)A new aeroelastic instability(distinct from regular panel flutter)arises from the complex interaction of the incident and reflected wave system with the panel flexural modes and thermal loads.For a relatively weak shock,the parameter space is still divided into four parts.For higher shock strength,the space will be divided into five parts but still four regions: flat and stable,transition region,buckling and basic flutter.Even for no thermal stress,the panel may buckle rather than only flutter.

(4)The dynamic pressure below the panel,the incident shock strength and the Mach number of the air flow all have a large effect on the aeroelastic stability boundaries of the panel.Further analysis will consider the effect of shock impingement location on the stability of the heated panel,as well as the three-dimensional panel.

Acknowledgements

The authors would like to thank Professor Earl Dowell for his help and suggestions regarding this paper and,in particular,for pointing out that by considering values of dynamic pressure and Mach number both ahead of and behind the shock and using these in turn in a standard flutter analysis that does not include the effects of a shock,one can obtain upper and lower bounds on the flutter boundary and limit cycle oscillations for the case with a shock.

The work is supported by the National Natural Science Foundation of China(No.11732013).

Appendix A

The non-dimensional variables are as follows: