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.
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 Section 2, 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 where and are the mass matrix and stiffness matrix, respectively. represents the dynamic values of three-directional displacement, and is the force vector. The size of is three times the number of element nodes.
According to the expansion theorem, the displacement can be written as
Hence, Eq. (1) can be written in generalized formulation as where vector represents the modal or generalized coordinates and is the generalized load vector. A FEM-based modal analysis leads to the truncated modal matrix with considered eigenmodes, as well as generalized mass, stiffness matrices, and generalized load vector given by
The modal matrix is the mass orthonormal, and is equal to unity, and is a diagonal matrix in which main diagonal elements are natural frequencies . The th generalized force vector element in the time domain can be written as where the pressure coefficient obtained using the third-order piston theory is integrated over the body surface and weighted with the th eigenmode. The third-order piston theory [10] can be expressed as where where is the free stream Mach number, is the speed of sound, is the specific heat ratio of air, and 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 and natural frequencies . 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 is in equilibrium with the radiation heat flux and the conductive heat flux : where is the adiabatic wall temperature and is the wall temperature. After ignoring the heat conduction flux into the structure, the radiation equilibrium temperature can be determined by solving a fourth-order algebraic equation: where is the Boltzmann constant, equal to , and is the emissivity of the surface. The adiabatic-wall temperature and the heat transfer coefficient 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 Section 3. 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. Let denote a specific configuration of the set of parameters, and is the number of sampling points selected on the considered parameter space. Thus, the corresponding modal matrix is calculated at flight condition. In this paper, multivariate interpolation in a tangent space to Grassmann manifold [20] is applied to generate a modal matrix at a new sampling point 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 all -dimensional subspaces of is called the Grassmann manifold which is denoted as . So each -dimensional subspace of can be considered as a point of the manifold . Here, , the dimension of the truncated modal matrix , is fixed independently from the operating parameters. The subspace spanned by belongs to the Grassmann manifold . Because the manifold is a Riemannian manifold [20], there is a tangent space that is a vector space at each point . Therefore, the representative matrix of a point in can be interpolated directly. A new subspace associated with the new operating parameters can be obtained by interpolating the known subspaces in four steps as follows. (1)Choose a reference point , which is treated as the origin point for the interpolation(2)Each point sufficiently close to is mapped to a matrix representing a point of using the logarithm map and based on thin singular value decomposition (SVD). This can be written aswhere is a diagonal matrix containing the singular values, and and are orthogonal matrices. (3)Compute the matrix for the target operating point by interpolating the corresponding entries of the matrix . For this purpose, use a multivariate interpolation scheme [21](4)Compute the exponential of to obtain the desired modal matrix . This can be written as
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 points , , in order to produce modal matrices and natural frequencies . In a multidimensional parametric space, a nonlinear functional relationship between natural frequencies and design variables is expressed by where is the natural frequency of the th mode shape and 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 of sampling points in the dimensional parameter space. The corresponding responses are denoted by . The Kriging approximation is expressed at an arbitrary point as where are polynomial functions, are the corresponding parameters, and 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 Kriging surfaces based on the sampled data set , . Therefore, the prediction of at an arbitrary point on parametric space can be shown to be where are the outputs from sampled values, is the vector of at the sampling points, is the correlation function matrix, is the correlation between unknown points and the known sampling points, and is the least squares predictor given by
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 Figure 1. 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.
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.
(a) Mode 1
(b) Mode 2
(c) Mode 3
(d) Mode 4
(e) Mode 5
The 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 Figure 4, 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 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.
6.2. Thermal Modal Analysis Using High-Fidelity Analysis Tools
In this section, two flight parameters (free stream Mach number and altitude ) 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) , 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 and . 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.
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 and . 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 (-value) 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.
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 , 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 and test point is the shortest, the sampling point 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 .
(a) Case 1: mode 1,
(b) Case 2: mode 1,
(c) Case 1: mode 2,
(d) Case 2: mode 2,
(e) Case 1: mode 3,
(f) Case 2: mode 3,
(g) Case 1: mode 4,
(h) Case 2: mode 4,
(i) Case 1: mode 5,
(j) Case 2: mode 5,
Here, the accuracy of the restructured thermal mode is validated at three arbitrary test points, , , and . 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 and , the MAC value is given as
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 Table 1 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 Figure 9. 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%.
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. Figures 10 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 , but the response of the unheated wing is convergent. In Figure 11, the response of the unheated wing is stable at , 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.
(a)
(b)
(a)
(b)
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
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no competing interests regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China under Grant 11732013 and 11472128.