Thermal-Aeroelastic Investigation of a Hypersonic Panel Vibration Based on a Developed MMC Method

Abstract

Hypersonic vehicles or engines usually employ complex thermal protecting shells. This sometimes brings multi-physics difficulties, e.g., thermal-aeroelastic problems like panel flutter etc. This paper aims to propose a novel optimization method versus thermal dynamic influence on panel vibration. A traditional panel structure was modelled and analyzed. After analyzing its dynamic characteristics of panel flutter, thermal effects were also included to propose thermal-aeroelastic analysis results of the present hypersonic panel. Then, a MMC (Moving Morphable Component) method was proposed to suggest dynamic optimization for such panel structures. The proposed method can provide arbitrary frequency control result in order to suggest a newly generated panel structure. Based on the optimal structures, dynamic analysis was presented again to verify the effectiveness of the optimization method. So aero-thermo-dynamic characteristics of the optimal panel structures could be investigated. It can be seen that the computational results presented significantly improved panel flutter results. The proposed dynamic optimization method can be employed for the design of panel structures versus high combustion temperatures or hypersonic aerodynamics.

1. Introduction

In recent years, hypersonic vehicles have attracted the attention of researchers for their extremely fast flight speeds, but hypersonic vehicles also face many problems during flight, such as the panel flutter problem. This problem may be more pronounced in the thermal environment generated by hypersonic flights. Panel flutter is a self-excited vibration behavior of the panel structure under the combined effect of aerodynamic, inertial and elastic forces. When panel flutter occurs, it can cause severe vibration of the panel leading to failure. Therefore, some significant works have presented investigations on the aeroelastic optimization of panels.

Wang et al. [1] proposed an energy method based on Galerkin method and a two-degree-of-freedom reduced-order model that can be used in the calculation of panel flutter. The aerodynamic forces are applied to structure by first-order piston theory. Zhang et al. [2] investigated the flutter characteristics of a nonlinear spring-supported composite Panel. The aerodynamic forces are calculated by the third piston theory, and the structural equations of motion are obtained by von Karmen’s nonlinear theory. Since the piston theory assumes that the perturbation at each point of the airfoil propagates in the direction of that point and ignores the interaction of the airfoil points, some works were presented to use other aerodynamic theories to calculate the panel flutter. Serafim et al. [3] used unsteady potential flow aerodynamics to calculate flutter characteristics of panel at subsonic, transonic and supersonic conditions.

As the flight velocity of the vehicle is increasing, the panel will be subjected to a high temperature environment; many studies [4,5] were presented on the panel flutter problem in the thermal environment. Ye et al. [6] studied the effect of temperature-dependent material properties on the buckling and flutter characteristics of a heated panel. The results show that the temperature-dependent material parameters have a very significant effect on the buckling and flutter characteristics of the panel. Chen et al. [7] investigated the flutter characteristics of angle tow composite curved panels under aerodynamic and thermal load. Similar conclusions were also reached as above. Abdollahi et al. [8] has given an investigation of aero-thermo-elastic flutter of functionally graded porous skew panel. The governing equations are obtained by Hamilton’s principle combined with the first-order shear theory and the first-order piston theory. The effects of constant, linear and nonlinear temperature fields on the flutter characteristics of structures are considered. Javadi et al. [9] investigated the aero-thermoelastic characteristics of porous 2D curved panels. Additionally, the effects of porosity distribution, yawed flow angle, curvature ratio and Mach number on the aeroelastic properties were investigated.

In order to improve the chattering characteristics of the panel, some studies have optimized and improved the panel; one of the optimization methods is to optimize the structural parameters of the panel. Bochkarev et al. [10] investigated the aeroelastic stability of shallow cylindrical shells stiffened with stringers; the aerodynamic forces were obtained by quasi-steady first-order piston theory, and the optimal location of stringers was achieved by varying the distance between stringers. Fazilati et al. [11] investigated the flutter characteristics of panels with fiber reinforcement and optimized the flutter of the panels by varying the parameters such as the angle of the fiber layer. Another approach to optimizing the panel is to use topological optimization methods such as the implicit optimization method SIMP. Stanford et al. [12] optimized the topological configuration of the reinforcing ribs with the SIMP method using the buckling and flutter characteristics of the panel as the optimization objective. The stiffeners are divided into finite elements and the topological configuration of the stiffeners is changed by varying the element density from 0 to 1. Stanford et al. [13] also used the SIMP method to optimize panels subjected to aerodynamic, elastic, inertial and thermal loads. Optimization is for the panel itself and not the ribs.

However, the abovementioned optimization studies may need some developments. For the optimization of reinforcement using parameter optimization [10,11], the optimization design domain is small, and the improvement of structural performance is limited. The boundaries of the structure obtained using implicit topology optimization [12,13] are not clear, which is not conducive to processing and manufacture.

