Fluid–Structure Interaction Dynamic Response of Rocket Fairing in Falling Phase

Abstract

A method based on fluid–structure coupling is used in this study to calculate the response of a rocket fairing as it is falling. Some cases of vibration divergence of the fairing were found, and the influence of some specific factors was analyzed. The aerodynamic forces are calculated by using computational fluid dynamics (CFD) software and the structural responses by the modal-superposition method. The data are then subjected to modal interpolation in the CFD solver for the next cycle of calculation. The dynamic pressure, Mach number, and angle of attack are fixed in this process. Given that the fairing has a fixed attitude during falling, its rotation is ignored in calculations for the simulation. The results are then used to propose a framework for the fluid–structure coupling-based analysis of a non-streamlined structure. The mechanism of the fairing is discussed based on this method, and the effects of the settings of the solver, Mach number, dynamic pressure, and structural stiffness on it are investigated. Dangerous and safe regions are identified as the fairing falls back to the ground. Three methods are then provided based on the above analysis to prevent damage to the fairing as it falls to ground, such as increasing structure rigidity, attitude control, and opening the parachute at high altitude. A comprehensive method was used to suppress the vibration of the fairing during the descent, which was proven to be effective.

1. Introduction

The fairing of a rocket is expensive, and accounts for a large part of the cost of launching. When the fairing is separated from the rocket at a high altitude, it falls back to the ground along with the remaining wreckage. This poses a safety hazard, and the fact that the fairing cannot be reused significantly increases the cost of the rocket. Figure 1 shows typical fairings. Recycling and reusing the fairing is a popular subject of investigation in commercial aerospace research and a new scheme for recovering the fairing has been proposed here. Following launch, the fairing detaches from the rocket as the air thins high in the atmosphere. A parachute is used to ease its fall as the fairing speed slows to subsonic speed and it is recovered when it reaches the ground. As the flexible shell of the fairing falls at a high speed, it is subjected to complex fluid–solid coupling under the action of airflow before the parachute opens, and thus may suffer damage. It is important to study this complex problem of fluid–solid coupling. Analyses of the fluid–structure interaction (FSI) that induce the failure of the non-streamlined structure during its descent are carried out to prevent accidents and provide a reference for the design of the fairing.

Methods to calculate aeroelasticity in the context of aerospace research can be divided into three categories: theoretical algorithms for the frequency domain of aeroelasticity [1], methods to calculate aeroelasticity based on a reduced-order model (ROM) [2], and coupled numerical calculations using computational fluid dynamics (CFD) and computational structural dynamics (CSD) [3]. The doublet-lattice method is an example of methods of calculation in the frequency domain [4], and reduced-order aerodynamic models include the Volterra series and the radial basis neural network [5,6,7]. With advances in CFD and CSD, numerical calculations based on the theory of FSI provide novel means for dealing with problems of unsteady aeroelasticity.

Researchers in aviation have developed mature methods to analyze the streamlines of wings and aircraft. Torregrosa used the ROM method to study non-linear aeroelastic phenomena and applied the results to beams [8]. Xie established a non-linear, geometric FSI model based on the “quasi-modal” method of a non-constant vortical lattice combined with structural dynamics, and used it to examine the problem of the gust response of large-scale wings with a flexible chord ratio [9]. Kenneth used numerical calculations based on CFD/CSD coupling to simulate the movement of a two-dimensional wing with gaps and analyzed its flutter characteristics. This work formed the foundation for the subsequent investigation of problems of aeroelasticity [10]. Hallissy and Cesnik developed a highly precise tool for the aeroelastic simulations of flexible wings through weak CFD/CSD coupling [11]. To study the linear and non-linear static aeroelasticity of flexible wings, Mian developed a numerical method that weakly couples an open-source structural dynamics software with the Reynolds-averaged Navier–Stokes (RANS) equations. The results of calculations of the model of an aeroelastic wing by using this method were in good agreement with experimental data [12]. Ilie developed a time-domain numerical model based on bidirectional CFD/CSD coupling to simulate aeroelastic responses of helicopter blades [13]. Aprovitola focused on the static pitching stability of the aircraft to examine the effects of variations in the Mach number on the shift in its aerodynamic center by using both a commercial (ANSYS-Fluent) tool and an open-source (SU2) code [14]. Franzmann used an experimental method to measure the coefficient of the damping moment of the pitch of a rocket at a certain angle of attack and obtained results that were consistent with those of CFD simulations [15]. Zhong studied a control system that can reliably separate the fairing under a high dynamic pressure and suppress interference in this process [16]. Seiji used a hybrid LES/RANS-based method to simulate the transonic flow-field around a rocket fairing [17].

