Fluid-structure interaction simulation for performance prediction and design optimization of parafoils

Abstract

Parachute design is challenging to achieve innovative progress if the dominant role of testing continues, as it will be an increasingly expensive and time-consuming work. The aim of this study is to establish a reliable and efficient design tool using existing advanced numerical modeling methods. This paper presents a numerical method based on two-way coupled fluid-structure interaction (FSI) strategies for predicting aerodynamic and flight performance for parafoil design optimization. The nonlinear finite element method was used for the canopy fabric model and flow field, and the fluid dynamics were solved by Reynolds-averaged Navier-Stokes with the Spalart-Allmaras turbulence model. The FSI simulations are performed to assess the aerodynamic performance and structural deformations of full-scale parafoil canopies. The equilibrium shape of the parafoil canopy under steady gliding states and the relevant flow field were analyzed to enhance confidence and understanding in the performance prediction of new parachutes. Three-dimensional FSI simulation results of parafoils show that the inflation caused flexible bulges of canopy cells, and the maximum lift coefficient increased more than 16% with a higher stall angle of attack than that of the rigid body model. A parafoil with a smaller leading edge inlet or a scaling down area can improve the aerodynamic performance, mainly manifested in a higher lift-to-drag ratio and better anti-stall performance. Finally, the prediction results of parafoil glide performance were verified by flight test data, and the prediction accuracy of the flexible model is more than 10% higher than that of the rigid model. This work makes the simulation tools a step closer to practical application.

Keywords:

• Fluid–structure interaction
• parafoil design
• aerodynamic characteristics
• performance prediction
• flexible deformation
• inflation

1. Introduction

Parafoil systems have been widely used in military delivery and aerospace recovery areas for their benefits of gliding performance and maneuverability, such as the Crew-Return Vehicles X-38 of NASA (Bennett & Fox, 2003), ALEX of the Institute of Flight Systems of the German Aerospace Center (DLR) (Jann, 2006), and fairing recovery program of SpaceX (Seedhouse, 2022). Failure of aerospace missions will lead to enormous economic and human losses. However, the mission success rate is still far from the requirements, which brings higher challenges to the performance prediction and design industry of parafoils (Bonetti et al., 2019; Cacan & Costello, 2018). Therefore, it is necessary to study the accurate prediction of flight performance, such as aerodynamic lift, drag and glide ratio, to improve the success rate of missions.

The aerodynamic performance prediction is critical to parafoil design optimization. Over the last decades, the development of parachute design industry has mainly relied on empirical or semi-empirical methods for the prediction of aerodynamics, based on data obtained from wind tunnel experiments and airdrop tests (Wang & He, 2018; Yang et al., 2017). Parachute design is challenging to achieve innovative progress if the dominant role of testing continues, as it will be an increasingly expensive and time-consuming work (Jiang et al., 2018). The aim of the present work is to establish a reliable and efficient design tool using existing advanced numerical modeling methods.

In terms of predicting the aerodynamic characteristics of parafoil systems in steady flight, many studies have used computational fluid dynamics (CFD) software to calculate the flow field around the parafoil, treating the canopy as a rigid body. The airflow around two-dimensional (2-D) airfoils (Balaji et al., 2005; Mittal et al., 2001; Mohammadi & Johari, 2010) and three-dimensional (3-D) parafoils (Cao & Zhu, 2013; Eslambolchi & Johari, 2013) has been analyzed using CFD software. In these studies, the Reynolds-Averaged Navier-Stokes (RANS) equations were used to solve the flow field around rigid parafoils and airfoil cross-sections. The above studies revealed that the sharp leading edge leads to separation of the airflow flowing through the upper and lower surfaces, and a parafoil with open leading edge produces less lift and more drag than that with smooth closed leading edge. To bridge the gap between flight dynamics simulation of parafoil systems and computational fluid dynamics, some scholars have combined aerodynamic data from CFD calculations into the dynamics model to improve the fidelity of model and prediction accuracy of flight performance for parafoil systems (Ghoreyshi et al., 2016; Wu et al., 2019).

Due to the highly flexible and reconfigurable canopy structure, parachute dynamics involve complex nonlinear fluid-structure interactions (FSI) during the working process. Some studies (Fairuz et al., 2013; Takizawa et al., 2016; Young et al., 2009) have demonstrated that the accurate representation of wing shapes is essential for reliable aerodynamics analysis. The above research based on CFD simulations ignored flexible deformations of the canopy, which limits the precise prediction of aerodynamic coefficients. However, accurate aerodynamic coefficients are essential to predict the dynamic flight performance of parafoil vehicles. To solve the above problems, with the development of computer technology, one of the solutions is to establish a high-fidelity dynamic model based on FSI simulations (Ghalandari et al., 2019; Liu et al., 2022).