Therefore, we used an explicit topology optimization method-MMC (Moving Morphable Component) [14] to optimize the aeroelastic properties of the panel. In this paper, dynamics and aeroelastic characteristics of a standard panel are firstly presented where thermal effects are included. Then, a developed MMC method with the optimization objective versus the difference between the first and the second natural frequencies is presented to control the critical flutter velocity of the panel. Meanwhile, for engineering applications, the constraints of the optimization process are proposed to obtain applicable panel optimization. Finally, dynamics and aeroelastic characteristics of three kinds of optimized structures are presented and compared.

2. Panel Structure Introduction

2.1. Geometry and Finite Element Modelling

A common reinforced panel structure with fixed support around its edge is shown in Figure 1 as a standard panel. The length and width of this panel are both 500 mm, the thickness of the ribs is 4 mm, and the thickness of the panel is 1 mm (see

Table 1. Panel geometry.

Parameter	Value
Length (m)	0.5
width (m)	0.5
Thickness (m)	0.004 (rib)
	0.001 (panel)

). The length and width of the panel are in the x and y directions, respectively. A temperature field (see Figure 2), varying linearly along the flow direction, was applied to the panel structure to include the thermal effect (i.e., the present temperature field is a hypothetical temperature field which was employed for method investigation, not a physical temperature field).

Material parameters [15] used in the panel structure are shown in

Table 2. Material parameters.

Temperature (K)	Elastic (Gpa)	Density (kg/m3)	Poisson’s Ratio	Coefficient of Thermal Expansion
293	72	2700	0.33	$2.31\times {10}^{-5}$
373	71			$2.46\times {10}^{-5}$
423	68			$2.56\times {10}^{-5}$
478	63			$2.69\times {10}^{-5}$

. Thermal effects are mainly found on elastic and thermal expansion.

Finite element modelling (FEM) was used to find dynamic characteristics (i.e., natural frequency and mode) (see Figure 3) which can also be employed for aeroelasticity analysis of the above structure. Kirchhoff thin plate elements were used to build the finite element model. In total, 7056 elements and 7225 nodes were modelled.

2.2. Dynamic Characteristics

The finite element governing equations for the motion of the panel can be written as Equation (1) [10].

$$\mathit{M}\ddot{\mathit{x}}+\mathit{C}\dot{\mathit{x}}+\mathit{K}\mathit{x}=\mathit{f}$$

(1)

Here, $\mathit{M}$,$\mathit{C}$,$\mathit{K}$ are the mass matrix, damping matrix and stiffness matric, respectively; $\mathit{x}$ denotes the node displacement vector; and $\mathit{f}$ represents the node load vector, which is the aerodynamic load calculated by supersonic lifting surface method in this paper. $\mathit{M}$,$\mathit{C}$,$\mathit{K}$ can be calculated by Equation (2).

$$\begin{gathered} \mathit{M}={\int }_{\mathrm{V}}\rho {\mathit{N}}^{\mathrm{T}}\mathit{N}\mathrm{d}\mathrm{V} \\ \mathit{C}={\int }_{\mathrm{V}}\mu {\mathit{N}}^{\mathrm{T}}\mathit{N}\mathrm{d}\mathrm{V} \\ \mathit{K}={\int }_{\mathrm{V}}{\mathit{B}}^{\mathrm{T}}\mathit{D}\mathit{B}\mathrm{d}\mathrm{V} \end{gathered}$$

(2)

Here, $\rho$ is material density; $\mu$ is the viscosity factor; $\mathit{N}$ is the shape function [16]; $\mathit{B}$ is the strain matrix; and $\mathit{D}$ is the elasticity matrix.

With the above finite element model, it is possible to obtain structural eigenvalues and eigenvectors by solving the generalized eigenvalue problem, which is written as:

$$\begin{gathered} \mathit{K}\mathit{\Phi }=\mathit{M}\mathit{\Phi }\mathit{\Lambda } \\ \mathit{\Lambda }=\left[\begin{array}{ccc}{\lambda }_1& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & {\lambda }_N\end{array}\right] \\ \mathit{\Phi }=\left[\begin{array}{ccc}{\mathit{\phi }}_1& \cdots & {\mathit{\phi }}_N\end{array}\right] \end{gathered}$$

(3)

Here, $\mathit{\Phi }$ represents the matrix of eigenvectors; $\mathit{\Lambda }$ represents the matrix of eigenvalues; ${\lambda }_j$ and ${\mathit{\phi }}_j$ $(j=\mathrm{1,2},\cdots ,N)$ are the $j\mathrm{t}\mathrm{h}$ order eigenvalue and eigenvector; and $N$ is the number of degrees of freedom of the structure.

Natural frequencies and modes of the panel are shown in

Table 3. Natural frequency of the standard panel (Hz).

