# Aerothermoelastic Analysis of a Hypersonic Vehicle Based on Thermal Modal Reconstruction

Authors: Zhiqiang Chen; Yonghui Zhao
Journal: International Journal of Aerospace Engineering
 (2019)

---

## Abstract

Hypersonic vehicles operate in a severe aerodynamic heating environment, which has a significant impact on their structural dynamic characteristics. Therefore, aerodynamic heating effects cannot be ignored when performing an aeroelastic analysis for a hypersonic vehicle. However, incorporating aerodynamic heating effects into the fluid-structural coupling analysis will result in extreme computational costs. Actually, after experiencing a sustained flight in a fixed state, the vehicle will eventually reach the thermodynamic equilibrium. Thus, the aeroelastic analysis can be efficiently performed by using the structural dynamic characteristics of the heated vehicle operating in each equilibrium state. The effects of aerodynamic heating show that the modal frequencies and modal shapes of the flexible structure are bound to change significantly in comparison with the unheated structure. In this paper, a method of thermal modal reconstruction is developed in order to directly generate the structural mode shapes and frequencies within the given parameter space without having to solve a high-fidelity thermal and structural problem. Once the modal data are available, the multivariate interpolation in a tangent space to Grassmann manifold is used to generate the modal matrix at the arbitrary selected parameter point. Besides, the Kriging interpolation method is used to establish the approximate relationships between natural frequencies and sampling points. Finally, an example of an aerodynamic heated control surface structure is used to validate the effectiveness of the proposed aerothermoelastic framework. It is demonstrated that the developed thermal modal reconstruction method has good robustness, very high computational efficiency, and sufficient accuracy over a wide parametric domain.


---

## Body

## 1. Introduction

Hypersonic vehicles generally refer to the flight of vehicles experiencing a sustained long-range maneuverable flight in the atmosphere layer or trans-atmosphere layer at a Mach number above 5. This flight condition makes the aerothermodynamic environment extremely complex around hypersonic vehicles. A hypersonic vehicle, such as the X-43, usually adopts a long, slender lifting body layout, and the body and control surfaces are flexible due to minimum-weight restrictions. Thus, the wings, rudders, and other components will be subjected to urgent aerothermoelasticity [1–3]. The accurate and reliable predictions of aerothermodynamic load, structural temperature distribution, thermal deformation, and thermal stress as well as vibration response of thermal structural are the most important and challenging tasks.Due to the high speed of the hypersonic vehicle, severe aerodynamic heating occurs at its surface, and the aerodynamic heat flux is conducted into the internal structure. Aerodynamic thermal loads lead to the complexity and difficulty to the aerothermodynamic analysis of hypersonic vehicles. McNamara et al. [4] performed a systematic fluid-solid coupling study of the hypersonic aeroelastic and aerothermoelastic behavior of a three-dimensional configuration and concluded that the aeroelastic behavior of a vehicle is sensitive to structural variations caused by heating. Tabiei and Sockalingam [5] developed a multiphysics framework based on a loosely coupled strategy in combination with the computational fluid dynamics (CFD) code “Fluent” and the material thermal and structural response code “LS-DYNA.” However, the high computational cost of CFD and computational thermostructural dynamics (CTSD) regretfully make these approaches impractical for use in hypersonic aerothermoelasticity. Furthermore, high-speed high-enthalpy tunnels are not suitable for the aerothermoelastic testing of hypersonic vehicles at the moment [6]. Therefore, developing an efficient and accurate multiphysics framework for aerothermoelastic analysis is an urgent task.There are two methodologies to use in the generation of efficient aerodynamic heating models. The first approach is based on approximate models according to simplified assumptions. Eckert’s reference enthalpy [7, 8] and reference temperature methods have been used extensively in simple and efficient approximations for the aerodynamic heating. Culler and McNamara [9] used simplified coupling procedures associated with comprehensive aerothermoelastic analysis in order to abate the computational effort. An alternative approach is to use reduced-order models (ROMs) [10–12] that are derived from high-fidelity analysis tools. Chen et al. [13] provided a computationally efficient and accurate ROM framework based on surrogate and proper orthogonal decomposition (POD) for incorporating aerodynamic heating into a coupled aerothermoelastic analysis. Based on the basis augmentation and surrogate modeling techniques, Falkiewicz and Cesnik [14] proposed a reduced-order aerothermoelastic simulation framework to improve the computational efficiency.The high-temperature and large temperature gradient could not only degrade the material property of the structure but also produce thermal deformation and thermal stress inside the structure. The above thermal effects have great influence on the natural vibration characteristics of the vehicle. With the development of finite element methods and computer technology, the accurate numerical method has been widely applied in the thermal mode analysis [15, 16]. It is known that the aerothermal environment is strongly dependent on the flight parameters, Mach number, angle of attack, altitude etc. However, the thermal modal analysis based on the finite element method cannot adjust the structural modes to the changes of flight parameters in real time. In order to improve the flight parameter adaptability of thermal modal analysis, the structural mode shape database can be precomputed for the sampling points in the user-defined parameter space, and then the interpolation operation can be invoked for generating the new structural mode shapes at the arbitrary selected point in the user-defined parameter space. Amsallem and Farhat [17] proposed an interpolation in a tangent space to Grassmann manifold to construct a new reduced-order basis associated with a new set of parameters by interpolating reduced-order bases that were precomputed for selected values of these parameters. In this paper, the modal matrix corresponding to the arbitrary selected point in the user-defined parameter space is obtained by using the interpolation method. Each modal frequency is obtained through a Kriging surface [18] constructed by individual frequencies obtained from the high-fidelity simulations at the sampling points in the user-defined parameter space.The aim of this paper is to develop an efficient and accurate aerothermoelastic analysis framework in which the aerodynamic heating effects are considered. The above thermal modal reconstruction method is used to update the natural vibration characteristics of the heated structure. The third-order piston theory is employed for predicting the unsteady aerodynamic loads. Based on the proposed aerothermoelastic analysis framework, the aerothermoelastic behavior of the hypersonic vehicles can be easily and efficiently predicted at the arbitrary selected point in the whole user-defined parameter space. This paper is organized as follows: in Section2, the formula of structural dynamics under aerothermal environment are presented; in Section 3, the aerodynamic heating and the thermal modal analysis are performed using the high-fidelity method; details on thermal modal reconstruction are given in Section 4; the complete aerothermoelastic framework is given in Section 5; in Section 6, an illustrative simulation is given to demonstrate the effectiveness of the proposed thermal modal reconstruction strategy; and finally, conclusions are drawn in Section 7.