In view of the complexity of initial geometry, previous studies on parachute FSI have been more about round parachutes (Qu et al., 2019; Stein et al., 2000; Takizawa & Tezduyar, 2012; Xue & Wen, 2021), while parafoil FSI has mainly focused on the inflation or equilibrium configurations (Altmann, 2015; Kalro & Tezduyar, 2000; Nie et al., 2018; Ortega et al., 2017; Zhang et al., 2021). The inflated shapes of parafoils with closed leading edge (Kalro & Tezduyar, 2000), and 3-D full-scale parafoils (Altmann, 2015; Nie et al., 2018; Ortega et al., 2017; Zhang et al., 2021) were simulated by combining fluid solver and structural dynamics solver. Recent studies have shown that the incompressible computational fluid dynamics techniques of LS-DYNA is a practical tool for precisely modeling FSI problems (Charles, 2017; Garrec et al., 2019; Peralta & Johari, 2015), which can transfer data between flow fluid and structural domains automatically to promote the solution efficiency. Although these studies improve the understanding of aerodynamic characteristics of parafoils, there is a lack of research on the effect of flexible canopy on flight performance predictions for the development and design optimization of parafoils.

The goal of this paper is to establish an accurate model to predict the aerodynamic and flight performance of highly nonlinear flexible parafoils based on advanced FSI technology, reducing design optimization costs. The LS-DYNA solver is used to solve highly flexible fabric structures and combined with the incompressible computational fluid dynamics solver to achieve a two-way FSI coupling. The fluid loads and structural displacement data are transferred through the FSI interface to solve the highly nonlinear problem. The aerodynamic characteristics and structural deformation results of different fully flexible 3-D parafoil models with rectangular planform are obtained under steady gliding conditions. Considering the influence of spanwise and chordwise reductions caused by flexible deformation of the parafoil, the visualizations of the 3-D flow field are analyzed. To the best of our knowledge, this is the first numerical prediction platform developed for the actual airdrop application of the parafoil system, which is used to solve the complex aerodynamics and structural dynamics for 3-D full-scale flexible parafoils. The main contributions of this paper are as follows:

1. We establish a general FSI coupling framework for the parafoil, which combines inflation deformations of the canopy with aerodynamic calculations to predict the parafoil system's aerodynamic characteristics more accurately.

2. Based on the established models, the influences of crucial geometric design parameters on aerodynamic performance are simulated and analyzed for parafoils.

3. Aerodynamic coefficients obtained by FSI simulations are integrated into flight dynamics equations to predict the flight characteristics of the parafoil system. The accuracy of the numerical prediction results is validated by airdrop experiments.

This paper is organized as follows. Section 2 describes the main issues that need to be addressed to predict the gliding performance of parafoil systems. Section 3 provides the governing equations used for fluid-structure interaction simulations. Section 4 presents different parafoil models constructed in the design stage and finite element models used for the following simulations. The FSI simulation results and discussions of FSI coupling are presented in Section 5. Predictions of parafoil glide performance are performed based on numerical simulations and verified by comparing with flight test data in Section 6. Finally, we conclude with some key findings in Section 7.

2. Problem statement

The schematic of the parafoil system in a steady gliding state is shown in Figure 1. The aerodynamic lift force $L_c$ and drag force $D_c$ of the parafoil system are obtained by decomposing the total aerodynamic force along the incoming flow direction $V_{\mathrm{\infty }}$ and its vertical direction (Knacke, 1987). According to the force balance of the system, we can obtain

Figure 1. Schematic diagram of parafoil system in steady gliding flight.

(1) $\begin{array}{rl}& L_c\sin \gamma =D_c\cos \gamma \end{array}$(1) (2) $\begin{array}{rl}& Mg=L_c\cos \gamma +D_c\sin \gamma \end{array}$(2) where $M=m_c+m_p$ is the total mass of the system (including ropes), g is the gravity acceleration, and γ is the glide angle. Equation (1(1) $\begin{array}{rl}& L_c\sin \gamma =D_c\cos \gamma \end{array}$(1) ) can be transformed into: (3) $\begin{array}{rl}& L_c/D_c=C_L/C_D=1/\tan \gamma \end{array}$(3) (4) $\begin{array}{rl}& C_L=L_c/0.5\rho V_{\mathrm{\infty }}^2{S}_{ref}\end{array}$(4) (5) $\begin{array}{rl}& C_D=D_c/0.5\rho V_{\mathrm{\infty }}^2{S}_{ref}\end{array}$(5) where Equation (3(3) $\begin{array}{rl}& L_c/D_c=C_L/C_D=1/\tan \gamma \end{array}$(3) ) is the formula for aerodynamic efficiency in gliding flight. The smaller the glide angle, the larger the lift-to-drag ratio and the better the gliding performance. $C_L$ and $C_D$ are the lift coefficient and drag coefficient, respectively. ρ is the air density, and $S_{ref}$ is the reference area of the parafoil canopy.

