Dynamic modelling of parafoil system based on aerodynamic coefficients identification

This study focuses on the identification of aerodynamic coefficients through airdrop data based on a known six-degree-of-freedom (6-DOF) model. Since the precise modelling of the parafoil system is a prerequisite for design of guidance, navigation and control systems, the accurate calculation of aerodynamic coefficients is vital. This paper first presents a rigid 6-DOF nonlinear model and simplifies it by linearizing the equations around the equilibrium point, which facilitates the acquisition of aerodynamic parameters next. Then, the recursive weighted least square method is applied to identify roll and yaw coefficients from airdrop test data. At last, simulations in the wind environment are implemented to analyse the dynamics of gliding and turning. The airdrop experiment also verifies the effectiveness of the dynamic model and the accuracy of identification.


R
parafoil rope length b, c canopy span and chord length d control line length that the trailing edge can be pulled down x, y, z components of position vector in the inertial coordinate system φ, θ, ψ Euler roll, pitch and yaw angles B d−s transformation matrix from the inertial to body coordinate system u, v, w components of velocity vector in the body coordinate system p, q, r components of angular velocity in the body coordinate system V c , W vector of velocity and angular velocity in the body coordinate system F W , F A gravity and aerodynamic force in the body coordinate system M A aerodynamic moment in the body coordinate system I T moment of inertial matrix of the parafoil system S W skew symmetric cross product operator of angular velocity, i.e.W× m mass of the parafoil and payload ρ density of air g acceleration of gravity S forced area of the parafoil system α angle of attack δ a asymmetric control deflection C L0 , C Lα , C Lδ a aerodynamic lift coefficients of the parafoil system C D0 , C Dα 2 , C Dδ a aerodynamic drag coefficients of the parafoil system C lφ , C lp , C lδ a aerodynamic roll coefficients of the parafoil system C m0 , C mα , C mq aerodynamic pitch coefficients of the parafoil system C nr , C nδ a aerodynamic yaw coefficients of the parafoil system δ bias constant deviation

