Mechanical Stability of Hybrid Corrugated Sandwich Plates under Fluid-Structure-Thermal Coupling for Novel Thermal Protection Systems

Hybrid corrugated sandwich (HCS) plates have become a promising candidate for novel thermal protection systems (TPS) due to their multi-functionality of load bearing and thermal protection. For hypersonic vehicles, the novel TPS that performs some structural functions is a potential method of saving weight, which is significant in reducing expensive design/manufacture cost. Considering the novel TPS exposed to severe thermal and aerodynamic environments, the mechanical stability of the HCS plates under fluid-structure-thermal coupling is crucial for preliminary design of the TPS. In this paper, an innovative layerwise finite element model of the HCS plates is presented, and coupled fluid-structure-thermal analysis is performed with a parameter study. The proposed method is validated to be accurate and efficient against commercial software simulation. Results have shown that the mechanical instability of the HCS plates can be induced by fluid-structure coupling and further accelerated by thermal effect. The influences of geometric parameters on thermal buckling and dynamic stability present opposite tendencies, indicating a tradeoff is required for the TPS design. The present analytical model and numerical results provide design guidance in the practical application of the novel TPS.


Introduction
Hybrid sandwich structures refer to metallic sandwich plates filled with stochastic foam [1]. Compared to traditional lightweight structures, the hybrid sandwich structures exhibit good mechanical performances as well as enhanced multi-functionality with respect to load bearing capacity, energy absorption and thermal protection [2]. Thus, they have demonstrated a wide range of applications in shipbuilding, train manufacturing and construction industry [3][4][5][6]. Among various core topologies of the sandwich structures (e.g., pin-reinforced [7], truss [8], honeycomb [9]), corrugated cores have received major research attention because of their easy fabrication and cost-efficient maintenance [10]. In recent years, there has been growing interest in applying hybrid corrugated sandwich (HCS) plates to the novel thermal protection systems (TPS) of vehicles operating in extreme temperature environments [11][12][13]. Typical representatives are hypersonic vehicles which are exposed to severe thermal and aerodynamic loading. Due to the versatility of the HCS plates, the novel load-bearing TPS offers more weight saving than traditional TPS, leading to performance improvements of the hypersonic vehicles. A significant issue emerged is the mechanical stability of the HCS plates under thermal and aerodynamic loading, while a thorough understanding of the fluid-structure-thermal coupling effects is crucial for the novel TPS design.

Layerwise Theory
Equivalent single theories and layerwise theories are two approaches that are mostly adopted to model the displacement fields for multi-layered laminated and sandwich plates. Compared to the most commonly used equivalent single theories [32][33][34], such as classic laminate theory, first-order or higher-order shear deformation theory, layerwise theories are capable of predicting more accurate static and dynamic behavior of soft-core sandwich structures by a separate description of displacement fields for each layer [35][36][37][38].
In this paper, a layerwise theory is employed to describe displacement fields for the HCS plates with three layers. The piecewise linear in-plane displacement field is assumed through the plate thickness, maintaining displacement continuity at layer interfaces and allowing different slopes in each layer. The displacement field of the kth layer (k = 1, 2, 3) is given as where U (k) , V (k) and U (k+1) , V (k+1) are the in-plane displacements at lower and upper surfaces of the kth layer and w 0 is the out-plane displacement of any point on the mid-plane of the sandwich plates.
where hk is the thickness and  

Layerwise Theory
Equivalent single theories and layerwise theories are two approaches that are mostly adopted to model the displacement fields for multi-layered laminated and sandwich plates. Compared to the most commonly used equivalent single theories [32][33][34], such as classic laminate theory, first-order or higher-order shear deformation theory, layerwise theories are capable of predicting more accurate static and dynamic behavior of soft-core sandwich structures by a separate description of displacement fields for each layer [35][36][37][38].
In this paper, a layerwise theory is employed to describe displacement fields for the HCS plates with three layers. The piecewise linear in-plane displacement field is assumed through the plate thickness, maintaining displacement continuity at layer interfaces and allowing different slopes in each layer. The displacement field of the kth layer (k = 1, 2, 3) is given as where U (k) , V (k) and U (k+1) , V (k+1) are the in-plane displacements at lower and upper surfaces of the kth layer and w 0 is the out-plane displacement of any point on the mid-plane of the sandwich plates. Ψ (k) 1 and Ψ (k) 2 are the linear interpolation functions through the thickness direction and ζ k is the local thickness coordinate, given as Appl. Sci. 2020, 10, 2790 4 of 22 where h k is the thickness and z (l) k and z (u) k are the z-axis coordinates of the lower and upper surfaces.

Strain-Displacement Relations
For a three-dimensional elasticity analysis, the in-plane strain is comprised of linear and nonlinear components. The in-plane strain vector ε (k) and transverse shear strain vector γ (k) of the kth layer (k = 1, 2, 3) in the corrugated hybrid sandwich plates can be given as The nonlinear in-plane strain vector ε (k) nl can also take the form