Substituting Equation (3(3) $\begin{array}{rl}& L_c/D_c=C_L/C_D=1/\tan \gamma \end{array}$(3) ) into Equation (2(2) $\begin{array}{rl}& Mg=L_c\cos \gamma +D_c\sin \gamma \end{array}$(2) ) yields (6) $W=Mg=0.5\rho V_{\mathrm{\infty }}^2{S}_{ref}{(C_L^2+C_D^2)}^{1/2}$(6) Therefore, the velocity of parafoil system in steady gliding state can be written as: (7) $V_{\mathrm{\infty }}=\frac{2W}{\rho S_{ref}(C_L^2+C_D^2)}$(7) The horizontal velocity $V_h$ and vertical velocity $V_d$ are obtained as: (8) $\begin{array}{rl}{V}_h& =V_{\mathrm{\infty }}\cos \gamma \\ V_d& =V_{\mathrm{\infty }}\sin \gamma \end{array}$(8) From Equations (7(7) $V_{\mathrm{\infty }}=\frac{2W}{\rho S_{ref}(C_L^2+C_D^2)}$(7) ) and (8(8) $\begin{array}{rl}{V}_h& =V_{\mathrm{\infty }}\cos \gamma \\ V_d& =V_{\mathrm{\infty }}\sin \gamma \end{array}$(8) ), it can be seen that the gliding velocity of parafoil system is mainly related to system weight W, air density ρ, and aerodynamic coefficients $C_L$ and $C_D$. Accurate aerodynamic coefficients are required for flight performance prediction in the design stage of parafoil systems. However, considering the highly nonlinear flexibility of parafoil canopy, previous methods based on lifting-line theory and CFD methods for calculating aerodynamic parameters of rigid vehicles are not applicable (Zhu et al., 2022). And there is no universal aerodynamic model suitable for various configurations to provide accurate aerodynamics for parafoil flight dynamics. To build a high-fidelity prediction dynamic model, FSI simulations in Section 5 were conducted to compute the equilibrium shape and aerodynamic forces simultaneously.

3. Governing equations of fluid-structure interaction

Consider a fluid domain ${\mathrm{\Omega }}_f\left(t\right)$ is occupied by incompressible, viscous air with a boundary ${\mathrm{\Gamma }}_f\left(t\right)$, and ${\mathrm{\Omega }}_s\left(t\right)$ is the structure domain with a boundary ${\mathrm{\Gamma }}_s\left(t\right)$. For the convenience of simplification, the form of time t is no longer shown below.

3.1. Structural mechanics

Due to the large structural displacements in inflating and trailing edge deflecting, the canopy deformation is regarded as a geometrical nonlinear problem. The equations of motion are written as (9) ${\rho }_s\frac{\mathrm{\partial }{\mathit{v}}_s}{\mathrm{\partial }t}={\rho }_s{\mathit{f}}_s+\mathrm{\nabla }\cdot {\mathit{\sigma }}_s\mathrm{i}\mathrm{n}{\mathrm{\Omega }}_s$(9) where ${\rho }_s$ is the material density, ${\mathit{v}}_s$ is the velocity of structure, ${\mathit{\sigma }}_s$ is the Cauchy stress tensor, and ${\mathit{f}}_s$ is the external body force acting on the structure. There is no stress and strain on the structure at the initial time.

The canopy material is simulated by a fabric with large deformation and nonlinear dynamic properties. The stress-strain relationship equations of the membrane element are given as (10) $\begin{array}{rl}& {ϵ}_a=\frac{1}{E_a}({\sigma }_a-{\nu }_a{\sigma }_b)\\ & {ϵ}_b=\frac{1}{E_b}({\sigma }_b-{\nu }_b{\sigma }_a)\\ & 2{ϵ}_{ab}=\frac{1}{G_{ab}}{\tau }_{ab}-\text{β}{\tau }_{ab}^3\end{array}$(10) where $ϵ,\nu$ and E represent the stress, Poisson's ratio and elastic modulus, respectively. In this notation subscripts a and b refer to the longitudinal direction and transverse direction, respectively. $G_{ab}$ is the shear modulus, ${\tau }_{ab}$ is the shear stress, and β is a nonlinear coefficient, which can be obtained by stress-strain tests.

A discrete cable element in the LS-DYNA material library is used for the suspension lines. The force generated by the suspension lines $F_c$ is non-zero only if the line is in tension, and the tensile force can be calculated by the cable material stiffness and offset length. Considering the damping and nonlinear characteristics of ropes, the dynamic governing equations are (11) $\begin{array}{rl}& F_c=f\left(ϵ\right)+\eta \dot{ϵ}(ϵ>0)\\ & ϵ=\frac{l-(l_0-\mathrm{\Delta }l)}{(l_0-\mathrm{\Delta }l)}\end{array}$(11) where $f\left(ϵ\right)$ denotes the nonlinear tensile function of ropes, η is the damping coefficient, ε is the strain of ropes, $\mathrm{\Delta }l$ denotes the element offset, $l_0$ and l represents the initial length and current length of the cable element, respectively.

The above equations are discrete in time. The discrete forms of the above structural equations are based on the linear time discontinuous Galerkin strategy. The specific details can be found in Bauchau (2011).

3.2. Fluid mechanics