Introduction
The parafoil system is a flexible aircraft constructed of the parafoil canopy, connecting ropes and the payload.The motion of turning or flare landing of the parafoil system can be achieved by pulling down the steering ropes, which connects on both sides of the trailing edges.Due to its stability and manoeuverability, the parafoil system is widely used in aerospace, military and civil fields [1][2][3].Therefore, to study the characteristics of the parafoil system, the first mission is to establish a precise mathematical model in accordance with the physical structure.
During the last several decades, a series of dynamics modelling studies were conducted, spanning from simplified three-degree-of-freedom (3-DOF) to 9-DOF models [4][5][6][7].In order to achieve the multi-objective global optimal homing trajectory, Zhang et al. [8] adopted the 3-DOF model involving the horizontal position and yaw angle.Prakash et al. [9] used a 4-DOF dynamic model to analyse the longitudinal motion's stability after considering the canopy and payload's velocities and pitch angle.These low DOF models may be useful to analyse the fundamental characteristics such as the position but not sufficiently accurate when there is an abrupt change in the roll of the system.Therefore, the 6-DOF model, regarding the connection of the canopy and payload as rigid and considering both the three-axis position and three Euler orientation angles, is utilized in many references [10][11][12][13].Li et al. [14] developed a linear dynamic model and modified it to a simplified model through data correlation analysis.Furthermore, an active model is established after consisted of the simplified model and the real-time model error.Sun et al. [15] combined the 6-DOF model and the Gaussian pseudo-spectrum method to achieve trajectory optimization in complex environments.Li et al. [16] designed lateral, longitudinal and speed controllers based on the 6-DOF model and used ecosystem particle swarm optimization to realize the proportionalintegral-derivative parameter tuning to accomplish trajectory tracking control.Adding two degrees, referring to the relative yaw and pitch between the canopy and payload on the basis of the 6-DOF model, Watanabe and Ochi [17] and Tan et al. [18] developed the 8-DOF model of the parafoil system with two suspending points.Lv et al. [19] improved a multibody dynamic model for a parafoil-UAV system that the suspension lines are modelled as a combination of several linear viscoelastic elements.Tao et al. [20,21] used the 8-DOF model to study the dynamic characteristics of parafoil systems in wind and rain environments.A higher DOF model compared to the 8-DOF, adding the relative roll between the canopy and the payload, is used to further realize accurate modelling, such as the 9-DOF [22][23][24].And the main 6/8/9-DOF are defined in detail as follows: • 6-DOF: The parafoil and the payload are regarded as a rigid body, which means imposing that the parafoil and the payload have the same attitude.The parafoil system includes three DOF for translational motion and three DOF for rotational motion, i.e. the threeaxis position and three Euler orientation angles.• 8-DOF: The parafoil and the payload are considered the relative motion.The model consists of six DOF of the parafoil and two DOF of the payload, which the former refers to the position and angles, and the latter refers to the relative pitch and yaw.• 9-DOF: On the basis of the 8-DOF, the relative roll motion is considered between the parafoil and the payload.
However, due to the strong nonlinear fluid-structure interaction (FSI) between the flexible canopy material and flow field, the dynamic models regarding the canopy structure as a rigid body may be lack of precision, so much literature has studied the dynamic modelling based on FSI simulations [25,26].Benney et al. [27] combined the deforming spatial domain stabilized space-time finite element algorithm with the Baldwin-Lomax turbulence model to propose a threedimensional FSI simulation model for a large parafoil.In order to solve the problem of deformation due to aerodynamics during planetary exploration flight of partially sealed parafoils, Ishida et al. [28] adopted the FSI method and used the precise code interaction coupling environment(preCICE) to analyse the instability mechanism caused by canopy deformation.And Zhu et al. [29] explored the effects of canopy inflation and trailing-edge deflections on aerodynamic performance based on FSI simulation.In addition, another crucial issue affecting the model precision is the determination of aerodynamic parameters.Owing to the high cost of airdrop experiments and the difficulty of wind tunnel testing, a method to calculate the aerodynamic coefficients based on the computational fluid dynamics (CFD) was adopted [30,31].Cao and Zhu [32] studied the effects of characteristic geometric parameters on parafoil lift and drag through simulations by CFD.Considering the effects of a leading-edge incision and trailing-edge deflection, Wu et al. [33] obtained the lift and drag coefficients by use of CFD and then applied the generalized predictive control based on a characteristic model to an 8-DOF mathematical model.Besides, Tao et al. [34] and Lv et al. [35] constructed the function relationship between the aerodynamic coefficients and the basic parameter, such as the angle of attack, deflection, etc.In view of the above work which has greater complexity and difficulty, the paper adopts the airdrop data to obtain aerodynamic coefficients, rather than the acquisition typically carried out through the simulation of CFD which cannot consider the effect of sudden uncertainty in real environment.
Once a suitable dynamic model is derived, it is particularly important to identify the aerodynamic parameters due to the difference in physical characteristics and the parafoil system's actual flight environment and then the system identification method plays an important role.Identification of nonlinear systems has become an active topic in recent years, such as Newton-Raphson [36,37] and Levenberg Marquardt [38,39].But these methods are very sensitive to initial values, parameters, etc., and may generate the problems of local optima.To avoid these problems, the recursive weighted least square (RWLS) method is adopted here.The LS method was proposed firstly by Gauss and vigorously promoted, which has been a common method used in the field of parameter identification in linear system.However, the LS method involves a large number of calculations and ignores the covariance matrix of the estimation error, which results in low identification accuracy.Then the RWLS is proposed gradually, which simplifies the difficulty of calculation and can get better effect of identification due to its strong nonlinearity of state equations.Besides, owing to the idea of recursion and weighting, the RWLS can make the most of the measured data so as to get fine prediction performance.Compared with the maximum likelihood (ML) method [40] and innovative learning machine-based method [41], which were adopted into the relatively higher DOF model that the coupling between motion parameters is strong, the RWLS method eases the work while getting great identification effect.Hence, the paper will introduce a 6-DOF model that can describe the parafoil system's basic motion and simplify it [42,43], then the RWLS method is directly used in the model to acquire aerodynamic coefficients.
In this paper, a rigid 6-DOF model is presented first.Then considering the motion characteristics when the parafoil system reaches a steady state, the detailed derivation of simplified model is elaborated.After that, the airdrop test data is adopted and processed to identify the aerodynamic coefficients using RWLS.Furthermore, numerical simulations are carried out to illustrate the characteristics of the parafoil system and compared with airdrop test.The contributions of this study are listed as follows: (1) A 6-DOF model of the parafoil system, which can provide a reference for the modelling, is established.The model can be used to analyse the motion characteristics and accumulate experience for parafoil airdrop engineering applications.(2) A method to simplify the 6-DOF is described in detail, which can overcome the nonlinearity and make further study easier.(3) The RWLS method is elaborated comprehensively and adopted into the identification of aerodynamic coefficients based on airdrop data, which can provide an idea for the precise dynamic modelling.
The remainder of this paper is organized as follows.In Section 2, the conventional 6-DOF dynamic model is established, and a simplifying assumption is proposed, which builds a linear function between the variables.In Section 3, given the features of the simplified model, the RWLS method is utilized based on airdrop data.In Section 4, the simulation results are analysed and compared with the experiment data.In Section 5, the conclusions obtained are presented.

