Study on Transient Dynamics and Aerodynamic Characteristics of a New Type of High-Damping Four-Winged Rotating Parachute Inflation Process

Abstract

To decrease the opening shock as well as improve the resistance coefficient and stability of the parachute of aviation weapon under subsonic conditions, a new kind of high-damping four-winged rotating (HFWR) parachute is investigated in this paper. The transient dynamic behavior and aerodynamic characteristics of the parachute during the inflation process are studied. Considering the permeability, the 3D folded finite element (FE) model of the HFWR parachute is established based on the direct folding modeling technique, and the inflation process of the parachute under subsonic flow is simulated using the multimaterial arbitrary Lagrange–Euler (ALE) method. A series of wind tunnel tests are conducted to verify the numerical results. Besides, the opening performances of the HFWR parachute and the round parachute, which includes the inflation process, the dynamic response of the swing angle, and the opening shock load varying with time, are compared under the same conditions. The results show that the opening performance of the HFWR parachute is superior to the round parachute under specific military background. The fluid-structure interaction (FSI) simulation results show good consistency with the wind tunnel tests, which indicates that the numerical modeling can effectively simulate and predict the opening performance and aerodynamic characteristics of the rotating parachute. The modeling method in this paper can help shorten the development cycle, improve the cost effectiveness, and optimize the design of the parachute.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Parachutes are widely used as devices for stability, posture adjustment, and deceleration in the air for stand-off aviation weapon systems. With the improvement of the tactical and technical performance and the operational requirements, it is necessary to design various new types of parachutes, which have a great impact on trajectory and attitude, to improve the stability and reliability of the parachute-payload system. To meet the tactical requirements of good stability, short endurance, and small opening shock load and volume and weight as well as low cost, the rotating parachute is employed in the aviation weapon system. Effective inflation of the parachute is not only the premise of maintaining the reliability of the weapon system but also the key to complete the combat mission. Therefore, it is necessary to investigate the parachute inflation process in depth. However, the inflation process of flexible parachute involves the transient large deformation and fluid-structure interaction (FSI) problems of structural dynamics with the aerodynamics of the surrounding flow field in which geometric nonlinearity and material nonlinearity coexist. Furthermore, since the inflation process of parachutes is generally within 1 second, it is very difficult to investigate the parachute inflation process by the experimental methods. Therefore, the numerical simulation method becomes an important alternative mean. At present, the research on the numerical simulation of the parachute inflation process is still in its infancy, which is mainly focusing on the round parachute, and there is still short of relevant research studies on the rotating parachute inflation process. Compared with the round drag parachute, the transient dynamic characteristics of the rotating parachute in the nonlinear, large deformation, and FSI inflation process will be more complicated due to the rotating characteristics. Therefore, it is necessary to study and explore the FSI process of parachute inflation utilizing numerical simulation.

In recent years, a great many advanced numerical simulation techniques and methods have been developed to solve parachute FSI coupling problems, such as the mass-spring-damper (MSD) model, deformable spatial domain/stabilized space-time (DSD/SST) technique, and arbitrary Lagrange–Euler (ALE) coupling method. The MSD model was widely used in the computational simulation of parachute dynamic coupling in the 1980s. Purvis coupled the kinematic equation of the canopy with that of the fluid in a tight coupling manner, but both the canopy structure model and the flow field model were simplified [1]. Then, Purvis improved the coupling model, in which the canopy was used as the MSD model and the two-dimensional Euler equation was used to calculate the flow field [2]. Later on, Stein and Benney used a simplified arbitrary Lagrangian–Eulerian (SALE) CFD program and the MSD model of the canopy to calculate the inflation process of the C-9 flat circular parachute in loose coupling manner [3, 4]. Using the MSD model and standard two-equation turbulence model, Yu and Ming simulated the main inflation process of a flat circle parachute in a loose coupling manner [5]. In addition to the MSD model, Tezduyar’s team simulated the fluid-structure interaction (FSI) modeling of parachutes based on the deformable spatial domain/stabilized space-time (DSD/SST) technique. They combined this format with the membrane FE model of the canopy to calculate the FSI process of the T-10 flat circular parachute in a loose coupling manner [6, 7]. To make the algorithm more stable and efficient for lighter structures, the coupling method was improved from loose coupling to quasidirect and direct coupling. Then, they proposed a new generation stabilized space-time fluid-structure interaction (SSTFSI) technique based on the DSD/SST technique, which has become the core technique of parachute numerical simulation of the team [8, 9]. Moreover, the team used a homogeneous geometric porosity model (HMGP) to solve the problems caused by structural permeability and simulated the steady descending state of a single ring-sail parachute and three ring-sail parachutes [10–13]. To solve the contact problem of the structural surface, they proposed the surface-edge-node contact tracking (SENCT) technique, which has been further developed into the basic technology of FSI simulation for parachute clusters. They used those techniques to perform the FSI simulation that focuses on parachute disreefing, including the disreefing of parachute clusters and the reefed stages of ring-sail parachutes [14–18]. Besides, a new explicit method based on the ALE coupling method in LS-DYNA code has been used to simulate and analyze the FSI problems of parachute inflation. Tutt and Taylor used the Euler–Lagrange penalty algorithm to simulate the inflation process of a flat circular parachute [19]. Then, they applied a new Euler–Lagrange coupling algorithm to introduce the permeability into the FSI simulation and simulated the inflation performance of the disk-gap-band parachute for supersonic Mars exploration under the influence of the wake of capsule [20–23]. Subsequently, they simulated the inflation process of a flat circular parachute in a finite mass scenario and conducted a series of indoor vertical parachute tests to verify the simulation results [24]. Cheng et al. simulated the inflation process of a C-9 flat circle parachute with this method and then transformed the shape of the canopy into a CFD model that uses a porous medium model to simulate the permeability [25]. Gao et al. used a multimaterial arbitrary Lagrange–Euler-coupled numerical method to simulate the inflation process of a slot parachute [26]. Then, they analyzed the FSI inflation phenomenon of a lifesaving slot parachute in low-speed finite mass inflation process by using the arbitrary Lagrange–Euler coupling penalty method, and they studied the coupling behaviors of a lifesaving parachute’s inflation process based on the multibody dynamic model and LS-DYNA nonlinear analysis code in low attitude [27, 28]. Moreover, they studied the transient dynamic behavior of a supersonic disk-gap-band parachute in a Mars entry environment involving FSI, and the numerical results were verified by the supersonic wind tunnel test data from NASA [29].

This paper focuses on the analysis of the dynamic coupling behaviors of a new kind of high-damping four-winged rotating (HFWR) parachute during inflation in an infinite mass scenario. The rest of the paper is organized as follows. Firstly, the structure layout and the folded model of the HFWR parachute are introduced in Section 2. Secondly, the governing equations of structure dynamics and fluid dynamics and the FSI coupling method are given in Section 3. Thirdly, the FE models of the HFWR parachute and fluid domain are established by utilizing the direct folding modeling technique and finite element method in Section 4. To verify the simulation results, the wind tunnel tests are conducted, and the test method and equipment of the wind tunnel test are introduced in Section 5. In Section 6, the results of transient dynamic behavior, inflation performance, and aerodynamic characteristics of the HFWR parachute are analyzed and compared with a round parachute. Finally, Section 7 concludes the paper.

2. Problem Statement and Modeling

2.1. Problem Statement

