Abstract
The inflatable reentry vehicle provides a new technical way in aerospace entry, descent, and landing. The structural failure of inflatable reentry vehicle experiment caused by thermal aeroelastic effect is serious, which needs to be further studied. A traditional numerical method about flexible vehicles separates the aeroheating and aeroelastic problems, resulting in poor matching with the actual test. In this paper, a thermal-fluid-solid coupling model considering inflation gas effect was established, which associates the aeroheating and aeroelastic modules and adopts the LES to improve the depicting ability of hypersonic flow. The model was used to solve the thermal aeroelastic characteristics under extreme aeroheating load. From aeroheating results, the large-scale vortex on windward generated by the interaction of the shock layer and boundary layer has great influence on aeroheating due to the heat dissipation, and the skin deformation also increases the surface friction and local heating near depressions. From aeroelastic analysis, the flexible structure performs violent forced vibration induced by the unsteady large-scale vortex on windward, and the aeroheating effect will significantly increase the thermal stress and natural vibration properties. The thermal-fluid-solid coupling method for the flexible structure proposed in this paper provides a reasonable reference for engineering.
1. Introduction
The traditional rigid aeroshell is not competent in complex return missions due to the large volume, low payload, and inconvenient operation. Aimed at these defects, the inflatable reentry vehicle experiment (IRVE) is becoming an emerging technology [1, 2]. However, the safety problems of the thermal protection system (TPS) and flexible inflation toroids under hypersonic flow are serious [3, 4]. The aerodynamic force will make large deformation and dynamic stress of flexible film [5], and the aeroheating will not only increase the surface temperature of TPS but also induce the thermal expansion of inflation gas [6]. For instance, the IRDT-2 and IRDT-2R launched by Russia were ablated under violent aeroheating [7, 8]; momentary structural instability appears on the inflatable reducer launched by JAXA due to the pressure difference between inflation gas and aerodynamic force [9]. Note that the large deformation of flexible film will in turn change the aerodynamic loads; thus, it is a complex thermal aeroelastic system [10].
Traditional thermal-fluid-solid coupling problems are mainly about conventional elastomers such as airfoils and turbine blades, but it becomes more complicated for flexible inflatable structures. Due to the unique properties of flexible materials, the structural dynamics shows strong geometric and material nonlinear features, which bring difficulty to the coupling of structural dynamics, hypersonic flow, and thermal dynamics [11]. In the past research about the IRVE system, the thermal aeroelasticity was divided into the aeroheating and aeroelastic problems based on the premise that the characteristic time of the flow field is much smaller than heat transfer [12]. Furthermore, these two problems were usually separated to, respectively, solve the surface temperature and structural dynamics under aerodynamic loads, and the one-way coupling method was most widely used. For example, McGuire et al. established the three-dimensional heat conduction model of TPS and studied the temperature response after aeroheating solution [13]; Xu researched the dynamic response of IRVE and checked the structural stress after obtaining the steady aerodynamic force [14]. For the aeroheating problem, the one-way coupling of aeroheating and heat transfer is proved to be in good agreement with the high-enthalpy wind tunnel test [15]. But for the aeroelastic problem, the large deformation effect of flexible film does have great influence on aerodynamic loads [16]; thus, some scholars attempted to develop the two-way fluid-solid coupling models about flexible structure. Wei et al. established the fluid-solid coupling model of folded membrane boom by the ALE method and solved the inflatable deployment process [17], and Karagiozis proposed an efficient iterative method to reveal the dynamic response of inflatable reducer under Ma 2 [18].
However, the existing coupling methods have some defects in describing the structural dynamics under extreme aeroheating loads, resulting in poor matching between numerical results and the actual flight test. Actually, the interaction between aeroheating and aeroelasticity is significant: the surface temperature and inflation gas expansion induced by aeroheating cannot be ignored in structural dynamics solution, and the large deformation effect of flexible film also has notable influence on aeroheating in turn. These effects have been already confirmed in some tests recently [19, 20]; thus, the independent solutions of aeroheating and aeroelastic problems are inappropriate. Besides, it is still not clear how the inflation gas participates in the coupling process in existing research, and the widely used Reynolds-average equation also has difficulty in depicting the details of hypersonic flow with strong disturbance. Therefore, it is vital to improve the multiple field coupling system and enhance the simulation capability of each field.
In this paper, a thermal-fluid-solid coupling model considering the thermal effect of inflation gas is established, which associates the aeroheating and aeroelastic problems, so as to reveal the actual physical process more closely than separated aeroheating and aeroelastic methods. Meanwhile, the ground vibration test is carried out to verify and correct the structural dynamics model, and the large eddy simulation (LES) is adopted to depict the details of hypersonic flow. Based on the model, the thermal aeroelastic characteristics under extreme aeroheating are investigated, and the mechanisms of external flow and inflation gas on aeroheating and structural dynamics are explored.
2. Methods and Numerical Models
2.1. Thermal-Fluid-Solid Coupling Scheme
In this paper, the aeroheating and aeroelastic problems are associated based on ANSYS software, as Figure 1 shows. The one-way coupling is adopted to solve the temperature of TPS and thermal expansion of inflation gas, and the two-way fluid-solid iterative coupling is used to solve the aeroelasticity. The coupling is specifically achieved by the time-marching method: in every physical time step, the transient CFD and structural dynamics, respectively, transmit the aerodynamic force and structural deformation to each other. The iteration is continuously carried out until the residual is within the specified range. In this process, the surface temperature and inflation pressure are set as the thermal boundary of solid, and the aeroheating obtained in the last time step is extracted as the input of heat conduction.
2.2. Research Object and Finite Element Model
The IRVE-3 aircraft is adopted in this study, and Figure 2 presents the schematic diagram [21]. Seven toroids inflate independently to provide structural stiffness, the toroids are connected with each other through six interlayers to increase the structural strength, and the skin is wrapped outside the toroids to form the thermal protection system (TPS). The maximum diameter of IRVE is 3 m, and the diameter of each inflation toroid is 0.22 m. The Kevlar fiber membrane is selected as the material [21], which is a kind of isotropic material.
(a)
(b)
The membrane static equation and vibration equation considering the effect of inflation pressure are deduced to establish the finite element model. A uniform film with a thickness is stretched at the horizontal boundary, and a pressure is exerted on it. At this point, the film will produce lateral displacement and the internal tension and due to deformation. is the function of and , and its governing equation is given by where is the modulus of elasticity. Set to be the mass per unit area, the differential equation for the free vibration is
The fiber membrane is a highly nonlinear material that has no bending resistance. It will produce strong large deformation and small-strain effects under the action of load; thus, the geometric nonlinearity must be considered in the study. Figure 3(a) shows the finite element model with a total grid of 100 thousand, in which the shell 41 element is applied. The shell 41 element is a kind of nonlinear element that can only bear tensile stress, and its stiffness is generated by the internal stress, which conforms to the characteristics of flexible film. Due to the circumferential symmetry of the structure, the swept method is used to divide the tetrahedral mesh to ensure orthogonality. The study on grid independence is shown in Figure 3(b), in which the maximum stress and natural frequency are chosen as the monitoring parameters.
(a)
(b)
The modal vibration test based on the hammering method is carried out to validate the structural dynamics method; the test object is a simplified IRVE system in the form of three-layer toroids (see Figure 4). The time-domain response is measured by acceleration sensors after applying an impulse excitation signal to the structure with a hammer, and the modal shape and frequency can be solved by modal parameter identification. In this test, the object is supported on the tire to simulate a free state, and there are ten measuring points on each layer of the toroid to measure the longitudinal stretching and torsion shapes. From the results, the first- and second-order frequencies vary little with the inflation pressure, and the deviations between numerical values and tests are about 9%. The third-order frequency varies approximately linear with the inflation pressure, and the deviation is about 7%. Considering that the pin, inflation valve, and other connecting devices will produce added mass and thereby reduce the natural frequency, thus the deviations are within the reasonable range.
(a)
(b)
(c)
2.3. CFD Model
These hypotheses are employed to simplify the CFD: (1) the boundary conditions remain constant in the simulation time and (2) the initial attack angle and sideslip angle are zero. Besides, the governing equations are described as the following: (1)Continuity equation(2)Momentum equation(3)Energy equation
Besides, the convective heat transfer between the high-enthalpy flow and the structure, as well as the radiation generated by the wall, should be considered when solving the aerodynamic heating. Set the blackbody radiation coefficient to be , and the emissivity is equal to 0.89; then, the thermal radiation equation can be stated as
The large eddy simulation (LES) based on the finite volume method is applied in CFD simulation. The WALE (Wall-Adapting Local Eddy Viscosity Model) is chosen as the SGS (subgrid scale) model, which ensures the simulation accuracy near the wall by considering the effects of momentum transfer of turbulent, and the ideal subgrid eddy viscosity distribution can be solved without adding a wall function or dynamically adjusting the model constant. The eddy viscosity coefficient in WALE is where is a model parameter and the value is generally 0.325. is the strain rate tensor, and in order to reflect the effect of the symmetric part and the antisymmetric part of the velocity gradient tensor, it is constructed as where is the generated term of buoyancy. The validation of the used WALE model and its comparison with other SGS models are carried out based on the double ellipsoid configuration. The double ellipsoid model is a typical hypersonic vehicle, where the three-dimensional bow shock, secondary shock, and three-dimensional separation flow exist around the vehicle. Many wind tunnel tests have been carried out and verified with each other, which can be used as a standard verification case of the hypersonic algorithm [22–24]. Figure 5 presents the double ellipsoid model and comparison on heat flux of different SGS models. The results show that obvious excessive dissipation exists in the WMLES (Algebraic Wall-Modeled LES Model) and SM (Smagorinsky) models, causing a large deviation in prediction. The simulation on dissipation of the subgrid scale by the WALE and KET (Dynamic Kinetic Energy Subgrid-Scale Model) is accurate before and after separation; thus, the heat flux predicted by these two models is in good agreement with the experiment.
(a)
(b)
(c)
The high-resolution scheme is adopted for spatial discretization, and the second-order Euler scheme is adopted for time discretization to ensure the conservation of time. Figures 6(a) and 6(b) present the sectional view of the fluid calculation domain, where the hybrid grids of both structured and unstructured are applied, and the total mesh number is 8 million. The mesh height of the first layer on the surface is 0.25 mm to ensure that the value of is approximately equal to 1. The total calculation time is 2 s, and the time step is 0.0002 s. In each time step, there are 4 virtual time steps of CSD, and in each CSD step, there are also 5 virtual time steps of CFD. The grid independence study is also carried out to validate the rationality of the grid number, as shown in Figure 6(c).
(a)
(b)
(c)
2.4. Aeroheating Method
The TPS system is divided into three parts from outside to inside: thermal protection layer, thermal insulation layer, and airtight layer [25], as Figure 7(a) shows. The aeroheating transfers through three function layers in turn and finally causes the thermal expansion of inflation gas. The heat conduction differential equation is expressed as follows:
(a)
(b)
In Equation (10), is temperature, is specific heat capacity, is thermal conductivity, and is the internal heat source. In respect of the process of thermal expansion, the average temperature of the innermost layer is extracted as the temperature of inflation gas, and then, the inflation pressure is solved by the state equation of ideal gas. Figure 7(b) gives the validation of the heat conduction model based on Ref. [25]; the temperature curve of the innermost layer is consistent with the test.
2.5. Boundary Conditions
In the previous research by the authors in Ref. [26], the heat flux during reentry is preliminarily solved based on flight dynamics and engineering algorithm. The heat flux increases first and then decreases, reaching the maximum at 60 km. In order to obtain the general laws of thermal aeroelasticity, eight conditions (points to ) from the hypersonic region are chosen for research, in which point corresponds to the extreme aeroheating condition (see Table 1). Their Mach numbers change from 5 to 12, and pressures reduce from 400 Pa to 24 Pa, representing different stages of hypersonic flight. According to Ref. [27], the Knudsen number () is less than 0.01 when the height is below 91 km, which indicates that it can be regarded as continuous fluid. In the hypersonic CFD simulation, the velocity, pressure, and temperature boundary conditions are set to the inlet, while no boundary condition is set for the outlet.
3. Hypersonic Flow Analysis
Take Ma 12 as the example, the flow field distributions are shown in Figure 8. A strong bow shock is formed at the leading of the structure due to strong compression, after which the flow deflects and travels along the surface with the Mach number decreasing rapidly, and a large amount of low-speed fluid accumulates near the stagnation point. The low-speed region also generates at the trailing edge due to the bluff body, forming a symmetric vortex. After crossing the shock wave, the temperature rises from 250 K to 1976 K and the pressure rises from 24 Pa to 4803 Pa near the stagnation point. In addition to the severe aeroheating on the windward, heating also occurs at the trailing edge due to the separation vortex. Besides, the high turbulent viscosity is mainly distributed near the shock wave, with the maximum value of 1.55 kg/(m·s).
(a)
(b)
(c)
(d)
In many research, the solved hypersonic flow around IRVE is circularly symmetric and no obvious pulsation is found, such as the study on the flow of Ma 12 using the -based SST turbulent model in Ref. [28]. This is inconsistent with the actual test phenomenon: the strong shock wave, viscous interference, and high-temperature effect existing in hypersonic flow will certainly lead to strong unsteady characteristics [29, 30]. Therefore, the LES WALE is adopted in this paper to make some attempts. As Figure 9 shows, a special physical image of Ma 12 can be obtained: the surface pressure is irregular, and the large-scale separated vortex on windward is discovered. However, the large vortex on windward is not obvious in all cases, only if the Mach number is greater than 6. For example, when the Mach number decreases to 5, the flow in front of windward is symmetrical without a separated vortex. To a certain extent, the difference between LES and RANS is little when the Mach number is below 6, especially for the simulation of aerodynamic loads on windward.
(a)
(b)
(c)
The generation mechanism of the vortex is analyzed by density gradient nephogram (see Figure 10). As we know, with the rise of the flight Mach number, the density of gas after shock wave increases, and therefore, the volume of fluid goes down. Therefore, the distance between shock wave and surface is small in hypersonic flow. On the other hand, the temperature of gas in the boundary layer will increase due to the viscous dissipation, leading to the drop of density; thus, the thickness of the boundary layer is large in hypersonic flow. Under these two factors, the shock wave and boundary layer have a tendency to interact with each other. However, when the Mach number decreases, the distance between shock wave and boundary layer becomes larger, causing little interaction and weakening vortex near the head cone, as Figure 10(b) shows.
(a)
(b)
The formation and development of the vortex are analyzed in detail in Figure 11. Originally, the interaction of shock wave and boundary layer causes the separation of the boundary layer, and a pair of large-scale vortices are generated from the bottom of the boundary layer, as observed in Figure 11(a). Affected by the vortex, the local pressure drops to a low level, and only at the flow reattachment region does the pressure recover. However, the vortex is unstable, and it will affect the inviscid fluid nearby. Take the vortex as a whole, the forces on the vortex are analyzed in Figure 11(c). There are three types of forces on the vortex: supporting force exerted by the IRVE system, driving force exerted by surrounding fluid (mainly driven by incoming flow), and viscous resistance exerted by the boundary layer. Owing to the angle between the incoming flow and surface, only the momentum in the direction parallel to the surface can become the effective driving force. If the driving force is greater than the resistance of viscous dissipation, the vortex will move downstream with the flow. Moreover, in the process of movement, the local pressure will reduce, and therefore, the flow faces with a severe adverse pressure gradient, which exacerbates the separation of the boundary layer. In the process of moving, the vortices on both sides of the head cone alternately generate, develop, and break. Affected by the vortex, the surface pressure shows strong pulsation, and the peak value of pressure is no longer located at the head cone (see Figure 11(b)).
(a)
(b)
(c)
For the purpose of quantitatively investigating the intensity and frequency of the vortex at different positions, six points (points to ) near skin depressions are established to monitor the static pressure, and the results are given in Figure 12. With the point moving from the head cone to the shoulder, the time-averaged pressure and maximum amplitude increase first and then decrease, reaching the peak value at point with these two values of 4210 Pa and 501 Pa, respectively. Different from the high frequency (105-107 Hz) and low amplitude of pressure pulsation in the turbulent boundary layer [31], the energy of the large-scale vortex on windward is large, but the frequencies are low (mainly under 400 Hz).
(a)
(b)
(c)
4. Aeroheating Analysis
The traditional aeroheating method about IRVE does not consider the effects of structural deformation and windward vortex captured by LES. In their research, the aeroheating mainly results from the shock wave, stagnation, and surface friction. Figure 13 reveals the difference when considering the aeroelasticity: the surface temperature is no longer circularly symmetric, and there is a growing trend on the whole skin. It is worth noting that the temperature in some regions is particularly high, while in other places, the temperature shows an obvious fluctuating distribution. This phenomenon is related to the large separation vortex on windward and the small vortex inside the depression, which will be specifically analyzed in the following.
(a)
(b)
The flow details near the skin are compared in Figure 14 to illustrate the effects of these two kinds of vortices. As for the large separation vortex, it is generated due to the interaction of shock wave and boundary layer near the head cone, which has been specifically analyzed above. The vortex will consume the kinetic energy and thereby produce huge heat dissipation, leading to a very high-temperature amplitude in the region where the large vortex propagates. However, the large vortex cannot cover all the skin, while in other regions without a large vortex, the effect of the small vortex inside depression is gradually reflected. Focusing on the flow inside the depressions, the fluid is blocked partly due to the impact on the front wall of depression, resulting in slight stagnation. Furthermore, with the growth in the amount of fluid stagnated on the front wall, the accumulated fluid will flow back in the depression, forming a vortex and circulating into the mainstream finally. This phenomenon will increase the thickness of the boundary layer and also exacerbate the heat transfer near depressions. However, the scale and dissipation effect of these two vortices are quite different; thus, the temperature curves also show significant differences (see Figure 13(b)).
Note that the large-scale vortex is accompanied by strong unsteady characteristics, so the surface heating will certainly be constantly changed. Choose three points from Figure 14 (named ) to display the temperature variations, in which point is on the head cone and points and are, respectively, within the influence of the large-scale vortex and disturbance by skin depression. As can be seen from Figure 15, the pulsation amplitudes of points (103 K) and (97 K) are remarkable, and the time-averaged temperature of point (1721 K) is 7.6% higher than the undeformed situation. These effects cause the maximum value of stagnation point (2031 K) to approach the allowable temperature (2100 K). Meanwhile, the peak region also changes due to the temperature pulsation, so it is not correct to think that the maximum heating always occurs at the head cone. It brings difficulties for thermal protection and indicates that the structural deformation and windward large-scale vortex solved by LES cannot be ignored in aeroheating analysis.
(a)
(b)
5. Nonlinear Structural Dynamics Analysis
Figure 16(a) presents the sectional views at a certain time under Ma 12. The deformation at the outer end is larger than that at the inner end overall, and the maximum value is located at the outermost skin with a value of 65 mm. In order to reflect the characteristics of deformation effectively, the meridians before and after deformation are compared in Figure 16(b). There are two kinds of deformation forms under aerodynamic loads: toroid bending and skin depression. Toroid bending refers to the backward deformation of the whole structure centered on the rigid cone, while skin depression is a kind of local deformation due to the gaps between toroids and skin.
(a)
(b)
The vibration of toroid reflects the global response of the IRVE system. From Figure 17(a), the vibration with large amplitude occurs around the balance position firstly, and then, the amplitude is maintained within a certain range. In the whole process, the time-averaged deformation is 48 mm and the maximum amplitude reaches 17 mm. According to the frequency spectrum (see Figure 17(b)), the dominant frequency is 16 Hz and the secondary frequency is 54 Hz, which correspond to the first two natural frequencies of the IRVE system, and the other frequencies are consistent with aerodynamic pulsation. Therefore, it performs a free attenuated vibration dominated by inertial force in the early stage but then becomes a forced vibration induced by unsteady aerodynamic force when the amplitude is stable. Combined with the results of other flight conditions (Ma 5-12), the structural vibrations are all maintained at a certain amplitude without divergence; thus, the global instability does not occur during the hypersonic region.
(a)
(b)
The depths of six depressions are monitored in Figure 18 to analyze the local vibration of skin. The closer to the shoulder, the larger the depression depth is. The time-averaged deformation of the outermost depression is 23.6 mm, and the maximum amplitude is 4.4 mm. Different from the global vibration of toroids, the free attenuated process is not significant; it mainly performs a forced response induced by aerodynamic force.
(a)
(b)
The difference in the free attenuated process after the initial incentive results from the stiffness and damping. The stiffness of toroid is provided by inflation gas; its stiffness and structural damping are commonly low, causing slow amplitude attenuation during free vibration. But for skin, its deformation is severely constrained when adhered to the toroid, resulting in higher structural stiffness. Therefore, the large vibration of toroid means that it is meaningful to study the natural vibration properties, especially when bearing the huge aeroheating. As we know, the surface temperature and inflation gas expansion caused by aeroheating have opposite effect on stiffness: high temperature will cause material softening and thereby reduce the elastic modulus and natural frequency, while inflation gas expansion will increase the stiffness. In order to distinguish their influence degrees, a comparison based on the control-variable method is conducted. As can be seen from Table 2, the first- and second-order frequencies are mainly affected by inflation gas expansion, whose rise, respectively, reaches 63% and 35%. Besides, the modal shapes do not change in thermal environment.
The huge difference of vibrations between toroid and skin indicates that the stress analysis also needs to separate them. Table 3 illustrates the independent influence of each factor on stress, and the results are compared with the state that only initial inflation pressure (20 kPa) is exerted. From the results, the aerodynamic force and surface temperature mainly affect skin stress, while the thermal expansion has a great influence on toroid stress.
6. Conclusions
The thermal-fluid-solid coupling model proposed in this paper can present the actual physical process more effectively than the traditional method of separated aeroheating and aeroelasticity. For the aeroheating problem, this model is competent in describing the heat conduction and inflation gas expansion, overcoming the defect that structural deformation is not considered in existing research. For the aeroelastic problem, this model considers the effect of aeroheating more than the traditional aeroelastic method, and it can clearly reveal the mechanism of surface temperature and inflation gas on nonlinear structural dynamics. Based on the model, some important phenomena are discovered: (1)In addition to the heating caused by shock wave, stagnation, and surface friction, the aeroheating in hypersonic flow is also affected by the large-scale separated vortex on windward, which is generated due to the interaction of shock wave and boundary layer near the head cone. The vortex shows strong pulsation characteristic and produces huge heat dissipation, resulting in a temperature pulsation amplitude of 103 K on the stagnation point(2)The skin depressions caused by structural deformation will hinder the flow near the surface, causing a thicker boundary layer and larger surface friction. This effect obviously intensifies the local heating and brings a fluctuating distribution of temperature, leading to a temperature rise of 7.6% near depression(3)From Ma 5 to Ma 12, the vibration amplitudes of toroids and skin are all maintained at a certain level without divergence; thus, both the global instability of toroid and local instability of skin do not occur under hypersonic conditions(4)Under aerodynamic force, the flexible structure performs violent vibration, which is essentially a forced response induced by a large-scale separated vortex on windward. The generation and breakdown of the vortex play the role of unsteady excitation, causing the maximum dynamic deformation of toroids and skin, respectively, to reach 65 mm and 28 mm(5)Under aeroheating effects, the high surface temperature and inflation gas expansion have great influence on thermal stress and natural vibration properties. On the one hand, thermal stress is mainly reflected on skin and has little effect on toroid. On the other hand, the inflation pressure determines the structural stiffness, and the first-order frequency increases by 63% after thermal expansion
Data Availability
The numerical data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
The research was supported by the National Natural Science Foundation of China (11602018) and the Special Funds of Foshan Technology Innovation Team of China (2016lt100123).