However, research on the fairing of the rocket has mainly focused on its acoustics and shape [18,19,20], and no references were found which studied issues of the rocket fairing recycling for engineering applications. The traditional doublet-lattice method cannot be applied to develop a flexible and non-streamlined configuration of the fairing. A CFD–CSD coupling method is required. Due to the lack of research on this object, the design department lacks cognition of influences of relevant paraments. It is worthwhile to study this problem.

In this paper, the motion of the fairing as it falls before the parachute opens is studied. Firstly, a framework for fluid–structure coupling-based analysis is proposed, which is based on CFD and the modal-superposition method to investigate the non-streamlined configuration of the rocket fairing and provide a dynamic response analysis for it reuse [21,22]. Unsteady aerodynamic forces and the structural response are calculated by a CFD software and the modal-superposition method. Data on the displacement are passed to the CFD solver through interpolation calculation for the next round of calculations. Then, simulations of the typical cylindrical vortex-induced vibrations that are consistent with past work are performed to verify the effectiveness of the proposed method. Furthermore, the dynamic response analysis in the recycling process of fairing is analyzed [23]. Typical working conditions are selected during the calculations to observe the stability of the response of the fairing while maintaining a constant dynamic pressure, Mach number, and angle of attack. The rotation of the fairing is ignored in the calculations. The influences of Mach number and dynamic pressure are discussed, and dangerous and safe regions for the descent of the fairing are identified. Three methods are then provided based on the above analysis to prevent damage to the fairing.

2. Basic Theories of FSI

Due to uncertainty in the movement of the fairing as it falls and limitations on its size, it is difficult to conduct wind tunnel tests. The motions of fluids and structures are often complex. Aerodynamic forces deform the structure of the fairing, and the displacement generated by its deformation changes the shape of fluids which cause structure deformation. The interaction between the fluid and the solid renders the problem unsteady in the time domain. Numerical methods are used to study this problem, which involve modeling the fluids and the structure as well as accurately calculating their responses and interactions.

The panel method is widely used to calculate the aerodynamic force due to its convenience of use, but it is not suitable for problems involving large deformation and non-linearity. Reduced-order models have received increasing attention in recent years, but their inadequate capability of generalization and the long training time needed by them limit their application. Unsteady aerodynamic forces were calculated by using CFD to accurately study the movement and mechanism of failure of the fairing as it falls.

Since the structure of the fairing is relatively simple, its main modes over several orders can reflect structural deformations in it. The modal-superposition method was used for structural calculations with several main modes to reduce the number of requisite calculations.

The numerical solution based on CFD/CSD coupling has advantages in terms of solving problems in aerospace research. Fluid–solid weak coupling is used to solve the problem, and the data are exchanged through the interface between the fluid and the solid.

Since structural displacement cannot be ignored during the descent of the fairing, bidirectional fluid–structure coupling is needed in calculations of the simulation.

2.1. Fluid Mechanics

A method based on the RANS equation is used to study the motion of the fairing in the atmosphere. The basic control equations of the fluid include the equation of continuity and the Navier–Stokes equation. They can be integrated into the following form [24,25]:

$$\frac{\partial }{\partial t}{\displaystyle {\iiint }_{V}W\mathrm{dV}+}{\displaystyle {\iint }_{\partial V}{F}_C}\mathrm{d}\mathrm{S}-{\displaystyle {\iint }_{\partial V}{F}_V}\mathrm{d}\mathrm{S}=0$$

(1)

where,

$$W=\left[\begin{array}{c}\rho \\ \rho \overrightarrow{v}\\ \rho E\end{array}\right]\hspace{1em}{F}_C=\left[\begin{array}{c}\rho V_r\\ \rho V\overrightarrow{v}+p\overrightarrow{n}\\ \left(\rho H+p\right)V_r+V_{t}p\end{array}\right]\hspace{1em}{F}_V=\left[\begin{array}{c}0\\ \overline{\tau }\cdot \overrightarrow{n}\\ k\left(\nabla T\cdot \overrightarrow{n}\right)+\left(\overline{\tau }\cdot \overrightarrow{v}\right)\cdot \overrightarrow{n}\end{array}\right]$$

(2)

In the above expressions, $W$ represents the vector of the conservative variables, $F_C$ and $F_V$ are the vectors of convective and viscous fluxes, respectively, $p$, $\rho$, $\overrightarrow{v}$, $T$, and $E$ denote the pressure, the density, the velocity vector, the temperature, and the total energy per unit mass, respectively, $V_t$ stands for the contravariant velocity of the face of a control volume and $V_r$ represents the contravariant velocity relative to the motion of the grid, $\overrightarrow{n}$ refers to the outward-facing unit normal vector of $\mathrm{dS}$, $k$ is the thermal conductivity coefficient, and $\overline{\tau }$ is the viscous stress tensor.