Stress-Strain Relations
For the HCS plates, the stress-strain relations of the kth layer (k = 1, 2, 3) in the thermal environment with respect to the global coordinate system are given as where Q (k) in and Q (k) tr are the in-plane and transverse shear stiffness matrices. ε (k) 0 is the initial strain vector caused by thermal expansion, given as where α (k) is the thermal expansion coefficient and ∆T(z) is the temperature change through the plate thickness, which could have uniform, linear or parabolic distribution. For the bottom and top face sheets (k = 1, 3), the stiffness matrices and thermal expansion coefficient in the global coordinate system are obtained from those in the material coordinate system with transformation matrices between the global and material coordinates, which are expressed in detail in [38].
For the foam-filled corrugated core (k = 2), the effective stiffness matrices and thermal expansion coefficient are obtained from a homogenization method [19], in which the coupling effects between the corrugated structure and the filled foam are concerned. When the MD of the foam-filled corrugated core coincides with the x-direction of the global coordinate system, the plane stress-reduced and transverse shear stiffness matrices the foam-filled corrugated core are given as where the detailed Q s ij , Q c ij and Q f m ij are given in Equations (A1)-(A3) (Appendix A). Moreover, the components of the effective thermal-stress tensor are given as Therefore, the in-plane thermal expansion coefficient of the foam-filled corrugated core is obtained by The effective density of the foam-filled corrugated core is given as

Finite Element Method
In the present work, an eight-noded C 0 isoparametric element with 9 degrees of freedom per node (4) , w 0 is used to develop a finite element model. For any kth layer (k = 1, 2, 3) of the HCS plates, the relationships between physical displacement vector u (k) , linear in-plane and transverse shear strain vectors ε (k) l and γ (k) l , vector g (k) that is related to nonlinear strain with element nodal displacement vector δ are given as where S (k) is the shape function matrix and B tr , G (k) are the differential operator matrices of the kth layer. These matrices can be divided into matrices H (k) i that are dependent on the thickness coordinate ζ k and matrices N (k) i that are independent of ζ k , where the subscript 'i' donates D, in, tr or G.

Kinetic Energy
The element kinetic energy of the HCS plates is given as where I (k) is the coefficient of the mass matrix for the kth layer and M e is the element consistent mass matrix, given as

Potential Energy
The element potential energy of the hybrid corrugated sandwich plates is given as Substituting Equation (14) into Equation (18), and neglecting terms with third and higher powers in displacement gradients, the total element potential energy can be expressed as U e in and U e tr are the element in-plane and transverse shear strain energy, given as where D (k) in and D (k) tr are the coefficients of the structural stiffness matrices for the kth layer and K e in and K e tr are the element in-plane and transverse shear stiffness matrices, given as V e T is the element potential energy due to the nonlinear strain and the initial stress, given as where S (k) G is the initial stress matrix of the kth layer, given as where E 3×3 denotes a 3×3 identity matrix and ⊗ denotes the Kronecker tensor product. The initial stress vector σ Appl. Sci. 2020, 10, 2790 7 of 22 By using Equation (15), the element potential energy V e T takes the form where D (k) G is the coefficient of the geometric stiffness matrix for the kth layer and K e G is the element geometric stiffness matrix, given as W e T is the work done by the equivalent external force caused by thermal expansion, given as where P

(k)
T is the coefficient of the equivalent nodal force vector for the kth layer and F e T is the element equivalent nodal force vector, given as

Work Done by External Force
Due to their structural application as the TPS of hypersonic vehicles, the hybrid corrugated sandwich plates expose their outer surface in free airflow and withstand transverse aerodynamic pressure. Piston theory is adopted to describe the unsteady aerodynamic load, given as where V, ρ a , and M a are the velocity, density and reference Mach number of free airflow, respectively. Dimensionless dynamic pressure is introduced by defining λ = λ/D, where D = Eh 3 /12 1 − ν 2 . The aerodynamic load evaluated by Piston theory is not only time variable, but also dependent to the gradient of out-of-plane displacement in x-direction.
The element work of the HCS plates is given as where K e a and C e a are the element aerodynamic stiffness and damping matrices and F e a is the element equivalent nodal load caused by aerodynamic pressure, given as where N wi = [0, 0, 0, 0, 0, 0, 0, 0, N i ], ( i = 1 ∼ 8) is the shape function matrix between the out-plane displacement w 0 and the element nodal displacement vector δ.

Governing Equations
By element assembling procedure, the global matrices of mass M, structural stiffness K, geometric stiffness K G , aerodynamic stiffness K a , aerodynamic damping C a and the global nodal force vector F T caused by thermal expansion are obtained according to Equations (17), (21), (26), (31) and (28). The total kinetic energy, potential energy and work done by external force are given in terms of the global nodal displacement vector d, as Hamilton's principle is employed to formulate governing equations, given as Substituting Equation (32) into Equation (33), the governing equation of coupled fluid-structure-thermal analysis for the HCS plates is obtained as The solution of Equation (34) is assumed to be a sum of a time-independent solution due to static deformation and a time-dependent solution due to dynamic deformation as d = d s + d d . Therefore, the static and dynamic equilibrium equations are expressed as