During the falling of the parachute-payload system, the HFWR parachute undergoes an inflation and rotation process and finally reaches a steady state, thereby playing a role of stabilization and deceleration. The HFWR parachute studied is suitable for stand-off aviation weapons under subsonic conditions, and all the computations are carried out using air properties at standard sea-level conditions. The air density is 1.226 kg/m³, and kinematic viscosity is 1.78 × 10⁻⁵ Pa·s.

2.2. Layout of the HFWR Parachute

The HFWR parachute, with a circular vent on the top, consists of canopy made of 4 pieces of nylon fabric gores and 28 suspension lines as well as some reinforcement belts. A stiffener connects the skirts of the canopy. The shape of parachute inflation is affected by the layout and length of the suspension lines, and thus the suspension lines are divided into 4 groups, which symmetrically distribute on the skirts of the canopy and intersect at the bottom. The suspension lines are different in length to ensure that the canopy forms air vent with the same direction of rotation. The top view of the parachute with expanding canopy is shown in Figure 1(a), and the layout of the suspension lines and canopy is depicted in Figure 1(b).

[figures omitted; refer to PDF]

2.3. Modeling of Folded Canopy

The parachute is folded and encapsulated in pack before inflating in practical applications. Modeling of the initial state of the parachute has a significant influence on the simulation results. The direct folding modeling technique is applied to establish the initial model of full-size parachute due to a compromise between the complexity of modeling flexible folded fabric canopy and the accuracy of calculation, as shown in Figure 2. The direct folding modeling technique is simple and effective, and more structure details, such as reinforcement belt and geometric porosity, can be described clearly. Furthermore, it is convenient to adjust the folded extent and the initial size of the air inlet as needed.

[figure omitted; refer to PDF]

The main structural parameters of the folded model of the parachute are calculated as follows. Each gore of the canopy is initially folded in the radial direction with an angle $\beta \approx 10$ degree, a top edge $L_{\text{arc}}=\left(\pi D_a/2N\right)$, edge of folded line $H=\left(D_c-D_a\right)/2$, and side edge $R=\sqrt{{\left(L-D_a/2\right)}^2+H^2}$, by which a profile of a single wing can be constructed. Here, D_c, D_a, and L are the structural diameter, aperture diameter, and half the length of the bottom edge of the canopy, respectively. Since the parachute is geometrically symmetrical, the folded model of the canopy shape with 4 gores can be obtained by placing 4 single wings evenly around the axis. The initial air inlet area of the folded model is 0.0305 m².

3. Mathematical Model

3.1. Special Considerations

Considering the complex dynamic characteristics of the HFWR parachute, the following hypotheses are made before the governing equations are given:

(1)

The influence of the interaction between canopy fabrics is neglected.

(2)

The effect of the initial folded model of the parachute on the opening performance is not considered.

(3)

The payload of the parachute-payload system is simplified as a massless point fixed at the inlet boundary. The aerodynamic forces of the payload and the suspension lines are ignored, and the tensions of the suspension lines are considered.

(4)

Regardless of the influence of crosswind, the attack angle of the inlet flow is zero.

(5)

The inlet flow is considered a steady incompressible ideal fluid, whose behavior is time-dependent. The velocity of flow at the inlet boundary is constant, and the nonreflective boundary is used for the remaining boundaries of the fluid domain to simulate infinite space.

(6)

The porosity of the canopy is assumed constant. The air density and dynamic viscosity are supposedly uniform.

3.2. Structure Dynamics

3.2.1. Membrane Shell of the Porous Canopy

Considering the canopy structure is made of flexible fabric, a two-dimensional membrane shell with large flexible and porous permeability can characterize the fabric of the canopy. Let ${\Omega }^s\in R^3$ denote the domain of the thin porous medium structure of canopy, which is discretized by Lagrangian shells based on the Belytschko–Lin–Tsay formulation [30], and $\partial {\Omega }^s$ denote its boundary, including the displacement boundary and the traction boundary. The movement of the thin porous medium structure of canopy ${\Omega }^s$ is described by $x_i\left(t\right),\left(i=\mathrm{1,2,3}\right)$, which can be expressed in terms of the reference coordinates $X_i\left(t\right),\left(i=\mathrm{1,2,3}\right)$ and time t, that is,$$\begin{array}{}\text{(1)}& x_i=x_i\left(X_{\alpha },t\right).\end{array}$$

The governing equations for the structure of porous media membrane are given as follows:$$\begin{array}{}\text{(2)}& {\rho }_s\frac{\text{d}\stackrel{⟶}{v}}{\text{d}t}=\stackrel{⟶}{\text{div}}\left(\overline{\overline{{\sigma }_s}}\right)+\stackrel{⟶}{f_s},\text{(3)}& {\rho }_s\frac{\text{d}e}{\text{d}t}=\overline{\overline{{\sigma }_s}}:\overline{\overline{\text{grad}}}\stackrel{⟶}{\left(v\right)}+{\stackrel{⟶}{f}}_s\cdot \stackrel{⟶}{v}.\end{array}$$

The boundary conditions are a specified displacement and a traction, where the exterior unit normal is denoted as ${\stackrel{⟶}{n}}_s$. The boundaries where the displacement and traction boundary conditions are applied are as follows:$$\begin{array}{}\text{(4)}& \stackrel{⟶}{x}\left(\stackrel{⟶}{X},t\right)=\stackrel{⟶}{D}\left(t\right)\mathrm{on}\partial {\Omega }_1^s,\text{(5)}& \overline{\overline{{\sigma }_s}}\cdot {\stackrel{⟶}{n}}_s=\stackrel{⟶}{T}\left(t\right)\mathrm{on}\partial {\Omega }_2^s,\end{array}$$where $\overline{\overline{{\sigma }_s}}$ is the Cauchy stress of the structure, ${\rho }_s$ is the density of the structure, $f_s$ is the force on the structure, e is the specific internal energy, $\text{d}\stackrel{⟶}{v}/\text{d}t$ is the acceleration, and $\stackrel{⟶}{D}$ and $\stackrel{⟶}{T}$ are the displacement vector and stress vector of the structure.

The solution of equations (2) and (3) is subject to equations (4) and (5) corresponding to the displacement boundary condition on boundary $\partial {\Omega }_1^s$ and the traction boundary condition on the boundary $\partial {\Omega }_2^s$, respectively. The traction boundary is the fluid-structure interaction surface, where the air pressure acts on the fabric canopy.

The Belytschko–Lin–Tsay membrane with a fabric material model (^∗MAT_FABRIC) is suited for a 3- or 4-node element, which is of nonlinear dynamic characteristics and large deformation. It can well simulate the thin porous medium structure of the canopy and the stress-strain relationship of the membrane, which is given as follows:$$\begin{array}{c}\text{(6)}& {\epsilon }_1=\frac{1}{E_1}\left({\sigma }_1-{\nu }_1{\sigma }_2\right),{\epsilon }_2=\frac{1}{E_2}\left({\sigma }_2-{\nu }_2{\sigma }_1\right),\\ 2{\epsilon }_{12}=\frac{1}{G_{12}}{\tau }_{12}+\alpha {\tau }_{12}^3,\end{array}$$where ${\sigma }_i$, ${\nu }_i$, and $E_i$$\left(i=\mathrm{1,2}\right)$ are, respectively, stress, Poisson’s ratio, and elastic modulus, in which $i=1$ means longitudinal and $i=2$ represents transverse; ${\tau }_{12}$ and ${\epsilon }_{12}$ represent shear stress and shear strain; $G_{12}$ is shear elasticity; and α denotes nonlinear coefficient which can be measured by stress-strain relationship.

