# Aerothermoelastic Load Calculation for Hypersonic Vehicles Based on Multiphysics Coupled Analysis

Authors: Hao Chen; Min Xu; Zihua Qiu; Dan Xie; Yabin Wang
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

To fulfill the design objective of a structure and thermal protection system, accurate load environment prediction is very important, so we present a high-fidelity aerothermoelastic load calculation method based on a partitioned computational fluid dynamics/computational structural dynamics/computational thermal dynamics (CFD/CSD/CTD) coupling analysis. For the data transformation between the CFD/CSD/CTD systems, finite element interpolation (FEI) is explored, and a shape-preserving grid deformation strategy is achieved via radical basis functions (RBFs). Numerical results are presented for validation of the proposed CFD/CSD/CTD coupling analysis. First, a simply supported panel in hypersonic flow is investigated for results comparison of the proposed coupling method and previous work. Second, a hypersonic forebody is investigated to explore the aerothermoelastic effects while considering the feedback between deformation and aerodynamic heating. The results show that the CFD/CSD/CTD coupling method is accurate for analysis of aerothermoelasticity. In addition, considering the aerothermoelastic effect, the shear force and bending movement increase with time before900s and decrease after 900s, and at 900s increased percentages of 5.7% and 4.1% are observed, respectively. Therefore, it is necessary to adopt high-fidelity CFD/CSD/CTD coupling in the design of a structure and thermal protection system for hypersonic vehicles.


---

## Body

## 1. Introduction

The air-breathing hypersonic vehicle (HSV) has drawn attention from the aerospace communities of several countries because of its possibilities to be fully reusable and to succeed the space shuttle to become the key component of next-generation earth/orbit transportation systems. However, the highly coupled nature of aerodynamics and structural dynamics requires analysis of the complex interactions between aerodynamic forces, aerodynamic heating, heat transfer, and structural deformation during the design and verification of an HSV. Generally, the forebody serves as both a compressor and a lift component; i.e., the forebody generates the proper inflow for the inlet and generates a considerable fraction of the lift. Thus, deformation of the forebody may affect the lift and thrust characteristics significantly. During the exposure to hypersonic flow, the forebody must withstand intense aerodynamic heating and severe aerodynamic forces. The aerodynamic heating degrades the material properties and induces high temperatures in the structure. The temperature gradients and structural constraints generate thermal stresses that alter the structural geometry significantly. In addition, the aerodynamic forces can cause the structural geometry to deform; such deformation redistributes the aerodynamic forces and heating. These complex interactions between aerodynamics, structural dynamics, and heat transfer constitute the aerothermoelastic problem shown in Figure1. Therefore, the design of the structure and thermal protection system (TPS) of an HSV is critical so that the aerothermoelastic problems can be analyzed properly and the load environment predicted accurately.Figure 1
Basic interactions of the aerothermoelastic problem.Considering the deformation induced by aerodynamic forces and thermal effects whose feedback would alter aerodynamic heating, the CFD/CSD/CTD coupling offers the highest fidelity to calculate the aerothermoelastic loads. In the 1980s, for the first time, Thornton [1] adopted this method, which was implemented using the Taylor-Galerkin algorithm and was used to obtain the transient temperature distributions and the corresponding thermal stress distributions of a leading edge. In [2], the aerothermoelastic responses of a metallic thermal protection system (MTPS) panel were calculated using the coupling and uncoupling method. The structural deformation and aerodynamic heating induced by deformation were negligible, whereas the change of surface temperature was significant. The lateral temperature gradient increases the thermal stress of the plate; however, the deformation does not dramatically alter the thermal load along the trajectory. In [3], two heat transfer mechanisms were investigated for the heated panel in hypersonic flow; the aerodynamic heating via fluid viscosity was found to dominate over the heating from unsteady local flow. Ref. [4] focused on the wall shear stress effect in the aerothermoelastic study of a skin panel, and the wall shear stress effect was found to be negligible in a linear analysis, whereas in a nonlinear analysis, the effect is significant, and the transient chaotic/dynamic instability is diminished. In [5], the aerothermoelastic characteristics of a low-aspect-ratio wing were studied using CFD/CSD/CTD coupling, and significant reductions in the flutter boundary of the heated wing were observed. The wing was also found to be susceptible to thermal buckling caused by thermal stresses. The research method adopted in [3, 4], which is based on a simplified aerodynamic model (piston theory) and a structural model (von Karman theory), is challenging to apply to a more general configuration, such as a forebody. Although [5] used the more advanced CFD/CSD/CTD coupling method, the feedback effect of the structural deformation to the aerodynamic heating was not considered.In the 21st century, Tran and Farhat [6] formulated the first nonlinear multiphysics aerothermoelastic framework based on CFD/CSD/CTD coupling completely by adding heat conduction equation and modifying the structural equation of the CFD/CSD coupling based aeroelastic framework established by Farhat and coworkers [7] and applied to an F-16 fighter in [8, 9]. Then the framework was extended to data-driven modeling and simulation of nonlinear degrading continuum systems by Michopoulos and Farhat [10, 11]. The works of [10, 11] were very meaningful as the near real-time material and structural health monitoring systems based on embedded sensor networks and associated sensor-based technologies could be implemented. However, the present paper focuses on the aerothermoelastic problem rather than the sensor related problem, so let us rewind to the aerothermoelastic governing equations. The governing equations of [6] are represented in mixed form. The fluid equation is a partial differential equation, while the structural equation and heat conduction equation are finite element discretized matrix equations. Although the structural equation considering thermal effects is derived in [10, 11], the heat conduction equation is just listed in [10–12] and one term in the governing equation is duplicated with the boundary condition which will be discussed later. In addition to these three governing equations, one each for a physical field, there is another pseudo-structural equation in the aerothermoelastic framework of [6, 10, 11] which is used for fluid grid deformation. Unlike the three irreplaceable governing equations for physical fields, the pseudo-structural equation can be replaced by other grid deformation methods. As the newly developed Radial Basis Functions (RBFs) methods maintain grid orthogonality well [13] which is very important for aerodynamic heating calculation and has been applied to aeroelastic problem successfully [14], the RBFs method with some modification other than the pseudo-structural method is used in this paper.Since the aerothermoelastic framework has only been used to several supersonic cases in [6, 10, 11], some new difficulties and phenomena may emerge when the framework is applied to hypersonic vehicles, especially for the air-breathing vehicles that can maintain hypersonic flight for more than one hour. Note that the three governing equations for physical fields are all time-dependent, suggesting that reasonable results of coupling solution can only be obtained by time marching with a physical time step smaller than the smallest characteristic time of the three physical fields. The characteristic time of fluid equation equates the characteristic length divided by the fluid speed and is in the order of 1e-3. Then the aerothermoelastic solution of a hypersonic vehicle with flight time more that one hour would need 1e6 iterations at least and that imposes high computational cost. To overcome the computational difficulty, this paper borrows the ideas of [1] by adopting the quasi-static assumption for the fluid and structure equation of which the characteristic time is small. Then the strict limitation is removed and time step could be much larger than the order of 1e-3, such as 30s in [1]. To investigate hypersonic vehicle with a long flight time for possible extreme load in the aerodynamic heating process, the heat conduction equation must be left in time-dependent form.The present study is based on the high-fidelity CFD/CSD/CTD coupling method, considering the feedback of the structural deformation. The CSD solver implements a second-order nonlinear finite element method, by which the shear locking can be avoided. In the CTD solver, the lumped heat capacity matrix is used to resolve the nonphysical oscillation problem. To preserve the accuracy of the displacements solved by the CSD solver, a finite element interpolation (FEI) is employed as the data transfer method. Moreover, for retaining the orthogonality of the CFD grid in the coupling computation, the shape-preserving grid deformation strategy based on RBFs is adopted. Next, the correctness of the CFD/CSD/CTD coupling method is verified in the case of a simply supported panel within hypersonic flow. Finally, the aerothermoelastic calculation is conducted for a forebody of an HSV using the high-fidelity coupling method; the structural deformation, temperature field, stress distribution, and loads are obtained. Thus, the influence of the introduction of aerodynamic heat coupling is determined by comparing the results of the CFD/CSD/CTD coupling method with the results of the CFD/CSD coupling method without the consideration of aerodynamic heating.


## 2. Solution Method

### 2.1. Aerothermoelastic Coupling Framework