Solution Procedure
In the present work, the aforementioned layerwise theory with the equivalent model of the foam-filled corrugated cores and the finite element formulation based on this theory are adopted to investigated the effects of the thermal and aerodynamic loads on the natural frequencies and structural stabilities of the HCS plates. The solution procedure of coupled fluid-structure-thermal analysis of the HCS plates is illustrated in Figure 2. Firstly, the static stabilities subjected to thermal loads are investigated by buckling analysis and evaluated by the critical temperature change. The critical temperature change limits the maximum temperature allowed in subsequent analysis for steady-state and uniform temperature distributions. Secondly, the dynamic characteristics subjected to thermal loads are investigated by pre-stressed modal analysis in thermal environments and evaluated by the natural frequencies and mode shapes. The geometric stiffnesses induced by the initial thermal stresses are obtained from static analysis. Thirdly, the dynamic stabilities subjected aerodynamic loads are investigated by coupled fluid-structure analysis based on the natural frequencies and modal shapes from modal analysis, which is also called aeroelastic analysis, and evaluated by the critical dynamic pressure. Finally, the dynamic stabilities subjected to both thermal and aerodynamic loads are investigated by coupled fluid-structure-thermal analysis, which is also called aerothermoelastic analysis, and evaluated by the critical dynamic pressure.
frequencies and modal shapes from modal analysis, which is also called aeroelastic analysis, and evaluated by the critical dynamic pressure. Finally, the dynamic stabilities subjected to both thermal and aerodynamic loads are investigated by coupled fluid-structure-thermal analysis, which is also called aerothermoelastic analysis, and evaluated by the critical dynamic pressure. The formulations for each analysis are given as Buckling analysis: Modal analysis: Pre-stressed modal analysis: Coupled fluid-structure analysis: Coupled fluid-structure-thermal analysis: The critical temperature change is defined as 3 10 cr cr TT      and the dimensionless natural frequency is defined as For the coupled fluid-structure-thermal analysis, a modal coordinate transformation is employed to improve the computation efficiency. A new vector of generalized modal coordinates is defined as dd = Φq, where Φ is assembled by normal eigenvectors obtained from Equation (38). Then Equation (40) can be written as The complex eigenvalue is obtained by solving Equation (42). The imaginary part of the eigenvalue is the frequencies and the real part is the damping terms. The dimensionless critical dynamic pressure The formulations for each analysis are given as Buckling analysis: Modal analysis: Pre-stressed modal analysis: Coupled fluid-structure analysis: Coupled fluid-structure-thermal analysis: The critical temperature change is defined as ∆T cr = ∆T cr α × 10 3 and the dimensionless natural frequency is defined as ω = ωa 2 ρh/D. For the coupled fluid-structure-thermal analysis, a modal coordinate transformation is employed to improve the computation efficiency. A new vector of generalized modal coordinates is defined as d d = Φq, where Φ is assembled by normal eigenvectors obtained from Equation (38). Then Equation (40) can be written as The complex eigenvalue is obtained by solving Equation (42). The imaginary part of the eigenvalue is the frequencies and the real part is the damping terms. The dimensionless critical dynamic pressure λ, corresponding to the critical airflow velocity, is determined for the HCS plates when the damping of a rth order mode changes from negative to positive or the coalescence of frequency curves of two adjacent mode occurs.

Numerical Examples and Results
A computer program developed by MATLAB programming language is used for coupled fluid-structure-thermal analysis of the HCS plates. A validation example is presented firstly. Then thermal buckling and pre-stressed modal analysis in thermal environments are conducted. Finally, coupled fluid-structure-thermal analysis is performed with parametric studies.

Validation Study
In order to show the accuracy of the present layerwise finite element model with a homogenization technique, a validation study is performed by comparing the dynamic characteristics of the HCS plates from the present method with those from an MSC.Nastran simulation of a detailed 3D finite element model under different boundary conditions and temperature change.
The HCS plates all made of Ti-6Al-4V with Rohacell 31 foam filled are taken into consideration. The material properties of the foam and base structures that are listed in Tables 1 and 2 refer to the values given by Han [19]. The number of the unit cells with triangular corrugated-core shape is n c = 25 and other dimensions are In the MSC.Nastran software, the corrugated sandwich plate and the filling foam are both simulated by CHEXA solid elements with totally 142,500 elements and 193,431 nodes as shown in Figure 3. This mesh refinement is considered to be sufficient enough that the results of the 3D finite element model can be the exact solutions.  , corresponding to the critical airflow velocity, is determined for the HCS plates when the damping of a rth order mode changes from negative to positive or the coalescence of frequency curves of two adjacent mode occurs.

Numerical examples and results
A computer program developed by MATLAB programming language is used for coupled fluidstructure-thermal analysis of the HCS plates. A validation example is presented firstly. Then thermal buckling and pre-stressed modal analysis in thermal environments are conducted. Finally, coupled fluid-structure-thermal analysis is performed with parametric studies.