Modes	Normal Modes	Thermal Modes
1	79.96	78.13
2	156.51	153.02
3	156.51	153.08
4	224.70	220.12
5	281.57	275.46
6	288.84	282.73

and Figure 4, respectively.

To describe the relative difference between the first two orders of natural frequency, we have defined a frequency ratio $R$ as Equation (4). The smaller the parameter R, the greater the difference in frequency between the first two orders, and the better the aeroelastic characteristics of the structure.

$$R=\frac{1}{\frac{{\omega }_2-{\omega }_1}{{\omega }_1}}=\frac{{\omega }_1}{{\omega }_2-{\omega }_1}$$

(4)

Here, $\omega$ is the is the circular frequency of the structure. If $R=1$, the second natural frequency is exactly twice the first natural frequency. If $R>1$, the second natural frequency is less than twice the first natural frequency. If $R<1$, the second natural frequency is more than twice the first natural frequency.

$R$ of standard structure without considering temperature effects is 1.0445 implies that the second natural frequency is less than twice the first natural frequency.

Next, we will consider thermal effects on the dynamic characteristics of the structure. Thermal effects on the panel can be found to influence the stiffness of the structure; the stiffness matrix of the panel under the thermal load can be written as Equation (5) [17]:

$$\mathit{K}={\mathit{K}}_{\mathit{T}}+{\mathit{K}}_{\mathsf{\sigma }}+{\mathit{K}}_{\mathit{N}\mathit{L}}$$

(5)

Here, ${\mathit{K}}_{\mathit{T}}$ denotes the stiffness matrix due to the change in the modulus of elasticity; ${\mathit{K}}_{\mathit{\sigma }}$ denotes the stiffness matrix due to thermal stress; and ${\mathit{K}}_{\mathit{N}\mathit{L}}$ denotes the nonlinear stiffness matrix due to large structural deformation.

Large deformation was not included here, so the above equation is simplified as:

$$\begin{gathered} \mathit{K}={\mathit{K}}_{\mathit{T}}+{\mathit{K}}_{\mathit{\sigma }} \\ {\mathit{K}}_{\mathit{T}}={\int }_{\mathrm{V}}\mathit{B}^{\mathrm{T}}{\mathit{D}}_{\mathit{T}}\mathit{B}\mathrm{d}\mathrm{V} \\ {\mathit{K}}_{\mathit{\sigma }}={\int }_{\mathrm{V}}\mathit{G}^{\mathrm{T}}\mathit{S}\mathit{G}\mathrm{d}\mathrm{V} \end{gathered}$$

(6)

Here, ${\mathit{D}}_{\mathit{T}}$ denotes elasticity matrix (i.e., a function of temperature); $\mathit{G}$ denotes the derivative of the shape function matrix; and $\mathit{S}$ denotes the thermal stress matrix due to thermal load.

By calculating the stiffness matrix of the structure with temperature load and solving the generalized eigenvalue problem as Equation (3), frequencies and modes of the panel structure subjected to temperature field can be obtained, as shown in Table 3 and Figure 5. Natural frequencies of the panel decrease as the temperature increases, the effect of temperature becomes more significant as frequency increasement. After considering the temperature effect, the first natural frequency decreases by less than 2 $\mathrm{H}\mathrm{z}$, while the 6th natural frequency decreases by more than 6 $\mathrm{H}\mathrm{z}$. The modes of the panel are not varied obviously.

2.3. Aero-Elastic Analysis

Before aeroelastic analysis of the panel structure is performed, aerodynamic loads on the panel should be evaluated. In this paper, aerodynamic load was computed with a supersonic lifting surface method [18,19]. The equation for the supersonic lifting surface method is specified as follows:

$$w_k=\frac{1}{4\pi {\rho }_{\mathrm{f}}{V}^2}\sum _j\Delta p_j{\iint }_{{\mathrm{S}}_j}{K}_{kj}\mathrm{d}\mathrm{S}\left(k=1,2,\cdots ,n;j=1,2,\cdots ,n\right)$$

(7)

Here, $w_k$ is the downwash speed at the control point of box $k$; ${\rho }_{\mathrm{f}}$ is the density of air; $V$ is the flight velocity; $\Delta p_j$ is the pressure applied on box $j$; ${\mathrm{S}}_j$ means the region where the inverted Mach cone from the $j$th downwash point intercepts the airfoil; $K_{kj}$ is the kernel function for aerodynamic calculations; and $n$ is the number of lifting surface boxes.

The pressure on the $j\mathrm{t}\mathrm{h}$ box can be written as follows:

$$\Delta p_j=\Delta c_{p_j}\left(\frac{1}{2}{\rho }_{\mathrm{f}}{V}^2\right)$$

(8)

Here $,\Delta p_j$ is the pressure applied on the $j\mathrm{t}\mathrm{h}$ panel, and $\Delta c_{p_j}$ is the pressure coefficient of the $j\mathrm{t}\mathrm{h}$ lifting surface box.