## 2. Formulation of Structural Dynamics under Aerothermal Environment

The full-order system of structural dynamic equations of motion in physical space is given by(1)Mu¨t+Kut=ft,where M and K are the mass matrix and stiffness matrix, respectively. ut represents the dynamic values of three-directional displacement, and ft is the force vector. The size of ut is three times the number of element nodes.According to the expansion theorem, the displacementut can be written as
(2)ut=∑i=1Nqitφi.Hence, Eq. (1) can be written in generalized formulation as
(3)Mgenq¨t+Kgenqt=fgent,where vector qt=q1t,q2t,⋯qNtT∈RN represents the modal or generalized coordinates and fgent is the generalized load vector. A FEM-based modal analysis leads to the truncated modal matrix Φ=φ1,φ2,⋯,φN with N considered eigenmodes, as well as generalized mass, stiffness matrices, and generalized load vector given by
(4)Mgen=ΦTMΦ,Kgen=ΦTKΦ,fgent=ΦTft.The modal matrixΦ is the mass orthonormal, and Mgen is equal to unity, and Kgen is a diagonal matrix in which main diagonal elements are natural frequencies ωi. The ith generalized force vector element in the time domain can be written as
(5)fgen,it=ρV∞2c22∫Scptφi⋅dS,where the pressure coefficient cpt obtained using the third-order piston theory is integrated over the body surface S and weighted with the ith eigenmode. The third-order piston theory [10] can be expressed as
(6)cp=2M∞2c1vna∞+c2vna∞2+c3vna∞3,where
(7)c1=M∞δ,c2=M∞4γ+1−4δ24δ4,c3=16M∞δ7a0M∞8+b0M∞6+c0M∞4+d0M∞2+e0,a0=3γ+142,b0=3γ2−12γ−74,c0=9γ+12,d0=−6,e0=4,δ2=M∞2−1,where M∞ is the free stream Mach number, a∞ is the speed of sound, γ is the specific heat ratio of air, and vn is the normal velocity of the surface.It is known that the hypersonic vehicle must survive the harsh thermal environment when experiencing a sustained long-range maneuverable flight. The structural stiffness matrix is significantly influenced by the material properties and thermal stresses resulting from the aerodynamic heating effects. Therefore, the aerodynamic heating effect will result in the significant change of the natural characteristics of the structure, such as mode shapeφi and natural frequencies ωi. Note that the effects of material property changes and thermal stresses are both included in the calculation of the mode shapes and natural frequencies of the heated structure [4].


## 3. Hypersonic Aerodynamic Heating and Thermal Modal Analysis

Aerodynamic heating in the boundary layer is obtained using the CFD method. According to the experience of previous studies [4, 6, 19], the surface temperatures are selected as the peak temperatures of the thermal equilibrium state, which occur when the aerodynamic heat flux q̇aero is in equilibrium with the radiation heat flux q̇rad and the conductive heat flux q̇cond:
(8)q̇aero=q̇rad+q̇cond,q̇aero=hhtTAW−TW,where TAW is the adiabatic wall temperature and TW is the wall temperature. After ignoring the heat conduction flux into the structure, the radiation equilibrium temperature TR can be determined by solving a fourth-order algebraic equation:
(9)q̇aero=q̇rad,hhtTAW−TR=σεTR4,where σ is the Boltzmann constant, equal to 1.3806505×10−23J/K, and ε is the emissivity of the surface. The adiabatic-wall temperature TAW and the heat transfer coefficient hht can be calculated separately using the CFD technique by specifying adiabatic wall and constant wall temperature as the boundary condition for the Navier-Stokes equations. These quantities are used to define the boundary conditions of the heat transfer problem in the FEM solver. Afterwards, a simple steady-state heat transfer analysis is performed using the “Steady State Analysis Solution” (Sol 153) in MSC Nastran with the selected “HEAT” option. The mode shapes and the natural frequencies of the heated structure are calculated using the “Nonlinear Statics Solution” (Sol 106) in MSC Nastran with the “NORMAL MODES” option selected.It is known that the changes in natural frequencies and mode shapes have a great impact on the aeroelastic behavior of the hypersonic vehicle. The overall goal of this paper is to directly generate the structural mode shapes and frequencies without having to solve a high-fidelity thermal and structural problem.


## 4. Thermal Modal Reconstruction

In order to obtain the modal information in the user-defined parameter space, the structural mode shapes and frequencies are calculated at the corresponding sample points using the method introduced in Section3. A novel method of thermal modal reconstruction is proposed using two advanced interpolation methods: mode shape vectors are obtained using multivariate interpolation in a tangent space to Grassmann manifold, and natural frequencies are obtained using Kriging interpolation at the arbitrary selected operation point in the user-defined parameter space.


### 4.1. Multivariate Interpolation in a Tangent Space to Grassmann Manifold

