2.1 Model

Four extents of the right wing damage are considered. The extent of damage is indicated by the missing wing span length. The maximum damage to the wing is 50%, which represents the worst-case damage reproducible during environmental drops performed by the Close Combat Weapons Systems Project Office and the US Marine Corp [Reference Lanham20]. These parameters may not be consistent with the actual damaged situations, but they help to obtain the relationship between the change of aerodynamic performance and the damage extent of missile wing.

The model used in this study is a simplified version of a typical supersonic cruise missile. In order to conveniently make a comparison before and after damage, the axisymmetric inlet and aerofoil of the model are simplified, and the case of one wing damage is discussed in this paper. The parameters of geometric model are shown in Fig. 1, and they are measured in millimeters. W represents the wing and R represents the rudder in Fig. 1(b). h represents the damaged width in Fig. 1(c).

Figure 1. Geometric model and damage diagram.

The wing of the missile is thin and light, so the change in the centre of gravity can be neglected. The parameters of the missile are shown in Table 1. The reference area is the maximum cross-sectional area of the missile body, and the reference length is the diameter of the missile. The moment reference point (MRP) is 2.5m from the missile head and it is also set as the point of rotation (POR) for dynamic simulation.

Table 1. Basic parameters of missile

2.2 Mesh

The research adopts structured meshes. High-speed compressible flows involve important fluid phenomena typically exhibiting shock waves, rarefactions and related complexities [Reference Struchkov, Kozelkov, Volkov, Kurkin, Zhuckov and Sarazov21]. Therefore, the mesh needs to be encrypted where the flow changes violently. For example, the mesh is encrypted when the geometric structure of the model changes. Figure 2 gives the computational mesh, including far field mesh and boundary layer mesh. The far-field length in the computational domain is 10 times the length and 20 times the diameter of the projectile.

Figure 2. Mesh generation of the missile.

Mesh independence is crucial to get the correct solution. An efficient mesh resolution must be determined to find correct calculation results and the appropriate amount of solution time [Reference Sumnu and Güzelbey22]. The solution in this study is implemented by generating different mesh resolution between about 1.2 and 6.2 million. Figure 3 shows mesh independence for the drag coefficient (C _d) versus the number of mesh resolution.

Figure 3. Mesh independence validation.

The calculation result of 6.2 million mesh resolution is used as reference values. Table 2 shows that the change rate of C _d is all less than 2%. Since the solution results are negligible changes above 2.85 million mesh resolution, considering the calculation time cost, 2.85 million is chosen as the efficient mesh resolution.

Table 2. Mesh independence error of Cd

3.1 Fluid governing equation

The continuity equation and energy equation are referred to the Ref. [Reference Elisov, Ishkov and Shakhov23]. The conservation form of the 3D dimensionless, unsteady, compressible, Reynolds-averaged Navier–Stokes equations in Cartesian coordinates is as follows [Reference Yang, Wang and Gao24]:

\begin{align*}\frac{{\partial {{\boldsymbol{{Q}}}}}}{{\partial t}} + \frac{{\partial {{\boldsymbol{{F}}}}}}{{\partial x}} + \frac{{\partial {{\boldsymbol{{G}}}}}}{{\partial y}} + \frac{{\partial {{\boldsymbol{{H}}}}}}{{\partial z}} = \left[ {\frac{{\partial {{{\boldsymbol{{F}}}}_v}}}{{\partial x}} + \frac{{\partial {{{\boldsymbol{{G}}}}_v}}}{{\partial y}} + \frac{{\partial {{{\boldsymbol{{H}}}}_v}}}{{\partial z}}} \right]\end{align*}

\begin{align*}{{\boldsymbol{{Q}}}} = \left[ \begin{array}{l}\rho \\[5pt] \rho u\\[5pt] \rho v\\[5pt] \rho w\\[5pt] \rho E\end{array} \right],{{\boldsymbol{{F}}}} = \left[ \begin{array}{l}\rho u\\[5pt] \rho {u^2} + p\\[5pt] \rho uv\\[5pt] \rho uw\\[5pt] \rho Hu\end{array} \right],{{\boldsymbol{{G}}}} = \left[ \begin{array}{l}\rho v\\[5pt] \rho uv\\[5pt] \rho {v^2} + p\\[5pt] \rho vw\\[5pt] \rho Hv\end{array} \right]\end{align*}

