Numerical visualization of drop and opening process for parachute-payload system adopting fluid–solid coupling simulation

In order to study the fluid–solid coupling dynamic characteristics of parachute-payload system during drop process and analyze the unsteady aerodynamic characteristics under finite mass opening conditions, an adaptive moving fluid mesh method is developed on the basis of the existing arbitrary Lagrangian–Eulerian (ALE) fluid–solid coupling method. The calculation results of open force and drop velocity on the C-9 parachute demonstrate the effectiveness of this method. On this basis, the effect of canopies with three different permeability on parachute-payload system motion characteristic including opening property, steady descent property and stability is studied. Comparative analysis is conducted for structures and characteristics of vortex with different canopy materials, and interference mechanism of unsteady flow for parachute-payload system in unsteady oscillation is revealed. The results show that the adaptive moving fluid mesh method can effectively eliminate restrictions of existing simulation methods for parachute-payload system and significantly reduce calculation time. For the lightweight parachute, permeability has significant effect on kinetic characteristic of parachute-payload system. Canopy with large permeability has small opening load and structural stress in opening stage. After opening, there are mainly small vortexes distributed evenly behind the canopy with good stability. However, canopy with small permeability has obvious breath behavior and oscillation in opening stage. The main vortexes periodically shed off after opening. With the change of permeability from small to large, Parachute-payload system eventually presents three steady descent modes: conical descent, gliding descent and stable vertical descent.

