3.1.1 Rotational stiffness K_θx around the folding axis

The schematic of the hinged–locked structure showed that contributions to the connection stiffness mainly consists of the following parts: (1) deformation caused by interfacial contacts of pins and holes and (2) bending deformation of the locking pin, as demonstrated in Fig. 5 [Reference Ren, Wang, Wang, Zhou, Chang, Cai and Lei17]. High deformation in each part corresponds to low rotational stiffness K _θx when the external load ΔM _x is constant. Therefore, each part of the deformation should be calculated separately and then combined to obtain the angle of outer fin and the connection stiffness.

i) Contact deformation

Clearances are reserved during the structural design to ensure the satisfactory operation state of the locking mechanism and avoid the pin from becoming stuck in the nonworking position of the hole. The preset clearance value between the locking pin and the hole wall of the folding mechanism is 1 mm. The vehicle is for one-shot use, thus, without considering the variation of clearance caused by friction and wear. Shape of the locking pin and the hole wall are sufficiently dissimilar, and the interface is nonconforming, thereby implying that the area of contact is remarkably smaller than characteristic dimensions of contacting bodies. The Hertz contact model [Reference Roark and Young18] can be used to calculate the contact deformation as follows:

(1) \begin{align}\frac{1}{{{E_{{\rm{eq}}}}}} = \frac{1}{2}\!\left( {\frac{{1 - v_1^2}}{{{E_1}}} + \frac{{1 - v_2^2}}{{{E_2}}}} \right)\end{align}

(2) \begin{align}\frac{1}{{{R_{{\rm{eq}}}}}} = \frac{1}{{{R_1}}} + \frac{1}{{{R_2}}}\end{align}

(3) \begin{align}a = \sqrt {\frac{{8P{R_{{\rm{eq}}}}}}{{\pi {E_{{\rm{eq}}}}}}} \end{align}

(4) \begin{align}\delta (t) = \frac{{4P\!\left( {1 - v_1^2} \right)}}{{\pi {E_{{\rm{eq}}}}}}\!\left( {\frac{1}{3} + \ln \frac{{4{R_{{\rm{eq}}}}}}{a}} \right)\end{align}

where E _eq and R _eq are the equivalent elastic modulus and equivalent radius, respectively; E ₁ and E ₂ are the elastic modulus of the curvature of the cylinder and plane, respectively; υ ₁ and υ ₂ are the Poisson’s ratio of the curvature of the cylinder and plane, respectively; R ₁ and R ₂ are the radius of the curvature of the cylinder and plane, respectively; a is the half-contact width; δ is the contact deformation; t is the time; and P is the vertical load.

3.1.2 Connection stiffness in other directions

The connection stiffness model presented in Fig. 6 is established to solve the translation stiffness values K _x , K _y and K _z as well as the rotational stiffness values K _θy and K _θz . The tensile load carried by the outer fin in the model is transferred through the cooperation of the shaft and hole. The quasi-rigidity assumption is adopted for the body and only the resistance of connections in contact areas is considered while analysing the connection stiffness because bodies provide no or very minimal additional stiffness to the system. The contact stiffness between the shaft and the outer fin hole are denoted k ₁ and k ₂, respectively, and the contact stiffness between shaft and the inner fin hole is denoted k ₃. The overall stiffness of the hinged–locked structure is then represented by a combination of these three parts, k ₁ – k ₃ as the outcome of the analysis.

Figure 6. Contact model of the shaft and hole connection.

The ESDU-78035 model [Reference Koshy, Flores and Lankarani20] is used to solve the contact stiffness k ₁ – k ₃ because the clearance in the rotation pair is small. The ESDU-78035 model is expressed as follows:

(5) \begin{align}\delta (t) = {F_N}\!\left( {\frac{{{\sigma _i} + {\sigma _j}}}{L}} \right)\left[ {\ln\! \left( {\frac{{4L\left( {{R_i} - {R_j}} \right)}}{{{F_N}\left( {{\sigma _i} + {\sigma _j}} \right)}}} \right) - 1} \right]\end{align}