The structure is subjected to fluid loading, which in turn is affected by structural deformations. This loading is obtained by solving RANS equations in arbitrary Lagrangian-Eulerian coordinate system (Lee & Xu, 2016; Winter & Sváček, 2021). The governing equations for an incompressible, viscous fluid flow are given as (12) $\begin{array}{rl}& {\rho }_f(\frac{\mathrm{\partial }{\overline{\mathit{v}}}_f}{\mathrm{\partial }t}+({\overline{\mathit{v}}}_f-{\overline{\mathit{v}}}_m)\cdot \mathrm{\nabla }{\overline{\mathit{v}}}_f)=\mathrm{\nabla }\cdot \overline{\mathit{\sigma }}+\mathrm{\nabla }\cdot {\mathit{\sigma }}^{\mathrm{T}}+{\rho }_f{\mathit{f}}_f\end{array}$(12) (13) $\begin{array}{rl}& \mathrm{\nabla }\cdot {\overline{\mathit{v}}}_f=0\mathrm{i}\mathrm{n}{\mathrm{\Omega }}_f\end{array}$(13) where ${\rho }_f$ is the fluid density, ${\mathit{f}}_f$ is the body force, ${\overline{\mathit{v}}}_f$ and ${\overline{\mathit{v}}}_m$ are the fluid velocity and mesh velocity respectively, ${\mathit{\sigma }}^T$ and $\overline{\mathit{\sigma }}$ are the turbulent stress term and Cauchy stress tensor respectively, which are given by (14) $\begin{array}{rl}\overline{\mathit{\sigma }}& =-\overline{p}\mathbf{I}+\mu (\mathrm{\nabla }{\overline{\mathit{v}}}_f+{\left(\mathrm{\nabla }{\overline{\mathit{v}}}_f\right)}^{\mathrm{T}})\end{array}$(14) (15) $\begin{array}{rl}{\mathit{\sigma }}^T& ={\mu }_t\overline{(\mathrm{\nabla }\cdot {\mathit{v}}_f^{\mathrm{\prime }}+{(\mathrm{\nabla }\cdot {\mathit{v}}_f^{\mathrm{\prime }})}^{\mathrm{T}})}\end{array}$(15) where $\overline{p}$ is the time-averaged pressure, $\mathbf{I}$ is the unit vector, ${\mathit{v}}_f^{\mathrm{\prime }}$ is the fluctuating component of velocity resulting from the turbulence, μ and ${\mu }_t$ are the fluid and turbulent dynamic viscosity respectively. The initial conditions and boundary conditions are required to supplement Equations (12(12) $\begin{array}{rl}& {\rho }_f(\frac{\mathrm{\partial }{\overline{\mathit{v}}}_f}{\mathrm{\partial }t}+({\overline{\mathit{v}}}_f-{\overline{\mathit{v}}}_m)\cdot \mathrm{\nabla }{\overline{\mathit{v}}}_f)=\mathrm{\nabla }\cdot \overline{\mathit{\sigma }}+\mathrm{\nabla }\cdot {\mathit{\sigma }}^{\mathrm{T}}+{\rho }_f{\mathit{f}}_f\end{array}$(12) ) and (13(13) $\begin{array}{rl}& \mathrm{\nabla }\cdot {\overline{\mathit{v}}}_f=0\mathrm{i}\mathrm{n}{\mathrm{\Omega }}_f\end{array}$(13) ). The normal velocity of slip-wall equals zero; the tangential velocity of nonslip-wall is zero; the tangential stress of the free surface must vanish, and the normal stress is balanced with the external applied normal force. According to these constraints, the boundary conditions are defined as (16) $\begin{array}{rl}{\overline{\mathit{v}}}_f& ={\mathit{v}}_f^D\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}_f^D\end{array}$(16) (17) $\begin{array}{rl}{\overline{\mathit{\sigma }}}_f\cdot \overrightarrow{n}& ={\mathit{t}}_f^N\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}_f^N\end{array}$(17) where ${\mathit{v}}_f^D$ denotes the imposed Dirichlet condition of velocity on ${\mathrm{\Gamma }}_f^D$, ${\overline{\mathit{\sigma }}}_f$ is the stress normal to the boundary surface, with $\overrightarrow{n}$ being the outward normal to the surface, and ${\mathit{t}}_f^N$ is the externally applied Neumann condition of normal stress on ${\mathrm{\Gamma }}_f^N$.

In this paper, the Spalart-Allmaras model (Spalart & Allmaras, 1992) is chosen as the turbulence model, as it is designed specifically for aerospace applications involving wall-bounded flows and has been shown to give good results for boundary layers of 2-D airfoils (Mohammadi & Johari, 2010; Winslow et al., 2018) and 3-D parafoils (Wu et al., 2019).

3.3. Coupling algorithm

To accurately capture the FSI phenomenon between flexible parafoils and flow field, it is necessary to couple the fluid and structure solvers running simultaneously. Owing to the FSI coupling capability of LS-DYNA can combine the incompressible CFD solver with its structural mechanics solver, the solution efficiency is greatly improved (Garrec et al., 2019).

Fluid-structure interactions are performed at the interface ${\mathrm{\Gamma }}_I$ until convergence is reached between fluid and structural domains. The ALE approach for mesh movement of the LS-DYNA code enables the coupling to occur with non-coincident fluid and structure meshes (Caldichoury et al., 2014). The convergence is achieved through the balance of dynamic and kinematic coupling conditions. Kinematic constraints guarantee identical velocities at the fluid-structure interface, while dynamic constraints assure stress stability of the interface. (18) $\begin{array}{rl}{\mathit{v}}_s& ={\overline{\mathit{v}}}_f\mathrm{a}\mathrm{t}{\mathrm{\Gamma }}_I\end{array}$(18) (19) $\begin{array}{rl}{\mathit{\tau }}_s\cdot \overrightarrow{n}& ={\overline{\mathit{\sigma }}}_f\cdot \overrightarrow{n}\mathrm{a}\mathrm{t}{\mathrm{\Gamma }}_I\end{array}$(19) where ${\mathit{\tau }}_s$ is the tensile force acting on structure along with the interface, and $\overrightarrow{n}$ is the outer normal of fluid-structure coupling interface. The force acting on the structure is determined by fluid loading on the interface. The above conditions are always met over time.