simulation. By developing efficient and reliable fluid-solid coupling method of parachute, it is helpful to better understand its dynamic characteristics, reveal the corresponding flow mechanism, and provide reference and guidance for the design of various parachute systems. Therefore, it has important theoretical and practical significance.
The traditional parachute analysis method mainly establishes the multi-body model of the parachutepayload system through empirical formula and theory, and then obtains the motion process of the parachutepayload system based on dynamics (Cockrell 1987;Knacke 1992). Due to mature related theories and high computational efficiency, this kind of method has been widely studied at home and abroad: Doherr et al. (Doherr and Schilling 1992), Tang et al. (Qian-gang et al. 2007) established the 9-DOF dynamic model of parachute; Cheng et al. (Wenke et al. 1998) established 11-DOF inelastic dynamic model of general parachute-payload system; Ke et al. (Ke et al. 2009) developed multi-body dynamic model suitable for complicated parachute-payload system based on Kane equation. Although the method based on multi-body dynamics can be used to calculate the whole motion of the parachute-payload system, the structural deformation and the surrounding flow field details in the actual fluid-solid coupling process cannot be considered in the calculation.
With the development of computational fluid dynamics (CFD) and computational structural mechanics (CSD), a variety of CFD/CSD coupled computing methods have appeared: Tezduyar et al. (Tezduyar and Sathe 2007) proposed DSD/SST method based on finite element and verified the calculation by taking a variety of parachutes; Peskin 2006,2009) separately conducted calculation and analysis for opening process of small parachute-payload system under two-dimensional and three-dimensional conditions based on immersed boundary method; Cao et al. (Cao et al. 2012) conducted coupling calculation for circle flat parachute and conical parachute based on SIMPLE algorithm and analyze properties of the flow field at different attack angles; Fan and Xia (Fan and Xia 2014) proposed fluid-solid coupling algorithm based on finite element method and finite volume method with pretreatment and conduct verification and calculation for opening process of C-9 parachute. Although this kind of method can obtain more detailed flow field details, it requires a large amount of calculation. At present, it is mainly used to analyze the flow field characteristics in the specific opening stage, and it is difficult to apply it to the calculation of the whole process of the Parachute-Payload system.
In recent years, with the popularization of ALE(Arbitrary Lagrangian-Eulerian) methods in the field of fluid-structure coupling and the development of relevant computing software and hardware, ALE-based fluid-structure coupling modeling and computing technology has become a research hotspot in the field of aerospace return and recovery simulation at home and abroad . In 2003, Tayler (Taylor 2003) used ALE method to calculate and analyze the fluid-solid coupling of cross parachute for the first time. In 2005 * 2009, Lingard (Lingard and Darley 1607), Tutt et al. (Tutt 2006;Tutt et al. 2010Tutt et al. ,2011 further developed the method and applied it to various kinds of parachutes and opening conditions , and good results have been obtained. Gao et al. (Gao et al. 2013 Ma et al. (XiaoDong et al. 2015) successfully applied ALE method in opening process computation of rotating parachute. However, the above research mainly focuses on calculation of inflation process in infinite mass under steady flow condition, and focus on the inflation performance and steady descent characteristics, while pay less attention to the kinetic characteristics of the payload and the overall stability of the parachute-payload system during the practical drop and opening (Benson and Hallquist 1986). At present, calculation for parachute-payload system opening in finite mass is still rare at home and abroad, and the fluid-solid coupling method under the condition of relevant finite mass is still in the initial stage of development. Coquet (Coquet et al. 2011) proposed a simulation method of fluid-solid coupling of the parachute in finite mass situation in 2011 and conducted verification and calculation, Cheng et al. (Cheng Han et al. 2014) conducted preliminary calculation and analysis by adopting similar method with C-9 parachute as computing object. However, the method needs to establish flow field grid for the whole motion area of the parachute-payload system, and the calculated amount is directly restricted by the motion range. Therefore, it is still very difficult to apply it to the analysis of widerange and long-term drop and opening process. In addition, it is difficult to visualize and accurately simulate the flow field vortex after opening the parachute. At present, experimental methods are used to study vortices in flow field (Lau et al. 2019;Ma et al. 2020).
In order to enable fluid-solid couple computation of parachute-payload system in finite mass situation more effectively, the main works of this thesis is as follows: Firstly, a numerical simulation method using adaptive moving fluid mesh is developed based on the research of existing ALE method. Then, the drop and opening process of the C-9 parachute is computed by using the method established; the reliability of computed results is also analyzed. Finally, on the basis of the C-9 parachute example, a detailed analysis of the drop and opening process of a certain kind parachute is carried out. Influence of permeability of canopy fabric on kinetic characteristic is deeply studied. The mechanism of unsteady flow disturbance during parachute oscillation is revealed, and some significant conclusions are obtained.
2 Numerical method 2.1 Model formation and adaptive control of the mesh Differing from the condition that parachute is fixed to one point and inflow speed is known in wind tunnel test(opening in infinite mass), practical drop process (opening in finite mass) of parachute-payload system usually involves a large spatial range of motion. Moreover, the velocity and attitude of the parachutepayload system change in real time with the fluid-solid coupling of the opening process, which cannot be predicted in advance. In the literature Cheng Han et al. 2014), computation for opening process in a finite mass situation is achieved through the way of expansion of computed field that covers the whole motion area. Therefore, calculated amount is directly restricted by the motion range of parachutepayload system, and usually only small range of motion can be calculated.
Since existing methods are difficult to meet the requirements of the simulation of wide-range and longterm drop process, this paper further develops an adaptive motion method of flow field calculation domain. Schematic of adaptive moving fluid mesh is shown in Fig. 1.
As shown in Fig. 1, due to adoption of moving fluid mesh, there is no need to generate mesh in the whole motion area. Instead, it is only necessary to establish a local flow field calculation domain of appropriate size around the parachute-payload system, and coverage for the whole opening motion area is achieved by controlling dynamic tracking of flow field grid for drop motion, which can significantly reduce grid amount required in drop motion process and eliminate limitations of existing methods on the time and range of motion.
Adaptive moving fluid mesh is realized based on the inherent free transformation characteristic in ALE method. The specific grid motion strategy is to take the body-axes system fixed to the payload mass center as the ALE dynamic reference system of mesh translation and rotation. The translation takes the position of Fig. 1 Schematic of adaptive moving fluid mesh mass center as the reference, while the rotation takes the velocity direction of mass center as the reference, that is, the flow domain always follows the motion of the payload and the axis direction is consistent with the motion direction of the payload. If the velocity vector of the payload mass center at some time is v d , the adaptive velocity vector w i of flow field grid node i at this moment can be expressed as: where r i is the radius vector of flow field grid node i in ALE dynamic reference system. Based on the above methods, computation model includes payload, parachute and self-adaptive fluid. For the computational domain of parachute-payload system consisting of payload and parachute, the modeling method is consistent with the infinite mass calculation in references 17, 19 and 23, the deformation and movement of the structure are tracked by the body-fitted grid described by Lagrange, and the appropriate mesh type is selected according to the characteristics of each structure component. Components of parachute-payload system are shown in Fig. 2; element formulation and mesh type of each component in parachute-payload system are shown in Table 1.
The computational domain of flow field is divided into hexahedral elements as spatial background grid embedded with structured mesh of parachute system. In order to better capture the flow field details around the area of canopy and tailing vortex, a local mesh encryption is adopted for background grid. In order to minimize the influence of boundary truncation on calculation accuracy, the mesh size and density of the flow field are determined by mesh independence verification on the basis of a trial calculation example (see Sect. 3). Under the condition of both calculation accuracy and calculation amount, the cross-section diameter of the flow field is 5D0, and the height is 7D0, where D0 is the nominal diameter of the parachute. The ratio of the size of the flow field element near the canopy to the size of the canopy surface element is controlled at about 2:1.
In the setting of boundary conditions, the adaptive flow field boundary follows the movement of the payload constantly changes in this paper. Therefore, different from the practice of setting the flow field boundary as the constant velocity inlet and wall or far field conditions in the calculation of infinite mass opening, this paper sets the outer surface of the entire flow field as the reflection-free boundary, and the boundary unsteady flux is determined by the velocity and direction of the background grid following the parachute-payload system.
Embedded mesh model for FSI simulation is as shown in Fig. 3.

Numerical model of the flow field
The governing equation is N-S equation based on the description of arbitrary moving reference system and its differential form is as follows Souli and Benson 2013): where q is fluid density; v is speed vector of fluid material; w is speed vector of flow field grid under (1) is that f is volume force; e is specific energy; r is Cauchy stress tensor, the equation is: where I is 3 9 3 unit matrix; p is fluid pressure. The article makes assumption of air as ideal gas, so p and q satisfy the ideal gas law; l is dynamic viscosity coefficient and expressed as the sum of laminar viscosity and turbulent viscosity. Prandtl Mixing-Length model is adopted for calculation of turbulent viscosity. This model has a certain ability to simulate turbulence and does not significantly increase the amount of fluidsolid coupling calculation (D. C. 1998). Eq.
(2) shows that ALE method is essentially superset of Lagrange method and Euler method: if w = v is always satisfied, ALE motion coordinates system coincides with material coordinates system, the above N-S equation is degraded into Lagrangian description; if grid motion speed w is always 0, ALE motion coordinates system coincides with Eulerian space coordinates system, the above N-S equation is degraded into Eulerian description. The solution of Eq. (2) is based on nonlinear finite element program LS-DYNA, which conducts spatial dispersion with finite element method and use explicit time marching with central difference scheme (Souli and Benson 2013). For both computational efficiency and accuracy, van Leer multidimensional upwind scheme (MUSCL) with second-order accuracy is adopted for the flow flux calculation  Numerical visualization of drop and opening process for parachute-payload system at the interface of flow field units. The research shows that (Gao et al. 2014) compared with first-order donor cell scheme, the scheme can obtain enough calculation accuracy with less mesh quantity.

