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

Cheng, H. ; Wei, L. ; Zhou, B. 等

1. 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.

2. Model development

2.1. 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):

[rho]Dv/Dt = [nabla] [sigma] + [rho]b, (1)

where v is the velocity vector, [sigma] is the stress vector, b is the body force vector, [rho] is the density.

The Eq. (1) is discretized in spatial and shown in Eq. (2):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2)

where [B.sub.Ij] is geometry matrix, [[sigma].sub.ji] is stress, [[bar.t].sub.i] is external force, [N.sub.I] and [N.sub.J] are shape functions, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is acceleration, [f.sup.int] is internal force matrix, [f.sup.ext] is external force matrix, M is mass matrix, a is acceleration matrix.

The central difference scheme is used in time do-main discretization as shown in Eq. (3):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3)

where d is displacement.

2.2. 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):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4)

where [S.sub.ij] is the second Piola-Kirchhoff stress, [E.sub.ij] is Green stress, E is elastic modulus, [upsilon] is Poisson ratio.

The lines and reinforcing belt are described by one-dimensional linear elastic constitutive equation as shown in Eq. (5):

[[sigma].sub.x] = E [[epsilon].sub.x]. (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 hydro-static part as shown in Eq. (6):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (6)

where [mu] is dynamic viscosity coefficient, and p is hydrostatic pressure. The Eq. (6) must be used together with the ideal gas state equation.

2.3. Coupling calculation based on contact algorithm

In order to reflect the fabric permeability, the Ergun equation is used to calculate the coupling force:

[delta]p = [avn + b[(vn).sup.2]]e, (7)

where [delta]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 [f.sub.couple] derived from the Eq. (7) is applied to both the fluid and fabric in opposite directions to satisfy force equilibrium, and the coupling force [f.sub.couple] are taken as a part of external force [f.sup.ext] in Eq. (2).

2.4. 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.

[FIGURE 1 OMITTED]

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).

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (8)

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):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (9)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the homogeneous coordinate after the displacement, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is that before the displacement.