3.2.2. Elastic Cable

The suspension lines are modeled by nonlinear discrete beam elements with an elastic cable material model (^∗MAT_CABLE_DISCRETE_BEAM); thus, no force will be produced in compression. The force, $F_t$, generated by the cable is nonzero if and only if the cable is under tension. The force is given by$$\begin{array}{}\text{(7)}& F_t=\max \left(F_0+K\cdot \Delta L,0\right),\end{array}$$where $F_0$ is the initial force and $\Delta L$ is the change in length, described as$$\begin{array}{}\text{(8)}& \Delta L=l_{\text{current}}-\left(l_0-l_{\text{offset}}\right),\end{array}$$and the stiffness K ($E\ge 0$ only) is defined as$$\begin{array}{}\text{(9)}& K=\frac{E\cdot A_c}{\left(l_0-l_{\text{offset}}\right)},\end{array}$$where $l_{\text{current}}$ is the current length of the element, $l_0$ is the initial length of the element, $l_{\text{offset}}$ denotes offset of the element length, E represents Young’s modulus, and $A_c$ denotes the cross-sectional area of the cable. For a slack cable, the offset should be input as a negative length. For an initial tensile force, the offset should be positive.

3.3. Fluid Dynamics

The Naiver–Stokes governing equations of incompressible flows are solved by the ALE method [31] (^∗SECTION_SOLID_ALE). Equations (10)–(12) give the equations of mass, momentum, and energy conservation for a Newtonian fluid in ALE formulation in the reference domain. Let ${\Omega }^f\in R^3$ represent the spatial domain occupied by air defined by ^∗MAT_NULL and $\partial {\Omega }^f$ denote its boundary.$$\begin{array}{}\text{(10)}& \frac{\partial {\rho }_f}{\partial t}+{\rho }_f\text{div}\left(\stackrel{⟶}{v}\right)+\left(\stackrel{⟶}{v}-\stackrel{⟶}{w}\right)\text{grad}\left({\rho }_f\right)=0,\text{(11)}& {\rho }_f\frac{\partial \stackrel{⟶}{v}}{\partial t}+{\rho }_f\left(\stackrel{⟶}{v}-\stackrel{⟶}{w}\right)\overline{\overline{\text{grad}}}\left(\stackrel{⟶}{v}\right)=\stackrel{⟶}{\text{div}}\left(\overline{\overline{{\sigma }_f}}\right)+\stackrel{⟶}{f},\text{(12)}& {\rho }_f\frac{\partial e}{\partial t}+{\rho }_f\left(\stackrel{⟶}{v}-\stackrel{⟶}{w}\right)\overline{\overline{\text{grad}}}\left(e\right)=\overline{\overline{{\sigma }_f}}:\overline{\overline{\text{grad}}}\left(\stackrel{⟶}{v}\right)+\stackrel{⟶}{f}\cdot \stackrel{⟶}{v},\end{array}$$where ${\rho }_f$ is the density of the fluid, $\stackrel{⟶}{f}$ is the body force, e is the specific internal energy, and $\stackrel{⟶}{v}$ and $\stackrel{⟶}{w}$ are the fluid and mesh velocity fields, respectively. In the Eulerian formulation $\stackrel{⟶}{w}=0$, this assumption eliminates the remeshing and smoothing process, but does not simplify the Navier–Stokes equations equations (10)–(12) which can be solved by the split approach [32, 33] and implemented in LS-DYNA®. $\overline{\overline{{\sigma }_f}}$ is the total Cauchy stress of the fluid, given by$$\begin{array}{}\text{(13)}& \overline{\overline{{\sigma }_f}}=-p\cdot \overline{\overline{Id}}+\mu \left(\overline{\overline{\text{grad}}}\left(\stackrel{⟶}{v}\right)+\overline{\overline{\text{grad}}}{\left(\stackrel{⟶}{v}\right)}^T\right),\end{array}$$where p is the pressure, $Id$ is a unit second order tensor, and μ is the dynamic viscosity of the fluid. The pressure is calculated by the ideal gas law (^∗EOS_IDEAL_GAS):$$\begin{array}{}\text{(14)}& p={\rho }_f\left(Cp-Cv\right)T_f,\end{array}$$where ${\rho }_f$ and $T_f$ are the air density and temperature and $Cp$ and $Cv$ are the specific heat capacities at constant pressure and volume, respectively.

Equations (10)–(12) are completed with appropriate boundary conditions. In this part, both the inflow boundary condition where the velocity is specified and the traction boundary condition which is the fluid-structure interaction surface $\left(\partial {\Omega }_2^f=\partial {\Omega }_2^s\right)$ are expressed as follows:$$\begin{array}{c}\text{(15)}& \stackrel{⟶}{v}=\stackrel{⟶}{g}\left(t\right)\mathrm{on}\partial {\Omega }_1^f,\overline{\overline{{\sigma }_f}}\cdot {\stackrel{⟶}{n}}_f=\stackrel{⟶}{h}\left(t\right)\mathrm{on}\partial {\Omega }_2^f,\end{array}$$where $\partial {\Omega }_1^f$ is the velocity boundary of the fluid field, $\partial {\Omega }_2^f$ is the traction boundary that is actually the parachute, $\stackrel{⟶}{g}$ and $\stackrel{⟶}{h}$ are the displacement vector and stress vector of the fluid, and ${\stackrel{⟶}{n}}_f$ is the exterior unit normal of the fluid boundary. The traction force is computed by the porous coupling method in the ALE formulation (^∗CONSTRAINED_LAGRANGE_IN_SOLID). This coupling method will be described in the next section.

3.4. FSI Coupling Method

The porous Euler–Lagrange coupling method is employed to solve the fast transient porous FSI problem because Eulerian formulation can simulate the large deformation of fluid and coupling can deal with the interaction between fluid and thin porous medium. The calculation of nodal forces is the main part of this process. The fluid and structural nodal forces and the coupling forces only affecting the nodes on the fluid-porous structure interface are calculated. The principle of this coupling is presented as follows. The slave node is a structure mesh node, while the master node is not a fluid mesh node. It can be regarded as a fluid particle within a fluid element, and the mass and velocity can be interpolated from the fluid element nodes by using finite element shape functions. Using the isoparametric coordinates of the fluid element, the location of the master node is calculated.

For each structure node, the penetration depth $\stackrel{⟶}{d}$ of each step is updated incrementally by the relative velocity ${\stackrel{⟶}{v}}_{\text{rel}}$ of the slave node and the master node in equation (16), that is,$$\begin{array}{}\text{(16)}& {\stackrel{⟶}{d}}^{n+1}={\stackrel{⟶}{d}}^n+{\stackrel{⟶}{v}}_{\text{rel}}^{\left(n+1\right)/2}\cdot \Delta t,\end{array}$$where ${\stackrel{⟶}{d}}^n$ is the penetration depth at time $t=t^n$ and the relative velocity ${\stackrel{⟶}{v}}_{\text{rel}}$ can be represented as follows:$$\begin{array}{}\text{(17)}& {\stackrel{⟶}{v}}_{\text{rel}}^{\left(n+1\right)/2}={\stackrel{⟶}{v}}_s^{\left(n+1\right)/2}-{\stackrel{⟶}{v}}_f^{\left(n+1\right)/2},\end{array}$$where the fluid velocity ${\stackrel{⟶}{v}}_f$ is the velocity at the master node and the structure velocity ${\stackrel{⟶}{v}}_s$ is the velocity at the slave node. The coupling acts only if penetration occurs, namely, ${\stackrel{⟶}{n}}_s\cdot {\stackrel{⟶}{d}}^n<0$, where ${\stackrel{⟶}{n}}_s$ is obtained by averaging the normal of structure elements connected to the structure node.