where σ _i(j) = (1−υ _i(j) ²)/E _i(j); υ _i and υ _j are the Poisson’s ratio of materials of the hole and shaft, respectively; E _i and E _j are the elastic modulus of materials of the hole and shaft, respectively; R _i and R _j are the outer hole and inner shaft radii, respectively; F _N is the total normal load in the contact area; L is the axial length of the contact zone; and t is the time. Dimensions and material properties of different subcomponents along the flight time are used into Equation 5 to obtain the contact stiffness k ₁ – k ₃.

Figure 7. FEM the deployable fin and diagram of the zone partition.

The structural connection stiffness is then calculated under the following assumptions:

1. 1) The translation stiffness K _y and K _z are treated the same, in which springs of k ₁ and k ₂ are parallel and connected in series with k ₃. Thus, the translation connection stiffness K _y (or K _z ) is expressed as follows:

(6) \begin{align}{K_y}\left( {or\,\,{\rm{ }}{K_z}} \right) = \frac{{\left( {{k_1} + {k_2}} \right){k_3}}}{{{k_1} + {k_2} + {k_3}}}\end{align}

1. 2) Assuming that the rotational stiffness K _θy and K _θz are the same and contributions of k ₁ and k ₂ to the torque are considered, then the rotational connection stiffness K _θy (or K _θz ) is expressed as follows:

(7) \begin{align}{K_{\theta y}}\left( {or\,\,{\rm{ }}{K_{\theta z}}} \right) = {k_1}{d_1} + {k_2}{d_2}\end{align}

where d ₁ (or d ₂) is the distance between the spring k ₁ (or k ₂) and the symmetrical plane of the inner fin.

1. 3) Determination of the translation stiffness K _x is indirectly and dominated by the interface between the inner and outer fins. The influence of the body stiffness is significantly greater than that of the contact stiffness due to the large contact area; hence, the contact stiffness K _x is set to 1 × 10¹⁰ N/m as a relatively large value.

3.2.1 FEM

Finite element analysis is conducted using ABAQUS/Standard to model the complex structure of the inner and outer fins, and quadratic tetrahedral thermal displacement coupling elements (C3D10MT) are assigned (Fig. 7). To verify the mesh irrelevance, a comparison between extremely coarse, coarser, fine and extremely fine is studied based on the variation of temperature of locking pin region versus time. The number of elements in each case is mentioned in Table 1. The variation of temperature versus time for different grid sizes is shown in Fig. 8. The fine mesh with 129,561 elements is chosen based on the study. The metal skeleton and composite heatshield of the outer fin model are established separately and then tied together regardless of the interface layer. Thus, displacement and temperature are considered continuous. Nodes near contact areas of inner and outer fins are rigidly constrained to the control point at the centre of the rotating shaft through MPC. A bush element is established between control points to simulate the connection relationship, and the nonlinear stiffness obtained using previous analytical methods is assigned to the bush element.

Table 1. Number of elements in each case

Table 2. Properties of superalloy steel

Figure 8. Variation of temperature of locking pin region versus time by using different meshes.

The following material properties are assigned to different subcomponents: superalloy steel for the inner fin and skeleton, hardened steel for pins and shaft, and high-performance composite material for heatshield. Degraded material properties are used in the model because of the evident rise of temperature in the entire fin. Tables 2–4 present material properties at elevated temperatures obtained through experimental measurement. The density of the superalloy steel, hardened steel and composite material is set to a constant value of 8,240, 7,780, and 1860 kg/m³, respectively, given that high temperature exerts a negligible influence on the behaviour of density. In addition, mechanical properties of the high-performance composite material in this study are basically unchanged in the analysis temperature range; hence, the constant value of elastic modulus, Poisson’s ratio and thermal expansion coefficient are utilised and emissivity is set to 0.8.

Table 3. Properties of hardened steel

Table 4. Properties of composite thermal insulation materials

3.2.2 Boundary conditions

The root of the all-moveable fin axis of a real aircraft is generally connected to the actuation and transmission system, such as rods, rocker arms and electric motor. The detailed simulation of these complex connections is beyond the scope of this study. Therefore, the simplified fixed displacement boundary condition is utilised through multipoint constraint (MPC) method, which is similar to the resolvent work in Reference [Reference Koshy, Flores and Lankarani20].