Then

$$w_k=\sum _j\frac{1}{8\pi }\Delta c_{p_j}{\iint }_{{\mathrm{S}}_j}{K}_{ij}\mathrm{d}\mathrm{S}$$

(9)

Let

$$D_{kj}=\frac{1}{8\pi }{\iint }_{{\mathrm{S}}_j}{K}_{kj}d\mathrm{S}$$

(10)

Then

$$w_k=\sum _j{D}_{kj}\Delta c_{p_j}$$

(11)

The above equation is written in the form of a matrix:

$$\begin{gathered} \mathit{w}=\mathit{D}∆{\mathit{c}}_{\mathit{p}} \\ \mathit{w}=\frac{1}{V}{\left[w_1,w_2,\cdots w\right]}^{\mathrm{T}} \\ ∆{\mathbf{c}}_{\mathit{p}}={\left[\Delta c_{p_1},\Delta c_{p_2},\cdots \Delta c_{p_n}\right]}^{\mathrm{T}} \\ \mathit{D}=\left[\begin{array}{ccc}{D}_{11}& \cdots & D_{1n}\\ ⋮& \ddots & ⋮\\ D_{n1}& \cdots & D_{nn}\end{array}\right] \end{gathered}$$

(12)

Here, $\mathit{w}$ is the downwash speed vector, $\mathit{D}$ is the matrix of aerodynamic influence coefficient, and $∆{\mathbf{c}}_{\mathit{p}}$ is the vector of pressure coefficient.

The downwash velocity and vibration displacement at the grid control points satisfy the following equations:

$$\begin{gathered} \mathit{w}=\left(\frac{\partial \mathit{z}}{\partial x}+\frac{1}{V}\frac{k}{b}\mathit{z}\right) \\ \mathit{z}={\left[z_1,z_2,\cdots ,z_n\right]}^{\mathrm{T}} \end{gathered}$$

(13)

Here, $k$=$\frac{\omega b}{V}$ is the reduced frequency, $b$ is the reference length, and $z_j$ $\left(j=1,2,\cdots ,n\right)$ is z-direction displacement at the control point of lifting surface box $j$.

Thus, the pressure vector can be expressed as:

$$\mathit{\Delta }\mathit{p}=\frac{1}{2}{\rho }_f{V}^2{\mathit{D}}^{-1}\left(\frac{\partial \mathit{z}}{\partial x}+\frac{1}{V}\frac{k}{b}\mathit{z}\right)$$

(14)

Here, $\mathbf{\Delta }\mathit{p}$ is a 1 × $n$-dimensional pressure vector.

After obtaining the pressure loads on the lifting surfaces, the pressure loads acting on the structure can be obtained by Infinite Plate Spline (IPS) [20].

Substituting the aerodynamic forces into the structural dynamics governing equations and applying a modal coordinate transformation to the equations, the governing equations for the flutter problem can be obtained as Equation (15) [21]:

$$\widehat{\mathit{M}}\ddot{\mathit{q}}\left(\mathit{t}\right)+\widehat{\mathit{C}}\dot{\mathit{q}}\left(\mathit{t}\right)+\widehat{\mathit{K}}\mathit{q}\left(\mathit{t}\right)=\frac{1}{2}{\rho }_f{V}^2\mathit{Q}\left(k,Ma\right)\mathit{q}\left(\mathit{t}\right)$$

(15)

Here, $V$ denotes air velocity; $\widehat{\mathit{M}}$, $\widehat{\mathit{C}}$ and $\widehat{\mathit{K}}$ are the modal mass matrix, modal damping matrix and modal stiffness matrix, respectively; $\mathit{Q}\left(k,Ma\right)$ is modal aerodynamic matrix.

The modal mass matrix, modal damping matrix, modal stiffness matrix and modal aerodynamic matrix can be computed by the following Equation (16):

$$\begin{gathered} \widehat{\mathit{M}}={\mathit{\Phi }}_{\mathit{s}\mathit{e}}^{\mathrm{T}}\mathit{M}{\mathit{\Phi }}_{\mathit{s}\mathit{e}} \\ \widehat{\mathit{C}}={\mathit{\Phi }}_{\mathit{s}\mathit{e}}^{\mathrm{T}}\mathit{C}{\mathit{\Phi }}_{\mathit{s}\mathit{e}} \\ \widehat{\mathit{K}}={\mathit{\Phi }}_{\mathit{s}\mathit{e}}^{\mathrm{T}}\mathit{K}{\mathit{\Phi }}_{\mathit{s}\mathit{e}} \\ \mathit{Q}\left(k,Ma\right)={\mathit{\Phi }}_{\mathit{s}\mathit{e}}^{\mathrm{T}}\left({\mathit{G}}_{\mathit{a}\mathit{s}}^{\mathrm{T}}-{\mathit{G}}_{\mathit{b}\mathit{s}}^{\mathrm{T}}\mathit{\Gamma }\right)\mathit{S}{\mathit{D}}^{-1}\left[{\mathit{G}}_{\mathit{b}\mathit{s}}+\frac{ik}{b}\left({\mathit{G}}_{\mathit{a}\mathit{s}}+\mathit{\Gamma }{\mathit{G}}_{\mathit{b}\mathit{s}}\right)\right]{\mathit{\Phi }}_{\mathit{s}\mathit{e}} \\ \mathit{\Gamma }=\left[\begin{array}{ccc}\Delta x_1& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & \Delta x_n\end{array}\right] \\ \mathit{S}=\left[\begin{array}{ccc}\Delta s_1& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & \Delta s_n\end{array}\right] \end{gathered}$$