The porous coupling forces are deduced by integrating the Ergun equation [34] of shell volume, that is,$$\begin{array}{}\text{(18)}& \frac{\text{d}p}{\text{d}r}=a\left(\mu ,\epsilon \right)\cdot {\stackrel{⟶}{v}}_{\text{rel}}\cdot {\stackrel{⟶}{n}}_s+b\left(\rho ,\epsilon \right)\cdot {\left({\stackrel{⟶}{v}}_{\text{rel}}\cdot {\stackrel{⟶}{n}}_s\right)}^2,\end{array}$$where r denotes normal direction of the shell element, ε represents porosity, $a\left(\mu ,\epsilon \right)$ is reciprocal permeability of the porous shell or viscous coefficient, and $b\left(\rho ,\epsilon \right)$ is inertia coefficient. Using the Ergun theory, the equations for calculating the coefficients of a and b can be obtained, namely,$$\begin{array}{c}\text{(19)}& a\left(\mu ,\epsilon \right)=\frac{150\mu {\left(1-\epsilon \right)}^2}{D^2{\epsilon }^3},b\left(\rho ,\epsilon \right)=\frac{1.75\rho \left(1-\epsilon \right)}{D{\epsilon }^3},\end{array}$$where D is the characteristic length defined by $D=6\left(1-\epsilon \right)V/S$, in which $V$ denotes the volume of the canopy and S is the permeated surface of the canopy.

Assume that air density and dynamic viscosity are constant, and then the viscous and inertia parameters in equation (18) are constant. By fitting the Ergun theoretical permeability with the experimental data from the research of Wang et al. [20], the values of the viscous and inertia parameters of MIL-C-7020 type III fabric, which is made up of the nylon canopy, are a = 1599174 kg/ms² and b = 48051 kg/m³.

Then, the porous coupling force F derived from equation (18) is applied to both the master and the slave nodes in opposite directions to satisfy the principle of action-reaction of forces at the interface coupling. Applying the porous coupling force to a structure coupling node gives$$\begin{array}{}\text{(20)}& F_s=-F.\end{array}$$

For the fluid, the porous coupling force is distributed to the fluid element nodes based on the shape functions at each node $i\left(i=1,\dots ,8\right)$, and the fluid force is scaled by the shape function $N_i$, described as$$\begin{array}{}\text{(21)}& F_f^i=N_i\cdot F,\end{array}$$where $N_i$ is the shape function at node i and the fluid force for a hexahedron element can be calculated as follows:$$\begin{array}{}\text{(22)}& {\displaystyle \sum _i^8{F}_f^i}=F.\end{array}$$

4. Numerical Simulation

4.1. Finite Element Model

The finite element (FE) model of the HFWR parachute and the fluid domain are established, respectively. Based on the folded model defined in Section 2, the FE model of the parachute is formed as shown in Figure 3. The canopy is meshed as a two-dimension (2D) tetrahedral membrane shell, and the suspension lines and reinforcement belts are meshed with one-dimensional discrete beam elements.

[figures omitted; refer to PDF]

To make a compromise between the calculation accuracy and efficiency, a gradual mesh of cubical fluid domain with a size of 5D_c $\times$ 5D_c $\times$ 8D_c is built up as shown in Figure 4. The fluid domain is meshed by three-dimensional (3D) hexahedron elements. The fluid grid size near the FE model of the parachute is gradually diminished, and the size of the fluid grid becomes larger gradually from the inside to the outside stage by stage. The fluid domain entrance adopts the velocity inlet boundary condition. To reduce the effect of wave reflection, a nonreflective boundary condition is employed on other walls. The parachute is placed at the center of the fluid domain near the inlet flow boundary. The statistical information for the FE model is presented in Table 1.

[figures omitted; refer to PDF]

Table 1

Statistical information of FE model.

Components	Elements		Materials	Parts
	Types	Numbers
Canopy	Tetrahedral shell	1,628	Fabric	Part 1
Reinforcement belt	Discrete beam	844	Cable	Part 2
Suspension line	Discrete beam	609	Cable	Part 3
Fluid	Hexahedral solid	93,600	Ideal gas	Part 4-5
Plant	Hexahedral solid	16	Rigid	Part 6
Total	—	100,883	5	6

4.2. Material Models

The canopy of the parachute is fabricated from nylon fabrics with low permeability, which is constructed using the FABRIC material model for airbag inflation in the LS-DYNA material library. Both edges of the canopy and the vent are reinforced with high-strength nylon fabric materials. The CABLE material model is used to simulate the mechanical behavior of the suspension lines and reinforcement belts. The material properties of the HFWR parachute are listed in Table 2.

Table 2

Material properties of the HFWR parachute.

Material property	Canopy	Reinforcement belt	Suspension line
	Membranes	Cable	Cable
Thickness/cross-sectional area	$1.0\times {10}^{-4}\text{m}$	$1.3\times {10}^{-8}{\text{m}}^2$	$7\text{.}07\times {10}^{10}{\text{m}}^2$
Density (kg/m3)	533.96	1194	1140
Young’s modulus (Pa)	$4\text{.}309\times {10}^8\text{m}$	$1.3\times {10}^{10}$	$1.63\times {10}^{10}$
Poisson’s ratio	0.14	—	—

5. Wind Tunnel Test

5.1. Test Model and Equipment

In order to verify the validity of the simulation results and evaluate the blowing resistance and strength of the HFWR parachute under subsonic airflow environment, low-speed wind tunnel tests and subsonic airflow blow tests were carried out. The low-speed wind tunnel tests were carried out in a single-reflux open low-speed wind tunnel. As shown in Figure 5(a), the test section of the wind tunnel is ϕ2.5 m $\times$ 5.6 m in size and the maximum test flow speed is 100 m/s. The low-speed wind tunnel can simultaneously complete the measurement of 15 force components, and the conventional pressure test can simultaneously complete the 64-channel test. The data acquisition module of the measurement and control system based on VXI bus technology has a high-speed, high-precision data sampling function, and the system accuracy is better than 0.1%. The main performance and flow field quality of low-speed wind tunnel are shown in Table 3. The subsonic tests were conducted in a transonic wind tunnel of an intermittent blow-down type. As shown in Figure 5(b), the cross section of the nozzle of the transonic wind tunnel is a 1.4 m $\times$ 1.0 m ellipse. It adopts two control modes: front chamber total pressure closed-loop control and pressure-regulating valve open-loop control, which can realize real-time continuous stepless regulation of transonic speed of 85 m/s–476 m/s. The accuracy of the velocity pressure control is less than 2%, and the fluctuation of the velocity pressure is less than 1.5%.

[figures omitted; refer to PDF]

Table 3

Main performance and flow field quality of low-speed wind tunnel.

Main performance		Flow field quality
Test section size	ϕ2.5 m $\times$ 5.6 m	Dynamic pressure field	$\le 0.5\text{\%}$
Maximum wind speed	100 m/s	Orientation field	$\overline{\Delta \alpha }\le {0.14}^{\circ }$
			$\overline{\Delta \beta }\le {0.21}^{\circ }$