Flow chart of the two-way coupling FSI strategy is presented in Figure 2. First, the fluid domain calculates initial aerodynamic forces based on initial conditions. The pressure data are then passed to the structural field, and deformations of the canopy under the fluid loading are calculated. And then, the structure domain transfers displacements of the updated canopy boundary to the flow field. After aerodynamic force and structural displacements converge at the current time step, the solution proceeds to the next step. The above process repeatedly iterates until the residuals are below a preset value (usually ${10}^{-3}$) and no longer vary.

Figure 2. Flow chart of computational algorithm.

4. Computational model

4.1. Parafoil models

As shown in Figure 3, the parafoil is a two-layer fabric structure consisting of a series of connected cells sewn together. The leading edge inlet is open, allowing air to enter the interior, and a ram-air wing shape is formed to obtain aerodynamic forces. The cells are separated by ribs, as shown in Figure 4. A rib has several holes to keep pressure balanced between cells. The leading edge inlet of the rib has a significant influence on the aerodynamic performance of parafoil. Based on different dimensions of the leading edge inlet and canopy area, we generated five parafoil models to analyze the influence of geometric parameters on aerodynamic performance. Detailed geometric parameters of parafoil are summarized in Table 1. Model B is obtained by proportionally scaling half of the area of Model A, which aims to analyze the influence of canopy area scaling on aerodynamic performance. Model C and Model D have different leading edge inlet sizes to investigate the effect of inlet length and cut angle on aerodynamic performance. In this paper, airdrop experiments were conducted using parafoil systems based on Model A and Model E to verify the numerical prediction results of the gliding performance.

Figure 3. Computer aided design model of parafoil system.

Figure 4. A rib.

Table 1. Dimension parameters of parafoil models.

Parameter	Model A	Model B	Model C	Model D	Model E
Chord length c (m)	0.82	0.6	0.6	0.6	3.18
Span length (m)	2.2	1.6	1.6	1.6	10.5
Aspect ratio	2.7	2.7	2.7	2.7	3.3
Leading inlet length (m)	0.1	0.1	0.03	0.06	0.36
Leading inlet angle (${}^{\circ }$)	157	157	147	147	135
Relative thickness	0.12	0.12	0.12	0.13	0.12
Average rope length (m)	1.5	1.1	1.1	1.1	6.8
Canopy area (m${}^2$)	2.0	1.0	1.0	1.0	34.0

4.2. Material parameters

All FSI simulations were conducted at standard sea-level conditions with a Reynolds number of $1.0\times {10}^6$. The incoming flow has a density of 1.225 kg/m${}^3$ and a viscosity of $1.7895\times {10}^{-5}\mathrm{P}\mathrm{a}\cdot \mathrm{s}$. The initial velocity of the incoming flow was 15 m/s, and the reference pressure was 101,325 Pa. The parafoil canopy and suspension lines were modeled by membrane and beam elements, respectively. The relevant material parameters are shown in Table 2.

Table 2. Material properties.

Item	Membrane	Cable
Young's modulus (Pa)	$4.31\times {10}^8$	$9.70\times {10}^{10}$
Density (kg/m${}^3$)	533.77	462.00
Poisson's ratio	0.14	–
Thickness (m)	$1.00\times {10}^{-4}$	–
Cross-section area (m${}^2$)	–	$4.91\times {10}^{-6}$

4.3. Computational mesh

As shown in Figure 5, a finite element mesh model with suspension lines of different lengths was established for the structural domain. One end of the suspension line converges to a fixed point, while the other end is connected to the lower surface of the canopy to restrict deformations effectively.

Figure 5. Initial mesh model of parafoil.

And a mesh model of the flow field was established as shown in Figure 6. The size of the flow field must be large enough to avoid boundary effects (Wu et al., 2019). The boundary conditions of velocity inlet, pressure outlet, and free-slip wall are represented in Figure 6, respectively, and the canopy surface was set as the non-slip boundary condition.

Figure 6. Mesh model of fluid domain (c denotes the chord length).

4.4. Mesh convergence

First, the grid sensitivity of parafoil models was studied. Three mesh models with different mesh sizes were used to calculate aerodynamic forces at an angle of attack (denoted by α) 6${}^{\circ }$. Each case was carried out in 20,000 iterations with a time step of $1.0\times {10}^{-4}$. Moreover, the values of the last 1000 steps were averaged. The obtained results are listed in Table 3. It can be seen that the difference between the results of medium mesh and fine mesh models was no more than 0.8%. Compromising accuracy and computational cost, we chose the medium grid sizes as a reference in the following simulations.

Table 3. Grid sensitivity examination.

Mesh	Number of grids	Number of nodes	$C_L$ (error)	$C_D$ (error)
Fine	4216197	876428	0.582	0.118
Coarse	1790759	448657	0.569 $(-2.2\mathrm{\%})$	0.106 $(-10.2\mathrm{\%})$
Medium	3346297	724289	0.584 $(+0.3\mathrm{\%})$	0.119 $(+0.8\mathrm{\%})$

5. Fluid-structure interaction simulation and discussion