(16)

Here, ${\mathit{\Phi }}_{\mathbf{s}\mathbf{e}}$ is the matrix of the selected modal vectors; ${\mathit{G}}_{\mathit{a}\mathit{s}}$ and ${\mathit{G}}_{\mathit{b}\mathit{s}}$ are the spline matrices representing the relationship between the normal displacement of the grid points of the lifting surfaces, the slopes along the air flow direction and the displacements of the structural nodes, respectively; $\Delta x_j$ and $\Delta s_j$ $(j=1,2,3\cdots n)$ are the chord length and the area of the $j\mathrm{t}\mathrm{h}$ panel.

Equation (15) can be solved by the P-K method [22,23]; the above equation can be rewritten as follows:

$$\begin{gathered} \mathit{Q}\left(k,Ma\right)=i{\mathit{Q}}_{\mathit{I}\mathit{m}}\left(k,Ma\right)+{\mathit{Q}}_{\mathit{R}\mathit{e}}\left(k,Ma\right) \\ p=\left(\gamma +i\right)k \\ \left[{\left(\frac{V}{b}\right)}^2\widehat{\mathit{M}}{p}^2+\frac{V}{b}\widehat{\mathit{C}}p+\widehat{\mathit{K}}-i\frac{1}{2}{\rho }_f{V}^2{\mathit{Q}}_{\mathit{I}\mathit{m}}\left(k,Ma\right)-\frac{1}{2}{\rho }_f{V}^2{\mathit{Q}}_{\mathit{R}\mathit{e}}\left(k,Ma\right)\right]\overline{\mathit{q}}=0 \end{gathered}$$

(17)

Here, $\overline{\mathit{q}}$ is the modal displacement amplitude vector.

Equation (17) can be transformed into the following canonical form:

$$\mathit{A}\left(k,Ma\right)\psi =p\mathit{B}\left(k,Ma\right)\psi$$

(18)

where

$$\begin{gathered} \mathit{A}\left(k,Ma\right)=\left[\begin{array}{cc}-\widehat{\mathit{K}}+\frac{1}{2}{\rho }_f{V}^2{\mathit{Q}}_{\mathit{R}\mathit{e}}\left(k,Ma\right)& \mathbf{0}\\ \mathbf{0}& {\left(\frac{V}{b}\right)}^2\widehat{\mathit{M}}\end{array}\right] \\ \mathit{B}\left(k,Ma\right)=\left[\begin{array}{cc}\frac{V}{b}\widehat{\mathit{C}}-\frac{1}{2}{\rho }_f{V}^2\frac{{\mathit{Q}}_{\mathit{I}\mathit{m}}\left(k,Ma\right)}{k}& {\left(\frac{V}{b}\right)}^2\widehat{\mathit{M}}\\ {\left(\frac{V}{b}\right)}^2\widehat{\mathit{M}}& \mathbf{0}\end{array}\right] \end{gathered}$$

(19)

The eigenvalues $p$ can be obtained by solving the above eigenvalue problem, which leads to the computation of $g$ and $f$. If the real part of eigenvalue $p$ is positive, the system is unstable, and if the real part of eigenvalue $p$ is negative, the system is stable.

$$\begin{gathered} g\approx 2\gamma =2\frac{Re\left(p\right)}{k} \\ f=\frac{kV}{2\pi b} \end{gathered}$$

(20)

$p$ can be calculated in the desired velocity range, then plot V–g, and the velocity at the point where $\gamma$ is 0 is the flutter velocity.

The lifting surface mesh of the standard panel has 20 lifting surface panels in both $i$ and $j$ directions. The flutter analysis was taken by the commercial software NASTRAN 2012. The V-g and V-f plots of panel flutter with or without thermal effect are shown in

Table 4. Critical flutter velocity of the standard panel.

Flutter	Velocity (m/s)
Natural modes	1050
Thermal modes	640