Geometric parameters
During the flight, the air chambers of the canopy are filled with gases, and the parafoil system will be in a completely unfolded and stable state that can complete actions next, such as turning.The front and side views are simplified, as shown in Figure 1.

Coordinate systems
In order to facilitate the analysis, two main coordinate systems are established, which are the inertial coordinate system O d x d y d z d and the body coordinate system O s x s y s z s , as shown in Figure 2. The origin O d of inertial coordinate system is chosen on the surface of the earth and the positive direction of the z d -axis is taken downward.The plane of O d x d y d is horizontal, where the x d -axis is chosen in the direction of north.The origin O s of the body coordinate system is chosen at the centre of gravity (CG) of the parafoil system.The positive direction of the z s -axis is pointed to the CG of the payload.The x s -axis is in the symmetry plane of the parafoil and is perpendicular to z s -axis.The y s -axis is chosen so that the system O s x s y s z s is a right-hand coordinate system.
The conversion from the O d x d y d z d to the O s x s y s z s can be achieved by Euler orientation angles and the transformation matrix with the sequence of "ZYX" i.e. rotating around the z-axis first, then y-axis and x-axis last, is expressed as where for arbitrary angle α, s α = sin(α), c α = cos(α) and t α = tan(α).
Then the relationship of velocity vector between the above systems can be expressed as where V c = uvw T .The relationship of angular velocity vector can be expressed as where W = p q r T .

Dynamic model
The parafoil system has translational and rotational motions when the aerodynamic force and moment act on it.For simplicity, the canopy and payload are regarded as a whole with the same velocity and Euler orientation angles.Before modelling, the following modelling assumptions [44] are used: (1) Apparent mass [45,46] and inertial effects are negligible in near steady-state conditions;  (2) Moments of aerodynamic forces on the centre of mass (CM) are much smaller than those of the aerodynamic rolling, pitching and yawing; (3) Payload drag is proportional to the surface area and therefore much smaller than the parafoil drag; (4) The contribution of the parafoil to the velocity due to different locations of the centre of parafoil pressure and CM is negligible (i.e. the velocity of the parafoil is considered the same as that of the CM).
Then the dynamic model in the body coordinate system is derived as follows: where the F W , F A , M A , S W , I T are expressed as And to facilitate the next analysis, Equations ( 4) and ( 5) can be rewritten as

Simplified model
Based on the established dynamic model, determining the values of the aerodynamic coefficients in Equations (7) and ( 8) is the next step so as to analyse the motion characteristics.However, the model shows problems of high nonlinearity and strong coupling, rendering more difficult estimating the parameters values.This work makes use of a model linearization step.When landing from a high altitude to the ground freely or under the condition of unilateral deflection, the parafoil system will reach a relatively steady state, especially the pitch angle.At the period, observation of the parafoil system shows that the velocity V c changes slightly, almost a constant [47], which implies that all forces are balanced, thus eliminating Equation (11).Then the linear equations are derived by calculating Equation ( 12) around the equilibrium point while eliminating the coupling terms.Last, a constant deviation δ bias is introduced to reduce the influence of the linearization process.
In this study, the parafoil system will take-off on the ground and be controlled to an initial position by humans.When the parafoil system is switched to the airdrop phase, there exists a buffer period, during which the data are not sufficiently efficient if used in the work.Besides, owing to the condition that the velocity of the parafoil system is constant, it is necessary to ensure that the drop is steady without sudden or frequent turning as much as possible.Considering both conditions, the state of the parafoil system is selected as the equilibrium point.Assume the values at the equilibrium point as follows: The simplified model is expressed as where

Analysis of simplified model
It can be seen from Equations ( 15)-( 21) that the relationship between the input and output variables are linear.Substituting Equations ( 16)-( 21) into Equation ( 15) extracts the common coefficients and Equation (22).
where ṗsim and ṙsim are output variables and φ sim , p sim , r sim and δ asim are input variables, which can be obtained from airdrop experiments; C lφ , C lp , C lδ a , C nr , C nδ a and δ bias are coefficients to be identified.