Numerical model of the structure
For the parachute canopy, constitutive model of fabric material based on Belytschko-Tsay membrane element is used as its governing equation. Pruett et al. (2009) tested the fabric constitutive model and proved that the model can accurately predict the geometric nonlinear behavior of folds and obtain reliable mechanical properties of materials such as stress and strain during large deformation of flexible fabric. The relationship between the stress increment Dr ij and strain increment De ij under the local coordinate system of the node 3 and node 4 membrane element is as follows (Hallquist 2006): where E ij , tt ij and t ij are separately Young modulus, shear modulus and Poisson's ratio. For suspension line, skirt tape, vent tape and riser, the linear elastic discrete cable unit is used for modeling because it only produces axial tension under the tension state and the axial tension T of the cable unit is expressed as: where E is the elastic modulus; A is the cross sectional area; L 0 and L t are separately initial length and the current time length of the rope, respectively. For payload, its own aerodynamic force is neglected in the present study, the shape of the payload is simplified as a cuboid with the same size as the actual payload. The kinetic equation of 6-DOF rigid body (Benson and Hallquist 1986) is used to describe its motion under the action of gravity and riser tension. Vector form of dynamic equation under body axis system of payload is as follows: where v b and x b are separately velocity vector and angular velocity vector of payload; G b and T b are separately gravity vector of the payload and tension velocity vector of the riser in the body-coordinate system; m and J are separately mass and inertia tensor of the payload; r T is radius vector where the riser hanging point is located. Eq. (6) are ordinary differential equations. Explicit time-marching method combined with the initial conditions can be used to solve the motion of the payload at each moment.