2.2. Basic Theory of Structural Dynamics

FSI was used to study the motion of the fairing during recovery. The equation of an n degree of freedom (DOF)-damped system’s forced vibration in the z direction is:

$$M\ddot{z}+C\dot{z}+Kz=f$$

(3)

where $z$ represents the vector formed by displacement of n degrees of freedom along the vibration direction, $M$, $C$, and $K$ denote the mass matrix, damping matrix, and stiffness matrix, respectively, and $f$ represents the load.

The natural modes of the model are ${\mathit{\varphi }}_1^{(n\times 1)} ~ {\mathit{\varphi }}_n^{(n\times 1)}$, and can be written as the following matrix:

$${\mathit{\varphi }}_1^{(n\times 1)} ~ {\mathit{\varphi }}_n^{(n\times 1)}$$

(4)

Multiplying the left side of Equation (3) by ${\Phi }^{\mathrm{T}}$ yields:

$${\Phi }^{T}M\ddot{z}+{\Phi }^{T}C\dot{z}+{\Phi }^{T}Kz={\Phi }^{T}f$$

(5)

The displacement vector $z$ can be represented by the generalized displacement $q$ corresponding to the mode of each order:

$$z=\Phi q={\displaystyle \sum _{i=1}^n{\mathit{\varphi }}_i{q}_i}$$

(6)

$${\Phi }^{T}M\Phi \ddot{q}+{\Phi }^{T}C\Phi \dot{q}+{\Phi }^{T}K\Phi q={\Phi }^{T}f$$

(7)

When the mode of each order is a mass-normalized mode, the damping matrix is ignored:

$$I\ddot{q}+C_{\mathrm{p}}\dot{q}+\Omega q={\Phi }^{T}f$$

(8)

where $I^{(n\times n)}$ is the identity matrix, $C_{\mathrm{p}}{}^{(n\times n)}=diag\left(2{\zeta }_i{\omega }_{\mathrm{i}}\right)$ is a diagonal matrix, and ${\Omega }^{(n\times n)}=diag\left({\omega }_i^2\right)$ is also a diagonal matrix. $M_p={\Phi }^{T}M\Phi$ represents the generalized mass matrix, and $\zeta$ is the damping ratio. When $\Phi$ consists of mass-normalized modes, $M_p$ is equal to $I$, $C_p={\Phi }^{T}C\Phi$ represents the generalized damping matrix, $K_p={\Phi }^{T}K\Phi$ represents the generalized stiffness matrix, and $K_p$ is equal to $\Omega$.

Generalized forces are involved in the calculation process. The program firstly reads the pressure of every aerodynamic grid at each time step, and then calculates the aerodynamic force of each aerodynamic grid. In the same time step, the mode shapes of the aerodynamic center are obtained from locations of nodes of the structure model and the mode shapes of the structure model through the Thin-Plate Spline (TPS) interpolation algorithm. The aerodynamic force is multiplied by the mode shape, and then by its area to obtain the generalized force for each aerodynamic mesh.

The generalized force is defined as:

$$F={\displaystyle \sum {\mathit{\varphi }}_i}{p}_i{S}_i$$

(9)

After the generalized aerodynamic force is obtained, it is substituted into Equation (8). Then, the generalized displacement is solved by the Runge–Kuta method. After the generalized displacement force is obtained, the modal-superposition method is used to solve the physical displacement of each node. Finally, the mesh in the CFD solver completes the deformation.

2.3. Thin-Plate Spline Interpolation Algorithm

When the modal method is used to deal with the problem of FSI, it usually involves using the grids of different nodes in the regions of the fluid and the solid. The modes at nodes of the fluid mesh are calculated by the TPS interpolation algorithm. A generalized dynamic equation is established at the central node of the grid and is solved to obtain the generalized displacement and generalized velocity. The nodal displacement of the boundary of FSI is obtained by modal value interpolation, and changes according to the generalized displacement once the motion of the fluid has been solved.

The TPS interpolation algorithm can be described by the following equation:

$$\left\{\begin{array}{l}{c}_{i1}=x_i\\ c_{i2}=y_i\\ c_{i3}=z_i\end{array}\right\},i\in [1\dots p]\equiv C_{p\times 3}=\left[\begin{array}{l}\begin{array}{ccc}{x}_1& y_1& z_1\end{array}\\ \begin{array}{cc}\begin{array}{cc}{x}_2& y_2\end{array}& z_2\end{array}\\ \begin{array}{ccc}& \dots & \end{array}\\ \begin{array}{ccc}{x}_p& y_p& z_p\end{array}\end{array}\right]$$

(10)