Determining the unsteady heat flux distribution is a fundamental task of heat transfer analysis. Different heat flux and radiation values are applied as input thermal loads or boundaries for various surface zones in FEM to simulate flight heating conditions reasonably. Specifically, the windward and leeward surfaces of the fin are denoted A and B, respectively, and then further divided into wedge (denoted 1) and plane (denoted 2) areas. Therefore, a total of four partitions (A1, A2, B1 and B2) are set, as shown in Fig. 7.

It should be noted that, for the wing or fin of high-speed vehicle, thermal protection structures of the leading-edge region, wedge region and planar region are usually designed separately. Although the distribution of heat flux is not uniform in each region, the same thickness of thermal protection structure is used through each region because of the processing characteristics. Generally, the location with highest heat flow in each region is utilised, ensuring that the designed thickness of heat protection system can meet the requirements along the whole flight profile. Selecting more points and dividing more regions can achieve better performance and obtain detailed thermal evaluation results. However, using existing conservative design methods, that is, only the highest heat flow point in each region is considered during the preliminary design process, is efficient and sufficient in engineering applications. The superiorities and disadvantages of the similar design methodology are specifically discussed in the literature [Reference D’Souza, McGuire, Torricelli, Visser and Hays21].

3.3.1 Cold-wall heat flux

Accurately predicting aerodynamic heating flux is the critically important task for the high-speed vehicle design. Generally, based on whether the shape is regular, the surface of the whole vehicle can be divided into large areas and local areas. Large areas include conical and cylindrical parts of aircraft, and the heat flux can be obtained by utilising engineering algorithms conveniently. Furthermore, the estimating formula can be combined with the trajectory to calculate the heat flow at different flight times. The local structural area includes the wing or fin structures. Due to the complex flow fields of local structures, it is often difficult to accurately predict heat flow using empirical formulas. Although wind tunnel test or CFD simulation can be used to obtain the heat flux conditions, obtaining the heat flow through the entire flight process faces the obstruction of computational complexity, time alocation and poor economy.

In this paper, a method (named interference factor method) that combines CFD numerical simulation and engineering algorithms is conducted. A more detailed description of the interference factor analysis method can be found in reference [Reference Liu, Tu, Sun, Tan and Yang22]. The analysis process and main calculation formulations are briefly listed as follows:

1. 1) Select a reference point in a large area, and select a point in a local surface, such as different region of fin. Define reference factors χ:

   (8) \begin{align}\chi = {\raise1.5ex\hbox{${{Q_f}}$}\!\! \bigg/\!\! \lower1.5ex\hbox{${{Q_b}}$}}\end{align}

Where: Q _f is the heat flux density of the fin, and Q _b is the heat flux of the reference point.

1. 2) By using STAR-CCM+ software [Reference Cinosi, Walker, Bluck and Issa23], reference factor interpolation tables under different flight conditions, such as flight altitude, Mach number, angle-of-attack and fin yaw can be obtained.

2. 3) The reference enthalpy method [Reference Huang, Wang and Ma24] is used to calculate the heat flux Q _b of the reference point through the whole flight trajectory. The formula is as follows:

Calculate the oblique wave angle β:

(9) \begin{align}\sin \beta = \sin {\theta _c}{\left( {\frac{{\gamma + 1}}{2} + \frac{1}{{M_\infty ^2{{\sin }^2}{\theta _c}}}} \right)^{\frac{1}{2}}}\end{align}

Where, θ _c is the half cone angle.

Calculate the physical parametres behind the shock wave:

(10) \begin{align}{H_e} = \frac{{\left( {7M_\infty ^2{{\sin }^2}\beta - 1} \right)\left( {M_\infty ^2{{\sin }^2}\beta + 5} \right)}}{{36 \cdot {{\left( {{M_\infty }\sin \beta } \right)}^2}}}{H_\infty }\end{align}

(11) \begin{align}{M_e} = \sqrt {5\left( {\frac{{{T_\infty }\left( {1 + 0.2M_\infty ^2} \right)}}{{{T_e}}} - 1} \right)} \end{align}

(12) \begin{align}{V_e} = {M_e} \cdot 20.0467 \cdot \sqrt {{T_e}} \end{align}

Where, H _e , M _e , and V _e are the static enthalpy, Mach number and velocity behind the oblique wave, respectively.

The Sutherland method is used to obtain viscosity coefficient on the outer edge of boundary layer:

(13) \begin{align}\mu _e^* = 1.458 \times {10^{ - 6}}\frac{{{{\left( {H_e^*} \right)}^{1.5}}}}{{H_e^* + 110.4}}\end{align}

Calculate the reference enthalpy based on the Eckert reference enthalpy method:

(14) \begin{align}H_e^* = 0.28{H_e} + 0.5{H_w} + 0.22{H_{re}}\end{align}

Where, H _w is the wall temperature enthalpy, H _re is the recovery enthalpy, and H _e is the boundary layer outer edge enthalpy.

Calculate the heat flux density at the reference point Q _b

(15) \begin{align}{Q_b} = 0.332\sqrt 3\Pr {^{ - \frac{2}{3}}}{\left( {\frac{{{\rho ^*}{\mu ^*}}}{{{\rho _e}{\mu _e}}}} \right)^{ - 0.5}} \times {\rho _e}{\mu _e}{\left( {0.5{R_{ex}}} \right)^{ - 0.5}} \times \left( {{H_r} - {H_w}} \right)\end{align}

1. 4) Predict the heat flux Q _f of the fin surface using the flowing equation

(16) \begin{align}{Q_f} = \chi \cdot {Q_b}\end{align}

3.3.2 Hot-wall heat flux

Because heat conduction is mainly along the thickness direction, one-dimensional heat transfer model can be used to simplify the analysis. The one-dimensional analysis program is integrated into cold-wall heat flow simulation program, and the outer surface temperature of the structure is iteratively calculated. Thus, the hot-wall heat flux can be obtained.

The one-dimensional heat conduction equation is shown as follows:

(17) \begin{align}{\rho _c}{C_{pc}}\frac{{\partial T}}{{\partial t}} = \frac{\partial }{{\partial y}}\left( {\lambda \frac{{\partial T}}{{\partial y}}} \right)\end{align}

Where, y is the coordinate along the thickness direction, ${\rho _c}$ is the density of the material, ${C_{pc}}$ is the specific heat, and T is the material temperature, and λ is the thermal conductivity.

For the metal-composite integrated structure, the governing equation is:

(18) \begin{align}{\rho _{c,i}}{C_{pc,i}}\frac{{\partial T}}{{\partial t}} = \frac{\partial }{{\partial y}}\left( {{\lambda _i}\frac{{\partial T}}{{\partial y}}} \right)\end{align}

Where, the subscript i represents the material parametres of the i th layer.

For the outer surface, the thermal radiation is considered, and heat transfer balance equation is:

(19) \begin{align}Q - \varepsilon \sigma T_w^4 = - \lambda {\left( {\frac{{\partial T}}{{\partial y}}} \right)_{y = 0}}\end{align}

The adiabatic condition is used for inner surface, and the equation is:

(20) \begin{align}0 = - \lambda {\left( {\frac{{\partial T}}{{\partial y}}} \right)_{y = \delta }}\end{align}

The heat flow equilibrium and temperature conditions of inter-surfaces satisfied:

(21) \begin{align} - {\lambda _1}\left( {\frac{{\partial {T_1}}}{{\partial y}}} \right) = - {\lambda _2}\left( {\frac{{\partial {T_2}}}{{\partial y}}} \right)\end{align}

(22) \begin{align}{T_1} = {T_2}\end{align}

5.1.1 Flight profile

The heat flux is closely related to flight speed, altitude, aircraft attitude and aircraft shape. The high-speed aircraft is regard as a particle, and a 3-degree of freedom aircraft dynamics model (Equation (23)) is used to solve the flight profile considering gravity, thrust and aerodynamic parametres. And then, the heat flux can be obtained by using the interference factor method as mentioned above.

(23) \begin{align}m\frac{{dv}}{{dt}} &= F\cos \alpha - D - mg\sin \theta \nonumber\\mv\frac{{d\theta }}{{dt}} &= F\sin \alpha + L - mg\cos \theta \nonumber\\\frac{{dx}}{{dt}} &= v\cos \theta \\\frac{{dy}}{{dt}} &= v\sin \theta \nonumber\\\frac{{dm}}{{dt}} &= - {m_f}\nonumber\end{align}

Where, v is the flight speed, α is the angle-of-attack, F is the thrust, m _f is the fuel consumption rate, m is the vehicle mass, θ is the trajectory inclination angle, L is the lift force and D is the drag force.