Stagnation pressure	$1.0\times {10}^5\text{Pa}$	Turbulence	0.11%
Dynamic pressure	5600 Pa	Dynamic pressure stability	$\le 0.004$
Kinetic energy ratio	3.5	Axial static pressure gradient	0.0036
Motor power	1800 kW	Maximum airflow temperature	$\le {40}^{\circ }\text{C}$

Considering the errors caused by the scaled model, a full-size test model of the HFWR parachute is employed, which is equal to the simulation model, as shown in Figures 5(c) and 5(d). The test content includes observing and recording the dynamic opening process of the HFWR parachute and measuring the inflation time, opening shock, and the steady rotation rate of the parachute. The layout of the test model of the parachute and equipment is shown in Figure 6. The test equipment includes a support device, force balance, high-speed camera, data acquisition system, and computer monitoring system. The support device is a single support rod device. The force balance is a beam strain-gauge balance. The high-speed camera is used to shoot the dynamic opening process of the parachute. The data acquisition system is used for data recording and real-time processing, and the computer monitoring system is for controlling the test system.

[figure omitted; refer to PDF]

5.2. Test Methods and Conditions

The parachute model is fixed on the single support rod device, and the beam strain-gauge balance is placed inside the support rod, which connects the parachute model with the force balance. When the airflow in the test section is stable at the predetermined dynamic pressure, the parachute-opening controller is activated to open the pack, and then the test model is pulled out by the action of the parachute-opening device. At the same time, the data acquisition system, the computer monitoring system, the high-speed camera, and the force balance start working. The data are recorded until the parachute is filled and reaches the steady state. The sensitive time of the high-speed camera and force balance is within 0.01 s, and the whole test time is about 5 s. The parachute-opening characteristic parameters, such as maximum opening shock load, resistance characteristic, inflation time, and steady rotation rate, can be determined by combining the opening shock results and time scale with high-speed camera recording.

The tests at an airflow speed of 50 m/s were conducted in a low-speed wind tunnel. Besides, to reduce the wall effect of the wind tunnel and test the blow resistance and strength of the HFWR parachute under high-speed airflow, the wind tunnel tests were carried out at subsonic airflow speeds of 200 m/s, 245 m/s, and 280 m/s, respectively. Table 4 gives the flow conditions of the test.

Table 4

Flow conditions of the test.

	$v_{\infty }$ (m/s)	$q_{\infty }$ (Pa)
Low speed	50	1506
	200	24100
Subsonic	245	36165
	280	47236

6. Results and Discussion

6.1. Parachute Inflation Dynamics

6.1.1. Opening Process and Stability

To study the inflation characteristics of the parachute, the inflation process of the HFWR parachute and round parachute with the same structural diameter and aperture diameter is simulated and compared at a flow speed of 50 m/s. The canopy deformation of the two kinds of parachute is presented. As illustrated in Figure 7, the airflow firstly acts on the fringe of the canopy skirt of the HFWR parachute, and then the canopy gradually inflates and gathers slightly towards the apex as the airflow enters. Under the action of the tension of suspension lines and aerodynamics, the canopy gradually tilts and starts to rotate, then gradually fills into a shape similar to “helicopter screws,” and finally reaches the steady state after fully inflated. The time to reach the steady state is about 0.3 s after inflating.

[figures omitted; refer to PDF]

Compared with the round parachute, the HFWR parachute has no obvious initial and main inflation stage and no obvious “breathing phenomenon.” This is due to the large gaps between the gores of the HFWR parachute, resulting in a relatively larger porosity and permeability of the canopy structure than the round parachute.

Figure 8 shows the comparison of the fully inflated shape of the HFWR parachute for the numerical simulation and wind tunnel test. It shows that the opening shape of the numerical model is consistent with the real parachute.

[figures omitted; refer to PDF]

Figure 9 shows the comparison of the swing angles around the axis of the parachute between the HFWR parachute and the round parachute under the same flow velocity. The swing of the HFWR parachute is slightly larger than that of the round parachute at the beginning. On the one hand, the suspension lines of the HFWR parachute initial numerical model are not straightened, while suspension lines of the round parachute numerical model are of equal length and straightened; on the other hand, the large structural permeability of the HFWR parachute results in turbulence of airflow, thus causing the parachute to swing. However, with constant inflation, the swing amplitude of the HFWR parachute gradually reaches a stable level and fluctuates slightly around 5 degrees, while the swing angle of the round parachute becomes larger and larger, showing irregular large fluctuations. Obviously, the HFWR parachute has better stability than the round parachute.

[figure omitted; refer to PDF]

6.1.2. Inflation Time and Projected Diameter

Table 5 shows the comparison of numerical simulation and wind tunnel tests for the fully inflated time of the HFWR parachute in different cases. Ignoring the straightening time of the suspension lines, the inflation time of numerical simulation is slightly shorter than that of the wind tunnel test. In the actual test situation, the asymmetry caused by the manufacturing process and installation of the test parachute model and the airflow out of sync in the wind tunnel will result in a certain difference between the test and simulation. However, the relative errors within a reasonable range will not affect the consistency of the test and simulation results.

Table 5

Calculated and measured inflation time of the HFWR parachute.

Inflow velocity (m/s)	Numerical inflation time (s)	Test inflation time (s)	Relative errors $\left(\text{\%}\right)$
50	0.0241	0.0290	20.33
280	0.0075	0.0086	14.67

Figure 10 shows the evolution over time of the projected diameter of the canopy under different airflow rates. The curves indicate a clear inflation behavior of the canopy. A large airflow speed would lead to an increase in the projected diameter of the canopy, but the increase is not large.

[figure omitted; refer to PDF]

6.1.3. Rotation Rate

To make the surface of the canopy inclined after fully inflated, the suspension lines are not equal in length, so the FE model of suspension lines cannot be established directly, which should be built up by the folded method. Therefore, in the initial stage of numerical simulation, the suspension lines are not all straightened, which causes the fluctuation of the rotation rate and the opening shock. After the inflation is completed, the parachute gradually reaches the steady state and rotates smoothly around the axis of the parachute. Figure 11 shows the curves of the rotation rate of the canopy with time under different airflow speeds. It shows that the larger the flow speed is, the larger the steady rotation rate of the canopy is and the faster the parachute reaches the steady state.

[figure omitted; refer to PDF]

To better understand the influence of the airflow speed on the steady rotation rate of canopy, a series of cases are calculated and analyzed with different inflow speeds and compared with the wind tunnel test results. There are four cases in all (depicted by A, B, C, and D), and the corresponding inflow speeds are 50 m/s, 200 m/s, 245 m/s, and 280 m/s, respectively. In Table 6, the comparison between numerical results and test results of the rotation rate is given. Both numerical and test results show that the steady rotation rate of canopy increases with the increase of the inflow speed and the relative errors are reasonable. The rate of $\omega /v$ for the HFWR parachute becomes smaller with a higher inflow speed slightly. This is because a large inflow speed would cause the mesh to deform more severely, leading to the canopy gradient decreases.

Table 6

Comparison between numerical results and test results of rotation rate.

Case	Numerical rotation rate (r/s)	Test result of rotation rate (r/s)	Relative errors $\left(\text{\%}\right)$	The rate of (ω/v)
A	8.81	9.21	4.54	0.176
B	30.63	35.5	15.89	0.153
C	34.73	37.33	7.48	0.142
D	37.46	40.83	8.99	0.134

6.2. Structure Dynamics Results