where p is the number of control points and $c_i\left(x_i,y_i,z_i\right),i=1,2,\dots ,p$ represents their location.

The regularization parameter, $\lambda$, can then be obtained, and the unknown TPS weights w and $a$ in the linear system of equations need to be solved.

$$\left[\begin{array}{cc}W& P\\ P^T& O\end{array}\right]\left[\begin{array}{c}{\overrightarrow{w}}^2\\ {\overrightarrow{a}}^2\end{array}\right]=\left[\begin{array}{c}\overrightarrow{v}\\ \overrightarrow{o}\end{array}\right]\equiv L_{\left(p+3\right)\times \left(p+3\right)}{\overrightarrow{x}}^2{}_{(p+3)}={\overrightarrow{b}}_{(p+3)}$$

(11)

where $W$, $P$, and $O$ are submatrices, and $\overrightarrow{\mathit{w}}$, $\overrightarrow{\mathit{a}}$, $\overrightarrow{\mathit{v}}$, and $\overrightarrow{\mathit{o}}$ are column vectors that are given by the following equations:

$$W_{ij}=U\left(\left||\begin{array}{cc}{c}_{i1}& c_{i2}\end{array}|-|\begin{array}{cc}{c}_{j1}& c_{j2}\end{array}|\right|\right)+I_{ij}\cdot {\alpha }^2\cdot \lambda \Vert i,j\in [1\dots p]\wedge \lambda \ge 0$$

(12)

$$U\left(r\right)=\{\begin{array}{l}{r}^2\cdot \log r\Vert r>0\\ \begin{array}{cc}0& \Vert r=0\end{array}\end{array}$$

(13)

$$\alpha =\frac{1}{p^2}{\displaystyle \sum _i^p{\displaystyle \sum _j^p\left(\left|\left|\begin{array}{cc}{c}_{i1}& c_{i2}\end{array}\right|-\left|\begin{array}{cc}{c}_{j1}& c_{j2}\end{array}\right|\right|\right)}}$$

(14)

$$P_{p\times 3}=\left[\begin{array}{ccc}1& c_{11}& c_{12}\\ 1& c_{21}& c_{22}\\ 1& ⋮& ⋮\\ 1& c_{p1}& c_{p1}\end{array}\right],O_{3\times 3}=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]$$

(15)

$$\begin{array}{cc}{P}^T{}_{ij}=P_{ij}& \Vert i\in [1\dots p]\wedge j\in [1\dots 3]\end{array}$$

(16)

$${\overrightarrow{v}}_{p\times 1}=\left[\begin{array}{c}{c}_{13}\\ c_{23}\\ ⋮\\ c_{p3}\end{array}\right],{\overrightarrow{o}}_{3\times 1}=\left[\begin{array}{c}0\\ 0\\ ⋮\\ 0\end{array}\right]{\overrightarrow{w}}_{p\times 1}=\left[\begin{array}{c}{w}_1\\ w_2\\ ⋮\\ w_p\end{array}\right],{\overrightarrow{a}}^2{}_{3\times 1}=\left[\begin{array}{c}{a}_1\\ a_2\\ ⋮\\ w_p\end{array}\right]$$

(17)

The value of z at any point can be interpolated by the following equation:

$$\tilde{z}\left(x,y\right)=a_1+a_{2}x+a_{3}y+{\displaystyle \sum _{i=1}^p{w}_{i}U\left(\left|\left|\begin{array}{cc}{c}_{i1}& c_{i2}\end{array}\right|-\left|\begin{array}{cc}x& y\end{array}\right|\right|\right)}$$

(18)

The generalized displacement of the fairing is analyzed to determine whether the vibrations are divergent.

2.4. Dynamic Mesh

Since structural deformation occurs in the process of fluid–structure coupling, many methods are available to update the mesh. Spring smoothing is used here to adjust the mesh to represent its deformation.

For fluid regions with triangular or quadrilateral meshes, an elastic deformation-based method of mesh adjustment can be used to smoothly adjust the positions of nodes in the domain of flow. This method can smoothly adjust the volume of mesh without changing the connectivity between meshes. The connection between any pair of nodes on the grid is idealized as springs connected to each other. A boundary node moves according to a given displacement, and causes all springs connected to this node to generate a force proportional to the displacement. In this way, the displacement of the nodes on the boundary is propagated through the volume of mesh in the fluid. From the perspective of balance, the resultant force of all springs at each node must be zero. This condition can be expressed by the following iterative equation:

$$\Delta {\overrightarrow{x}}^{m+1}{}_i=\frac{{\displaystyle {\sum }_j^{n_i}{k}_{ij}\Delta {\overrightarrow{x}}^m{}_j}}{{\displaystyle {\sum }_j^{n_i}{k}_{ij}}}$$

(19)