The first step to solve the aerothermoelastic problem is to model it. Michopoulos and Farhat [12] provided the following governing equations:(1)∂Jw∂tξ+J∇x·Fw-∂x∂tξw=J∇x·Rw+JSw,(2)ρS∂2uS∂t2-∇xσSεS,TS=b,(3)ρSCpS∂TS∂t+∇x·-λS∇xTS-σS∂εS∂t-Qi-QΓ=0.Eq. (1) represents an alternative form of the Navier-Stokes (NS) equations describing the fluid dynamics, derived in an arbitrary Lagrangian-Eulerian (ALE) reference frame, and (2) and (3) are the structural dynamics and heat conduction equations in the structural domain, respectively.In (1), w is the conservative fluid state vector, F and R represent the convective and viscous terms, S(w) is a source term, x(ξ,t) is the position vector of fluid points, and J=∂x/∂ξ refers to the Jacobian determinant of the deformation gradient. In (2) and (3), ρS represents the material density; σS and εS refer to the structural stress and strain tensors, respectively; uS is the displacement vector of structural points; b is the body force vector; CpS and λS are the specific heat and thermal conductivity of the material, respectively; TS is the temperature field in structure; and Qi and QΓ are the internal heat flux and the heat flux at the fluid-structure interface Γ [12], respectively. The internal heat flux Qi expressed in [12] stands for the effects such as phase transformation, electromagnetic heating, and electrical resistance heating. None of the effects are considered in this paper, so Qi is neglected.Note that the third term in (3), which is an internal heat source term proportional to the stress and strain rate, expressed as σS∂εS/∂t, reflects the internal heating effects induced by thermoelastic coupling within the structure. Regarding the fluid-structure coupling, the following boundary conditions should be satisfied [12]:(4)σS·n=-pn+σF·n,(5)∂x∂t=∂uS∂t=vF,(6)xt-xtt=0=uF,(7)TS=TF,(8)λS∇TS·n+ε^σ^TS4-Te4=λF∇TF·n=QΓ,F.Note that the first term λS∇TS·n of (8) is just QΓ of (3), suggesting that the first term in boundary condition (8) is duplicated with QΓ in governing equation (3). So QΓ of (3) is removed to correct the error.In this study, considering the large differences in the characteristic time of each physical field of the aerothermoelastic problem (the characteristic time of heat conduction is much longer than the characteristic time of aerodynamic and structural dynamic response), a quasi-static solution is used; i.e., the time-dependent term is retained in the heat conduction equation to consider the transient effect, and for the aerodynamic and structural dynamics equations, only the steady characteristics are considered. Therefore, the grid velocity in (1) is ignored, but the deformation must be updated in the coupling iterations. For the same reason, the time-dependent thermoelastic coupling term expressed as σS∂εS/∂t is also ignored. Moreover, gravity is the only body force considered in this paper. The detailed solutions of the aerodynamic, structural, and heat transfer equations are described in Sections 2.2, 2.3, and 2.4, respectively.The iterative solution of the aerothermoelastic coupling problem is shown in Figure2 [6]. When the iterative solution is obtained, the quasi-static aerodynamic equation is first solved, and the obtained heat flux on surface is passed to the heat conduction solver through the two-dimensional interpolation algorithm. The obtained surface pressure is also passed to the structural response solver by the same algorithm. Next, the unsteady heat conduction equation is solved under the new heat flux boundary condition, and then the solved temperature field is passed to the structural response solver through the three-dimensional interpolation algorithm. The obtained temperature distribution on the surface is passed to the aerodynamic solver through the two-dimensional interpolation algorithm. Finally, the structural equilibrium equation is solved in a quasi-static manner under the new surface pressure and temperature field, and then the obtained deformation is passed to the aerodynamic solver through the dynamic mesh method to complete the iteration of a coupling step. In the process, both two-dimensional and three-dimensional interpolation are achieved by FEI uniformly. The dynamic mesh method is implemented using a shape-preserving grid deformation strategy based on RBFs to maintain the orthogonality and surface accuracy of the original CFD grid. The FEI and the RBF-based dynamic mesh method used in this study are described in Sections 2.5 and 2.6.Figure 2
Aerothermoelastic coupling iterations [6].


### 2.2. Solution of the Aerodynamic Equation

The NS equation is solved to obtain the aerodynamic forces and heating using the finite volume method (FVM). The NS equation used is given by (9), in which the time-dependent term is retained for pseudo-time advancing:(9)∂∂t∭ΩwdV+∬∂ΩF-R•ndS=0.In (9), w is the conservation variable vector, F is the convective term, R is the viscous term, and Ω and ∂Ω are the control body and its boundary, respectively. When (9) is solved, the convective term is calculated using the ausmpw+ scheme because of its accuracy of aerodynamic forces and heating in the hypersonic regime, and the minmod limiter is chosen because of its robustness. The viscous term is calculated using the central difference scheme. The LUSGS scheme is used as the pseudo-time advancing method, and the convergence is accelerated by the adoption of a local time step. The boundary conditions in the far field are determined according to the direction of local flow velocity; i.e., the variables of ghost cells at the boundary are assigned to the inflow values when air flows inward, and the assignment is performed by extrapolation of the inner cells when air flows outward. On the wall boundary, the non-slip condition is applied to the velocity, the normal gradient of the pressure in the surface is enforced to be zero, the Dirichlet condition is used for the temperature, and the wall temperature is given by the interpolation algorithm described below.


### 2.3. Solution of the Structural Equation

The structural equilibrium equation is solved by the finite element method. The static equilibrium equation without the time-dependent term used in this study is(10)∇xσSεS,TS=g,where g stands for gravity.Because the effect of thermal effects must be considered in solving the aerothermoelastic coupling problem, incorporation of the nonlinearity caused by the thermal effect is necessary. The thermal nonlinear effects can be divided into two types: (1) changes in the material properties caused by temperature changes and (2) changes in the stiffness caused by thermal stress. The combination of the two can be written as follows:(11)K=KE+Kσ.K E is derived from the first-order strain without changing the linear property of the geometric equation. Kσ is introduced by the second-order strain, similar to the geometric nonlinearity. Considering that the structural deformation may be larger under the combination of aerodynamic forces and aerodynamic heating, the updated Lagrangian (U.L.) frame with geometric nonlinearity is used to solve the equilibrium equation in which the thermal effect is introduced. The final incremental equilibrium equation at time t+1 from the moment t is obtained as follows:(12)KttL+KttNLΔu=Ptt+Δt-Itt,where KttL is the linear stiffness matrix calculated by the material properties with temperature changes and is equivalent to KE in (11); KttNL is a nonlinear stiffness matrix introduced by the second-order strain, where the thermal effects are taken into account by subtracting the thermal strain from the whole strain; Δu is the displacement increment vector; Ptt+Δt is the load vector integrated from the external force; Itt is the load vector integrated from the internal force. The (KttL+KttNL) on the left side is the tangent stiffness matrix at time t, and Ptt+Δt-Itt on the right side is the load vector of the unbalanced forces. The pressure and stress distributions required to calculate Ptt+Δt are given by the interpolation algorithm described below. In this paper, the Almansi strain is used in assembling the nonlinear stiffness matrix KttNL and (12) is solved by using the LU decomposition for sparse matrix.To avoid shear locking, the second-order hexahedral element with 20 nodes is adopted in this study for structural modeling. The shape function of this kind of hexahedral element is expressed as follows:(13)Ni=1+ξiξ1+ηiη1+ζiζξiξ+ηiη+ζiζ-2ξi2ηi2ζi28+1-ξ21+ηiη1+ζiζ1-ξi2ηi2ζi24+1+ξiξ1-η21+ζiζ1-ηi2ξi2ζi24+1+ξiξ1+ηiη1-ζ21-ζi2ξi2ηi24,i=1,2,…,20


### 2.4. Solution of the Heat Conduction Equation

The heat conduction problem is solved by the finite element method; that is, the temperature distribution of the internal and the surface structure is obtained by solving the Fourier equation. The Fourier equation without the internal heat source used in this study is represented by(14)ρSCpS∂TS∂t+∇x·-λS∇xTS=0.The first step to solve this equation is spatial discretization via the finite element method; the obtained semidiscrete equation is expressed as follows:(15)CTT˙+KTT=PT,where CT and KT are the heat capacity matrix and the heat conduction matrix, respectively, and PT is the thermal load vector. The heat flux distributions required to calculate PT are given by the interpolation algorithm described below.Between timestn and tn+1, a numerical difference is introduced for the time-dependent term of (15); for the other terms, linear interpolation is introduced. The resulting expression is written as(16)CTTn+1-TnΔt=θPT,n+1+1-θPT,n-KTθTn+1+1-θTn,and it is ultimately represented as(17)CT+θΔtKTTn+1=θPT,n+1+1-θPT,n+CT-1-θΔtKTTn.To ensure the stability in time,θ is assigned to 1. Eq. (17) is then solved by the preconditioned conjugate gradient method based on the element-by-element (EBE) technique to obtain the temperature distribution at the next moment. To avoid nonphysical oscillation, the lumped heat capacity matrix is adopted.The 8-node, 3-dimensional solid hexahedron is used in this study for heat conduction calculation and its shape function is represented by(18)Ni=1+ξiξ1+ηiη1+ζiζ8,i=1,2,…,8.


### 2.5. Data Transfer by FEI

