Accurate calculation of aerodynamic coefficients of parafoil airdrop system based on computational fluid dynamic

Abstract

Accurate calculation of canopy aerodynamic parameters is a prerequisite for precise modeling of a parafoil airdrop system. This investigation analyses the aerodynamic performance of the canopy in airdrop testing combining the leading-edge incision and the trailing-edge deflection. Aerodynamic parameters of the canopy are obtained using the computational fluid dynamic simulations, and then, the output data are used to estimate the deflection and incision factors. The estimated lift and drag coefficients instead of the traditional parameters based on lifting-line theory are incorporated into the eight degrees of freedom dynamic model of an airdrop system and make some simulations. The effectiveness of the proposed method for calculating aerodynamic coefficients is verified by actual airdrop testing.

Introduction

In recent years, the parafoil airdrop system (PADS) gradually becomes active in the field of aerospace. It has the abilities of patrolling, military reconnaissance, surveillance, advertising, and so on. Compared with the fixed-wing aircraft and rotor unmanned aerial vehicle (UAV), the PADS has the advantages of high loading ratio and low cost. ¹ The Aviation Industry Corporation of China developed a kind of PADS,with the the thrust-to-weight ratio reaching 1:3, whereas the thrust-to-weight ratio of the other types of UAV was just 1:1. Moreover, the PADS can fly at low speed and low altitude, ² the flight speed of the parafoil is 30–70 km/h, and the flight height ranges from several meters to several thousand meters. In addition, the PADS is easy to be carried by aircraft and missile and can be delivered remotely. ³

The canopy of the parafoil is entirely constructed from fabric with no rigid members, which allows it to be packed and deployed. When inflated, it is just like a low-aspect-ratio wing, containing the upper and lower surfaces, a cross-section, and a rectangular platform. The ribs are sewn in the chord direction at a number of spanwise intervals forming a series of cells. The leading edge of the wing is open over its length so that the ram-air pressure maintains the wing shape. ⁴ The aerodynamic characteristics of the parafoil have two distinctive features: one is the strong nonlinearity and coupling due to the flexibility of the canopy and another is that the leading-edge incision and trailing-edge deflection (brake deflection) of the canopy are hard to incorporate into analytical models. These two elements make it difficult to accurately model for parafoils, which pose a new challenge to the control of the PADS.

The modeling and control are the most significant problems in the airdrop sysems. Xu and Chen ⁵ made a review of modeling and control during the airdrop process of transport aircraft, in which the airdrop systems were classified into two scenarios according to the airdrop height. The PADSs belong to the low altitude airdrop, which mainly focus on the control and aerodynamic characteristics to ensure safety; therefore, it is important to compute the aerodynamics of the canopy accurately. In the past decades, most investigations have been conducted on the mechanism modeling of parafoils. Goodrick ⁶ and Barrows ⁷ established asix-degree-of-freedom (DOF) of PADS, wherein the parafoil and payload were regarded as rigid connection. Müller et al. ⁸ proposed an eight-DOF model of parafoil systems, in which the relative motions between airfoil and payload were taken into account. Zhu et al. ⁹ yielded an eight-DOF dynamic model of the powered parafoil according to the Kirchhoff motion equation. Meanwhile, the nine-DOF models of parafoils were also investigated by some researchers. ^{10 –13} It is valuable to mention that all the abovementioned studies computed the aerodynamic parameters of the canopy based on the lifting-line theory, ¹⁴ which involves a series of assumptions to account for the effects of the leading-edge incision and the brake deflection of the canopy.

To control the airdrop systems, Xu ¹⁵ proposed a dynamic surface control of nonlinear systems, which can deal with the uncertainty. There are also many others control methods that can deal with the uncertainty of the PADSs, for example, neural control ^16,17 and learning control. ¹⁸ However, too much uncertainty in the system always makes the control algorithm unsatisfactory. In order to control a parafoil system more accurately, a precise mathematical model that includes the kinematics and aerodynamics is more desirable. Nowadays, three methods can be used to accurately compute the aerodynamic force, which are the airdrop testing, the wind–tunnel experiment, and the computational fluid dynamic (CFD). ¹⁹ Airdrop testing has high cost, and the measured data are limited, which is often used to verify experiments. Owing to large dispersion and poor regularity, it is difficult to utilize the wind–tunnel testing data for simulation analysis. With the development of computers, CFD is an attractive method due to its high reproducibility and low cost.

The existing studies based on CFD mostly focused on the qualitative analysis, and the investigations on the quantitative analysis of the effect of leading-edge incision and trailing-edge deflection were rarely carried out. Ghoreyshi et al. ²⁰ investigated the aerodynamic characteristics of ram-air parachute with open and closed round inlets under the assumption of rigid canopy and compared the results with available wind–tunnel test data. Based on CFD, two estimation methods were used to calculate stability derivatives. Ghoreyshi et al. ²¹ conducted a numerical simulation on parafoils that had plate form. The airflow regime on canopy was also analyzed in detail. The presence of a complex vortex system inside the air chamber and in the proximity of the inlet was revealed by Andreutti et al. ²² using CFD solvers. Moreover, it pointed out that the Euler simulation was able to describe the macroscopic flow features. In order to study the unique aerodynamic characteristics of the PADS, Altmann ²³ combined loads with flexible deformations to investigate the effects of spanwise compression of the wing, cell billowing, and S-shapening of the profiles on aerodynamic characteristics. He also proposed some positive proposals on the design of parafoil systems. Zhu and Cao ²⁴ investigated the influence of arc-angle, leading-edge cut and the plate form of the parafoil on the aerodynamic performance. The abovementioned simulations were carried out from the perspective of the parafoil design, and the modeling and control of the PADS were not considered. Liang et al. ²⁵ studied the modeling of a parafoil in the wind and rain environments. The flap deflection was also investigated, and the modified aerodynamic function was obtained using least squares method. ²⁶ Tao et al., ²⁷ by combing the CFD and the Kirchhoff motion equation, proposed a more accurate model of the parafoil. However, in these models, the leading-edge incision and the ribs were neglected.