where $\Delta {\overrightarrow{x}}_i$ is the displacement of the node, $m$ is the number of iterations, $n_i$ is the number of nodes adjacent to node $i$, and $k_i{}_j$ is the spring constant between the given node and adjacent nodes, where the spring constant can be defined as:

$$k_{ij}=\frac{k_{fac}}{\sqrt{\left|{\overrightarrow{x}}_i-{\overrightarrow{x}}_j\right|}}$$

(20)

where $k_{fac}$ is the spring constant factor.

Since the displacement along the boundary is known, the equation is solved for sweeping all nodes inside the fluid domain through a Jacobi matrix. During the solution process, the updated positions of the node can be expressed as:

$${\overrightarrow{x}}_i^{n+1}={\overrightarrow{x}}_i^n+\Delta {\overrightarrow{x}}_i^{k,converged}$$

(21)

2.5. CFD–CSD Coupling

Figure 2 shows a weak coupling strategy to integrate the CFD and CSD methods.

The CFD method outputs the pressure at each moment, and then calculates the aerodynamic force at each aerodynamic face of the fairing. The mode shapes of the centers of the aerodynamic faces are obtained by using the TPS interpolation algorithm based on the locations of the structural nodes. The mode shapes are multiplied by the aerodynamic forces to obtain the generalized forces for each aerodynamic mesh. The generalized displacement is then obtained by using the generalized force in the fourth-order Runge–Kuta method to solve the equation of structural dynamics. The generalized displacement is multiplied by the mode shape to obtain the deformation of the mesh. Curves of the temporal response of the generalized force and the generalized displacement can be obtained by repeating the above process at each time step.

3. Research Object

3.1. Research Object and Working Conditions

Figure 3 shows the shape of the rocket fairing and provides the definitions of the coordinate and the angle of attack. The radius of the shell was approximately 1.5 m, the radius of the head was approximately 1 m, the length of the shell was 5 m, and the length of the fairing was 9 m. The head and fuselage were made of different materials.

Table 1. Material properties of each part of the fairing.

Fairing	Young’s Modulus (Pa)	Poisson’s Ratio	Density (kg/m3)
Head	$6.9\times {10}^{10}$	0.33	3200
Shell	$1.5\times {10}^{10}$	0.33	2500

lists the materials of each part of the fairing and their properties.

According to the typical state of a falling fairing, the Mach numbers $Ma=0.6$ and $Ma=0.5$ were selected as the typical working conditions for detailed calculation and analysis. The Mach number ranged from 0.2 to 0.85, and five values of the angle of attack were considered: 0°, 20°, 45°, 70°, and 90°. The effects of dynamic pressure, Mach number, angle of attack, and other factors on vibrations induced by FSI were explored by calculating for different conditions, as shown in

Table 2. Working conditions.

Working Condition	Dynamic Pressure (Pa)	Mach Number	Temperature (K)
1	700	0.5	218.6
2	1000	0.6	220

.

3.2. Grid for the Structure and Fluid

The fairing of the rocket was modeled by shells with stiffeners. Figure 4 shows a schematic diagram of the structural grid, and the red star in it represents the displacement-monitoring point.

A modal analysis of the above structure was carried out. The first four modes were selected, which may cause failure owing to structural vibrations. The first-order mode was the torsional mode, and the second-order mode was the breathing mode. Figure 5 shows the mode shapes and the frequencies of the fourth mode. The analysis yielded vibrational frequencies of the fourth-order modes of 2.2, 5.8, 9.1, and 11.5 Hz.

The coordinate system of the CFD model was consistent with that of the structural model. The area of fluid calculation was a spherical area around the fairing with a radius of 50 m. The spherical surface of the calculation area was set as the far-field boundary, and surface of FSI of the fairing was set as the boundary of the wall. The computational grids were tetrahedral mesh. The surface meshes of the fairing were triangular elements, and the total number of meshes was 1.88 million. The mesh was refined in the area near the fairing to better capture the flow-field details in the areas near the wall and in the wake area of the fairing, with drastic changes in velocity and pressure. A schematic diagram of the grid is presented in Figure 6.

The CFD solver was Ansys Fluent and generalized displacement was calculated by the modal-superposition method, which has been mentioned in Section 2.2. The modal-superposition method is a self-programmed C language program in UDF. The first four modes were selected for the response calculation with a time step of 0.005 s and a turbulence model of $k-\omega SST$.

4. Verification of Methods

4.1. Grid Independence

The above model of the fairing was used with a Mach number of $Ma=0.6$, dynamic pressure of $\overline{p}=800 \mathrm{Pa}$, and angle of attack of $\alpha ={90}^{°}$. Figure 7 compares the generalized displacement of three grid systems. Their number of grid elements and minimum volume are listed in