Based on the FSI frameworks established in Section 3, and taking the uninflated smooth parafoil shape in Figure 3 as the initial model, we conducted FSI simulations under different angles of attack in the range of $0^{\circ }$–${24}^{\circ }$. The calculations were performed on TianHe-1 at National Supercomputer Center in Tianjin and on 28 cores belonging to a single compute node. Each case was carried out in 20,000 iterations and took 24 h using the medium mesh size in Table 3.

5.1. Structural deformation

Aerodynamic performance is closely related to canopy shape, so precise aerodynamic force predictions depend on accurate canopy geometries. Figure 7 displays the shapes of the Model A canopy before and after inflating. In detail, Table 4 presents geometric dimension changes of Model A before and after inflation, where $T_{max}$ means the maximum thickness of cells. It can be seen that the bulge of cells caused by inflation led to a 30% increase in thickness. As such, the canopy size decreased along the chordwise and spanwise directions. We can conclude that the dimension changes of the canopy before and after inflating are apparent and cannot be ignored. In addition, stabilizers and the front of the ribs also appeared to have slight deformations (see green arrows), which may be caused by the lack of constraints on the canopy in the numerical model.

Figure 7. Isometric views of the canopy (a) before and (b) after inflation.

Table 4. Geometric dimension changes of parafoil canopy caused by inflation.

Model A	Span	Chord	Area	Aspect ratio	$T_{max}$
Initial geometry	2.20 m	0.82 m	2.00 m${}^2$	2.70	0.10 m
Inflated canopy	1.90 m	0.75 m	1.50 m${}^2$	2.50	0.13 m
Relative shrinkage	−13.60%	−8.50%	−25.00%	−7.40%	+30.00%

5.2. Flow field simulation analysis

This section presents detailed fluid field results from fluid-structure interaction simulations of parafoil models with different characteristic geometries, leading edge inlet sizes, and canopy area scalings. CFD simulations (without coupling with the structural domain) were also run and compared with FSI simulations to analyze the influence of inflation deformations of the structural domain. The aerodynamic coefficients predicted by the flow field of FSI simulations supplement the flight performance predictions.

5.2.1. Effect of parafoil inflation

From Section 5.1, we know that the flexibility of canopy structure leads to deformations along spanwise and chordwise directions, which also affects the aerodynamic performance of parafoil. This section then evaluates the differences in aerodynamic response by comparing the flexible deformation model with the rigid body model.

Figure 8 presents external pressure contours and streamlines on the upper surfaces of parafoils. The flexible and rigid models showed different pressure distributions in the spanwise direction. Before the chordwise maximum thickness line from the leading edge, the negative pressure area in the middle cells decreased gradually and disappeared toward the end-tip cells. This may be due to the arc dihedral angle of the canopy along the spanwise direction. After flow passed through the maximum thickness point, pressure increased gradually toward the trailing edge, consistent with the NASA model simulation results (Cao & Zhu, 2013). In addition, the high-pressure area near the trailing edge of the upper surface of the flexible parafoil model was smaller than that of the rigid model, which may be caused by the change of the arc dihedral angle due to the inflation deformations.

Figure 8. Streamlines around upper surface of (a) Inflated flexible model and (b) Rigid model with pressure contour (Direction from bottom to top).

The flow behaviour on the upper surface in Figure 8 shows that streamlines around end-tip cells began to tilt from outside toward the center of the canopy. In contrast, the streamlines showed the opposite direction when flowing through the lower surface in Figure 9. This phenomenon resulted from the sharp-edge effects, which further created wake vortex cores at the trailing edge of the finite wing (Ghoreyshi et al., 2016).

Figure 9. Streamlines around lower surface of (a) Inflated flexible model and (b) Rigid model with pressure contour.

Figure 10 compares the lift and drag coefficient curves of the rigid and flexible models, with angles of attack ranging from 0${}^{\circ }$ to 20${}^{\circ }$. Lift coefficients increased approximately linearly in the range of small angles of attack before the stall, where the stall implies a sudden decrease in the lift coefficients and an increase in the drag coefficients. The maximum lift coefficient of the rigid model was 0.73 at a stall angle of attack 12${}^{\circ }$. In comparison, the maximum lift coefficient of the flexible model was 0.85 at a stall angle of attack 14${}^{\circ }$, which is significantly higher (by 16.4%) than that of the rigid model. Moreover, the lift coefficients of the flexible model were slightly larger than that of the rigid model at small angles of attack, and the differences became evident after the stall. After reaching stall conditions, lift coefficients suddenly dropped and changed nonlinearly. The drag coefficients consistently increased nonlinearly. There were few differences between rigid and flexible models under small angles of attack. The growth slope of the rigid model was greater, resulting in differences between the two models becoming more evident.

Figure 10. Comparison of lift and drag coefficients between flexible and rigid model.

5.2.2. Effect of airfoil inlet

The basic airfoil inlet is a crucial factor affecting the aerodynamic performance of parafoil (Cao & Zhu, 2013). To investigate the influence of airfoil shape on the aerodynamic performance of a parafoil system, we constructed two canopy models composed of airfoils with different inlet lengths, named Model C (inlet length 0.05c, c denotes the chord length) and Model D (inlet length 0.1c), respectively. For specific geometric parameters, please refer to Table 1.

