# Parachute-Payload System Flight Dynamics and Trajectory Simulation

Authors: Giorgio Guglieri
Journal: International Journal of Aerospace Engineering
 (2012)

---

## Abstract

The work traces a general procedure for the design of a flight simulation tool still representative of the major flight physics of a parachute-payload system along decelerated trajectories. An example of limited complexity simulation models for a payload decelerated by one or more parachutes is given, including details and implementation features usually omitted as the focus of the research in this field is typically on the investigation of mission design issues, rather than addressing general implementation guidelines for the development of a reconfigurable simulation tool. The dynamics of the system are modeled through a simple multibody model that represents the expected behavior of an entry vehicle during the terminal deceleration phase. The simulators are designed according to a comprehensive vision that enforces the simplification of the coupling mechanism between the payload and the parachute, with an adequate level of physical insight still available. The results presented for a realistic case study define the sensitivity of the simulation outputs to the functional complexity of the mathematical model. Far from being an absolute address for the software designer, this paper tries to contribute to the area of interest with some technical considerations and clarifications.


---

## Body

## 1. Introduction

The purpose of a parachute is to decelerate and provide stability to a payload in flight. The aerodynamic and stability characteristics of the parachute system are governed by the geometry of the parachute as such careful consideration is paid to this in the design process. The effects of deployment and opening force are critical in the safe operation of the parachute and the integrity of the payload. The opening characteristics also feature heavily in the selection of geometry and other parameters in the design process.Parachutes for aerospace applications [1–4] are in general symmetric about the canopy axis. This axis passes through the center of the canopy and the confluence point of the suspension lines. The canopy is the cloth surface that inflates to provide the desired lift, drag, and stability. The suspension lines transmit the retarding force from the canopy to the payload either directly or through a riser attached below the confluence point of the suspension lines. The deceleration force may be distributed on the payload over more than one mechanical joint linked to the riser by a set of short hardly extensible strips (bridles).There are a number of different kinds of parachutes that have been designed for various applications. The different applications parachutes are typically used for pilot, drogue, deceleration, descent, extraction, supersonic drogue and stabilization, flight termination, and landing.The dynamics of parachutes are complex and difficult to model accurately. During both the inflation process and the terminal descent stage, the dynamics of a parachute are governed by a coupling between the structural dynamics of the parachute system and the surrounding fluid flow. Both of these dynamic systems must be addressed as a coupled system to gain a proper representation of the dynamic system as a whole.When the parachute is in a steady state, the air flowing around the decelerator will separate at some location on the canopy. The shedding of the vortices from the canopy can affect the stability and cause a periodic motion of both parachute and payload. The wake from a porous parachute consists of air that flowed around the canopy and air that flowed through the canopy. A payload body in the speed range of parachute usage sheds a very turbulent wake. Part of the flow that is entering the parachute is therefore of a disturbed nature and should be considered regarding the aerodynamic performance of the parachute. For many types of parachutes, this change in oncoming airflow can be quite significant during the time required for the parachute to inflate. The implication of a rapid deceleration is that second-order effects are likely to be present.To summarize, calculation of parachute deployment, inflation, and deceleration requires the numerical solution to the equations of motion for a viscous, turbulent, separated airflow. The parachute is also a flexible body having dynamic behavior coupled with the behavior of the flow, which passes through and around it. From the above description it is obvious that a full-time dependent solution of this system is far from being easily feasible. To make a mathematical model that is feasible, simplifications must be made, as long as the model can be validated satisfactorily by experiment or by comparison with reference data.The overall behavior of parachutes is related to various parameters: added masses, filling time, parachute shape (inflated canopy elongation), porosity, suspension line length, reefing, clustering, snatch loads at deployment, and aero-mechanical and inflation instability. In the past, most of these effects could be generally modeled in an imprecise way by simulation tools. A comprehensive computational technique is presented in [5–7] for carrying out three-dimensional simulations of parachute fluid-structure interactions, and this technique is applied to simulations of airdrop performance and control phenomena in terminal descent. The technique uses a stabilized space-time formulation of the time-dependent, three-dimensional Navier-Stokes equations of incompressible flows for the fluid dynamics part. A finite-element formulation derived from the principle of virtual work is used for the parachute structural dynamics. The parachute is represented as a cable-membrane tension structure. Coupling of the fluid dynamics with the structural dynamics is implemented over the fluid-structure interface, which is the parachute canopy surface.According to the different missions, several types of payloads have been used in combination with aerodynamic decelerators: paratroops, equipment, hardware, materiel, weapons, missiles, aircraft, unmanned aerial vehicles, aerospace lifting, and nonlifting spacecraft. The present analysis is focused on aerospace applications for planetary and atmospheric entry vehicles, where the payload is typically a blunt body. The purpose of a blunt body is primarily to provide a large source of drag to facilitate a deceleration. Applications of blunt bodies can be seen with both manned reentry and planetary exploration missions. An outline of aerodynamic decelerators for robotic planetary exploration missions is given in [8]. In the development programs of such blunt bodies and also in subsequent studies, it has been shown that they may be dynamically unstable in all or part of the sub-, trans-, super-, and in some cases hypersonic speed regimes.Blunt bodies for which there exist examples can be classified into two categories, large angle cones and capsules [3]. An example of large angle cone that has performed actual missions is the Viking probe. It should be noted that in the two Viking missions confirmed the dynamic instability upon entry into the Martian atmosphere. The examples of capsules are well-known ones due to the intense space activities from the former USSR and the USA. The capsules that were developed for manned flights in this period were Mercury, Gemini, Apollo and Soyuz, which are still considered a reference for performing missions to this day. As was the case with the large angle cones, there were found to be speed ranges over which the capsules were unstable [9].The flow field that is associated with the capsule configuration is highly complex. For almost the entire speed range that the capsules operate over, the flow remains attached on the forward face. At the point of maximum diameter, the flow is accelerated such that the boundary layer rapidly grows to the point of separation. After the maximum diameter, the flow then remains separated and turbulent. This flow is unstable and coupled with an unsteady near wake the dynamic instability associated with capsules is produced [11].As dynamic instability exists for blunt bodies over various speed ranges, to complete missions successfully, there is then a requirement for some kind of accurate stability assessment with a potential impact on both stability augmentation (if any reaction control system is implemented) and mission design (parachute deployment sequence). This issue applies for probes and capsules as thrusters and parachutes are the unique available sources of additional damping.The measurement of stability derivatives for probes and capsules was a concern of designers since the origin of space flight [9], as analytical methods did not provide (at least in the past) adequate estimation of these relevant parameters. The wind tunnel experimental techniques considered are free oscillation, forced oscillation, free flight tests, and even ballistic range experiments. It should be noted that wind tunnel experiments are not without limitations, so there remains a strong desire to make available numerical data, for the purposes of both finding solutions and comparing with experimental results.