Algorithm of interface coupling and canopy permeability
Ergun theory is adopted to describe the permeability of canopy fabrics. Permeability under pressure difference can be calculated through the following Ergun equation (Ergun 1952): where DP is the pressure difference between two sides of fabric; d n is the fabric thickness; V n is the average normal velocity of air passing through the fabric surface, that is, permeability; a and b are, respectively, the viscosity coefficient representing the viscous dissipation and the inertia coefficient representing the kinetic energy loss (hereinafter referred to as the Ergun coefficient), which are generally obtained through the permeability test of the fabric. A lot of research (Tutt 2006;Aquelet et al. 2006;Wang et al. 2006) shows that, although Eq. (7) is an empirical formula derived from statistics, the above equation can accurately simulate the impact of fabric permeability on parachute performance on a macro-scale if a reliable Ergun coefficient can be obtained.
Interface coupling force between canopy element and flow field is calculated through penalty function method. In every time step Dt, the contact algorithm of coupling of flow field node and structure node is as shown in Fig. 4 (Souli and Benson 2013).
The coupling force between the canopy element and the flow field element is calculated by the penalty function in the contact algorithm. Schematic of fluid-structure interface penalty coupling is shown in Fig. 4 (Souli and Benson 2013).
As shown in Fig. 4, contact force F f on flow field node and contact force F s on structure node are equal and opposite: Fs ¼ kd À C d , k is the contact stiffness determined by the material properties, C is the damping coefficient, d is the contact displacement which is calculated by the relative velocity between the structure node and the flow field node, that is d where v s and v f are separately velocity vectors of structure node and flow field node.

Method validation
In order to verify the effectiveness of the methods in the paper, C-9 parachute dropping tests conducted by Calvin Lee et al. (Lee 1984;Lee et al. 1997)  In this paper, the starting point is the deploying moment of C-9 parachute in the experiment. According to the data given in the original 36 , the corresponding time is t = 1.6 s and the deployment velocity is 13.4 m/s. In order to get closer to the instantaneous natural form of the canopy, the method proposed in the literature ) was used to establish the initial folding shape of the canopy. Other specific model dimensions and material parameters are taken from the literature (Lee 1984;Lee et al. 1997) and will not be given here.
Compared with conventional CFD calculation, due to the additional calculation of structural deformation and coupling force of interface, the computational load of FSI is obviously greater than that of CFD under the same flow field grid. Therefore, it is necessary to choose reasonable mesh size and density to achieve the balance between calculation accuracy and calculation amount. In order to simultaneously investigate the influence of the relative size of the computational domain and grid quantity in the moving fluid mesh method presented in this paper on the computational accuracy, three sets of flow field grids with different ranges of size and quantities of elements are established for verification, as shown in Fig. 5, D0 is the nominal diameter of the of parachute canopy, the number of flow field elements corresponding to mesh model A, B and C is 63972, 98,700 and 157,920, respectively. Figure 6 shows the comparison between the calculated results (Mesh C) and the shape change of the C-9 parachute in practical drop process (Lee 1984;Potvin and McQuilling 2011). It can be seen that inflation process of the simulated parachute is in good agreement with practical situation, which shows that FSI method adopted in this paper can simulate the details of canopy inflation shape in practical drop tests. Comparison of the open force results of C-9 canopy validation case (Lee et al. 1997) is shown in Fig. 7. It can be seen from Fig. 7 that size of flow field has obvious effect on calculation results, the variation process and trend of open force calculated by the mesh configuration of each flow field are basically consistent with the experimental data. With increases of size of flow field and amount of mesh, peak open force and its occurrence time are gradually close to test results. When mesh C is adopted, the peak value of open force and the time of emergence are in good agreement with the test results. Meanwhile, calculation results of open force increases faster than test data in initial inflation stage before 2.0 s, which may be caused by the difference of initial canopy shape used in this paper and the actual shape at the moment of deployment; compared to the test data, the calculation result of open force fluctuates obviously in steady descent stage after 3.0 s, which may be due to the fact that the suspension line and the damping effect  Comparison of simulated C-9 canopy inflation shape with drop test (Potvin and McQuilling 2011) caused by the payload motion are ignored in the current calculation. Considering the possible error and uncertainty in the actual test, the deviation of the calculation result of the open force using mesh C has been within the acceptable range. Figure 8 shows comparison of payload descent velocity results. As can be seen from the figure, the calculated descent velocity is smaller than test result before 2.5 s, but maximum difference is not greater than 1.5 m/s. In combination with Fig. 7, we can see that fast descent velocity of payload is caused by the large open force in initial stage. Variation trend of descent velocity and final steady descent velocity are very close to test result, which indicates that calculation results has already accurately reflected dynamic properties of real drop and opening process under local flow field size (5D0 of radial direction, 7D0 of longitudinal direction).
In terms of computation, when mesh C is adopted, the calculation time of this example on a common PC (CPU main frequency 3.5 GHz, quad-core SMP parallel) is 40276 s (11 h 11 min 16 s). Compared with the calculation time of 400 h required by the global flow field method on the same example in reference (Cheng Han et al. 2014), the adaptive moving fluid mesh method adopted in this paper has obvious advantages.