For further comparison, Figures 11–13 show the pressure contours of slice sections near the central ribs of two models at different angles of attack. The negative pressure region (blue) near maximum thickness along the chordwise direction corresponded to accelerated airflow flowing through upper surfaces. The longer the incision, the larger the separation bubble at the leading edge. As the increase of angle of attack, high-pressure regions on the lower surfaces of the two models expanded gradually. Therefore, pressure differences between upper and lower surfaces increased, and so did lift forces. As shown in Figure 12, compared to Model C, Model D had a larger positive pressure area (red) on the lower surface, so the pressure difference between upper and lower surfaces, i.e. the lift, was greater. It proves that a parafoil with a bigger leading edge inlet has a higher lift. When the angles of attack increased to ${14}^{\circ }$, as shown in Figure 13, there was a slight pressure difference between the front and rear of the interior cells, especially for Model C. However, severe fluid separation has occurred on the upper surface of Model D, which may seriously deteriorate the aerodynamic performance. It reveals that the aerodynamic performance of Model D deteriorated earlier.

Figure 11. Pressure contours along chordwise of (a) Model C and (b) Model D at $\alpha =2^{\circ }$.

Figure 12. Pressure contours along chordwise of (a) Model C and (b) Model D at $\alpha =6^{\circ }$.

Figure 13. Pressure contour along chordwise of (a) Model C and (b) Model D at $\alpha ={14}^{\circ }$.

The vorticity contours in Figures 14–16 present the variation of boundary layer behaviours with the angle of attack for two parafoil models. When $\alpha =4^{\circ }$, a thin upper boundary layer increased gradually after the maximum thickness point of the rib as airflow flowed toward the trailing edge, while the lower boundary layer increased from the leading edge inlet. With the increase of angles of attack, the thickness of the lower boundary layer decreased, indicating that fluid fitted more closely to the airfoil surface and improved the aerodynamic performance. When the angle of attack increased to ${10}^{\circ }$, the boundary layer started earlier from the upper leading edge. Comparing the boundary layers in Figures 15 and 16, the upper boundary layer thickness increased with the angle of attack, while the lower layer thickness reduced. At $\alpha ={16}^{\circ }$, the upper boundary layer separation of Model D was more severe than Model C, which confirmed that the parafoil with a smaller inlet has a better anti-stall performance.

Figure 14. Vorticity magnitude contours along central cross-sections of Model C (left) and Model D (right) at $\alpha =4^{\circ }$.

Figure 15. Vorticity magnitude contours along central cross-sections of Model C (left) and Model D (right) at $\alpha ={10}^{\circ }$.

Figure 16. Vorticity magnitude contours along central cross-sections of Model C (left) and Model D (right) at $\alpha ={16}^{\circ }$.

To quantitatively evaluate the flow behaviour observed above, Figure 17 shows the aerodynamic parameters of two models in the range of $\alpha =0^{\circ }$ to 24${}^{\circ }$.

Figure 17. Comparison of (a) Lift coefficients, (b) Drag coefficients and (c) Lift to drag ratios between Model C and Model D.

The lift coefficients $C_L$ and drag coefficients ${\mathrm{C}}_{\mathrm{D}}$ of Model D were larger than Model C in the whole valid range of α, and Model D stalled earlier (about $\alpha ={14}^{\circ }$). This result can be explained by more severe flow separations in Model D observed in Figure 16. On the contrary, the lift-to-drag ratio L/D of Model C was larger than Model D in Figure 17(c). The maximum L/D of Model C occurred around $\alpha =6^{\circ }$, which was about 22% higher than that of Model D at $\alpha =8^{\circ }$. It shows that the effect of drag is more evident than lift, so the improvement of parafoil gliding performance should be mainly focused on reducing drag of all types.

5.2.3. Effect of canopy scaling

This section studies the effect of canopy area scaling on aerodynamic performances by comparing Model A of 2.0-m${}^2$ canopy and Model B of 1.0-m${}^2$ canopy. The results are plotted in Figure 18. With the canopy area reduced by half, the growth slope of lift coefficient $C_L$ increased, and the stall angle of attack was delayed from 14${}^{\circ }$ to 16${}^{\circ }$. The growth rate of drag coefficient $C_D$ changed slightly, which shows that scaling down the canopy area has a more significant impact on lift than drag. From Figure 18(b), the maximum L/D increased from 5.1 (at $\alpha ={10}^{\circ }$) to 5.4 (at $\alpha =8^{\circ }$), an increase of about 6%. It reveals that reduction of canopy area could improve gliding performance in the range of affordable payload weight of parafoil.

Figure 18. Comparison of (a) Lift and drag coefficients and (b) Lift to drag ratio between full-scale Model A and scaled-down Model B.

From the above analysis in Section 5.2, the effects of inflating deformation, basic airfoil inlet, and canopy area scaling on the aerodynamic performance of parafoil were analyzed. The error of lift coefficient between the flexible inflated and rigid body models was up to 16%, proving that the inflation deformation has a significant influence on the aerodynamics of parafoil. Reducing the inlet length of parafoil can improve the gliding performance and increase the range of stall angles of attack. The canopy area scaling down could be beneficial to improve the lift-to-drag ratio. These simulation results provide some references for the parafoil design industry, where the geometric shape of parafoil can be adjusted in the design stage to meet application requirements.

6. Glide performance prediction and airdrop tests validation

6.1. Glide performance prediction