## 2. Background on Parachute-Payload Modeling and Simulation

Several studies analyze the descent and landing trajectories of parachute-payload systems. Generally, these analyses consist of performing a simulation (typically a 3–6 DOFs rigid body model is adopted) of the atmospheric entry phase to predict deceleration peaks, descent attitudes, and terminal conditions. In addition, a stochastic dispersion analysis (Monte Carlo or similar) is usually performed to assess the impact of off-nominal conditions that may arise to determine the robustness of the mission design.A two-dimensional parachute model is presented in [12] to compute the various characteristics of the steady descent of a parachute system. A three degree-of-freedom analysis is presented and validated in [13] giving the longitudinal motion of a typical vehicle during the recovery phase. The parachute and the payload are supposed to be rigid and interconnected by an elastic riser. Aerodynamic loads acting on the two subsystems are considered. Computer results showed good agreement with test results in terms of oscillation amplitude and frequency, riser force, and parachute wrapup about the vehicle for the simulation of a pad-abort situation. The three-dimensional motion of a freely descending parachute is studied in [14] with a five degree-of-freedom analysis (roll motion is neglected). Exact expressions are given for the longitudinal and lateral small disturbance stability of the gliding motion of parachutes. The analysis confirms that large longitudinal disturbance of most parachutes will result in a large pitching motion, whereas a large lateral disturbance will usually cause a large angle vertical coning motion (coning mode). The longitudinal mode damps out very quickly in the stable case. The three-dimensional motion of a nonrigid parachute and payload system is studied in [15]. Both the parachute and the payload are assumed to have five degrees of freedom (roll about axis of symmetry is again neglected). They are coupled together by a fixed-length connector. The general nonlinear equations are linearized using small perturbation theory. The evaluation of the stability of an unstable payload decelerated by a parachute is performed. The authors observed that increasing riser length and parachute weight promotes system instability. A nine degree-of-freedom computer program was developed in [16] for the simulation of the trajectory and the dynamic behavior of a rotating parachute system. An accurate mathematical model of the joint between the load and the parachute was found to be necessary to predict the dynamic behavior of a rotating decelerated system. A computer model based on a six degree-of-freedom analysis is described in [17] and compared with drop test data. The payload is rigidly connected, the aerodynamic forces on canopy and payload are determined by the instantaneous angle of attack of the impinging airflow, the apparent masses are constant, but they depend on the direction of the acceleration. Full nonlinear equations of motion for the axisymmetric parachute have been obtained in [18]. In particular, the correct form of the added mass tensor for a rigid axisymmetric parachute in ideal flow has been implemented in a six degree-of-freedom computer model [19], and the results indicate that added mass effects are significant. In particular, the component of added mass along the axis of symmetry has a strong effect on parachute dynamic stability. However, design and testing experience shows that dynamic stability of the parachute is a second order design problem [3] for high-performance decelerators that usually have both high static and dynamic stability due to the porosity of the canopy. Many works, including recent models developed and validated by means of drop tests, consider the bridles as rigid elements [20, 21].