Figures 12–14 show the pressure, the von Mises stress, and the effective plastic strain distribution contours of the fully inflated canopy at a flow rate of 50 m/s, respectively. Because of the permeability and obstruction effect of the fabric canopy and the flow passing the slotted area, the flow speed of outer canopy is significantly higher than that in the inner area, which results in a differential pressure on the surface of the canopy. Therefore, it is obvious that the negative pressure and stress concentration occur at the fringe of canopy connected with the suspension lines, and the maximum stress occurs at the corners of the slotted area.

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

Figure 15 illustrates the curves of stress and strain varying with time for the maximum stress element (10999) at the corner of the gores at different flow speeds. The stress and strain of the element have the same trend and fluctuate slightly with the rotation of the parachute. As the airflow speed rises, the maximum stress level increases. Therefore, in the design process of the parachute, it should be considered to strengthen the fabric in the area with relatively large stress levels by appropriately arranging the reinforcement belts.

[figures omitted; refer to PDF]

6.3. Aerodynamic Analysis

6.3.1. Fluid Domain Characteristics

The entire dynamic process of the fluid domain can be obtained from the computation. Figure 16 shows the pressure contours of the fluid domain at a flow rate of 50 m/s. The variation of the fluid domain pressure with the inflation process of the canopy is given in detail. From the pressure distribution, it is obvious that the airflow gradually gathers under the canopy with the obstruction effect of the canopy, resulting in the rarefaction of airflow near the upper surface of the canopy. Therefore, the pressure level under the canopy is higher than that above it. The pressure difference of the airflow between the canopy surfaces is about 2500 Pa when the canopy is fully inflated at 0.0241 s. After that, the pressure difference decreases gradually until the parachute reaches the steady state and reduces to about 1000 Pa. The pressure difference provides resistance to the parachute.

[figures omitted; refer to PDF]

To better understand the fluid domain characteristics and the evolvement mechanism of the vortexes, the FSI results of the fluid domain are further analyzed. Velocity contours and streamline of the fluid domain at different times are shown in Figure 17. In this case, the inflow speed is 50 m/s. As the airflow passes the canopy, the canopy expands from the initial folded state to resist the fluid under the action of aerodynamic pressures, which makes the airflow velocity at the fringe of the canopy significantly higher than that near the surface of canopy. Besides, due to the presence of the vent on the top of the canopy, the velocity of the airflow around the aperture is larger than that of the surrounding airflow. Then, symmetrical vortexes gradually develop above the canopy. With the rise of airflow and the rotation of the parachute, the vortexes grow and move separately upward in the airflow direction. After the parachute reaches the steady state, the fluid around the parachute exhibits periodic vortex shedding. A large number of vortexes produce much vortex resistance, which makes the HFWR parachute have greater resistance than the round parachute in the same nominal diameter, as shown in Figure18(b).

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

6.3.2. Opening Shock and Resistance Characteristics

To analyze the influence of the airflow velocity on the opening load of the parachute, the opening load varying with time at low speed (50 m/s and 60 m/s) and subsonic speed (245 m/s and 280 m/s) is simulated. In Figure 19, the maximum opening load occurs when the canopy is completely open, which is consistent with the variation rule of the inflation process in an infinite mass scenario. Subsequently, the opening load gradually decreases and eventually reaches a steady value with fluctuating slightly as canopy shrinks slightly and rotates. The opening load increases with the rise of airflow speed, which is more obvious with a larger airflow speed.

The numerical maximum opening shock of the HFWR parachute for all cases is transformed into the resistance characteristic by (23), that is,$$\begin{array}{}\text{(23)}& F_{\max }=\frac{1}{2}\rho v^2\left(CA\right)K_D,\end{array}$$where ρ denotes the fluid density, $v$ signifies flow speed, $\left(CA\right)$ denotes resistance characteristic of the fully open parachute, and $K_D$ represents opening shock coefficient at infinite mass [35].

To verify the validity of the simulation results, the results of the resistance characteristic and the maximum opening shock are compared with the wind tunnel test results, as shown in Table 7. Assume that the characteristic length of the Reynolds number is the nominal diameter of the canopy.

Table 7

Comparison of maximum opening shock and resistance characteristic results.

Inflow velocity (m/s)	Reynolds number ($\times {10}^6$)	Numerical maximum opening shock (N)	Test maximum opening shock (N)	Numerical resistance characteristic (m2)	Test resistance characteristic (m2)	Relative errors ($\text{\%}$)
50	2.58	1483.5	1226	0.88	0.80	9.1
200	10.33	15173	14257	0.56	0.53	5.4
245	12.65	20824	18265	0.51	0.45	11.8
280	14.46	26200	21358	0.50	0.44	12.0

Table 7 shows that the numerical results are consistent with the wind tunnel test results. The numerical results are slightly larger than that of the test, which is caused by the fact that the numerical model of the HFWR parachute does not fully consider the influence of the damping force and the friction between the canopy and suspension lines; however, the relative errors are acceptable. Moreover, both the numerical results and the test results decrease with increasing flow rate, and the variation with the Reynolds number is small. Therefore, the simulation results have high reliability and can serve as a reference.

In the process of parachute development, to ensure the reliability and effectiveness of the parachute-payload system, it is necessary to reduce the value of the maximum opening shock to avoid damage to the structure of the parachute. Furthermore, it is necessary to increase the resistance coefficient of parachute as much as possible to ensure that the parachute has a smaller size. Therefore, the numerical results of the opening load of the HFWR parachute and round parachute are compared. As shown in Figure 18, the maximum opening shock and steady resistance of the round parachute are much larger than those of the HFWR parachute with the same structural diameter, while the maximum opening shock and steady resistance of the round parachute are smaller than those of the HFWR parachute with the same nominal diameter.

[figures omitted; refer to PDF]

Details of the comparison of the maximum opening shock and steady resistance of the HFWR parachute and the round parachute at the airflow velocity of 50 m/s and 280 m/s are listed in Table 8. With the same structural diameter, the maximum opening shock and steady resistance of the round parachute are 2.08 and 1.56 times that of the HFWR parachute at the airflow speed of 50 m/s and 2.67 and 2.40 times that of the HFWR parachute at the airflow speed of 280 m/s, respectively. Consequently, in the case of higher speed, the use of a HFWR parachute can greatly reduce the opening shock. Besides, with the same nominal diameter, the maximum opening shock ratio and steady resistance ratio of the round parachute and the HFWR parachute are about 0.66 and 0.83 at the airflow speed of 50 m/s and 0.79 and 0.76 at the airflow speed of 280 m/s, respectively. Therefore, the HFWR parachute has a higher resistance coefficient than the round parachute.

Table 8

Comparison of opening shock and steady resistance of two types of parachutes.

Flow velocity (m/s)	Same structural diameter		Same nominal diameter
	50	280	50	280
Max of opening shock of round parachute (N)	3083.1	69900	976.39	20640
Max of opening shock of HFWR parachute (N)	1483.5	26200	1483.5	26200
Max of opening shock ratio	2.08	2.67	0.66	0.79
Steady resistance of round parachute (N)	877.7	38700	464.94	12242
Steady resistance of HFWR parachute (N)	561.70	16135	561.70	16135
Steady resistance ratio	1.56	2.40	0.83	0.76

7. Conclusion