The FEI used in this study for all of the data transfer incorporates the technique of numerical inverse isoparametric mapping from [15–17] and is implemented using forward FEI after local coordinates obtained by an iterative algorithm. To avoid duplicate computation, the local coordinates are calculated and stored in a preprocessing stage, and only forward FEI is run in the time marching simulation.In a 3D quadrilateral element which consists of curved surface, the isoparametric mapping is defined as mapping from local coordinate vectorξ=(ξ,η)T to global coordinates x=(x,y,z)T; that is, (19)x=∑i=1ncellNiξxi,where suffix i is the node index of the element. Thus, a general variable w in this element can be expressed as(20)w=∑i=1ncellNiξwi.Although (19) is not complex, the inverse mapping is not straightforward. Thus, an iterative algorithm is adopted to obtain (ξ,η) from (x,y,z). Note that (19) is formed by three equations with only two independent variables to solve. The equations are solved using the Newton iteration method with the least square method.After the local coordinates vectorξ is calculated, the FEI of w is done by (20).


### 2.6. Shape-Preserving Grid Deformation Strategy Based on RBFs

Because the CFD method is used to calculate the aerodynamic heating, the orthogonality and accuracy of the deformed grid should be ensured. To maintain grid orthogonality, the RBFs [14, 18] method is used, where the radial basis function is Wendland’s C2 [13] function. For the precise transfer of the surface deformation from structure to fluid, the point selection process is not implemented. Instead, all of the points on the fluid-structure interface are used in the RBF method. After RBF deformation is performed, a shape-preserving grid deformation strategy is applied to aerodynamic surface to maintain the consistency of interfaces on both sides.Take the transfer of deformation inX-direction as an example; that is,(21)Δxs=ΦXs,XsWx,(22)Δxd=ΦXd,XsWx,where Δxd and Δxs are the x deformation of the fluid and the structural grid, respectively, Wx refers to the RBF interpolation coefficient vectors, and Xs and Xd indicate the original surface coordinates of the structural mesh and the fluid grid, respectively. Φ is a matrix formed by functions and is given by(23)ΦXd,Xs=ϕd11d0ϕd21d0⋯ϕdn1d0ϕd12d0ϕd22d0⋯ϕdn2d0⋯ϕd1md0ϕd2md0⋯ϕdnmd0m×n,where dij indicates the distance between the i-th structural point and j-th fluid point.Finally, the RBF interpolation can be summarized as follows: first, a linear algebraic system is obtained by solving (21), and then the solved Wx is substituted into (22) to obtain the deformation of the interpolated fluid grid in the X-direction. The deformations in the Y- and Z-directions are similar, whereas the interpolation coefficients are different.Afterwards, all the nodes of aerodynamic surface mesh are projected to the deformed structural surface by the FEI method described in Section2.5. Next, the updated deformation obtained by FEI is distributed into the complete fluid grid using the TFI [19] method.


## 2.1. Aerothermoelastic Coupling Framework

The first step to solve the aerothermoelastic problem is to model it. Michopoulos and Farhat [12] provided the following governing equations:(1)∂Jw∂tξ+J∇x·Fw-∂x∂tξw=J∇x·Rw+JSw,(2)ρS∂2uS∂t2-∇xσSεS,TS=b,(3)ρSCpS∂TS∂t+∇x·-λS∇xTS-σS∂εS∂t-Qi-QΓ=0.Eq. (1) represents an alternative form of the Navier-Stokes (NS) equations describing the fluid dynamics, derived in an arbitrary Lagrangian-Eulerian (ALE) reference frame, and (2) and (3) are the structural dynamics and heat conduction equations in the structural domain, respectively.In (1), w is the conservative fluid state vector, F and R represent the convective and viscous terms, S(w) is a source term, x(ξ,t) is the position vector of fluid points, and J=∂x/∂ξ refers to the Jacobian determinant of the deformation gradient. In (2) and (3), ρS represents the material density; σS and εS refer to the structural stress and strain tensors, respectively; uS is the displacement vector of structural points; b is the body force vector; CpS and λS are the specific heat and thermal conductivity of the material, respectively; TS is the temperature field in structure; and Qi and QΓ are the internal heat flux and the heat flux at the fluid-structure interface Γ [12], respectively. The internal heat flux Qi expressed in [12] stands for the effects such as phase transformation, electromagnetic heating, and electrical resistance heating. None of the effects are considered in this paper, so Qi is neglected.Note that the third term in (3), which is an internal heat source term proportional to the stress and strain rate, expressed as σS∂εS/∂t, reflects the internal heating effects induced by thermoelastic coupling within the structure. Regarding the fluid-structure coupling, the following boundary conditions should be satisfied [12]:(4)σS·n=-pn+σF·n,(5)∂x∂t=∂uS∂t=vF,(6)xt-xtt=0=uF,(7)TS=TF,(8)λS∇TS·n+ε^σ^TS4-Te4=λF∇TF·n=QΓ,F.Note that the first term λS∇TS·n of (8) is just QΓ of (3), suggesting that the first term in boundary condition (8) is duplicated with QΓ in governing equation (3). So QΓ of (3) is removed to correct the error.In this study, considering the large differences in the characteristic time of each physical field of the aerothermoelastic problem (the characteristic time of heat conduction is much longer than the characteristic time of aerodynamic and structural dynamic response), a quasi-static solution is used; i.e., the time-dependent term is retained in the heat conduction equation to consider the transient effect, and for the aerodynamic and structural dynamics equations, only the steady characteristics are considered. Therefore, the grid velocity in (1) is ignored, but the deformation must be updated in the coupling iterations. For the same reason, the time-dependent thermoelastic coupling term expressed as σS∂εS/∂t is also ignored. Moreover, gravity is the only body force considered in this paper. The detailed solutions of the aerodynamic, structural, and heat transfer equations are described in Sections 2.2, 2.3, and 2.4, respectively.The iterative solution of the aerothermoelastic coupling problem is shown in Figure2 [6]. When the iterative solution is obtained, the quasi-static aerodynamic equation is first solved, and the obtained heat flux on surface is passed to the heat conduction solver through the two-dimensional interpolation algorithm. The obtained surface pressure is also passed to the structural response solver by the same algorithm. Next, the unsteady heat conduction equation is solved under the new heat flux boundary condition, and then the solved temperature field is passed to the structural response solver through the three-dimensional interpolation algorithm. The obtained temperature distribution on the surface is passed to the aerodynamic solver through the two-dimensional interpolation algorithm. Finally, the structural equilibrium equation is solved in a quasi-static manner under the new surface pressure and temperature field, and then the obtained deformation is passed to the aerodynamic solver through the dynamic mesh method to complete the iteration of a coupling step. In the process, both two-dimensional and three-dimensional interpolation are achieved by FEI uniformly. The dynamic mesh method is implemented using a shape-preserving grid deformation strategy based on RBFs to maintain the orthogonality and surface accuracy of the original CFD grid. The FEI and the RBF-based dynamic mesh method used in this study are described in Sections 2.5 and 2.6.Figure 2
Aerothermoelastic coupling iterations [6].


## 2.2. Solution of the Aerodynamic Equation

The NS equation is solved to obtain the aerodynamic forces and heating using the finite volume method (FVM). The NS equation used is given by (9), in which the time-dependent term is retained for pseudo-time advancing:(9)∂∂t∭ΩwdV+∬∂ΩF-R•ndS=0.In (9), w is the conservation variable vector, F is the convective term, R is the viscous term, and Ω and ∂Ω are the control body and its boundary, respectively. When (9) is solved, the convective term is calculated using the ausmpw+ scheme because of its accuracy of aerodynamic forces and heating in the hypersonic regime, and the minmod limiter is chosen because of its robustness. The viscous term is calculated using the central difference scheme. The LUSGS scheme is used as the pseudo-time advancing method, and the convergence is accelerated by the adoption of a local time step. The boundary conditions in the far field are determined according to the direction of local flow velocity; i.e., the variables of ghost cells at the boundary are assigned to the inflow values when air flows inward, and the assignment is performed by extrapolation of the inner cells when air flows outward. On the wall boundary, the non-slip condition is applied to the velocity, the normal gradient of the pressure in the surface is enforced to be zero, the Dirichlet condition is used for the temperature, and the wall temperature is given by the interpolation algorithm described below.


## 2.3. Solution of the Structural Equation