## 3. The Mathematical Model of the Payload-Parachute System

A first mathematical representation of a parachute-payload system (designated as SM1) is defined according to a very comprehensive approach that enforces the simplification of the coupling mechanism between the payload and the parachute, with an adequate level of physical insight still available in terms of sensitivity of system performance metrics to design parametric changes. In modeling a payload-parachute system, two bodies and one device have been used. They are the parachute canopy and the payload, which are then connected by a single riser. In reality the parachute canopy is connected to the payload by suspension lines, a riser section, and then a set of bridles. In the simulation only the riser and the bridles are modeled in terms of a dynamic response (suspension lines are neglected). However, in terms of connection points only the riser has been modeled. The riser is assumed to have flexible connections at both ends and provides stiffness above a certain threshold distance and zero stiffness below this distance (slack conditions).The payload is considered as a rigid body with six degrees of freedom. The forces and the moments that act on the body are provided by its aerodynamics, by the weight, by the inertial actions, and, finally, by the force applied by the riser in the suspension point. The parachute canopy acts as a rigid body and, by the action of aerodynamic drag, strains the riser which then transmits a force to the payload. Figure1 shows the general philosophy used in modeling the system.Figure 1
Flow chart of the simulation model.In terms of a dynamic response, the riser is modeled as having linear stiffness and damping, where the force at any given point in time is given by(1)FR=klrε+cε̇.
Note that k is stiffness, lr is riser length, ε and ε̇ are strain and strain rate, respectively, and c is the damping coefficient. The above equation is implemented in the code to calculate the force present in the riser at any point. The equivalent stiffness of the bridles and riser is given by the following equation:(2)k=nb⋅krkbkr+nb⋅kb,
where nb is the number of bridles, kr is riser stiffness and kb is single bridle stiffness. After calculation of the equivalent stiffness, the damping coefficient c can be calculated by the following formulation (mP is the parachute mass and ξ is the damping ratio):(3)c=2ξmPkmP.
The values for strain and strain rate are computed by the difference in displacement and velocity of each end of the riser, that is, as the displacement and the velocity of the payload and that of the parachute canopy. Also note that where this calculation occurs a check is in place for the condition where the riser is slack:(4)ε=sc-sp-l0lr,ε̇=Vc-Vp,
where ε and ε̇ are, respectively, strain and strain rate, sc and sp are, respectively, payload and parachute displacement, lr is the riser length, and Vc and Vp are, respectively, the payload and parachute velocity. From these two calculations, the force in the riser can be derived by application of (1).The performance of the parachute has been modeled using an approximation for the added massma. When a parachute passes through the air, it drags some of the surrounding with it this means that this air must be accelerated from rest up to the speed of the parachute. The parachute will also accelerate some of the air immediately in front and behind. The energy that is required to facilitate this acceleration is taken from the kinetic energy of the parachute-payload system. These fluid inertia effects are taken into account in the equation of motion. The idea is that the mass has a component added to take into account the extra mass associated with the accelerated air. This can be seen below, where the change in momentum of the parachute system is equal to the sum of the drag force and gravitational force component:(5)ddt[M*V]=-12ρV2CDSP+mPsin⁡γ,(6)M*=ma+mP.
After combining (5) and (6) and then rearranging, the following expression for system acceleration (or deceleration) can be used:(7)(mP+ma)dVdt=-12ρV2CDSP-Vdmadt+mPg⋅sin⁡γ.
However, for (7) to be useful, an expression or knowledge of ma(t) must be available. Shown below is an expression for added mass, similar to that for apparent mass, that has been used widely:(8)ma=kaρ43πR3.
The value of ka, a constant, can be found in various ways [1, 3]. Klimas calculated values of ka for porous hemispherical ribbon parachutes. In the present case, added mass is calculated by using a slight variation. In the previously presented method, the volume of a hemisphere is used, whereas in the model used in the simulation the volume of an ellipsoid is used. Hence, added mass has been calculated using the following expression:(9)ma=nPkaρ43πR2h,
where nP is the number of parachutes, ka is the added mass coefficient, R is the inflated radius, and h is the inflated height. The integer nP for the number of parachutes is necessary because a cluster of multiple parachutes is being modeled as one parachute of equivalent size. The added mass remains constant throughout inflation, so the values for parachute radius and height are constant, whereas in reality they are of course changing quite significantly through the inflation process. The added mass coefficient is given by:(10)ka=1.068(1-1.465p-0.25975p2+1.2626p3),
where p is the parachute porosity (see [3, 22] for further details concerning the effect of porosity on added masses). This coefficient remains constant through inflation. The porosity for both drogue and main parachutes is set at a value of 20% (p=0.2).The significance of fluid inertial effects is relevant. As an example, opening times and loads increase with altitude. The reason for this can be seen when it is considered that at increasing altitude the true airspeed will increase for constant indicated airspeed, as such the inertial effects of the air are greater. After comparison with flight test data, the added mass approximation tends to overestimate the deceleration that will occur. From a design point of view, this means that results obtained from the added mass approximation can be viewed as conservative. So while the concept of added mass is only an approximation, the effects of fluid inertia in transient parachute aerodynamics are significant enough to be included in modeling. At present this method of accounting for the inertial effects of air is the most practical tool available to be used in the transient phase of parachute operation.The parachute as previously stated is modeled as a separate body (rigid canopy). The deceleration of this body is therefore calculated by extension of (7) including the riser tension force:(11)dVdt=V̇P=Fr-DP-mPgsin⁡γ(mP+ma),
where a is acceleration, Fr is the riser force, DP is the parachute drag, mP is the mass of the parachute, ma is the added mass, and γ is the flight path angle. The drag force provided by the parachute is calculated using the following simple expression:(12)DP=12ρVP2CDSP⋅η,
where ρ is air density, VP is parachute velocity, CD is the canopy drag coefficient, SP is the canopy projected area, and η is an efficiency factor used to account for payload wake effects (η=0.7÷0.75 for the drogue parachute and η=1 for the main parachute). The product CDSP (drag area) of a parachute depends on its type, inflated shape, size, Reynolds’ number, Mach’s number, Froude’s number, material elasticity, and porosity.In (11) the third term on the numerator is the gravitation force component. The denominator is then the addition of the mass of the parachute and the added air mass.Due to the complexity involved with the filling process, filling time cannot be calculated by a purely analytical method. The alternative once again lies with an empirical approach that is validated by flight test. Knacke defined the following empirical relation for filling time:(13)tf=8D0VS0.9,
where D0 is nominal diameter and VS is velocity at line stretch. It should be noted here that (13) is dimensionally incorrect. This means that the actual physics of the system are not properly represented. While (13) has been shown to provide reasonable answers, care should always be taken in using such an equation so that it is not being used outside the intended range of conditions.An underlying assumption used in the modeling of the parachute canopy is that it remains aligned with the velocity vector at all times. This corresponds to neglecting the lift and moment coefficients of the parachute.The second major assumption made used in the parachute canopy model is that the added mass remains constant throughout inflation and parachute operation.Whilst in the model the value of added mass is kept constant at the value for the fully inflated parachute (conservative approach for the estimation of deceleration peak), the value for drag area is ramped up as the parachute inflates. It is this feature of the model that provides a simulation of the inflation process. Inflation modeling is done by ramping the value of the drag area,CDSP, this being the multiplication of the drag coefficient and the projected area. The effective drag area is time-scheduled according to the planned reefing stages and can be described by the following expressions:(14)(CDSP)t=(CDSP)0⋅[∑k=1i-1τk+τi(t-titfi)ni],fort∈[ti÷tfi],(CDSP)t=(CDSP)0⋅∑k=1iτk,fort∈[tfi÷ti+1],
where reefing starting time is ti, duration tfi, and ratio τk. The power for the growth, ni, has been set at 2 for the initial stages and at 2.5 for the final stage, for both drogue and main parachutes (the exponents are selected to fit the available flight test data). Selection of these values is justified by a better prediction of peak loads. This value of CDSP is then substituted into the equation for parachute drag.The final assumption made is with the clustered parachutes. In the simulation one parachute of equivalent size has been used, with a cluster efficiency factor applied.In Figure2, the forces acting on both the parachute and payload can be seen (with the exception of the body aerodynamic forces and moments).Figure 2
Layout of the suspension system (parachute-payload configuration).The forces and moments acting on the payload are calculated as the relevant contributions are added. The equations are(15)X=12ρVS(VCX+CXqqd)-FRcos⁡αcos⁡β,Y=12ρVS(VCY+CYrrd)-FRsin⁡β,Z=12ρVS(VCZ+CZqqd)-FRsin⁡αcos⁡β,L=12ρVSd(VCl+Clppd)-FRsin⁡β(zcg-zp)d-FRsin⁡αcos⁡β(ycg-yp)d,M=12ρVSd(VCM+CMqqd)+FRcos⁡αcos⁡β(zcg-zp)d+FRsin⁡αcos⁡β(xcg-xp)d,N=12ρVSd(VCN+CNrrd)-FRsin⁡β(xcg-xp)d-FRcos⁡αcos⁡β(ycg-yp)d,
where X, Y, and Z are the body forces, FR is the riser force (tension), L, M, and N are the body moments, V is velocity, d is payload diameter, S is reference surface area, p is roll rate, q is pitch rate, r is yaw rate, α  is angle of attack, β  is angle of sideslip, and xp, yp, and zp are the parachute-payload attachment coordinates. All coordinates are normalized with a reference length.In the above-presented model, there are a number of featured that need to be discussed. The aerodynamic data for the payload—which the simulation uses as input—are given in the longitudinal plane only. In order to calculate the correct coefficient, the center of gravity must be shifted accordingly. This is done using the following equations:(16)CM=CM-CZ(xcg-xref)-CX(zcg-zref),CMq=CMq-CZq(xcg-xref)-CXq(zcg-zref),
where CM is the moment coefficient, CZ is the force coefficient in the Z direction, CX is the force coefficient in the X direction, CMq is the moment coefficient due to pitch rate, CZq is the force coefficient in the Z direction due to pitch rate, CXq is the force coefficient in the X direction due to pitch rate, xcg and zcg are the x and z locations of the payload center of gravity, and finally xref and zref are the x and z locations of the reference data center of gravity. The coordinates are normalized with a reference length, that is, for axis-symmetric bodies the payload diameter d.The lateral-directional coefficients are obtained assuming the geometric symmetry of the payload. In taking this approach, the aerodynamic response to angle of attack and angle of sideslip have been decoupled. The implications of this decoupling are that there are no aerodynamic roll moments (L=0), no sideslip effects in the longitudinal plane, and finally no angle of attack effects in the lateral plane. It should also be noted that, since the body is axis-symmetric, total angle of attack would need to be considered for aerodynamic coefficients interpolation, as implemented in [23].The body axis system has been adopted, and this system can be seen in Figure3. The system of differential equations that describes the parachute-payload system can be broken into sections, six rigid body equations of motion, four quaternion equations, three position equations, and the velocity and displacement of the parachute.Figure 3
The reference frame and the aerodynamic angles for the payload.The system written in residual form can be seen below:(17)r1=u̇-((X+mgCBE(1,3))m-qw+rv),r2=v̇-((Y+mgCBE(2,3))m-ru+pw),r3=ẇ-((Z+mgCBE(3,3))m+qu-pv),r4=IXXṗ-IXYq̇-IXZṙ-IXZpq+(IZZ-IYY)qr+IXYpr+(r2-q2)IYZ-L,r5=IYYq̇-IXYṗ-IYZṙ-IYZpq+(IXX-IZZ)pr-IXYqr+(p2-r2)IXZ-M,r6=IZZṙ-IXZṗ-IYZq̇+IXZqr+(IYY-IXX)pq-IYZpr+(q2-p2)IXY-N,r7=ė0+12(e1p+e2q+e3r),r8=ė1-12(e0p-e3q+e2r),r9=ė2-12(e3p+e0q-e1r),r10=ė3+12(e2p-e1q-e0r),r11=ẋ-CBE(1,1)u-CBE(1,2)v-CBE(1,3)w,r12=ẏ-CBE(2,1)u-CBE(2,2)v-CBE(2,3)w,r13=ḣ+CBE(3,1)u+CBE(3,2)v+CBE(3,3)w,r14=ṡp-Vpr15=V̇p-(Fr-Dp-mpgsin⁡γ)(mp+ma),
where u,v, and w are the body axis velocities in the X,Y, and Z directions, respectively, p,q, and r are the roll rate, pitch rate, and yaw rate in body axis, e0 through to e3 are the quaternion values for orientation, γ is the flight path angle, x,y,and h give the position in the Earth fixed axis, sp is the absolute displacement of the parachute canopy, and finally Vp is the parachute canopy total velocity.The magnitude of the velocity vector, that is, the airspeed, is found from the three body-axis velocities:(18)V=u2+v2+w2.
The derivative of the velocity vector is then calculated by(19)V̇=(uu̇+vv̇+wẇ)V.
The angle of attack is calculated by(20)α=tan-1⁡(|wu|).
The angle of sideslip is calculated by(21)β=sin⁡-1⁡(|vV|).
Note that the angles of attack and sideslip are put into the correct phase after this calculation, since the inverse trigonometric function will only ever return values between ±90°.The conversion from the Euler angles to quaternion values is computed by the following set of equations:(22)e0=cos⁡(ψ2)⋅cos⁡(θ2)⋅cos⁡(ϕ2)+sin⁡(ψ2)⋅sin⁡(θ2)⋅sin⁡(ϕ2),e1=cos⁡(ψ2)⋅cos⁡(θ2)⋅sin⁡(ϕ2)-sin⁡(ψ2)⋅sin⁡(θ2)⋅cos⁡(ϕ2),e2=cos⁡(ψ2)⋅sin⁡(θ2)⋅cos⁡(ϕ2)+sin⁡(ψ2)⋅cos⁡(θ2)⋅sin⁡(ϕ2),e3=sin⁡(ψ2)⋅cos⁡(θ2)⋅cos⁡(ϕ2)+cos⁡(ψ2)⋅sin⁡(θ2)⋅sin⁡(ϕ2).
The transformation matrix that relates from Earth fixed to body axis is calculated using quaternion values:
(23)CBE=[e02+e12-e22-e322(e1e2+e0e3)2(e1e3-e0e2)2(e1e2-e0e3)e02-e12+e22-e322(e2e3+e0e1)2(e0e2+e1e3)2(e2e3-e0e1)e02-e12-e22+e32].
The flight path angle is calculated using a sequence of C transformation matrices. There exists a problem in extracting the flight path angle as it approaches 90°, in that there exists a singularity in the inverse sine function. To get around this problem, a third axis system has been introduced that is rotated 90° in pitch with respect to the wind axis. The flight path angle is then given by the following:(24)CWB=C1(0)⋅C2(α)⋅C3(-β),CBW=CWB-1,CGW=C1(0)⋅C2(90°)⋅C3(0),CWE=CWB⋅CBE,CGE=CGW⋅GWE,γ=sin⁡-1⁡(-CGE(1,3))-90°.
The coordinates of the payload along the trajectory are calculated from(25)s=x2+y2+(h-h0)2,
where h0 is the initial altitude.The atmospheric profile (density and speed of sound) is approximated with a cubic curve fit of the reference atmosphere.A second simulation model (designated as SM2) was defined, implemented, and validated in [23]. The major differences with respect to the less complex SM1 model are(i)
the aerodynamics of the bodies are defined in terms of total angle of attack due to the geometrical symmetry of payload and parachute;(ii)
the dynamics of the parachute are modeled with a full-state rigid body six degree-of-freedom representation;(iii)
the added masses of the parachute are defined by a tensor including rotational inertial properties [18, 19];(iv)
the aerodynamics of the parachute include lift and damping coefficients;(v)
the suspension system is represented by a more realistic layout with distributed elements and suspension links whose strain is estimated with an iterative numerical method;(vi)
the effects of atmospheric turbulence and asymmetries (either geometrical or inertial unbalance) can be included.