Validation study
In order to show the accuracy of the present layerwise finite element model with a homogenization technique, a validation study is performed by comparing the dynamic characteristics of the HCS plates from the present method with those from an MSC.Nastran simulation of a detailed 3D finite element model under different boundary conditions and temperature change.
The HCS plates all made of Ti-6Al-4V with Rohacell 31 foam filled are taken into consideration. The material properties of the foam and base structures that are listed in Table 1 and Table 2 refer to the values given by Han [19]. The number of the unit cells with triangular corrugated-core shape is nc = 25 and other dimensions are a/h = 50, b/a = 1, h1 = h3 = hf = 0.1 h, t = 0.1 h. In the MSC.Nastran software, the corrugated sandwich plate and the filling foam are both simulated by CHEXA solid elements with totally 142,500 elements and 193,431 nodes as shown in Figure 3. This mesh refinement is considered to be sufficient enough that the results of the 3D finite element model can be the exact solutions.     Figure 4 shows the convergence study of the presented layerwise model. The dimensionless fundamental frequency of the HCS plates under SSSS and CCCC boundary conditions with different mesh sizes is illustrated in Figure 4. It is observed that the solution converges with the mesh size of 10×10. Therefore, the following examples are calculated with 10 × 10 mesh size.  Figure 4 shows the convergence study of the presented layerwise model. The dimensionless fundamental frequency of the HCS plates under SSSS and CCCC boundary conditions with different mesh sizes is illustrated in Figure 4. It is observed that the solution converges with the mesh size of 10×10. Therefore, the following examples are calculated with 10 × 10 mesh size.  The first six dimensionless natural frequencies under different boundary conditions obtained from the present method and MSC.Nastran simulation are shown in Figure 5. It is noted that the results calculated by the present method are in good agreement with those by MSC.Nastran, and the relative error is less than 5% for all cases within the frequency range of interest.  The first six dimensionless natural frequencies under different boundary conditions obtained from the present method and MSC.Nastran simulation are shown in Figure 5. It is noted that the results calculated by the present method are in good agreement with those by MSC.Nastran, and the relative error is less than 5% for all cases within the frequency range of interest.  Figure 4 shows the convergence study of the presented layerwise model. The dimensionless fundamental frequency of the HCS plates under SSSS and CCCC boundary conditions with different mesh sizes is illustrated in Figure 4. It is observed that the solution converges with the mesh size of 10×10. Therefore, the following examples are calculated with 10 × 10 mesh size.  The first six dimensionless natural frequencies under different boundary conditions obtained from the present method and MSC.Nastran simulation are shown in Figure 5. It is noted that the results calculated by the present method are in good agreement with those by MSC.Nastran, and the relative error is less than 5% for all cases within the frequency range of interest. Pre-stressed modal analysis in thermal environments of the clamped HCS plates is performed then. The first six dimensionless natural frequencies against temperature change ΔT are shown in Figure 6. It is found that for the first three modes, the results from the present method fit well with those from the MSC.Nastran simulation; meanwhile, for higher-order modes, the relative error of the dimensionless frequencies is a little larger as temperature change increases, but still within an acceptable range. Figure 7 shows the first six mode shapes of the clamped HCS plate with temperature change ΔT = 100. The mode shapes from the present method are as same as those from the MSC.Nastran simulation.
The present method not only has enough calculation accuracy but also has less computational cost compared to commercial software. All the calculations are performed using an Intel Core i7-4790K CPU with a 4GHz main frequency and a 16G memory. It takes about 240 seconds for MSC.Nastran to finish the calculation for one thermal condition; however, the average calculating time for the present method is only 1.5 seconds. Therefore, within an acceptable error range, the present method can save an enormous amount of computing time, much more efficient than MSC.Nastran simulation.  Pre-stressed modal analysis in thermal environments of the clamped HCS plates is performed then. The first six dimensionless natural frequencies against temperature change ∆T are shown in Figure 6. It is found that for the first three modes, the results from the present method fit well with those from the MSC.Nastran simulation; meanwhile, for higher-order modes, the relative error of the dimensionless frequencies is a little larger as temperature change increases, but still within an acceptable range. Figure 7 shows the first six mode shapes of the clamped HCS plate with temperature change ∆T = 100. The mode shapes from the present method are as same as those from the MSC.Nastran simulation.
The present method not only has enough calculation accuracy but also has less computational cost compared to commercial software. All the calculations are performed using an Intel Core i7-4790K CPU with a 4GHz main frequency and a 16G memory. It takes about 240 seconds for MSC.Nastran to finish the calculation for one thermal condition; however, the average calculating time for the present method is only 1.5 seconds. Therefore, within an acceptable error range, the present method can save an enormous amount of computing time, much more efficient than MSC.Nastran simulation. Pre-stressed modal analysis in thermal environments of the clamped HCS plates is performed then. The first six dimensionless natural frequencies against temperature change ΔT are shown in Figure 6. It is found that for the first three modes, the results from the present method fit well with those from the MSC.Nastran simulation; meanwhile, for higher-order modes, the relative error of the dimensionless frequencies is a little larger as temperature change increases, but still within an acceptable range. Figure 7 shows the first six mode shapes of the clamped HCS plate with temperature change ΔT = 100. The mode shapes from the present method are as same as those from the MSC.Nastran simulation.
The present method not only has enough calculation accuracy but also has less computational cost compared to commercial software. All the calculations are performed using an Intel Core i7-4790K CPU with a 4GHz main frequency and a 16G memory. It takes about 240 seconds for MSC.Nastran to finish the calculation for one thermal condition; however, the average calculating time for the present method is only 1.5 seconds. Therefore, within an acceptable error range, the present method can save an enormous amount of computing time, much more efficient than MSC.Nastran simulation.

Buckling amalysis and pre-stressed modal analysis
For the clamped HCS plates subjected to the thermal loads, the static stabilities and dynamic characteristics are evaluated by thermal buckling analysis and pre-stressed modal analysis in thermal environments. The HCS plates are made of Ti-6Al-4V and foam-filled with Rohacell 31. The reference geometric parameters are nc = 25, a/h = 50, b/a = 1, h1 = h3 = hf = 0.1h, t = 0.1h. Figure 8 shows the critical temperature and dimensionless fundamental frequency with ΔT = 100 against geometric parameters change.
The effects of the geometric parameters on the critical temperature change are discussed firstly. The number of the unit cells nc varies from 5 to 50. With more unit cells, the critical temperature change firstly increases and then decreases with a maximum value at about nc = 15. The inclined angle θ takes the values from 35° to 80°. The maximum critical temperature change is achieved at the minimum angle, corresponding to a triangular corrugated shape. The thickness of the face sheets hf and the thickness of the corrugated members t both vary from 0.025h to 0.2h. The critical temperature change increases and then decreases as the facesheet thickness increases, while it decreases as the corrugation thickness increases.
The effects of geometric parameters on the dimensionless fundamental frequency with temperature change ΔT = 100 have a similar tendency as those on the critical temperature change.