Calculation configuration
Parachute constructed profile is shown in Figure 9. Canopy consists of eight gores, tile diameter of the canopy D0 = 1000mm, vent diameter Dv = 0.1D0; the length of suspension line Ls = D0, eight suspension lines are seized and connected to riser with payload, the length of riser Lr = 0.08Ls.Skirt tape and vent tape are separately arranged on the bottom and vent of the canopy, four radial reinforcing belts run through the vent and are connected with the suspension line at both ends.
Canopy permeability (including permeability of structural geometry and fabric material) is the main factor affecting parachute aerodynamic performance (Knacke 1992). Similar parachutes are adopted to conduct study on inflation process under steady inflow condition in the literature 20, and the effect of geometric permeability is also analyzed. Because practical light payloads are mostly drop in batches, it is more likely to use solid parachute canopy with simple manufacturing and low cost (such as US army C-9, G-11, G-12 and other types of parachutes). In order to investigate the effect of permeability on various performance in drop process, three typical canopy fabrics with different permeability are selected for simulation analysis. Material properties of canopy fabric used are shown in Table 2. 544 kam silk is the fabric with micro-permeability, 59,225 kam silk is the fabric with medium permeability and 509 kam silk is the fabric with large permeability (for sequent calculation results, mat544, mat59225 and mat509 used separately refer to three fabric materials of canopies). In the table, a and b are Ergun coefficient mentioned above, which are obtained by least square fitting of multiple sets of test data based on the method in reference , The permeability test results within the range of 0 * 800 Pa pressure difference are shown in Fig. 10, and the range of error lines in the figure is determined by the standard deviations of the three repeated tests under each pressure difference. Fig. 7 Open force results of C-9 canopy validation case

Computational condition
The above example analysis shows that initial shape of the parachute at deployment time has large uncertainty and it is difficult to model accurately. In order to facilitate the subsequent comparative verification of the corresponding drop test and reduce the uncertain interference factors under the actual drop condition as far as possible, the gravity drop method similar to the test in reference (Desabrais et al. 2007) is adopted in this paper, and atmospheric wind interference is not considered. Initial conditions and main physical conditions of simulation are as shown in Table 3.