Analysis and processing of airdrop test data
In order to make the values identified more practical to the actual environment, the paper adopts the airdrop test data.Our group has developed a parafoil system in which the payload is a small powered vehicle as shown in Figure 3, and the details of embedded control can refer to in the description of Ref. [48] in Section 4.
An experiment was carried out in Weifang City, Shandong Province, China in November 2020, in which the parafoil system was lifted up by thrust and released.Before using the airdrop data for identification, three steps should be done.First, because the airdrop data have a great impact on the identification effect, the paper simulates the flight trajectory in near 10 experiments and determines one in accordance with the conditions described in Section 2.4.Then, since the parafoil system used in the experiment is powered, the work intercepts the phase without the control of power after visualizing the input change curve.Finally, owing to the measurement noise of sensors, Kalman filtering is performed on airdrop data in order to be further analysed.For example, the comparison of roll angle is shown in Figure 4.
Determining one feasible test and processing the data in MATLAB, the flight trajectory is shown in Figure 5.It can be seen that the parafoil system descended from a height of nearly 124 m, turning first and then gliding.The dynamic change curves of the input and output variables in Section 3.1 are shown in Figure 6.We can know that the flap deflection changes during the flight, and the value is between [−1,0].The roll angle fluctuates between [2,18] degrees.Angular velocities and accelerations are acquired by solving Equation (3) in MATALB on the basis of airdrop data, and their values are shown in Figure 6(c,d).

Identification of aerodynamic coefficients
In order to identify the aerodynamic parameters in Equation ( 22), the objective function is established as follows:  where z i is the model output and can be expressed as x is the parameter to be identified and is expressed as n i is the zero average measurement noise.The measurement noise is introduced to account for the effect of unknown disturbances in the measured signals.
h i is the state matrix, which maps the input variables to output variables, and the expression is: Due to its superiority described in Section 1, the RWLS method is adopted here.The basic thought of the RWLS is that the current estimated value is modified from the last estimated value and weighted, where the weighting factor is Q.So, the estimated value based on the previous i sets measurement data is: When a new set of input and output data is added to the former measurement data, the model equation is: where Then according to the optimal thought, the estimated value of the sequence i + 1 can be expressed as Let Then through a series of derivations, the recurrence formulas can be expressed as P is initialized as a 6×6 diagonal matrix, with each diagonal entry set to 0.05.The 2×2 matrix Q consists of the diagonal entries Q 1,1 and Q 2,2 which take values 0.00475 and 0.0005, respectively [11].And the x is initialized as −0.1, −0.16, 0.002, −0.14, −0.02, 0. Import the input and output data into MATLAB and then get the values of those six aerodynamic coefficients, shown in Table 1.

Validation
The model validation is conducted to analyse the motion characteristics.The geometric parameters of the parafoil system are listed in    first.Based on the above properties, the work analyses the difference and adjust the values correspondingly.
The specific values are shown in Table 3.

Simulation of motion performances
Simulation is an efficient and essential way to study the characteristics of the plant before practical experienment.Ward and Costello [49] presented the development of a control law to implement glide slope control on an autonomous airdrop system, and demonstrated an improvement in landing accuracy by simulation and flight test.Tao et al. [50] designed the guidance law based on the hybrid approach that combines the crosstrack error and the line of sight, which was simulated in a semi-physical simulation platform.Chen et al. [51] developed the guidance by modern multi-objective evolutionary algorithms without assuming any predefined trajectories, and adopted the software developed with Kotlin, Java and JavaFX to implement solution.In addition, Zhang et al. [10] carried out the simulation of motion during the flight, and Farì and Grande [52] implemented the dynamic models and G&C algorithms in the MATLAB/Simulink and Dymola environments.
In this paper, a MATLAB code has been utilized to simulate the 6-DOF model, and the flowchart of the algorithm is presented in Figure 7.

Gliding
The parafoil system glides with no manipulation.Then the corresponding trajectory, speed and angle characteristics are shown in Figure 8.
It can be seen that the trajectory in the horizontal trajectory is a straight line, and thus the velocity in the y-direction is 0. The velocity in the x-direction is 6.2 m/s and in the z-direction is 1.92 m/s, with the gliding ratio, defined as the ratio between forward and descending distances, is 3.25.The pitch angle stabilizes at 18.6 degrees, and there is no roll and yaw motion.