The structural equilibrium equation is solved by the finite element method. The static equilibrium equation without the time-dependent term used in this study is(10)∇xσSεS,TS=g,where g stands for gravity.Because the effect of thermal effects must be considered in solving the aerothermoelastic coupling problem, incorporation of the nonlinearity caused by the thermal effect is necessary. The thermal nonlinear effects can be divided into two types: (1) changes in the material properties caused by temperature changes and (2) changes in the stiffness caused by thermal stress. The combination of the two can be written as follows:(11)K=KE+Kσ.K E is derived from the first-order strain without changing the linear property of the geometric equation. Kσ is introduced by the second-order strain, similar to the geometric nonlinearity. Considering that the structural deformation may be larger under the combination of aerodynamic forces and aerodynamic heating, the updated Lagrangian (U.L.) frame with geometric nonlinearity is used to solve the equilibrium equation in which the thermal effect is introduced. The final incremental equilibrium equation at time t+1 from the moment t is obtained as follows:(12)KttL+KttNLΔu=Ptt+Δt-Itt,where KttL is the linear stiffness matrix calculated by the material properties with temperature changes and is equivalent to KE in (11); KttNL is a nonlinear stiffness matrix introduced by the second-order strain, where the thermal effects are taken into account by subtracting the thermal strain from the whole strain; Δu is the displacement increment vector; Ptt+Δt is the load vector integrated from the external force; Itt is the load vector integrated from the internal force. The (KttL+KttNL) on the left side is the tangent stiffness matrix at time t, and Ptt+Δt-Itt on the right side is the load vector of the unbalanced forces. The pressure and stress distributions required to calculate Ptt+Δt are given by the interpolation algorithm described below. In this paper, the Almansi strain is used in assembling the nonlinear stiffness matrix KttNL and (12) is solved by using the LU decomposition for sparse matrix.To avoid shear locking, the second-order hexahedral element with 20 nodes is adopted in this study for structural modeling. The shape function of this kind of hexahedral element is expressed as follows:(13)Ni=1+ξiξ1+ηiη1+ζiζξiξ+ηiη+ζiζ-2ξi2ηi2ζi28+1-ξ21+ηiη1+ζiζ1-ξi2ηi2ζi24+1+ξiξ1-η21+ζiζ1-ηi2ξi2ζi24+1+ξiξ1+ηiη1-ζ21-ζi2ξi2ηi24,i=1,2,…,20


## 2.4. Solution of the Heat Conduction Equation

The heat conduction problem is solved by the finite element method; that is, the temperature distribution of the internal and the surface structure is obtained by solving the Fourier equation. The Fourier equation without the internal heat source used in this study is represented by(14)ρSCpS∂TS∂t+∇x·-λS∇xTS=0.The first step to solve this equation is spatial discretization via the finite element method; the obtained semidiscrete equation is expressed as follows:(15)CTT˙+KTT=PT,where CT and KT are the heat capacity matrix and the heat conduction matrix, respectively, and PT is the thermal load vector. The heat flux distributions required to calculate PT are given by the interpolation algorithm described below.Between timestn and tn+1, a numerical difference is introduced for the time-dependent term of (15); for the other terms, linear interpolation is introduced. The resulting expression is written as(16)CTTn+1-TnΔt=θPT,n+1+1-θPT,n-KTθTn+1+1-θTn,and it is ultimately represented as(17)CT+θΔtKTTn+1=θPT,n+1+1-θPT,n+CT-1-θΔtKTTn.To ensure the stability in time,θ is assigned to 1. Eq. (17) is then solved by the preconditioned conjugate gradient method based on the element-by-element (EBE) technique to obtain the temperature distribution at the next moment. To avoid nonphysical oscillation, the lumped heat capacity matrix is adopted.The 8-node, 3-dimensional solid hexahedron is used in this study for heat conduction calculation and its shape function is represented by(18)Ni=1+ξiξ1+ηiη1+ζiζ8,i=1,2,…,8.


## 2.5. Data Transfer by FEI

The FEI used in this study for all of the data transfer incorporates the technique of numerical inverse isoparametric mapping from [15–17] and is implemented using forward FEI after local coordinates obtained by an iterative algorithm. To avoid duplicate computation, the local coordinates are calculated and stored in a preprocessing stage, and only forward FEI is run in the time marching simulation.In a 3D quadrilateral element which consists of curved surface, the isoparametric mapping is defined as mapping from local coordinate vectorξ=(ξ,η)T to global coordinates x=(x,y,z)T; that is, (19)x=∑i=1ncellNiξxi,where suffix i is the node index of the element. Thus, a general variable w in this element can be expressed as(20)w=∑i=1ncellNiξwi.Although (19) is not complex, the inverse mapping is not straightforward. Thus, an iterative algorithm is adopted to obtain (ξ,η) from (x,y,z). Note that (19) is formed by three equations with only two independent variables to solve. The equations are solved using the Newton iteration method with the least square method.After the local coordinates vectorξ is calculated, the FEI of w is done by (20).


## 2.6. Shape-Preserving Grid Deformation Strategy Based on RBFs

Because the CFD method is used to calculate the aerodynamic heating, the orthogonality and accuracy of the deformed grid should be ensured. To maintain grid orthogonality, the RBFs [14, 18] method is used, where the radial basis function is Wendland’s C2 [13] function. For the precise transfer of the surface deformation from structure to fluid, the point selection process is not implemented. Instead, all of the points on the fluid-structure interface are used in the RBF method. After RBF deformation is performed, a shape-preserving grid deformation strategy is applied to aerodynamic surface to maintain the consistency of interfaces on both sides.Take the transfer of deformation inX-direction as an example; that is,(21)Δxs=ΦXs,XsWx,(22)Δxd=ΦXd,XsWx,where Δxd and Δxs are the x deformation of the fluid and the structural grid, respectively, Wx refers to the RBF interpolation coefficient vectors, and Xs and Xd indicate the original surface coordinates of the structural mesh and the fluid grid, respectively. Φ is a matrix formed by functions and is given by(23)ΦXd,Xs=ϕd11d0ϕd21d0⋯ϕdn1d0ϕd12d0ϕd22d0⋯ϕdn2d0⋯ϕd1md0ϕd2md0⋯ϕdnmd0m×n,where dij indicates the distance between the i-th structural point and j-th fluid point.Finally, the RBF interpolation can be summarized as follows: first, a linear algebraic system is obtained by solving (21), and then the solved Wx is substituted into (22) to obtain the deformation of the interpolated fluid grid in the X-direction. The deformations in the Y- and Z-directions are similar, whereas the interpolation coefficients are different.Afterwards, all the nodes of aerodynamic surface mesh are projected to the deformed structural surface by the FEI method described in Section2.5. Next, the updated deformation obtained by FEI is distributed into the complete fluid grid using the TFI [19] method.


## 3. Results and Discussions

### 3.1. Aerodynamic Coefficients and Pressure Verification

The blunt biconic model in [20] is used to verify the accuracy of the CFD solver in aerodynamic coefficients and pressure calculation. The aerodynamic grid of the blunt biconic model is shown in Figure 3. The parameters used for the current verification study are listed in Table 1.Table 1
The computational parameters for the blunt biconic.
M∞ p ∞ / P a T ∞ / K α / ∘ 6 2200 61.4 5Figure 3
The aerodynamic grid of the blunt biconic.The results of the normal-force coefficient (CN), axial-force coefficient (CA), and pitching-moment coefficient (Cm) obtained by the present CFD solver and the ones from the experiment of [20] are listed and compared in Table 2. The maximum error is 2.04% for the pitching-moment coefficient which means the accuracy of the CFD solver in aerodynamic coefficients is good.Table 2
The comparison of the aerodynamic coefficients of the blunt biconic.
C N C A C m CFD 0.156 0.100 -0.0049 Exp 0.153 0.102 -0.005 Error 1.96% 1.96% 2.04%The pressure distribution (Ps) in the longitudinal symmetric plane obtained by the CFD solver and the one from the experiment of [20] are compared in Figure 4. Figure 4 shows the CFD solutions match well with the experimental results.Figure 4
The comparison of the pressure distribution in the longitudinal symmetric plane.


### 3.2. Aerodynamic Heating Verification

The compressible laminar boundary layer of the plate (1.8m) is used to verify the accuracy of the CFD solver in aerodynamic heating calculation. The aerodynamic grid of the plate is shown in Figure 5. The parameters used for the current verification study are listed in Table 3.Table 3
The computational parameters for the blunt biconic.
M ∞ H ∞ / k m T w / K 5 20 288.15Figure 5
The aerodynamic grid of the plate.The aerodynamic heating distribution by the CFD solvers is compared with the theoretical result of [21] and the solution of Eckert’s reference method of [22] in Figure 6. The figure shows the CFD solver has a better accuracy in aerodynamic heating calculation than Eckert’s reference method compared to the theoretical method.Figure 6
The comparison of aerodynamic heating distribution by different methods.


### 3.3. Heat Conduction Verification

The heated panel model (1m×1m×0.1m) is used to verify the accuracy of the CTD solver in transient temperature calculation. The panel is initially at a uniform temperature of 30K. The upper and right surfaces are isotherm boundary at 100K, while the lower and left surfaces are adiabatic. The grid for heat conduction calculation is shown in Figure 7. The material parameters used for the current verification study are listed in Table 4.Table 4
The material parameters for the heated panel.
λ ρ C p /W/(m⋅K) / k g / m 3 / J / ( k g ⋅ K ) 0.5 3 2Figure 7
The thermal grid of the heated panel.The temperature distributions at1s by analytic method and the present CTD solver are compared in Figure 8. Note that the CTD solver uses a time step of Δt=0.01s and the time matching is carried out for 100 steps. The results of Figure 8 shows that the temperature distribution by the CTD solver agrees well with the analytic solution.Figure 8
The comparison of the temperature distributions at1s.
(a)
By analytic method (b)
By CTD solver