Buckling Amalysis and Pre-Stressed Modal Analysis
For the clamped HCS plates subjected to the thermal loads, the static stabilities and dynamic characteristics are evaluated by thermal buckling analysis and pre-stressed modal analysis in thermal environments. The HCS plates are made of Ti-6Al-4V and foam-filled with Rohacell 31. The reference geometric parameters are n c = 25, a/h = 50, b/a = 1, h 1 = h 3 = h f = 0.1h, t = 0.1h. Figure 8 shows the critical temperature and dimensionless fundamental frequency with ∆T = 100 against geometric parameters change.
The effects of the geometric parameters on the critical temperature change are discussed firstly. The number of the unit cells n c varies from 5 to 50. With more unit cells, the critical temperature change firstly increases and then decreases with a maximum value at about n c = 15. The inclined angle θ takes the values from 35 • to 80 • . The maximum critical temperature change is achieved at the minimum angle, corresponding to a triangular corrugated shape. The thickness of the face sheets h f and the thickness of the corrugated members t both vary from 0.025h to 0.2h. The critical temperature change increases and then decreases as the facesheet thickness increases, while it decreases as the corrugation thickness increases.
The effects of geometric parameters on the dimensionless fundamental frequency with temperature change ∆T = 100 have a similar tendency as those on the critical temperature change. Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 21

Coupled fluid-structure-thermal analysis
Coupled fluid-structure-thermal analysis is conducted for the clamped HCS plates at sea level with air density ρa = 1.225kg/m 3 and reference Mach number Ma = 3.0. Figure 9 illustrates the effects of temperature change (ΔT = 100) and aerodynamic pressure on the first two dimensionless natural frequencies.

Coupled Fluid-Structure-Thermal Analysis
Coupled fluid-structure-thermal analysis is conducted for the clamped HCS plates at sea level with air density ρ a = 1.225 kg/m 3 and reference Mach number M a = 3.0. Figure 9 illustrates the effects of temperature change (∆T = 100) and aerodynamic pressure on the first two dimensionless natural frequencies.

Coupled fluid-structure-thermal analysis
Coupled fluid-structure-thermal analysis is conducted for the clamped HCS plates at sea level with air density ρa = 1.225kg/m 3 and reference Mach number Ma = 3.0. Figure 9 illustrates the effects of temperature change (ΔT = 100) and aerodynamic pressure on the first two dimensionless natural frequencies.  For the HCS plates only subjected to thermal loads, the natural frequencies are affected only by the temperature changes. As the temperature increases, the natural frequencies decrease due to the reduced overall stiffnesses caused by the thermal-induced initial stresses. When the sandwich plates are only subjected to aerodynamic loads, the unsteady aerodynamic pressure generates the damping and stiffness that are related to the aerodynamic pressure and with the change of the dynamic pressure, the first two natural frequencies gradually approach to each other. The dynamic pressure corresponding to the coincidence of the first two consecutively varying frequencies is the critical dynamic pressure. For the HCS plate subjected to both thermal and aerodynamic loads, on the one hand, the increasing dynamic pressure causes the first two natural frequencies to coincide with each other and on the other hand, the temperature change causes the decrease in the natural frequencies, resulting in a reduction in the critical dynamic pressure.
For a given temperature field, the dynamic stabilities of the HCS plate under fluid-structure-thermal coupling mainly depend on the imposing aerodynamic pressure, which are measured by the critical dynamic pressure. It is seen from Figure 9 that the critical dynamic pressure reduces with the additional thermal effects.

Parameter Study under Uniform Temperature Field
The effects of geometric and material parameters on critical dynamic pressure of the HCS plates in a uniform temperature field are investigated. The HCS plates are made of Ti-6Al-4V filled with Rohacell 31 foam. The reference geometric parameters are nc = 25, a/h = 50, b/a = 1, Figure 10 illustrates the effects of aspect ratio b/a and airflow direction on the dimensionless critical dynamic pressure. The aspect ratio takes the value of 0.5, 1.0, 1.5 and the airflow direction is parallel or perpendicular to machine direction (MD) of the corrugated plate. It is shown that with an increasing aspect ratio, the critical dynamic pressure decreases regardless of the airflow direction. Meanwhile, the streamwise has a limited influence on the critical dynamic pressure and the critical dynamic pressure is slightly higher when it is parallel to MD as opposed to perpendicular to MD. For the HCS plates only subjected to thermal loads, the natural frequencies are affected only by the temperature changes. As the temperature increases, the natural frequencies decrease due to the reduced overall stiffnesses caused by the thermal-induced initial stresses. When the sandwich plates are only subjected to aerodynamic loads, the unsteady aerodynamic pressure generates the damping and stiffness that are related to the aerodynamic pressure and with the change of the dynamic pressure, the first two natural frequencies gradually approach to each other. The dynamic pressure corresponding to the coincidence of the first two consecutively varying frequencies is the critical dynamic pressure. For the HCS plate subjected to both thermal and aerodynamic loads, on the one hand, the increasing dynamic pressure causes the first two natural frequencies to coincide with each other and on the other hand, the temperature change causes the decrease in the natural frequencies, resulting in a reduction in the critical dynamic pressure.