To meet the demand of PADS control, accurate modeling is essential. Leading-edge incision and trailing-edge deflection have a significant impact on the precision of the modeling of the PADSs, which were neglected in the traditional method of modeling. This problem highlights the motivation for this investigation, and the ribs of the parafoil are also taken into account in this article. Moreover, the quantitative analysis of the influence of the leading-edge incision is conducted. Furthermore, trailing-edge deflection (the only control input of PADS) is investigated. By solving the Navier–Stokes equation, the aerodynamic parameters of the canopy can be obtained and then is used to identify the aerodynamic function. Finally, the aerodynamic functions are verified by the airdrop testing.

The remainder of this article is organized as follows. “Problem formulation” section describes the origin and formulates the problem. In “Computational approach of aerodynamic coefficients” section, the process of numerical simulations based on CFD is introduced in detail. The next section is to ensure the effectiveness of numerical simulations. “Simulations and analysis” and “Airdrop testing” sections provide the simulation results to verify the proposed approach.

Problem formulation

To calculate the aerodynamic force, the traditional method is dividing the canopy of the parafoil into eight segments evenly along the span direction. ²⁸ The lift coefficient on every segment times a factor from outside to inside (0.6, 1.0, 1.16, 1.24). According to the local velocity and angle of attack (AOA), the aerodynamic force and moment on each segment are computed. The entire aerodynamic force and moment are obtained through the addition of each segment. The aerodynamic on every segment of a parafoil is computed as follows

