- Email: [email protected]
Thermal flutter prediction at trajectory points of a hypersonic vehicle based on aerothermal synchronization algorithm
Aerospace Science and Technology 94 (2019) 105381
Contents lists available at ScienceDirect
Aerospace Science and Technology www.elsevier.com/locate/aescte
Thermal flutter prediction at trajectory points of a hypersonic vehicle based on aerothermal synchronization algorithm Tongqing Guo a , Ennan Shen b , Zhiliang Lu a,∗ , Di Zhou a , Jiangpeng Wu b a
Key Laboratory of Unsteady Aerodynamics and Flow Control, Ministry of Industry and Information Technology, Nanjing University of Aeronautics and Astronautics, Yudao Street 29, Nanjing, Jiangsu 210016, China b Shenyang Aircraft Design and Research Institute, Shenyang, Liaoning 110035, China
a r t i c l e
i n f o
Article history: Received 6 June 2019 Received in revised form 30 July 2019 Accepted 3 September 2019 Available online 5 September 2019 Keywords: Hypersonic vehicle Trajectory Thermal flutter Aerothermal Synchronization algorithm CFD/CSD
a b s t r a c t Due to orders of magnitude differences in time scale between structural heat transfer and aeroelastic responses, one-way aerothermal-aeroelastic coupling is adopted to develop a thermal flutter prediction method for a hypersonic vehicle operating along a desired trajectory. In view of the strong dependency of the heat transfer process on the unsteady hypersonic trajectory, an aerothermal synchronization algorithm is established in a non-inertial frame of reference by formulating the governing equations of fluid flow and heat transfer into a unified form. Then the heated free-vibration frequencies and mode shapes are calculated at each trajectory point by using a finite-element analysis. Consequently, the flutter computations are performed on the transiently heated structure at each trajectory point by utilizing a coupled computational fluid dynamics (CFD)/computational structural dynamics (CSD) method. Because of the mass dissimilarity caused by directly increasing the dynamic pressure of a compressible flow, the technique of variable stiffness is introduced to evaluate the flutter dynamic pressure at the point of mass similarity and the stiffness margin of flutter. The present method is applied to the thermal flutter computations of a hypersonic all-movable rudder operating along a given trajectory. The computed temperature differences between the synchronization and conventional partitioned methods, and the significant effects of aerodynamic heating on the structural modes and the flutter characteristics are analyzed in detail. © 2019 Elsevier Masson SAS. All rights reserved.
1. Introduction Hypersonic vehicles operate over a broad range of Mach numbers and must fly within the atmosphere for sustained periods of time to meet the needs of an airbreathing propulsion system [1–3]. The body, surface panels, and aerodynamic control surfaces of a hypersonic vehicle are generally flexible due to minimum-weight restrictions. The combined extreme aerodynamic heating and loading acting on the system produce complex fluid-thermal-structural interactions [4]. The testing of aerothermoelastically scaled windtunnel models in hypersonic flow is not feasible [3,5,6]. Thus, computational aerothermoelasticity is essential and has been an active area of research. Fig. 1 illustrates that the aerothermoelastic problem can be conceptually divided into an aerothermal problem and an aeroelastic problem [2,7].
*
Corresponding author. E-mail address: [email protected] (Z. Lu).
https://doi.org/10.1016/j.ast.2019.105381 1270-9638/© 2019 Elsevier Masson SAS. All rights reserved.
Due to high computational efficiency and ease of implementation, approximate theories for the unsteady aerodynamics are extensively employed in the aerothermoelastic analysis, such as piston theory [8,9], Van Dyke second-order theory [10], and unsteady Newtonian impact theory [11]. These theories are based on the inviscid and quasi-steady flow assumptions and neglect the real gas effects. Because of significant viscosity interactions and intense aerodynamic heating, the exact solution to the hypersonic aerothermodynamic problem can only be obtained by solving the unsteady Navier-Stokes (N-S) equations [2]. With rapid increases in computational power, high-fidelity CFD-based aeroelastic analysis has emerged as a viable approach for aerothermoelasticity. Meanwhile, both the heat transfer and structural dynamics are usually analyzed by the finite-element method (FEM). Recent research has also focused on developing computationally efficient reduced-order models (ROMs) from CFD and FEM solvers [12,13]. As shown in Fig. 1, prediction of aerothermoelasticity requires simultaneous solutions of the fluid, structure, and heat transfer problems. In addition to numerical issues, a realistic tightly coupled analysis is computationally expensive owing to the disparate time scales. In view of orders of magnitude differences in time
2
T. Guo et al. / Aerospace Science and Technology 94 (2019) 105381
chronization algorithm will be developed in a non-inertial frame in order to accurately solve the aerothermal problem along a hypersonic trajectory. A CFD/CSD-based time-domain method with variable stiffness will be then extended to the flutter computations of the heated vehicle at each trajectory point. Finally, the proposed thermal flutter prediction method will be applied to a hypersonic all-movable rudder along a given flight trajectory. 2. An aerothermal synchronization algorithm along a hypersonic trajectory
Fig. 1. Basic structure of the aerothermoelastic problem [2].
In an inertial frame of reference, we have developed a synchronization method on a static mesh [18] for the fluid-thermal interactions. In Ref. [18], this method has been well verified through its applications to the aerodynamic heating on a cylindrical leading edge at a Mach number 6.47 and on a biconic model with increasing freestream Mach number from 0.37 to 9.86 within 0.35 s. Due to the significant variations of flight altitude, Mach number and posture along a hypersonic trajectory, the above synchronization method is to be extended to a non-inertial frame with the dynamic grid methodology in this section. 2.1. Unified governing equations
Fig. 2. Cycle of the serial staggered approach for partitioned aerothermal coupling.
scale between the heat transfer and aeroelastic responses, historically, the aerothermoelastic problem has been decoupled into an aerothermal problem and an aeroelastic problem on the heated structure [2,6,7]. In the case of small static aeroelastic coupling, this one-way coupling approach neglecting path 2 in Fig. 1 has been universally applied to the CFD-based aerothermoelastic analysis [3,6,14,15]. Otherwise, two-way coupling should be considered [2,16]. For a hypersonic vehicle flying within the atmosphere, the flight altitude, Mach number and posture vary with time significantly along the trajectory. Therefore, the structural heat transfer process strongly depends on the unsteady hypersonic trajectory. As reviewed in Ref. [2], the majority of aerothermoelastic analyses were conducted at a fixed flight state. In the limited number of studies considering the flight trajectory [3,17], the unsteady trajectory was decomposed into a series of quasi-steady flight states, and then a one-way coupling analysis was implemented at each trajectory point. In one-way aerothermal-aeroelastic coupling, an accurate aerothermal solution is vital to the subsequent aeroelastic prediction. There are two basic categories of aerothermal coupling [2]: partitioned and synchronized. As shown in Fig. 2, in a partitioned approach, the fluid and the heat transfer are solved using separate solvers in a serial staggered sequence. Hence, two issues arise: the wall temperature T w and heat flux q w transfer at the fluid-structural interface with mismatched grids and the time lag between two solvers. Note that, when a quasi-steady flow assumption is used, the path 2 in Fig. 2 denotes a steady computation under the operating conditions at t = t 1 . To eliminate these drawbacks, a synchronization method combines the governing equations of two fields into a consistent scheme and marches them forward in time simultaneously. However, existing synchronization methods are all established in an inertial frame of reference, which are not suitable for a hypersonic vehicle operating along a trajectory with a wide range of Mach numbers. Because the transient heating of a hypersonic vehicle is inherently linked to its unsteady flight trajectory, in this work, a syn-
The N-S equations are applied to simulate the unsteady fluid flow along a hypersonic trajectory. The variations of flight altitude and attack angle are modeled by changing the freestream parameters and using the dynamic grid methodology, respectively. For efficiency and convenience, a translating frame of reference is introduced to consider the variation of Mach number. In the translating frame moving at the transient flight Mach number, a mathematical derivation leads to the following integral form of the unsteady N-S equations on a dynamic grid
∂ ∂t
ˆ
˛
W d +
( F c − F v )dS = 0
(1)
∂
where denotes a dynamic control volume with a surface element dS, W , F c and F v are the vectors of the conservative variables, the convective fluxes and the viscous fluxes defined as
⎤
⎡
ρ W = ⎣ ρv ⎦, ρE ⎡
0
⎤
⎡
ρ (V r − V g ) F c = ⎣ ρ ( V r − V g ) v + pn ⎦ , ρ E (V r − V g ) + p V ⎤
(2)
⎦ F v = ⎣ τ¯¯ · n ¯ k(∇ T · n) + (τ¯ · v ) · n In the above expressions, ρ , p, T , E, v denote the density, pressure, temperature, absolute total energy per unit mass and the absolute velocity vector, respectively. k is the thermal conductivity coefficient and τ¯¯ is the viscous stress tensor. n refers to the outward facing unit normal vector of dS. V and V r respectively stand for the absolute and relative contravariant velocities on dS, i.e., V = v · n and V r = v r · n with v r being the relative velocity vector. V g represents the contravariant velocity of the face of the control volume. Hence, V g = v g · n with v g being the grid velocity vector. In a translating frame of reference, the absolute velocity v is used as an independent variable in Eq. (1). In this manner, the introduction of other relative variables is avoided and the computational formula for pressure has the same form as that in an inertial frame. Also, the related numerical scheme construction and boundary condition specification can be easily implemented. For an undeformed structure, the heat transfer equation without internal heat source can be obtained by only remaining the
T. Guo et al. / Aerospace Science and Technology 94 (2019) 105381
terms of time derivative and temperature gradient in the fluid flow energy equation, that is
∂ ∂t
ˆ
˛
ρs E s d +
−ks (∇ T s · n)dS = 0
(3)
∂
where the subscript s denotes the solid domain and E s represents the internal energy per unit mass. E s = c s T s with c s being the structural specific heat coefficient. Now, the governing equations for fluid flow and structural heat transfer are formulated into a unified integral form in a translating frame of reference. Additionally, the following Geometric Conservation Law (GCL) has to be fulfilled in order to avoid the errors induced by the deformations of the control volumes [19]
∂ ∂t
¨
˛ d −
v g · ndS = 0
(4)
∂
2.2. Unified numerical approaches The cell-centered structured finite volume method is applied to the spatial discretization of Eqs. (1) and (3). The convective and viscous fluxes are evaluated by using the improved advection upstream splitting method (AUSM+ ) [20] and the central scheme, respectively. The Spalart-Allmaras one-equation turbulence model [21] is adopted in the turbulence modeling. For a steady problem, the implicit lower-upper symmetric Gauss-Seidel (LU-SGS) scheme [22] with the local time step is employed in the temporal discretization. The dual time-stepping approach [19] is used in the unsteady simulation. A parallel algorithm is utilized for higher computational efficiency. 2.3. Boundary conditions
k f ∇ T f · n = −k s ∇ T s · n
3. CFD/CSD-based time-domain method for thermal flutter prediction The structural transient temperature distributions along a trajectory have been determined in the above section. The freevibration frequencies and mode shapes of the transiently heated structure are then calculated at each desired trajectory point in time using a finite-element analysis. Both the material-property degradations and the thermal-stress effects are included in this calculation. Consequently, the thermal flutter at each trajectory point can be predicted by employing an existing CFD/CSD-based time-domain method as well as the corresponding heated modal data. 3.1. CFD and CSD solvers The unsteady CFD solver in section 2 is used to calculate the unsteady aerodynamic forces, and the RBF-TFI dynamic grid method is employed to model the structural deformation. At each trajectory point, the transient solution from the aerothermal simulation supplies the initial flow field and the thermal boundary conditions for the coupled CFD/CSD computation. For a flexible vehicle undergoing small structural deformations, the deformation at any point of the structure can be described by a finite series of modes as follows
r (x, y , z, t ) =
m
h i (x, y , z)qi (t )
(5)
where the subscripts f and s refer to the fluid and structural domains, respectively. 2.4. Dynamic grid generation An efficient hybrid radial basis functions-transfinite interpolation (RBF-TFI) method [23] is applied to the multi-block dynamic grid generation along a flight trajectory. In RBF-TFI, all block vertexes with known deformations are selected as the center points and the grid deformations of all block edges are accurately calculated by RBF. Then, the arc-length-based 2-D and 3-D TFIs are used to fast evaluate the grid deformations on all block faces and inside all blocks, respectively. In a region with complex deformation, the original RBF-TFI in Ref. [23] may fail if the center points of RBF are too few to reflect the real deformation. Instead of only choosing block vertexes as the center points, an effective improvement is to increase the number of center points properly. In this work, the initial multi-block grid is coarsened in a specified ratio to generate a new coarse grid, and all grid points with known deformations are then selected as the center points.
(6)
i =1
where r is the structural deformation vector, h i is the vector of the i-th mode shape, q i is the generalized coordinate of mode i, and m is the number of modes. Then, the following structural equations of motion can be derived from Lagrange’s equations
Mqtt + Gqt + K q = A
At the fluid-structural interface with matching grids, the boundary conditions of noslip and zero normal pressure gradient are imposed. The interface temperature is evaluated from the following temperature continuity and heat flux equilibrium conditions
Ts = T f
3
(7)
where q = [q1 , q2 , · · · , qm ] , M , G and K are respectively the generalized mass, damping and stiffness matrices, and A is the generalized aerodynamic force vector. M , G and K can be obtained from the structural finite-element analysis, and A is calculated from the following surface integral over the aerodynamic surfaces T
‹
Ai =
− p (x, y , z, t )n + τ¯¯ (x, y , z, t ) · n hi (x, y , z)dS
(8)
where p (x, y , z, t ) and τ¯¯ (x, y , z, t ) are the transient pressure and viscous stress tensor, respectively, which are determined from the unsteady CFD results. A second-order hybrid linear multi-step scheme using a predictor-corrector procedure [24,25] is employed to solve Eq. (7), in which the flow field is solved only once at each physical time step by introducing a polynomial extrapolation of the generalized aerodynamic force to the corrector step. With respect to the data transfer at the fluid-structure interface, at first, the mode shapes are interpolated from the structural points to the CFD grid on the vehicle surfaces by RBF. They are further interpolated to the whole CFD grid by RBF-TFI. Consequently, the coupled CFD/CSD iterations are performed only on the CFD grid and the dynamic grid is efficiently generated by using Eq. (6) together with the interpolated spatial modal data. 3.2. CFD/CSD-based time-domain method with variable stiffness A coupled CFD/CSD method is applied to the flutter prediction of the heated vehicle at each trajectory point, in which the aerodynamic equations are advanced simultaneously with the
4
T. Guo et al. / Aerospace Science and Technology 94 (2019) 105381
Fig. 3. Solution procedure of the present thermal flutter prediction method.
structural equations in the time domain. After enough CFD/CSD coupling iterations, the flutter characteristics can be assessed by analyzing the histories of generalized coordinates. More details about our CFD/CSD-based time-domain method can be found in Refs. [25–27]. In Refs. [2] and [3], the aerothermoelastic behavior of a wing along the trajectory was computed by holding the Mach number and altitude constant for each trajectory point and only increasing the dynamic pressure until flutter, and a flutter margin was then determined from the ratio of the virtual-flutter dynamic pressure to the freestream dynamic pressure. As we pointed out in Ref. [27], increasing the dynamic pressure would result in mass dissimilarity, since it means to increase the fluid density in the case of compressible flow at specified flight altitude and Mach number. In this work, the technique of variable stiffness presented in Ref. [27] is employed to calculate the flutter dynamic pressure at the point of mass similarity, in which the structural stiffness is gradually decreased until flutter is achieved. Theoretical analysis indicates that the flutter dynamic pressure Q F and the flutter frequency f F at the point of mass similarity are finally evaluated as
(2) Simulate the unsteady fluid-thermal fields along the trajectory by using the aerothermal synchronization algorithm. (3) At each trajectory point, interpolate the transient structural temperature distribution on the finite-volume model for aerothermal analysis to the finite-element model for structural dynamics analysis. (4) At each trajectory point, calculate the free-vibration frequencies and mode shapes of the heated structure by using the finite-element analysis including material-property degradations and thermal-stress effects. (5) At each trajectory point, interpolate the mode shapes from the structural points to the CFD grid on the aerodynamic surfaces by RBF, and then interpolate them to the whole CFD grid by RBF-TFI. (6) At each trajectory point, utilize the transient fluid-thermal fields as the initial flow flied and the thermal boundary conditions, and conduct the flutter computation of the heated vehicle by employing the CFD/CSD-based time-domain method with variable stiffness. Resultantly, the thermal flutter characteristics are obtained for all trajectory points. 4. Numerical examples and discussions
Q Fa = Q Fm / N F
(9)
f Fa = f Fm / N F
(10)
where the superscripts a and m refer to the actual vehicle and the computational model with decreased stiffness, respectively, and N F denotes the critical stiffness coefficient of flutter relative to the actual value. Hence, the stiffness margin of flutter equals (1 − N F ). 3.3. Solution procedure Up to now, the thermal flutter computational method for a hypersonic vehicle along its flight trajectory has been established by combining the aerothermal synchronization algorithm with the separate CFD/CSD-based time-domain flutter calculations. Its main solution procedure illustrated in Fig. 3 can be outlined below. (1) Initialize the fluid flow and structural temperature fields according to the operation condition at the first point on the trajectory, and compute the corresponding steady-state flow solution under an isothermal boundary condition.
As depicted in Fig. 4, in this section, the present method is applied to the thermal flutter analysis of a hypersonic all-movable rudder operating along a trajectory. Its leading-edge radius is 2.5 × 10−3 m. The stiffness of the torsional spring attached to the hinge joint at the half-root-chord location is 1 × 107 (N m)/deg. The structural density, specific heat coefficient and thermal conductivity coefficient are 7900 kg/m3 , 500 J/(kg K) and 16.3 W/(m K), respectively. The initial structural temperature equals 293.15 K. Fig. 5 displays the multi-block structured grid of the rudder. In the fluid and structural domains with matching interface grid, the cell numbers are about 1.5 and 0.2 million, respectively. The non-dimensional wall coordinate of the first off-wall grid node approximately equals one. The hypersonic trajectory of the rudder is represented by five discretized points on the trajectory. The operating conditions at each point are listed in Table 1. The operating parameters at any time are linearly interpolated from two neighboring trajectory points.
T. Guo et al. / Aerospace Science and Technology 94 (2019) 105381
5
Table 1 The operating conditions at each point on the trajectory. Trajectory point No.
Time (s)
Mach number
Attack angle (deg.)
Altitude (km)
Pressure (Pa)
Temperature (K)
Density (kg/m3 )
1 2 3 4 5
0.0 4.0 8.0 10.0 12.0
1.5 2.5 4.0 5.0 6.0
12.0 8.0 4.0 2.0 2.0
4 5.676 7.36 7.925 8
61660.4 49345.5 39069.4 36038.9 35651.6
262.166 251.285 240.366 236.701 236.215
0.819346 0.684097 0.566242 0.530406 0.525786
Fig. 4. Planform and cross-sectional views of the rudder (unit: m). Fig. 5. Multi-block structured grid of the rudder.
4.1. Aerothermal computations The physical time step for aerothermal synchronization computation is set as 5 × 10−4 s. For comparison, both the synchronization algorithm and the conventional partitioned method in Fig. 2 using a quasi-steady flow approximation [17] are employed to compute the aerothermal problem of current rudder along its trajectory. Fig. 6 compares the computed time histories of heat flux and temperature at the leading-edge of the root section between two methods. For the synchronization method, the time histories well reflect the discontinuity of operating conditions at each trajectory point in Table 1. For the partitioned method, a quasi-steady flow computation without structural heat transfer leads to a higher heat flux. The maximum temperature difference is up to 8% at trajectory point 4 (t = 10 s). 4.2. Heated and unheated structural modes The transient structural temperature distributions from the previous aerothermal synchronization computation are interpolated
Table 2 Temperature-dependent material properties. Temperature (K)
Modulus of elasticity (Pa)
Thermal expansion coefficient (1/K)
293.15 500 800 1500
7.84 × 1011 2 × 1010 1 × 1010 5 × 109
1.66 × 10−5 1.72 × 10−5 1.81 × 10−5 2 × 10−5
from the finite-volume model to the finite-element model. Fig. 7 displays the interpolated temperature distributions on the finiteelement model at trajectory point 4 and point 5. Table 2 lists the temperature-dependent material properties. Then the free-vibration frequencies and mode shapes of the heated structure are calculated by using the finite-element analysis including material-property degradations and thermal-stress effects. For comparison, Table 3 provides the first ten modal fre-
Fig. 6. Time histories of heat flux (left) and temperature (right) at the leading-edge of the root section.
6
T. Guo et al. / Aerospace Science and Technology 94 (2019) 105381
Fig. 7. Transient temperature distributions on the finite-element model at two trajectory points: 4 (left) and 5 (right). (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)
Fig. 8. Comparison of the first four mode shapes between the cold (left) and hot (right) models at trajectory point 5.
quencies of the cold and hot models. Here, the cold and hot models refer to excluding and including aerodynamic heating effects, respectively. At trajectory points 4 and 5, aerodynamic heating reduces the first ten natural frequencies by nearly 20% and 50% except the first torsional mode. In this work, the heating effects on the hinge joint are neglected. Fig. 8 further compares the first four mode shapes between two models at trajectory point 5. As can be seen from Fig. 8, the first two modes are the first torsional mode and the first bending mode, respectively; there is an apparent difference in each mode shape between two models, but the mode order is unchanged.
4.3. Flutter computations of the cold and hot models To analyze the aerodynamic heating effects on the flutter characteristics, flutter computations are performed on both the cold and hot models. Fig. 9 gives the time histories of generalized coordinates of the hot model at trajectory points 4 and 5. As can be seen from Fig. 9a to Fig. 9c, at trajectory point 4, the time histories of generalized coordinates vary from convergence to divergence with the decrease of structural stiffness coefficient N, and the critical stiffness coefficient of flutter is about 0.18. Fig. 10 shows the time histories of generalized coordinates of the cold model at the
T. Guo et al. / Aerospace Science and Technology 94 (2019) 105381
7
Table 3 Comparison of the first ten modal frequencies between the cold and hot models. Mode order
Cold model
Hot model Trajectory point 4
1 (first torsion) 2 (first bending) 3 4 5 6 7 8 9 10
Trajectory point 5
Frequency (Hz)
Frequency (Hz)
% decrease
Frequency (Hz)
% decrease
14.527 30.954 76.406 101.44 143.7 158.88 164.81 204.57 220.28 245.96
13.883 23.916 56.694 78.24 111.97 118.28 127.95 164.0 177.15 186.67
4.43% 22.7% 25.8% 22.9% 22.1% 25.5% 22.3% 19.8% 19.6% 24.1%
12.617 16.285 32.461 39.997 54.421 57.776 61.372 82.135 84.576 92.612
13.1% 47.4% 57.5% 60.5% 62.1% 63.6% 62.5% 59.8% 61.6% 62.3%
Fig. 9. Time histories of the generalized coordinates of the hot model.
Table 4 Comparison of the flutter characteristics between the cold and hot models. Case
Critical stiffness coefficient NF Stiffness margin of flutter Flutter dynamic pressure (Pa) Flutter frequency (Hz)
Trajectory point 4
Trajectory point 5
Cold model
Hot model
% decrease
Cold model
Hot model
% decrease
0.066 93.4% 9.56 × 106 24.42
0.18 82% 3.51 × 106 18.01
– 12.2% 63.3% 26.2%
0.082 91.8% 1.1 × 107 21.60
>1 <0
– – – 33.3%
– 14.41
8
T. Guo et al. / Aerospace Science and Technology 94 (2019) 105381
Fig. 10. Time histories of the generalized coordinates of the cold model.
critical flutter point. In all cases including the cold and hot models, a coupled bending-torsion flutter is obtained. Table 4 outlines the flutter characteristics of the cold and hot models. For the hot model, the large reduction in the natural frequencies results in a significant decrease in the flutter margin. At trajectory point 4, the reductions in the stiffness margin, the flutter dynamic pressure and the flutter frequency are up to 12.2%, 63.3% and 26.2%, respectively. At trajectory point 5, the stiffness margin of the cold model is 91.8%, but the heated rudder has fluttered at the original stiffness (N = 1), which can be seen from Fig. 9d.
This work was supported by the National Natural Science Foundation of China (Grant No. 11872212) and a project funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions.
5. Conclusions
References
For the hypersonic aerothermal problem along a flight trajectory, a synchronization algorithm is developed in a non-inertial frame of reference by formulating the governing equations of fluid flow and structural heat transfer into a unified form. The synchronization algorithm is then combined with the CFD/CSD-based time-domain method with variable stiffness. Consequently, the thermal flutter prediction method is achieved. The present method has been applied to the thermal flutter computations of a hypersonic all-movable rudder along a given flight trajectory. From our development process and numerical analyses, we conclude as follows.
[1] J.J. Bertin, Hypersonic Aerothermodynamics, AIAA, Reston, VA, 1994. [2] J.J. McNamara, P.P. Friedmann, Aeroelastic and aerothermoelastic analysis in hypersonic flow: past, present, and future, AIAA J. 49 (2011) 1089–1122. [3] J.J. McNamara, P.P. Friedmann, K.G. Powell, B.J. Thuruthimattam, R.E. Bartels, Aeroelastic and aerothermoelastic behavior in hypersonic flow, AIAA J. 46 (2008) 2591–2610. [4] J.J. Bertin, R.M. Cummings, Fifty years of hypersonics: where we’ve been, where we’re going, Prog. Aerosp. Sci. 39 (2003) 511–536. [5] J. Dugundji, J. Calligeros, Similarity laws for aerothermoelastic testing, J. Aerosp. Sci. 29 (1962) 935–950. [6] N. Lamorte, P.P. Friedmann, Hypersonic aeroelastic and aerothermoelastic studies using computational fluid dynamics, AIAA J. 52 (2014) 2062–2078. [7] R.L. Bisplinghoff, Some structural and aeroelastic considerations of high-speed flight, J. Aeronaut. Sci. 23 (1956) 289–329. [8] H. Ashley, G. Zartarian, Piston theory: a new aerodynamic tool for the aeroelastician, J. Aeronaut. Sci. 23 (1956) 1109–1118. [9] W.W. Zhang, Z.Y. Ye, C.A. Zhang, F. Liu, Supersonic flutter analysis based on a local piston theory, AIAA J. 47 (2009) 2321–2328. [10] M. Van Dyke, A Study of Second-Order Supersonic Flow Theory, NACATR 1081, 1951. [11] D. Zhang, S. Tang, Q.J. Zhu, R.G. Wang, Analysis of dynamic characteristics of the rigid body/elastic body coupling of air-breathing hypersonic vehicles, Aerosp. Sci. Technol. 48 (2016) 328–341. [12] Z.Q. Chen, Y.H. Zhao, R. Huang, Parametric reduced-order modeling of unsteady aerodynamics for hypersonic vehicles, Aerosp. Sci. Technol. 87 (2019) 1–14. [13] C.Q. Gao, W.W. Zhang, Z.Y. Ye, A new viewpoint on the mechanism of transonic single-degree-of-freedom flutter, Aerosp. Sci. Technol. 52 (2016) 144–156. [14] K. Ye, Z.Y. Ye, C.N. Li, J. Wu, Effects of the aerothermoelastic deformation on the performance of the three-dimensional hypersonic inlet, Aerosp. Sci. Technol. 84 (2019) 747–762. [15] F. Chen, H. Liu, S.T. Zhang, Time-adaptive loosely coupled analysis on fluidthermal-structural behaviors of hypersonic wing structures under sustained aeroheating, Aerosp. Sci. Technol. 78 (2018) 620–636. [16] C. Yang, G.S. Li, Z.Q. Wan, Aerothermal-aeroelastic two-way coupling method for hypersonic curved panel flutter, Sci. China, Technol. Sci. 55 (2012) 831–840. [17] S. Zhang, F. Chen, H. Liu, Integrated Fluid-Thermal-Structural Analysis for Predicting Aerothermal Environment of Hypersonic Vehicles, AIAA paper 2014-1394, 2014.
(1) The synchronization algorithm can accurately predict the hypersonic aerothermal problem along a flight trajectory, including the unsteady effects of the variations of flight altitude, Mach number and attack angle. For current rudder along its flight trajectory, a maximum temperature difference at the monitored point is up to 8% between the synchronization and conventional partitioned algorithms. (2) The present thermal flutter computational method can quantitatively evaluate the flutter margin at each point on the flight trajectory. For current rudder at trajectory points 4 and 5, aerodynamic heating reduces the first ten natural frequencies by nearly 20% and 50% except the first torsional mode partly due to neglecting the heating effects on the hinge joint. Consequently, the flutter margin decreases significantly: at trajectory point 4, the reductions in the stiffness margin and the flutter dynamic pressure are up to 12.2% and 63.3%, respectively; at trajectory point 5, the stiffness margin of the cold model is 91.8%, but the heated rudder has fluttered at the original stiffness.
Declaration of competing interest There is no competing interest. Acknowledgements
T. Guo et al. / Aerospace Science and Technology 94 (2019) 105381
[18] E.N. Shen, Z.L. Lu, D. Zhou, T.Q. Guo, Development of a synchronization method for fluid-thermal study of hypersonic flow, Trans. Nanjing Univ. Aeronaut. Astronaut. 35 (2018) 973–985. [19] A.L. Gaitonde, A Dual-Time Method for the Solution of the 2d Unsteady NavierStokes Equations on Structured Moving Meshes, AIAA Paper 1995-1877, 1995. [20] M.S. Liou, A sequel to AUSM: AUSM+ , J. Comput. Phys. 129 (1996) 364–382. [21] P.R. Spalart, S.R. Allmaras, A One-Equation Turbulence Model for Aerodynamic Flows, AIAA paper 1992-0439, 1992. [22] H. Rieger, A. Jameson, Solution of Steady 3-d Compressible Euler and NavierStokes Equations by an Implicit LU Scheme, AIAA Paper 88-0619, 1988. [23] L. Ding, Z.L. Lu, T.Q. Guo, An efficient dynamic mesh generation method for complex multi-block structured grid, Adv. Appl. Math. Mech. 6 (2014) 120–134.
9
[24] W.W. Zhang, Y.W. Jiang, Z.Y. Ye, Two better loosely coupled solution algorithms of CFD based aeroelastic simulation, Eng. Appl. Comput. Fluid Mech. 1 (2007) 253–262. [25] T.Q. Guo, Z.L. Lu, D. Tang, T.G. Wang, L. Dong, A CFD/CSD model for aeroelastic calculations of large-scale wind turbines, Sci. China, Technol. Sci. 56 (2013) 205–211. [26] T.Q. Guo, Z.L. Lu, Y.J. Wu, A time-domain method for transonic flutter analysis with multidirectional coupled vibrations, Mod. Phys. Lett. B 23 (2009) 453–456. [27] T.Q. Guo, Z.L. Lu, A CFD/CSD model for transonic flutter, Comput. Mater. Continua 2 (2005) 105–111.
Contents lists available at ScienceDirect
Aerospace Science and Technology www.elsevier.com/locate/aescte
Thermal flutter prediction at trajectory points of a hypersonic vehicle based on aerothermal synchronization algorithm Tongqing Guo a , Ennan Shen b , Zhiliang Lu a,∗ , Di Zhou a , Jiangpeng Wu b a
Key Laboratory of Unsteady Aerodynamics and Flow Control, Ministry of Industry and Information Technology, Nanjing University of Aeronautics and Astronautics, Yudao Street 29, Nanjing, Jiangsu 210016, China b Shenyang Aircraft Design and Research Institute, Shenyang, Liaoning 110035, China
a r t i c l e
i n f o
Article history: Received 6 June 2019 Received in revised form 30 July 2019 Accepted 3 September 2019 Available online 5 September 2019 Keywords: Hypersonic vehicle Trajectory Thermal flutter Aerothermal Synchronization algorithm CFD/CSD
a b s t r a c t Due to orders of magnitude differences in time scale between structural heat transfer and aeroelastic responses, one-way aerothermal-aeroelastic coupling is adopted to develop a thermal flutter prediction method for a hypersonic vehicle operating along a desired trajectory. In view of the strong dependency of the heat transfer process on the unsteady hypersonic trajectory, an aerothermal synchronization algorithm is established in a non-inertial frame of reference by formulating the governing equations of fluid flow and heat transfer into a unified form. Then the heated free-vibration frequencies and mode shapes are calculated at each trajectory point by using a finite-element analysis. Consequently, the flutter computations are performed on the transiently heated structure at each trajectory point by utilizing a coupled computational fluid dynamics (CFD)/computational structural dynamics (CSD) method. Because of the mass dissimilarity caused by directly increasing the dynamic pressure of a compressible flow, the technique of variable stiffness is introduced to evaluate the flutter dynamic pressure at the point of mass similarity and the stiffness margin of flutter. The present method is applied to the thermal flutter computations of a hypersonic all-movable rudder operating along a given trajectory. The computed temperature differences between the synchronization and conventional partitioned methods, and the significant effects of aerodynamic heating on the structural modes and the flutter characteristics are analyzed in detail. © 2019 Elsevier Masson SAS. All rights reserved.
1. Introduction Hypersonic vehicles operate over a broad range of Mach numbers and must fly within the atmosphere for sustained periods of time to meet the needs of an airbreathing propulsion system [1–3]. The body, surface panels, and aerodynamic control surfaces of a hypersonic vehicle are generally flexible due to minimum-weight restrictions. The combined extreme aerodynamic heating and loading acting on the system produce complex fluid-thermal-structural interactions [4]. The testing of aerothermoelastically scaled windtunnel models in hypersonic flow is not feasible [3,5,6]. Thus, computational aerothermoelasticity is essential and has been an active area of research. Fig. 1 illustrates that the aerothermoelastic problem can be conceptually divided into an aerothermal problem and an aeroelastic problem [2,7].
*
Corresponding author. E-mail address: [email protected] (Z. Lu).
https://doi.org/10.1016/j.ast.2019.105381 1270-9638/© 2019 Elsevier Masson SAS. All rights reserved.
Due to high computational efficiency and ease of implementation, approximate theories for the unsteady aerodynamics are extensively employed in the aerothermoelastic analysis, such as piston theory [8,9], Van Dyke second-order theory [10], and unsteady Newtonian impact theory [11]. These theories are based on the inviscid and quasi-steady flow assumptions and neglect the real gas effects. Because of significant viscosity interactions and intense aerodynamic heating, the exact solution to the hypersonic aerothermodynamic problem can only be obtained by solving the unsteady Navier-Stokes (N-S) equations [2]. With rapid increases in computational power, high-fidelity CFD-based aeroelastic analysis has emerged as a viable approach for aerothermoelasticity. Meanwhile, both the heat transfer and structural dynamics are usually analyzed by the finite-element method (FEM). Recent research has also focused on developing computationally efficient reduced-order models (ROMs) from CFD and FEM solvers [12,13]. As shown in Fig. 1, prediction of aerothermoelasticity requires simultaneous solutions of the fluid, structure, and heat transfer problems. In addition to numerical issues, a realistic tightly coupled analysis is computationally expensive owing to the disparate time scales. In view of orders of magnitude differences in time
2
T. Guo et al. / Aerospace Science and Technology 94 (2019) 105381
chronization algorithm will be developed in a non-inertial frame in order to accurately solve the aerothermal problem along a hypersonic trajectory. A CFD/CSD-based time-domain method with variable stiffness will be then extended to the flutter computations of the heated vehicle at each trajectory point. Finally, the proposed thermal flutter prediction method will be applied to a hypersonic all-movable rudder along a given flight trajectory. 2. An aerothermal synchronization algorithm along a hypersonic trajectory
Fig. 1. Basic structure of the aerothermoelastic problem [2].
In an inertial frame of reference, we have developed a synchronization method on a static mesh [18] for the fluid-thermal interactions. In Ref. [18], this method has been well verified through its applications to the aerodynamic heating on a cylindrical leading edge at a Mach number 6.47 and on a biconic model with increasing freestream Mach number from 0.37 to 9.86 within 0.35 s. Due to the significant variations of flight altitude, Mach number and posture along a hypersonic trajectory, the above synchronization method is to be extended to a non-inertial frame with the dynamic grid methodology in this section. 2.1. Unified governing equations
Fig. 2. Cycle of the serial staggered approach for partitioned aerothermal coupling.
scale between the heat transfer and aeroelastic responses, historically, the aerothermoelastic problem has been decoupled into an aerothermal problem and an aeroelastic problem on the heated structure [2,6,7]. In the case of small static aeroelastic coupling, this one-way coupling approach neglecting path 2 in Fig. 1 has been universally applied to the CFD-based aerothermoelastic analysis [3,6,14,15]. Otherwise, two-way coupling should be considered [2,16]. For a hypersonic vehicle flying within the atmosphere, the flight altitude, Mach number and posture vary with time significantly along the trajectory. Therefore, the structural heat transfer process strongly depends on the unsteady hypersonic trajectory. As reviewed in Ref. [2], the majority of aerothermoelastic analyses were conducted at a fixed flight state. In the limited number of studies considering the flight trajectory [3,17], the unsteady trajectory was decomposed into a series of quasi-steady flight states, and then a one-way coupling analysis was implemented at each trajectory point. In one-way aerothermal-aeroelastic coupling, an accurate aerothermal solution is vital to the subsequent aeroelastic prediction. There are two basic categories of aerothermal coupling [2]: partitioned and synchronized. As shown in Fig. 2, in a partitioned approach, the fluid and the heat transfer are solved using separate solvers in a serial staggered sequence. Hence, two issues arise: the wall temperature T w and heat flux q w transfer at the fluid-structural interface with mismatched grids and the time lag between two solvers. Note that, when a quasi-steady flow assumption is used, the path 2 in Fig. 2 denotes a steady computation under the operating conditions at t = t 1 . To eliminate these drawbacks, a synchronization method combines the governing equations of two fields into a consistent scheme and marches them forward in time simultaneously. However, existing synchronization methods are all established in an inertial frame of reference, which are not suitable for a hypersonic vehicle operating along a trajectory with a wide range of Mach numbers. Because the transient heating of a hypersonic vehicle is inherently linked to its unsteady flight trajectory, in this work, a syn-
The N-S equations are applied to simulate the unsteady fluid flow along a hypersonic trajectory. The variations of flight altitude and attack angle are modeled by changing the freestream parameters and using the dynamic grid methodology, respectively. For efficiency and convenience, a translating frame of reference is introduced to consider the variation of Mach number. In the translating frame moving at the transient flight Mach number, a mathematical derivation leads to the following integral form of the unsteady N-S equations on a dynamic grid
∂ ∂t
ˆ
˛
W d +
( F c − F v )dS = 0
(1)
∂
where denotes a dynamic control volume with a surface element dS, W , F c and F v are the vectors of the conservative variables, the convective fluxes and the viscous fluxes defined as
⎤
⎡
ρ W = ⎣ ρv ⎦, ρE ⎡
0
⎤
⎡
ρ (V r − V g ) F c = ⎣ ρ ( V r − V g ) v + pn ⎦ , ρ E (V r − V g ) + p V ⎤
(2)
⎦ F v = ⎣ τ¯¯ · n ¯ k(∇ T · n) + (τ¯ · v ) · n In the above expressions, ρ , p, T , E, v denote the density, pressure, temperature, absolute total energy per unit mass and the absolute velocity vector, respectively. k is the thermal conductivity coefficient and τ¯¯ is the viscous stress tensor. n refers to the outward facing unit normal vector of dS. V and V r respectively stand for the absolute and relative contravariant velocities on dS, i.e., V = v · n and V r = v r · n with v r being the relative velocity vector. V g represents the contravariant velocity of the face of the control volume. Hence, V g = v g · n with v g being the grid velocity vector. In a translating frame of reference, the absolute velocity v is used as an independent variable in Eq. (1). In this manner, the introduction of other relative variables is avoided and the computational formula for pressure has the same form as that in an inertial frame. Also, the related numerical scheme construction and boundary condition specification can be easily implemented. For an undeformed structure, the heat transfer equation without internal heat source can be obtained by only remaining the
T. Guo et al. / Aerospace Science and Technology 94 (2019) 105381
terms of time derivative and temperature gradient in the fluid flow energy equation, that is
∂ ∂t
ˆ
˛
ρs E s d +
−ks (∇ T s · n)dS = 0
(3)
∂
where the subscript s denotes the solid domain and E s represents the internal energy per unit mass. E s = c s T s with c s being the structural specific heat coefficient. Now, the governing equations for fluid flow and structural heat transfer are formulated into a unified integral form in a translating frame of reference. Additionally, the following Geometric Conservation Law (GCL) has to be fulfilled in order to avoid the errors induced by the deformations of the control volumes [19]
∂ ∂t
¨
˛ d −
v g · ndS = 0
(4)
∂
2.2. Unified numerical approaches The cell-centered structured finite volume method is applied to the spatial discretization of Eqs. (1) and (3). The convective and viscous fluxes are evaluated by using the improved advection upstream splitting method (AUSM+ ) [20] and the central scheme, respectively. The Spalart-Allmaras one-equation turbulence model [21] is adopted in the turbulence modeling. For a steady problem, the implicit lower-upper symmetric Gauss-Seidel (LU-SGS) scheme [22] with the local time step is employed in the temporal discretization. The dual time-stepping approach [19] is used in the unsteady simulation. A parallel algorithm is utilized for higher computational efficiency. 2.3. Boundary conditions
k f ∇ T f · n = −k s ∇ T s · n
3. CFD/CSD-based time-domain method for thermal flutter prediction The structural transient temperature distributions along a trajectory have been determined in the above section. The freevibration frequencies and mode shapes of the transiently heated structure are then calculated at each desired trajectory point in time using a finite-element analysis. Both the material-property degradations and the thermal-stress effects are included in this calculation. Consequently, the thermal flutter at each trajectory point can be predicted by employing an existing CFD/CSD-based time-domain method as well as the corresponding heated modal data. 3.1. CFD and CSD solvers The unsteady CFD solver in section 2 is used to calculate the unsteady aerodynamic forces, and the RBF-TFI dynamic grid method is employed to model the structural deformation. At each trajectory point, the transient solution from the aerothermal simulation supplies the initial flow field and the thermal boundary conditions for the coupled CFD/CSD computation. For a flexible vehicle undergoing small structural deformations, the deformation at any point of the structure can be described by a finite series of modes as follows
r (x, y , z, t ) =
m
h i (x, y , z)qi (t )
(5)
where the subscripts f and s refer to the fluid and structural domains, respectively. 2.4. Dynamic grid generation An efficient hybrid radial basis functions-transfinite interpolation (RBF-TFI) method [23] is applied to the multi-block dynamic grid generation along a flight trajectory. In RBF-TFI, all block vertexes with known deformations are selected as the center points and the grid deformations of all block edges are accurately calculated by RBF. Then, the arc-length-based 2-D and 3-D TFIs are used to fast evaluate the grid deformations on all block faces and inside all blocks, respectively. In a region with complex deformation, the original RBF-TFI in Ref. [23] may fail if the center points of RBF are too few to reflect the real deformation. Instead of only choosing block vertexes as the center points, an effective improvement is to increase the number of center points properly. In this work, the initial multi-block grid is coarsened in a specified ratio to generate a new coarse grid, and all grid points with known deformations are then selected as the center points.
(6)
i =1
where r is the structural deformation vector, h i is the vector of the i-th mode shape, q i is the generalized coordinate of mode i, and m is the number of modes. Then, the following structural equations of motion can be derived from Lagrange’s equations
Mqtt + Gqt + K q = A
At the fluid-structural interface with matching grids, the boundary conditions of noslip and zero normal pressure gradient are imposed. The interface temperature is evaluated from the following temperature continuity and heat flux equilibrium conditions
Ts = T f
3
(7)
where q = [q1 , q2 , · · · , qm ] , M , G and K are respectively the generalized mass, damping and stiffness matrices, and A is the generalized aerodynamic force vector. M , G and K can be obtained from the structural finite-element analysis, and A is calculated from the following surface integral over the aerodynamic surfaces T
‹
Ai =
− p (x, y , z, t )n + τ¯¯ (x, y , z, t ) · n hi (x, y , z)dS
(8)
where p (x, y , z, t ) and τ¯¯ (x, y , z, t ) are the transient pressure and viscous stress tensor, respectively, which are determined from the unsteady CFD results. A second-order hybrid linear multi-step scheme using a predictor-corrector procedure [24,25] is employed to solve Eq. (7), in which the flow field is solved only once at each physical time step by introducing a polynomial extrapolation of the generalized aerodynamic force to the corrector step. With respect to the data transfer at the fluid-structure interface, at first, the mode shapes are interpolated from the structural points to the CFD grid on the vehicle surfaces by RBF. They are further interpolated to the whole CFD grid by RBF-TFI. Consequently, the coupled CFD/CSD iterations are performed only on the CFD grid and the dynamic grid is efficiently generated by using Eq. (6) together with the interpolated spatial modal data. 3.2. CFD/CSD-based time-domain method with variable stiffness A coupled CFD/CSD method is applied to the flutter prediction of the heated vehicle at each trajectory point, in which the aerodynamic equations are advanced simultaneously with the
4
T. Guo et al. / Aerospace Science and Technology 94 (2019) 105381
Fig. 3. Solution procedure of the present thermal flutter prediction method.
structural equations in the time domain. After enough CFD/CSD coupling iterations, the flutter characteristics can be assessed by analyzing the histories of generalized coordinates. More details about our CFD/CSD-based time-domain method can be found in Refs. [25–27]. In Refs. [2] and [3], the aerothermoelastic behavior of a wing along the trajectory was computed by holding the Mach number and altitude constant for each trajectory point and only increasing the dynamic pressure until flutter, and a flutter margin was then determined from the ratio of the virtual-flutter dynamic pressure to the freestream dynamic pressure. As we pointed out in Ref. [27], increasing the dynamic pressure would result in mass dissimilarity, since it means to increase the fluid density in the case of compressible flow at specified flight altitude and Mach number. In this work, the technique of variable stiffness presented in Ref. [27] is employed to calculate the flutter dynamic pressure at the point of mass similarity, in which the structural stiffness is gradually decreased until flutter is achieved. Theoretical analysis indicates that the flutter dynamic pressure Q F and the flutter frequency f F at the point of mass similarity are finally evaluated as
(2) Simulate the unsteady fluid-thermal fields along the trajectory by using the aerothermal synchronization algorithm. (3) At each trajectory point, interpolate the transient structural temperature distribution on the finite-volume model for aerothermal analysis to the finite-element model for structural dynamics analysis. (4) At each trajectory point, calculate the free-vibration frequencies and mode shapes of the heated structure by using the finite-element analysis including material-property degradations and thermal-stress effects. (5) At each trajectory point, interpolate the mode shapes from the structural points to the CFD grid on the aerodynamic surfaces by RBF, and then interpolate them to the whole CFD grid by RBF-TFI. (6) At each trajectory point, utilize the transient fluid-thermal fields as the initial flow flied and the thermal boundary conditions, and conduct the flutter computation of the heated vehicle by employing the CFD/CSD-based time-domain method with variable stiffness. Resultantly, the thermal flutter characteristics are obtained for all trajectory points. 4. Numerical examples and discussions
Q Fa = Q Fm / N F
(9)
f Fa = f Fm / N F
(10)
where the superscripts a and m refer to the actual vehicle and the computational model with decreased stiffness, respectively, and N F denotes the critical stiffness coefficient of flutter relative to the actual value. Hence, the stiffness margin of flutter equals (1 − N F ). 3.3. Solution procedure Up to now, the thermal flutter computational method for a hypersonic vehicle along its flight trajectory has been established by combining the aerothermal synchronization algorithm with the separate CFD/CSD-based time-domain flutter calculations. Its main solution procedure illustrated in Fig. 3 can be outlined below. (1) Initialize the fluid flow and structural temperature fields according to the operation condition at the first point on the trajectory, and compute the corresponding steady-state flow solution under an isothermal boundary condition.
As depicted in Fig. 4, in this section, the present method is applied to the thermal flutter analysis of a hypersonic all-movable rudder operating along a trajectory. Its leading-edge radius is 2.5 × 10−3 m. The stiffness of the torsional spring attached to the hinge joint at the half-root-chord location is 1 × 107 (N m)/deg. The structural density, specific heat coefficient and thermal conductivity coefficient are 7900 kg/m3 , 500 J/(kg K) and 16.3 W/(m K), respectively. The initial structural temperature equals 293.15 K. Fig. 5 displays the multi-block structured grid of the rudder. In the fluid and structural domains with matching interface grid, the cell numbers are about 1.5 and 0.2 million, respectively. The non-dimensional wall coordinate of the first off-wall grid node approximately equals one. The hypersonic trajectory of the rudder is represented by five discretized points on the trajectory. The operating conditions at each point are listed in Table 1. The operating parameters at any time are linearly interpolated from two neighboring trajectory points.
T. Guo et al. / Aerospace Science and Technology 94 (2019) 105381
5
Table 1 The operating conditions at each point on the trajectory. Trajectory point No.
Time (s)
Mach number
Attack angle (deg.)
Altitude (km)
Pressure (Pa)
Temperature (K)
Density (kg/m3 )
1 2 3 4 5
0.0 4.0 8.0 10.0 12.0
1.5 2.5 4.0 5.0 6.0
12.0 8.0 4.0 2.0 2.0
4 5.676 7.36 7.925 8
61660.4 49345.5 39069.4 36038.9 35651.6
262.166 251.285 240.366 236.701 236.215
0.819346 0.684097 0.566242 0.530406 0.525786
Fig. 4. Planform and cross-sectional views of the rudder (unit: m). Fig. 5. Multi-block structured grid of the rudder.
4.1. Aerothermal computations The physical time step for aerothermal synchronization computation is set as 5 × 10−4 s. For comparison, both the synchronization algorithm and the conventional partitioned method in Fig. 2 using a quasi-steady flow approximation [17] are employed to compute the aerothermal problem of current rudder along its trajectory. Fig. 6 compares the computed time histories of heat flux and temperature at the leading-edge of the root section between two methods. For the synchronization method, the time histories well reflect the discontinuity of operating conditions at each trajectory point in Table 1. For the partitioned method, a quasi-steady flow computation without structural heat transfer leads to a higher heat flux. The maximum temperature difference is up to 8% at trajectory point 4 (t = 10 s). 4.2. Heated and unheated structural modes The transient structural temperature distributions from the previous aerothermal synchronization computation are interpolated
Table 2 Temperature-dependent material properties. Temperature (K)
Modulus of elasticity (Pa)
Thermal expansion coefficient (1/K)
293.15 500 800 1500
7.84 × 1011 2 × 1010 1 × 1010 5 × 109
1.66 × 10−5 1.72 × 10−5 1.81 × 10−5 2 × 10−5
from the finite-volume model to the finite-element model. Fig. 7 displays the interpolated temperature distributions on the finiteelement model at trajectory point 4 and point 5. Table 2 lists the temperature-dependent material properties. Then the free-vibration frequencies and mode shapes of the heated structure are calculated by using the finite-element analysis including material-property degradations and thermal-stress effects. For comparison, Table 3 provides the first ten modal fre-
Fig. 6. Time histories of heat flux (left) and temperature (right) at the leading-edge of the root section.
6
T. Guo et al. / Aerospace Science and Technology 94 (2019) 105381
Fig. 7. Transient temperature distributions on the finite-element model at two trajectory points: 4 (left) and 5 (right). (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)
Fig. 8. Comparison of the first four mode shapes between the cold (left) and hot (right) models at trajectory point 5.
quencies of the cold and hot models. Here, the cold and hot models refer to excluding and including aerodynamic heating effects, respectively. At trajectory points 4 and 5, aerodynamic heating reduces the first ten natural frequencies by nearly 20% and 50% except the first torsional mode. In this work, the heating effects on the hinge joint are neglected. Fig. 8 further compares the first four mode shapes between two models at trajectory point 5. As can be seen from Fig. 8, the first two modes are the first torsional mode and the first bending mode, respectively; there is an apparent difference in each mode shape between two models, but the mode order is unchanged.
4.3. Flutter computations of the cold and hot models To analyze the aerodynamic heating effects on the flutter characteristics, flutter computations are performed on both the cold and hot models. Fig. 9 gives the time histories of generalized coordinates of the hot model at trajectory points 4 and 5. As can be seen from Fig. 9a to Fig. 9c, at trajectory point 4, the time histories of generalized coordinates vary from convergence to divergence with the decrease of structural stiffness coefficient N, and the critical stiffness coefficient of flutter is about 0.18. Fig. 10 shows the time histories of generalized coordinates of the cold model at the
T. Guo et al. / Aerospace Science and Technology 94 (2019) 105381
7
Table 3 Comparison of the first ten modal frequencies between the cold and hot models. Mode order
Cold model
Hot model Trajectory point 4
1 (first torsion) 2 (first bending) 3 4 5 6 7 8 9 10
Trajectory point 5
Frequency (Hz)
Frequency (Hz)
% decrease
Frequency (Hz)
% decrease
14.527 30.954 76.406 101.44 143.7 158.88 164.81 204.57 220.28 245.96
13.883 23.916 56.694 78.24 111.97 118.28 127.95 164.0 177.15 186.67
4.43% 22.7% 25.8% 22.9% 22.1% 25.5% 22.3% 19.8% 19.6% 24.1%
12.617 16.285 32.461 39.997 54.421 57.776 61.372 82.135 84.576 92.612
13.1% 47.4% 57.5% 60.5% 62.1% 63.6% 62.5% 59.8% 61.6% 62.3%
Fig. 9. Time histories of the generalized coordinates of the hot model.
Table 4 Comparison of the flutter characteristics between the cold and hot models. Case
Critical stiffness coefficient NF Stiffness margin of flutter Flutter dynamic pressure (Pa) Flutter frequency (Hz)
Trajectory point 4
Trajectory point 5
Cold model
Hot model
% decrease
Cold model
Hot model
% decrease
0.066 93.4% 9.56 × 106 24.42
0.18 82% 3.51 × 106 18.01
– 12.2% 63.3% 26.2%
0.082 91.8% 1.1 × 107 21.60
>1 <0
– – – 33.3%
– 14.41
8
T. Guo et al. / Aerospace Science and Technology 94 (2019) 105381
Fig. 10. Time histories of the generalized coordinates of the cold model.
critical flutter point. In all cases including the cold and hot models, a coupled bending-torsion flutter is obtained. Table 4 outlines the flutter characteristics of the cold and hot models. For the hot model, the large reduction in the natural frequencies results in a significant decrease in the flutter margin. At trajectory point 4, the reductions in the stiffness margin, the flutter dynamic pressure and the flutter frequency are up to 12.2%, 63.3% and 26.2%, respectively. At trajectory point 5, the stiffness margin of the cold model is 91.8%, but the heated rudder has fluttered at the original stiffness (N = 1), which can be seen from Fig. 9d.
This work was supported by the National Natural Science Foundation of China (Grant No. 11872212) and a project funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions.
5. Conclusions
References
For the hypersonic aerothermal problem along a flight trajectory, a synchronization algorithm is developed in a non-inertial frame of reference by formulating the governing equations of fluid flow and structural heat transfer into a unified form. The synchronization algorithm is then combined with the CFD/CSD-based time-domain method with variable stiffness. Consequently, the thermal flutter prediction method is achieved. The present method has been applied to the thermal flutter computations of a hypersonic all-movable rudder along a given flight trajectory. From our development process and numerical analyses, we conclude as follows.
[1] J.J. Bertin, Hypersonic Aerothermodynamics, AIAA, Reston, VA, 1994. [2] J.J. McNamara, P.P. Friedmann, Aeroelastic and aerothermoelastic analysis in hypersonic flow: past, present, and future, AIAA J. 49 (2011) 1089–1122. [3] J.J. McNamara, P.P. Friedmann, K.G. Powell, B.J. Thuruthimattam, R.E. Bartels, Aeroelastic and aerothermoelastic behavior in hypersonic flow, AIAA J. 46 (2008) 2591–2610. [4] J.J. Bertin, R.M. Cummings, Fifty years of hypersonics: where we’ve been, where we’re going, Prog. Aerosp. Sci. 39 (2003) 511–536. [5] J. Dugundji, J. Calligeros, Similarity laws for aerothermoelastic testing, J. Aerosp. Sci. 29 (1962) 935–950. [6] N. Lamorte, P.P. Friedmann, Hypersonic aeroelastic and aerothermoelastic studies using computational fluid dynamics, AIAA J. 52 (2014) 2062–2078. [7] R.L. Bisplinghoff, Some structural and aeroelastic considerations of high-speed flight, J. Aeronaut. Sci. 23 (1956) 289–329. [8] H. Ashley, G. Zartarian, Piston theory: a new aerodynamic tool for the aeroelastician, J. Aeronaut. Sci. 23 (1956) 1109–1118. [9] W.W. Zhang, Z.Y. Ye, C.A. Zhang, F. Liu, Supersonic flutter analysis based on a local piston theory, AIAA J. 47 (2009) 2321–2328. [10] M. Van Dyke, A Study of Second-Order Supersonic Flow Theory, NACATR 1081, 1951. [11] D. Zhang, S. Tang, Q.J. Zhu, R.G. Wang, Analysis of dynamic characteristics of the rigid body/elastic body coupling of air-breathing hypersonic vehicles, Aerosp. Sci. Technol. 48 (2016) 328–341. [12] Z.Q. Chen, Y.H. Zhao, R. Huang, Parametric reduced-order modeling of unsteady aerodynamics for hypersonic vehicles, Aerosp. Sci. Technol. 87 (2019) 1–14. [13] C.Q. Gao, W.W. Zhang, Z.Y. Ye, A new viewpoint on the mechanism of transonic single-degree-of-freedom flutter, Aerosp. Sci. Technol. 52 (2016) 144–156. [14] K. Ye, Z.Y. Ye, C.N. Li, J. Wu, Effects of the aerothermoelastic deformation on the performance of the three-dimensional hypersonic inlet, Aerosp. Sci. Technol. 84 (2019) 747–762. [15] F. Chen, H. Liu, S.T. Zhang, Time-adaptive loosely coupled analysis on fluidthermal-structural behaviors of hypersonic wing structures under sustained aeroheating, Aerosp. Sci. Technol. 78 (2018) 620–636. [16] C. Yang, G.S. Li, Z.Q. Wan, Aerothermal-aeroelastic two-way coupling method for hypersonic curved panel flutter, Sci. China, Technol. Sci. 55 (2012) 831–840. [17] S. Zhang, F. Chen, H. Liu, Integrated Fluid-Thermal-Structural Analysis for Predicting Aerothermal Environment of Hypersonic Vehicles, AIAA paper 2014-1394, 2014.
(1) The synchronization algorithm can accurately predict the hypersonic aerothermal problem along a flight trajectory, including the unsteady effects of the variations of flight altitude, Mach number and attack angle. For current rudder along its flight trajectory, a maximum temperature difference at the monitored point is up to 8% between the synchronization and conventional partitioned algorithms. (2) The present thermal flutter computational method can quantitatively evaluate the flutter margin at each point on the flight trajectory. For current rudder at trajectory points 4 and 5, aerodynamic heating reduces the first ten natural frequencies by nearly 20% and 50% except the first torsional mode partly due to neglecting the heating effects on the hinge joint. Consequently, the flutter margin decreases significantly: at trajectory point 4, the reductions in the stiffness margin and the flutter dynamic pressure are up to 12.2% and 63.3%, respectively; at trajectory point 5, the stiffness margin of the cold model is 91.8%, but the heated rudder has fluttered at the original stiffness.
Declaration of competing interest There is no competing interest. Acknowledgements
T. Guo et al. / Aerospace Science and Technology 94 (2019) 105381
[18] E.N. Shen, Z.L. Lu, D. Zhou, T.Q. Guo, Development of a synchronization method for fluid-thermal study of hypersonic flow, Trans. Nanjing Univ. Aeronaut. Astronaut. 35 (2018) 973–985. [19] A.L. Gaitonde, A Dual-Time Method for the Solution of the 2d Unsteady NavierStokes Equations on Structured Moving Meshes, AIAA Paper 1995-1877, 1995. [20] M.S. Liou, A sequel to AUSM: AUSM+ , J. Comput. Phys. 129 (1996) 364–382. [21] P.R. Spalart, S.R. Allmaras, A One-Equation Turbulence Model for Aerodynamic Flows, AIAA paper 1992-0439, 1992. [22] H. Rieger, A. Jameson, Solution of Steady 3-d Compressible Euler and NavierStokes Equations by an Implicit LU Scheme, AIAA Paper 88-0619, 1988. [23] L. Ding, Z.L. Lu, T.Q. Guo, An efficient dynamic mesh generation method for complex multi-block structured grid, Adv. Appl. Math. Mech. 6 (2014) 120–134.
9
[24] W.W. Zhang, Y.W. Jiang, Z.Y. Ye, Two better loosely coupled solution algorithms of CFD based aeroelastic simulation, Eng. Appl. Comput. Fluid Mech. 1 (2007) 253–262. [25] T.Q. Guo, Z.L. Lu, D. Tang, T.G. Wang, L. Dong, A CFD/CSD model for aeroelastic calculations of large-scale wind turbines, Sci. China, Technol. Sci. 56 (2013) 205–211. [26] T.Q. Guo, Z.L. Lu, Y.J. Wu, A time-domain method for transonic flutter analysis with multidirectional coupled vibrations, Mod. Phys. Lett. B 23 (2009) 453–456. [27] T.Q. Guo, Z.L. Lu, A CFD/CSD model for transonic flutter, Comput. Mater. Continua 2 (2005) 105–111.