Effects of aspect ratio and airflow direction
For a given temperature field, the dynamic stabilities of the HCS plate under fluid-structurethermal coupling mainly depend on the imposing aerodynamic pressure, which are measured by the critical dynamic pressure. It is seen from Figure 9 that the critical dynamic pressure reduces with the additional thermal effects.

Parameter study under uniform temperature field
The effects of geometric and material parameters on critical dynamic pressure of the HCS plates in a uniform temperature field are investigated. The HCS plates are made of Ti-6Al-4V filled with Rohacell 31 foam. The reference geometric parameters are nc = 25, a/h = 50, b/a = 1, h1 = h3 = hf = 0.1h, t = 0.1h. Figure 10 illustrates the effects of aspect ratio b/a and airflow direction on the dimensionless critical dynamic pressure. The aspect ratio takes the value of 0.5, 1.0, 1.5 and the airflow direction is parallel or perpendicular to machine direction (MD) of the corrugated plate. It is shown that with an increasing aspect ratio, the critical dynamic pressure decreases regardless of the airflow direction. Meanwhile, the streamwise has a limited influence on the critical dynamic pressure and the critical dynamic pressure is slightly higher when it is parallel to MD as opposed to perpendicular to MD.

Effects of aspect ratio and airflow direction
(a) (b) Figure 10. Effects of (a) aspect ratio and (b) airflow direction on dimensionless critical dynamic pressure against temperature change.
2. Effects of number of unit cell and inclined angle Figure 11 illustrates the effects of number of unit cells nc and inclined angle θ on the dimensionless critical dynamic pressure. The number of unit cells varies from 15 to 45 and the inclined angle varies from 30° to 75°. The critical dynamic pressure increases with the increase in the number of unit cells or the inclined angle. Within a smaller range of temperature change, the discrepancy of the critical dynamic pressure between different numbers of unit cells or different 2. Effects of number of unit cell and inclined angle Figure 11 illustrates the effects of number of unit cells n c and inclined angle θ on the dimensionless critical dynamic pressure. The number of unit cells varies from 15 to 45 and the inclined angle varies from 30 • to 75 • . The critical dynamic pressure increases with the increase in the number of unit cells or the inclined angle. Within a smaller range of temperature change, the discrepancy of the critical dynamic pressure between different numbers of unit cells or different corrugation angles is quite evident. However, for a larger temperature change, the critical dynamic pressure does not change much with these two changing parameters.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 15 of 21 corrugation angles is quite evident. However, for a larger temperature change, the critical dynamic pressure does not change much with these two changing parameters.
(a) (b) Figure 11. Effects of (a) number of unit cell and (b) inclined angle on dimensionless critical dynamic pressure against temperature change.
3. Effects of thickness of face sheets and corrugated members Figure 12 illustrates the effects of the thickness of the face sheets hf and the corrugated members t on the dimensionless critical dynamic pressure with temperature change ΔT = 100. The thickness of face sheets takes the value of 0.05, 0.10, 0.20 and the thickness of corrugated members varies from 0.02 to 0.20. It is clearly seen that increasing either the facesheet or corrugation thickness can both improve the critical dynamic pressure. Moreover, the stability enhancement for the HCS plates is more obvious when thickening the facesheet thickness as opposed to the corrugation thickness.

Effects of Material properties
Material properties of the foam and corrugated sandwich plates are also taken into consideration. Table 3 shows the effects of the material combination of the foam and base structures on the critical temperature change and dimensionless critical dynamic pressure with temperature change ΔT = 100. Ti, Steel and Al denote Ti-6Al-4V, 304 stainless steel and Aluminum 1050 for the base structures, respectively. RC31, RC51, RC71, AC and empty denote the corrugated sandwich plate filled with Rohacell 31, Rohacell 51, Rohacell 71, aluminum foam (Alporas) and no filler, respectively. The material properties are given in Table 1 and Table 2.
It is shown that the base materials of the face sheets and the corrugated members have a much greater influence on the critical dynamic pressure than the foam material. The filling foam increases 3. Effects of thickness of face sheets and corrugated members Figure 12 illustrates the effects of the thickness of the face sheets h f and the corrugated members t on the dimensionless critical dynamic pressure with temperature change ∆T = 100. The thickness of face sheets takes the value of 0.05, 0.10, 0.20 and the thickness of corrugated members varies from 0.02 to 0.20. It is clearly seen that increasing either the facesheet or corrugation thickness can both improve the critical dynamic pressure. Moreover, the stability enhancement for the HCS plates is more obvious when thickening the facesheet thickness as opposed to the corrugation thickness.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 15 of 21 corrugation angles is quite evident. However, for a larger temperature change, the critical dynamic pressure does not change much with these two changing parameters.
(a) (b) Figure 11. Effects of (a) number of unit cell and (b) inclined angle on dimensionless critical dynamic pressure against temperature change.
3. Effects of thickness of face sheets and corrugated members Figure 12 illustrates the effects of the thickness of the face sheets hf and the corrugated members t on the dimensionless critical dynamic pressure with temperature change ΔT = 100. The thickness of face sheets takes the value of 0.05, 0.10, 0.20 and the thickness of corrugated members varies from 0.02 to 0.20. It is clearly seen that increasing either the facesheet or corrugation thickness can both improve the critical dynamic pressure. Moreover, the stability enhancement for the HCS plates is more obvious when thickening the facesheet thickness as opposed to the corrugation thickness.