and Figure 6. The critical flutter velocity of the panel without considering thermal effect is 1050 m/s, and this value becomes 640 m/s when thermal effect is included. It can be found that the critical flutter velocity of the structure becomes lower and the aeroelastic stability of the structure becomes worse due to thermal influence.

3. Optimization via Both Thermal-Aeroelastic and Manufacture Consideration

3.1. Optimization Method

In order to improve the aeroelastic stability of the panel structure, a MMC (Moving Morphable Component) method was proposed to optimize the design of the panel structure. The MMC method was initially proposed by Guo et al. [14]. The basic idea of this method is to use the deformation and movement of components to achieve topology optimization of the structure. Areas covered with components are filled with material and areas without components are not filled with material. Compared to traditional implicit topology optimization methods such as the SIMP method, the MMC method has the advantage that the optimized structure is expressed explicitly. It can avoid the intermediate density element problem that occurs with implicit methods. The final optimization result can be explicitly expressed as a function and can be directly imported into some commercial CAD software.

The topological optimization equation based on MMC can be specified as Equation (21) [24].

$$\begin{gathered} \mathrm{Find}\mathit{d}={\left({\mathit{d}}_1,\dots ,{\mathit{d}}_{nc}\right)}^{\mathrm{T}} \\ \mathrm{Minimize}I=I(\mathit{d}) \\ \mathrm{s}.\mathrm{t}. \\ {g}_j(\mathit{d})\le 0,j=1,\dots ,l \\ \mathit{d}\subset {\mathcal{U}}_d \end{gathered}$$

(21)

Here, $nc$ represents the total number of components; $\mathit{d}={\left({\mathit{d}}_1,\dots ,{\mathit{d}}_{nc}\right)}^{\mathrm{T}}$ represents the vector of the design variables; $I$ represents the objective equation; $g_j(\mathit{d})\le 0,j=1,\dots ,l$ represents the constraint equations; and ${\mathcal{U}}_{\mathit{d}}$ is the admissible sets that $\mathit{d}={\left({\mathit{d}}_1,\dots ,{\mathit{d}}_{nc}\right)}^{\mathrm{T}}$ belongs to.

In this paper, the components are treated as reinforcement, and the area not covered by the components is the panel (see Figure 7). We would like to use the variation in the reinforcement to adjust the frequency characteristics of the panel structure. Mostly, panel flutter occurs mostly due to the coupling of the first and the second natural modes. Therefore, in order to increase the flutter velocity of the panel, we increase the difference between the first and second natural frequencies of the panel structure.

The proposed topology optimization in Case A is to make the variable $R$ as small as possible, the number of components was set to eight, the problem formulation can be specified as Equation (22).

$$\begin{gathered} \mathrm{Find}\mathit{d}={\left({\mathit{d}}_1,\dots ,{\mathit{d}}_8\right)}^{\mathrm{T}} \\ \mathrm{Minimize}R \\ \mathrm{s}.\mathrm{t}. \\ \mathit{K}{\mathit{\phi }}_j={\omega }_j^2\mathit{M}{\mathit{\phi }}_j,j=1,\dots n,\dots ,N \\ \mathit{d}\subset {\mathcal{U}}_d \\ {V}_t\le {\overline{V}}_t \end{gathered}$$

(22)

Here, $\mathit{K}$ and $\mathit{M}$ are the stiffness matrix and mass matrix, respectively, ${\overline{V}}_t$ is the upper bound of the available materials volume.

In Case B, in order to make the optimized structural reinforcement form simple, we add the compliance as one of the constraints in this optimization example. The optimization equation can be written as Equation (23).

$$\begin{gathered} \mathrm{Find}\mathit{d}={\left({\mathit{d}}_1,\dots ,{\mathit{d}}_8\right)}^{\mathrm{T}} \\ \mathrm{Minimize}R \\ \mathrm{s}.\mathrm{t}. \\ \mathit{K}{\mathit{\phi }}_j={\omega }_j^2\mathit{M}{\mathit{\phi }}_j,j=1,\dots n,\dots ,N \\ \mathit{d}\subset {\mathcal{U}}_d \\ {V}_t\le {\overline{V}}_t \\ C\le \overline{C} \end{gathered}$$

(23)

Here, $C$ is the compliance of the structure, and $\overline{C}$ is the upper bound of compliance.

In order to make the reinforcement form closer to the standard structure and further simplify the structure, we change the number of components to four in Case C, but the objective function and constraints are not changed. The optimization column of this problem is given by Equation (24).

$$\begin{gathered} \mathrm{Find}\mathit{d}={\left({\mathit{d}}_1,\dots ,{\mathit{d}}_4\right)}^{\mathrm{T}} \\ \mathrm{Minimize}R \\ \mathrm{s}.\mathrm{t}. \\ \mathit{K}{\mathit{\phi }}_j={\omega }_j^2\mathit{M}{\mathit{\phi }}_j,j=1,\dots n,\dots ,N \\ \mathit{d}\subset {\mathcal{U}}_d \\ {V}_t\le {\overline{V}}_t \\ C\le \overline{C} \end{gathered}$$