## 4. Software Implementation

The simulation software is written in Fortran language. This program simulates the dynamics of the parachute-payload system (time domain integration). The user is offered predefined initial conditions from reference flight conditions or the option of trimming the payload at any desired altitude and using those generated initial conditions. The inertial properties of the payload (mass, moments of inertia, and position of the center of gravity) and the parachute characteristics (size, mass, opening, and staging sequences) are defined for a selected altitude ranges (mission table lookup).The solver used in the simulation program is DASSL (differential algebraic system solver) originally developed by Petzold [24–26] for the solution of systems of differential algebraic equations (DAEs). This routine is based on backward differentiation formulas (BDFs). The solution algorithm then attempts to make the residual r(t) equal to zero at each time step:(26)r(t)=f(ẋ,x,t)=ẋ-f′(x,t).
Equation (26) also shows that it is very simple to modify existing explicit ODE (ordinary differential equation) systems to be used with an implicit formulation. With DAE solvers [27], the equations of motion can be implemented directly in the form of residuals. Therefore, no symbolic expansions are needed to identify acceleration terms, and no ad-hoc algorithms need to be used to determine the vector of derivatives at every time step. The DAE solver may automatically adjust the order of the integration formula and the integration step size to achieve the desired accuracy, whereas typical solvers have fixed order, fixed step size, and no automatic accuracy control.The initial conditions for the payload (with or without decelerator deployed) may be computed for a given initial altitudeh0, returning the quasi-trim values for V velocity, α angle of attack, β angle of sideslip, θ pitch angle, and γ flight path angle (stabilized fall). Note that it has been decided to leave roll as the untrimmed equation. The purpose is to solve a system of n simultaneous nonlinear equations in n unknowns. It solves the problem f(x)=0 where x is a vector with components x(1),…,x(n) and f is a vector of nonlinear functions. Each equation is of the form:(27)fk(x(1),…,x(n))=0,fork=1,…,n.
The algorithm used for the solution is based on an iterative method, which is a variation of Newton’s method using Gaussian elimination in a manner similar to the Gauss-Seidel process. Convergence is roughly quadratic. All partial derivatives required by the algorithm are approximated by first difference quotients. The convergence behavior of this code is affected by the ordering of the equations, and it is advantageous to place linear and mildly nonlinear equations first in the sequence. The convergence is started from statically stable conditions (initial payload tumbling is avoided). Numerical integration was performed assuming a time step of 0.001 s with an absolute and relative error tolerance of 10−7. Longitudinal and directional planes are considered separately. The system of equations for the longitudinal plane is defined as follows:(28)f1=M,f2=X-mgsin⁡θ,f3=Z+mgcos⁡θ,f4=θ-γ-α,
where(29)M=12ρV2SdCM+FRcos⁡α(zcg-zp)d+FRsin⁡α(xcg-xp)d,X=12ρV2SCX-FRcos⁡α,Z=12ρV2SCZ-FR⋅sin⁡α.