### 3.4. Structural Displacement Verification

The simply supported panel model (1m×1m×0.1m) used to verify the accuracy of the CSD solver in displacement calculation is shown in Figure 9. The panel is initially at a uniform temperature of 300K. The material parameters used for the current verification study are listed in Table 5. The model is uniformly discrete by 120 (60×1×2) second-order hexahedral elements with 20 nodes.Table 5
The material parameters for the simply supported panel.
E ν α ρ /MPa  / K - 1 / k g / m 3 1 0.3 1e-6 2700Figure 9
The simple supported panel model.The mid-plate displacements obtained by the CSD solver based on nonlinear FEM and von Karman theory in [23] are compared in Figure 10. The displacement in Figure 10(a) is calculated under the pressure qa=1Pa, while the displacement in Figure 10(b) is calculated under the pressure qa=1Pa and the heated temperature T=1300K. No matter under which condition, the FEM results agree well with the von Karman theory.Figure 10
The comparison of the mid-plate displacements.
(a)
q a = 1 P a (b)
q a = 1 P a and T=1300K


### 3.5. Aerothermoelastic Panel Verification

To verify the CFD/CSD/CTD coupling method, the multilayer flat panel of [23] is used as the testing model. There are three layers in the panel, namely, the radiation shield, the thermal insulation, and the main structure. The material and thickness of each layer are listed in Table 6, and the thermophysical properties for a temperature of 300K are listed in Table 7. The emissivity of the radiation shield is assumed to be constant at 0.7. In this model, only the stiffness of the main structure is considered, and the stiffnesses of the other two layers are neglected. The panel is supported by immovable, simple supports.Table 6
Composition of the panel.
Layer Thickness/mm Material Radiation shield 7.4 PM-2000 Thermal insulation 10 IMI Main structure 5 Ti-6Al-2Sn-4Zr-2MoTable 7
Properties of the panel at 300 K.
Layer Density Heat capacity Thermal conductivity /kg/m3 /J/(kg⋅K) /W/(m⋅K) Radiation shield 359 465 0.250 Thermal insulation 73.0 729 0.0258 Main structure 4540 463 6.89The panel is located on a5∘ inclined surface of a wedge-shaped forebody. The length of the panel is 1m. The test condition is at altitude 30km and Mach number 8. The initial temperature is assumed to be 300K.The results of the temperature distribution of the upper surface and the displacement at1200s after 60,000 steps of aerothermoelastic coupling iteration are compared with the results of [23] in Figures 11 and 12, respectively. From the comparison, it is evident that the aerothermoelastic coupling method of this study is sufficiently accurate to investigate effects of the introduction of aerodynamic heat coupling.Figure 11
Comparison of temperature distribution of the upper surface (t=1200s).Figure 12
Comparison of displacement.


### 3.6. Aerothermoelastic Forebody Study

After the above verification, the aerothermoelastic load calculation method based on CFD/CSD/CTD coupling is applied to the two-dimensional forebody as shown in Figure13. This model is in the XY plane, and the total length in the X-direction is 1m. The upper surface is parallel to the X-direction, and the angle between lower surface and the X-direction is 10∘. For the sake of thermal protection, the leading edge between the upper surface and the lower surface is rounded by a circle with a 10mm radius. The fluid grid of the forebody is shown in Figure 14 and there are 33320 cells in this grid. Note that cells are clustered near the surface of the forebody.Figure 13
Geometry of the forebody.Figure 14
Aerodynamic grid of the forebody.The finite element model for heat conduction using 8-node, 3-dimensional solid hexahedron elements is shown in Figure15. The total solid element number is 6900 and the total number of nodes is 13662. As shown in Figure 16, there are 1662 cells and 12612 nodes in the structural model formed by 20-node, 3-dimensional solid hexahedron elements. The constraints are imposed at the rear of the structural grid. The material used is a type of high-temperature titanium alloy designated Ta19. The emissivity of the forebody is assumed to be constant at 0.8.Figure 15
Heat conduction grid of the forebody.Figure 16
Structural grid of the forebody.There are two types of loads studied in this work, namely, the shear force and bending moment. The positive direction of the shear force is defined as the positiveY-direction. The bending moment is integrated in the negative Z-direction (leading edge upward), and the integral point is at the rear of the forebody with the same Y-coordinate as the center of circle at leading edge. As the model used is two-dimensional, the loads are integrated by setting the length of the third direction to be 5mm.


#### 3.6.1. Cases Description

Table8 lists three cases (A-1, A-2, A-3) considered in this study. The flight condition of the three cases is Mach number 5, altitude 20km, and angle of attack 0∘. Case A-1 uses CFD/CSD/CTD coupling to calculate aerothermoelastic results at this condition. These results are discussed and compared to the results calculated at the same baseline condition by cases A-2 and A-3 which use CFD/CSD coupling and CFD, respectively, in Section 3.6.2. As the quasi-static assumption is adopted in CFD/CSD/CTD coupling, relatively large time steps are used in case A-1. The time step at 0s is set as 0.01s. Then the time step is amplified by the ratio 1.14 at every time step until it reaches 25.95s. After that the time step stands as a constant until the final time 4000s is reached.Table 8
Cases for aerothermoelastic analysis.
Case Calculation method A-1 CFD/CSD/CTD coupling A-2 CFD/CSD coupling A-3 CFD


#### 3.6.2. Aerothermoelastic Analysis

In this subsection, the results of case A-1 by CFD/CSD/CTD coupling method are discussed successively in the order of temperature, stress, deformation, and loads (shear force and bending moment). In the discussion of stress, deformation, and loads, the effect of aerodynamic heating analysis is carried out by comparing the results of CFD/CSD/CTD coupling with the ones of case A-2 by CFD/CSD coupling without aerodynamic heating. Finally, the aerothermoelastic loads of case A-1 are listed and compared with the aeroelastic loads of case A-2 and the aerodynamic loads of case A-3.Figure17 illustrates the transient temperature distributions given by the aerothermoelastic coupling method at various time instants. There are two main directions for the heat conduction, i.e., the tangential direction from the heated leading edge to the cool rear and the normal direction from the surface inward. Moreover, the area with large temperature gradients moves backward over time. When the final time (4000s) is reached, the temperature distributions become steady, and the thermal saturation status with a stagnation temperature of 1117K at leading edge is obtained. Although temperature gradients remain at 4000s because of radiation, the magnitude is much smaller than that at other time instants.Figure 17
Temperature distribution from CFD/CSD/CTD coupling: (a)t=3.6s, (b) t=50s, (c) t=440s, (d) t=4000s.
(a) (b) (c) (d)Figure18 illustrates the stress distributions obtained using the CFD/CSD/CTD coupling method at various time instants. Figure 19 shows the mechanical stress distributions obtained using the CFD/CSD coupling method without consideration of the thermal effects. The results in Figure 18 show that initially the stress maximum is obtained under the thermal shock. Subsequently, the stress concentration area moves backward over time, in agreement with the area with large temperature gradients. In addition, the magnitude of stress decreases, and the concentrated stress in the solid head nearly vanishes when thermal saturation is reached at 4000s. This effect can be explained by the temperature gradients producing a thermal stress that represents the majority of the stress in the solid head. Figure 18 also shows that the stresses on the upper and lower panel of the forebody nearly do not change over time because the stress is produced primarily by the mechanical stress caused by the aerodynamic forces. The similar stress distributions of the upper and lower panel shown in Figure 18(d) and Figure 19 also support this explanation.Figure 18
Stress distribution from CFD/CSD/CTD coupling: (a)t=3.6s, (b) t=50s, (c) t=440s, (d) t=4000s.
(a) (b) (c) (d)Figure 19
Stress distribution from CFD/CSD coupling.The displacements of the stagnation point in theX- and Y-direction are illustrated in Figures 20 and 21, respectively. Figure 20 shows that the magnitude of the X displacement increases monotonically as time goes on. Before 1000s, the X displacement increases rapidly as the structural temperature increases, whereas the displacement increases only minimally after 2000s. This behavior arises because the displacement is contributed predominantly by thermal dilation. Figure 21 shows that the displacement of the stagnation point in the Y-direction first increases and then decreases with time increasing, with an extreme point around 900s. This phenomenon can also be explained by the thermal dilation.Figure 20
Displacement of stagnation in theX-direction from CFD/CSD/CTD coupling.Figure 21
Displacement of stagnation in theY-direction from CFD/CSD/CTD coupling.The shear force and the bending moment of the forebody are shown in Figures22 and 23, respectively. Note that the tendency of the shear force and the bending moment with time is consistent with that of the displacement in the Y-direction because the bend of structure illustrated by the Y displacement magnifies the local inclination initially and then increases the pressure difference between the upper and lower surfaces. The maximum of the shear force and the bending moment appears around 900s rather than at 4000s, when the thermal saturation is reached, reflecting the importance of considering thermal effects via transient aerothermal coupling.Figure 22
Shear force from CFD/CSD/CTD coupling.Figure 23
Bending moment from CFD/CSD/CTD coupling.Figure24 shows a comparison between the deformation obtained using CFD/CSD/CTD coupling when thermal saturation is reached (case A-1, shown in red) and the one obtained using CFD/CSD coupling (case A-2, shown in blue) without considering the thermal effects. The original configuration used by case A-3 is drawn in black to make the comparison clearer. The deformation without thermal effects (blue) is small, and its configuration nearly coincides with the original one, whereas the deformation from aerothermoelastic coupling is significant. This result implies that, in an aerothermoelastic analysis, it is important to consider the thermal effects using CFD/CSD/CTD coupling.Figure 24
Deformation comparison by different methods around leading edge.Table9 compares the loads obtained by CFD/CSD/CTD coupling when the displacement in the Y-direction reaches its maximum with the ones obtained by CFD/CSD coupling without the consideration of the thermal effects. To make the comparison clearer, the loads calculated by CFD only on the original configuration without any deformation are also listed in the table. Relative to the CFD results, the shear force and bending moment obtained using CFD/CSD/CTD coupling increase by 5.7% and 4.1%, respectively, whereas the loads obtained using CFD/CSD coupling vary by only 1.5% and 0.7%, respectively. The results emphasize the importance of the proper introduction of the thermal effects by CFD/CSD/CTD coupling in the load computation.Table 9
Comparison of loads.
Case Calculation method Shear force/N Bending moment/N⋅m A-1 CFD/CSD/CTD coupling 109.45 49.53 A-2 CFD/CSD coupling 105.12 47.24 A-3 CFD 103.54 47.59