Analysis of computational results
5.1 Analysis of opening performance Figure 11 shows canopy opening diameter time history results of three different canopy materials. Drop process can be divided into three stages. Taking canopy with mat544 for example, the shape change process of the canopy corresponding to its three stages is given: i) In the first stage of 0 * 0.42 s, the airflow enters from the bottom and reaches the top of the canopy. In this process, the diameter of the bottom opening increases relatively slowly ii) In the second stage of 0.42-0.63 s, as the top of the canopy is saturated, the airflow that continues to enter rapidly expands the bottom edge of the canopy. In this process, the bottom opening diameter increases rapidly. At the same time, it can be seen that the bottom edge diameter of the canopy does not maintain a stable state after it reaches the stable opening diameter for the first time, but continues to increase under the action of inertial force and aerodynamic coupling, resulting in obvious top collapse and excessive inflation iii) In the third stage after 0.63 s, the canopy begins to contract and breathe. During this process, the diameter of the bottom opening fluctuates back and forth with the breath and gradually tends to be stable.
This process is consistent with the analysis conclusion of the opening process of the solid parachute canopy given in reference (Potvin and McQuilling 2011) and accords with the actual law of the opening movement of this kind of parachutes. Figure 12 shows open force time history results in the drop process. Combining the both Fig. 11 and Fig. 12, it can be found that fluctuation frequency and intensity of opening load are approximately synchronous with fluctuation of the canopy shape, which indicates that fluctuation of open force is essentially caused by the change of canopy in inflation process; meanwhile, the frequency and intensity of breath fluctuation of 509 kam silk with large permeability are obviously lower than the two others, and relevant peaking and vibration values of open force are significantly reduced. This indicates that the fabric with high permeability can achieve the effect similar to the slit or hole in the structure of the canopy, which is beneficial to slow down the inflation speed, reduce the opening load and improve the opening stability. Table 4 shows comparison of canopy opening performance parameters. t f is the time of full inflation which refers to the time that canopy is fully inflated as it achieves steady descent for the first time 1 , t mo is the time of maximum opening of the canopy, t p is the time of peak opening load, F p is the maximum open force.
It can be seen from Table 4 that peak open force occurs before achieving maximum opening of the canopy. The less permeability, the faster inflation speed and the earlier maximum dynamic opening load  Figure 13 shows canopy surface effective stress contour at peak opening load instant. It can be seen that large dynamic opening load of the canopy with low permeability is corresponding to large stress of the canopy; stress distribution and peak stress value of mat544 with micro-permeability are obviously greater than those of mat59225 with medium permeability and mat509 with large permeability.  Meanwhile, it can be seen that the canopy stress is mainly distributed around the vent and the bottom of the canopy, which indicates that slotting or holing will inevitably cause stress concentration. Therefore, vent loop, bottom band and other reinforcing components are needed to avoid tearing. Figure 14 shows payload descent velocity time history results. In the initial descent stage, the parachute is approximately in free fall that the descent speed is increasing rapidly. This is because the external force of payload suffered is mainly contributed by its own gravity in this stage, parachute deployment is not fully and the air drag of parachute is small. Since the moment of full inflation, the significant difference in descent speed emerges: parachute-payload system using mat509 with large permeability has the fastest descent speed, the variation of descent velocity which stabilizes at 7.9m/s after 3.3 seconds is relatively smooth; descent speed of mat59225 with medium permeability which stabilizes at 7.6m/s after 4 seconds is between mat509 and mat544; Parachute-payload system using mat544 with micro-permeability has the slowest descent speed with obvious fluctuations between 7.15m/s and 7.45m/s after 5 seconds. According to calculation equation of drag coefficient during terminal descent :C D = 2W T /(qv e 2 S 0 ) (W T is the total weight of parachute-payload system, v e is the terminal descent velocity and S 0 is nominal area of the canopy), drag coefficient variation during terminal descent is shown in Fig. 15.

Analysis of motion characteristics of parachute-payload system
Drag coefficient of mat544 fluctuates at 0.8, which coincides with the range of drag coefficient 0.75 * 0.80 given in empirical data (Knacke 1992). Increasing of permeability will cause decreasing of drag coefficient, drag coefficient of mat509 is 0.68 which decreases by 15% compared with mat544. Figure 16 and Fig. 17 show payload flight trajectory and path angle time history results. However, wind disturbance and other uncertainty factors are not considered in calculation, the motion of drop still exists uncertainty after parachute opening, which is because circle flat parachute is a typical unsteady parachute. Numerical visualization of drop and opening process for parachute-payload system Unstable attitude of parachute-payload system will cause continuous variation on motion speed and direction.
To further indicate the effect of canopy permeability on stability of parachute-payload system. As shown in Fig. 18, attitude angle time history results of parachute and payload is presented, and Ux and Ux represent the rotation about the X-axis and Y-axis in a fixed spatial coordinate system, respectively. Compare Fig. 18a with Fig. 18b, it can be found that, in essence, the fluctuation of the attitude angle of the payload is the result of the superposition of the large swing of the parachute with long period and the relative swing of the payload with small amplitude in short period. Meanwhile, it can be seen that attitude changes of mat59225 with medium permeability and mat509 with large permeability are both smooth and the fluctuation range of attitude angle is less than ± 15 AE during terminal descent; there is obvious unstable oscillation in the whole descent process of mat544 with micro-permeability and the fluctuation range of attitude angle is close to ± 40.
It can be seen from the changes of flight path angle in Fig. 17, parachute-payload system using mat509 finally descends at small attitude angle and about 90 AE flight path angle; parachute-payload system using mat59225 finally glides down at small flight path angle; parachute-payload system using mat544 finally descends with coning motion unstably. The above three motion modes are consistent with practical airdrop (Knacke 1992).