According to the coordinates of fluid meshes before and after the displacement, the fluid meshes moving velocity [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] can be obtained. Then the convection velocity c ([MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], 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.

[rho][v.sub.i,t] + [rho][v.sub.i,j][c.sub.j]=[[sigma].sub.ji,j] + [rho][b.sub.i]. (10)

The distorted fluid meshes are reconstructed by solving the Laplace Eq. (11):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (11)

After the meshes reconstruction, the convection terms [phi] (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.

3. 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.

[FIGURE 2 OMITTED]

While the flow field domain was divided by hexahedral meshes (1,150,000). The fluid meshes and structure meshes interspersed with each other. The specific size is shown in Fig. 3.

[FIGURE 3 OMITTED]

The C9 parachute dropping in windless environment (Model A) and lateral wind (3 m/s) condition (Model B) were calculated respectively. The other working conditions are shown in Table 1.

4. 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[degrees] while the test angle was less than 90[degrees] in actual dropping. It is proved that the dynamic load of 90[degrees] is greater than other angles based on a large number of tests.

[FIGURE 4 OMITTED]

5. Analysis of numerical results

5.1. 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/s2, and the velocity maintained at about 6 m/s gradually. The dropping velocities were consistent to the actual stable dropping velocity.

[FIGURE 5 OMITTED]

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. It can be found in Fig. 6 that the lateral displacement was about 8 m.

[FIGURE 6 OMITTED]

5.3. Analysis of structure and flow field change

Fig. 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).

[FIGURE 7 OMITTED]

6. 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.

Acknowledgements

The investigations described in this paper are a part of the research project No F2015KF03 realized at Civil Aviation Flight Technology and Flight Safety Research Base.

References

(1.) Xia, G.; Cheng, W.K.; Qin Z.Z. 2002. Case study of main parachute malfunction in aerospace recovery, Spacecraft Recovery & Remote Sensing 23(4): 4-8 (in Chinese).

(2.) Purvis, J.W. 1981. Theoretical analysis of parachute inflation including fluid kinetics, AIAA Report 1981-1925.

(3.) Stein, K.R.; Benney, R.J.; Steeves, E.C. 1993. A computational model that couples aerodynamic structural dynamic behavior of parachutes during the opening process, NASA-ADA264115.

(4.) Yu L.; Ming, X. 2007. Study on transient aerodynamic characteristics of parachute opening process, Acta Mechanica Sinica 23(6): 627-633. http://dx.doi.org/10.1016/j.cja.2014.02.021.

(5.) Kim, Y.S; Peskin, C.S. 2009. 3-D Parachute simulation by the immersed boundary method, Computers and Fluids 38: 1080-1090. http://dx.doi.org/10.1016/j.compfluid.2008.11.002.

(6.) Tutt, B.A; Taylor, A.P. 2005. The use of LS-DYNA to simulate the inflation of a parachute canopy, AIAA-2005-1608.

(7.) Takizawa, K.; Spielman, T; Moorman, C.; Tezduyar, T.E. 2012. Fluid structure interaction modeling of spacecraft parachutes for simulation-based design, Journal of Applied Mechanics 79: 1-9. http://dx.doi.org/10.1115/1.4005070.

(8.) Bao, H.J.; Jin, X.G.; Pen, Q.S. 2000. Basic Algorithms of Computer Animation. Hangzhou, Zhejiang University publishing house, 20 p. (in Chinese).

(9.) Souli, M.; Ouahsine, A.; Lewin, L. 2000. ALE formulation for fluid-structure interaction problems, Computer Methods in Applied Mechanics and Engineering 190: 659-675. http://dx.doi.org/10.1016/S0045-7825(99)00432-6.

(10.) Cheng, H.; Yu, L.; Rong, W.; Jia, H. 2012. A numerical study of parachute inflation based on a mixed method, Aviation 16(4): 115-123. http://dx.doi.org/10.3846/16487788.2012.753676.

(11.) Ewing, E.G; Knacke, T.W; Bixby, H.W. 1988. Recovery Systems Design Guide, Beijing. Aviation industry press, 238 p. (in Chinese).

Summary

The extra wind field is an important factor affecting on parachute dropping. However, the existing test or numerical methods are difficult to investigate this effect. In order to investigate parachute airdropping in lateral wind environment, the explicit finite element method was used to calculate this complicated coupling process. The structure and flow field were discretized by Lagrangian meshes, and the coupling between canopy and air was realized by penalty function. To reduce the calculation amount, the fluid meshes were moved based on local coordinate transformation and the calculation domain was limited in a finite space. While, the boundary conditions were not changed with fluid meshes moving, and a stable lateral wind environment could be simulated. In this paper, the C9 parachute dropping process was calculated based on the above numerical method. This method can also be a reference for parachute design.

Keywords: fluid structure interaction, finite element method, extra wind field, graphical deformation.

Received December 03, 2015

Accepted February 06, 2017

H. Cheng (*), L. Wei (**), B. Zhou (***), Q. He (****), P. Xiao (*****)

(*) Civil Aviation Flight University of China, Guanghan 618307, China, E-mail: [email protected]

(**) Civil Aviation Flight University of China, Guanghan 618307, China, E-mail: [email protected]

(***) Civil Aviation Flight University of China, Guanghan 618307, China, E-mail: [email protected]

(****) Civil Aviation Flight University of China, Guanghan 618307, China, E-mail: [email protected]

(*****) Civil Aviation Flight University of China, Guanghan 618307, China, E-mail: [email protected] Table 1 Working conditions of calculation Condition Opening altitude, m 1830 Atmospheric pressure, Pa 8.12e4 Atmospheric density, kg/[m.sup.3] 1.023 Payload, N 980 Snatch velocity, m/s 20 Contrail declining angle, [degrees] 90