Flight parameters such as Mach number, altitude, and angle of attack have great influence on the magnitude and distribution of temperature due to aerodynamic heating. Letxi=x1i,x2i,⋯,xNpii=1,2,…,Ns denote a specific configuration of the set of Np parameters, and Ns is the number of sampling points selected on the considered parameter space. Thus, the corresponding modal matrix Φxi is calculated at xi flight condition. In this paper, multivariate interpolation in a tangent space to Grassmann manifold [20] is applied to generate a modal matrix Φxn at a new sampling point xn by the precomputed modal matrix. For a detailed discussion of the multivariate interpolation method, see reference [17]. Here, only the most relevant properties and operating steps are discussed in the following.Grassmann manifold is a generalized concept of the projective space, and the set of allN-dimensional subspaces of RMN<M is called the Grassmann manifold which is denoted as GN,M. So each N-dimensional subspace S of can be considered as a point of the manifold GN,M. Here, N, the dimension of the truncated modal matrix Φxi∈RM×N, is fixed independently from the operating parameters. The subspace Si spanned by Φxi belongs to the Grassmann manifold GN,M. Because the manifold GN,M is a Riemannian manifold [20], there is a tangent space TS that is a vector space at each point Si. Therefore, the representative matrix Γ of a point in TS can be interpolated directly. A new subspace Sn associated with the new operating parameters xn can be obtained by interpolating the known subspaces Sii=1n−1 in four steps as follows.
(1)
Choose a reference pointSxi0∈GM,N,i0∈1,2,⋯,n−1, which is treated as the origin point for the interpolation
(2)
Each pointSxi sufficiently close to Sxi0 is mapped to a matrix Γi representing a point of TSxi0 using the logarithm map and based on thin singular value decomposition (SVD). This can be written as
(10)I−Φxi0ΦTxi0ΦxiΦTxi0Φxi−1=UiΣiViTthinSVD,Γi=Uitan−1ΣiViT,where Σi is a diagonal matrix containing the singular values, and Ui and Vi are orthogonal matrices.
(3)
Compute the matrixΓn for the target operating point xn by interpolating the corresponding entries of the matrix Γi. For this purpose, use a multivariate interpolation scheme [21]
(4)
Compute the exponential ofΓn to obtain the desired modal matrix Φxn. This can be written as
(11)Γn=UnΣnVnT,Φxi=Φxi0VncosΣn+UnsinΣn.


### 4.2. Kriging Interpolation of Natural Frequencies

Based on the high-fidelity system, a series of numerical simulations are performed at a set of sampling pointsxi=x1i,x2i,⋯,xNpi, i=1,2,…,Ns, in order to produce Ns modal matrices and natural frequencies Fi=f1i,f2i,⋯,fNi. In a multidimensional parametric space, a nonlinear functional relationship between natural frequencies and design variables is expressed by
(12)fr=yrx1,x1,⋯,xNp,r=1,2,…,N,where fr is the natural frequency of the rth mode shape and N is the total number of modes retained after truncation. In this paper, the Kriging interpolation method, due to its outstanding ability to model the local behavior of the function, is used to establish the approximate relationships between natural frequencies and sampling points.Consider the case ofNs sampling points xi,i=1,2,⋯,Ns in the Np dimensional parameter space. The corresponding responses are denoted by fri,i=1,2,⋯,Ns. The Kriging approximation is expressed at an arbitrary point x as
(13)yrx=hTxβ+Zx,where hx are polynomial functions, β are the corresponding parameters, and Zx is a stochastic process. The reader is referred to reference [10] for the details of the Kriging interpolation method. Then, we can follow the reference to build N Kriging surfaces based on the sampled data set xi,fri, i=1,2,…,Ns. Therefore, the prediction of yr at an arbitrary point x on parametric space can be shown to be
(14)yrx=hTxβ^+rTxR−1y−Hβ^,where y are the outputs from sampled values, H is the vector of hx at the sampling points, R is the correlation function matrix, rx is the correlation between unknown points x and the known sampling points, and β^ is the least squares predictor given by
(15)β^=HR−1H−1HR−1y.


## 4.1. Multivariate Interpolation in a Tangent Space to Grassmann Manifold

Flight parameters such as Mach number, altitude, and angle of attack have great influence on the magnitude and distribution of temperature due to aerodynamic heating. Letxi=x1i,x2i,⋯,xNpii=1,2,…,Ns denote a specific configuration of the set of Np parameters, and Ns is the number of sampling points selected on the considered parameter space. Thus, the corresponding modal matrix Φxi is calculated at xi flight condition. In this paper, multivariate interpolation in a tangent space to Grassmann manifold [20] is applied to generate a modal matrix Φxn at a new sampling point xn by the precomputed modal matrix. For a detailed discussion of the multivariate interpolation method, see reference [17]. Here, only the most relevant properties and operating steps are discussed in the following.Grassmann manifold is a generalized concept of the projective space, and the set of allN-dimensional subspaces of RMN<M is called the Grassmann manifold which is denoted as GN,M. So each N-dimensional subspace S of can be considered as a point of the manifold GN,M. Here, N, the dimension of the truncated modal matrix Φxi∈RM×N, is fixed independently from the operating parameters. The subspace Si spanned by Φxi belongs to the Grassmann manifold GN,M. Because the manifold GN,M is a Riemannian manifold [20], there is a tangent space TS that is a vector space at each point Si. Therefore, the representative matrix Γ of a point in TS can be interpolated directly. A new subspace Sn associated with the new operating parameters xn can be obtained by interpolating the known subspaces Sii=1n−1 in four steps as follows.
(1)
Choose a reference pointSxi0∈GM,N,i0∈1,2,⋯,n−1, which is treated as the origin point for the interpolation
(2)
Each pointSxi sufficiently close to Sxi0 is mapped to a matrix Γi representing a point of TSxi0 using the logarithm map and based on thin singular value decomposition (SVD). This can be written as
(10)I−Φxi0ΦTxi0ΦxiΦTxi0Φxi−1=UiΣiViTthinSVD,Γi=Uitan−1ΣiViT,where Σi is a diagonal matrix containing the singular values, and Ui and Vi are orthogonal matrices.
(3)
Compute the matrixΓn for the target operating point xn by interpolating the corresponding entries of the matrix Γi. For this purpose, use a multivariate interpolation scheme [21]
(4)
Compute the exponential ofΓn to obtain the desired modal matrix Φxn. This can be written as
(11)Γn=UnΣnVnT,Φxi=Φxi0VncosΣn+UnsinΣn.