Effects of Material properties
Material properties of the foam and corrugated sandwich plates are also taken into consideration. Table 3 shows the effects of the material combination of the foam and base structures on the critical temperature change and dimensionless critical dynamic pressure with temperature change ΔT = 100. Ti, Steel and Al denote Ti-6Al-4V, 304 stainless steel and Aluminum 1050 for the base structures, respectively. RC31, RC51, RC71, AC and empty denote the corrugated sandwich plate filled with Rohacell 31, Rohacell 51, Rohacell 71, aluminum foam (Alporas) and no filler, respectively. The material properties are given in Table 1 and Table 2.
It is shown that the base materials of the face sheets and the corrugated members have a much greater influence on the critical dynamic pressure than the foam material. The filling foam increases

Effects of Material properties
Material properties of the foam and corrugated sandwich plates are also taken into consideration. Table 3 shows the effects of the material combination of the foam and base structures on the critical temperature change and dimensionless critical dynamic pressure with temperature change ∆T = 100. Ti, Steel and Al denote Ti-6Al-4V, 304 stainless steel and Aluminum 1050 for the base structures, respectively. RC31, RC51, RC71, AC and empty denote the corrugated sandwich plate filled with Rohacell 31, Rohacell 51, Rohacell 71, aluminum foam (Alporas) and no filler, respectively. The material properties are given in Tables 1 and 2. It is shown that the base materials of the face sheets and the corrugated members have a much greater influence on the critical dynamic pressure than the foam material. The filling foam increases both the critical temperature change and dynamic pressure to a certain degree compared to the unfilled cases. The HCS plates made of Ti-6Al-4V and filled with Rohacell 51 foam have the best performance for both static and dynamic stabilities.

Parameter Study under Uneven Temperature Field
An uneven temperature field is considered for the coupled fluid-structure-thermal analysis of the HCS plates. For the bottom and top face sheets, the temperature fields are treated as uniform due to their relatively thin thickness compared with the core layer. For the foam-filled corrugated core, the temperature field can be a linear or parabolic distribution. Therefore, the temperature field for the bottom face sheet is assumed to be ∆T 1 = 50, the temperature field for the top face sheet is assumed to be ∆T 3 = 300 and the temperature field for the foam-filled corrugated core is given.
In linear distribution: In parabolic distribution: For the HCS plates under this uneven temperature field, the structural dynamic stabilities are mainly determined by the geometric parameters of the foam-filled corrugated core. Therefore, the effects of the inclined angle θ and the corrugation thickness t on the dimensionless critical dynamic pressure under the linear or parabolic temperature field are investigated Figure 13 shows the critical dynamic pressure and the equivalent density against the inclined angle and corrugation thickness in a linear or parabolic temperature field. With the increase in the inclined angle, the equivalent density decreases, and the critical dynamic pressure increases for both linear and parabolic temperature distributions. It seems that increasing the inclined angle is a clever way to improve structural efficiency; however, as shown in Figure 8, the corrugated sandwich plate with a larger inclined angle has a smaller critical temperature change, which is more prone to occur thermal buckling. As for the corrugation thickness, with its increasing, the critical dynamic pressure improves, accompanied by the deterioration of mass and buckling characteristics. Meanwhile, for the given temperature field of bottom and top face sheets, the critical dynamic pressure obtained from the linear temperature distribution is smaller than that from the parabolic distribution. This means that the assumption of linear temperature distribution will lead to more conservative results for the HCS plates under aerodynamic pressure and uneven temperature fields.

Discussion
The validation example in Section 3.1 has proved that the present layerwise finite element model possesses the advantages of both computational accuracy and efficiency compared to MSC.Nastran simulation. This model is thus conducive to the efficient parameter studies of the HCS plates.
The frequency characteristics presented in Figure 9 demonstrate that the unsteady aerodynamic pressure causes two adjacent frequency to gradually coalesce with each other and the coalescence point indicates a critical stable state. The temperature change causes the decrease in natural frequencies which leads to an earlier occurrence of the structural instability.
The geometric and material parameters, including number of unit cells, inclined angle, facesheet thickness, corrugation thickness as well as materials of foam and base structures, have been taken into consideration in the thermal buckling, pre-stressed modal and coupled fluid-structure-thermal analysis of the HCS plates.
In general, the effects of the geometric parameters on the critical temperature change are similar to those on the natural frequency, but quite different from those on the critical dynamic pressure. As shown in Figure 8(a) and Figure 11(a), with the increase in the number of unit cells, the critical temperature change and the fundamental frequency in thermal environments firstly go up and then down, while the critical dynamic pressure keeps growing. As shown in Figure 8(b) and Figure 11(b), with the increase in the inclined angle, the critical temperature change and the fundamental frequency falls, and the critical dynamic rises. The improvement of the critical dynamic pressure due to the increase in inclined angle is obvious under no thermal loads or lower temperature changes. As shown in Figure 8(c) and Figure 12, as the facesheet thickness increases, the critical temperature change and the fundamental frequency reaches a maximum when it is equal to the corrugation thickness, while the critical dynamic presents a great improvement. As shown in Figure 8(d) and Figure 12, as the corrugation thickness increases, the critical temperature change and the fundamental frequency increase, while the critical dynamic decreases. It can be observed that for most cases, parameter changes improve static stabilities while deteriorates dynamic stabilities in coupled fluidstructure-thermal analysis. Due to the different influences of the geometric parameters on the critical temperature change and dynamic pressure, the final structural efficiency of the HCS plates is a tradeoff between static and dynamic stabilities under fluid-structure-thermal coupling. Table 3 displays the effects of material properties of foam and base structures. The critical temperature change is closely related to both the foam and base materials and the critical dynamic Meanwhile, for the given temperature field of bottom and top face sheets, the critical dynamic pressure obtained from the linear temperature distribution is smaller than that from the parabolic distribution. This means that the assumption of linear temperature distribution will lead to more conservative results for the HCS plates under aerodynamic pressure and uneven temperature fields.