## 5. Results

The reference mission described in [23] was used to evaluate the two simulation models. The data available from the reference flight tests [28, 29] are used for the validation of the simulation of capsule terminal reentry dynamics. The capsule is substantially a 70%-scaled version (diameter is 2.8 m) of the original Apollo Command Module decelerated by a single conical ribbon drogue parachute (nominal diameter 5.8 m) inflated in two reefing stages and a cluster of three main polyconical slotted parachutes (nominal diameter 22.9 m) inflated in three reefing stages. The sequence of inflation (see Table 1) and the parachute drag profile are time-scheduled according to the planned reefing stages and inflation times. Note that staging is obtained through parachute reefing (i.e., through restricted canopy deployment) with the purpose of limiting peak loads and deployment shocks. A detailed description of parachute deployment sequences for aerospace applications is given in [2, 3].Table 1
The parachute deployment sequence.
SequenceTime (s)Altitude (m)Mortar firing (Pilot Jettison)0.0013994Back cover separation—13666Drogue snatch/1st stage inflation3.3413430Drogue 2nd stage inflation/disreef9.3412546Drogue bridles cut—6461Main snatch/1st stage inflation82.306364Main 2nd stage inflation88.306092Main 3rd stage inflation/disreef94.555950Main bridle 1 cut—5360The aerodynamic data for the payload—which the simulation uses as input—were obtained with wind tunnel static and forced oscillation tests performed at Politecnico di Torino [10]. The wind tunnel is a closed-loop tunnel with a cylindrical working section of 3 m diameter by 5 m length. The model reproduces the scaled geometry of the NASA Apollo Command Module (model diameter is 340 mm). The longitudinal loads on the scaled vehicle (axial force, normal force, and pitching moment) are measured by an internal 3-component balance fit within the capsule model. Body axes are adopted for the reduction of experimental results centered in the actual reference center of gravity. The angle of attack is set to 180° when the spherical base (thermal shield) is exposed to the wind. The model is supported (see Figure 4) by a vertical strut that can rotate over the complete range of angle of attack (α=±180°). Positioning is performed by a step-motor that is controlled by a driving unit interfaced to the data acquisition PC. The stability derivatives are evaluated according to the small amplitude direct forced oscillation technique [30]. The oscillation is generated by the step-motor (Δθ=1°-5°) and f=0.5–5 Hz. The position of the driving shaft is acquired by means of a digital encoder. The measurement repeatability for the averaged static coefficients is σ≈±0.1% estimated over the full range while the measurement repeatability for the longitudinal stability derivatives (measured from frequency response) is σ≈±1%.The experimental setup [10].
(a)(b)The reentry profile is correctly reproduced by the two simulators—as presented in Figure5—and the trajectory, shaped by the parachute opening sequence, is matched by both SM1 and SM2 models. The trend of the load factor exhibits a set of marked peaks given by the staging of both main and drogue parachute. The SM1 results exhibit an oscillatory behavior in the second part of the parachute deployment phase, induced by the oscillatory behavior of the payload, not found for SM2 model results. Note that low-altitude flight test data are affected by the asymmetry of the suspension system induced by the cut of one of the bridles of the main parachutes (see Table 1).The simulation of the trajectory profile for the parachute-payload system.
(a)(b)(c)(d)The ability of the SM1 model to estimate the riser and parachute loading is outlined in Figure6. The trend of forces is coherent with the deceleration profile presented in Figure 5. The offset between the force acting on the canopy and the riser is well marked for the main clustered parachute, modeled as an equivalent single decelerator. This offset is due to the mutual influence of the added masses (inertia-induced delays) and the elongation of the riser. The effect of the elasticity of the main cable is also visible in terms of damped oscillatory strain, triggered by the bouncing of the payload (and by its attitude dynamics) after each inflation phase (staged reefing).The loads acting on the riser and the decelerator (drogue and main parachute).
(a)(b)(c)(d)The projected trajectory is compared in Figure7. Accurate trajectory coordinates were not available from flight test data. The profiles (mirrored for the purpose of checking their symmetry) are very similar in the vertical plane. Nevertheless, a slightly different path is found comparing the traces for the lateral-directional plane. This is the impact of the aerodynamic model adopted for the payload, which for the SM1 case decouples and superimposes the effects for the longitudinal and the lateral-directional planes. The use of total angle of attack for the interpolation and reconstruction of payload aerodynamic coefficients (as in SM2 model) provides a more accurate fit.The simulation of the trajectory profile for the parachute-payload system.
(a)(b)The aerodynamic angles of the payload are plotted in Figure8. The major discrepancy between the two models is the level of dynamic stability of the modal response (short-period dynamics). This can be explained with the extreme sensitivity of the aerodynamic coefficients to center of gravity location and angle of attack (see Figures 9 and 10). Remind that the center of gravity location is updated during the descent, accounting for covers, deployment bags, and parachutes release. As a matter of fact, minor parametric changes (α, xcg, and zcg) shift the attitude for the rotational equilibrium of the payload and alter its dynamic stability. This point is also a concern for the validation of simulation models by comparison with flight test data [23], in which a transition between different equilibrium states and levels of dynamic stability can be externally forced by the atmospheric perturbations, exciting the combined parachute-payload system. Another contributing element is the dynamics of the parachute that in the higher fidelity model (SM2) allows for attitude changes decoupled from the suspended body, providing a more realistic behavior of the parachute-payload system. Differently, the parachute trajectory in the SM1 model is aligned instantly with the velocity vector (flight path) of the payload. Furthermore, a multilink suspension system (as modeled in SM2 implementation) induces a stronger steering effect on the payload.The simulation of the angles of attack and sideslip for the payload.
(a)(b)(c)(d)Figure 9
The aerodynamic coefficients of the payload (pitching moment static stability) [10].Figure 10
The aerodynamic coefficients of the payload (pitching moment dynamic stability) [10].In order to verify the modal response, at least for the longitudinal short-period natural frequency, the pitch rate spectral response for the two simulation models is compared with the elaboration of available flight logs in Figures11 and 12 (note that from flight test reports gyro signal measurements are very noisy after drogue release). The comparison of the spectral data is based on the crosscheck of the index of similarity (see Table 2), defined according to (30), where Ai is the amplitude for the given frequency ωi: (30)s0=∑i=1NAi,s1=∑i=1NAi⋅ωi⟶σ1=s0Nσ2=s1s0σ3=s2s0⟶I1=σ1aσ1bI2=σ2aσ2bI3=σ3aσ3b︸a→bs2=∑i=1NAi⋅ωi2.
The results show that the range for natural frequency of the short period response is matched by both simulation models as demonstrated by the fact that I2 and I3 are close to unity. Other than that, the averaged amplitude of the spectral response is less precisely reproduced by the simulation models, as experimental data represent a modal response that is overexcited by atmospheric disturbances and suspension system asymmetries, mainly at parachute deployment.Table 2
The index of similarity for the pitch rate spectral response.
Flight test⇒ SM1Flight test⇒ SM2Drogue parachuteDrogue parachuteI1I2I3I1I2I31.3921.0050.9890.6451.1451.110Main parachuteMain parachuteI1I2I3I1I2I30.5781.0221.0290.4631.0881.098The effect of model complexity on pitch rate spectral response (drogue parachute).
(a)(b)The effect of model complexity on pitch rate spectral response (main parachute).
(a)(b)