The angle-of-attack is designed according to the desired goal. The entire flight process of the trajectory can be divided into three stages:

1. 1) After the vehicle is launched, the engine ignites, and the missile undergoes a climb process in which the fuel burns out. Thereafter, the trajectory is a parabolic trajectory similar to a ballistic missile.

2. 2) After re-entering the atmosphere, the projectile does not directly dive to hit the target. The vehicle is pulled up again, achieving a jumping manoeuver. Relying on the lift force without power, the second parabola is completed.

3. 3) Due to the fact that the trajectory undergoes a second takeoff, the speed of the projectile during reentry is somewhat lower than that of the first, and the trajectory tends to flatten under the effect of dense atmospheric lift. Therefore, it is necessary to make the projectile bow and dive after the second reentry.

The variation equation of AOA is:

(24) \begin{align}\alpha \left( t \right) &= - 4{\alpha _{\max ,j}}Z\left( {1 - Z} \right) \nonumber\\[6pt] Z &= {e^{ - {a_j}\left( {t - {t_j}} \right)}},\,j = 1,2,3\end{align}

Where, subscript j represents the number of jumps; α _max,j is the maximum allowable angle-of-attack for the j th jump; and the design value a _i can adjust the speed of the turn; t _i is the jump start time.

The calculated flight profile, AOA, cold-wall heat flux and hot-wall heat flux are shown in the following Fig. 13.

Figure 13. Flight profile and heat flux results during the flight.

Results of heat flux (Figs13(c) and 13(d)) show that: in the initial stage, the heat flux reaches the first peak due to the low altitude and the large AOA. Then the heat flux is relatively small due to the thin air at a higher altitude. When the vehicle performs the second pull-up manoeuver with large AOA, the heat flux increases obviously again. However, when the vehicle performs the final dive, although the altitude is relatively low, and the angle-of-attack attitude is relatively large, the velocity of the end of the flight is relatively small. The heat flux is therefore not as obvious as the first two peaks.

5.1.2 Response of thermal displacement coupling

The heat flux of A1, A2 and B1 is shown in Fig.13(d), and the flux of B2 is 75% of that of B1. The temperature variation of the metal skeleton near the locking pin region versus time is shown in Fig. 14. The temperature rises faster in the initial state due to the peak heat flux and decreases slightly because of the high-temperature outer surface’s thermal radiation at the end of the flight. Two nodes near the locking hole are selected, and the distance change (Fig. 15) in the upper and lower surfaces of the locking hole obtained by subtracting the displacement of the two nodes reveals the clearance variation.

Figure 14. Temperature of metal skeleton near the hole of locking pin versus time.

Figure 15. Clearance change of the locking hole versus time.

5.2.1 Rotational stiffness K_θx

The clearance and temperature obtained by previous analyses are assigned into the semianalytical model of the hinged–locked structure.

First, the reference temperature of the contact zone is determined via the average temperature of the node. Temperature-dependent material properties are used in the Hertz model (Equation 1) to calculate contact deformation.

Second, the temperature rise of the whole locking pin is assumed unclear. Therefore, properties at ambient temperature are used in the beam model to obtain the bending deformation of pin while excluding thermal stress.

Finally, the obtained K _θx values are illustrated in Fig. 16.

The (i) clearance variation due to thermal expansion and (ii) material property degradation, which considers only one factor, are compared to determine the influence of two key factors on the connection stiffness (Fig. 16). The connection stiffness considering both factors K _v, the result of considering the material factor alone K _vm, the result of considering the clearance factor alone K _vc, the reference value K _v, relative error of K _vm (Re _vm) and relative error of K _vc (Re _vc) are listed in Table 8.

Table 8. Values of connection stiffness K_θx at typical moments

Figure 16. Connection stiffness K _θx versus time.

The results in Fig. 16 and Table 8 present the following:

1. 1) The connection stiffness changes nonlinearly with time and shows an inverse correlation with the trend of temperature. The initial state of K _v (1.278 × 10⁵ Nm/rad) is considered the reference, and the value at the end time of heating is 1.233 × 10⁵ Nm/rad, which is 3.52% lower than the initial value.