Discussion
The validation example in Section 3.1 has proved that the present layerwise finite element model possesses the advantages of both computational accuracy and efficiency compared to MSC.Nastran simulation. This model is thus conducive to the efficient parameter studies of the HCS plates.
The frequency characteristics presented in Figure 9 demonstrate that the unsteady aerodynamic pressure causes two adjacent frequency to gradually coalesce with each other and the coalescence point indicates a critical stable state. The temperature change causes the decrease in natural frequencies which leads to an earlier occurrence of the structural instability.
The geometric and material parameters, including number of unit cells, inclined angle, facesheet thickness, corrugation thickness as well as materials of foam and base structures, have been taken into consideration in the thermal buckling, pre-stressed modal and coupled fluid-structure-thermal analysis of the HCS plates.
In general, the effects of the geometric parameters on the critical temperature change are similar to those on the natural frequency, but quite different from those on the critical dynamic pressure. As shown in Figures 8a and 11a, with the increase in the number of unit cells, the critical temperature change and the fundamental frequency in thermal environments firstly go up and then down, while the critical dynamic pressure keeps growing. As shown in Figures 8b and 11b, with the increase in the inclined angle, the critical temperature change and the fundamental frequency falls, and the critical dynamic rises. The improvement of the critical dynamic pressure due to the increase in inclined angle is obvious under no thermal loads or lower temperature changes. As shown in Figures 8c  and 12, as the facesheet thickness increases, the critical temperature change and the fundamental frequency reaches a maximum when it is equal to the corrugation thickness, while the critical dynamic presents a great improvement. As shown in Figures 8d and 12, as the corrugation thickness increases, the critical temperature change and the fundamental frequency increase, while the critical dynamic decreases. It can be observed that for most cases, parameter changes improve static stabilities while deteriorates dynamic stabilities in coupled fluid-structure-thermal analysis. Due to the different influences of the geometric parameters on the critical temperature change and dynamic pressure, the final structural efficiency of the HCS plates is a trade-off between static and dynamic stabilities under fluid-structure-thermal coupling. Table 3 displays the effects of material properties of foam and base structures. The critical temperature change is closely related to both the foam and base materials and the critical dynamic pressure is mainly determined only by the base materials. Compared to the unfilled sandwich plates, the foam filling improves the structural efficiency to a certain degree. Among the materials adopted in the analysis, the best combination is a Ti-6Al-4V-made corrugated sandwich plate fill with Rohacell 51 foam.
Furthermore, additional investigations about the effects of airflow direction and aspect ratio on the critical dynamic pressure are conducted. Figure 10 indicates that regardless of the aspect ratio, the critical dynamic pressure is always slightly higher when the fiber or machine direction of the HCS plates is parallel to the streamwise than perpendicular to the streamwise.
Finally, uneven temperature fields of linear and parabolic distributions are taken into the consideration. As can be seen in Figure 13, the critical dynamic pressure is smaller in the linear assumption than in the parabolic assumption, indicating that a linear distribution between the top and bottom face sheets brings about a conservative solution.
The future research will focus on the designing, manufacturing and testing of the HCS plates for the novel TPS. Further studies may include 1) optimization design aimed at a balanced structural efficiency; 2) introducing heat transfer into the present fluid-structure-thermal coupling analysis.

Conclusions
In this paper, a C 0 layerwise finite element model with a thermo-elastic homogenization procedure is presented for the coupled fluid-structure-thermal analysis of the HCS plates under thermal and aerodynamic loading. The thermal effect is evaluated by steady temperature change and the unsteady aerodynamic pressure is calculated by the first-order Piston theory. Parameter studies are performed to evaluate the mechanical stability. The main conclusions are summarized as follows: 1.
The present method for the coupled fluid-structure-thermal analysis of the HCS plates possesses both computational accuracy and efficiency.

2.
The dynamic instability of the HCS plates can be induced by fluid-structure coupling and further accelerated by involved thermal effects.

3.
As for frequency characteristics of the HCS plates, unsteady aerodynamic pressure causes two adjacent frequency curves to gradually coalesce with each other, meanwhile, temperature change causes natural frequencies to decrease, resulting in an earlier occurrence of the frequency coalescence. 4.
The influences of design variables on dynamic instability and thermal buckling instability show opposite tendencies, and therefore a tradeoff is required for the TPS design.

5.
A linear assumption of temperature fields through plate thickness brings out conservative solutions compared to parabolic distribution, providing safety margins in the preliminary design. 6. The present analytical model and numerical results for the HCS plates provide design guidance in terms of the practical application of the novel TPS.
Author Contributions: Conceptualization, W.Z., C.Y. and Z.W.; methodology, W.Z. and Z.W.; formal analysis, W.Z.; writing-original draft preparation, W.Z.; writing-review and editing, C.Y. and Z.W. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.