## 4.2. Kriging Interpolation of Natural Frequencies

Based on the high-fidelity system, a series of numerical simulations are performed at a set of sampling pointsxi=x1i,x2i,⋯,xNpi, i=1,2,…,Ns, in order to produce Ns modal matrices and natural frequencies Fi=f1i,f2i,⋯,fNi. In a multidimensional parametric space, a nonlinear functional relationship between natural frequencies and design variables is expressed by
(12)fr=yrx1,x1,⋯,xNp,r=1,2,…,N,where fr is the natural frequency of the rth mode shape and N is the total number of modes retained after truncation. In this paper, the Kriging interpolation method, due to its outstanding ability to model the local behavior of the function, is used to establish the approximate relationships between natural frequencies and sampling points.Consider the case ofNs sampling points xi,i=1,2,⋯,Ns in the Np dimensional parameter space. The corresponding responses are denoted by fri,i=1,2,⋯,Ns. The Kriging approximation is expressed at an arbitrary point x as
(13)yrx=hTxβ+Zx,where hx are polynomial functions, β are the corresponding parameters, and Zx is a stochastic process. The reader is referred to reference [10] for the details of the Kriging interpolation method. Then, we can follow the reference to build N Kriging surfaces based on the sampled data set xi,fri, i=1,2,…,Ns. Therefore, the prediction of yr at an arbitrary point x on parametric space can be shown to be
(14)yrx=hTxβ^+rTxR−1y−Hβ^,where y are the outputs from sampled values, H is the vector of hx at the sampling points, R is the correlation function matrix, rx is the correlation between unknown points x and the known sampling points, and β^ is the least squares predictor given by
(15)β^=HR−1H−1HR−1y.


## 5. Overview of the Aerothermoelastic Framework

In this section, the theoretical approaches described above are combined in order to obtain an efficient and accurate aerothermoelastic framework. The complete aerothermoelastic framework is illustrated by the flowchart in Figure1. Firstly, given the control surface configuration and the bounds of the parameter space, the uniformly distributed sampling points are produced in the parameter space using the successive local enumeration algorithm [10], which is an enhanced algorithm of the Latin hypercube sampling. At each sampling point, the aerodynamic heating solutions on the surface of structure obtained by the CFD solver are computed. Secondly, a steady-state heat transfer analysis is performed using MSC Nastran (Sol 153). The mode shapes and frequencies of the heated structure are determined using MSC Nastran (Sol 106). Thirdly, after collecting the modal matrix data, multivariate interpolation in a tangent space to Grassmann manifold is applied to generate a modal matrix at the new parameter point. Fifthly, the Kriging interpolation method is used to establish the approximate relationships between natural frequencies and sampling points. Finally, unsteady aerodynamic forces obtained by the third-order piston theory are coupled with structural modes to predict the aerothermoelastic behavior of hypersonic vehicles in the considered parameter space.Figure 1
Schematic overview of the aerothermoelastic framework.


## 6. Results and Discussion

### 6.1. Structural and Aerodynamic Models of the Control Surface

To validate the developed aerothermoelastic framework, a representative control surface structure flying in the hypersonic environment is considered. A detailed description of the structure is given in reference [22, 23]. The planform geometric and cross-sectional views of the wing are shown in Figure 2. The top and bottom skin layers, respectively, include two 3.8 mm-thick thermal protection system layers, totally 15.2 mm of the thermal protection system material. The thickness of all stiffeners is 25.4 mm. The thermal protection system consists of an outer heat shield and middle insulation layer on top of the skin, as shown in Figure 3 of reference [23]. The material of the heat shield is chosen to be René 41, which is efficient in terms of mechanical properties at high temperatures. The emissivity of the heat shield is set to be 0.85. For the insulation layer, Min-K is chosen due to its relatively low thermal diffusivity. For both skin and stiffeners of the structure, the titanium alloy TIMETAL 834 is chosen. The main thermal and mechanical property indexes of the three materials are shown in Table 1 of reference [23], where “T-dep.” indicates that the property index is a temperature-dependent parameter [24–26], and “Neglect” indicates that the property is neglected in the structural model. The finite element model used for the thermal and structural modeling aspects of the study is shown in Figure 3. The model contains 2812 degrees of freedom and 6886 elements. The heat shield and insulation are each modeled using one layer of six-node solid wedge elements. The top and bottom skins and stiffeners are modeled using three-node 2D triangular elements. There are 3456 solid elements and 3430 triangular elements in the model. In addition, rigid elements (RBE2) are used between each skin node and the corresponding node at the outer surface of the insulation. Finally, it is important to point out that the rigid-body modes are suppressed by fixing the wing at the root.Figure 2
Geometric model of the control surface.Table 1
MAC value associated with each mode shape.
Mode
MAC forxa5.8,6
MAC forxb7.4,10
MAC forxc9.2,14
1
0.9987
0.9992
0.9986
2
0.9976
0.9981
0.9971
3
0.9867
0.9908
0.9875
4
0.9824
0.9854
0.9813
5
0.9715
0.9726
0.9676Figure 3
Kriging surfaces of modal frequencies.
(a)
Mode 1
(b)
Mode 2
(c)
Mode 3
(d)
Mode 4
(e)
Mode 5The CFD solver uses an implicit finite volume algorithm, and the upwind-biased spatial differencing scheme is used to solve the time-dependent Reynolds-averaged Navier-Stokes (RANS) equations. The Spalart-Allmaras turbulence model is used for closure of the RANS equations. As illustrated in Figure4, the CFD mesh is a vertically symmetric H-H grid with 53 points spanwise, 153 points chordwise, and 43 points extending vertically from the surface (approximately 6.1×105 cells). We can see that the majority of grid points are clustered near the wing surface, leading edge, and mid-chord since these locations correspond to the maximum flow gradients and the boundary layers in hypersonic flows are relatively thick. This grid does not include any sections upstream of the wing surface since the considered flow is hypersonic and disturbances cannot propagate upstream.Figure 4
Finite element model of control surface.