(1) \begin{align}{{\boldsymbol{{H}}}} = \left[ \begin{array}{l}\rho w\\[5pt] \rho uw\\[5pt] \rho vw\\[5pt] \rho {w^2} + p\\[5pt] \rho Hw\end{array} \right],{{{\boldsymbol{{F}}}}_{{\boldsymbol{{v}}}}} = \left[ \begin{array}{c}0\\[5pt] {\tau _{xx}}\\[5pt] {\tau _{xy}}\\[5pt] {\tau _{xz}}\\[5pt] {\tau _{xx}}u + {\tau _{xy}}v + {\tau _{xz}}w + {q_x}\end{array} \right]\end{align}

\begin{align*}{{{\boldsymbol{{G}}}}_{{\boldsymbol{{v}}}}} = \left[ \begin{array}{c}0\\[5pt] {\tau _{yx}}\\[5pt] {\tau _{yy}}\\[5pt] {\tau _{yz}}\\[5pt] {\tau _{yx}}u + {\tau _{yy}}v + {\tau _{yz}}w + {q_y}\end{array} \right]\end{align*}

\begin{align*}{{{\boldsymbol{{H}}}}_{{\boldsymbol{{v}}}}} = \left[ \begin{array}{c}0\\[5pt] {\tau _{zx}}\\[5pt] {\tau _{xy}}\\[5pt] {\tau _{zz}}\\[5pt] {\tau _{xx}}u + {\tau _{zy}}v + {\tau _{zz}}w + {q_z}\end{array} \right]\end{align*}

where p is pressure, $\rho$ is density, (u, v, w) are the Cartesian velocity components, E is total specific energy, H is total specific enthalpy, H = E + p/ $\rho$ , Q is the conserved variable, F , G , H are the convective terms, F _v , G _v and H _v are the viscous terms.

3.2 Discrete scheme and solution method

The discrete scheme used in this study is Roe-FDS. The spatial discrete scheme is the second order upwind. The solution method is to adopt the density-based solution, velocity component and density are taken as basic variables, and the pressure is obtained by state equation. The density based solution is solved by coupling algorithm. The continuity equation, momentum equation and energy equation are solved simultaneously, and then other scalar equations are solved sequentially.

3.3 Turbulence model

Spalart-Allmaras (S-A) turbulence model is a simple single equation model for a modified turbulent kinematic viscosity [Reference Spalart and Allmaras25]. It is specially designed for aerospace applications. And it has a good application in the calculation of missile external flow field [Reference Meng, Yan, Huang and Tong26]. Compared with the two-equation model, it has lower encryption requirements for the computing grid on the wall. And it has the advantage of good stability. Previous studies [Reference Render, Samaad-Suhaeb, Yang and Mani27] on battle-damaged aerofoils had shown this model produced reasonable result, and outperformed the k- $\varepsilon$ model slightly. The transport equation of the Spalart-Allmaras model is adopted as the governing equation, and its form can be obtained from the Refs [Reference Jinbum and Seungsoo28, Reference Spalart and Allmaras29].

3.4 Boundary and initial conditions

Pressure-far-field condition is adopted at the entrance and radial boundary of the computational domain. Pressure-out condition is adopted at the outlet, the temperature and pressure are the same as the far-field. The surface of missile adopts non-slip solid wall boundary condition. The initial value of the computational flow field is set as the free flow value. The fluid types in the computational domain is ideal gas, and the viscosity is calculated using Sutherland’s law.

3.5 Dynamic derivative identification method

Forced motion offers a good alternative to predict the aerodynamic coefficients of an aircraft [Reference Allen and Ghoreyshi30]. The classic method for calculating dynamic derivatives is the small amplitude oscillation technique. It evolved from the wind tunnel tests and can calculate subsonic, transonic and supersonic states [Reference Wu, Lv and Cao31]. The Fourier expansion of unsteady aerodynamic moment on rigid aircraft with small amplitude forced pitch vibration is given as follows:

(2) \begin{align}{M_z} = {M_{{z_0}}} + {\bar M_z}\sin\!(\omega t + \lambda ) + u(t)\end{align}

According to the concept of aerodynamic derivative, the unsteady aerodynamic moment of the rigid aircraft with small amplitude forced pitch vibration can be expressed by Taylor expansion as follows:

(3) \begin{align}{M_Z} = {M_{{z_0}}} + M_z^\alpha {{\Delta }}\alpha + M_z^{\dot \alpha }{{\Delta }}\dot \alpha + M_z^{{\omega _z}}{\omega _z} + M_z^{{{\dot \omega }_z}}{\dot \omega _z} + {{\hat{\Delta} }}\left( {{{\Delta }}\alpha, {\omega _z}} \right)\end{align}

Where, $\;{{{M}}_{{z}}}$ is the transient unsteady aerodynamic moment; ${{{M}}_{{{{z}}_{{0}}}}}$ is the aerodynamic moment at the equilibrium position. ${{M}}_{{z}}^{{\alpha }}$ is the static derivative of aerodynamic moment; ${{M}}_{{z}}^{{{\dot \alpha }}}$ is the first dynamic derivative of aerodynamic moment to AoA; ${{M}}_{{z}}^{{{{\omega }}_{{z}}}}$ is the zero-order dynamic derivative of aerodynamic moment to pitch angular velocity; ${{M}}_{{z}}^{{{{{\dot \omega }}}_{{z}}}}$ is the first dynamic derivative of aerodynamic moment to pitch angular velocity; $\hat{\Delta}$ is a higher-order derivative term. When the rigid aircraft oscillates in small amplitude, its motion equation can be simplified as:

(4) \begin{align}\left\{ {\begin{array}{*{20}{l}}{\theta = {\theta _m}\sin\!(\omega t)}\\[5pt] {\dot \theta = \omega {\theta _m}\cos\!(\omega t) = {\omega _z}}\\[5pt] {\ddot \theta = - {\omega ^2}{\theta _m}\sin\!(\omega t) = {{\dot \omega }_z}}\\[5pt] {{{\Delta }}\alpha = \theta = {\theta _m}\sin\!(\omega t)}\\[5pt] {{{\Delta }}\dot \alpha = \dot \theta = \omega {\theta _m}\cos\!(\omega t)}\end{array}} \right.\end{align}

Substituting Equation (4) into Equation (3) and omitting high-order component $\hat{\Delta}$ :

(5) \begin{align}\eqalign{ & {M_Z} = (M_Z^\alpha - {\omega ^2}M_Z^{{{\dot \omega }_Z}}){\theta _m}\sin\!(\omega t) \cr & + (M_Z^{\dot \alpha } + M_Z^{{\omega _Z}})\omega {\theta _m}\cos\!(\omega t) + {M_{{Z_0}}} \cr} \end{align}

When the calculation time of the unsteady problem is long enough, the aerodynamic moment changes periodically, and the single point method can be adopted [Reference Oktay and Akay32], making ωt = 2nπ to remove the influence of the initial effect, so Equation (5) can be written as follows:

(6) \begin{align}M_z^{\dot \alpha } + M_z^{{\omega _z}} = \frac{{{{\bar M}_z}\sin \lambda }}{{\omega {\theta _m}}} = \frac{{{M_{\omega t = 2n\pi }} - {M_{{z_0}}}}}{{\omega {\theta _m}}}\end{align}

Make dimensionless treatment on the above formula. The reduced frequency is k = ωl/2V _∞ , where l is the reference length and V _∞ is the far-field free inflow velocity, substituting and obtaining:

(7) \begin{align}{C_{mq}} + {C_{m\dot \alpha }} = \frac{{{M_{\omega t = 2n\pi }} - {M_{{z_0}}}}}{{kqsl{\theta _m}}} = \frac{{{C_{{M_{\omega t}}_{ = 2n\pi }}} - {C_{{M_{{z_0}}}}}}}{{k{\theta _m}}}\end{align}

where ${{{C}}_{{{{M}}_{{{\omega t = 2n\pi }}}}}}$ is the instantaneous pitching moment coefficient at ωt = 2nπ, and ${{{C}}_{{{{M}}_{{{z0}}}}}}$ is the steady moment coefficient at the equilibrium position. k is the reduced frequency, q is the dynamic pressure, s and l are the reference area and length respectively, and ${{{\theta }}_{{m}}}$ is the amplitude.

3.6 Method validation

The Basic Finner missile model is the international standard model to validate the aerodynamic calculation. The model of the missile is shown in Fig. 4. The structural parameters are consistent with those in Ref. [Reference Bhagwandin and Sahu33], that is, d = 30mm. The computational mesh of the missile is shown in Fig. 5, and the number of grids is 2.6m. The length of the far field in the calculation domain is 10 times the length of the missile and the diameter is 20 times the diameter of the missile. The calculation conditions in this study are consistent with the Ref. [Reference Mi34], that is, Ma = 1.58, initial attack angle θ ₀ = 0°, frequency reduction rate k = 0.0158226, amplitude θ _m = 1°, so the oscillation law can be written as $\alpha$ = θ ₀ + θ _m sin(ωt) = 1°sin(571.7t).

Figure 4. Basic Finner missile model.

Figure 5. Mesh generation of Basic Finner missile.

Experimental data was obtained from a combination of free-flight tests in a ballistic range and wind tunnel measurements at different test facilities [Reference Dupuis and Hathaway35, Reference Dupuis36].

The hysteresis loops of pitch moment coefficient obtained from the calculation in this study is shown in Fig. 6. The curve basically coincides with the CFD calculation value in the Ref. [Reference Mi34], indicating that the error is small. The errors of $C_{{m}{\dot{\alpha}}} + C_{mq}$ are shown in Table 3. Error-C represents the error with the calculated value [Reference Green, Spence and Murphy37], and Error-E represents the error with the experimental value [Reference Bhagwandin and Sahu33]. Compared with the experimental and calculated values, the error of the S-A model obtained in this study is larger, but it is still less than 8%. One aspect of the analysis that contributes to this error may be due to insufficient number of grids and insufficient precision of numerical methods. On the other hand, the edge of the numerical calculation model is sharp, but the test piece usually has a certain chamfer, not a completely sharp edge, which may also cause errors between the numerical calculation and the test results. This error is acceptable, indicating that the calculation method has high reliability for calculating the aerodynamic characteristics of the missile model in this study.

Figure 6. Hysteresis loops of Cmy (Cmy represents the pitch moment coefficient).

Table 3. ${{{C}}_{{{m}{\dot \alpha }}}}{{\; + \;}}{{{C}}_{{{mq}}}}$ comparison of Basic Finner missile

4.1 Static aerodynamic characteristics

First, the static aerodynamic characteristics of all models with different damage extent are simulated. The calculation conditions are listed in Table 4.

Table 4. Static calculation conditions

The differences before and after damage are analysed in terms of the flow field. Due to large number of calculation conditions and limited space, only some representative flow field contours of the damaged models are selected for display, i.e., the models with the right wing of 0.3l and 0.5l damage at $\alpha$ = 4° and $\alpha$ = 8°.

Figures 7 and 8 show the Mach and local pressure contour around the models. The oblique shock waves formed in front of the head, wings, rudders and fins of the missile can be clearly seen in the figures. The low Mach number region at the tail of the missile is also obvious. By comparison, it is found that when the AoA changes from $\alpha$ = 4° to $\alpha$ = 8°, a local high Mach number region appears behind the wing on the undamaged side (left side). On the side of the damaged wing, the region is small. When the damage extent reaches 0.5l, the above phenomenon hardly appears, which makes the flow field on both sides of the missile asymmetric.

Figure 7. Mach and local pressure contour at $\alpha$ = 4°.

Figure 8. Mach and local pressure contour at $\alpha$ = 8°.

Figures 9 and 10 show the streamlines and vortices around the models. After being damaged, the flow field asymmetry becomes larger on both sides of the missile. This asymmetry is more obvious when the AoA increases. The vortices on both sides of the missile differs significantly. The vortices formed by the damaged wing are longer but smaller in diameter compared to the undamaged left side. The streamlines show that the vortices on the damaged side become more concentrated. Compared with the left wings, the vortices formed by the right wings of the two damaged models are closer to the body of the missile and extend to the rear side of the fin at $\alpha$ = 8°, which affects the flow field around the rudders and fins behind the wings.

Figure 9. Streamlines and vortices at $\alpha$ = 4°.

Figure 10. Streamlines and vortices at $\alpha$ = 8°.

To further explore the effect of wing damage on the flow field around the missile, the pressure distribution around the wings, rudders and fins of the missile are shown separately to compare their differences (the missile body is hidden for better visibility). Due to the limitation of space, the representative pressure contour distribution at $\alpha$ = 8° is selected for display. As shown in Fig. 11, the wingtip vortices of the damaged wings are smaller than that of the undamaged wing and appear mainly in part of the wing parallel to the body. However, the pressure distribution at the leading edge of the wing is almost unaffected by the damage. As the damage extent increases, the wingtip vortices become smaller and the pressure gradient changes more drastically. This leads to a decrease in the lift coefficient and a non-linear relationship with the damage extent at the same AoA. The wingtip vortices effect of damaged models is weakened, leading to reduction in lift force along with a reduction in lift-induced drag force.

Figure 11. Pressure contour around wings at $\alpha$ = 8°.

Figures 12 and 13 show that the damaged wing also affects the flow field around the rudders and fins at the rear of the body. Due to the partial loss of right wing, the airflow through the right side of the missile body is faster and the pressure is lower, and a local low pressure region is formed between the upper and lower rudders and fins on the right side. The pressure distribution on the right rudders and fins, especially on the right upper rudder and the fin surface, produces more obvious changes, this is, the pressure distribution on the upper and lower surfaces is not symmetrical. The increased damage extent also has a great influence on the flow field around the rudder and fin at the right upper position. As a result, a counterclockwise roll moment is generated by the rudders and fins of the missile, respectively, to counteract the clockwise roll moment that occurs due to the partial missing of the right wing.

Figure 12. Pressure contour around rudders at $\alpha$ = 8°.

Figure 13. Pressure contour around fins at $\alpha$ = 8°.

Figure 14 shows the static aerodynamic force and moment coefficient. Compared to the undamaged model, both the drag and lift coefficient of the damaged model are smaller due to the loss of both windward and lift surfaces. With the increase of damage extent, both C _d and C _l are sensitive to the change with AoA. The rate of change is listed in Tables 5 and 6. The change rate of the drag coefficient increases with increasing damage extent as well as with increasing AoA, while the change rate of the lift coefficient increases with increasing damage extent and decreases with increasing AoA. Meanwhile, Table 6 also shows the nonlinear relationship between the reduction of lift coefficient and the damaged extent at the same AoA.

Figure 14. Static aerodynamic forces and moments.

Table 5. C_d change at different AoA

Table 6. C_l change at different AoA

The pitch moment C _my of the missile is negative, which indicates that it can restrain the pitch motion and the missile is statically stable. The figure shows that as the AoA increases, the pitch moment values of several damaged models still basically obey linear changes. The more the missile wing is damaged, the larger the static pitch moment is, but the absolute value is decreasing, and the suppression of pitch motion decreases. That is, the static stability of the missile becomes worse. Due to the front of the missile wing is triangular, the surface area decreases more behind the MRP than in the front when the missile loses partial right wing.

The roll moment C _mx of the undamaged model is nearly 0 at all the angles of attack. As the area loss on the right wing increases, the absolute value of the roll moment increases and gradually changes from linear to non-nonlinear. This is due to the missing lifting surface area of the missile’s right wing is not linearly related to the missing wingspan.

4.2 Dynamic aerodynamic characteristics

Based on the static solutions, dynamic grid method is used to simulate unsteady small amplitude forced vibration to calculate dynamic aerodynamic parameters. The pitch motion is defined as $\alpha$ = θ ₀ + θ _m sin(ωt) = θ ₀ + 0.25° sin(823.2t), and the corresponding reduced frequency is k = ωl/2V _∞ = 0.24. The roll motion is defined as $\alpha$ = θ ₀ + θ _m sin(ωt) = θ ₀ + 0.25° sin(343t), and the corresponding k is 0.1. Other parameters are the same as those in the static calculation.

Due to the unsteady flow involved, the independence of the time step must be validated before the dynamic characteristics calculation. The time step t _y = 1×10⁻⁵s, 2×10^-5s and 4×10^-5s are used to calculate the pitch moment coefficient of the undamaged model at $\alpha$ = 0°. And the time step t _x = 1×10^-4s, 5×10^-5s and 2.5×10^-5s are used to calculate the roll moment coefficient of the undamaged model at $\alpha$ = 0°. The results in the Fig. 15 show that there is almost no difference in the calculation results. Considering the calculation time cost and calculation accuracy, the time step for calculating the pitch motion is t _y = 2×10^-5s, and the time step for calculating the roll motion is t _x = 5×10^-5s in this study.

Figure 15. Validation of moment coefficient with time during pith/roll motion (time step independence validation).

4.2.1 Pitch dynamic derivative

The hysteresis loops of the pitch moment coefficient are presented in Fig. 16. The hysteresis loops is obviously elliptical in shape, and the angle between its major axis and the x-axis is defined as the azimuth, as shown in Fig. 17. The circumferential direction of all curves is counterclockwise, which indicates the missile obtained damping characteristics of pitch moment at all attack angles. Although the different extents of the missing right wing parts, the azimuths and sizes of the hysteresis loops hardly change at the same AoA, only causing different areas of the hysteresis loops: the hysteresis loops area increases slightly with increasing of the damage extent.

Figure 16. Hysteresis loops of pitch moment coefficient at different AoA.

Figure 17. Diagram of the azimuth.

Comparing the hysteresis loops at different AoA, the azimuths change. The azimuths of $\alpha$ = 4° are opposite to that of $\alpha$ = 8° and $\alpha$ = 12°. Under the condition of the $\alpha$ = 4° and the same absolute value of oscillation instantaneous AoA, the longitudinal damping value of positive instantaneous AoA is larger than that of the negative instantaneous AoA. On the contrary, the longitudinal damping value at the negative instantaneous AoA is larger than that at the positive instantaneous AoA under the condition of $\alpha$ = 8°. The results of $\alpha$ = 12° is the same as that of $\alpha$ = 8°.

The pitch combined dynamic derivatives are shown in Fig. 18. With the increase of the AoA, the pitch combined dynamic derivatives ${{{C}}_{{{m}{\dot \alpha }}}}{{\; + \;}}{{{C}}_{{{mq}}}}$ of the damaged model and the undamaged models all increase at first and then decrease. This phenomenon is consistent with the changing trend of azimuth with the AoA mentioned earlier. However, the specific change of ${{{C}}_{{{m}{\dot \alpha }}}}{{\; + \;}}{{{C}}_{{{mq}}}}$ among individual damage models is different. At $\alpha$ = 0°, the values of ${{{C}}_{{{m}{\dot \alpha }}}}{{\; + \;}}{{{C}}_{{{mq}}}}$ decrease with the increase of damage extent. At $\alpha$ = 4°, the values of ${{{C}}_{{{m}{\dot \alpha }}}}{{\; + \;}}{{{C}}_{{{mq}}}}$ are basically the same for different damage extent. When AoA is $\alpha$ = 8°, the ${{{C}}_{{{m}{\dot \alpha }}}}{{\; + \;}}{{{C}}_{{{mq}}}}$ of the damaged model becomes larger than that of the undamaged model, and this phenomenon is more obvious at $\alpha$ = 12°. The combined dynamic derivative calculation indicates that the differences among damaged and undamaged models are mainly caused by the different missing areas of effective lifting surface before and after the moment reference point.

Figure 18. Combined dynamic derivative of pitch moment.

With the increase of AoA, the absolute values of the pitch combined dynamic derivatives first increases and then decreases with the increase of damage extent, which indicates the pitch damping characteristics don’t get worse in all cases, even slightly increase at small AoA. For the missile in this study, the dynamic performance is the worst at $\alpha$ = 4°. Compared with the undamaged model, the dynamic performance of the damaged models is even better when the AoA is within 5°. However, when the AoA exceeds 5°, it becomes worse rapidly. Generally speaking, at small angles of attack, the damage to one side of wing has little effect on the pitch dynamic performance of the missile.

4.2.2 Roll dynamic derivative

The hysteresis loops of the roll moment coefficient are presented in Fig. 19. After the wings being damaged, the hysteresis loop areas of their roll moment coefficients decrease, and the greater the damage extent, the smaller the hysteresis loops area, which is obvious at $\alpha$ = 0°. Different from the pitch direction, damage causes a great influence on the dynamic performances in the roll direction. The roll combined dynamic derivatives ${{{C}}_{{{mx\dot \alpha }}}}{{\; + \;}}{{{C}}_{{{mxq}}}}$ of several models at different AoA are shown in Fig. 20. For all models, the absolute values of the roll combined dynamic derivatives first decrease and then increase with the increase of the AoA. Compared to the undamaged model, the values of dynamic derivative of damaged models decrease a lot. When the damage is 0.5l, the absolute value of roll combined dynamic derivative decreases even more than 20%. Considering that the missile has X-shaped rudders and fins, this value is quite large. And the absolute value of the roll combined dynamic derivatives first decreases and then increases with the increase of the AoA, and this trend is the same as that of pitch combined dynamic derivatives. The turning point appears between $\alpha$ = 4° and $\alpha$ = 8°. According to Fig. 18, it is guessed that it may appear around $\alpha$ = 5°. However, different from the pitch direction, the roll performance of damaged models deteriorates at all angles of attack, which is also consistent with reality.

Figure 19. Hysteresis loops of roll moment coefficient at different AoA.

Figure 20. Combined dynamic derivative of roll moment.

4.2.3 Pitch-roll coupling dynamic derivative

The unsteady roll moment coefficient is shown in Fig. 21. The pitch-roll coupling moment caused by pitch oscillation is in the same order of magnitude as the roll moment, so it cannot be neglected. At $\alpha$ = 0°, the pitch-roll moments of all models differ by several orders of magnitude from other angles of attack, and it is nearly 0, so the hysteresis loop is not shown. When AoA is more than 0°, the force on the left and right sides of the missile is unbalanced in pitch oscillation, which in turn leads to roll moment. Meanwhile, the pressure centre is away from the centre of gravity in the longitudinal and transverse directions respectively. It leads to an increase in the areas of the hysteresis loops, which makes the dynamic coupling characteristics between the pitch and roll directions increasing apparently. With the increase of the initial AoA, the azimuths of hysteresis loops also decrease, and the pitch-roll hysteresis loops tend to be horizontal when the instantaneous AoA changes.

Figure 21. Unsteady roll moment coefficient caused by the pitch oscillation.

4.3 Damage limit of missile wing

After the unilateral wing of the missile is damaged, roll dispersion results from lift non-averaging produced by body-fixed vehicle asymmetries will appear [Reference Fuess38, Reference Crenshaw39]. Some aerodynamic performances can be saved by controlling the deflection of the rudders. The damage limit of the missile wing is judged since the missile does not roll. That is, the static balance of the missile after the wing damage is maintained by four rudder deflections. The full roll moment of the missile can be expressed as:

(8) \begin{align}{C_m} = {C_{{m_0}}} + {C_m}(\alpha ) \cdot \alpha + {C_m}(\beta ) \cdot \beta \end{align}

where, C _m is the full roll moment, C _m0 is the zero-liter moment, C _m ( $\alpha$ ) $\cdot$ $\alpha$ is the roll moment caused by the AoA $\alpha$ , and C _m ( $\beta$ ) $\cdot$ $\beta$ is the trim moment caused by the rudder deflection angle $\beta$ .

Figure 22 shows the roll moment provided by one rudder at different deflection angles. In the figure, $\beta$ represents the deflection angle of a single rudder. It can be seen from the figure that with the increase of the AoA, the static roll moment provided by the smaller rudder deflection angles first increase and then decrease. The roll moments provided by the larger rudder deflection angles decrease continuously.

Figure 22. Roll moment coefficient caused by one rudder deflection.

The deflection efficiency of the remaining three rudders is regarded as the superposition of one rudder. Through this simplified processing, the roll moment provided by the four rudders at different angles of attack can be obtained. The maximum deflection angle of each rudder is limited to 15°. Figure 23 shows the static roll moments that can be provided when all four rudders are deflected to 15°, as well as the roll moment required to maintain balance for four models with different extents of damage at different angles of attack.

Figure 23. The static roll moments provided by four rudders, and the roll moments required of different damaged models.

As can be seen in Fig. 23, even at a maximum angle-of-attack of $\alpha$ = 12°, all damaged models in this study can achieve roll moment balancing. Next, the damage limit of missile wing at $\alpha$ = 12° will be studied, which is the worst flight condition. According to the relationship between the change of roll moment coefficient and the damage extent of missile wing, it is found that it basically conforms to the linear change law: C _mx = −0.2941(h/l). When the four rudders are all deflected to 15°, the maximum roll moment available is C _mxr = −0.1829. Then, the maximum allowable wing damage extent is h/l = 0.625 when the missile does not roll.

Finally, the correctness of the damage limit needs to be validated by numerical simulation. The results are shown in Table 7. C _mxt represents the roll moment coefficient of the missile with right wing h/l = 0.625 damage trimmed by rudders deflection. The relative error C _mxt /C _mxr is −2.023%. This is a relatively small value. This shows that the limit damage extent of missile wing is h/l = 0.625 with high credibility.

Table 7. Trim validation of the h/l = 0.625 damaged model at $\alpha$ = 12°