Table 3. The number of elements in the three grid systems.

Grid	Elements	Minimum Volume (m3)
Fine	2,362,564	$4.8\times {10}^{-8}$
Baseline	1,876,113	$5.1\times {10}^{-8}$
Coarse	1,517,761	$5.5\times {10}^{-8}$

. The mesh number was changed by adjusting the number of nodes in the circumferential and y directions on the fairing surface. In this process, the number of meshes in the inner spherical region is mainly adjusted to verify the influence of different grids. The results of the baseline grid and the fine grid were similar, and the amplitude of vibrations of the baseline grid was larger than that of the coarse grid. The generalized displacement was adequately captured by the baseline grid. Thus, the baseline grid system was used in subsequent calculations.

4.2. The Coupling Method

Since the object of calculation of the fairing was special, a single cylinder with a similar aerodynamic shape to the fairing was used for simulations of vortex-induced vibrations with a single degree of freedom. The aim was to verify the correctness of the method used to calculate the aerodynamic force and the deformation of the fairing under different working conditions. The result of the simulations was compared with results from past work to verify the correctness of the user-defined function (UDF) program [26,27,28]. In one time step, UDF was used to read the force from each grid and calculate the generalized force. Then, the UDF calculated the generalized displacement through the Runge–Kuta method and calculated the displacement through the modal-superposition method. Finally, the UDF updated the deformation of each node in turn. The CSD calculation was finished and the CFD solver started to run after this process.

Figure 8 shows the physical model of the vortex-induced vibrations of a single cylinder, simplified as a spring-mass-damper system. The relevant parameters were as follows: flow speed $U=0.146 \mathrm{m}/\mathrm{s}$, flow density $\rho =1.225 \mathrm{kg}/\mathrm{s}$, diameter $D=0.01 \mathrm{m}$, length $L=0.02 \mathrm{m}$, and kinematic viscosity $\mu =1.48\times {10}^{-6} {\mathrm{m}}^2/\mathrm{s}$. The non-dimensional mass of the cylinder was ${\mathrm{m}}^{*}=10$. To encourage high-amplitude oscillations, the structural damping coefficient was set to zero. The Reynolds number was $Re=100$ and the reduced velocity was $Ur=2.4~12$.

$$\begin{array}{l}{m}^{*}=\frac{m}{\frac{1}{4}\rho \pi D^2}\\ Ur=\frac{U}{f_{n}D}\end{array}$$

(22)

where $m$ is the mass of a cylinder with unit length and $f_n$ is its natural vibrational frequency.

Figure 9 shows variations in the maximum transverse displacement of the cylinder with the reduced velocity. The results of the simulation were close to those reported by Singh and Mittal [26,27,28]. Due to differences in mesh division and deformation methods, the maximum amplitude ratio was slightly higher than the reference values. However, its frequency locking range and the overall trend change were consistent. The reasons for the difference in the results of the other two articles can be found in their articles. A comparison of the variations in the maximum amplitude of the vortex-induced vibrations of a single cylinder with the reduced velocity, $Ur$, verified the correctness of the UDF program in this paper, which was used to read the force and adjust the mesh deformation.

5. Numerical Results

To analyze the influence of the flight of the fairing as it detaches from the rocket and falls to the ground, different dynamic pressures and angles of attack, in the Mach number range of 0.2 to 0.85, were used to calculate its characteristics of FSI.

A number of cases were calculated to determine the dangerous and safe regions for the fairing as it fell and are used here to propose several methods to assess the vibrational damage to it.

5.1. Typical Working Conditions

The structural responses of the fairing under an unsteady aerodynamic force with several angles of attack were calculated to assess its stability in different states as it fell to the ground.

Figure 10, Figure 11, Figure 12 and Figure 13 show the generalized displacement–time responses and fourth-order generalized force–time responses of the first four modes under two typical working conditions with different angles of attack ($\alpha =0^{°}, {70}^{°}, {90}^{°}$).

Figure 12 shows that the generalized displacement of the fairing gradually converged when $\alpha =0$, which shows that the vibrations were stable. When $\alpha ={70}^{°}$, the response of the first-order generalized displacement in the time domain was close to a constant oscillation in amplitude, indicating that the structure was in a critically stable state. When $\alpha ={90}^{°}$, the amplitude of the first-order generalized displacement gradually expanded and the vibrations of the structure diverged.

To further analyze the condition of divergence at $\alpha ={90}^{°}$, the displacement of the monitoring point at the end of the fairing was selected to draw an image in the time domain and its power spectral density was analyzed. Figure 14 shows that the deformation diverged in the X direction, which led to the destruction of the structure. This occurred mainly due to the influences of the first-order mode.