### 6.2. Thermal Modal Analysis Using High-Fidelity Analysis Tools

In this section, two flight parameters (free stream Mach number5≤M∞≤10 and altitude 5km≤H≤15km) are considered, which have great influence on the magnitude and distribution of temperature due to aerodynamic heating. The angle of attack is taken as 0.0 deg for all simulations. The 36 flight conditions (sampling points) xi=M∞i,Hi, i=1,2,…,36 are artificially chosen over the range of the 2D parameter space, as shown in Figure 5. Under sustained aerothermal environment, the leading edge of the hypersonic wing will be subjected to the higher heat fluxes; hence, the stronger fluid-thermal-structural coupling phenomena will be expected. This is due to the fact that the flow is brought to rest at the leading edge, thus the kinetic energy of the flow is converted into heat. The external aeroheating flux gradually transfers to the interior of the wing, and then the structural temperature rises increasingly until a steady state is reached. In this paper, the aerodynamic heating of the wing is calculated using an aero-thermo coupling technique (see Section 3). When the aerothermodynamic loop has converged, the temperature distribution of the wing is obtained, as shown in Figure 6 for M∞=7 and H=10km. It can be seen that the peak temperatures occur at the leading edge with a maximum temperature of 1, 276 K, and the temperature gradient of the leading edge is obviously larger than any other part of the wing. Spain et al. [27] have pointed out the negative impact of high temperatures due to the material degradation and thermal stresses created by temperature gradients.Figure 5
Grids for full-order CFD computations.Figure 6
Operating domain (the red circles represent sampling points).In reference [4], the aeroelastic and aerothermoelastic analysis of the hypersonic wing was carried out using 5 structural modes. In Figure 7, case 1 presents the first five vibration modes and frequencies of the wing without considering aerodynamic heating effect. Case 2 presents the first five vibration modes and frequencies of the heated wing at M∞=7 and H=10km. It can be seen that the first four modal shapes are almost unchanged, while the fifth changes significantly in modal shapes with a trend towards bending-torsion coupling. It can be inferred that aerodynamic heating affects the modal shapes of higher order more easily for the wings. Table 2 shows the difference value (D-value) Δf of the first five natural frequencies between case 1 and case 2. From Table 1, we can see that the natural frequency of different modes tends to decrease, and the higher the mode order is, the more rapidly its frequency falls. Obviously, aerodynamic heating has a great effect on the natural vibration characteristics of the structure.Figure 7
Temperature distribution of the radiation equilibrium at the wall (K) atM∞=7 and H=10km.Table 2
The difference of the first five natural frequencies between case 1 and case 2.
D-value
Mode 1
Mode 2
Mode 3
Mode 4
Mode 5
Δf(Hz)
1.51
3.26
6.74
9.67
13.53


### 6.3. Thermal Modal Reconstruction

The purpose of this part is to examine the accuracy of the structural mode shapes and natural frequencies obtained using multivariate interpolation in a tangent space to Grassmann manifold and Kriging interpolation method. Firstly, a truncated modal matrixΦ5M∞i,Hi, i=1,2,…36 is constructed by the first five mode shape vectors at each sampling point. For the modal matrix at any point of the parameter space, the interpolation algorithm is applied in the tangent space to the Grassmann manifold at the chosen reference point. Since the dimensionless linear distance between M∞1,H1 and test point M∞∗,H∗ is the shortest, the sampling point M∞1,H1 is chosen as a reference and origin point, as shown in Figure 8. Here, we use the standard 5-point bivariate interpolation scheme to compute each entry of the matrix Γ∗ at five cross-shaped points M∞i,Hi,i=1,2,⋯,5.Figure 8
The first six vibration modes and frequencies of the wing, (1) case 1: without the thermal effects, (2) case 2: with the thermal effects.
(a)
Case 1: mode 1,f1=23.83Hz
(b)
Case 2: mode 1,f1=22.32Hz
(c)
Case 1: mode 2,f2=40.72Hz
(d)
Case 2: mode 2,f2=37.46Hz
(e)
Case 1: mode 3,f3=71.45Hz
(f)
Case 2: mode 3,f3=64.71Hz
(g)
Case 1: mode 4,f4=89.82Hz
(h)
Case 2: mode 4,f4=80.15Hz
(i)
Case 1: mode 5,f5=107.27Hz
(j)
Case 2: mode 5,f5=93.74HzHere, the accuracy of the restructured thermal mode is validated at three arbitrary test points,xa5.8,10, xb7.4,6, and xc9.2,14. Note that locations of these test points are not coincident with those of the sampling points. The accuracy of thermal mode restructured is assessed by comparing the natural frequencies and mode shapes with those computed using the high-fidelity method. For the mode shape vectors, the comparison is performed using the modal assurance criterion (MAC) [28]. For two modes φi and φj, the MAC value is given as
(16)MCA=φjTφj2φjTφiφjTφj.The value of MAC is limited between zero (representing inconsistent mode shapes) and one (representing fully consistent mode shapes); it can be used to compare the modal vectors quantitatively. The results listed in Table1 show that the mode shape vectors obtained by the interpolation method are in good agreement with those computed using the high-fidelity method at the three arbitrary test points.After the natural frequencies are obtained at the sampling points, Kriging surfaces of the first five modal frequencies are calculated, as shown in Figure9. It is shown that the natural frequencies of different modes of the wing tend to decrease with the increasing Mach number and the decreasing altitude, and the natural frequencies of the higher order fall more rapidly. Natural frequencies of high-order modes get nearer one another more easily, and that leads to bending-torsion coupling. Table 3 shows the comparison of the first five modal frequencies at the three arbitrary test points obtained by Kriging interpolation with their counterparts obtained from the high-fidelity method. These results illustrate that the error of the proposed thermal modal reconstruction method is below 4%.Figure 9
Selection of interpolation points.Table 3
Comparisons of the first five modal frequencies at the three arbitrary test points.
Test pointxa5.8,6
Test pointxb7.4,10
Test pointxc9.2,14
Mode
High-fidelity method
Kriging interpolation
Relative error (%)
High-fidelity method
Kriging interpolation
Relative error (%)
High-fidelity method
Kriging interpolation
Relative error (%)
1
22.74
22.82
3.5
22.47
22.53
2.7
22.31
22.23
3.6
2
39.42
39.34
2.1
38.5
38.41
2.3
36.42
36.31
3
3
68.51
68.42
1.3
66.67
66.54
1.9
63.88
63.72
2.4
4
84.49
84.65
1.9
82.34
82.19
1.8
78.29
78.12
2.1
5
98.63
98.47
1.6
95.37
95.15
2.3
91.12
90.94
1.9