Analysis of flow field characteristics
In order to further analyze the cause of the instability of the parachute-payload system with low permeability, Fig. 19 shows velocity contour map during descent. Figure 19a shows vortex development during parachute inflation at initial drop stage (0.4 * 0.8 s). It can been seen that vortex structures behind the canopy is symmetric and no obvious oscillation of the parachute occurs at the stage. The low pressure area caused by the formation of the symmetric vortex system causes the canopy to expand outward under the pressure difference. When the canopy expands to a certain degree, the vortex system begins to shed off the canopy. In this process, the canopy shrinks under the action of inflow pressure. It is the periodic generation and separation of the vortex system that causes the breath phenomenon of the canopy. Fig. 19b shows vortex development during parachute oscillation (2.0 * 2.2 s). It can be seen that a complex asymmetric separation vortex system dominates the rear of the canopy. Parachute oscillates due to lateral force caused by separated vortex. Main separation vortexes on both sides fall off alternately and move downstream relative to the canopy. Figure 20 shows comparison of canopy surrounding flow field at t = 5.0 s. It can be seen from Fig. 20a that, when the air permeability of the canopy is large, the vortex structure generated by the airflow in the rear trailing vortex area after being disturbed by the canopy is more broken, and the vortex distribution is more uniform, that is because the air passing through the canopy surface accelerates the separation of the rear airflow and reduces the backflow area, thus inhibiting the formation of large-scale separation vortexes. Compared with mat544 with micro-permeability, there are no obvious main separation vortexes behind the Numerical visualization of drop and opening process for parachute-payload system canopy 509 with large permeability, and the main vortexes are small ones with uniform distribution. Therefore, lateral force is small and the whole parachute-payload system is steady.
The vorticity given in Fig. 20b is the y-direction (perpendicular to the paper surface inward) component, so a positive vorticity represents a clockwise rotation and a negative vorticity means counterclockwise rotation. It can be seen that the separation vortexes on both sides of the central axis of each parachute are opposite, and there are large vortexes at the bottom edge and the edge of the vent, which indicates that the flow moves downstream in the form of separation vortexes after the separation occurs at the bottom edge and the edge of vent, and the viscous dissipation vortexes gradually decrease with the movement. At the same time, it can be seen that the vorticity distribution on both sides of the rear of the large permeability Canopy 509 with small vortexes is symmetric and continuous, while the vorticity asymmetry and discontinuity in the rear of the small permeability Canopy 544 with large separation vortexes are obvious.

Conclusion
1. An adaptive moving fluid mesh method is developed on the existing fluid-solid coupling method in this paper. Validation of the method is carried out based on the C-9 parachute. The results show that the method can effectively eliminate restrictions of existing simulation methods for parachute-payload  system. The results of the calculation are in good agreement with the test data which demonstrates the reliability of this method. 2. For the three parachute-payload systems with different air permeability of canopy, calculation results under the finite mass of drop and opening process shows that canopy fabric with large permeability is beneficial to slow down inflation speed, reduce opening load and structural stress of the canopy structure and improve opening stability, where the effect is similar to slotting and vent of the canopy structure. 3. In allusion to the solid parachute canopy studied in this paper, parachute-payload system has obvious fluctuation of attitude angle during terminal descent. The larger permeability, the slower attitude angle fluctuation. Parachute-payload systems of Canopy 544 with micro-permeability, Canopy 59,225 with medium permeability and Canopy 509 with large permeability present three steady descent modes after fully opening: conical descent, gliding descent and stable vertical descent. 4. The analysis of flow field characteristics of each parachute-payload system shows that the canopy breath behavior is caused by periodical generation and separation of vortex; the low pressure area generated during the formation of the symmetric vortex system will make the canopy expand outward. When the vortexes shed off, the canopy will contract under the action of inflow dynamic pressure. After the canopy is fully opened, asymmetrically large-scale vortexes dominate the rear of the low-permeability canopy, and the main separation vortexes on both sides alternately shed off, causing the oscillation of the parachute-payload system, However, the small-scale vortexes evenly distributed behind the highpermeability canopy are dominant, so the parachute suffers less lateral force and the whole parachutepayload system is relatively stable.