Turning under unilateral deflection
The parafoil system will turn if pulling down one of the steering ropes.The working conditions are set as follows: After 50 s gliding, the left flap is set to 20%, and the corresponding motion characteristics are shown in Figure 9.
It can be seen that the parafoil system is in the state of steady flight in the first 50 s, then it begins to turn left, and the trajectory in the horizontal plane is a circle with a radius of 102 m.And if pulling down the right flap by 20%, the parafoil system will turn right, shown in Figure 6(a), with the same radius.The vertical velocity stabilizes at 1.99 m/s, almost unchanged.But the velocities in the x and y directions change abruptly, showing a sinusoidal wave.When the velocity V x reaches the maximum, the V y is zero, and vice versa.The roll and pitch angles happen to change, from 0 to 0.57 degrees and −18.6 degrees to −19.17 degrees, respectively.The yaw angle increases almost as a linear function of the time when the turning motion begins.

Wind responses
Considering the actual environment, due to the unique nature of the parafoil material whose mean density is close to the air, the wind will affect the flight trajectory to some extent.In this section, the wind speed is given relative to the inertial coordinate system and let it be V w = 0 3 0 T m/s, which is along the direction of y d .Then the trajectories of gliding and turning are shown in Figure 10.
When gliding in the wind, the V w is added to the model at the time of 50 s, then the vehicle moves no longer just along the x-axis but towards the y-direction, which is consistent with the direction of wind speed, as shown in Figure 10(a).When turning in the wind, the condition is as follows: After 50 s, the left steering rope is pulled down by 20%, and after 80 s, when the vehicle reaches a steady-state, the wind V w is added to the simulation environment.The trajectory in the horizontal plane is an upward spiral curve, just shown in Figure 10(b).Both the motions last 200 s.

Airdrop experiment
In order to verify the validity of the dynamic model and the accuracy of the aerodynamic coefficients identified, it is meaningful to compare the airdrop data with the model simulation data.As analysed in Section 4.2, the parafoil system is sensitive to the wind.Thus, the paper assumes that the wind is constant and added into the simulation environment.The result is shown in Figure 11.It can be seen that the trajectory of simulation almost overlaps with the flight track of the actual parafoil system, while they still exhibit small numerical errors.It may be caused by the following reasons: (1) The dynamic model established is 6-DOF and very simple without considering potential factors that affect the parafoil system's motion, such as the apparent mass and the payload drag, etc. (2) The airdrop data collected by sensors exhibits some measurement noise.After processing, the noise is amplified and may lead to the inaccuracy of parameters identification.(3) The wind speed is assumed to be constant in the simulation environment.However, the wind in airdrop experiment is time-varying, which may cause errors in comparison to simulation and airdrop data.

Conclusion
In this paper, a dynamic 6-DOF rigid model is constructed, and the RWLS is used to identify the aerodynamic coefficients with flight data after simplifying the established model.Then the following conclusions are mainly obtained: (1) The 6-DOF model is simplified through linearization to solve the problems of strong nonlinearity and high coupling, which lays the foundation for the subsequent identification of aerodynamic coefficients.
(2) Using experimental airdrop data, the aerodynamic coefficients are identified by RWLS, which the result can quickly converge to be stable and the numerical values are within the empirical range so as to get preliminary verification of accuracy.(3) By simulating, the paper analyses the motion characteristics of the parafoil system when gliding and turning with and without wind, which conforms to the theory.Finally, the airdrop experiment is conducted to verify the effectiveness of the dynamic model and the accuracy of the values of roll and yaw coefficients.
However, the model needs further improvement.The apparent mass, the difference in characteristics between the canopy and payload, and the real-time wind field should be considered in the model in the next step so as to build a foundation for future research.

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

Figure 1 .
Figure 1.The schematics of the parafoil canopy.(a) The front view.(b) The side view.

Figure 4 .
Figure 4. Comparison with or not Kalman filtering.

Figure 5 .
Figure 5. Flight path.(a) Trajectory in the altitude.(b) Trajectory in the horizontal plane.

Figure 7 .
Figure 7.The flowchart of the algorithm.

Figure 10 .
Figure 10.Trajectory in the horizontal plane in wind.(a) Gliding in wind.(b) Turning in wind.

Figure 11 .
Figure 11.Trajectories in the horizontal plane.

Table 1 .
The estimated model coefficients.

Table 2 .
Geometric parameters of the parafoil system.

Table 2 .
And the aerodynamic coefficient values inTable 1 are substituted into the dynamic model derived in Section 2.3.As for the determination of the remaining

Table 3 .
Aerodynamic parameters of the parafoil system.