(24)

The sensitivity of the objective function to the design variables is derived as follows:

$$\frac{\partial R}{\partial \alpha }=\frac{\partial \frac{{\omega }_1}{{\omega }_2-{\omega }_1}}{\partial \alpha }=\frac{\partial {\omega }_1}{\partial \alpha }\frac{1}{\left({\omega }_2-{\omega }_1\right)}-\left(\frac{\partial {\omega }_2}{\partial \alpha }-\frac{\partial {\omega }_1}{\partial \alpha }\right)\frac{{\omega }_1}{{\left({\omega }_2-{\omega }_1\right)}^2}$$

(25)

The derivation of the sensitivity of the eigenvalues to the design variables can be seen as follows [25].

$$\mathit{K}{\mathit{\phi }}_j-{\omega }_j^2\mathit{M}{\mathit{\phi }}_j=\mathbf{0}$$

(26)

Then

$$\left(\frac{\partial \mathit{K}}{\partial \alpha }-2{\omega }_j\frac{\partial {\omega }_j}{\partial \alpha }\mathit{M}-{\omega }_j^2\frac{\partial \mathit{M}}{\partial \alpha }\right){\mathit{\phi }}_j+\left(\mathit{K}-{\omega }_j^2\mathit{M}\right)\frac{\partial {\mathit{\phi }}_j}{\partial \alpha }=\mathbf{0}$$

(27)

If the eigenvector is normalized with respect to $\mathit{M}$, then