Figure 15 shows a diagram of structural deformation in one cycle when the motion of the fairing was unstable. Due to the divergence of first-order vibrations, the torsion of the junction between the cone and the cylinder became increasingly prominent and eventually led to vibration-induced damage to the fairing.

5.2. Influence of Mach Number

Since many variables influence the working conditions, a single variable was considered for analysis. The operational dynamic pressure and the angle of attack were fixed, and a series of Mach numbers were selected for calculation. The influence of different Mach numbers on the characteristics of FSI of the fairing was analyzed and the rules of their influences were summarized. Two work conditions were considered: $\overline{p}=500 \mathrm{Pa}$ and $\alpha ={70}^{°}$, and $\overline{p}=700 \mathrm{Pa}$ and $\alpha ={90}^{°}$.

To discuss the stability of the generalized displacement of the fairing, the definition of its logarithmic decay rate, $\delta$, with reference to damping was provided. Figure 16 shows the vibrational amplitude as a function of time. A value of $\delta$ greater than zero shows that the motion diverged. The larger the value of $\delta$ was, the more severe the vibrational divergence, and vice versa. The definition is as follows:

$$\delta =\frac{\ln \left(\frac{x_n}{x_0}\right)}{n}$$

(23)

where $x_n$ is the generalized displacement of the nth spike and $x_0$ is the generalized displacement at the selected initial spike.

The settings of the CFD solver were $\overline{p}=600 \mathrm{Pa}$ and $\alpha ={90}^{°}$. The FSI-induced vibrations of the fairing under different Mach numbers were then calculated. Figure 17 shows that the logarithmic decay rate of the fourth-order generalized displacement continued to increase with the Mach number. In the vicinity of $Ma=0.6$ and $Ma=0.7$, the first-order vibrations diverged, and FSI-induced failure occurred. However, as the Mach number continued to increase to around 0.8, the vibrations converged.

The operational condition of $\overline{p}=800 \mathrm{Pa}$ and $\alpha ={45}^{°}$ was similarly considered. Figure 18 shows that the first-order generalized displacement diverged in the Mach number range of 0.3 to 0.5. The first-order vibrations converged as the Mach number continued to increase. A comparison of the two results of the calculations shows that it is necessary to pay attention to different intervals of the Mach number in case of flutter divergence in case of first-order vibrations under different flight conditions. Since the non-linear factors are obvious in the whole calculation process, the mechanism is very complicated. According to the analysis, the primary reason for this phenomenon may be that the position and intensity of the shock wave generated on the aircraft surface are different at different angles of attack.

5.3. Influence of Dynamic Pressure

The influence of different flight pressures on the FSI-induced vibration divergences of the fairing was analyzed. Two working conditions, with $Ma=0.4$ and $\alpha ={45}^{°}$, and $Ma=0.6$ and $\alpha ={90}^{°}$, were selected for the calculation. Figure 19 and Figure 20 show that the problem of divergence of the first-order modal vibrations was encountered. As the dynamic pressure increased, the logarithmic decay rate corresponding to the first-order generalized displacement continued to increase, which means that the degree of divergence was more severe. The fairing was damaged due to the divergence of the first-order vibrations, which is the first-order torsional mode-induced damage of concern in this study.

5.4. Summary of the Results of Calculations

The fall of the fairing was calculated in different states, and the results were summarized. Figure 21 shows the hazard and safety zones at different angles of attack (green points represent safe states, the orange points represent the critical stable states, and the red points represent the states of structural failure).

When $\alpha =0$, the dangerous zone was concentrated in the range of $Ma=0.4$. As the angle of attack increased, this zone changed to the range of $Ma=0.6$. As the angle was further increased, the zone of hazard continued to expand.

6. Improvement Methods

6.1. Increasing Structural Rigidity

An analysis of the structural displacement of the fairing shows that damage to it had been caused mainly by the first-order torsion mode. Thus, its structural stiffness was enhanced and considered to suppress the damage caused by vibrations. Given that a variety of variables affect structural stiffness and that it is inconvenient to try to control all of them, a relatively simple method of changing the thickness of the shell was chosen to change its structural stiffness.

The working condition of $Ma=0.6$, $\overline{p}=600 \mathrm{Pa}$, and $\alpha ={90}^{°}$ was selected for calculation and comparison. The settings of the CFD solver were the same as before.

To increase the structural rigidity, the fairing shell was thickened to 1.2, 1.5, and 1.8 times its original thickness. The mass and thickness of the shell exhibited a quadratic relationship, and its stiffness and thickness had an approximately cubic relation. The mode shapes did not change in the meantime. Taking 1.5 times the thickness as an example for illustration, the modal frequency was increased to 1.5 times that of the original, and the vibration frequencies of the fourth-order modal were 3.3, 8.7, 13.7, and 17.3 Hz, respectively. Therefore, the unit mass matrix was multiplied by the double of 1.5 and the stiffness matrix by the cube of 1.5. Other calculation parameters remained unchanged. The generalized displacement of the first-order mode was then observed to determine whether the vibrations diverged.