$\{\begin{array}{c}{F}_{Li}=k_i{c}_{Li}\left(0.5\rho S_i\sqrt{u_i^2+w_i^2}\right)\left[\begin{array}{c}{w}_i\\ 0\\ -u_i\end{array}\right]\\ F_{Di}=-c_{Di}\left(0.5\rho S_i{V}_i\right)\left[\begin{array}{c}{u}_i\\ 0\\ w_i\end{array}\right]\end{array}$

1

where the subscript $i\left(i=1,2,\dots \mathrm{,8}\right)$ denotes the number of segment; k is the product factor, S represents the area of each segment; ρ is the air density; and u, v, and w denote the components of the velocity vector in its coordinate system at the pressure center, respectively. The lift and drag coefficients are often computed based on the lifting-line theory and given as follows

$\{\begin{array}{c}cl=C_{\text{L}\alpha }\left(\alpha -{\alpha }_0\right){{\displaystyle \text{cos}}}^2\left(\varphi \right)+j{{\displaystyle \text{sin}}}^2\left(\alpha -{\alpha }_0\right)\text{\hspace{0.17em}cos}\left(\alpha -{\alpha }_0\right)\\ cd=C_{\text{D}0}+\frac{C_{LC}^2{\left(\alpha -{\alpha }_0\right)}^2\left(1+\delta \right)}{\pi \text{AR}}+j{{\displaystyle \text{sin}}}^3\left(\alpha -{\alpha }_0\right)\end{array}$

2

In the determination of $cl$ , $C_{\text{L}\alpha }$ denotes the lift-curve slope, α is the AOA, ${\alpha }_0$ is the zero lift angle, and ϕ is the arc anhedral angle. The lift-curve slope has the following form

$C_{\text{L}\alpha }=\frac{\pi C_{\text{L}\alpha }^a\text{AR}}{\pi \text{AR}+C_{\text{L}\alpha }^{\alpha }\left(1+\tau \right)}$

3

where $C_{\text{L}\alpha }^a$ is the two-dimensional lift-curve slope of a parafoil with small aspect ratio, AR is the aspect ratio, and τ is a small positive factor.

$C_{L\alpha }^a=2\pi \text{ARtanh}\left(\frac{C_{L\alpha }^{a_0}}{2\pi \text{AR}}\right)$

4

where $C_{L\alpha }^{a_0}$ is the two-dimensional lift-curve slope and the theoretical value is $2\pi$ , $\Delta C_{\text{L}}=j{{\displaystyle \text{sin}}}^2\left(\alpha -{\alpha }_0\right)\text{cos}\left(\alpha -{\alpha }_0\right)$ increases with AOA over and above the basic slope. This increase is a function of aspect ratio, the shape of the wings lateral edges, and the component of velocity normal to the wing. The nonlinear component has already been investigated but is still not well understood, where j is a function of AR

$j=\{\begin{array}{ll}3.33-1.33\text{AR},\hfill & 1<\text{AR}<2.5\hfill \\ 0,\hfill & \text{AR}>2.5\hfill \end{array}$

5

In the determination of $cd$ , $C_{\text{D}0}$ is the profile drag and the second term represents the induced drag, δ denotes a small factor, and $\Delta D=j{{\displaystyle \text{sin}}}^3\left(\alpha -{\alpha }_0\right)$ is a further drag component.

From the above analysis, although the lifting-line theory can be used to estimate the aerodynamic characteristics of a canopy, the effect of leading-edge cut and the flexibility of the parafoil are neglected, and it is not explicit how to deal with the trailing-edge deflection. For illustrating the influence of the incision, we take three 2-D models as an example. The sketches of the parafoil are shown in Figure 1, r is the rope length, b is the spanwise length, ϕ denotes the arc anhedral angle, c denotes the chord length, e represents the relative thickness, h is the length of the leading edge cut, ε is spanwise radian, and the detailed parameters of three models are presented in Table 1. For three models in Table 1, we compute the lift and drag coefficients based on CFD, and the results are shown in Figure 2.

Figure 1.

Parafoil aerodynamic shape.

Figure 2.

Lift and drag coefficients of three models.

Table 1.

Detailed parameters of 2-D parafoils.

Models	Length (m)	Angle (°)	Relative thickness	Chord length (m)
A	0.36	135	0.12	3.18
B	0.36	150	0.12	3.18
C	0.18	135	0.12	3.18

The effects of the leading-edge cut on the aerodynamic performance are relevant to the length and angle. Figure 3 shows the velocity distribution on the wing surface. The longer the incision, the more serious the speed loss as the airflow passes through the leading edge as well as the more front position of peak velocity in upper wing surface. Under these circumstances, the pressure of the upper surface increases with the deceleration of the airflow; therefore, the aerodynamic performance of the parafoil becomes severe. Different from the incision length effect, the incoming velocity does not change with the increase in the slope of the leading edge, and it will result in airflow separation at the trailing edge and deteriorate the aerodynamic performance. Figures 2 and 3 show the significant effect of the leading-edge incision. Qualitative analysis of the effect of the leading-edge incision is far from enough, and it is crucial for accurate modeling of parafoil systems.

Figure 3.

Velocity distribution of different models: (a) model A, (b) model C, (c) model A, and (d) model B.

Computational approach of aerodynamic coefficients

To calculate the aerodynamic parameters accurately, it is necessary take into account the leading-edge incision and trailing-edge deflection. The primary purpose of this section is the quantitative analysis of these effects based on CFD. Through comparing the pressure and velocity distributions, the flow regime is revealed, and then, the least squares method is used to identify the aerodynamic parameters. The estimated results are compared with the theoretical results to illustrate the effect of the incision and deflection.

Before computing the aerodynamic parameters, some reasonable assumptions are made.

Because we mainly focus on the effect of the trailing-edge deflection and leading-edge incision, the suspension lines are not taken into account in this article.

The main purpose of this investigation is to study the stable flight performance; the effect of the canopy flexibility is neglected.

Governing equation

The flow field of the parafoil satisfies the following control equation

$\frac{\partial \left(\rho \phi \right)}{\partial t}+\text{div}\left(\rho u\phi \right)=\text{div}\left(\Gamma \text{grad}\phi \right)+S$

6

where φ is the generalized variable, ρ denotes the density, u represents the velocity of the free-stream flow, Γ denotes the diffusion coefficient, and S is the generalized source term. In equation (6), from the left to right are the transient term, the convective term, the diffusion term, and the source term, respectively. φ, Γ, and S have particular forms in different flow equations. Because of the parafoil working conditions, it is not necessary to solve the energy function. We give the forms of φ, Γ, and S in continuity equation and momentum conservation equation in Table 2.

Table 2.

Detailed parameters of governing equations.

Function	φ	Γ	S
Continuity equation	1	0	0
Momentum equation	$u_i$	μ	$-\frac{\partial p}{\partial x_i}+S_i$

Parafoil models and flow field settings

The CFD models are created by DesignModeler (version 17.2) in Ansys Workbench. As shown in Figure 4, it has 15 cells, and there are holes on the rib to ensure the pressure balance of each air chamber. Specific parameters of the parafoils are given in Table 3. The boundary of the flow field should be far away enough from the disturbance source; ¹⁴ however, it is not appropriate if the boundary is too far away because too many grids can cause error accumulation. In practice, the dimension of the flow field is set to 10c, 5c, and 10c, as shown in Figure 5.

Figure 4.

Physical model in the airdrop testing: (a) CFD model of parafoil and (b) rib model of parafoil. CFD: computational fluid dynamics.

Figure 5.

Flow field dimension.

Table 3.

Detailed parameters of parafoils.

Models	Length (mm)	Angle (°)	Relative thickness	Chord length (mm)
1	282.73	150.3	0.185c	2011.18
2	282.73	135	0.185c	2011.18
3	100.56	150.3	0.185c	2011.18
4	100.56	135	0.185c	2011.18

The $Re$ number of the parafoil is about $1.4\times {10}^6$ ; therefore, the turbulence flow has to be solved. To ensure the accuracy of simulation results, two turbulence models are utilized here to compute dynamic parameters.

Spalart Allmaras model

The Spalart–Allmaras (SA) model ²⁹ is a model that solves a modeled transport equations for the turbulent viscosity. It is designed specifically for aerospace applications involving wall-bounded flows and has been shown to give good results for boundary layers subjected to adverse pressure gradients. The transport variable in the SA model is given as follows

$\frac{\partial }{\partial t}\left(\rho \tilde{v}\right)+\frac{\partial }{\partial x_i}\left(\rho \tilde{v}{u}_i\right)=G_v+\frac{1}{{\sigma }_{\tilde{v}}} \frac{\partial }{\partial x_j}\left\{\left(\mu +\rho \tilde{v}\right)\frac{\partial \tilde{v}}{\partial x_j}\right\}+C_{b2}\rho {\left(\frac{\partial \tilde{v}}{\partial x_j}\right)}^2]-Y_v+S_{\tilde{v}}$

7

All details of these parameters can be referred in the previous work. ³⁰ Mohammadi and Johari ³¹ combined the incompressible Navier–Stokes equations with the SA turbulence model, examined the flowfield and forces on the central cross section of a parafoil canopy using the finite volume method, and obtained some positive results, which illustrated the efficiency of the SA models.

Renormalization group (RNG) $\text{k}-\epsilon$ model

From the equation of turbulent kinetic energy k, the equation of turbulent dissipation rate ε is introduced. The RNG $k-\epsilon$ model ^32,33 has a similar form with the standard $k-\epsilon$ model, and there are some improvements based on the standard $k-\epsilon$ model. For a class of flows, the RNG $k-\epsilon$ model is more accurate than standard $k-\epsilon$ model. The form of the RNG $k-\epsilon$ model is given as follows

$\frac{\partial }{\partial t}\left(\rho k\right)+\frac{\partial }{\partial x_i}\left(\rho k{u}_i\right)=\frac{\partial }{\partial x_j}\left({\alpha }_k{\mu }_{\text{eff}}\frac{\partial k}{\partial x_j}\right)+G_k+G_b-\rho \epsilon -Y_M+S_k$

8

and

$\begin{array}{l}\frac{\partial }{\partial t}\left(\rho \epsilon \right)+\frac{\partial }{\partial x_i}\left(\rho \epsilon u_i\right)=\frac{\partial }{\partial x_j}\left({\alpha }_{\epsilon }{\mu }_{\text{eff}}\frac{\partial \epsilon }{\partial x_j}\right)\\ +C_{1\epsilon }\frac{\epsilon }{k}\left(G_k+C_{3\epsilon }{G}_b\right)\\ -C_{2\epsilon }\rho \frac{{\epsilon }^2}{k}-R_{\epsilon }+S_{\epsilon }\end{array}$

9

The RNG $k-\epsilon$ model also can give an accurate prediction on the variety of wall reverse pressure gradient. All of the simulations are completed through FLUENT 17.2. Taking the left and bottom of the flow field as the velocity inlets, the fluid is regarded as incompressible air with the incoming velocity of 10 m/s. The front and back sides are set as no-slip wall, and the reminders are pressure outlet. Governing equations have a second-order precision. In addition, the SIMPLIC algorithm is used in this article to deal with the coupling between the velocity and pressure.

Mesh generation

To solve the governing equation, the flow field should be discretized. The mesh generation is finished in integrated computer engineering and manufacturing (ICEM)-CFD, and the air chamber is also taken into account. It is indisputable that structured meshes have higher accuracy and less memory requirement compared with unstructured meshes. Owing to the complexity of parafoil models, tetrahedral meshes are chosen in this investigation.

However, only tetrahedral meshes are not sufficient to capture the high gradients of velocity and pressure near the canopy. There are three layers of the near-wall zone:

Viscous sublayer: It is the closest layer to the wall; in this segment, the viscous force is much higher than the turbulent shear stress. The flow of fluid is similar to laminar flow motion.

Transition layer: Its viscous force and the turbulent shear stress are within the same range. The fluid motion is between turbulent and laminar motion. The transition layer is very thin and is always contributed to logarithmic layer.

Logarithmic layer: It is the outermost layer of the near-wall zone. The viscous force is much smaller than the turbulent shear stress. The flow is in an entirely developed tumultuous state. Velocity distribution is close to logarithmic law.

The SA and RNG $k-\epsilon$ turbulence models are only valid for fully developed turbulence, unable to solve the near wall-zone problem. Up to now, this problem can be solved using two methods, which include wall functions and the low-Re number turbulence model.

Wall function is a semiempirical formula. It is utilized to combine quantities on the wall with the quantities that are unknown in the turbulent core region. To describe the flow in the viscous sublayer and the logarithmic layer, two dimensionless parameters, $u^{+}$ and $y^{+}$ , are introduced, where $u^{+}$ denotes the velocity and $y^{+}$ denotes the distance

$u^{+}=\frac{u}{u_{\tau }}$

10

$y^{+}=\frac{\Delta y\rho u_{\tau }}{\mu }$

11

where u is the time average velocity of the fluid, $u_{\tau }$ is the wall friction velocity, ${\mu }_{\tau }={\tau }_w/\rho$ in which ${\tau }_w$ denotes the wall shear stress, and $\Delta y$ denotes the distance from the first node to the wall surface. It is essential to control the height of the first layer for the desirable results. Many investigations suggested that $y^{+}$ is in the range of 60–300 for standard wall function and $y^{+}=1$ or $y^{+}<4\sim 5$ for enhanced wall function. As for the SA model, either the standard wall function $\left(y^{+}\ge 30\right)$ or enhanced wall function (y⁺ ≥ 1) is appropriate. In this investigation, in SA model, $y^{+}\ge 30$ .

Prism grids are used to treat the boundary layers. The first step of creating computational meshes is to generate the tetrahedral mesh in the whole domain, and then, prism cells are typically extruded from the canopy layer by layer. The first node should be placed at the viscous sublayer, and the height must meet the requirements of wall conditions. The mesh of the flow field and parafoil are shown in Figure 6.

Figure 6.

Mesh generation for flow field and parafoil model: (a) meshes of entire fluid and (b) prism mesh near canopy.

Aerodynamic parameters identification

Taking the leading-edge incision and the trailing-edge deflection into account, the aerodynamic function can be modified as follows

$\begin{array}{ll}cl\hfill & =C_{\text{L}\mathrm{,0}}+C_{\text{L},\alpha }\alpha +C_{\text{L},\delta e}{\delta }_e+C_{\text{L},\gamma }\gamma +C_{\text{L},h}h\hfill \\ cd\hfill & =C_{\text{D}\mathrm{,0}}+C_{\text{D},\alpha }\alpha +C_{\text{D},\delta e}{\delta }_e+C_{\text{D},\gamma }\gamma +C_{\text{D},h}h\hfill \end{array}$

12

where γ and h are the angle and length of the leading edge incision, respectively, and ${\delta }_e$ is brake deflection.

From Figure 2, the leading-edge cut coefficients $C_{\text{L},\gamma }$ and $C_{\text{L},h}$ are nonlinear with respect to the AOA. The least squares method is used here to estimate the factors. We assume that all of the influence factors concerning the AOA have the following form, where $X,Y\in \{L,D\}$

$\begin{array}{ll} \hfill & C_{X,\delta e}=C_{X,\delta e0}+C_{X,\delta e\alpha }\alpha +C_{X,\delta e{\alpha }^2}{\alpha }^2\hfill \\ \hfill & C_{X,\gamma }=C_{X,\gamma 0}+C_{X,\gamma \alpha }\alpha +C_{X,\gamma {\alpha }^2}{\alpha }^2\hfill \\ \hfill & C_{X,h}=C_{X,h0}+C_{X,h\alpha }\alpha +C_{X,h{\alpha }^2}{\alpha }^2\hfill \end{array}$

13

The aerodynamic data required is off-line. According to the least squares method, ³⁴ the unknown parameters are estimated as follows

$\theta ={\left({\Phi }_N^T{\Phi }_N\right)}^{-1}{\Phi }_N^T{Y}_N$

14

where

$\begin{array}{l}\theta =\left[{\theta }_1,{\theta }_2,{\theta }_3\right]\hfill \end{array}$

15

$\begin{array}{ll}{\theta }_1\hfill & =\left[C_{X,\delta e0},C_{X,\delta e\alpha },C_{X,\delta e{\alpha }^2}\right]\hfill \\ {\theta }_2\hfill & =\left[C_{X,\gamma 0},C_{X,\gamma \alpha },C_{X,\gamma {\alpha }^2}\right]\hfill \\ {\theta }_3\hfill & =\left[C_{X,h0},C_{X,h\alpha },C_{X,h{\alpha }^2}\right]\hfill \end{array}$

16

${\Phi }_N=\left[\begin{array}{ccc}{\Phi }_{N11}& {\Phi }_{N12}& {\Phi }_{N13}\\ {\Phi }_{N21}& {\Phi }_{N22}& {\Phi }_{N23}\\ {\Phi }_{N31}& {\Phi }_{N32}& {\Phi }_{N33}\end{array}\right]$

17

${\Phi }_{Nj1}=\left[\begin{array}{ccc}{u}_{Nj11}\left(0\right)& u_{Nj12}\left(0\right)& u_{Nj13}\left(0\right)\\ u_{Nj11}\left(0\right)& u_{Nj12}\left(1\right)& u_{Nj13}\left(1\right)\\ ⋮& ⋮& ⋮\\ u_{Nj11}\left(39\right)& u_{Nj12}\left(39\right)& u_{Nj13}\left(39\right)\end{array}\right]$

18

in which

$\{\begin{array}{l}{u}_{Nj11}\left(i\right)={\delta }_e\left(\text{floor}(i\mathrm{,10})/3\right)\\ u_{Nj12}\left(i\right)=\alpha \left(\text{mod\hspace{0.17em}}(i\mathrm{,10})\times 2+2\right){\delta }_e\left(\text{floor}(i/10)/3\right)\\ u_{Nj13}\left(i\right)={\alpha }^2\left(\text{mod}(i\mathrm{,10})\times 2+2\right){\delta }_e\left(\text{floor\hspace{0.17em}(}i/10\text{)}/3\right)\end{array}$

19

${\Phi }_{Nj2}=\left[\begin{array}{ccc}{u}_{Nj21}\left(0\right)& u_{Nj22}\left(0\right)& u_{Nj23}\left(0\right)\\ u_{Nj21}\left(0\right)& u_{Nj22}\left(1\right)& u_{Nj23}\left(1\right)\\ ⋮& ⋮& ⋮\\ u_{Nj21}\left(39\right)& u_{Nj22}\left(39\right)& u_{Nj23}\left(39\right)\end{array}\right]$

20

also, we have

$\{\begin{array}{l}{u}_{Nj21}\left(i\right)=\gamma \left(j\right)\\ u_{Nj22}\left(i\right)=\alpha \left(\text{mod\hspace{0.05em}}\left(i\mathrm{,10}\right)\times 2+2\right)\gamma \left(j\right)\\ u_{Nj23}\left(i\right)={\alpha }^2\left(\text{mod\hspace{0.05em}}\left(i\mathrm{,10}\right)\times 2+2\right)\gamma \left(j\right)\end{array}$

21

${\Phi }_{Nj3}=\left[\begin{array}{ccc}{u}_{Nj31}(0)& u_{Nj32}(0)& u_{Nj33}(0)\\ u_{Nj31}(0)& u_{Nj32}(1)& u_{Nj33}(1)\\ ⋮& ⋮& ⋮\\ u_{Nj31}(39)& u_{Nj32}(39)& u_{Nj33}(39)\end{array}\right]$

22

$\{\begin{array}{l}{u}_{Nj31}\left(i\right)=h\left(i\right)\\ u_{Nj32}\left(i\right)=\alpha \left(\text{mod\hspace{0.05em}\hspace{0.17em}}\left(i\mathrm{,10}\right)\times 2+2\right)h\left(j\right)\\ u_{Nj33}\left(i\right)={\alpha }^2\left(\text{mod\hspace{0.05em}}\left(i\mathrm{,10}\right)\times 2+2\right)h\left(j\right)\end{array}$

23

In these equations, $i\in \left(0,1,2\dots 39\right)$ , $j\in \left(1,2,3\right)$ is an aerodynamic data identifier, ${\Phi }_{Nij}$ are inputs, mod is modulo operation, and floor is integral function.

Numerical method validation

To verify the effectiveness of the numerical method, it is necessary to compare the simulation results with the available wind–tunnel test data. In addition, mesh dimensions can influence the calculations. Simulation results are meaningful only when the results do not change with the number of grids. Mesh-independent validation is to illustrate that the simulation results are independent with the mesh.

Comparison of wind–tunnel test and CFD simulations

It is helpless that model 1 used in airdrop testing has no available wind–tunnel test data. We adopt model 1 (not the model in this article) in wind–tunnel test ³⁵ as the validated model. The results are shown in Figures 7 and 8. The calculation results of CFD are consistent with the data in wind–tunnel experiment at the small AOA. With the increase of AOA, the simulation results are larger than those of wind–tunnel tests. The reason is that the cell billowing phenomenon of the canopy can result in the deterioration of aerodynamic performance, ²³ which is neglected in this article. Overall, the calculation results of CFD have a certain credibility in a reasonable range.

Figure 7.

Comparison of lift and drag coefficients between wind–tunnel test and CFD simulations. CFD: computational fluid dynamics.

Figure 8.

Comparison of lift and drag coefficients between wind–tunnel test and CFD simulations under deflection. CFD: computational fluid dynamics.

Mesh-independent validation

The simulation results based on CFD are related to the density of the mesh; in general, the results are credible only when there is enough mesh in the flow field. To guarantee the accuracy of CFD simulation, the mesh independent is performed first. The mesh density is presented in Table 4. We change the mesh density by adjusting the mesh size. It is worth noting that the specific number of points, faces, and cells may be little different from each other because of the discrepancy of molded dimensions. Table 4 presents the specific mesh parameters of model 1.

Table 4.

Detailed parameters of parafoils.

Mesh name	Number of points	Number of faces	Number of cells
Coarse	451,291	4,502,223	2,175,852
Medium	1,100,539	8,342,030	3,811,382
Fine	1,389,433	13,212,513	6,334,365

The lift and drag coefficients predicted by SA turbulence model are shown in Figure 9. To select a suitable turbulence model, four models with no deflection are also computed by the RNG $k-\epsilon$ model. The comparison result between SA and RNG is shown in Figure 10. From Figure 9, it is apparent that solutions of the medium mesh and the fine mesh match everywhere, the coarse mesh can predict correctly only within a specific range. The larger the leading-edge cut, the more sensitive the simulation results concerning the mesh. As shown in Figure 10, the parameters predicted by the SA model are larger than the RNG model, the overall error is less than 10%. Aerodynamic coefficients of model 2 are slightly different in various turbulence models. SA turbulence model and RNG $k-\epsilon$ model can give an accurate prediction under the condition of small size incision on the canopy. Through the comparison with the airdrop testing, it is founded that the parameters predicted by the SA model are more suitable than the RNG $k-\epsilon$ model; therefore, the SA model is used for computing the parameters of four models with deflection. Because the medium mesh has a good consistency with the fine mesh, in all the subsequent simulations, the medium mesh is used to save time.

Figure 9.

Mesh-independent verification: (a) model 1, (b) model 2, (c) model 3, and (d) model 4.

Figure 10.

Comparison between the SA and RNG $k-\epsilon$ turbulence models: (a) model 1, (b) model 2, (c) model 3, and (d) model 4. SA: Spalart–Allmaras.

Simulations and analysis

Effects of leading edge incision

In the traditional method, the influence of leading-edge cut was neglected. However, the incision always makes the stall angle smaller and deteriorates the aerodynamic performance of the parafoil. Figure 11 shows the cl and cd of four models with no brake deflection based on CFD. From Figure 11, the lift coefficient decreases while the drag coefficient increases under the influence of leading-edge cut. The larger the AOA is, the more significant the impact on the aerodynamic performance. Comparing four models, with the decrease of leading-edge angle and the increase of length, the canopy obtains smaller lift coefficient and larger drag coefficient at the same AOA. The weight of incision length is higher than the angle.

Figure 11.

Lift and drag coefficients based on CFD. CFD: computational fluid dynamics.

The lift-to-drag ratio is shown in Figure 12. It indicates that the leading-edge incision has a small influence on the stall angle. In addition, lift-curve slopes are basically consistent in the case of equal incision length. The lift-to-drag ratio of model 1 is 25% greater than model 2. From models 1 and 3, the length of leading-edge cut has an adverse effect on the aerodynamic performance. In addition, the impact of the leading edge length is larger than the angle. Furthermore, the smaller the leading-edge incision size is, the more slowly the lift-to-drag ratio changes, it is also beneficial to the stability of parafoil system.

Figure 12.

Lift-to-drag ratio of four models.

In order to further analyze the influence of leading-edge incision on the aerodynamic performance, the pressure distribution of four models at $\alpha =10$ is given in Figure 13. This figure offers the truth that the leading-edge cut can affect the pressure distribution of the canopy, and with the increase of incision size, the low-pressure zone on the upper surface moves forward as well as the area of the low-pressure zone becomes smaller. Figure 14 shows the bottom surface pressure distribution, once again, it is affected by the leading-edge cut, and the high-pressure region is reduced by the incision. Moreover, we take out the velocity distribution of the middle plane in the chord direction of three models. In addition, the velocity distribution is shown in Figure 15. The results are the same as the previous work, the incision is longer, and the speed loss more severe when the airflow passes through the leading edge. Reducing the angle will lead to airflow separation at the trailing edge and deteriorate the aerodynamic performance of a parafoil.

Figure 13.

Pressure distribution on upper surface: (a) model 1, (b) model 2, (c) model 3, and (d) model 4.

Figure 14.

Pressure distribution on bottom surface: (a) model 1, (b) model 2, (c) model 3, and (d) model 4.

Figure 15.

Velocity distribution on middle cross section: (a) model 1, (b) model 2, and (c) model 3.

The analysis of the leading-edge effects gives some instructions on the control of the airdrop systems. Sometimes increasing the length of incision does not influence the stability of the system, it only decreases the speed of airflow passing through the parafoil canopy. However, the angle always leads to the airflow separation at trailing edge; then, the stall angle becomes smaller, although it has little influence on lift and drag coefficients, and the range of effective AOA of PADS also becomes smaller. Therefore, if the parafoil has a relatively long incision, the control inputs can be large to implement the lift lose; while the incision angle is relatively large, the brake deflection should be small in the case of the stall of a parafoil.

Effects of trailing-edge deflection

PADSs can realize the turning motion through pulling a single steering line, and the soft landing can be achieved by pulling the left and right steering lines on the canopy trailing edge symmetrically. The pressure distribution on the canopy is not symmetrical in the case of a single steering line is pulling down. The spanwise circulation of the canopy does not satisfy the elliptic distribution. It is another drawback of traditional models. NASA wind tunnel test ³⁵ and the airdrop testing showed that during brake deflection, only a quarter from the rear of the canopy deflects, as shown in Figure 16. The whole deflection is 75°. In this article, we investigate three deflection of four models, that is, $\delta =1/3\left(25°\right)$ , $\delta =2/3\left(50°\right)$ , and $\delta =1\left(70°\right)$ , respectively. δ is the brake deflection.

Figure 16.

CFD model with the deflection of model 2. (a) Airdrop deflection, (b) δ = 1/3, (c) δ = 2/3, and (d) δ = 1. CFD: computational fluid dynamics.

To analyze the effect of different trailing edge deflections, the aerodynamic parameters of four models within three kinds of deflections are presented in Figure 17. No matter what the size of the leading-edge incision is, the brake deflection can increase both the lift and drag coefficients. At the small AOA, the lift coefficient increases with the increase in δ. However, as the AOA gradually increases, the lift coefficients reduce instead. In these four models, the drag coefficients increase all the time. Moreover, at the same brake deflection, different models have different aerodynamic characteristics.

Figure 17.

Lift and drag coefficients of four models under deflection: (a) model 1, (b) model 2, (c) model 3, and (d) model 4.

The pressure contours of three models at $\delta =1$ are presented in Figure 18. In the case of unilateral brake deflection, the pressure distribution on the canopy is not symmetrical. From Figure 18(a) and (b), with the increase of incision length, the high-pressure region on the bottom surface decreases, and the low-pressure region on the upper surface increases. Due to the length of the incision, the airflow needs to go through a longer deceleration process. Hence, the upper surface of the low-pressure zone will becomes larger. With the increase in deflection, the effects become more and more apparent. As shown in Figure 18, the pressure distribution on the upper and bottom surfaces have little difference between models 3 and 4. The small angle of the incision leads to the airflow separation and then deteriorates the aerodynamic performance.

Figure 18.

Pressure distribution when $\delta =1$ : (a) model 1, (b) model 3, (c) model 4, (d) model 1, (e) model 3, and (f) model 4.

The lift-to-drag ratio increment of different models under the same brake deflection is shown in Figure 19. In the case of small deflection, the lift-to-drag ratio of models 1 and 3 change gently. The aerodynamic characteristics of models 2 and 4 strongly vary with AOA. These results indicate that in terms of aerodynamic parameters, the negative effect of increasing length is greater than that of decreasing angle, as for the stability (it refers to the degree of change of aerodynamic parameters) of the parafoil, the effect of angle is much greater than length. As the brake deflection increases, although the increment of the lift-to-drag ratio in models 3 and 4 is larger, the stability is getting worse. Models 1 and 2 change slowly in this case. In the previous study, ³⁶ the brake deflection is converted to the variation of the airflow AOA. From this point of view, when the deflection increases, the AOA is more suitable for models 1 and 2. Overall, the smaller the leading-edge cut is, the better the aerodynamic performance of the parafoil will be.

Figure 19.

Increment of lift-to-drag ratio of four models under deflection: (a) δ = 1/3, (b) δ = 2/3, and (c) δ = 1.

Parameters identification results

Because the brake deflection coefficients and the leading-edge cut coefficients are static parameters, parameters identification can be carried out off-line. The results are obtained by the least squares method. Specific parameter values are listed in Tables 5 and 6.

Table 5.

Identification parameters related to cl.

cl parameters	Identification
$C_{\text{L},\delta a0}$	0.6094
$C_{\text{L},\delta a\alpha }$	0.038556
$C_{\text{L},\delta a{\alpha }^2}$	0.0029
$C_{\text{L},\gamma a0}$	0.5304
$C_{\text{L},\gamma a\alpha }$	0.0249
$C_{\text{L},\gamma a{\alpha }^2}$	0.00013
$C_{\text{L},ha0}$	0.8637
$C_{\text{L},ha\alpha }$	0.1251
$C_{\text{L},\gamma a{\alpha }^2}$	0.00102

Table 6.

Identification parameters related to cd.

cd parameters	Identification
$C_{\text{D},\delta a0}$	0.09819
$C_{\text{D},\delta a\alpha }$	0.026212
$C_{\text{D},\delta a{\alpha }^2}$	−0.00102
$C_{\text{D},\gamma a0}$	0.2169
$C_{\text{D},\gamma a\alpha }$	0.0293
$C_{\text{D},\gamma a{\alpha }^2}$	0.00021
$C_{\text{D},ha0}$	0.3081
$C_{\text{D},ha\alpha }$	0.0745
$C_{\text{D},\gamma a{\alpha }^2}$	−0.00082

From Tables 5 and 6, it is known that the result is the same as that of numerical calculation, the influence of incision length is greater than the angle of the incision, and the coefficient of the quadratic term is very small. To evaluate the result of parameter identification results, we compare it with the theoretical data, and the results are shown in Figure 20. The parameter identification results are consistent with the theory, the lift coefficient calculated by CFD is smaller than that calculated by theory, and the drag coefficient is larger, because CFD takes the effects of incisions and turbulence into account. With the increase in attack angle and brake deflection, the deviation between them is growing. From the analysis above, CFD has the ability to accurately predict the aerodynamic performance of the canopy in a certain range.

Figure 20.

Comparison between CFD simulations and theory results. CFD: computational fluid dynamics.

Airdrop testing

The effectiveness of the proposed approach based on CFD can be validated by airdrop testing. Our research group conducted an airdrop testing using model 1 in this article at Xu Chang city, Henan Province. The parafoil system was carried off to 400 m above sea level and released from a hot air balloon, as shown in Figure 21. The elevation of the airdrop ground was 90 m, the elevation of the releasing site was 403 m, and the altitude loss of parachute opening was 40 m. The left steering line is pulled down by $\delta =1/3$ .

Figure 21.

Airdrop testing: (a) released from balloon and (b) soft landing.

The flight data collected by the GPS module were in the form of latitude and longitude. After conversion and processing, comparison results of simulations and airdrop test are shown in Figures 22 and 23.

Figure 22.

Comparison of airdrop testing and simulation: (a) horizontal velocity and (b) vertical velocity.

Figure 23.

Horizontal trajectory in airdrop testing and simulation.

From Figure 21, the horizontal velocity and the vertical velocity model show good agreement with that of the airdrop testing, although there are still some deviations. We all know that the horizontal motion of atmosphere has a greater influence than the vertical. ³⁷ That is why the error in the horizontal direction is larger than that in the vertical direction. The wind speed measured by the anemometer is very small on the day of airdrop testing, and the trailing-edge deflection is set as constant, and this is the reason that the velocity of the vertical direction changes little. In addition, the cell billowing always appears on the canopy of the parafoil under the action of airflow. ²³ The cell billowing also affects the aerodynamic force which will be the future work. Another source of errors is the complicated noise in the flight process. From the above analysis, the aerodynamic coefficients calculation based on CFD can accurately describe the characteristics of the actual parafoil system.

Conclusions

In this article, the effects on the aerodynamic characteristics of the leading-edge incision and the trailing-edge deflection are addressed. By comparing the pressure and velocity distribution, following results are obtained.

Increasing the length and angle has adverse effect on the aerodynamic performance. The length has a great effect on the lift and drag coefficients, while the angle has an influence on the stability of the parafoil.

In the case of the unilateral brake deflection, the pressure distribution is not symmetric. The increment of the lift to drag ratio is relevant with the incision size.

The incision and deflection factors are identified using the CFD method and least squares method. By comparing with the theoretical value, the necessity for accurate calculation of aerodynamic parameters is explained.

The estimated parameters are incorporated into an eight-DOF model; simulation results of basic motions of the parafoil system and comparisons with results from airdrop testing illustrate that the proposed method in this article can accurately predict the aerodynamic parameters. This work can lay the foundation for the accurate modeling of the PADS.

The future work will include as follows:

Taking the fluid–structure interaction between canopy and airflow into account, the accurate modeling of the PADS will be realized.

Based on the accurate model, we will realize the integrated design of the parafoils, control, and trajectory optimization of the PADS.