### 6.4. Aerothermoelastic Behavior of the Wing

In this section, the aeroelastic and aerothermoelastic analysis of the wing is carried out using the unsteady aerodynamics of the third-order piston theory in hypersonic flow. At the constant altitude, the flutter point is found by increasing the Mach number until the net damping of the system becomes zero. Figures10 and 11 show time histories and the FFT spectrum of wing tip displacement at two different Mach numbers. In Figure 10, the response of the heated wing exhibits harmonic oscillations at M∞=6.13, but the response of the unheated wing is convergent. In Figure 11, the response of the unheated wing is stable at M∞=6.31, while the response of the heated wing is divergent. It can be seen that the flutter frequency decreases from 34.3 Hz to 32.5 Hz due to the aerodynamic heating effect. It is concluded that the aerodynamic heating effect can result in earlier transition to flutter. Besides, the effect of increasing altitude on the flutter boundary is also examined, as shown in Figure 12, where the results of aeroelastic and aerothermoelastic simulations are presented. We can see that the flutter boundaries of the heated wing are smaller than those of the unheated wing at the same altitude; the higher the altitude is, the larger the flutter speed is. Obviously, the aerodynamic heating effect leads to the difference of the flutter boundary. These simulation results indicate that the severe aerodynamic heating effect significantly influences the aeroelastic stability of the hypersonic vehicle through a combination of the degradation of material properties and the induced thermal stresses.Figure 10
Time histories and FFT spectrum of wing tip displacement atM∞=6.13.
(a)
(b)Figure 11
Time histories and FFT spectrum of wing tip displacement atM∞=6.31.
(a)
(b)Figure 12
Flutter boundaries of the unheated and heated wing.A comparison is made for the computational cost to demonstrate the efficiency of the proposed aerothermoelastic framework. Firstly, aerothermal environment is predicted using the CFD method, and the computation time is found to be approximately 4.3 h. Then, the corresponding CFD results are transferred to the FEM solver to obtain the thermal modes. Finally, the aerothermoelastic responses are predicted by combining the third-order piston theory with the structural modes obtained from the thermal modal reconstruction. In this paper, when the simulation time-step size is set to be 2000, the total computation time is about 4.4 h. However, when the flight parameters change, we have to repeat the above process. Obviously, it is very tedious and time-consuming. In our developed method, once thermal modal reconstruction is constructed, the aerothermoelastic analysis of the hypersonic vehicle can be efficiently performed at the arbitrary point in the considered parameter space. For the same case, the computational time using the developed aerothermoelastic framework is found to be less than 2 minutes. Thus, a very high computational efficiency is achieved.


## 6.1. Structural and Aerodynamic Models of the Control Surface

To validate the developed aerothermoelastic framework, a representative control surface structure flying in the hypersonic environment is considered. A detailed description of the structure is given in reference [22, 23]. The planform geometric and cross-sectional views of the wing are shown in Figure 2. The top and bottom skin layers, respectively, include two 3.8 mm-thick thermal protection system layers, totally 15.2 mm of the thermal protection system material. The thickness of all stiffeners is 25.4 mm. The thermal protection system consists of an outer heat shield and middle insulation layer on top of the skin, as shown in Figure 3 of reference [23]. The material of the heat shield is chosen to be René 41, which is efficient in terms of mechanical properties at high temperatures. The emissivity of the heat shield is set to be 0.85. For the insulation layer, Min-K is chosen due to its relatively low thermal diffusivity. For both skin and stiffeners of the structure, the titanium alloy TIMETAL 834 is chosen. The main thermal and mechanical property indexes of the three materials are shown in Table 1 of reference [23], where “T-dep.” indicates that the property index is a temperature-dependent parameter [24–26], and “Neglect” indicates that the property is neglected in the structural model. The finite element model used for the thermal and structural modeling aspects of the study is shown in Figure 3. The model contains 2812 degrees of freedom and 6886 elements. The heat shield and insulation are each modeled using one layer of six-node solid wedge elements. The top and bottom skins and stiffeners are modeled using three-node 2D triangular elements. There are 3456 solid elements and 3430 triangular elements in the model. In addition, rigid elements (RBE2) are used between each skin node and the corresponding node at the outer surface of the insulation. Finally, it is important to point out that the rigid-body modes are suppressed by fixing the wing at the root.Figure 2
Geometric model of the control surface.Table 1
MAC value associated with each mode shape.
Mode
MAC forxa5.8,6
MAC forxb7.4,10
MAC forxc9.2,14
1
0.9987
0.9992
0.9986
2
0.9976
0.9981
0.9971
3
0.9867
0.9908
0.9875
4
0.9824
0.9854
0.9813
5
0.9715
0.9726
0.9676Figure 3
Kriging surfaces of modal frequencies.
(a)
Mode 1
(b)
Mode 2
(c)
Mode 3
(d)
Mode 4
(e)
Mode 5The CFD solver uses an implicit finite volume algorithm, and the upwind-biased spatial differencing scheme is used to solve the time-dependent Reynolds-averaged Navier-Stokes (RANS) equations. The Spalart-Allmaras turbulence model is used for closure of the RANS equations. As illustrated in Figure4, the CFD mesh is a vertically symmetric H-H grid with 53 points spanwise, 153 points chordwise, and 43 points extending vertically from the surface (approximately 6.1×105 cells). We can see that the majority of grid points are clustered near the wing surface, leading edge, and mid-chord since these locations correspond to the maximum flow gradients and the boundary layers in hypersonic flows are relatively thick. This grid does not include any sections upstream of the wing surface since the considered flow is hypersonic and disturbances cannot propagate upstream.Figure 4
Finite element model of control surface.