#### 3.6.3. Computational Costs

The computational costs of the submodules of CFD/CSD/CTD coupling in one time step are listed in Table10. All calculations are performed on a local PC with one Intel(R) Core(TM)2 Duo processor, 2.13 GHz and 2048 MB memory by serial computing. Table 10 shows that submodules CFD and CSD are the most time-consuming parts of the CFD/CSD/CTD coupling method. To accelerate the CFD/CSD/CTD coupling method, more efforts could be used to reduce the CPU time of the CFD and CSD submodules, such as replacing the LUSGS and LU decomposition method by more advanced generalized minimal residual method (GMRES) or implementing parallel computing.Table 10
Comparison of computational costs between submodules of CFD/CSD/CTD coupling.
Calculation submodule CPU time/s CPU ratio CFD 141.24 53.26% CSD 119.79 45.17% CTD 3.71 1.40% Data transfer 0.17 0.06% Grid deformation 0.3 0.11%


## 3.1. Aerodynamic Coefficients and Pressure Verification

The blunt biconic model in [20] is used to verify the accuracy of the CFD solver in aerodynamic coefficients and pressure calculation. The aerodynamic grid of the blunt biconic model is shown in Figure 3. The parameters used for the current verification study are listed in Table 1.Table 1
The computational parameters for the blunt biconic.
M∞ p ∞ / P a T ∞ / K α / ∘ 6 2200 61.4 5Figure 3
The aerodynamic grid of the blunt biconic.The results of the normal-force coefficient (CN), axial-force coefficient (CA), and pitching-moment coefficient (Cm) obtained by the present CFD solver and the ones from the experiment of [20] are listed and compared in Table 2. The maximum error is 2.04% for the pitching-moment coefficient which means the accuracy of the CFD solver in aerodynamic coefficients is good.Table 2
The comparison of the aerodynamic coefficients of the blunt biconic.
C N C A C m CFD 0.156 0.100 -0.0049 Exp 0.153 0.102 -0.005 Error 1.96% 1.96% 2.04%The pressure distribution (Ps) in the longitudinal symmetric plane obtained by the CFD solver and the one from the experiment of [20] are compared in Figure 4. Figure 4 shows the CFD solutions match well with the experimental results.Figure 4
The comparison of the pressure distribution in the longitudinal symmetric plane.


## 3.2. Aerodynamic Heating Verification

The compressible laminar boundary layer of the plate (1.8m) is used to verify the accuracy of the CFD solver in aerodynamic heating calculation. The aerodynamic grid of the plate is shown in Figure 5. The parameters used for the current verification study are listed in Table 3.Table 3
The computational parameters for the blunt biconic.
M ∞ H ∞ / k m T w / K 5 20 288.15Figure 5
The aerodynamic grid of the plate.The aerodynamic heating distribution by the CFD solvers is compared with the theoretical result of [21] and the solution of Eckert’s reference method of [22] in Figure 6. The figure shows the CFD solver has a better accuracy in aerodynamic heating calculation than Eckert’s reference method compared to the theoretical method.Figure 6
The comparison of aerodynamic heating distribution by different methods.


## 3.3. Heat Conduction Verification

The heated panel model (1m×1m×0.1m) is used to verify the accuracy of the CTD solver in transient temperature calculation. The panel is initially at a uniform temperature of 30K. The upper and right surfaces are isotherm boundary at 100K, while the lower and left surfaces are adiabatic. The grid for heat conduction calculation is shown in Figure 7. The material parameters used for the current verification study are listed in Table 4.Table 4
The material parameters for the heated panel.
λ ρ C p /W/(m⋅K) / k g / m 3 / J / ( k g ⋅ K ) 0.5 3 2Figure 7
The thermal grid of the heated panel.The temperature distributions at1s by analytic method and the present CTD solver are compared in Figure 8. Note that the CTD solver uses a time step of Δt=0.01s and the time matching is carried out for 100 steps. The results of Figure 8 shows that the temperature distribution by the CTD solver agrees well with the analytic solution.Figure 8
The comparison of the temperature distributions at1s.
(a)
By analytic method (b)
By CTD solver


## 3.4. Structural Displacement Verification

The simply supported panel model (1m×1m×0.1m) used to verify the accuracy of the CSD solver in displacement calculation is shown in Figure 9. The panel is initially at a uniform temperature of 300K. The material parameters used for the current verification study are listed in Table 5. The model is uniformly discrete by 120 (60×1×2) second-order hexahedral elements with 20 nodes.Table 5
The material parameters for the simply supported panel.
E ν α ρ /MPa  / K - 1 / k g / m 3 1 0.3 1e-6 2700Figure 9
The simple supported panel model.The mid-plate displacements obtained by the CSD solver based on nonlinear FEM and von Karman theory in [23] are compared in Figure 10. The displacement in Figure 10(a) is calculated under the pressure qa=1Pa, while the displacement in Figure 10(b) is calculated under the pressure qa=1Pa and the heated temperature T=1300K. No matter under which condition, the FEM results agree well with the von Karman theory.Figure 10
The comparison of the mid-plate displacements.
(a)
q a = 1 P a (b)
q a = 1 P a and T=1300K


## 3.5. Aerothermoelastic Panel Verification

To verify the CFD/CSD/CTD coupling method, the multilayer flat panel of [23] is used as the testing model. There are three layers in the panel, namely, the radiation shield, the thermal insulation, and the main structure. The material and thickness of each layer are listed in Table 6, and the thermophysical properties for a temperature of 300K are listed in Table 7. The emissivity of the radiation shield is assumed to be constant at 0.7. In this model, only the stiffness of the main structure is considered, and the stiffnesses of the other two layers are neglected. The panel is supported by immovable, simple supports.Table 6
Composition of the panel.
Layer Thickness/mm Material Radiation shield 7.4 PM-2000 Thermal insulation 10 IMI Main structure 5 Ti-6Al-2Sn-4Zr-2MoTable 7
Properties of the panel at 300 K.
Layer Density Heat capacity Thermal conductivity /kg/m3 /J/(kg⋅K) /W/(m⋅K) Radiation shield 359 465 0.250 Thermal insulation 73.0 729 0.0258 Main structure 4540 463 6.89The panel is located on a5∘ inclined surface of a wedge-shaped forebody. The length of the panel is 1m. The test condition is at altitude 30km and Mach number 8. The initial temperature is assumed to be 300K.The results of the temperature distribution of the upper surface and the displacement at1200s after 60,000 steps of aerothermoelastic coupling iteration are compared with the results of [23] in Figures 11 and 12, respectively. From the comparison, it is evident that the aerothermoelastic coupling method of this study is sufficiently accurate to investigate effects of the introduction of aerodynamic heat coupling.Figure 11
Comparison of temperature distribution of the upper surface (t=1200s).Figure 12
Comparison of displacement.


## 3.6. Aerothermoelastic Forebody Study