During actual flight, the installation angle can be adjusted by the length of suspension ropes so that the parafoil system can be steady in a gliding state. Figure 19 shows the predicted velocities of Model A parafoil system with different suspended weights. Wing loading is the ratio of suspended weight to wing area. Horizontal and descent velocities increased with wing loading–the higher the attack angle, the smaller the velocities. The horizontal velocity can reach 8 m/s to 10 m/s with a payload in the range of 10 kg to 20 kg. With a 13 kg payload, horizontal and descent velocities at a steady angle of attack $\alpha =6^{\circ }$ were about 9.9 m/s and 4.7 m/s, respectively.

Figure 19. (a) Horizontal and (b) Vertical velocities vs payload weight.

6.2. Airdrop tests validation

To verify the flight performance prediction capability of established parafoil models in this paper, we conducted airdrop tests on parafoil systems of Model A with a 13 kg payload and Model E with a 100 kg payload, respectively.

6.2.1. Airdrop test for model A

The airdrop system as shown in Figure 20 began to glide at time t = 65 s. Figure 21(a) plots velocity change curves of the experimental system within $t=80\mathrm{s}$, and the corresponding trajectory is plotted in Figure 21(b). In almost windless experimental conditions, the gliding velocities of the parafoil system were about 10.1 m/s in the horizontal direction and 5.0 m/s in the vertical direction.

Figure 20. Airdrop experiment of Model A.

Figure 21. (a) Velocity curves and (b) Trajectory of airdrop test for Model A.

Based on aerodynamic coefficients calculated by FSI simulations for the flexible parafoil model and by CFD simulations for the rigid body model, we can use the numerical model to predict the gliding velocities of the parafoil system under the same experimental conditions. The specific results are listed in Table 5. By comparing numerical results with the experimental data, we can see that prediction errors for velocity components of the flexible model were not exceed 6%, while the errors of the rigid model reached 11%. Furthermore, the prediction error for the glide ratio of the flexible model was up to 10.5% lower than that of the rigid model. These results demonstrate the feasibility and accuracy of the proposed FSI model.

Table 5. Flight performance comparisons between simulation and airdrop testing of Model A.

Parameters	Airdrop test	Flexible model (error)	Rigid model (error)
Horizontal velocity (m/s)	10.10	9.94 (1.6%)	11.20 (11.0%)
Vertical velocity (m/s)	5.00	4.72 (5.6%)	4.82 (3.6%)
Glide ratio	2.02	2.11 (4.5%)	2.32 (15.0%)

6.2.2. Airdrop test for model E

As shown in Figure 22, the airdrop system of Model E began to glide at $t=40\mathrm{s}$. Figure 23(a) plots the velocity change curves of the experimental system within $t=130\mathrm{s}$, and the corresponding trajectory is plotted in Figure 23(b). In almost windless conditions, the gliding velocities of the parafoil system were about 9.69 m/s in the horizontal direction and 3.24 m/s in the vertical direction. The specific velocity comparisons and errors between the airdrop test data and simulation results of the flexible model are listed in Table 6, and the errors are also within a reasonable range. This demonstrates the universality of the FSI model in this paper for the aerodynamic predictions of different parafoil geometries.

Figure 22. Airdrop experiment of Model E.

Figure 23. (a) Velocity curves and (b) Trajectory of airdrop test for Model E.

Table 6. Velocity comparisons between simulation and airdrop testing of Model E.

Velocities	Airdrop test	Flexible model	Errors
Horizontal velocity (m/s)	9.69	10.19	5.10%
Vertical velocity (m/s)	3.24	3.11	4.00%
Glide ratio	3.00	3.27	9.00%

7. Conclusion

This paper presents a fluid-structure interaction simulation method for predicting and optimizing the aerodynamic and flight performance of parafoils in the design stage. The visualization and analysis of various parafoil geometry shapes and flow behaviour, combined with material properties, make flight performance predictions more reliable by addressing the challenges associated with accurate representation of the parafoil geometry and aerodynamic characteristics. The main findings of this paper are as follows:

1. There were significant differences between the simulation results of flexible and rigid parafoil models. The main features of canopy inflation shape are bulge of cells and slight deformation of ribs. The results show a 16.4% increase in the maximum lift coefficient compared to the rigid model.

2. The parafoil with a smaller leading edge inlet or scaling down area could improve the aerodynamic performance, mainly manifested by higher lift-to-drag ratio and delayed stalling angle of attack. Improving the flight speed and the range of manipulable angles of attack can be achieved by appropriately reducing the inlet length or area of the parafoil, considering the inflating feasibility of the parafoil and the payload weight.

3. The prediction results of flight performance in steady gliding were consistent with airdrop test data, and the accuracy of the flexible model was more than 10% higher than that of the rigid model–this proved applicability of the proposed method for predicting the actual performance of parafoil systems.

In future applications, the proposed method can be used to predict the dynamic response of parafoil under control line retraction maneuvers, aiming to close the gap between simulation and experiment further.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the National Natural Science Foundation of China (Grant NoS. 61973172, 61973175, 62003175 and 62003177), the Key Technologies Research and Development Program of Tianjin (Grant No. 19JCZD-JC32800), and Research and Innovation Project for Postgraduates in Tianjin (Grant No. 2021YJSO2B02), this project is also funded by the China Scholarship Council.