Numerical study of lateral wind effect on parachute dropping based on finite element method

H. Cheng*, L. Wei**, B. Zhou***, Q. He****, P. Xiao***** *Civil Aviation Flight University of China, Guanghan 618307, China, E-mail: chenghanstorm@sina.com **Civil Aviation Flight University of China, Guanghan 618307, China, E-mail: 597614657@qq.com ***Civil Aviation Flight University of China, Guanghan 618307, China, E-mail: zhoubinzhoubin@163.com ****Civil Aviation Flight University of China, Guanghan 618307, China, E-mail: hqcq@mail.nwpu.edu.cn *****Civil Aviation Flight University of China, Guanghan 618307, China, E-mail: poaix@163.com


Introduction
Parachute is widely used in aerospace, aviation, weapon, and other fields.Whether the parachute works properly or not would affect the safety of personnel or system.There were many serious accidents occurred in the history of aerospace and aviation because of the parachute failure [1].However, there are many factors affect the parachute working, such as canopy material, fabric permeability, working environment and so on.In the above factors, the extra wind is complicated and changeable which affected by topography, altitude, temperature and other factors.Therefore parachute dropping has high randomicity and difficulty in data collection, and it is difficult to investigate the extra wind effect on parachute working by airdropping experiments.While the numerical method with high efficiency and repeatability gradually becomes an important approach for parachute design and research.
The parachute working is a complicated transient and nonlinear process.Although it is very difficult to investigate this coupling process based on numerical method, many scholars applied different methods or models to research.Purvis [2] simplified the parachute structure and flow field and realized two-dimensional coupling simulation.Stein, Benny and et al. [3] proposed CFD/MSD coupling model, Yu [4] and other scholars applied this model to simulate parachute opening in infinite mass situation.Kim [5] used IB (Immersed Boundary) method, and calculated parachute opening under low Reynolds number.Tutt [6] used ALE (Arbitrary Eulerian Lagrangian) method realized parachute working process, and this method was the most widely used in engineering.Takizawa [7] used SSTFSI (Stabilized Space Time Fluid Structure Interaction) calculated the structure and flow field of a fully inflated parachute.
However, the parachute flight characteristics were not considered and the effect of extra wind was ignored in most researches above, which has limited reference value for engineering application.The relevant research works about the extra wind field effect on parachute dropping have not yet been reported therefore investigating this effect has important research value.In this paper, the lateral wind effect on parachute dropping was investigated.The flow field and canopy structure were discretized by finite elements.That the velocity of lateral wind was defined as boundary condition was used to simulate the lateral wind field.The coupling between canopy and air was realized based on penalty function.To reduce the calculation amount, the graphical deformation technology was used to realize the fluid meshes motion.While the boundary condition was not changed, the stable wind field environment could be maintained.

Governing and discretization equations
In this paper, the explicit finite element method based on Updated Lagrangian (UL) scheme is used to calculate the entire parachute dropping.Both the canopy structure and flow field are discretized by Lagrangian meshes.It is unnecessary to solve the mass conservation equation as finite volume method because the Lagrangian meshes will deform with the deformation of materials.In addition, the temperature field is not considered in calculation.Therefore, only the momentum equation needs to be solved, and this equation based on UL scheme described is shown in Eq. (1): where v is the velocity vector, σ is the stress vector, b is the body force vector, ρ is the density.The Eq. ( 1) is discretized in spatial and shown in Eq. (2): 0, where BIj is geometry matrix, σji is stress, The central difference scheme is used in time domain discretization as shown in Eq. ( 3): where d is displacement.

Constitutive equations of materials
The deformation features of fabrics are small strain and large rotation during parachute dropping.Therefore the canopy is described by Kirchhoff material constitutive equation in this paper.Meanwhile, the canopy is divided by three-node triangular shell elements.The constitutive equation of canopy in Voigt notation can be simplified as Eq. ( 4):


  where Sij is the second Piola-Kirchhoff stress, Eij is Green stress, E is elastic modulus, υ is Poisson ratio.The lines and reinforcing belt are described by one-dimensional linear elastic constitutive equation as shown in Eq. ( 5): In this paper, both the initial dropping velocity and lateral wind velocity are less than 0.3 Mach number.Therefore, the air could be defined as uncompressible material.The volume and density are constant during air mass moving, which is called as equal volume motion.The stress description is divided into deviatoric part and hydrostatic part as shown in Eq. ( 6): where μ is dynamic viscosity coefficient, and p is hydrostatic pressure.The Eq. ( 6) must be used together with the ideal gas state equation.