In this paper, the inflation process of a new type of HFWR parachute for the stand-off aviation weapon is studied. The FSI finite element model of the parachute and the surrounding fluid domain is established. The transient dynamic behavior and fluid domain evolution of the HFWR parachute in an infinite mass scenario are investigated. The variation of rotation rate and opening shock of the HFWR parachute at different flow speeds is analyzed. The FSI simulation results are verified by the wind tunnel test data. Besides, the comparison of the opening performances of the HFWR parachute and round parachute is conducted. Based on the results obtained in this paper, the following conclusions are summarized:

(1)

The HFWR parachute can quickly reach the steady state with good inflation performance, and there is no obvious “breathing phenomenon” during the inflation process.

(2)

The stress concentration appears at the connection of the fringe of canopy and the suspension lines, so it is necessary to strengthen the fabric in the area with relatively large stress levels.

(3)

A large number of persistent vortexes distributing symmetrically along the parachute axis are generated above the HFWR parachute, which periodically sheds off, making the HFWR parachute possess greater resistance coefficient.

(4)

With the same structural diameter, the maximum opening shock of the HFWR parachute is much smaller than that of the round parachute. With the same nominal diameter, the HFWR parachute has a higher resistance coefficient than the round parachute. Moreover, the stability of the HFWR parachute is better than that of the round parachute.

To sum up, the HFWR parachute has a good dynamic performance with many excellent characteristics, such as small size and opening shock, large resistance coefficient, and good stability, making the HFWR parachute more suitable for the application of the stand-off aviation weapon system. The numerical model and simulation methods used in this study can effectively guide the design and engineering applications of the rotating parachute. Nevertheless, the environmental factors considered in the study are still relatively simple. In future work, more advanced numerical simulation methods can be explored to study the influence of more complex environmental factors on the stability of rotating parachute.

Authors’ Contributions

Y.L. was responsible for conceptualization, data curation, formal analysis, investigation, methodology, software, visualization, and original draft preparation. C.J. acquired funding and supervised the study. C.J. and M.L. were responsible for project administration. Y.L. and L.M. validated the study. Y.L. and L.M. were responsible for the wind tunnel test. Y.L., L.M., and C.J. reviewed and edited the manuscript.

Acknowledgments

The authors are grateful to the support of the National Natural Science Foundation for this research. This research was funded by the National Natural Science Foundation of China (grant no. 11672042).

Glossary

Abbreviations

2D:Two-dimensional

3D:Three-dimensional

ALE:Arbitrary Lagrange–Euler

FSI:Fluid-structure interaction

FE:Finite element

HFWR parachute:High-damping four-winged rotating parachute

Max:Maximum.

$\beta$:Angle between the folded gore

$L_{\text{arc}}$:Top edge of a piece of gore

R:Side edge of the canopy

H:Edge of folded line

D_a:Aperture diameter

D_c:Structural diameter

L:Half the length of the bottom edge of the canopy

${\Omega }^s$:Domain of the structure of canopy

$\partial {\Omega }^s$:Boundary of the structure

${\rho }_s$:Density of the structure

$\overline{\overline{{\sigma }_s}}$:Cauchy stress of the structure

e:Specific internal energy

$\left(\text{d}\stackrel{⟶}{v}/\text{d}t\right)$:Acceleration

$f_s$:Force on the structure

${\stackrel{⟶}{n}}_s$:Exterior unit normal of the structure boundary

$\stackrel{⟶}{D}$:Displacement vector of the structure

$\stackrel{⟶}{T}$:Stress vector of the structure

$\partial {\Omega }_1^s$:Displacement boundary of the structure

$\partial {\Omega }_2^s$:Traction boundary of the structure

σ:Stress

$\nu$:Poisson’s ratio

E:Elastic modulus

${\tau }_{12}$:Shear stress

${\epsilon }_{12}$:Shear strain

$G_{12}$:Shear elasticity

α:Nonlinear coefficient

$F_t$:Force of cable

$F_0$:Initial force of cable

$\Delta L$:The change in length of cable

K:Stiffness

$l_{\text{current}}$:Current length of the element

$l_0$:Initial length of the element

$l_{\text{offset}}$:Offset of the element length

$A_c$:Cross-sectional area of the cable

${\Omega }^f$:Spatial domain

$\partial {\Omega }^f$:Boundary of the spatial domain

${\rho }_f$:Density of the fluid

f:Body force

$\stackrel{⟶}{v}$:Fluid velocity

$\stackrel{⟶}{w}$:Mesh velocity

$\overline{\overline{{\sigma }_f}}$:Total Cauchy stress of the fluid

p:Pressure

$Id$:Unit second order tensor

μ:Dynamic viscosity of the fluid

$T_f$:Air temperature

$Cp$:Specific heat capacities at constant pressure

$Cv$:Specific heat capacities at constant volume

$\stackrel{⟶}{g}$:Displacement vector of the fluid

$\stackrel{⟶}{h}$:Stress vector of the fluid

${\stackrel{⟶}{n}}_f$:Exterior unit normal of the fluid boundary

$\partial {\Omega }_1^f$:Displacement boundary of the fluid

$\partial {\Omega }_2^f$:Traction boundary of the fluid

$\stackrel{⟶}{d}$:Penetration depth

${\stackrel{⟶}{v}}_{\text{rel}}$:Relative velocity of the slave node and master node

${\stackrel{⟶}{v}}_s$:Structure velocity

${\stackrel{⟶}{v}}_f$:Fluid velocity

r:Normal direction of the shell element

ε:Porosity

$a\left(\mu ,\epsilon \right)$:Viscous coefficient

$b\left(\rho ,\epsilon \right)$:Inertia coefficient

D:Characteristic length

$V$:Volume of the canopy

S:Permeated surface of the canopy

F:Porous coupling force

$F_s$:Coupling force to a structure node

$F_f$:Fluid force

$N_i$:Shape function

$v$:Flow speed

ρ:Flow density

CA:Resistance characteristic

$K_D$:Opening shock coefficient at infinite mass.

References

[1] J. W. Purvis, "Theoretical analysis of parachute inflation including fluid kinetics," Journal of Aircraft, vol. 19 no. 4, pp. 290-296, DOI: 10.2514/3.57392, 1982.

[2] J. W. Purvis, "Numerical prediction of deployment, initial fill, and inflation of parachute canopies," Proceedings of the 8th Aerodynamic Decelerator and Balloon Technology Conference,DOI: 10.2514/6.1984-787, .

[3] K. R. Stein, R. J. Benney, "Parachute inflation: a problem in aeroelasticity," 1994. Technical report

[4] R. J. Benney, K. R. Stein, "Computational fluid-structure interaction model for parachute inflation," Journal of Aircraft, vol. 33 no. 4, pp. 730-736, DOI: 10.2514/3.47008, 1996.

[5] L. Yu, X. Ming, "Study on transient aerodynamic characteristics of parachute opening process," Acta Mechanica Sinica, vol. 23 no. 6, pp. 627-633, DOI: 10.1007/s10409-007-0112-3, 2007.

[6] K. Stein, R. Benney, V. Kalro, T. E. Tezduyar, J. Leonard, M. Accorsi, "Parachute fluid–structure interactions: 3-D computation," Computer Methods in Applied Mechanics and Engineering, vol. 190 no. 3-4, pp. 373-386, DOI: 10.1016/s0045-7825(00)00208-5, 2000.

[7] V. Kalro, T. E. Tezduyar, "A parallel 3D computational method for fluid–structure interactions in parachute systems," Computer Methods in Applied Mechanics and Engineering, vol. 190 no. 3-4, pp. 321-332, DOI: 10.1016/s0045-7825(00)00204-8, 2000.