## 6.2. Thermal Modal Analysis Using High-Fidelity Analysis Tools

In this section, two flight parameters (free stream Mach number5≤M∞≤10 and altitude 5km≤H≤15km) are considered, which have great influence on the magnitude and distribution of temperature due to aerodynamic heating. The angle of attack is taken as 0.0 deg for all simulations. The 36 flight conditions (sampling points) xi=M∞i,Hi, i=1,2,…,36 are artificially chosen over the range of the 2D parameter space, as shown in Figure 5. Under sustained aerothermal environment, the leading edge of the hypersonic wing will be subjected to the higher heat fluxes; hence, the stronger fluid-thermal-structural coupling phenomena will be expected. This is due to the fact that the flow is brought to rest at the leading edge, thus the kinetic energy of the flow is converted into heat. The external aeroheating flux gradually transfers to the interior of the wing, and then the structural temperature rises increasingly until a steady state is reached. In this paper, the aerodynamic heating of the wing is calculated using an aero-thermo coupling technique (see Section 3). When the aerothermodynamic loop has converged, the temperature distribution of the wing is obtained, as shown in Figure 6 for M∞=7 and H=10km. It can be seen that the peak temperatures occur at the leading edge with a maximum temperature of 1, 276 K, and the temperature gradient of the leading edge is obviously larger than any other part of the wing. Spain et al. [27] have pointed out the negative impact of high temperatures due to the material degradation and thermal stresses created by temperature gradients.Figure 5
Grids for full-order CFD computations.Figure 6
Operating domain (the red circles represent sampling points).In reference [4], the aeroelastic and aerothermoelastic analysis of the hypersonic wing was carried out using 5 structural modes. In Figure 7, case 1 presents the first five vibration modes and frequencies of the wing without considering aerodynamic heating effect. Case 2 presents the first five vibration modes and frequencies of the heated wing at M∞=7 and H=10km. It can be seen that the first four modal shapes are almost unchanged, while the fifth changes significantly in modal shapes with a trend towards bending-torsion coupling. It can be inferred that aerodynamic heating affects the modal shapes of higher order more easily for the wings. Table 2 shows the difference value (D-value) Δf of the first five natural frequencies between case 1 and case 2. From Table 1, we can see that the natural frequency of different modes tends to decrease, and the higher the mode order is, the more rapidly its frequency falls. Obviously, aerodynamic heating has a great effect on the natural vibration characteristics of the structure.Figure 7
Temperature distribution of the radiation equilibrium at the wall (K) atM∞=7 and H=10km.Table 2
The difference of the first five natural frequencies between case 1 and case 2.
D-value
Mode 1
Mode 2
Mode 3
Mode 4
Mode 5
Δf(Hz)
1.51
3.26
6.74
9.67
13.53


## 6.3. Thermal Modal Reconstruction

The purpose of this part is to examine the accuracy of the structural mode shapes and natural frequencies obtained using multivariate interpolation in a tangent space to Grassmann manifold and Kriging interpolation method. Firstly, a truncated modal matrixΦ5M∞i,Hi, i=1,2,…36 is constructed by the first five mode shape vectors at each sampling point. For the modal matrix at any point of the parameter space, the interpolation algorithm is applied in the tangent space to the Grassmann manifold at the chosen reference point. Since the dimensionless linear distance between M∞1,H1 and test point M∞∗,H∗ is the shortest, the sampling point M∞1,H1 is chosen as a reference and origin point, as shown in Figure 8. Here, we use the standard 5-point bivariate interpolation scheme to compute each entry of the matrix Γ∗ at five cross-shaped points M∞i,Hi,i=1,2,⋯,5.Figure 8
The first six vibration modes and frequencies of the wing, (1) case 1: without the thermal effects, (2) case 2: with the thermal effects.
(a)
Case 1: mode 1,f1=23.83Hz
(b)
Case 2: mode 1,f1=22.32Hz
(c)
Case 1: mode 2,f2=40.72Hz
(d)
Case 2: mode 2,f2=37.46Hz
(e)
Case 1: mode 3,f3=71.45Hz
(f)
Case 2: mode 3,f3=64.71Hz
(g)
Case 1: mode 4,f4=89.82Hz
(h)
Case 2: mode 4,f4=80.15Hz
(i)
Case 1: mode 5,f5=107.27Hz
(j)
Case 2: mode 5,f5=93.74HzHere, the accuracy of the restructured thermal mode is validated at three arbitrary test points,xa5.8,10, xb7.4,6, and xc9.2,14. Note that locations of these test points are not coincident with those of the sampling points. The accuracy of thermal mode restructured is assessed by comparing the natural frequencies and mode shapes with those computed using the high-fidelity method. For the mode shape vectors, the comparison is performed using the modal assurance criterion (MAC) [28]. For two modes φi and φj, the MAC value is given as
(16)MCA=φjTφj2φjTφiφjTφj.The value of MAC is limited between zero (representing inconsistent mode shapes) and one (representing fully consistent mode shapes); it can be used to compare the modal vectors quantitatively. The results listed in Table1 show that the mode shape vectors obtained by the interpolation method are in good agreement with those computed using the high-fidelity method at the three arbitrary test points.After the natural frequencies are obtained at the sampling points, Kriging surfaces of the first five modal frequencies are calculated, as shown in Figure9. It is shown that the natural frequencies of different modes of the wing tend to decrease with the increasing Mach number and the decreasing altitude, and the natural frequencies of the higher order fall more rapidly. Natural frequencies of high-order modes get nearer one another more easily, and that leads to bending-torsion coupling. Table 3 shows the comparison of the first five modal frequencies at the three arbitrary test points obtained by Kriging interpolation with their counterparts obtained from the high-fidelity method. These results illustrate that the error of the proposed thermal modal reconstruction method is below 4%.Figure 9
Selection of interpolation points.Table 3
Comparisons of the first five modal frequencies at the three arbitrary test points.
Test pointxa5.8,6
Test pointxb7.4,10
Test pointxc9.2,14
Mode
High-fidelity method
Kriging interpolation
Relative error (%)
High-fidelity method
Kriging interpolation
Relative error (%)
High-fidelity method
Kriging interpolation
Relative error (%)
1
22.74
22.82
3.5
22.47
22.53
2.7
22.31
22.23
3.6
2
39.42
39.34
2.1
38.5
38.41
2.3
36.42
36.31
3
3
68.51
68.42
1.3
66.67
66.54
1.9
63.88
63.72
2.4
4
84.49
84.65
1.9
82.34
82.19
1.8
78.29
78.12
2.1
5
98.63
98.47
1.6
95.37
95.15
2.3
91.12
90.94
1.9