Coupling calculation based on contact algorithm
In order to reflect the fabric permeability, the Ergun equation is used to calculate the coupling force: where Δp is differential pressure, a is linear resistance coefficient, b is quadratic resistance coefficient, e is fabric thickness; n is normal vector; v is fluid velocity vector through the fabric.
Then the coupling force couple f derived from the Eq. ( 7) is applied to both the fluid and fabric in opposite directions to satisfy force equilibrium, and the coupling force

Motion and reconstruction of fluid meshes
The fluid meshes described by Lagrangian meshes would be distorted after each time step.In order to reduce the computation amount at the same time, the fluid meshes need to be moved and reconstructed after each time step.Here, Fig. 1 is used to illustrate the fluid meshes processing.; The local coordinate system will move after each time step, and a transformation matric T can be obtained according to the displacement [8].Therefore, the new homogeneous coordinate of each node on fluid meshes can be calculated based on Eq. ( 9): where  is the homogeneous coordinate after the displacement,   is that before the displacement.
According to the coordinates of fluid meshes before and after the displacement, the fluid meshes moving velocity v ( ˆΔ / Δt  vx ) can be obtained.Then the convection velocity c ( =c v v , where, v is materials moving velocity) which based on flow field meshes as a reference can be calculated.Therefore the materials velocity v in Eq. ( 1) is replaced by convection velocity c.
The distorted fluid meshes are reconstructed by solving the Laplace Eq. ( 11): .
After the meshes reconstruction, the convection terms ϕ (such as mass, momentum and other flow field information) are updated based on MUSCL (Monotone Upwind Schemes for Conservation Laws) scheme [9] with second order accuracy in this paper.

Example
The typical C9 parachute was taken as the calculation object, and the structure and material parameters of C9 parachute could be found in the reference [10].The finite element model of parachute was established according to those parameters.The canopy was discretized by three-node triangular elements (1,4000).The lines and reinforcing belts were by bar elements (1,932).As shown in Fig. 2, a local coordinate was established on payload model.

Comparison between Model A and test results
It can be found that the calculation results were basically consistent with the dropping test [11] (Fig. 4).The entire process can be divided into three phases.
The pre-inflation phase (0.0-0.41): with the opening of canopy, the aerodynamic deceleration area and dynamic load increased.When the vent of the parachute fully opened, the first peak appeared in both numerical and experimental results.
The fully inflation phase (0.41-1.0): after the vent was fully opened, the decelerating effect was significant.Therefore the dynamic load began to decline after the first peak.When the lower part of the canopy expanded, the deceleration area reached the maximum value and the second peak appeared.
The stable dropping phase (after 1.0): when the deceleration shape was stable, the dynamic load decreased followed the velocity's decreasing.
We also found that the dimensionless load curve of calculation was slightly higher than the test curve, because the contrail declining angle in this work was 90° while the test angle was less than 90° in actual dropping.It is proved that the dynamic load of 90° is greater than other angles based on a large number of tests.

Analysis of deceleration characteristics
The acceleration and velocity changing of payload in windless environment (Model A) and lateral wind environment (Model B) are shown in Fig. 5.The changing trend is almost the same.In stable dropping phase, both Model A and Model B trended to 0 m/s 2 , and the velocity maintained at about 6 m/s gradually.The dropping velocities were consistent to the actual stable dropping velocity.In windless environment, Model A would have a small amount of displacement both in x and y directions in stable dropping phase due to 'breathing' phenomenon.The displacement trajectory was approximately a line.However, Model B displaced in x direction in pre-inflation phase, which is not obvious due to the canopy was not fully expanded.With the canopy expanding gradually, the x direction projection area of canopy increased, and the lateral displacement was obvious under the action of lateral wind.
It can be found in Fig. 6 that the lateral displacement was about 8 m.   7 shows the equivalent stress contour and velocity vector.Both Model A and B had the same tendency in shape and stress changing.But Model B had a certain deflection angle caused by lateral wind effect.In preinflation phase, canopies bottom opened firstly, and then expanded uniformly (Fig. 7, a).After the vents opened fully, the canopies appeared a 'squid' shape, and the canopies began to inflate from top to bottom.In this moment, the maximum stress value concentrated on the top part, which might be the weakest part during the entire dropping process (Fig. 7, b-c).The acceleration curves appeared the first peak and the deceleration effect was obvious (Fig. 5).That the canopies fully expanded indicated the end of fully inflation phase (Fig. 7, d).The canopies provided the maximum aerodynamic drag, which caused the second acceleration peak (Fig. 5).Then, the deceleration areas would continue to increase because of inertia (Fig. 7, e).At last, the canopies appeared specific 'breathing' phenomenon in stable dropping phase, and the acceleration and velocity trended to be stable (Fig. 5).
However, the flow field changes were different.In windless environment, the air flowed out from inside of canopy and formed a symmetric vortex at the canopy bottom.With the canopy inflating, the outside vortex gradually rose to the canopy top and maintained a stable vortex structure in stable dropping phase.While, the vortex structure change in lateral wind environment was complicated and irregular.In pre-inflation phase, the vortex structure was not obvious (Fig. 7, a).After the vent fully opened, a vortex appeared at leeward side of canopy top and the vortex on upwind side was at the canopy bottom (Fig. 7, b-c).When the canopy fully expanded, the vortex on upwind side rose to the canopy top and was weaker than the leeward side's vortex (Fig. 7, d-e).