After the above verification, the aerothermoelastic load calculation method based on CFD/CSD/CTD coupling is applied to the two-dimensional forebody as shown in Figure13. This model is in the XY plane, and the total length in the X-direction is 1m. The upper surface is parallel to the X-direction, and the angle between lower surface and the X-direction is 10∘. For the sake of thermal protection, the leading edge between the upper surface and the lower surface is rounded by a circle with a 10mm radius. The fluid grid of the forebody is shown in Figure 14 and there are 33320 cells in this grid. Note that cells are clustered near the surface of the forebody.Figure 13
Geometry of the forebody.Figure 14
Aerodynamic grid of the forebody.The finite element model for heat conduction using 8-node, 3-dimensional solid hexahedron elements is shown in Figure15. The total solid element number is 6900 and the total number of nodes is 13662. As shown in Figure 16, there are 1662 cells and 12612 nodes in the structural model formed by 20-node, 3-dimensional solid hexahedron elements. The constraints are imposed at the rear of the structural grid. The material used is a type of high-temperature titanium alloy designated Ta19. The emissivity of the forebody is assumed to be constant at 0.8.Figure 15
Heat conduction grid of the forebody.Figure 16
Structural grid of the forebody.There are two types of loads studied in this work, namely, the shear force and bending moment. The positive direction of the shear force is defined as the positiveY-direction. The bending moment is integrated in the negative Z-direction (leading edge upward), and the integral point is at the rear of the forebody with the same Y-coordinate as the center of circle at leading edge. As the model used is two-dimensional, the loads are integrated by setting the length of the third direction to be 5mm.


### 3.6.1. Cases Description

Table8 lists three cases (A-1, A-2, A-3) considered in this study. The flight condition of the three cases is Mach number 5, altitude 20km, and angle of attack 0∘. Case A-1 uses CFD/CSD/CTD coupling to calculate aerothermoelastic results at this condition. These results are discussed and compared to the results calculated at the same baseline condition by cases A-2 and A-3 which use CFD/CSD coupling and CFD, respectively, in Section 3.6.2. As the quasi-static assumption is adopted in CFD/CSD/CTD coupling, relatively large time steps are used in case A-1. The time step at 0s is set as 0.01s. Then the time step is amplified by the ratio 1.14 at every time step until it reaches 25.95s. After that the time step stands as a constant until the final time 4000s is reached.Table 8
Cases for aerothermoelastic analysis.
Case Calculation method A-1 CFD/CSD/CTD coupling A-2 CFD/CSD coupling A-3 CFD


### 3.6.2. Aerothermoelastic Analysis

In this subsection, the results of case A-1 by CFD/CSD/CTD coupling method are discussed successively in the order of temperature, stress, deformation, and loads (shear force and bending moment). In the discussion of stress, deformation, and loads, the effect of aerodynamic heating analysis is carried out by comparing the results of CFD/CSD/CTD coupling with the ones of case A-2 by CFD/CSD coupling without aerodynamic heating. Finally, the aerothermoelastic loads of case A-1 are listed and compared with the aeroelastic loads of case A-2 and the aerodynamic loads of case A-3.Figure17 illustrates the transient temperature distributions given by the aerothermoelastic coupling method at various time instants. There are two main directions for the heat conduction, i.e., the tangential direction from the heated leading edge to the cool rear and the normal direction from the surface inward. Moreover, the area with large temperature gradients moves backward over time. When the final time (4000s) is reached, the temperature distributions become steady, and the thermal saturation status with a stagnation temperature of 1117K at leading edge is obtained. Although temperature gradients remain at 4000s because of radiation, the magnitude is much smaller than that at other time instants.Figure 17
Temperature distribution from CFD/CSD/CTD coupling: (a)t=3.6s, (b) t=50s, (c) t=440s, (d) t=4000s.
(a) (b) (c) (d)Figure18 illustrates the stress distributions obtained using the CFD/CSD/CTD coupling method at various time instants. Figure 19 shows the mechanical stress distributions obtained using the CFD/CSD coupling method without consideration of the thermal effects. The results in Figure 18 show that initially the stress maximum is obtained under the thermal shock. Subsequently, the stress concentration area moves backward over time, in agreement with the area with large temperature gradients. In addition, the magnitude of stress decreases, and the concentrated stress in the solid head nearly vanishes when thermal saturation is reached at 4000s. This effect can be explained by the temperature gradients producing a thermal stress that represents the majority of the stress in the solid head. Figure 18 also shows that the stresses on the upper and lower panel of the forebody nearly do not change over time because the stress is produced primarily by the mechanical stress caused by the aerodynamic forces. The similar stress distributions of the upper and lower panel shown in Figure 18(d) and Figure 19 also support this explanation.Figure 18
Stress distribution from CFD/CSD/CTD coupling: (a)t=3.6s, (b) t=50s, (c) t=440s, (d) t=4000s.
(a) (b) (c) (d)Figure 19
Stress distribution from CFD/CSD coupling.The displacements of the stagnation point in theX- and Y-direction are illustrated in Figures 20 and 21, respectively. Figure 20 shows that the magnitude of the X displacement increases monotonically as time goes on. Before 1000s, the X displacement increases rapidly as the structural temperature increases, whereas the displacement increases only minimally after 2000s. This behavior arises because the displacement is contributed predominantly by thermal dilation. Figure 21 shows that the displacement of the stagnation point in the Y-direction first increases and then decreases with time increasing, with an extreme point around 900s. This phenomenon can also be explained by the thermal dilation.Figure 20
Displacement of stagnation in theX-direction from CFD/CSD/CTD coupling.Figure 21
Displacement of stagnation in theY-direction from CFD/CSD/CTD coupling.The shear force and the bending moment of the forebody are shown in Figures22 and 23, respectively. Note that the tendency of the shear force and the bending moment with time is consistent with that of the displacement in the Y-direction because the bend of structure illustrated by the Y displacement magnifies the local inclination initially and then increases the pressure difference between the upper and lower surfaces. The maximum of the shear force and the bending moment appears around 900s rather than at 4000s, when the thermal saturation is reached, reflecting the importance of considering thermal effects via transient aerothermal coupling.Figure 22
Shear force from CFD/CSD/CTD coupling.Figure 23
Bending moment from CFD/CSD/CTD coupling.Figure24 shows a comparison between the deformation obtained using CFD/CSD/CTD coupling when thermal saturation is reached (case A-1, shown in red) and the one obtained using CFD/CSD coupling (case A-2, shown in blue) without considering the thermal effects. The original configuration used by case A-3 is drawn in black to make the comparison clearer. The deformation without thermal effects (blue) is small, and its configuration nearly coincides with the original one, whereas the deformation from aerothermoelastic coupling is significant. This result implies that, in an aerothermoelastic analysis, it is important to consider the thermal effects using CFD/CSD/CTD coupling.Figure 24
Deformation comparison by different methods around leading edge.Table9 compares the loads obtained by CFD/CSD/CTD coupling when the displacement in the Y-direction reaches its maximum with the ones obtained by CFD/CSD coupling without the consideration of the thermal effects. To make the comparison clearer, the loads calculated by CFD only on the original configuration without any deformation are also listed in the table. Relative to the CFD results, the shear force and bending moment obtained using CFD/CSD/CTD coupling increase by 5.7% and 4.1%, respectively, whereas the loads obtained using CFD/CSD coupling vary by only 1.5% and 0.7%, respectively. The results emphasize the importance of the proper introduction of the thermal effects by CFD/CSD/CTD coupling in the load computation.Table 9
Comparison of loads.
Case Calculation method Shear force/N Bending moment/N⋅m A-1 CFD/CSD/CTD coupling 109.45 49.53 A-2 CFD/CSD coupling 105.12 47.24 A-3 CFD 103.54 47.59


### 3.6.3. Computational Costs

The computational costs of the submodules of CFD/CSD/CTD coupling in one time step are listed in Table10. All calculations are performed on a local PC with one Intel(R) Core(TM)2 Duo processor, 2.13 GHz and 2048 MB memory by serial computing. Table 10 shows that submodules CFD and CSD are the most time-consuming parts of the CFD/CSD/CTD coupling method. To accelerate the CFD/CSD/CTD coupling method, more efforts could be used to reduce the CPU time of the CFD and CSD submodules, such as replacing the LUSGS and LU decomposition method by more advanced generalized minimal residual method (GMRES) or implementing parallel computing.Table 10
Comparison of computational costs between submodules of CFD/CSD/CTD coupling.
Calculation submodule CPU time/s CPU ratio CFD 141.24 53.26% CSD 119.79 45.17% CTD 3.71 1.40% Data transfer 0.17 0.06% Grid deformation 0.3 0.11%


## 3.6.1. Cases Description