## 6.4. Aerothermoelastic Behavior of the Wing

In this section, the aeroelastic and aerothermoelastic analysis of the wing is carried out using the unsteady aerodynamics of the third-order piston theory in hypersonic flow. At the constant altitude, the flutter point is found by increasing the Mach number until the net damping of the system becomes zero. Figures10 and 11 show time histories and the FFT spectrum of wing tip displacement at two different Mach numbers. In Figure 10, the response of the heated wing exhibits harmonic oscillations at M∞=6.13, but the response of the unheated wing is convergent. In Figure 11, the response of the unheated wing is stable at M∞=6.31, while the response of the heated wing is divergent. It can be seen that the flutter frequency decreases from 34.3 Hz to 32.5 Hz due to the aerodynamic heating effect. It is concluded that the aerodynamic heating effect can result in earlier transition to flutter. Besides, the effect of increasing altitude on the flutter boundary is also examined, as shown in Figure 12, where the results of aeroelastic and aerothermoelastic simulations are presented. We can see that the flutter boundaries of the heated wing are smaller than those of the unheated wing at the same altitude; the higher the altitude is, the larger the flutter speed is. Obviously, the aerodynamic heating effect leads to the difference of the flutter boundary. These simulation results indicate that the severe aerodynamic heating effect significantly influences the aeroelastic stability of the hypersonic vehicle through a combination of the degradation of material properties and the induced thermal stresses.Figure 10
Time histories and FFT spectrum of wing tip displacement atM∞=6.13.
(a)
(b)Figure 11
Time histories and FFT spectrum of wing tip displacement atM∞=6.31.
(a)
(b)Figure 12
Flutter boundaries of the unheated and heated wing.A comparison is made for the computational cost to demonstrate the efficiency of the proposed aerothermoelastic framework. Firstly, aerothermal environment is predicted using the CFD method, and the computation time is found to be approximately 4.3 h. Then, the corresponding CFD results are transferred to the FEM solver to obtain the thermal modes. Finally, the aerothermoelastic responses are predicted by combining the third-order piston theory with the structural modes obtained from the thermal modal reconstruction. In this paper, when the simulation time-step size is set to be 2000, the total computation time is about 4.4 h. However, when the flight parameters change, we have to repeat the above process. Obviously, it is very tedious and time-consuming. In our developed method, once thermal modal reconstruction is constructed, the aerothermoelastic analysis of the hypersonic vehicle can be efficiently performed at the arbitrary point in the considered parameter space. For the same case, the computational time using the developed aerothermoelastic framework is found to be less than 2 minutes. Thus, a very high computational efficiency is achieved.


## 7. Conclusion

In order to efficiently predict aerothermoelastic behavior of the hypersonic vehicle, we develop an aerothermoelastic framework based on thermal modal reconstruction and piston theory. First, the thermal modal data are obtained beforehand at each sampling point. Then, new mode shapes and natural frequencies are generated based on the proposed thermal modal reconstruction method in the prescribed parameter space. Finally, the aerothermoelastic analysis is carried out using the third-order piston theory aerodynamics in hypersonic flow. The principal conclusions are summarized as follows:(1)
The natural frequencies of different orders for the heated wing tend to decrease with increasing Mach number and decreasing altitude, the natural frequencies of the higher order fall more rapidly. The first four modal shapes remain almost unchanged, while the fifth modal shape changes significantly with the trend towards bending-torsion coupling. It can be inferred that modal shapes of higher order are more sensitive to the aerodynamic heating effects. Therefore, aerodynamic heating has great impact on the natural vibration characteristics of the hypersonic vehicle
(2)
Comparisons of the first five mode shapes and modal frequencies for the three arbitrary test points illustrate the accuracy of the proposed thermal modal reconstruction method on the considered parameter space. The flutter boundaries of the heated wing were predicted by using the developed thermal modal reconstruction method. Simulation results demonstrate that aerodynamic heating has an important impact on the flutter characteristics of the hypersonic vehicle. The developed thermal modal reconstruction method is very efficient since the natural vibration characteristics of the heated structure at an arbitrary selected operating point in the parameter space can be easily obtained without resorting to the high-fidelity thermal and structural simulations. Therefore, it is very suitable for aeroelastic analysis of the hypersonic vehicle undergoing large parameter variations



---

*Source: 8384639-2019-03-13.xml*