Conclusions
In this paper, the C9 parachute working process in windless and lateral wind environments were calculated.The structure change and flow field change were obtained.
The experimental data was used to verify the accuracy of numerical method.According to results analysis, the deceleration characteristics were not affected by lateral wind but the trajectory would be affected greatly, which was consistent with practical engineering.The results could reflect the general working laws of parachute in practice.
The numerical method used in this paper could be a reference for parachute design.
force, NI and NJ are shape functions, Ji v is acceleration, int f is internal force matrix, ext f is external force matrix, M is mass matrix, a is acceleration matrix.
couple f are taken as a part of external force ext f in Eq. (2).

Fig. 1
Fig. 1 Fluid meshes' motion (the cyan meshes represent structure elements, the wireframe meshes represent fluid elements, and red arrows are velocity boundary) Three noncollinear nodes are chosen on structure, and the coordinates of these nodes are xA, xB and xC.A local coordinate system (red) can be defined according to these three nodes, and the axis's vectors are shown in Eq. (8).      .

Fig. 2
Fig. 2 Parachute model and local coordinate system While the flow field domain was divided by hexahedral meshes (1,150,000).The fluid meshes and struc-

Fig. 4
Fig. 4 Comparison of dynamic load (where, t is time, tf is the time of parachute fully inflated, F is dynamic load, (CDS)o qs is stable resistance)

Fig. 5
Fig. 5 Acceleration and velocity of payload 5.2.Analysis of payload displacement Fig. 6 shows the displacements of Model A and B.In windless environment, Model A would have a small amount of displacement both in x and y directions in stable dropping phase due to 'breathing' phenomenon.The displacement trajectory was approximately a line.However, Model B displaced in x direction in pre-inflation phase, which is not obvious due to the canopy was not fully expanded.With the canopy expanding gradually, the x direction projection area of canopy increased, and the lateral displacement was obvious under the action of lateral wind.

Fig. 6
Fig. 6 Displacement of payload 5.3.Analysis of structure and flow field change

Fig.
Fig.7shows the equivalent stress contour and velocity vector.Both Model A and B had the same tendency in shape and stress changing.But Model B had a certain deflection angle caused by lateral wind effect.In preinflation phase, canopies bottom opened firstly, and then expanded uniformly (Fig.7, a).After the vents opened fully, the canopies appeared a 'squid' shape, and the canopies began to inflate from top to bottom.In this moment, the maximum stress value concentrated on the top part, which might be the weakest part during the entire dropping process (Fig.7, b-c).The acceleration curves appeared the first peak and the deceleration effect was obvious (Fig.5).That the canopies fully expanded indicated the end of fully inflation phase (Fig.7, d).The canopies provided the maximum aerodynamic drag, which caused the second acceleration peak (Fig.5).Then, the deceleration areas would continue to increase because of inertia (Fig.7, e).At last, the canopies appeared specific 'breathing' phenomenon in stable dropping phase, and the acceleration and velocity trended to be stable (Fig.5).However, the flow field changes were different.In windless environment, the air flowed out from inside of canopy and formed a symmetric vortex at the canopy bottom.With the canopy inflating, the outside vortex gradually rose to the canopy top and maintained a stable vortex structure in stable dropping phase.While, the vortex structure change in lateral wind environment was complicated and irregular.In pre-inflation phase, the vortex structure was not obvious (Fig.7, a).After the vent fully opened, a vortex appeared at leeward side of canopy top and the vortex on upwind side was at the canopy bottom (Fig.7, b-c).When the canopy fully expanded, the vortex on upwind side rose to the canopy top and was weaker than the leeward side's vortex (Fig.7, d-e).