## 6. Conclusions

The present work outlines a comparative analysis of two simulation models (SM1 and SM2) of a parachute-payload system with different levels of complexity. An example of limited complexity reconfigurable simulation models for a payload decelerated by one or more parachutes is given, including details and implementation features usually omitted as the focus of the research in this field is typically on the investigation of mission design issues, rather than addressing general implementation guidelines for the development of a reconfigurable simulation tool. The dynamics of the system are modeled through a simple multibody model that represents the expected behavior of an entry vehicle during the final deceleration phase by means of a system of aerodynamic decelerators. The simulators are designed according to a comprehensive vision that enforces the simplification of the coupling mechanism between the payload and the parachute, with an adequate level of physical insight still available in terms of sensitivity of system performance metrics to design parametric changes.A reference mission described was used to evaluate the two simulation models. The aerodynamic data for the payload—which the simulation uses as input—were obtained with wind tunnel static and forced oscillation tests performed at the Politecnico di Torino.The trajectory is reproduced by the two simulators as the profiles, mainly shaped by the parachute opening sequence, are matched by both SM1 and SM2 models, even if the use of total angle of attack for the interpolation and reconstruction of payload aerodynamic coefficients (as in SM2 model) provides a more accurate fit. The riser and parachute loading are estimated in a way that is coherent with the deceleration profiles.The major discrepancy between the two models is the level of dynamic stability of the modal response (short-period dynamics). This can be explained with the extreme sensitivity of the aerodynamic coefficients to center of gravity location and angle of attack. Another contributing element is the dynamics of the parachute which in the higher fidelity model (SM2) allows for attitude changes decoupled from the suspended body, providing a more realistic behavior of the parachute-payload system. Differently, the parachute trajectory in the SM1 model is aligned instantly with the velocity vector of the payload, probably a too crude approximation if the purpose of the analysis is the attitude dynamics of the suspended vehicle. Furthermore, a multilink suspension system (as modeled in SM2 implementation) induces a steering effect on the payload, neglected by the single-point suspension of SM1 model.The spectral results show that the range for the natural frequency of the short-period response—as measured in flight—is matched by both simulation models as demonstrated by the fact that the indices of similarityI2 and I3 are close to unity. Differently, the averaged amplitude of the spectral response is less precisely reproduced by the simulation models, as the experimental data represent a modal response that is over-excited by atmospheric disturbances and suspension system asymmetries, at parachute deployment mainly.As a general conclusion, the lower fidelity simulation tool SM1 exhibits a limited advantage in terms of computational workload, mainly providing a simplified approach for preliminary trajectory estimation and parachute sizing.



---

*Source: 182907-2012-05-29.xml*