$${\mathit{\phi }}_j^{\mathrm{T}}\mathit{M}{\mathit{\phi }}_k={\delta }_{kj}=\left\{\begin{array}{c}0k\ne j\\ 1k=j\end{array}\right.$$

(28)

Left-multiplying Equation (27) by ${\mathit{\phi }}_j^{\mathrm{T}}$, we have

$$\frac{\partial {\omega }_j}{\partial \alpha }=\frac{1}{2{\omega }_j}\left({\mathit{\phi }}_j^{\mathrm{T}}\frac{\partial \mathit{K}}{\partial \alpha }{\mathit{\phi }}_j-{\omega }_j^2{\mathit{\phi }}_j^{\mathrm{T}}\frac{\partial \mathit{M}}{\partial \alpha }{\mathit{\phi }}_j\right)$$

(29)

Calculations of $\frac{\partial \mathit{K}}{\partial \alpha }$ and $\frac{\partial \mathit{M}}{\partial \alpha }$ can be referred to as follows [24].

$$\begin{gathered} \frac{\partial \mathit{K}}{\partial \alpha }=\left(\frac{E}{4}\left(\sum _{k=1}^{NE}\sum _{j=1}^{4}q{\left(H\left({\chi }_j^e\right)\right)}^{q-1}\frac{\partial H\left({\chi }_j^e\right)}{\partial a}\right){\mathit{k}}^s\right. \\ \frac{\partial \mathit{M}}{\partial \alpha }=\left(\frac{\rho }{4}\sum _{k=1}^{NE}\sum _{j=1}^4\frac{\partial H\left({\chi }_j^e\right)}{\partial a}\right){\mathit{M}}^s \end{gathered}$$

(30)

Here, E and $\rho$ are the material elasticity and density, $NE$ is number of elements, $q$ is an integer, ${\chi }_j^e$ denotes topology description function, $H$ is Heaviside function, ${\mathit{k}}^s$ and ${\mathit{M}}^s$ represent the stiffness matrix and mass matrix of the element when $E=1$, $\rho =1$ and fully filled with material, respectively.

The sensitivity of the volume constraint function and compliance constraint function to the design variables is given by the following equation [24].

$$\begin{gathered} \frac{\partial \mathit{V}}{\partial \alpha }=\frac{1}{4}\sum _{k=1}^{NE}\sum _{i=1}^4\frac{\partial H\left({\chi }_j^e\right)}{\partial \alpha } \\ \frac{\partial \mathit{C}}{\partial {\alpha }_{}}={\mathit{U}}^{\mathrm{T}}\frac{\partial \mathit{K}}{\partial {\alpha }_{}}\mathit{U} \end{gathered}$$

(31)

Here, $\mathit{U}$ is the displacement vector of the structure.

3.2. Optimization Result of the Panel

The optimization process based on MMC method [19] is shown in Figure 8. Firstly, there is an initial design, then the structural topological description function is updated, then a finite element analysis is performed to calculate the structural dynamics eigenvalues and eigenvectors, followed by the calculation of the value and sensitivity of the constraint function and the objective function, and the design is updated. The above steps will be repeated until the convergence requirement is satisfied.

The initial design before the optimization is shown in Figure 9. The geometric parameters are the same as those of the standard structure, as shown in Table 1. The optimization results for the three constraint cases are shown in Figure 10. The optimization results can be easily imported into CAD software to generate geometric models (see Figure 11). The optimization result of Case A has the most complex form of reinforcement, which is not conducive to engineering manufacture. As the constraints are increased and strengthened, the reinforcement forms are simplified. The reinforcement form of Case C is a very simple reinforcement form and very similar to the standard structure, which can be manufactured easily.

4. Results and Discussion

In Section 3, the optimization results for the three constraint cases are obtained, here dynamics and aeroelastic characteristics of the three optimized panel structures will be simulated and discussed.

FEM modelling of the three optimization cases is proposed for the calculation of dynamics and aeroelastic characteristics (see Figure 12). The nodes and element information of the three models are shown in

Table 5. FEM modelling of three cases.

Case A		Case B		Case C
Nodes	Elements	Nodes	Elements	Nodes	Elements
7289	7300	7535	7524	7233	7082

.

Natural frequencies of the three optimized panel structures are shown in

Table 6. Natural frequencies of three optimized panel structures ($\mathrm{H}\mathrm{z}$).

Modes	Case A		Case B		Case C
	Normal Modes	Thermal Modes	Normal Modes	Thermal Modes	Normal Modes	Thermal Modes
1	78.92	78.34	79.08	78.46	84.30	82.24
2	165.34	163.77	160.44	158.88	169.86	165.68
3	169.89	168.31	161.78	160.26	170.50	166.32
4	215.02	213.16	201.37	199.48	212.01	206.96
5	247.95	246.29	232.20	230.92	240.79	234.71
6	278.23	276.96	246.05	244.66	262.11	255.57

, and the modes are shown in Figure 13, Figure 14, Figure 15, Figure 16, Figure 17 and Figure 18. For the optimized panel structure, the increase in temperature similarly causes natural frequencies of the structure to decrease while the mode remains unchanged. The frequency ratio $R$ of standard panel and optimized panels is shown in

Table 7. Frequency ratio $R$ of standard panel and optimized panels.

Frequency Ratio	Standard Structure	Case A	Case B	Case C
$\mathit{R}$	1.0445	0.9132	0.9720	0.9853

. Compared with the standard structure, the frequency ratios of the three optimized panel structures are all lower and less than 1, which means that the difference between the first and the second natural frequency of the optimized panel structure increases. The second frequency is greater than twice of the first natural frequency. It is proved that the frequency of the panel structure can be controlled by the proposed optimization design method in this paper. It also can be seen that as the constraints are increased or strengthened, the forms of panel reinforcement become simpler; however, the adjustment range of the frequency of the panel structure becomes smaller as well. So, a good balance between the optimization of aeroelasticity and manufacture is found.

The aeroelastic characteristics of the optimized panel are then evaluated. With the flutter calculation method described in Section 2, the critical flutter velocity of the optimized panel structure can be obtained, as shown in

Table 8. Critical flutter velocity of standard and optimized panels.

Flutter Velocity (m/s)	Standard Structure	Case A	Case B	Case C
Normal modes	1050	1200	1250	1110
Thermal modes	640	1140	1040	1080

. V-g and V-f plots of the three optimized panel structures are shown in Figure 19, Figure 20 and Figure 21.

By comparing the panel structures of the four forms of reinforcement, the panel structure of Case A has the smallest frequency ratio versus the largest frequency spacing, and its critical flutter speed with or without temperature load are relatively higher. The panel structure of Case C has the simplest form of reinforcement; its critical flutter speed is slightly lower than Case A, but higher than the standard structure. This indicates that as the constraints are increased and enhanced, the frequency range of the optimized structure becomes smaller and the increasement in critical flutter speed is also limited. Overall, the aeroelastic stability of the optimized panel structure is significantly improved compared to the standard structure.

5. Conclusions

In this paper, optimization of a surrounding fixed-support panel with reinforcing ribs is studied. The MMC method is proposed after the structural dynamics and flutter characteristics of the structure with or without considering thermal effects. From different results of three optimized cases, it could be concluded that:

• Increased temperature reduces natural frequencies and critical flutter velocity of the surrounding fixed-support square panel, which means that its aeroelastic stability should be reconsidered.
• The MMC optimization method can change the dynamic characteristics of the panel in both natural frequencies and modes. However, the reinforcement forms obtained under unconstrained conditions are complicated; thus, multiple constraints should be considered in engineering applications.
• As the number of constraints increases, the frequency adjustment range of the panel structure becomes smaller. However, compared to the standard structure, the critical flutter speed of optimized panel structure can be improved by adjusting the frequency difference.

The developed MMC method proposed in this paper can be used for dynamic and aeroelasticity optimization of hypersonic structures, which are significantly influenced by aero-heating effects.