Table8 lists three cases (A-1, A-2, A-3) considered in this study. The flight condition of the three cases is Mach number 5, altitude 20km, and angle of attack 0∘. Case A-1 uses CFD/CSD/CTD coupling to calculate aerothermoelastic results at this condition. These results are discussed and compared to the results calculated at the same baseline condition by cases A-2 and A-3 which use CFD/CSD coupling and CFD, respectively, in Section 3.6.2. As the quasi-static assumption is adopted in CFD/CSD/CTD coupling, relatively large time steps are used in case A-1. The time step at 0s is set as 0.01s. Then the time step is amplified by the ratio 1.14 at every time step until it reaches 25.95s. After that the time step stands as a constant until the final time 4000s is reached.Table 8
Cases for aerothermoelastic analysis.
Case Calculation method A-1 CFD/CSD/CTD coupling A-2 CFD/CSD coupling A-3 CFD


## 3.6.2. Aerothermoelastic Analysis

In this subsection, the results of case A-1 by CFD/CSD/CTD coupling method are discussed successively in the order of temperature, stress, deformation, and loads (shear force and bending moment). In the discussion of stress, deformation, and loads, the effect of aerodynamic heating analysis is carried out by comparing the results of CFD/CSD/CTD coupling with the ones of case A-2 by CFD/CSD coupling without aerodynamic heating. Finally, the aerothermoelastic loads of case A-1 are listed and compared with the aeroelastic loads of case A-2 and the aerodynamic loads of case A-3.Figure17 illustrates the transient temperature distributions given by the aerothermoelastic coupling method at various time instants. There are two main directions for the heat conduction, i.e., the tangential direction from the heated leading edge to the cool rear and the normal direction from the surface inward. Moreover, the area with large temperature gradients moves backward over time. When the final time (4000s) is reached, the temperature distributions become steady, and the thermal saturation status with a stagnation temperature of 1117K at leading edge is obtained. Although temperature gradients remain at 4000s because of radiation, the magnitude is much smaller than that at other time instants.Figure 17
Temperature distribution from CFD/CSD/CTD coupling: (a)t=3.6s, (b) t=50s, (c) t=440s, (d) t=4000s.
(a) (b) (c) (d)Figure18 illustrates the stress distributions obtained using the CFD/CSD/CTD coupling method at various time instants. Figure 19 shows the mechanical stress distributions obtained using the CFD/CSD coupling method without consideration of the thermal effects. The results in Figure 18 show that initially the stress maximum is obtained under the thermal shock. Subsequently, the stress concentration area moves backward over time, in agreement with the area with large temperature gradients. In addition, the magnitude of stress decreases, and the concentrated stress in the solid head nearly vanishes when thermal saturation is reached at 4000s. This effect can be explained by the temperature gradients producing a thermal stress that represents the majority of the stress in the solid head. Figure 18 also shows that the stresses on the upper and lower panel of the forebody nearly do not change over time because the stress is produced primarily by the mechanical stress caused by the aerodynamic forces. The similar stress distributions of the upper and lower panel shown in Figure 18(d) and Figure 19 also support this explanation.Figure 18
Stress distribution from CFD/CSD/CTD coupling: (a)t=3.6s, (b) t=50s, (c) t=440s, (d) t=4000s.
(a) (b) (c) (d)Figure 19
Stress distribution from CFD/CSD coupling.The displacements of the stagnation point in theX- and Y-direction are illustrated in Figures 20 and 21, respectively. Figure 20 shows that the magnitude of the X displacement increases monotonically as time goes on. Before 1000s, the X displacement increases rapidly as the structural temperature increases, whereas the displacement increases only minimally after 2000s. This behavior arises because the displacement is contributed predominantly by thermal dilation. Figure 21 shows that the displacement of the stagnation point in the Y-direction first increases and then decreases with time increasing, with an extreme point around 900s. This phenomenon can also be explained by the thermal dilation.Figure 20
Displacement of stagnation in theX-direction from CFD/CSD/CTD coupling.Figure 21
Displacement of stagnation in theY-direction from CFD/CSD/CTD coupling.The shear force and the bending moment of the forebody are shown in Figures22 and 23, respectively. Note that the tendency of the shear force and the bending moment with time is consistent with that of the displacement in the Y-direction because the bend of structure illustrated by the Y displacement magnifies the local inclination initially and then increases the pressure difference between the upper and lower surfaces. The maximum of the shear force and the bending moment appears around 900s rather than at 4000s, when the thermal saturation is reached, reflecting the importance of considering thermal effects via transient aerothermal coupling.Figure 22
Shear force from CFD/CSD/CTD coupling.Figure 23
Bending moment from CFD/CSD/CTD coupling.Figure24 shows a comparison between the deformation obtained using CFD/CSD/CTD coupling when thermal saturation is reached (case A-1, shown in red) and the one obtained using CFD/CSD coupling (case A-2, shown in blue) without considering the thermal effects. The original configuration used by case A-3 is drawn in black to make the comparison clearer. The deformation without thermal effects (blue) is small, and its configuration nearly coincides with the original one, whereas the deformation from aerothermoelastic coupling is significant. This result implies that, in an aerothermoelastic analysis, it is important to consider the thermal effects using CFD/CSD/CTD coupling.Figure 24
Deformation comparison by different methods around leading edge.Table9 compares the loads obtained by CFD/CSD/CTD coupling when the displacement in the Y-direction reaches its maximum with the ones obtained by CFD/CSD coupling without the consideration of the thermal effects. To make the comparison clearer, the loads calculated by CFD only on the original configuration without any deformation are also listed in the table. Relative to the CFD results, the shear force and bending moment obtained using CFD/CSD/CTD coupling increase by 5.7% and 4.1%, respectively, whereas the loads obtained using CFD/CSD coupling vary by only 1.5% and 0.7%, respectively. The results emphasize the importance of the proper introduction of the thermal effects by CFD/CSD/CTD coupling in the load computation.Table 9
Comparison of loads.
Case Calculation method Shear force/N Bending moment/N⋅m A-1 CFD/CSD/CTD coupling 109.45 49.53 A-2 CFD/CSD coupling 105.12 47.24 A-3 CFD 103.54 47.59


## 3.6.3. Computational Costs

The computational costs of the submodules of CFD/CSD/CTD coupling in one time step are listed in Table10. All calculations are performed on a local PC with one Intel(R) Core(TM)2 Duo processor, 2.13 GHz and 2048 MB memory by serial computing. Table 10 shows that submodules CFD and CSD are the most time-consuming parts of the CFD/CSD/CTD coupling method. To accelerate the CFD/CSD/CTD coupling method, more efforts could be used to reduce the CPU time of the CFD and CSD submodules, such as replacing the LUSGS and LU decomposition method by more advanced generalized minimal residual method (GMRES) or implementing parallel computing.Table 10
Comparison of computational costs between submodules of CFD/CSD/CTD coupling.
Calculation submodule CPU time/s CPU ratio CFD 141.24 53.26% CSD 119.79 45.17% CTD 3.71 1.40% Data transfer 0.17 0.06% Grid deformation 0.3 0.11%


## 4. Conclusions

In this paper, the history of aerothermoelastic solution based on multiphysics coupling is reviewed and the framework developed by Farhat and coworkers is discussed and modified to perform an aerothermoelastic analysis method based on CFD/CSD/CTD coupling for load calculation. The heat conduction equation is corrected by removing the duplicated term. The pseudo-structural method for grid deformation is replaced by a shape-preserving grid deformation strategy based on RBFs to maintain the orthogonality and surface accuracy of the original CFD grid for aerodynamic heating calculation. The fluid and structure equation are solved in a quasi-static manner to accelerate the coupling calculation by using large time step.Then aerothermoelastic analysis method based on CFD/CSD/CTD coupling is performed and verified to investigate the aerothermoelastic effects on the loads of an HSV, such as the shear force and bending moment for the forebody. Using the coupling method, the temperature, stress, displacement, shear force, and bending moment of the forebody in hypersonic flow for a period over4000s are obtained. Based on the numerical results, some conclusions are summarized as follows:(i)
The accuracy of the CFD/CSD/CTD coupling method is verified comprehensively: the key solvers such as CFD, CSD, and CTD are verified individually as well as the coupling method. All these verification works laid a solid foundation for the aerothermoelastic study of the forebody by the CFD/CSD/CTD coupling method.(ii)
The largest increases of the shear force and bending moment via CFD/CSD/CTD coupling are 5.7% and 4.1%, respectively. In comparison, the increases via CFD/CSD coupling are only 1.5% and 0.7%, respectively. This phenomenon emphasizes the importance of adoption of the CFD/CSD/CTD coupling for load computation.(iii)
The stress reaches its maximum concentrating at the leading edge area in the beginning. As the time increases, the magnitude of the stress decreases, and the concentrated stress in the solid head due to temperature gradients nearly vanishes, while the mechanical stress is sustained on the upper and lower panels due to the aerodynamic forces.(iv)
Under the combined action of aerodynamic forces and aerodynamic heating, bending deformation occurs. By comparing the deformation with the one gained by CFD/CSD coupling without consideration of the thermal effects, it is found that the displacement is contributed predominantly by thermal dilation, and the importance of the thermal effects is emphasized.



---

*Source: 6051924-2018-09-18.xml*