[8] T. E. Tezduyar, S. Sathe, "Modelling of fluid–structure interactions with the space–time finite elements: solution techniques," International Journal for Numerical Methods in Fluids, vol. 54 no. 6-8, pp. 855-900, DOI: 10.1002/fld.1430, 2007.

[9] T. E. Tezduyar, S. Sathe, T. Cragin, "Modelling of fluid–structure interactions with the space–time finite elements: arterial fluid mechanics," International Journal for Numerical Methods in Fluids, vol. 54 no. 6-8, pp. 901-922, DOI: 10.1002/fld.1443, 2007.

[10] T. E. Tezduyar, S. Sathe, M. Schwaab, J. Pausewang, J. Christopher, J. Crabtree, "Fluid-structure interaction modeling of ringsail parachutes," Computational Mechanics, vol. 43 no. 1, pp. 133-142, DOI: 10.1007/s00466-008-0260-8, 2008.

[11] T. E. Tezduyar, S. Sathe, J. Pausewang, M. Schwaab, J. Christopher, J. Crabtree, "Interface projection techniques for fluid-structure interaction modeling with moving-mesh methods," Computational Mechanics, vol. 43 no. 1, pp. 39-49, DOI: 10.1007/s00466-008-0261-7, 2008.

[12] K. Takizawa, S. Wright, C. Moorman, T. E. Tezduyar, "Fluid-structure interaction modeling of parachute clusters," International Journal for Numerical Methods in Fluids, vol. 65 no. 1–3, pp. 286-307, DOI: 10.1002/fld.2359, 2011.

[13] K. Takizawa, T. Spielman, T. E. Tezduyar, "Space-time FSI modeling and dynamical analysis of spacecraft parachutes and parachute clusters," Computational Mechanics, vol. 48 no. 3, pp. 345-364, DOI: 10.1007/s00466-011-0590-9, 2011.

[14] K. Takizawa, M. Fritze, D. Montes, T. Spielman, T. E. Tezduyar, "Fluid-structure interaction modeling of ringsail parachutes with disreefing and modified geometric porosity," Computational Mechanics, vol. 50 no. 6, pp. 835-854, DOI: 10.1007/s00466-012-0761-3, 2012.

[15] K. Takizawa, T. E. Tezduyar, "Computational methods for parachute fluid-structure interactions," Archives of Computational Methods in Engineering, vol. 19 no. 1, pp. 125-169, DOI: 10.1007/s11831-012-9070-4, 2012.

[16] K. Takizawa, T. Spielman, C. Moorman, T. E. Tezduyar, "Fluid-structure interaction modeling of spacecraft parachutes for simulation-based design," Journal of Applied Mechanics, vol. 79 no. 1,DOI: 10.1115/1.4005070, 2012.

[17] K. Takizawa, T. E. Tezduyar, J. Boben, N. Kostov, C. Boswell, A. Buscher, "Fluid-structure interaction modeling of clusters of spacecraft parachutes with modified geometric porosity," Computational Mechanics, vol. 52 no. 6, pp. 1351-1364, DOI: 10.1007/s00466-013-0880-5, 2013.

[18] K. Takizawa, T. E. Tezduyar, C. Boswell, R. Kolesar, K. Montel, "FSI modeling of the reefed stages and disreefing of the Orion spacecraft parachutes," Computational Mechanics, vol. 54 no. 5, pp. 1203-1220, DOI: 10.1007/s00466-014-1052-y, 2014.

[19] B. Tutt, A. Taylor, "The use of LS-DYNA to simulate the inflation of a parachute canopy," Proceedings of the 18th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar, pp. 56-64, DOI: 10.2514/6.2005-1608, .

[20] J. Wang, N. Aquelet, B. Tutt, I. Do, M. Souli, "Porous Euler-Lagrange coupling: application to parachute dynamics," Proceedings of the 9th International LS-DYNA Users Conference, .

[21] N. Aquelet, B. Tutt, "Euler-Lagrange coupling for porous parachute canopy analysis," The International Journal of Multiphysics, vol. 1 no. 1, pp. 53-68, DOI: 10.1260/175095407780130517, 2007.

[22] J. Lingard, M. Darley, "Simulation of parachute fluid structure interaction in supersonic flow," Proceedings of the 18th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar, pp. 47-55, DOI: 10.2514/6.2005-1607, .

[23] J. S. Lingard, M. Darley, J. C. Underwood, "Simulation of mars supersonic parachute performance and dynamics," Proceedings of the 19th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar, pp. 67-77, .

[24] B. Tutt, R. Scott, G. Noetscher, R. Charles, "Finite mass simulation techniques in LS-DYNA," Proceedings of the 21st AIAA Aerodynamic Decelerator Systems Technology conference and Seminar,DOI: 10.2514/6.2011-2592, .

[25] H. Cheng, L. Yu, W. Rong, H. Jia, "A numerical study of parachute inflation based on a mixed method," Aviation, vol. 16 no. 4, pp. 115-123, DOI: 10.3846/16487788.2012.753676, 2012.

[26] X. Gao, Q. Zhang, Q. Tang, "Transient dynamic modeling and analysis of complex parachute inflation with fixed payload," Journal of Aerospace Engineering, vol. 28 no. 4,DOI: 10.1061/(asce)as.1943-5525.0000430, 2014.

[27] X. Gao, Q. Zhang, Q. Tang, "Fluid-structure interaction analysis of parachute finite mass inflation," International Journal of Aerospace Engineering, vol. 2016,DOI: 10.1155/2016/1438727, 2016.

[28] X. Gao, Q. Zhang, Q. Tang, "Fluid-structure interaction and flight dynamics analysis of parachute-payload system during uncontrolled airdrop process," Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, vol. 232 no. 13, pp. 2499-2512, DOI: 10.1177/0954410017718565, 2018.

[29] X. Gao, Q. Zhang, Q. Tang, "Numerical modelling of mars supersonic disk-gap-band parachute inflation," Advances in Space Research, vol. 57 no. 11, pp. 2259-2272, DOI: 10.1016/j.asr.2016.03.010, 2016.

[30] T. Belytschko, J. I. Lin, T. Chen-Shyh, "Explicit algorithms for the nonlinear dynamics of shells," Computer Methods in Applied Mechanics and Engineering, vol. 42 no. 2, pp. 225-251, DOI: 10.1016/0045-7825(84)90026-4, 1984.

[31] N. Aquelet, J. Wang, B. A. Tutt, I. Do, H. Chen, "Euler-Lagrange coupling with deformable porous shells," Proceedings of the ASME 2006 Pressure Vessels and Piping/ICPVT-11 Conference, pp. 285-290, .

[32] D. J. Benson, "Computational methods in Lagrangian and eulerian hydrocodes," Computer Methods in Applied Mechanics and Engineering, vol. 99 no. 2-3, pp. 235-394, DOI: 10.1016/0045-7825(92)90042-i, 1992.

[33] T. J. R. Hughes, W. K. Liu, T. K. Zimmermann, "Lagrangian-eulerian finite element formulation for incompressible viscous flows," Computer Methods in Applied Mechanics and Engineering, vol. 29 no. 3, pp. 329-349, DOI: 10.1016/0045-7825(81)90049-9, 1981.

[34] S. Ergun, "Fluid flow through packed columns," Chemical Engineering Progress, vol. 48, pp. 89-94, 1952.

[35] L. R. Wang, Parachute Theory and Application, 1997.