Figure 22 and Figure 23 show the time-dependent curves of the generalized displacement corresponding to the fourth-order mode at different thicknesses.

As the thickness of the fairing increased, the first-order vibrations gradually transitioned from divergence to convergence. This verifies that enhancing the structural stiffness can somewhat suppress vibrational divergences, where this can be used to suppress flutter.

6.2. Attitude Control

Figure 21 shows the dangerous zone of fairing, and then it is necessary to find a feasible path to safely recycle the fairing. However, the dangerous zone was not the same under different angles of attack, and it was impossible to safely land the fairing while maintaining a certain angle of attack. A strategy of controlling the angle of attack within different angles at different Mach numbers was proposed to avoid dangerous zones that can cause vibrational damage to the fairing as it falls.

SpaceX adjusts the attitude and speed of the fairings of its rockets through nitrogen propulsion and parachutes built into the fairing [29]. In this same way, the movement of the separated fairing is controlled so that it does not disintegrate. Its landing position is accurately determined at the same time.

Figure 24 shows several control schemes for the angle of attack during the process of falling of the fairing. The green and yellow lines represent the states of change of the fairing during descent. The orange dots represent its critical states and the red dots the dangerous states. However, this method encounters technical complications, such as the need for an extra control device and the complex transformation of the fairing.

6.3. Opening the Parachute at High Altitude

Normally, the parachute will open at a certain altitude. However, the fairing disintegrates before the parachute opens due to damage caused to it by the vibrations. The fairing cannot be recycled. Figure 24 shows that the damage in our simulations was mainly concentrated in the Mach number range of 0.4 to 0.7. The parachute of the fairing can be opened before its speed drops to $Ma=0.7$, such that it can avoid the entire dangerous zone, and thus, disintegration.

This measure does not require complex modifications and redesign. However, the structural properties of the fairing change after the parachute opens, and this may lead to new problems related to its damage or destruction. It then becomes necessary to identify the novel states of the fairing and zones of hazard for it.

6.4. Combination of the Three Methods

Each of the above three solutions has advantages and disadvantages, which are listed in

Table 4. Advantages and disadvantages of the three methods.

Method	Advantages	Disadvantages
Increasing structural rigidity	Easy to alter	Increased launch cost Increased weight
Attitude control	Light in added weight	Complicated to alter
Opening the parachute at high altitude	Relatively light in added weight	Increased certain launch cost Relatively complicated to alter

.

The effectiveness of a single method is limited. In practice, we can use a combination of different methods. To prove the validity of the method, increasing structural rigidity and attitude control were chosen to ensure the stability of the fairing during the falling process. Path one (green) in Figure 24 and 1.5 times the original thickness were considered.

Figure 25 shows the detailed situation of path one during the falling process, such as Mach number, dynamic pressure, and angle of attack. Figure 26 shows that the vibration of the first-order generalized displacement converged in the whole path. This verifies that this comprehensive method can suppress vibrational divergences and could be used to ensure the safety of the fairing in the falling process.

7. Conclusions

This study used the coupled CFD–CSD method to simulate the fall of the fairing of a rocket. The effects of the angle of attack, Mach number, and dynamic pressure during its descent were investigated. Safe and dangerous zone for it were identified during its descent, and three improvements were proposed to avoid its destruction. The conclusions of this study can be summarized as follows.

• A framework of non-streamlined configurations with fluid–structure interactions was established. Several examples were used to verify that the proposed method can be used to calculate the FSI of the fairing and confirmed that the theoretical results corresponded to the actual situation. The work here provides ideas for future research on FSI involving objects with similar non-streamlined configurations.
• Through the analysis and calculation of multiple working conditions, the dangerous zone and safe zone as the fairing fell were obtained, which were consistent with the actual falling situation. When $\alpha \le {45}^{°}$, the hazardous zone occurred at $Ma=0.4$; when $\alpha \ge {45}^{°}$, it was concentrated in the Mach number range of 0.6 to 0.8. At the same time, the influences according to dynamic pressure and Mach number were also analyzed.
• According to the analysis, there is a risk of structural damage to the fairing as it falls. To suppress the vibration, a variety of possible methods were presented, such as enhancing the structural rigidity, flight attitude control, and opening the parachute at a high altitude. To verify the effectiveness of the method, a comprehensive method was used to calculate the vibration of the fairing during the descent. The fairing can land safely and avoid disintegration during the descent.