2. 2) The final K _v is 1.270 × 10⁵ Nm/rad when only the clearance effect is considered. This finding showed a decrease of only 0.59% compared with the initial value. The effect of clearance is insignificant in this analysis because the corresponding expansion of the metal hole at the end is only 0.245 mm, which is smaller than the initial preset clearance of 1 mm.

3. 3) Table 7 also shows that Re _vm and Re _vc are less than 0.73% and 3.41%, respectively, during the heating process. For the connection model with a preset clearance ensures the normal operation of the mechanism, material property degradation due to high temperature is the main cause of the nonlinear variation of the connection stiffness, and the effect of the clearance change caused by thermal expansion is insignificant.

5.2.2 Connection stiffness in other directions

Translational (K _y and K _z ) and rotational (K _θy and K _θz ) stiffnesses are mainly determined by contact properties between the shaft and holes. The temperature of the wall surface and the clearance between the hole and the shaft are obtained and then assigned into Equation (5) to solve the connection stiffness.

Figure 17 shows the variation of K _y , K _z , K _θy , and K _θz and the downward trend in each curve because of the material property degradation. The stiffness of the translational connection is 2.288 × 10⁸ and 2.117 × 10⁸ N/m at the initial and final moments, respectively. If the initial value is used as a reference value, then the final translational stiffness decreases by 7.47%. In addition, the rotational stiffness is approximately two orders of magnitude smaller than the translational stiffness, as determined by d ₁ and d ₂.

Figure 17. Connection stiffness versus time (K _y , K _z , K _θy and K _θz ).

5.3.1 Thermal mode

Numerical test cases are developed to investigate the contribution of connection stiffness and material property degradation to the natural frequency of the fin. Two cases are notable, and material property degradation is considered in both cases. The effect of connection stiffness and temperature-dependent material properties are preserved in the first case because the time-varying connection stiffness is used in model. Meanwhile, the constant connection stiffness is used in the second case of interest during the whole heating process (the connection stiffness is maintained equal to the initial value), and the material effect of the fin body is included as the only variable parametre. Figure 19 shows the variation of the first three natural frequencies of the deployable fin. Values at typical moments, including the frequency under variable connection stiffness f _v, frequency under constant connection stiffness f _c and relative error between two stiffness model Re, are listed in Table 9.

Table 9. First three modal frequencies of the deployable fin at typical moments

The results in Fig. 19 and Table 9 indicate the following:

1. 1) Natural frequencies of each order show a downward trend during the heating process. The first three orders of natural frequencies corresponding to the thermal mode are 1.67, 7.75 and 16.28 Hz lower than that of the ambient temperature for the time-varying stiffness model. The first-order mode shape of the deployable fin is an approximate rigid body mode shape with constraints and the second- and third-order mode shapes are self-elastic deformation.

2. 2) The frequency of the time-varying stiffness model is significantly lower than that of the fixed-value spring model because the connection stiffness of the spring decreases at high temperatures.

3. 3) The first-order mode is rigid, the connection stiffness plays a dominant role, the third-order frequency is the elastic deformation of the fin surface and the connection stiffness exerts a marginal effect on the deformation area far from the connection area. Hence, the relative error of the first frequency between the fixed and time-varying stiffness models is 1.46%, while the relative error of the third frequency is only 0.46% at the end of time.

5.3.2 Contribution rate of connection stiffness to frequency

The fixed connection stiffness results only reflect the material property degradation of the fin body, while the time-varying connection stiffness results comprehensively reflect the combined effect of the weakening of the connection stiffness and the variation of material properties. Frequency changes Δf _v and Δf _c relative to their respective initial frequencies are calculated, and the contribution rate index of the connection stiffness change to the frequency is defined as follows:

(25) \begin{align}C = \frac{{\Delta {f_{\rm{v}}} - \Delta {f_{\rm{c}}}}}{{\Delta {f_{\rm{v}}}}} \times 100\% \end{align}

where a large contribution rate index indicates a strong influence of the connection stiffness.

Figure 20 shows the index C at different times. The contribution rate of the connection stiffness to the first-order frequency fluctuates slightly but is maintained at 50%–60% as the heating progresses. The contribution rate to the second- and third-order frequencies is about 30% and 10%, respectively. The results revealed that the change of connection stiffness affects the first-order frequency the most evidently, followed by the second- and third-order frequencies.

Figure 20. Contribution rate of the connection stiffness on frequency at different times.