Review Article

Proc IMechE Part G:

J Aerospace Engineering
Numerical simulation of parachute 0(0) 1–31
! IMechE 2017

inflation: A methodological review Reprints and permissions:

[Link]/[Link]
DOI: 10.1177/0954410017705900
[Link]/home/pig

Yihua Cao1, Shuai Nie1 and Zhenlong Wu1,2

Abstract
The parachute inflation process involves fluid–structure interaction problems posing several mathematical and engin-
eering challenges, e.g. accurate aerodynamics calculations for bluff-body geometries involving with moving boundary,
appropriate structural models in predicting the behavior of canopy, and realization of the coupling between the fluid and
structure. These challenges attract the attention of scholars worldwide, and considerable achievements have been
obtained in applying numerical methods and simulations to design multifarious parachutes. In this paper, the authors
highlight the advances in the following fields: the methods suitable for time-dependent flow around bluff-body geome-
tries, the accurate structural models in consideration of the under-constrained and no-compression nature of the
canopy, and the advantages and disadvantages of different coupling algorithms in terms of numerical stability and com-
putational economics. Moreover, in order to simulate the parachute inflation more realistically, we focus on accurate
representation of three physical phenomena, as follows: an appropriate model of the flow through porous media, an
accurate treatment of the wrinkling phenomenon of the canopy, and a consistent representation of the impact-contact
problem associated with the inflation process. Finally, based on a review of existing literature, we offer recommendations
for future research on the application of numerical methods for simulating the inflation process.

Keywords
Parachute inflation, VEM and grid-based CFD, membrane model, porosity, wrinkling, impact-contact

Date received: 21 November 2016; accepted: 28 March 2017

Introduction experimentation. Numerical methods for simulating

The parachute is a preferred aerodynamic decelerator parachute inflation have attracted significant atten-
for addressing airborne scenarios such as airdrops, tion from the parachute community because para-
ejection escapes, lifesaving, and aircraft recovery. It chute designers believe that numerical simulation
has also become clear that the parachute descent tech- methods can assist with designing parachutes of arbi-
nique is indispensable for space exploration, as trary geometrical shapes and structures, particularly
demonstrated, for example, by the deceleration sub- under complex working conditions. After completing
system of the Mars Science Laboratory (MSL)1 cre- the design process with numerical methods, the
ated by NASA. The extensive application of related wind tunnel and airdrop tests can only serve
parachutes likely stems from their small folding verification purposes. Numerical simulation methods
volume, lightweight and outstanding decelerating have considerably shortened the parachute design
effects. During an aeronautics or astronautics mis- cycle-time, and the cost can also be reduced by elim-
sion, the shape and velocity of a parachute change inating unnecessary experiments. Furthermore, wind
significantly. This process of change can be divided tunnel tests cannot realistically simulate the perform-
into the following three stages: deployment, inflation ance of a parachute under actual Mars atmospheric
and terminal descent.2 Of the three stages, inflation is
of particular engineering significance and particularly 1
School of Aeronautic Science and Engineering, Beihang University,
subject to technical difficulties. Beijing, China
2
Research on the parachute inflation process has National Laboratory of Aeronautics and Astronautics, Beihang
been conducted for nearly a century, during which University, Beijing, China
time the methodology has changed significantly,
Corresponding author:
from an experimental method to a semi-experimental Zhenlong Wu, Beijing University of Aeronautics and Astronautics,
and semi-empirical method to finally, the current Xueyuan Road, 37th Haidian District, Beijing 100191 China.
numerical simulation associated with theory and Email: jackilongwu@[Link]
2 Proc IMechE Part G: J Aerospace Engineering 0(0)

conditions. Nonetheless, numerical simulation meth- and F and R represent the convective and diffusive
ods can help the simulations better approximate real fluxes, respectively. Because the Navier–Stokes equa-
conditions if the Mars atmospheric parameters are tions are strongly nonlinear equations, calculating
known, and these parameters are typically easily them numerically is a herculean task. A detailed
obtainable. description of solving the Navier–Stokes equations
The numerical method for simulating parachute (consisting of discretization scheme, simplification of
inflation seems to offer unique advantages. the equations, algorithm and turbulence model
Nonetheless, applying numerical method also encoun- according to different conditions) is described.4–11,27,80
ters with formidable challenges. Firstly, the parachute Solving Navier–Stokes equations for different airflow
inflation is intrinsically a nonlinear, large-displacement conditions for parachute inflation is discussed in
fluid–structure interaction phenomenon: the sur- Chapter 3.
rounding airflow field that generates aerodynamics
forces acting on the canopy depends on the shape of
The structure subproblem
the canopy, while the shape of the parachute depends
on the aerodynamic forces.3 The irregular movement Typically, the wetted canopy (and cable) of the system
of canopy increases difficulties in predicting the aero- is considered the structure component, whose motion
dynamics of it. Secondly, the difficulties also arise from is governed by Newton’s Second Law, given certain
that the mechanical behaviour of the canopy is unpre- internal constraints. Because the canopy (or cable) is a
dictable due to the under-constrained and non- flexible structure and a continuum system, the motion
compressive nature of the flexible fabric used to of the moving canopy can be described using
make the parachute. Thirdly, the choice of coupling Lagrangian form with general equations, assuming a
algorithms remains challenging in terms of the preci- spatial structural mechanics domain st with bound-
sion and stability of the data transfer. Finally, to simu- ary st at time t 2 ð0, TÞ
late this dynamic process more realistically requires
full consideration of various details, e.g. the porosity d2 xs dxs
of fabric, the wrinkling phenomenon of canopy and the  2
¼ f c þr ð2Þ
dt dt
impact-contact problems. Therefore, it is crucial to
identify the numerical methods employed in simulating where at position xs ,  denotes the density of fabric;
the parachute inflation process. A large volume of f denotes the external load exerted on the canopy that
related research has been conducted, resulting in sev- corresponds to the gravity of structure, the fluid trac-
eral significant insights. tion, the proportional pressure and the contact force
Considering the lack of studies that integrate vari- (ubiquitous in a folded inflatable structure); c denotes
ous methods to perform a numerical simulation of the the damping coefficient; and  denotes the Cauchy
parachute inflation process, this paper aims to collect stress tensor. The derivatives in equation (2) should
research in this field as a reference for future work. generally be considered partial derivatives.
Based on existing research on the numerical methods Additionally, the considerable displacement and rota-
used to simulate parachute inflation, the deficiencies tions involved make this problem nonlinear. Thus, to
of the proposed simulation models are summarized, solve the equation of structural motion precisely,
and future research directions are suggested towards equation discretization is necessary. Many studies
devising more detailed and more accurate numerical simplify the structure of canopy (and cable) into a
simulation models for parachute inflation. different force capacity model and adopt certain
assumptions based on engineering requirements. The
constitutive equation for determining the structural
Problem statement component also differs based on material. The discret-
ization scheme of the generalized structure equation
Formulation of the fluid subproblem
and the different structural model of the structural
The time-dependent Navier–Stokes equations of component are described in Chapter 4.
mass, momentum and energy can be calculated for a
domain  with a nearly continuous boundary @ as
The transmission conditions

@ ðJWÞ The fluid and the structure influence each other at
þJrx  FðWÞ ¼ Jrx  RðWÞ ð1Þ
@t  their coterminous interface through so-called trans-
mission conditions. The fluid–structure coupling is
where equation (1) is non-dimensional conservative accomplished by exchanging nodal information
form of the Navier–Stokes equations. Here, t denotes along the mutual boundary. The course of this cou-
time; x denotes the position in Cartesian coordinates pling occurs through a data transferring process
associated with fluids;  denotes the position in refer- between nodes in the fluid dynamic (FD) model and
ence configuration; J ¼ detðdx=dÞ; W is the fluid the structure dynamic (SD) model. Generally, the
state vector comprising the conservative variables; transmission conditions can be divided into the
Cao et al. 3

following two classes of conditions: dynamic and with the ensuing complex geometries, numerous
kinematic.12 The dynamic condition describes the experiments were conducted on either scaled-down
continuity of the load, whereas the kinematic condi- parachutes or full-scale parachutes. Most of the
tion accounts for the continuity of motion. To address experiments focused on rigid bodies to observe the
the parachute inflation problem, the coupling equa- wake behind axisymmetric geometries. The structure
tion for the dynamic condition can be written as of the ‘‘steady-state’’ wake behind axisymmetric
Given the arbitrary wetted boundary t bodies such as discs, cups, and even spheres has
been extensively studied in cases of moderate to
S  n ¼ p  n þ F  n ð3Þ high Reynolds numbers. It is summarized that the
wake of a bluff body is inherently unstable due to
where p denotes the pressure field; S and F denote the presence of a recirculation bubble and the accom-
the structure stress tensor and the fluid viscous stress panying reverse flow and that the flows beyond indi-
tensor, respectively; and n denotes the normal vector vidual canopies or between ribbons generate complex
at the boundary. time-dependent interactions.3
The coupling equation of the kinematic condition The near wake of an inflating flexible canopy was
can be written as also investigated.13 The temporal evolution of the vel-
ocity field immediately downstream from the canopy
@xs was measured. It is found that the inflation of the
¼ uF ð4Þ
@t canopy occurred in the following three stages: a
shear layer fully attached to the surface of the
where xs and uF denote, respectively, the displacement canopy; a shear layer separated from the canopy sur-
field of the structure and the velocity field of the fluid. face from the apex of the parachute; and the shear
This equation illuminates the compatibility between layer completely detached from the canopy surface
the velocity fields of the structure and the fluid at at the skirt edge. The near wake of an inflating para-
the solid interface. For an inviscid flow, equation (4) chute is shown in Figure 1.
is transformed to apply to a slip wall boundary

@xs The VEM

uF  n ¼ n ð5Þ
@t The numerical simulation of a flow field involves the
vortex element method (VEM) and the conventional
The kinematic condition of the boundary interface grid-based computational fluid dynamics (CFD)
can be interpreted as the Neumann boundary condi- method. Vortex element methods place vortex elem-
tion for solids and the Dirichlet boundary condition ents in regions where vorticity is present; then, a diver-
for fluids. The stabilization of the boundary interfaces gence free flow field can be substituted by the vorticity
coupling process, the stability analysis of the coupled field. The VEM was selected for several reasons. First,
FSI computations, data transmission and the strategy it is a grid-less method for solving the unsteady
of moving mesh are detailed in Chapter 5. Navier–Stokes equations. Compared with the time-
A penetrable boundary, a common instance consuming grid-based CFD, the VEM is more com-
referred to when considering the porosity of the putationally efficient because significant effort is not
canopy, should be treated distinctly. This distinct required to reconstruct the grid of the fluid field at
treatment will be detailed in Chapter 6. Chapter 6 every time step due to changes in radical geometries
also highlights certain tasks, including addressing during the parachute inflation. Additionally, the
the wrinkling phenomenon of the canopy and the VEM is more stable because of less numerical diffu-
impact-contact problem encountered during para- sion associated with the transport of the vorticity.
chute inflation. The paper concludes and provides rec- Nonetheless, the VEM also poses disadvantages.
ommendations for future research in Chapter 7. The greatest disadvantage is the limit of the applica-
tion condition. Traditionally, VEM is used to model
the incompressible flow with either a moderate or a
Parachute aerodynamics relatively high Reynolds number. In transonic and
Parachute aerodynamics is an issue for in-depth con- supersonic flow conditions, where certain parachutes
sideration for parachute designers. An important work, the feasibility of applying the VEM in a com-
issue within parachute aerodynamics concerns the pressible regime must be extended in the future.
flow around complicated shapes. Parachute designers To understand the distinct mechanism of the VEM,
prefer bluff-body geometries to model the parachute it is helpful to consider the vorticity transport
aerodynamics, given the characteristics of parachutes. equation
To model the surrounding flow field properly, the
parachute designer should clarify the flow regime D! @!
¼ þ u  r! ¼ !  ru þ r2 ! ð6Þ
around the parachute. To determine the massively Dt @t
separated unsteady flow characteristics associated
4 Proc IMechE Part G: J Aerospace Engineering 0(0)

the shear stress near the canopy surface with the

velocity gradient for the no-slip boundary condition,
the VEM must not do so, though it can be used to
resolve the boundary layer velocity gradient.
However, when introducing the boundary condition
into the VEM, it must be considered that considerable
insight and experience are required even for a 2D
flow.14,15
The 2D vortex-discrete method has been utilized
extensively to solve the flow around a parachute-like
geometry: from a finite wake flow16 to a complete flow
field description,17 from the assumption of a fixed
vortex shedding offset point18 to the implementation
of a dynamic separation point,19 from deliberately
ignoring the no-slip condition to fully considering
the viscosity influence,20–24 and did really provide
valuable insight into the aerodynamic characteristics
of a parachute system and improved the understand-
ing toward the aerodynamic performance of para-
chute. However, the ultimate goal is to apply VEM
to establish a guideline for designing a parachute prior
to experiments by understanding the behaviour of 3D
canopies.25 Accordingly, 3D vortex methods have
attracted researchers’ interest. Dating back to the vor-
ticity transport equation (6), the term !  ru, which
represents the change in vorticity due to vortex
stretching, remains for three-dimensional cases. In
fact, the stretching of vortex lines concentrates vorti-
city, prompting an increase in velocity fluctuations
and a decrease in the minimum length scale in the
flow. In other words, the turbulence effects contribut-
ing to the vortex lines stretching will also be respon-
sible for the turbulence. It is the presence of the
stretching term that results in the complexity of
three-dimensional flow and increases difficulties in
simulating the flow with a vortex method such as
vortex merging, a change of vorticity and a geometric
configuration variation of vortex from vortex stretch-
ing. Despite the manifold difficulties encountered in
this field, pronounced advances have also been made.
According to the different discretization schemes of
vorticity, many corresponding approaches emerge.
Currently, the most frequently used 3D vortex meth-
Figure 1. Vorticity fields of an inflating canopy.13 ods are the vortex filaments, the vortex particle and
the vortex-in-cell (VIC) methods.
As is well known, many unsteady three-dimen-
sional flows around a bluff-body can be characterized
Note that it is feasible to analyse the two-dimen- by high-Re flow and the presence of high vorticity
sional incompressible flow by the VEM regarding the embedded in an irrotational flow. Wakes, free-shear
transport equation of vorticity because the term layers, vortex rings, and trailing vortices are the main
‘‘vortex stretching’’ in the velocity gradient (!  ru flow patterns typifying this flow. Generally, the vorti-
only exists in three-dimensional space) vanishes, and city appears through viscous interaction between
the pressure does not appear explicitly in the govern- fluids with the solid surface and then ‘‘grows’’ to be
ing equation, thus indicating the velocity field can be the shear-layer or the vortex sheet. Under normal con-
determined without calculating the pressure field. ditions, this shear-layer or sheet will either subdivide
Hence, the velocity field can be determined from the into particles resembling thin tubes or sticks or roll up
vorticity field by integrating the influence of the vor- into filaments. The vortex filaments method for 3D
ticity field on any point in the domain. Compared flows were likely first proposed by Lenard.26 An arbi-
with the grid-based CFD method for modelling trary curve in space can be approximated as a
Cao et al. 5

sequence of linear segments between successive nodes. accurate treatment of boundary conditions33; and
Thus, the continuous filaments represented by splines the capability of performing an accurate particle
with integration in space and time enable simulating redistribution in the presence of a boundary.34 Due
the evolution of the shear layer. It is widely believed to the reduced computational cost of grid properties,
that the strength of a slender vortex filament depends many researchers tend to use temporary grids for the
on its core structure and that viscous stresses and the vortex particle method (also known as the vortex-in-
stretching of the filament length prompt the changes cell method) in addressing flow computation without
of the core structure. A detailed description of the regarding the grid-free nature of VEM. The VIC
vortex filaments method can be found in Leonard.26,27 method is inherently a hybrid Lagrangian–Eulerian
Generally, the continuous-vortex-filaments repre- algorithm, which uses a grid to calculate the velocity
sentation of vorticity can be replaced by use of field, but particles are used to track the vorticity. The
vortex balls such as vortons, vortex arrows or motivation for using the particles method is reducing
vortex sticks. Inspired by the inference, the vortex the dimensionality of the phase space, and elegantly
particle method for 3D flows of homogenous incom- and robustly addressing the transport terms in the
pressible inviscid fluids was proposed.28 These seg- equations, whereas grid solvers can rely on the Fast
ments (regarded as particles) characterized by a Fourier transform (FFT) to evaluate the fields at a
function of variables, including the core radius, geo- minimal numerical cost.35 Great strides have been
metric configuration, circulation and vector elements, made by implementing the VIC method for a three-
could be used to describe the flow regime in the wake. dimensional flow in aspects of celerity, accuracy and
The particle is required to account for the diffusion robustness, including an appropriate boundary treat-
effects. Leonard29 proposed a method to address the ment,33 a truncation error correcting,36 an application
vortex stretching analogue of core-spreading tech- of VIC to a direct numerical simulation,37 a coupling
nique. He assumed that the local differences in method between a vortex method with a penalization
vortex stretching created local variations in the core method38 and a higher computational efficiency for
radius and the geometric configuration. Because of implementing VIC algorithms.39
the feasibility of various flows, the vortex particle The accuracy of different three-dimensional vortex
method gained rapid development, for example, in element methods was validated by a series of numer-
the following respects: a rapid algorithm with multi- ical experiments.40–42 The simulation results revealed
pole expansions30 and an active control error31; an that different methods tended to coincide with the
exact approach of viscous effects by adopting a par- observed image.43 Figure 2 shows the results of differ-
ticle strength exchange scheme (PSE32) and an ent methods and experiments. The results of different

Figure 2. (a) Flow regime structure generated by the vortex filaments method (by Summers (Re ¼ 300))40; (b) Calculated with
the vortex particles method (by Lavi (Re ¼ 100))41; (c) Calculated with the vortex-in-cell method (by Kim (Re ¼ 100))42; (d) Flow
visualization (by Marcelo (Re ¼ 297)).43
6 Proc IMechE Part G: J Aerospace Engineering 0(0)

methods all revealed a ‘‘bubble-like’’ vorticity for a

long period behind the sphere at a moderate
Reynolds number. This result is consistent with the
experiments except subtle discrepancies due to the tur-
bulent random effects.
Shriyama44 pioneered the application of the 3D
vortex method for simulating flows around
‘‘parachute-like’’ objects. He proposed a series of
vortex sticks to use as the elements to simulate the
wake flow field of the parachute. A set of vortex
blobs was placed along the boundary layer, and the
vorticity convection from the boundary layer to the
near-wake could be conceptually replaced by the tran-
sition from vortex blobs to vortex filaments. These
vortex filaments disconnect from each other and sep-
arate into vortex sticks when their strengths reach a
given cut-off value. A rigid circular disc and part of a
sphere were modelled with his method, and the results
turned to be inspiring. Unfortunately, the vortex
method employed by Shriyama and his colleagues
was neither sufficiently detailed nor definitively con-
clusive, and the description of their numerical experi-
ments was imprecise and vague. Nonetheless, the
vortex sticks method, as a hybrid method of the
vortex filaments and vortex particles methods, could
be considered as a good attempt.
Additionally, the Sandia Laboratories extended the
2D vortex method code VAPRA (Vortex method for
PARachutes in Axisymmetric flow) to a 3D flow
field.45 Similar to the 2D VAPRA vortex distribution
method, a vortex sheet was attached to the canopy
surface to produce a normal velocity at the surface.
Then, the vortex sheet was split into two sheets inside
and outside the surface to satisfy the no-slip condi- Figure 3. (a) Wake behind a three ring ringslot; (b) Wake
tion, while the sum of the two new sheet strengths was behind a rotating cross parachute.45
set equal to the strength of the original sheet. The
sheets then diffused away from the surface and con- VEM is limited in low-speed flow conditions where
verted into particles (vortons) to adhere to the flow certain parachutes work, the feasibility of applying
regime in the near-wake of converting, stretching and the VEM in a compressible regime must be extended
diffusing to satisfy the vorticity transport equation. in the future. Second, all these existing methods are
This code for simulating 3D flow around several para- devised to address those issues with the presumption
chute types tends to be remarkable. However, this that the parachute is rigid. The interaction between
method also needs to be examined with more complex the flexible fabric and the vortex needs to be given full
geometries and more realistic circumstances. consideration. Moreover, the VEM codes are far from
Recently, the improved VIC method is frequently mature and the generality of the VEM is also pursued
applied in finding a solution of the flow field around in the future.
bluff bodies such as an airliner, a parachute,46–48
including a modified numerical scheme of the VIC
method that allows an effective application of the
Conventional grid-based CFD methods
mosaic-skeleton approximation of matrices for the Grid-based CFD methods were conceived following
evolution of a 3D vorticity field in an unbounded the availability of advances in computer science and
domain and the implementation of a parallel algo- remain tightly coupled with the evolution of computer
rithm for the 3D vortex method on a multiprocessor science. With advances in the development of numer-
platform. The latest improved algorithm can be ical methodologies, significant progress was made in
applied to solve aerodynamics problems around clus- grid generation methods with respect to such aspects
ters of parachutes. Although the VEM is beginning to as the robustness and accuracy of discretization
be used as a tool to predict the flow-field of the para- schemes and efficient solution techniques for viscous
chute, this method is far from mature and needs to be flows. There is a strong interest in linking CFD with
improved in the future. First, the application of the other disciplines (such as structural mechanics or heat
Cao et al. 7

conduction). Compared with the VEM, the grid-based

CFD method seems more mature, and a large number
of commercially available codes were developed to
solve compressible or incompressible flows. Based
on different principles of the governing equation spa-
tial discretization, CFD methods can be divided into
the following three branches: finite difference meth-
ods, finite element methods and finite volume meth-
ods. This paper is concerned with the application of
different CFD methods to problems of parachute
aerodynamics.
A large number of mature CFD codes based upon
different discretization schemes were applied in pre-
dicting the parachute aerodynamics, and numerous
achievements were made.49–62 Among them, Cao
et al. have been endeavouring to employ the finite
volume CFD method to generate solutions for the
flow field in and around different parachutes. Their
attempted methods include the semi-implicit method
for pressure-linked equations (SIMPLE, a code for
incompressible Navier–Stokes equations) algorithm
with the k  " two-equation turbulent model for a
flat circular parachute, a conical parachute46 and a
parafoil54; the semi-implicit method for pressure-
linked equations consistent (SIMPLEC, an improved
code relative to SIMPLE for incompressible Navier–
Stokes equations) algorithm to solve shear stress
transport (SST) k  ! two-equation turbulence
Navier–Stokes equations for a hemispherical para-
chute55; and mixed Jameson/total-variation-diminish-
ing (TVD) schemes to solve Spalart-Allmaras one-
equation Navier–Stokes equations for a flat circular
parachute.56 Reasonable results were obtained. All
the methods proved effective in predicting the drag
coefficient and pressure distribution of parachutes
according to engineering practice. Figure 4 shows
some results from Cao et al. applying their methods Figure 4. Some results for different parachutes: (a) Flat cir-
to different parachutes. cular parachute, (b) hemispherical parachute, (c) parafoil.
Note that the aforementioned CFD simulations are
based mainly on the assumption that the common _ to trace the tran-
term, the fluid grid point velocity x,
boundary separating the fluid and the structure is sta- sient fluid grid position. Consequently, equation (1)
tionary. Considering the nature of the fabric, this can be rewritten
assumption deviates from the realistic life with the 
time-variation effect ignored, thus rendering the simu- @ ðJWÞ
_
þJrx ðFðWÞ  xWÞ ¼ Jrx RðWÞ ð7Þ
lation results unconvincing, though they provide valu- @t 
able insight into the aerodynamic characteristics of a
parachute. To simulate the parachute inflation pro- where x is a function of time and is used to denote the
cess accurately, free surface flow problems must be time-dependent position or displacement of each fluid
considered. A classical method for analysing free sur- mesh point. Clearly, equation (7) and mesh moving
face flow problems involving large free surface equations (particularized in ‘Interface data transfer’
motions is based on an arbitrary Lagrangian- section) are directly coupled. A detailed solution of
Eulerian (ALE) kinematical description of the fluid the ALE equations of fluid is available in Hirt
domain in which the nodal points are displaced inde- et al.64 The ALE method inherently describes the
pendently from the fluid motion. Combining the grid motion independent of space coordinates based
merits63 of the Eulerian method to describe fluid on an Euler description and material coordinates
motion and the Lagrangian method to describe par- based on a Lagrange description and thus provides
ticle motion can be greatly adaptive for addressing a novel solving strategy for hyperbolic problems
such problems. Specifically, the modified governing involving a moving boundary. Given this philosophy,
equations of fluid motion introduce an additional many relevant versions are developed, including the
8 Proc IMechE Part G: J Aerospace Engineering 0(0)

finite difference-based simplified arbitrary

Lagrangian–Eulerian Navier–Stokes code,65,66 the
finite element-based ALE algorithm67–70 and the
finite volume-based ALE algorithm.71,72 However,
the ALE finite volume method has not yet been
employed in numerical simulations involving para-
chutes, which differs from the former two. Although
the ALE method exhibits advantages in addressing
highly non-linear FSI issues, we found that the ALE
method has some inherent flaws: the stability of this
algorithm is heavily dependent of the mesh quality
and the choice of time-step, and the time accuracy
of ALE method is first-order, which probably dissa-
tisfies the precision requirement.
Another distinct method for determining the flow
field around a parachute is the deforming-spatial-
domain/stabilized space-time (DSD/SST) method.73
The DSD/SST method, a finite element method, exhib-
ited favourable applicability for solving flow problems
involving complex geometries and moving boundaries
and interfaces, particularly for addressing such chal-
lenges as an unsteady flow with interfaces; fluid–struc-
ture interaction; and the aerodynamics of complex
shapes. The example of a parachute in an airdrop
system can be considered with the DSD/SST method.
With a space-time method, the finite element formula-
tion of the above problem is devised over the para-
chute’s space-time domain, while automatically
considering the motion of boundaries and interfaces. Figure 5. The space-time slab as an abstract representation
(top) and in a 2D context (bottom).78
At each time step, the new locations of boundaries and
interfaces are calculated as a part of an overall solution
and are used as the boundary and interface conditions linear in time. Moreover, stabilization techniques
for the next time step. The boundary and interface such as the Streamline-Upwind/Petrov-Galerkin
function requires calculation to capture the motion (SUPG), the Pressure-Stabilizing/Petrov-Galerkin
of boundary and interface. Inherently, the space-time (PSPG), and the Galerkin/Least-Squares (GLS) for-
formulation is an interface-tracking method. mulations are used to prevent numerical oscillations
Specifically, the DSD/SST finite element formula- that might be caused by the advections term and cer-
tion uses interpolation functions, which are piecewise tain combinations of interpolations functions for
bilinear and continuous in space and piecewise linear determining velocity and pressure. The more detailed
and discontinuous in time. To solve the problem given solution of incompressible and compressible flow is
the discontinuity in time, fully discretized equations of specified in Aliabadi and Tezduyar.74
one space-time slab at a time step must be con- The DSD/SST algorithm is extensively exploited in
structed. Figure 5 shows the space-time slab for the studying the aerodynamics of several types of para-
DSD/SST formulation. The domain is divided into a chutes, including a round parachute,75 a ram-air par-
sequence of N space-time slabs Qn, where Qn repre- afoil,76 a cross parachute,77 a ringsail parachute78 and
sents a slice of the space-time domain between the others.79 Outstanding achievements have been
time levels tn and tn þ 1. The integrations involved obtained, and remarkable progress has been made
the weak form are then conducted over Qn. with respect to determining the flow field around a
Introducing different trial solutions and weight func- parachute. Figure 6 shows the fruitful results yielded
tions corresponding to appropriate finite-dimensional by applying the DSD/SST algorithm. Currently, the
space-time function spaces, the different formulations T*AFSM is used for accurately and quickly modelling
of the compressible or incompressible Navier–Stokes parachutes of complex geometry, for sophisticated
equations can be written in integral form to satisfy the applications and in a complicated environment.
Neumann-type boundary condition. Then, the However, the DSD/SST method, generating large
GMRES (Generalized Minimal RESidual) algorithm numbers of finite element equations, requires a high
is used to iteratively solve the implicit system arising performance parallel computing platform, namely
from the space-time discretization scheme to eliminate mainframe computers with a multi-processor.
the highly computational cost associated with space- The immersed boundary (IB) method80 can also be
time formulations based on interpolation functions a viable alternative to simulate flow around moving
Cao et al. 9

Figure 6. Fruitful results obtained with the DSD/SST algorithm. (a) Flow field around a round parachute.75 (b) Pressure distribution
across a ram-air parafoil.76 (c) Flow structure around a cross parachute.77 (d) Velocity vectors around rings and sails.78

boundaries. A distinct difference between the conven- dependent geometry of a boundary distinctly.
tional approach and the IB method exists is that the Instead, the immersed boundary is tracked as a
former engages structured or unstructured grids that sharp Lagrangian entity, and its motion can be inter-
conform to the boundary. Generally, constructing preted as the convection and diffusion of the fluid cell
body-conformal grids is cumbersome, and particu- attached to it. The IB method has been used to simu-
larly in cases involving moving meshes, a new body- late the parachute inflation process by Kim and
conformal grid must be constructed at each time step. Peskin.81,82 Their numerical computations for a para-
Additionally, one seeks a grid that offers adequate chute indicated the potential of the IB method for
local resolution with a minimum of grid points. assisting with parachute design. A comparison of
These conflicting requirements would deteriorate the some of the aforementioned methods is shown in
grid quality and then negatively affect the accuracy Table 1. However, the existing versions of the IB
and convergence of the solver. By contrast, for a method have some inadequacies. For example, the
simulation performed using the IB method, the grid solution of the fluid equations is not efficient enough
complexity and quality are independent of the com- and the numerical simulations presented in their
plexity of the geometry, and the treatment with paper are limited in generalizability by the low and
moving boundaries proceeds relatively simply by moderate Reynolds numbers. The IB method might
employing a stationary, non-deforming Cartesian not possess the generality to address all of the simu-
grid, thus significantly reducing body-conformal lation scenarios even with various amendments.
gird-generating operations. Considering that the IB method’s fluid equations
Conventional methods proceed by imposing the must be solved in a regular Cartesian grid, the lattice
boundary conditions upon the discretization of equa- Boltzmann (LB) method,11 a regular lattice-based
tions (1) on a body-conformal grid directly. With the scheme for fluid flow simulation, can be incorporated
IB method, the immersed boundary moves at the local into the IB method to make an IB-LB coupling. The
fluid velocity and slightly resembles the familiar no- LB method is particularly successful in fluid flow
slip condition of a viscous fluid. The boundary con- applications involving interfacial dynamics and com-
dition is accounted for by the governing equations by plex boundaries, and its simplicity, efficiency, and par-
adding a source item (momentum or force), which allelism have been extensively verified. Combining the
reproduces the effect of a boundary. Thus, the fluid merits of the IB and LB methods, the IB-LB coupling
solver does not determine the complicated time- scheme might be promising for simulating FSI and
10 Proc IMechE Part G: J Aerospace Engineering 0(0)

Table 1. Comparison of CFD methods.

Methods Advantages Disadvantages

VEM Vortex filaments method Computationally efficient; Limit of the application condition;
Vortex particles method stable; Neglecting the flexibility of canopy;
VIC method Pre-maturation;
Conventional grid-based Commercial CFD codes Readily available Rigid-canopy assumption
CFD methods ALE algorithms Good applicability Poor stability and time-precision
DSD/SST algorithms High precision Computationally expensive
IB method Non-deforming grid Lack of the generality

Figure 7. A plain weave structure (cross-sectional image).89 Figure 8. A multi-element canopy model of six plates, sym-
metric wedge.91

moving boundary problems. Some preliminary but should be favoured. Parachute designers prefer to
successful attempts have been conducted83–88 and value the macroscopic structural behaviour of
the application of IB-LB coupling scheme on the canopy and to select sufficiently strong materials
parachute is also anticipated. that minimize the total canopy weight. Therefore, an
accurate macro-scale model may provide an excellent
method for creating parachute structural dynamic
Structure modeling of canopy fabric simulation scenarios.
The parachute canopy is composed of textile materials
such as cloth, nylon, Dacron or other fibres.
Simplified model of a parachute canopy
Generally, these fibres intersect to act as warp and
weft yarns,89 as shown in Figure 7. Several studies An accurate structural model of a parachute canopy
have attempted to determine the mechanism of has long been pursued by many researchers to pre-
fabric structures, and many different models have cisely predict the dynamic behaviour of the canopy.
also been constructed corresponding to the different The first attempt was a simplified 2D parachute
modelling philosophies. A broad classification of model,16 and whose canopy is idealized as a set of
these models corresponds to micro-scale, meso-scale massless rigid plate elements hinged to one another,
and macro-scale.90 The micro-scale model facilitates as shown in Figure 8. The motion of these plate elem-
yarn and fabric modelling from their basic constitu- ents was governed by the dynamic equilibrium condi-
ents and is likely highly computationally intensive for tion that the total moment about any of these hinges
predicting fabric properties. The meso-scale model encompassing unsteady aerodynamic pressure on the
considers the correlative structural properties of plate structure and the rigging line tension must equal
yarn, their interaction with each other and the mech- zero. Unfortunately, the details of their work appear
anical properties associated with the fabric. The vague and require elaboration, though their results
macro scale model focuses on macroscopic mechan- were fairly consistent across the conducted experi-
ical properties, starting with the properties of a por- ments. Subsequently, Robert91 developed a shell
tion of a unit cell and establishing laws for type parachute model involving a parabolic shell of
extrapolating the overall properties of the fabric. revolution with a constant length that could replace a
During the inflation process (when the peak force symmetric canopy without an apex vent. Geometric
often appears), the canopy undergoes a large displace- constraints and the zero total moment criterions
ment and an arbitrary deformation in an extremely about the apex were applied to govern the motion
brief time, when the aerodynamic loads should be sus- of an impervious canopy. Noticeable discrepancies
tained, and proper sufficiently strong fabric materials between the assumptions of rigid multi-plates or
Cao et al. 11

fixed length shell and the characteristics of a real

canopy make their methods counterintuitive and
useful only for qualitative analysis.
The stress applied to a canopy depends on the shape
of an inflated parachute, and this shape relies in turn on
the applied stress. Given this interactional relationship
between the stress and shape, the two variables can be
solved simultaneously. During the inflation process,
the filling time is divided into several time steps. In
each tiny time step, the equilibrium condition of
forces is pursued: the summation of the aerodynamic
forces acting on the canopy is reacted with the riser
force, and the aerodynamic force on each canopy
unit is reacted with forces from the adjacent units.
The D’Alembert principle is then applied to solve a
sequence of equilibrium equations, as is customary in
structural analysis for solving a static problem. Based
on the described procedure, a pressure-stress equilib-
rium method can be devised by studying the shape and
stress distribution of an inflated parachute. This is how
the CANO92 and CALA,93 famous codes for the struc-
tural analysis of parachutes, works.
However, both the CANO and CALA models neg-
lect the inertial term of a parachute in hypothesizing
the quasi-static condition. Generally, the aerodynamic
or the structural dynamic behaviour of a parachute
opening problem cannot be decoupled to consider the
time-variant characteristic of a parachute. Increasing
complexity and difficulty in the coupling aerodynam- Figure 9. (a) Mass spring damper model; (b) force analysis of
ics to structure model prompted several attempts to a user-defined node66
develop an appropriate structural model to represent
parachute inflation. Noteworthy, among these An improved multi-node model94 was developed
attempts is the representative mass-spring-damper based on the same modelling philosophy as a MSD
(MSD) model.66 model by considering the latitudinal direction effects.
The MSD model is established by linking a series In the multi-node model, a force caused by the latitu-
of lumped mass nodes by the parallel connection of dinal strain of a canopy is added to the force analysis
springs and dampers, as shown in Figure 9. The force of each node, and the canopy is simplified into a series
analysis is performed to specify four major force con- of discrete elastic strings satisfying generalized
tributions experienced by each user-defined node, Hooke’s law. Then, a set of first-order ODEs that
including a normal contribution from pressure differ- are nonlinear in space and first order in time is derived
entials on the canopy, two meridional spring and from the governing equations of node motion using
damper contributions from the fabric stretching, and Newton’s Second Law. An explicit Runge-Kutta tech-
a gravitational contribution. Then, Newton’s Second nique is used to iteratively solve the multi-node sub-
Law is applied to obtain a set of nonlinear differential routine. The input parameters and output variables of
equations to govern the motion of nodes. The MSD the multi-node mode are identical with those in the
subroutines require comprehensible input such as a MSD model. The results of a numerical experiment of
pressure distribution along the meridional length of a conical parachute implied that the multi-node model
the canopy and a time step. The subroutines return considerably practicable for predicting the behaviour
the position, velocity and acceleration of each interior of the canopy during the inflation process based on
node as output to the CFD program, which dynam- the approximations used in this model. However, this
ically realizes the coupling process between aero- method is still over-simplified.
dynamics and structural model. A numerical Recently, a three-dimensional spring system model
experiment of a half scale C-9 solid cloth canopy was proposed.95 In this model, the canopy fabric is
was conducted to show that the change in the discretized by a homogeneously triangulated mesh
canopy shape and the payload force vs. time curve and modelled as a series of lumped mass nodes
predicted using the MSD model fit the experimental located at the endpoints of every triangle (shown in
predictions well. Notably, however, the MSD model is Figure 10). The stress tensor on the plane is realized
based on the axisymmetric assumption, and certain by a set of springs connecting to the neighbouring
three-dimensional factors are not considered. points. The spring system supplies restoring forces
12 Proc IMechE Part G: J Aerospace Engineering 0(0)

membrane structure of the canopy, is anchored in

classical elasticity. For tractability, mostly frequently
used assumptions are adopted: the fabric material of
the canopy is treated as a linear elastic continuum
material with no bending stiffness, and the ratio of
thickness to the smallest radius of curvature is signifi-
cantly less than 1. The constitutive equation in the
deformed configuration (i.e. the current configur-
ation) accounting for the kinematic nonlinearities
must be rewritten for the undeformed configuration
(i.e. the reference configuration) based on curvilinear
coordinates (embedded in the reference surface) to
link the second Piola-Kirchhoff stress tensor S with
the Green-Lagrange strain tensor E. (The interested
Figure 10. Fabric surface modelled by Kim et al.95 reader can consult Ref.97 for a full description of the
transformation of the stress and strain tensors from
the current configuration to the reference configur-
against stretching and compression. The similar geo- ation.) The membrane element for a canopy in a
metric shape and good agreement of certain relevant state of planar stress can be expressed as follows
results with experimental results seem to indicate the
effectiveness of the spring model. This approximation Sij ¼ Cij E ð8Þ
cannot precisely predict the structure behaviour
because a spring system cannot reflect the real struc- For isotropic material, C is the forth-order elasti-
tural properties of fabric accurately, nor could the city tensor (associated with the material property)
stress tensor from bending and twisting be expressed. such that

Membrane finite element model Cij ¼ lm Gij G þ m ðGi Gj þ Gi Gj Þ

Researchers have been dedicated to finding better where lm and m denote the Lame parameters and
solutions for highly flexible fabric membrane model- constants associated, respectively, with the young
ling, and the technology of the finite element method modulus and Poisson coefficient of the material; Gij
has turned these solutions into reality. The term designates the components of the contravariant
‘‘membrane’’ may signify different meanings for dif- metric tensor; and the superscripts ij assume the
ferent practitioners. For example, it could represent a values 1 and 2.
transport medium for a chemist, or an idealized model Correspondingly, the cable element for the suspen-
of a shell structure for a structural engineer. For a sion line in a state of uniaxial stress can be expressed as
parachute designer, the membrane is a surface (of
small thickness) with zero or near-zero bending stiff- S11 ¼ Ec G11 G11 E11 ð9Þ
ness to accord with the no-compressive characteristic
of a flexible fabric. Due to the nonlinear trait of the where Ec is the young modulus of the cable material.
canopy, the degree of nonlinearity should be incorpo- For anisotropic materials to manufacture the
rated into the membrane field equations based on the canopy, the angles between the x-axis and the two
formulation philosophy. Theories ranging from small fibre directions (indicated by a dashed line) are 
strain/finite rotation (geometrically nonlinear only) to and , respectively, as shown in Figure 11. The
fully nonlinear in terms of both geometry and mater- Young modulus of fibres is oriented in the first and
ial are being developed by many highly qualified second directions as E1 and E2, respectively.
investigators. Assuming that both sets of fibres undergo deform-
The majority of research on the membrane problem ation in parallel, and the shear effects among the
has been investigated for static loading configurations fibres is negligible, the material properties C can be
or loading regimes with the quasi-static assumption for written in the following form
analysing the governing field equations. Currently, the
mature nonlinear theories for the membrane model C ¼ C1 ðÞE1 þ C2 ð ÞE2 ð10Þ
include the Foppl-von Karman theory for a small
strain with moderate rotation, the continuum mech- where C1 and C2 refer to the material property
anics theory for a large deformation and energy/vari- matrixes associated with the fibre braided angles 
ation principles. These theories are well documented and , respectively. Noticeably, the material property
and detailed in the existing literature.96 of the intersecting fibres depends on the fibre orienta-
The equation formulation of the combined initial- tion. Details of the fibre orientation are specified in
boundary value, which models the response of the Zhang et al.98
Cao et al. 13

where Pcur and Pref denote the external force in the

covariant coordinate system of the current configur-
ation and in the Cartesian coordinate system of refer-
ence configuration, respectively; g and ei are the
covariant base vector and the normal base vector,
respectively; and S and S0 denote the small region
of the current configuration and the reference config-
uration, respectively.
The corresponding virtual external work in the ref-
erence configuration can be expressed as
ZZ
ext
@W ¼ Pref @ui dS0 ð13Þ
@0

where @ui is the virtual displacement in the i direction.

According to the principle of virtual work, the gov-
erning equation (2) (neglecting the damping term for
the linear elastic continuum material) can be rewritten
as
ZZZ ZZZ
Figure 11. The weave pattern of the anisotropic membrane u€ i @ui d0 þ @Eij Sij d0
network.
0 0
ZZ ð14Þ
Coated fabrics, which are composed of certain  Pref @ui dS0 ¼ 0
membrane materials such as polymers or composites,
@0
are also widely used in manufacturing the canopy.
Hence, the viscoelastic effects of these materials
should be considered for improved accuracy. A linear and further simplified as
viscoelastic relation could be sufficient to model the
membrane materials with only a small strain. MU€ þ R  F ¼ 0 ð15Þ
Viscoelastic media can be considered a combination
of the elastic element and the viscous element on the where M represents the consistent mass matrix; R rep-
one hand and the spring and the dashpot for the mech- resents the internal stress vector; and F represents
anical model on the other hand. The standard linear external resultant force vector. For the above equa-
solid model, which is a mixture of the famous Kelvin tions, the nonlinearities of the displacement field make
viscoelastic model and the Maxwell viscoelastic model, impossible closed form solutions for the governing
is generally exploited to model the stress–strain rela- equation (15). Hence, numerical methods can be
tions of viscoelastic material. Details of viscoelastic used to solve the partial difference equations. In
relations are given in René de and Crisfield.99 developing the mathematical formulation of the non-
Finally, the virtual internal work is linear dynamic problem, a tangent stiffness matrix K
ZZZ should be added to the governing equations to char-
int acterize the dynamic material properties. An artificial
@W ¼ @Eij Sij d0 ð11Þ
structural damping matrix C should also be used
0 towards determining the steady solution.
Considering the displacement field Ut at time t and
Importantly, the pressure is always normal to the Utþt at time t þ t, the standard governing equa-
current configuration. The aerodynamic pressure tions can be written as
changes the magnitude and direction of the configur-
ation that ensues from the dynamic change in the MU€ tþt þ CU_ tþt þ Ktþt  U ¼ Ft  Rt ð16Þ
parachute shape and velocity and belongs to the cat-
egory of non-conservative loads. A transformation of where K ¼ @ ðRFÞ
@U , C ¼ cM M þ cK K, U ¼ U
tþt
 Ut
the external force from the current configuration to The governing equation (16) is typically solved by a
the reference configuration must be realized to match spatial discretization with the versatile finite element
the internal force. The transformational relation can method and a temporal discretization with either an
be written as explicit or an implicit finite difference method.100 For
spatial discretization, with appropriate function
dS spaces for different methods, diverse element types
Pref ¼ Pcur ð e  g Þ ð12Þ
*i * dS0 (e.g., a triangular space or a quadrilateral element)
would be created. For temporal discretization, in the
14 Proc IMechE Part G: J Aerospace Engineering 0(0)

case of explicit methods, the current state is computed

based on the forward difference schemes. Hence, the
stability of explicit methods depends on the size of the
time step: a finer spatial mesh thus requires a smaller
time step. Furthermore, provided that the mass
matrix is diagonal, the explicit schemes will have an
overwhelming advantage because the matrix inversion
is no longer required. In the case of implicit methods,
simultaneous equations are required to be solved at
each time step, and the matrix inversion is involved,
which is typically computationally expensive.
However, the implicit schemes are unconditionally
stable, and the relatively larger time step is allowed. Figure 12. Belystchko-Lin-Tsay shell 4-element geometry.102
By implementing the preceding theory, Accorsi and
Leonard presented an SD code, called TENSION,
capable of predicting the nonlinear dynamic response
of highly flexible structures consisting of membranes combining with the local nodal forces and moments,
and cables. In their method, a structured or an the nodal forces-stresses relations are acquired by
unstructured grid with low-order or high-order inter- applying the principle of virtue work (for further
polation functions is employed to perform spatial dis- detail, consult Wang et al.102). Many numerical
cretization. For the time marching method, the SD experiments67,68,102 were conducted to validate that
model includes a central difference scheme (an explicit the shell theory is considerably effective. NASTRAN
method), a Newmark scheme (an implicit method) and ABAQUS also provide similar shell formulation
and a modified Newmark scheme (Hilber-Hughes- for membrane analysis.
Taylor, also an implicit method). The code is
exploited by coupling with the aforementioned
DSD/SST algorithm to simulate the dynamic behav-
Fluid–structure interaction methodology
iour of the parachute canopy and the suspension lines The parachute inflation process is an issue of large-
during the inflation and terminal descent processes by displacement fluid–structure interaction. Modelling
T*AFSM. Numerous practical applications demon- large-displacement fluid–structure problems typically
strate the outstanding serviceability of this SD requires finding a solution to the Navier–Stokes equa-
model in predicting the structural behaviour of tions in a moving domain by coupling with the equa-
highly flexible membranes and cables, particularly tions of large-displacement structural motion. The
an unstructured grid with an HHT scheme for solving development of robust and efficient solution tech-
the governing equations. The interested reader should niques for such problems presents one of the greatest
consult Accorsi and Leonard101 for further detail. challenges in computational mechanics.
Moreover, many commercial codes prefer shell for- Computational methods for solving FSI problems
mulation to model the solid domain occupied by fab- can be classified as monolithic coupled solutions (or
rics. Consider, for example, the Belystchko-Lin-Tsay fully coupled solutions) and staggered coupled solu-
shell formulation (used in LS-DYNA for membrane tions (or partitioned coupled solutions).
analysis). The Belystchko-Lin-Tsay shell 4-node elem-
ent (shown in Figure 12) is based on a co-rotational
coordinate system and a constitutive computation
Coupling method
using a deformation rate. The embedded curvilinear In a monolithic approach, the governing equations
local coordinates that deform with the element are are recast by integrating the fluid and structure equa-
defined according to the four corner nodes. The tions of motion and then solved simultaneously. The
angle between the fibre direction and the unit new single governing equations can be expressed as
normal to the local coordinate experiences changes
 " #  
as the element deforms. The magnitude of the angle Aff Afs Xkf Bf
is limited to maintaining the plane shell geometry. ¼ ð17Þ
Asf Ass Xfs Bs
With the current coordinates, the velocity of a given
point in the shell is determined with the known vel-
ocity of the mid-surface and the rotation of the elem- where the subscripts f and s represent fluid and solid,
ent according to Reissner–Mindlin theory. Then, the respectively; k is the iterative time step; and A, X
deformation rates are computed at the centre of the and B denote the system matrix, the solution matrix
element. Integrating with the material model, the new and the external force matrix, respectively. This equa-
Cauchy stresses are obtained in every element. The tion indicates that determining the fluid–structure
stresses then are integrated to act as an internal interaction relation is crucial. For simple and small-
force and a moment at the element centre. Then, scale structural problems (e.g., involving one or two
Cao et al. 15

degrees of freedom), a fluid–structure interaction

matrix can definitely be written.103 For typically
highly nonlinear fluid problems and structural prob-
lems with multi-degrees of freedom, however, the
interaction relation is highly complex. Commonly,
the fluid and solid equations are discretized and
solved simultaneously,73–77,104 yielding a large
number of coupled nonlinear algebraic equations.
Consequently, the solution of such problems by a
monolithic algorithm is computationally challenging,
mathematically and economically suboptimal and
unmanageable in terms of software.
Alternatively, the FSI problems can be solved with
a staggered procedure. A staggered approach solves
the structure and the fluid by two separate solvers and
typically obtains the coupled solution by exchanging
information only at synchronization points. The stag- Figure 13. Schematic diagram of the CSS (left) and CPS
gered couple method is easy to understand and con- (right) procedures.105
venient to use because the specialized solvers for the
fluid and structural equations are readily available.
An elementary yet frequently-used staggered proced- (structure displacement possesses the same import-
ure for solving FSI problems goes as follows: First, an ance as does the fluid) makes the alternative unrealis-
initial guess for the wall shape is assumed, and the tic, as noted earlier; and (2) the computational cost
fluid equations are solved; then, according to the involved in devising closely coupled procedures coun-
dynamic conditions mentioned in ‘The transmission teracts the advantages of the relatively larger time-
conditions’ section, the traction that fluid exerts on step they may afford.
the wall is evaluated; after solving the structural equa- The CSS procedure allows intra-field parallelism
tions with the new traction, new wall shape is updated with any computational fluid dynamics or computa-
which is taken as the initial boundary for the next tional structural mechanics algorithm. However, the
time step. Clearly, such a partitioned sub-cycle is principle of the CSS procedure prohibits inter-field
known as a loosely coupled solution algorithm. parallelism: these two systems cannot be advanced
However, for improved accuracy or robustness, an simultaneously. To reduce the calculating time,
additional corrector iteration that can be incorpo- advancing the fluid and structural solvers simultan-
rated into each sub-cycle transforms the loosely eously and in a loosely coupled way is expected.
coupled procedure into a strongly/closed coupled The modified partitioned procedure for implement-
solution procedure.105 ing inter-field parallelism is designed as shown in
The basic staggered algorithm is referred to as a Figure 13. As the schematic diagram shows, the
conventional serial staggered (CSS) procedure, a fluid and structure paths are computed in parallel
type of sequential solution method, as graphically along the timeline. The inter-field I/O transfer is exe-
depicted in Figure 13. The simplicity and rationality cuted only at the beginning of each time interval. The
of the CSS method are considerably attractive and simple parallel staggered algorithm is referred to as
make it the most popular procedure among parti- the conventional parallel staggered (CPS) procedure.
tioned procedures developed for time-domain FSI Compared to the CSS method, the CPS method
computation. In the parachute inflation problem, requires relatively small time-steps, while ensuring
the fluid flow typically requires a higher temporal higher numerical stability and accuracy. Some
resolution than does the structural computation; improved parallel versions are outlined in Farhat
thus, for guaranteeing a certain accuracy in the flow et al.107 Additionally, thorough mathematical analysis
solution, the time step involved in overall computa- of different schemes, whether explicit or implicit, is
tion should treat the smallest time step used in the discussed in Fernández.108 Currently, most numerical
flow solution as a benchmark. The CSS procedure is simulations of parachute inflation are performed with
only first-order time-accurate, even when the under- the partitioned loosely coupled method, either serial
lying separate solver is second-order time-accurate, as or parallel.
mathematically demonstrated in Piperno et al.106
Many improved alternatives correct the low accuracy Stability of coupled fluid–structure
of the CSS procedure, yet some of them, such
interaction computations
as monolithic schemes and fully sub-iterations
(closely coupled methods), are restricted in the follow- As explained, the partitioned approach is attractive
ing two ways: (1) the difficulty in developing mono- because software modularity permits the selection of
lithic schemes for a highly complex FSI problem an appropriate solver among the well-established
16 Proc IMechE Part G: J Aerospace Engineering 0(0)

solvers for each domain. Nevertheless, its efficiency is

inferior to its monolithic counterparts because of an
‘‘added-mass effect,’’ a source of numerical instabil-
ities commonly inherent in problems involving large
deformation and lightweight structure. Therefore, the
study of stability and the convergence analysis of par-
titioned coupling algorithms are extremely important.
This section aims to integrate aspects of research in
this field as a reference to provide sufficient informa-
tion for conducting a realistic simulation. When con-
ducting a FSI calculation, numerical instability often Figure 14. Schematic representation of a simplified
parachute.
occurs under the following conditions:

1. The density ratio of fluid to solid is excessively high; boundary. A closed expression for max of a generic
2. An increase in fluid viscosity prompts an earlier geometry can hardly be found. The procedure pro-
occurrence of instabilities, whereas an increase in posed in Brenner110 can be followed to determine the
structural stiffness of parachute improves stability; range of max , adapted to a single subdomain problem.
3. Different temporal discretization schemes influ- A detailed stability analysis of different algorithms
ence the stability of the FSI solver; of FSI problems is presented in Causin et al.109 They
4. The time step size influences the stability of the found that the explicit algorithms were uncondition-
FSI calculation. ally unstable if

In Chapter 2, we discussed fluid governing equa- S hS

51 ð19Þ
tions, the structure governing equations and the fun- f max
damental FSI coupling conditions. To investigate the
stability of the FSI calculation, it is necessary to first Indeed, the instability condition of the explicit
provide an effective FSI governing equation. For our algorithms becomes increasingly restrictive as
mathematical analysis, we employ a simplified 2D S hS =f increases and max decreases. Namely, the
canopy of thickness hS , density S , Young’s modulus numerical instabilities may occur when the density
E, and Poisson’s coefficient . In this case, the inter- of the structure is lower than a certain threshold for
action of the air and the canopy occurs at the a given geometrical configuration or when the length
common boundary , and the motion of the canopy of the domain exceeds a certain threshold for a given
is allowed only in the radial direction. The schematic structure. They also found that even when
representation of the simplified model is shown in S hS =f max 4 1, the instabilities may occur if the
Figure 14, where the dashed lines represent the axis Young’s modulus of the structure is relatively large.
of symmetry for the parachute. In such cases, decreasing time step appropriately can
We refer to the governing equation of FSI in stabilize the scheme.
Causin et al.109 as follows In practice, the explicit algorithms might not be
effective for FSI problems under certain conditions.
@2 xS @2 xfS One option is to use an implicit scheme instead. Here,
S hS þ  M
f a þ a0 xS ¼ Pext
S
ð18Þ
we introduce a frequently used Dirichlet/Neumann
@t2 @t2
iterative method. The prototype of the D/N implicit
EhS
where a0 ¼ r2 ð1 2 Þ, and Ma denotes the added-mass coupling algorithms can be programmed as follows
operator matrix, which represents the interaction of
fluid with the structure and acts as an extra mass. In 1. Initialize xnþ1
0 ;
particular, the presence of the fluid changes the nat- 2. For k50 until the convergence of xnþ1
k ;
ural vibration frequencies of the structure.
Mathematical verification indicates that the operator (a) Update the fluid interior mesh and solve the fluid:
Ma is compact, self-adjoint and positive. (Note that
because we emphasized analysing the added-mass P~ nþ1 ~nþ1
kþ1 ¼ Fðdkþ1 Þ:
instabilities, some nonlinearity terms of the material
in equation (18) are omitted.) (b) Solve the structure:
Regardless of whether we employ an explicit or an
implicit time-marching scheme to complete the tem- d~nþ1 ~ nþ1
kþ1 ¼ SðPkþ1 Þ:
poral discretization, an estimate of the maximum
eigenvalue of Ma, denoted by max , must be deter- (c) Relaxation:
mined before performing the mathematical analysis.
Note that max is a purely geometric quantity asso- dnþ1 ~nþ1 nþ1
kþ1 ¼ !dkþ1 þ ð1  !Þdk :
ciated with size of the fluid domain and the common
Cao et al. 17

(d) Check the convergence: Although the D/N iterative method is widely used
for solving implicit coupling in a partitioned manner,
Check whether kdnþ1 nþ1
kþ1  dk k 5 tol. this method may result in poor convergence behav-
iour with increasing domain length. This limitation
(3) Proceed to the next time step. prompts the development of new variants, for
instance, based on the use of transpiration tech-
If we choose the Implicit Euler time integration niques,116 reduced order models,117 vector extrapola-
scheme for the fluid and the first order backward tion,118 interface artificial compressibility,119 and
difference time integration scheme for the structure Robin-Neumann coupling.120 Interested readers
(IE-BDF scheme), then the D-N iterative method con- should refer to the corresponding references for add-
verges to the solution to equation (18) stably if and itional detail. Recently, a novel partitioned, loosely
only if coupling scheme was introduced.121 Known as the
kinematically coupled -scheme, this scheme does
2ðS hS þa0 t2 Þ not suffer from the added-mass effect and is uncondi-
05!5 ð20Þ
tionally stable for all parameters in the problem. The
S hS þ f max þ a0 t2
interested reader can consult Canic et al.122 for an in-
Furthermore, if the time step is sufficiently small, depth discussion of this scheme.
then equation (20) becomes

2S hS
Interface data transfer
05!5 ð21Þ
S hS þ f max The realization of the coupling process is performed
by data exchange via coupling nodes. Specifically, the
Finally, when S hS ¼ f max , the relaxation factor coupling is enforced by transferring the velocity and
should be strictly less than 1 to maintain convergence displacements from the structure to the fluid and in
of the iterative process. From equation (20), we can turn, transferring the traction from the fluid to the
see that the material and geometrical properties (such structure (from the viewpoint of the direction of
as the fluid and solid densities, the Young’s modulus transfer, the coupling process of parachute inflation
of the structure, and the maximum eigenvalue of the involves two-way or bi-directional coupling). To
added-mass matrix) influence not only the allowable accomplish the coupling, one must connect two
relaxation factor but also the time step size of the meshes at the shared boundary.69 A situation in
calculation. which the meshes of subsystem form a one-to-one
With regard to the stability, the choice of the corresponding relation is expected (shown in
appropriate time evolution schemes is crucial for Figure 15). This situation requires conformity
large-scale computations due to possible severe between two meshes in the same manner at the
time-step restriction. The contradiction between sta- shared interface (identical number of nodes, identical
bility and the computational cost has been frequently nodal co-ordinates and identical interpolation func-
verified.111–114 The fully implicit schemes are uncon- tions for finite element analysis). When these require-
ditionally stable, while the explicit and semi-explicit ments are satisfied, the two meshes can be easily
ones are subjected to a CFL condition. The nature connected by equating the degrees of freedom of the
of such a CFL condition depends on the fluid space corresponding nodes at the shared boundary.
dimension and on the elastic body codimension. On However, time and effort are required in generating
the contrary, the fully implicit method is probably the mesh to ensure consistency.
too expensive for practical application sometimes if Considering that the two meshes are generated
without emendation. The combination of an adap- independently in different subsystems, the dissimilar
tively chosen time step and a predictor-corrector mesh seems more practical (shown in Figure 15). For
strategy would be a good attempt. Moreover, the the nodal information exchange of dissimilar meshes,
influence of the time integration schemes on the
fluid and solid accelerations used in the FSI calcula-
tion is also investigated.115 The results suggest that
the higher order time integration schemes can estab-
lish a more stringent condition to maintain stability
of the FSI, and the critical relaxation factor is higher
when using a structural time integration scheme that
does not involve numerical damping. Another
important parameter is the time step size. Smaller
time step size results in a more restrictive condition,
and as the time step size approaches zero, the value
of the critical relaxation factor converges to a spe- Figure 15. Pairing of a fluid grid point and a wet structural
cific value. element: matching (left) and dissimilar (right) meshes.
18 Proc IMechE Part G: J Aerospace Engineering 0(0)

the interpolation operation before data transmission following aspects of data transfer: (1) the global con-
is essential. In practice, the fluid mesh must be signifi- servation of certain parameters (e.g., loads and
cantly more refined compared with the structure energy) over the interface; (2) the local conservation
mesh. Assume bf t is the fluid boundary interfacing of certain parameters (e.g., displacement and tem-
with the structure domain, and st is the solid bound- perature) over the interface; (3) efficiency and accur-
ary interfacing with fluid domain (in fact, these two acy; or (4) conformity on the order of solvers.
interfaces tend to be identical). When transferring dis-
placements and velocity data from the solid interface
to the fluid interface, bf
Dynamic mesh technique
t is considered the ‘‘receiving
terminal,’’ whereas st is the ‘‘transmitting terminal.’’ When the fluid domain boundaries undergo a rela-
Thus, the data transfer is realized by a node projec- tively larger motion, moving or likely deforming
tion and a plane interpolation or a numerical approxi- boundaries should be considered to solve the asso-
mation in a certain manner. The fluid forces are ciated fluid equations. Methods for solving such prob-
transferred in an identical way. A bi-directional data lems can be classified as either dynamic mesh methods
transfer between the fluid domain and structure or fixed grid methods.72 There are many reasons to
domain is conducted at least once during a physical favour the fixed-mesh methods, particularly when the
time step. Thus, errors introduced by mathematical interface geometric complexity appears unaffordable
methods will accumulate.123 Without measures, or unmanageable for a dynamic mesh method. Kim
these cumulative errors may influence the stability of et al. successfully simulated the process of parachute
the overall algorithms and add incentives for reducing inflation with a local index coating algorithm using a
the frequency of the data transfer and improving the front-tracking method (a fixed method, see Kim
precision of the interpolation scheme. et al.95 for further detail). However, the accuracy of
The existing literature specifies different a non-moving method in terms of tracking the inter-
approaches to address issues of data transfer between face is inferior to that of a dynamic mesh method,
incompatible meshes, such as nearest neighbour inter- because that the resolution of the boundary layer is
polation, the projection method and the spline inter- independent of the complexity of the interface geom-
polation method. If the two grids nearly match, the etry, but limited by the resolution of the mesh inside.
nearest interpolation method is a good choice.124 Generally, for slight deformations, a deforming
Otherwise, the weighted residual method can be mesh will be appropriate for ensuring precision.
applied, based more specifically on an orthogonal With rational modifications, the deforming capacity
projection from the transmitting terminal to the of certain deforming mesh methods can be enhanced,
receiving terminal.125 The starting point is the conser- and the mesh quality can be preserved even in the case
vation of displacement over the interface, as equation of large deformations. In the case of severe deform-
(22) indicates. A weighted residual method and a ations, a mesh reconstruction is indispensable for
series of mathematical transformation are used to sat- ensuring accuracy. Manual re-meshing may assure
isfy this equation as the quality of mesh, if sufficient patience is observed.
However, when the iteration process is lengthy, the
MF x ¼ MS xS ð22Þ time step selected is relatively smaller to improve
accuracy, and the workload becomes time-consuming.
where MF and MS denote the matrix-form expres- Therefore, it is crucial to design an efficient and rapid
sion of the integral over the interface of the basis automatic mesh generator. In parachute fluid–
function multiplied by the Galerkin-based weighted structure interactions, the overall geometry exhibits
function. Two methods based on orthogonal projec- rigid motions as well as structural distortions; a spe-
tion, Gauss interpolation and the Intersection cial-purpose mesh generator may not afford such
method, are applied to evaluate these integrals (con- movement and shapes. Thus, a mesh-moving tech-
sult Boer et al.126 for additional detail). nique should also be incorporated into the automatic
The spline-based method can also be used to mesh generator for such problems.
exchange data. The representation of the displace- Based on the established models, the existing
ment at both the fluid and structure interfaces dynamic mesh methods can be classified as either a
depends on the different basis functions of interpol- physical model method or a mathematical interpol-
ation. The selected basis functions can be diverse, ation method.129 A mathematical interpolation
such as radial basis functions (RBFs), Multi-quadric method, a direct algebraic re-meshing approach, can
bi-harmonic splines (MQ) and Thin-plate splines be used to construct boundary-conforming grids
(TPS).127 In FSI computations, certain radial basis based on the concept of TransFinite (multivariate)
functions with a compact support tended to be more Interpolation (TFI)30 or other mathematical interpol-
advantageous than other basis functions for data ation methods.131,132 In the mathematical interpol-
transfer.128 ation method, the connection mode between
Any of the named data transfer methods may adjacent nodes is no longer required, and the nodal
necessarily involve a compromise with respect to the position is the only concern. Therefore, the data
Cao et al. 19

structure is relatively simple, and the corresponding

data storage space decreases. However, this mathem-
atical method is only employed to generate structured
grid and is unsuitable for unstructured/hybrid grids
widely used for complex shapes. This method is also
ill-suited for a multi-block domain, which seriously
restricts its universality. Therefore, the mathematical
interpolation method is commonly used together with
physical model methods, thus as a part of hybrid
approach.
The next subsection discusses the feasibility of
physical model methods for large-deformation prob-
lems such as the parachute inflation process. For the
most general case, the dynamic mesh is represented by
a discrete or continuous pseudo-structural system,
corresponding to a discrete or continuous model. A
brief introduction of each method is provided with Figure 16. Spring model by Batina.134
reference to theoretical methods and empirical
research.

Discrete model. The behaviour of such a mechanical and thereby improve the robustness for enlarging the
system is governed by the semi-discrete equations method’s range of application.133 Unfortunately, this
for dynamic equilibrium133 improved method requires a kinematic formulation
and conversions of forces and moments based on a
Mx€ þ Dx_ þ Kx ¼ gðxS Þ on  ð23Þ predefined condition of sufficiently small deform-
ations and rotations. When applying this model in
where x is the displacement vector of the aforemen- three-dimensional situation, rather complex formula-
tioned grid point; xS is the displacement vector on the tions and appreciable computational overheads are
moving boundary determined by simultaneously sol- moreover involved. Zeng and Ethier135 proposed a
ving equations (1) and (2); M, D, and K denote the semi-torsional spring model to avoid misshapen elem-
fictitious mass, damping and stiffness matrices, ents and modified the model for application to realis-
respectively; and the mesh motion triggered by the tic 3D problems. Extensive numerical experiments
displacement of the structure is represented by gðxS Þ. showed that the algorithms maintain the mesh quality
For most research on aerodynamic meshes docu- even under severe deformations. Based on those
mented in the aerodynamic computation literature, results, the modified semi-torsional spring model can
regardless of whether continuous or discrete methods be applied for solving the parachute FSI process.
are employed, the quasi-static version of the pseudo
structural system is of particular concern. Thus, equa- Continuous model. The equations governing the dis-
tion (23) becomes placement of internal nodes can also be written in a
continuous form,133 namely as
Kx ¼ 0 ð24Þ
@2 x
The grid points attached to the boundary possess   divðE : "ðxÞÞ ¼ gðxS Þ on  ð25Þ
@t2
an identical location and velocity to the solid surface
to satisfy the continuity condition of kinematics. where E is the fictitious stress tensor. To maintain
Therefore isotropic linear elasticity, E is defined as

xðtÞ ¼ xS ðtÞ on St E ¼ ltrð"ðxÞÞI þ 2"ðxÞ ð26Þ

@x @xS
¼ on St where
@t @t
In a discrete approach, the tension/compression 1 
"ðxÞ ¼ rx þ ðrxÞT
spring (shown in Figure 16) is often adopted as an 2
analogy to construct the left-hand side of equation E
(24).134 However, this method was limited to small ¼
2ð1 þ Þ
deformations. When used for large deformations, it
E
is of ineffective, and the adverse crossover of gridlines l¼
appears. An additional torsional spring can be used to ð1 þ Þð1  2Þ
inhibit the interpenetration of neighbouring triangles
20 Proc IMechE Part G: J Aerospace Engineering 0(0)

where, trðÞ is the trace operator, and I denotes the solid surface and assigned higher rigidity in the gov-
identity tensor. erning equations of grid motion than other fluid elem-
The boundary  is composed of F and S , which ents. The two following different methods of solving
are the far field boundary and the moving wall bound- the inconsistent equations are recommended: a separ-
ary, respectively, and can be considered the ate solution option known as ‘‘SEMMT-Multiple
Neumann-type and Dirichlet boundary conditions. Domain’’ and a uniform solution option known as
The Neumann-type and Dirichlet boundary condi- ‘‘SEMMT-Single Domain’’ (for additional detail,
tions are represented as consult Tezduyar137).

n   ¼ Ff on F Comparison among different coupling strategies

x ¼ xS on S using a two-dimensional parafoil
Numerical calculations on a two-dimensional parafoil
The fictitious young modulus and the Poisson coef- with a cut at 10% chord length from the leading edge
ficient of materials determined the element stiffness and at an angle of 45 with respect to the free-stream
matrix and the compressibility of grids.136 direction are carried out to investigate the advantages
Generally, the boundary layers of solid surfaces and disadvantages of different coupling schemes, such
must be fully controlled and require higher reso- as the explicit Drichlet/Neumann (EDN) method,
lutions for accurately measuring the aerodynamic implicit Drichlet/Neumann (IDN) method and semi-
forces involved. Empirically, the mesh size of the implicit Robin-Neumann (SIRN) method.138 The sche-
first layer in the boundary layer is usually set to be matic diagram of parafoil is illustrated in Figure 18 and
of the fifth order of the characteristic length. the detailed parameters are shown in Table 2.
Consequently, the smaller elements tend to exhibit Here, the semi-implicit method for pressure-linked
significantly larger distortions even if the boundary equations (SIMPLE) algorithm with a k  " two-
experiences tiny motions, thus leading to poor mesh equation turbulence model is employed in solving
quality. To overcome it, selective treatment based on the flow field around and inside the parafoil.
element size is made by altering the conversion of
Jacobian from the element domain to the physical
domain. This Jacobian-based stiffening (JBS)
method renders the smaller elements ‘‘harder’’ by
introducing a stiffening power (the effect is shown in
Figure 17). Additionally, to address body-fitted grids,
a solid-extension mesh moving technique (SEMMT)
is proposed (see Figure 17). This method tends to
treat thin fluid elements as an extension of solid elem-
ents, and the thin elements are ‘‘pinned’’ onto the Figure 18. Schematic diagram of a 2D parafoil section.138

Figure 17. (a) Deformed mesh with a stiffening power of 0.0, 1.0, and 2.0 (first row).119; (b) deformed meshes for the SEMMT Single
Domain (second row).137
Cao et al. 21

Table 2. Geometric parameter of parafoil section. Table 3. Computational parameter of fluid and canopy.
?
Parameter/unit c/m h/m
/ Airfoil Parameter Value

Value 1.00 0.1 135 NASA LS1-0417 Density of fluid (kg/m3) 1.225
Velocity (m/s) 10
Angle of attack ( ) 0
Density of structure (kg/m3) 533.8
Computational domain is constructed for the 2D Young’s modulus (Pa) 4.309  108
incompressible simulations by extruding the mesh
Poisson’s ratio 0.14
from the first 105  c high mesh layer around the
Thickness of canopy (cm) 0.1
airfoil, as shown in Figure 19. An inner mesh with
totally 40 layers is firstly extruded with a lower
growth rate of 1.05 to guarantee an accurate capture
of the boundary-layer turbulent flow field informa-
tion. Then a second extrusion of 91-layer mesh is Table 4. Number of iterations for the three schemes with an
made with a higher growth rate of 1.2 to finally optimal choice of the relaxation parameter.
form a sufficiently large domain (15  c in radius) to
EDN SIRN IDN
avoid the effects of the outlet boundary condition but
also to enhance the calculation efficiency. We choose !opt 0.15 1 1
an under-inflated parafoil as the initial model. In this Number of iterations 68 8 82
model, the entire canopy is released except that the Computational time 3 h 46 min 2 h 14 min 6 h 54 min
point A, B and C are omnidirectionally constrained.
The segment AC is set to be horizontal with no pre-
stress, while the segment BC ‘‘sags’’ naturally with its
shape satisfying the catenary equations. Besides, for
the initial configuration, an extremely small gap
Parachute permeability
between segment AC and BC, located in the trailing The air-flow around a parachute is irrevocably asso-
edge, is pre-reserved on the premise of not influencing ciated with separation and vortex shedding. Vortex
the flow regime around and inside the section. The shedding typically occurs and affects the canopy
membrane finite element method is employed in mod- motion asymmetrically, consequently prompting
eling the structural deforming. The structural material parachute instability. Certain ventilated measures
of the canopy is nylon and the relevant physical par- are adopted to suppress the asymmetric wake oscilla-
ameters are exhibited in Table 3. We design a mod- tions at the cost of reducing the aerodynamic drag.139
ified spring-TFI dynamic mesh method to detect the Therefore, a ventilation effect should be considered
grid moving, which is described in Nie et al.138 We for more accurately simulating parachute behaviour.
compare the performance of EDN, IDN and the There are two broad classifications of ‘‘leakages’’
SIRN schemes. In all these tests, we have employed in a parachute system when the parachute performs
an optimal relaxation factor and the computational its mission. The first type, venting, is a geometric ven-
time-step is set to be 0.0005 s. Table 4 exhibits the tilation and refers to the geometrical shape of the
optimal relaxation parameter of different schemes, parachute. Generally, the parachute canopy is
the corresponding number of iterations and computa- equipped with a vent at its apex, which helps to con-
tional time. Figures 20 and 21 plot the normalized trol the inflation process and to stabilize the near
residual versus the number of iterations and the dis- wake region. For some featured-purpose parachutes,
placement of the upper surface along the chord direc- different types of gaps are typically designed between
tion, respectively. The results highlight the SIRN ribbons or ringslots on their structural configurations,
scheme for its excellent efficiency without compromis- thus suppressing the large-scale vortex motion behind
ing the computational accuracy. the canopy, such as a ring-sail parachute, a disc-gap-
band parachute, and a ring slot parachute. The
Physical phenomena accompanying the numerical method solves these large ‘‘holes’’ through
the geometry model construction.140 The second type,
inflation process porosity, attaching to material permeability, stems
The main framework of the FSI problems involving from the porous nature of the canopy fabric.
the parachute opening process was emphasized in the The ‘‘leakage’’ of this section is generally evaluated
above-mentioned discussion. Next, we focus on three with experiments or a mathematical model. Thickly
main physical phenomena, including the embodiment dotted holes distributed on the canopy make it
of porosity in the numerical method for parachute porous (shown in Figure 22). It is computationally
permeability, the reflection of the wrinkling phenom- unmanageable to resolve these holes individually
enon in the membrane structure, and the flexible because doing so requires that the fluid mesh-width
structure impact-contact phenomenon. placed in a pore should be smaller than the diameter
22 Proc IMechE Part G: J Aerospace Engineering 0(0)

of the pore. Only if this requirement is satisfied, might ALE algorithm is employed in assessing a TP8 low
air flow through these holes be captured by the refined attitude troop parachute.142 A qualitative assessment
mesh.81 Rather than tracing each pore individually of the permeability effects is shown in Figure 23. The
(the microscopic scale), quantities of interest are mea- behaviour of the parachute wake changes completely
sured over areas covering many pores, and such while air flow passes through the porous canopy.
space-averaged quantities (the macroscopic scale) Specifically, the large recirculation of air behind the
change in a relatively regular manner and hence are canopy is blown further downstream so that the influ-
amenable to theoretical treatment. ence on the induced instability from vortex-shedding
Assume the canopy is a thin shell of porous mater- weakens.
ial with pores oriented perpendicular to the tangent
plane of their location and exhibiting uniform, iso-
tropic porosity. The term porosity ’ refers to the per-
Wrinkling mechanics of the membrane structure
meability of the fabric and is defined as the fraction of The computational challenge associated with pre-
the total volume of the canopy that is occupied by the dicting the behaviour of the parachute canopy
pores and that only depends on weft/warp density of (membrane) derives from the under-constrained and
the fabric. The diameter of each pore is also assumed no-compression nature of the canopy. The undercon-
to be small compared with the thickness of the canopy strained characteristic causes large displacements
h so that the air flow rushes out of each pore in a without an accompanying strain energy, while no-
direction normal to the surface. The porous coupling compression causes an out-of-plane response, named
force can be obtained with the Ergun equation102 wrinkling (shown in Figure 24). The membrane has

p
¼ að, ’Þvp þ bð, ’Þv2p ð27Þ
h 1
EDN
SIRN
where vp ¼ uðxS Þ  x_ S denotes the penetrating vel- IDN
normalized residual

ocity or seepage velocity;  is the dynamic viscosity;

p is the differential pressure between the inside and 0.01
the outside of the canopy; að, ’Þ is the reciprocal
permeability of the canopy or a viscous coefficient;
and bð, ’Þ is the inertial coefficient. The deterministic
representation and value range of each parameter is
1E-4
noted in Donald and Adrian.141 For parachute appli- 0 30 60 90
cation, the inertial term should be dominant. This number of iterations

correlation can be considered an extension of

Darcy’s law and reflects the coherence of the fluid Figure 20. Number of iterations for the three schemes with
momentum pattern. The Ergun equation with the the counterpart optimal choice of the relaxation parameter.138

Figure 19. (a) Extruded circle fluid domain. (b) Near mesh around the airfoil wall. (c) Partial enlarged drawing of the near mesh
around the trailing edge.138
Cao et al. 23

identically zero bending rigidity, and when the mem- Therefore, constraints should be enforced to prevent
brane sustains uniaxial in-plane tension, it tends to the interpenetration of those contact elements.
respond with an in-plane contraction. If the contrac- Original finite element formulations cannot accom-
tion exceeds Poisson’s effect, the parachute will modate the treatment of contact, so ad hoc modifica-
buckle, and many narrow wrinkles will appear with tions should be added to extend the variational
crests and troughs along the tensile direction.143 formulation based on by the finite element method.
During inflation, the canopy in some regions of the Then, the contact elements are formulated148 as
structure may be in uniaxial tension. Over the past
two decades, analysing the wrinkling phenomenon MU€ þ CU_ þ ðK þ KC ÞU ¼ F  R  RC ð28Þ
in membranes has been extensively and profoundly
studied from computational perspectives. Before we where KC is the tangent stiffness matrix corresponding
develop a mathematical description of this degenerate to the constraining force RC . The constraint condition
state, the concrete judgement criterion should be sum- is imposed on the contact elements with the form of
marized. At any point on the parachute surface, the constraining force. Typically, the active contact con-
canopy must be in one of the following three states: a straints are determined by a penalty function method
‘‘slack’’ state, characterized by an absence of stretch- (normal stiffness, tangential stiffness and displacement
ing in any direction; a ‘‘taut’’ state, characterized by are employed to determine the normal contact force
tension in all directions; or a ‘‘wrinkled’’ state. Several and the tangential contact force), and the estimated
wrinkling criteria were categorized into three types contact forces are determined with Lagrange multipli-
and the mixed criteria (based on both principal stres- ers to enforce contact conditions.
ses and strains) facilitate the most accurate character- We assume the impact between gores of parachutes
ization of a real membrane state144 as pertains to the completely inelastic collision. The con-
If min 4 0 taut tact between gores can be classified into the following
If min 40 and max 4 0, wrinkled three types: the frictionless, static and dynamic fric-
If max 40, slack tional. Hughes et al.149 presented a detailed analysis
Then, three approaches, based on the premise of of impact frictionless problems in which the tangent
membrane theory, aim to develop an enhanced model constraints were ignored. Parisch150 provided consist-
to eliminate the spurious compressive stresses allowed ent tangent stiffness matrices when addressing
by the classical model.145 The first approach tries to quasi-static large deformation. Laursen and Simo151
modify the deformation gradient tensor based on the developed an approach for solving dynamic frictional
wrinkling direction, thereby correcting the stress and contact problems. They employed the penalty and
strain. A similar treatment occurs in Fan and Xia59 augmented Lagrangian regularization methods to
and Lu et al.146 The second approach attempts to create kinematic contact conditions. Notably, the
directly modify the total stiffness matrix to achieve use of the variational formulations to solve the
the absence of compressive stresses (see Liu et al.144 quasi-static or dynamic contact problems is not math-
and Rossi et al.147). The third approach applies the ematically rigorous because the non-differentiability
relaxed energy functional method in adjusting the of Coulomb’s friction law limits properly expressing
structural equations to eliminate the influence of the variational formulations.152 Furthermore, the per-
the compressive stresses (see Gil and Bonet 145 for a formance of those contact elements, on account of
detailed discussion). Validation examples of the accuracy and convergence, depends on the user
former two approaches involving parachute behavior defined parameters. Due to the mathematical incon-
are provided, though alternatives can also be sistences associated with the variational equations,
implemented. much endeavour was made to seek accurate expres-
sion of the contact governing formulation.
Impact-contact phenomenon The exact variational representation of frictional
contact problems produces a variational inequality
of the flexible structure
because the imposed restrictions on admissible dis-
The physical impact-contact phenomenon is another placement, velocity and stress field are not equalities,
problem posed by a highly folded parachute in the but inequalities. The contact governing equations of
incipient state of the inflation process. This situation the general form together with the displacement and
becomes particularly apparent in the case of studying traction boundary conditions are discussed in
the motion of clusters of parachutes. During the infla- Belytschko and Neal,153 which provide a detailed
tion of a parachute or clusters of parachutes, different solution of the penalty and Lagrangian methods.
gores of parachute may collide and subsequently stick Meguid and Czekanski154 outlined the basic steps of
together in the normal direction, a behavior known as deriving the variational inequality formulation for
impact-contact in engineering structures. However, in elasto-plastic material undergoing large displacement
numerical computation, the large displacement of the and strains.
structure may generate nonphysical solutions in which Another important point in addressing impact-
two canopies occupy the same place at the same time. contact problems is determining the elements that
24 Proc IMechE Part G: J Aerospace Engineering 0(0)

should be involved in the contact calculation, i.e. col-

5
EDN lision detection. Commonly, static interference check-
4
SIRN ing is performed repeatedly while the structure is
displacement (mm)

IDN moving, i.e. ‘‘by step check.’’154 Contact detection

3 searchers were developed to monitor whether the
interference among multi-bodies exists. These search-
2 ers typically monitor the geometrical relationship of
the target surface to the contact surface.155 Sathe and
1 Tezduyar developed a surface-edge-node contact
tracking (SENTC) algorithm. A topologically hier-
0 archical search is a major aspect of SENCT (for
0.0 0.2 0.4 0.6 0.8 1.0
x (m) details of this search, consult156). The contacting sur-
faces are separated by a prescribed minimum distance
Figure 21. Displacement of the upper surface along the to guarantee the quality of fluid mesh between them.
chord direction for the three schemes.138 Another equally important aspect of the simulation
of multibody systems with contact-impact is detecting

Figure 24. Wrinkled (solid), pseudo and lengthened (dashed

Figure 22. Porous boundary.81 line) surface.143

Figure 23. Visualization of the permeability effect.142

Cao et al. 25

Figure 25. Flow chart of the computational contact procedure.157

the precise instant of impact. Because the governing

help with parachute design, as follows:
equations of structural motion are described by a set
of differential equations, and the error of response 1. At present, the three-dimensional VIC method has
system is particularly sensitive to constraints viola- been successfully applied in predicting aero-
tion, selecting a time step to integrate the time deriva- dynamics of parachute, while the application of
tives of the state variables is crucial to account for the the vortex filaments method and the vortex par-
dynamics of multibody systems. A conservative ticles method is still in exploration.
approach is to employ smaller time steps throughout 2. Grid-based CFD method has been extensively
the computational process. This approach, however, applied in computing the aerodynamic of para-
results in poor computational efficiency. Flores and chute, and the ALE, DSD/SST and IB method
Ambrosio157 outlined a computational contact pro- exhibit the outstanding advantages in handling
cedure for the dynamic analysis of a constrained mul- with free-surface flow problems.
tibody system (shown in Figure 25) and presented a 3. The membrane finite element model can certainly
general and comprehensive approach to automatically be the future trend in modeling the structure of
adjusting the time step, using variable time-step inte- canopy fabric.
gration algorithms, in the vicinity of contact of multi- 4. Different coupling strategies have advantages and
body systems. Ganter and Isarankura158 proposed a disadvantages in terms of precision, numerical
‘‘space partitioning’’ technique to detect dynamic col- instability, and computational economic.
lisions during computer simulations, and thus dra- 5. The existing dynamic mesh methods can be classi-
matically reducing the required computation time. fied as a physical model method or a mathematical
interpolation method. The former shows a wider
range of the application than the latter, while exhi-
Conclusions and recommendations
bits worse computational efficiency than the latter.
Numerical simulations of the parachute inflation pro- 6. Currently, the space-averaged method is the
cess must account for a complex, modular system common practice in treating the porosity of
accompanied by flow around complex geometries, canopy. The frequently-used model to evaluate
intricate structural mechanics and the interaction the porous coupling force is based on the Ergun
between fluid and structure. Such simulations also equation.
address many other issues, seemingly small yet 7. The accurate criteria to characterize a real mem-
important, including the flow through porous media, brane state is concluded and three approaches to
wrinkling of flexible fabric and impact-contact. mathematical description of the wrinkling phenom-
Tremendous challenges arising from the behaviour enon based upon different principles are outlined.
of flexible fabric in fluids have inspired researchers 8. The mechanism of the impact-contact problems is
worldwide to engage in this field, and great advances clarified and the common practice in handling
have recently been made to use numerical methods to with the impact-contact problems is listed.
26 Proc IMechE Part G: J Aerospace Engineering 0(0)

Nonetheless, the achievements obtained are far involved in parachute inflation are not yet perfect.
from adequate, and further research on certain The elaborate representation of flows through
issues is necessary. To satisfy the requirements of porous media in relation to the permeability of
engineering practice, we propose the following prob- fabric, the membrane response of the wrinkling
lems as worthy of greater attention in future research. phenomenon and the flexible impact-contact prob-
lem between fabrics depends on the development
1. Wake recontact occurs when axial velocities in the of methods or/and models in related fields to
wake region exceed that of the parachute, thus simulate the parachute inflation process.
allowing the wake to overtake the parachute
during inflation. Inherently, the interaction In order to improve the study in this field, future
between the trailing vortex and the canopy leads challenges lie in developing more accurate and efficient
to the parachute collapsing. Previously, the fluid and structural models with the inclusion of the
canopy was considered a rigid wall, regardless of more realistic turbulence effect, and closer scrutiny of
the interaction between the flexible fabric and the the effect of fabric properties, such as stress–strain
vortex. Recent experiments have demonstrated behaviour and anisotropy. One of the increasing inter-
that the vortex core breakdown and the disintegra- ests is gained in the application of the IB-LB coupling
tion cause unsteady, asymmetric deformations on method, which has been successfully employed in other
the canopy surface. Therefore, a more accurate fields. Besides, the more accurate and robust algo-
representation of the flexible wall condition rithms in handling with FSI calculations are also pur-
should be incorporated into the VEM. sued. This will serve as a motivation to develop the
2. The development of parachute inflation models monolithic coupling strategy and design more efficient
has advanced considerably in recent decades. implementations of the fully implicit schemes, e.g. the
These models have been used to aid the design adaptive Newton continuation strategy. Another gain-
of parachutes. Nonetheless, it is necessary to ing interest is the Robin–Neumann coupling method,
develop more refined grids throughout the flow which exhibits excellent stability in FSI calculations.
field, more convenient re-gridding techniques, Furthermore, elaborate treatments of flow through
improved models and methods emerging in the porous media in relation to the permeability of fabric
fluid or/and structure disciplines, and more and the wrinkling of membrane are desired to simulate
mature numerical stability techniques that can the parafoil inflation process for better accuracy in
eliminate artificial oscillation efficiently. practical applications.
3. Existing information exchange techniques are typ-
ically based on mathematical interpolation or Acknowledgement
numerical approximation, in which methods the The authors would like to acknowledge Mr. Junsen at
accumulative error may influence the accuracy Beihang University for technical discussions and Mr.
and stability of the overall numerical algorithms. Shufeng Liu at Beihang University for linguistic editing.
To avoid this influence, practices such as reducing
the frequency of data transfer and improving the Declaration of Conflicting Interests
precision of interpolation schemes should be pur- The author(s) declared no potential conflicts of interest with
sued. The partitioning procedure employed in respect to the research, authorship, and/or publication of
coupled solutions of FSI problems should also be this article.
improved by simplifying explicit/implicit treatment,
sub-cycling, load balancing and software modular- Funding
ity. The most fundamental solutions for such prob- The author(s) received no financial support for the research,
lems can be determined with monolithic algorithms authorship, and/or publication of this article.
if the computational, mathematical and economical
challenges mentioned above are well addressed. References
4. The most profound difficulty in dynamic mesh 1. Way DW, et al. Mars science laboratory: entry, descent,
techniques entails realizing the quick tracking of and landing system performance. IEEE Aerosp Conf
grid points, while maintaining high mesh quality. Proc 2007; 1467: 1–19.
The single dynamic mesh technique barely main- 2. Knacke TW. Parachute recovery systems: design manual.
tains large deformability and high quality simul- Santa Barbara, CA, USA: Para Pub, 1992.
3. Peterson CW, Strickland JH and Higuchi H. The fluid
taneously. Thus, attention should be paid to the
dynamics of parachute inflation. Ann Rev Fluid Mech
hybrid method of physical model and mathemat-
1996; 28: 361–387.
ical interpolation methods with the combination 4. Anderson JD. Computational fluid dynamics: the basics with
of multi-block grids techniques. applications. New York, NY, USA: McGraw-Hill, 1995.
5. Additional details are discussed in Chapter 6 to 5. Wilcox DC. Turbulence modeling for CFD. La Canada,
more accurately simulate parachute behaviour. CA, USA: DCW Industries, 2006, pp.363–367.
Nevertheless, methods for addressing the afore- 6. Hu HH, Patankar NA and Zhu MY. Direct numerical
mentioned three main physical phenomena simulation of fluid-solid systems using the arbitrary
Cao et al. 27

Lagrangian-Euler technique. J Comput Phys 2001; 169: Systems Technology Conference, San Diego, CA,
427–462. USA, 9–11 April 1991.
7. Karimian SMH and Schneider GE. Pressure-based con- 26. Leonard A. Numerical simulation of interacting, three-
trol-volume finite element method for flow at all speeds. dimensional vortex filaments. Lecture Notes Phys 1975;
AIAA J 2012; 33: 1611–1618. 35: 245–250.
8. Jenny P, Lee SH and Tchelep HA. Adaptive multiscale 27. Leonard A. Computing three-dimensional incompress-
finite-volume method for multiphase flow and transport ible flows with vortex elements. Ann Rev Fluid Mech
in porous media. Siam J Multiscale Model Simul 2004; 1985; 17: 523–559.
3: 50–64. 28. Winckelmans G and Leonard A. Improved vortex
9. Hinze M and Kunisch K. Three control methods for methods for three-dimensional flows. Mathematical
time-dependent fluid flow. Flow Turbul Combust 2000; aspects of vortex dynamics. Math Aspects Vortex Dyn
65: 273–298. 1989; 22: 25–35.
10. Tezduyar T, Aliabadi S, Behr M, et al. Flow simulation 29. Leonard A. Vortex methods for flow simulation.
and high performance computing. Comput Mech 1996; J Comput Phys 1980; 37: 289–335.
18: 397–412. 30. Greengard L and Rokhlin V. A fast algorithm for par-
11. Martys NS and Chen H. Simulation of multicomponent ticle simulations. J Comput Phys 1997; 135: 280–292.
fluids in complex three-dimensional geometries by the 31. Salmon JK and Warren M S. Skeletons from the tree-
lattice Boltzmann method. Phys Rev Stat Phys Plasmas code closet. J Comput Phys 1994; 111: 136–155.
Fluids Relat Interdis Topic 1996; 53: 743. 32. Degond P and Mas-Gallic S. The weighted particle
12. Opstal TM Van and Brummelen EHV. A finite-ele- method for convection-diffusion equations. I. The case of
ment/boundary-element method for large-displacement an isotropic viscosity. Math Comput 1989; 188: 485–507.
fluid–structure interaction with potential flow. Comput 33. Cottet GH and Koumoutsakos PD. Vortex methods:
Meth Appl Mech Eng 2013; 266: 57–69. theory and practice. The Edinburgh Building, Cambridge,
13. Desabrais KJ and Johari H. The flow in the near wake UK: Cambridge University Press, 2000, p.354.
of an inflating parachute canopy. In: 16th AIAA aero- 34. Ploumhans P and Winckelmans GS. Vortex methods
dynamic decelerator systems technology conference and for high-resolution simulations of viscous flow past
seminar, Boston, MA, USA, 21–24 May 2001. bluff bodies of general geometry. J Comput Phys
14. Sarpkaya T. Computational methods with vortices-The 2000; 165: 354–406.
1988 Freeman scholar lecture. J Fluid Eng 1988; 111: 35. Ould-Salihi ML, Cottet GH and El Hamraoui M.
5–52. Blending finite difference and vortex methods for
15. Leonard A. Vortex methods for flow simulation. incompressible flow computations. Siam J Sci Comput
J Comput Phys 1980; 37: 289–335. 2000; 22: 1655–1674.
16. Reddy KR and Roberts BW. Inflation of a multi-element 36. Cottet GH. Artificial viscosity models for vortex and
parachute structure. In: 5th AIAA aerodynamic deceler- particle methods. J Comput Phys 1996; 127: 299–308.
ator systems technology conference, AIAA-75-1353, 37. Cottet GH and Poncet P. Advances in direct numerical
Albuquerque, New Mexico, 17–19 November 1975. simulations of 3D wall-bounded flows by Vortex-in-
17. Klimas P C. Inflating parachute canopy differential Cell methods. J Comput Phys 2004; 193: 136–158.
pressures. J Aircraft 1979; 16: 861–862. 38. Cottet GH, Gallizio F, Magni A, et al. A vortex penal-
18. Frucht Y. A discrete free vortex method of analysis for ization method for flows with moving immersed obs-
inviscid axisymmetric flows around parachute canopies. tacle. In: ASME-JSME-KSME joint fluids engineering
In: 11th Aerodynamic Decelerator Systems Technology conference, Hamamatsu, Japan, 24–29 July 2011,
Conference, San Diego, CA, USA, 9–11 April 1991. pp.3703–3708. New York, NY, USA: American
19. Mccoy H and Werme T. Axisymmetric vortex lattice Society of Mechanical Engineers.
method applied to parachute shapes. In: 9th 39. Rasmussen JT, Cottet GH and Walther JH. A multi-
Aerodynamic Decelerator and Balloon Technology resolution remeshed Vortex-In-Cell algorithm using
Conference, Albuquerque, NM, USA, 7–9 October patches. J Comput Phys 2011; 230: 6742–6755.
1986. 40. Summers DM. Complementary modes of impulse gen-
20. Chorin AJ. Numerical study of slightly viscous flow. eration for flow past a three-dimensional obstacle.
J Fluid Mech 1973; 57: 785–796. Theor Comput Fluid Dyn 2007; 21: 15–38.
21. Chorin AJ. Vortex sheet approximation of boundary 41. Zuhal LR, Dung DV, Sepnov AJ, et al. Core spreading
layers. J Comput Phys 1978; 27: 428–442. vortex method for simulating 3D flow around bluff
22. Cheer AY. Numerical study of incompressible slightly bodies. J Eng Tech Sci 2014; 46: 436–454.
viscous flow past blunt bodies and airfoils. SIAM J Sci 42. Yoo-Chul K, Suh JC and Kyung-Jun L. Vortex-in-cell
Stat Comput 1983; 4: 685–705. method combined with a boundary element method for
23. Strickland JH. Prediction method for unsteady axisym- incompressible viscous flow analysis. Int J Num Meth
metric flow over parachutes. J Aircraft 1994; 31: Fluids 2012; 69: 1567–1583.
637–643. 43. Pelegrini MF and Vieira EDR. Flow past a sphere at
24. Behr V, Wolfe W, Peterson C, et al. An overview of the moderate Reynolds numbers. In: 17th international con-
development of a vortex based inflation code for para- gress of mechanical engineering, COBEM, Sao Paulo,
chute simulation (VIPAR). In: 15th AIAA aerodynamic Spanish, 10–14 November 2003.
decelerator systems technology conference, Toulouse, 44. Shirayama S and Kuwahara K. Computation of flow
France, 8–11 June 1999. past a parachute by a three-dimensional vortex method.
25. Sarpkaya T. Methods of analysis for flow around para- In: 24th Aerospace Sciences Meeting, Reno, NV, USA,
chute canopies. In: 11th Aerodynamic Decelerator 6–9 January 1986.
28 Proc IMechE Part G: J Aerospace Engineering 0(0)

45. Strickland J. On the development of a gridless inflation Meeting & Exhibit, Reno, NV, USA, 12–15 January
code for parachute simulations. Office of Scientific & 1998.
Technical Information Technical Reports, 2000. 61. Luo H, Baum JD and Lohner R. Extension of Harten-
46. Aparinov AA, Setukha AV and Zhelannikov AI. Lax-van Leer scheme for flows at all speeds. AIAA J
Numerical simulation of separated flow over three- 2005; 43: 1160–1166.
dimensional complex shape bodies with some vortex 62. Weiss JM and Smith WA. Preconditioning applied to
method. Am Inst Phys Conf Ser 2014; 1629: 69–76. variable and constant density time-accurate flows on
47. Aparinov AA and Setukha AV. Application of mosaic- unstructured meshes. In: 25th AIAA Fluid Dynamics
skeleton approximations in the simulation of three- Conference, Colorado Springs, CO, USA, 20–23 June
dimensional vortex flows by vortex segments. Comput 1994.
Math Math Phys 2010; 50: 937–948. 63. Hirt CW, Amsden AA and Cook JL. An arbitrary
48. Aparinov AA and Setukha AV. Parallelization in the Lagrangian-Eulerian computing method for all flow
vortex method for solving aerodynamic problems. speeds. J Comput Phys 1974; 14: 227–253.
Moscow, Russia: [Link] Moscow State 64. Hirt CW, Amsden AA and Cook JL. An arbitrary
University, 2013, pp.406–418. Lagrangian-Eulerian computing method for all flow
49. Nelsen J. Computational fluid dynamic studies of a speeds. J Comput Phys 1974; 14: 227–253.
solid and ribbon 12-gore parachute canopy in subsonic 65. Amsden AA, Ruppel HM and Hire CW. SALE: a sim-
and supersonic flow. Office of Scientific & Technical plified ALE computer program for fluid flow at all
Information Technical Reports, 1995. speeds. Los Alamos Scientific Lab. Report, LA-8095,
50. Lafarge RA, Nelsen JM and Gwinn KW. A novel CFD/ 1980.
structural analysis of a cross parachute. 99 general and 66. Benney RJ and Stein KR. Computational fluid–
miscellaneous. In: 32nd AlAA Aerospace Sciences Meeting structure interaction model for parachute inflation.
and Exhibit, Reno, NV, USA, 10–13 January 1994. J Aircraft 1996; 33: 730–736.
51. Sahu J. 3-D parachute descent analysis using coupled 67. Liu WK, Herman C, Jiun-Shyan C, et al. Arbitrary
CFD and structural codes. In: 13th Aerodynamic Lagrangian-Eulerian Petrov-Galerkin finite elements
Decelerator Systems Technology Conference, for nonlinear continua. Comput Meth Appl Mech Eng
Clearwater Beach, FL, USA, 15–18 May 1995. 1988; 68: 259–310.
52. Sahu J and Benney R. Prediction of terminal descent 68. Ibos C, Lacroix C, Chuzet L, et al. SINPA, a full 3D
characteristics of parachute clusters using CFD. In: fluid–structure software for parachute simulation. Am
14th Aerodynamic Decelerator Systems Technology Inst Aeronaut Astronaut 1997; 91: 313–323.
Conference, San Francisco, CA, USA, 3–5 June 1997. 69. Tutt B and Taylor AP. The use of LS-DYNA to simu-
53. Cao Y and Wang K. Numerical simulation of para- late the inflation of a parachute canopy. In: 18th AIAA
chute fluid–structure interaction in terminal descent. aerodynamic decelerator systems technology conference
Sci China 2012; 55: 3131–3141. and seminar, Munich, Germany, 23–26 May 2005.
54. Cao Y and Zhu X. Effects of characteristic geometric 70. Tutt B, Taylor AP, Berland JC, et al. The use
parameters on parafoil lift and drag. Aircraft Eng of LS-DYNA to assess the performance of
Aerosp Tech 2013; 85: 280–292. airborne systems North America candidate ATPS
55. Cao Y, Song Q, Wu Z, et al. Flow field and topological main parachutes. In: 18th AIAA aerodynamic deceler-
analysis of hemispherical parachute in low angles of ator systems technology conference and seminar,
attack. Modern Phys Lett B 2011; 24: 1707–1725. Munich, Germany, 23–26 May 2005.
56. Cao Y and Jiang C. Numerical simulation of the flow 71. Koschdon K and Schäfer M. A Lagrangian-Eulerian
field around parachute during terminal descent. finite-volume method for simulating free surface flows
Aircraft Eng Aerosp Tech 2007; 79: 268–272. of granular avalanches. In: Dynamic response of granu-
57. Izadi MJ and Dawoodian M. CFD analysis of drag lar and porous materials under large and catastrophic
coefficient of a parachute in a steady and turbulent con- deformations. Berlin, Germany: Springer Berlin
dition in various Reynolds numbers. In: ASME fluids Heidelberg, 2003, pp.83–108.
engineering division summer meeting, Vail, Colorado, 72. Li J, Hesse M, Ziegler J, et al. An arbitrary Lagrangian
USA, 2–6 August 2009, pp.2285–2293. New York, Eulerian method for moving-boundary problems and
NY: American Society of Mechanical Engineers. its application to jumping over water. J Comput Phys
58. Izadi MJ and Razzaz RB. 3D Numerical simulation of 2005; 208: 289–314.
a parachute with two air vented canopies in a top-to-top 73. Tezduyar TE. CFD methods for three-dimensional
formation. In: ASME fluids engineering division summer computation of complex flow problems. J Wind Eng
meeting collocated with the heat transfer, energy sustain- Ind Aerodyn 1999; 81: 97–116.
ability, and 3rd energy nanotechnology conferences, 74. Aliabadi SK and Tezduyar TE. Space-time finite elem-
Jacksonville, Florida, USA, 10–14 August 2008, ent computation of compressible flows involving
pp.435–443. New York, NY: American Society of moving boundaries and interfaces. Comput Meth Appl
Mechanical Engineers. Mech Eng 1993; 107: 209–223.
59. Fan Y and Xia J. Simulation of 3D parachute fluid– 75. Stein K, Benney R, Kalro V, et al. Parachute fluid–
structure interaction based on nonlinear finite element structure interactions: 3-D computation. Comput Meth
method and preconditioning finite volume method. Appl Mech Eng 2000; 190: 373–386.
Chin J Aeronaut 2014; 27: 1373–1383. 76. Kalro V and Tezduyar TE. A parallel 3D computa-
60. Sharov D and Nakahashi K. Low speed precondition- tional method for fluid–structure interactions in para-
ing and LU-SGS scheme for 3-D viscous flow computa- chute systems. Comput Meth Appl Mech Eng 2000; 190:
tions on unstructured grids. In: 36th Aerospace Sciences 321–332.
Cao et al. 29

77. Sathe S, Benney R, Charles R, et al. Fluid–structure 97. Valdés JG. Nonlinear analysis of orthotropic membrane
interaction modeling of complex parachute designs and shell structures including fluid–structure interaction.
with the space–time finite element techniques. Comput PhD Thesis, Universitat Politècnica De Catalunya,
Fluids 2007; 36: 127–135. Spain, 2007.
78. Takizawa K and Tezduyar TE. Computational methods 98. Zhang W, Leonard JL and Accorsi ML. Analysis of
for parachute fluid–structure interactions. Arch Comput geometrically nonlinear anisotropic membranes:
Meth Eng 2012; 19: 125–169. theory and verification. Finite Element Anal Des
79. Takizawa K, Montes D, Fritze M, et al. Methods for 2005; 41: 963–988.
FSI modeling of spacecraft parachute dynamics and 99. René de B and Crisfield MA. Nonlinear finite element
cover separation. Math Models Meth Appl Sci 2013; analysis of solids and structures. Chichester, UK: Wiley
23: 375–381. & Sons, 2012.
80. Mittal R and Iaccarino G. Immersed boundary meth- 100. Jenkins CH. Nonlinear dynamic response of mem-
ods. Ann Rev Fluid Mech 2005; 37: 239–261. branes: state of the art – update. Appl Mech Rev
81. Kim Y and Peskin CS. 2–D parachute simulation by the 1996; 49: S41.
immersed boundary method. Siam J Sci Comput 2006; 101. Accorsi M and Leonard J. Structural modeling of
28: 2294–2312. parachute dynamics. AIAA J 2000; 38: 139–146.
82. Kim Y and Peskin CS. 3-D parachute simulation by the 102. Wang J, Aquelet N, Tutt B, et al. Porous Euler-
immersed boundary method. Comput Fluids 2009; 38: Lagrange coupling: application to parachute dynam-
1080–1090. ics. In: 9th international LS-DYNA users conference,
83. Cheng Y and Zhang H. Immersed boundary method Dearborn, Michigan USA, 4–6 June 2006.
and lattice Boltzmann method coupled FSI simulation 103. Morton SA, Melville RB and Visbal MR. Accuracy
of mitral leaflet flow. Comput Fluids 2010; 39: 871–881. and coupling issues of aeroelastic Navier-Stokes solu-
84. Zhang J, Johnson PC and Popel AS. An immersed tions on deforming meshes. J Aircraft 1998; 35:
boundary lattice Boltzmann approach to simulate 798–805.
deformable liquid capsules and its application to micro- 104. Heil M. An efficient solver for the fully coupled solution
scopic blood flows. Phys Biol 2007; 4: 285–295. of large-displacement fluid–structure interaction prob-
85. Penga Y, Shua C, Chewa YT, et al. Application of lems. Comput Meth Appl Mech Eng 2004; 193: 1–23.
multi-block approach in the immersed boundary-lattice 105. Farhat C and Lesoinne M. Two efficient staggered
Boltzmann method for viscous fluid flows. J Comput algorithms for the serial and parallel solution of
Phys 2006; 218: 460–478. three-dimensional nonlinear transient aeroelastic
86. Gan Y, Xu A, Zhang G, et al. Two-dimensional lattice problems. Comput Meth Appl Mech Eng 2000; 182:
Boltzmann model for compressible flows with high 499–515.
Mach number. Phys A Stat Mech Appl 2008; 387: 106. Piperno S, Farhat C and Larrouturou B. Partitioned
1721–1732. procedures for the transient solution of coupled aero-
87. Shu C, Liu N and Chew YT. A novel immersed bound- elastic problems. Comput Meth Appl Mech Eng 2001;
ary velocity correction–lattice Boltzmann method and 190: 3147–3170.
its application to simulate flow past a circular cylinder. 107. Farhat C and Lesoinne M. Improved staggered algo-
J Comput Phys 2007; 226: 1607–1622. rithms for the serial and parallel solution of three-
88. Pan C, Luo LS and Miller C T. An evaluation of lattice dimensional nonlinear transient aeroelastic problems.
Boltzmann schemes for porous medium flow simula- Comput Meth Appl Mech Eng 2000; 182: 499–515.
tion. Comput Fluids 2006; 35: 898–909. 108. Fernández MA. Coupling schemes for incompressible
89. Gong RH, Ozgen B and Soleimani M. Modeling of fluid–structure interaction: implicit, semi-implicit and
yarn cross-section in plain woven fabric. Text Res J explicit. SeMa J 2011; 55: 59–108.
2009; 79: 1014–1020. 109. Causin P, Gerbeau JF and Nobile F. Added-mass
90. Chichani S and Guha A. A method of modeling fabric effect in the design of partitioned algorithms for
shear using finite element analysis. J Inst Eng 2014; 96: fluid–structure problems. Comput Meth Appl Mech
1–7. Eng 2005; 194: 4506–4527.
91. Roberts BW. Aerodynamic inflation of shell type, para- 110. Brenner S C. The condition number of the Schur com-
chute structures. J Aircraft 1974; 11: 390–397. plement in domain decomposition. Numer Math 1998;
92. Reynolds D. Stress analysis of ribbon parachutes. In: 83: 187–203.
5th Aerodynamic Deceleration Systems Conference, 111. Tu C and Peskin C S. Stability and instability in the
Albuquerque, NM, USA, 17–19 November 1975. computation of flows with moving immersed bound-
93. Sundberg WD. New solution method for steady-state aries: a comparison of three methodsJ. Siam J Sci Stat
canopy structural loads. J Aircraft 2012; 25: 1045–1051. Comput 1992; 13: 1361–1376.
94. Yihua C, Zhuo W, Qianfu S, et al. Numerical simula- 112. Boffi D, Gastaldi L and Heltai L. On the CFL condi-
tion of fluid–structure interaction in the opening pro- tion for the finite element immersed boundary method.
cess of conical parachute. Aeronaut J 2009; 113: Comput Struct 2007; 85: 775–783.
191–200. 113. Hoppe RHW and Linsenmann C. An adaptive
95. Kim JD, Li Y and Li X. Simulation of parachute FSI Newton continuation strategy for the fully implicit
using the front tracking method. J Fluids Struct 2013; finite element immersed boundary method. J Comput
37: 100–119. Phys 2012; 231: 4676–4693.
96. Jenkins CH and Leonard JW. Nonlinear dynamic 114. Brezzi F and Fortin M. Mixed and hybrid finite elem-
response of membranes: state of the art. Appl Mech ent methods. New York: Springer, Berlin–Heidelberg,
Rev 1991; 44: 319–328. 1991.
30 Proc IMechE Part G: J Aerospace Engineering 0(0)

115. Wong KK, Thavornpattanapong P, Cheung SC, et al. 132. Boer AD, Schoot MSVD and Bijl H. Mesh deform-
Numerical stability of partitioned approach in ation based on radial basis function interpolation.
fluid–structure interaction for a deformable thin- Comput Struct 2007; 85: 784–795.
walled vessel. Comput Math Meth Med 2013; 2013: 133. Farhat C, Degand C, Koobus B, et al. Torsional
437–452. springs for two-dimensional dynamic unstructured
116. Deparis S, Fernández MAF, Luca FMA, et al. fluid meshes. Comput Meth Appl Mech Eng 1998;
Acceleration of a fixed point algorithm for fluid–struc- 163: 231–245.
ture interaction using transpiration conditions. 134. Batina JT. Unsteady Euler algorithm with unstruc-
ESAIM Math Model Num Anal 2003; 37: 1307–1311. tured dynamic mesh for complex-aircraft aerodynamic
117. Vierendeels J, Lanoye L, Degroote J, et al. analysis. AIAA J 1991; 29: 327–333.
implicit coupling of partitioned fluid–structure inter- 135. Zeng D and Ethier CR. A semi-torsional spring ana-
action solvers using reduced-order models: implicit logy model for updating unstructured meshes in 3D
coupling of partitioned fluid–structure interaction sol- moving domains. Finite Elem Anal Des 2005; 41:
vers using a reduced order model. Fluid-Structure 1118–1139.
Interaction 2005; 53: 1–18. 136. Stein K, Tezduyar T and Benney R. Mesh moving
118. Küttler U and Wall WA. Vector extrapolation for techniques for fluid–structure interactions with large
strong coupling fluid–structure interaction solvers. displacements. J Appl Mech 2003; 70: 58–63.
J Appl Mech 2009; 76: 353–360. 137. Benney R, Tezduyar T and Stein K. Computational
119. Järvinen E, Råback P, Lyly M, et al. A method methods for modeling parachute systems. Comput Sci
for partitioned fluid–structure interaction computa- Eng 2003; 5: 39–46.
tion of flow in arteries. Med Eng Phys 2008; 30: 138. Nie S, Cao Y and Wu Z. Numerical simulation of
917–923. parafoil inflation via a Robin–Neumann transmis-
120. Badia S, Nobile F and Vergara C. Fluid–structure par- sion-based approach. Proc IMechE, Part G: J
titioned procedures based on Robin transmission con- Aerospace Engineering 2017; 1–14.
ditions. Int J Account Educ Res 2008; 27: 333–335. 139. Higuchi H and Takahashi F. Flow past two-dimen-
121. Bukač M, Čanić S, Glowinski R, et al. Fluid–structure sional parachute models. J Aircraft 1989; 26: 641–649.
interaction in blood flow capturing non-zero longitu- 140. Tezduyar TE, Sathe S, Schwaab M, et al. Fluid-struc-
dinal structure displacement. J Comput Phys 2012; ture interaction modeling of ringsail parachutes.
235: 515–541. Comput Mech 2008; 43: 133–142.
122. Canic S, Muha B and Bukac M. Stability of the kine- 141. Donald AN and Adrian B. Convection in porous media.
matically coupled b-scheme for fluid–structure inter- Berlin, Germany: Springer Berlin, 2013, pp.371–377.
action problems in hemodynamics. Mathematics 142. Aquelet N and Tutt B. Euler-Lagrange coupling for
2012; 12: 54–80. porous parachute canopy analysis. Int J Multiphys
123. Tezduyar TE and Sathe S. Modelling of fluid–struc- 2007; 1: 53–68.
ture interactions with the space–time finite elements: 143. Roddeman DG, Drukker J, Oomens CWJ, et al. The
solution techniques. Int J Num Meth Fluids 2007; 54: wrinkling of thin membranes: part I – theory. J Appl
855–900. Mech 1987; 54: 884–887.
124. Thévenaz P, Blu T and Unser M. Interpolation revis- 144. Liu X, Jenkins CH and Schur WW. Large deflection
ited. IEEE Transact Med Imag 2000; 19: 739–758. analysis of pneumatic envelopes using a penalty par-
125. Cebral J and Lohner R. Conservative load projection ameter modified material model. Finite Elem Anal Des
and tracking for fluid–structure problems. AIAA J 2001; 37: 233–251.
1997; 35: 687–692. 145. Gil AJ and Bonet J. Finite element analysis of partly
126. Boer AD, Zuijlen AHV and Bijl H. Review of coupling wrinkled reinforced prestressed membranes. Comput
methods for non-matching meshes. Comput Meth Appl Mech 2007; 40: 595–615.
Mech Eng 2007; 196: 1515–1525. 146. Lu K, Accorsi M and Leonard J. Finite element ana-
127. Smith MJ, Cesnik CES and Hodges DH. Evaluation lysis of membrane wrinkling. Int J Num Meth Eng
of some data transfer algorithms for noncontiguous 2001; 50: 1017–1038.
meshes. J Aerosp Eng 2014; 13: 52–58. 147. Rossi R, Lazzari M, Vitaliani R, et al. Simulation of
128. Beckert A and Wendland H. Multivariate interpol- light-weight membrane structures by wrinkling model.
ation for fluid–structure-interaction problems using Int J Num Meth Eng 2005; 62: 2127–2153.
radial basis functions. Aerosp Sci Technol 2001; 5: 148. Laursen TA. Computational contact and impact mech-
125–134. anics: fundamentals of modeling interfacial phenomena
129. Zhang W, Gao C and Ye Z. Research progress on in nonlinear finite element analysis. New York, NY,
mesh deformation method in computational aeroelas- USA: Springer, 2002.
ticity. Acta Aeronaut Astronaut Sin 2014; 35: 303–319. 149. Hughes TJR, Taylor RL, Sackman JL, et al. A finite
130. Eriksson LE. Generation of boundary-conforming element method for a class of contact impact pro-
grids around wing-body configurations using trans- blems. Comput Meth Appl Mech Eng 1976; 8:
finite interpolation. AIAA J 1982; 20: 1313–1320. 249–276.
131. Liu X, Qin N and Xia H. Fast dynamic grid deform- 150. Parisch H. A consistent tangent stiffness matrix for
ation based on Delaunay graph mapping. J Comput three-dimensional non-linear contact analysis. Int J
Phys 2006; 211: 405–423. Num Meth Eng 1989; 28: 1803–1812.
Cao et al. 31

151. Laursen TA and Simo JC. A continuum-based finite 155. Konyukhov A and Schweizerhof K. Computational
element formulation for the implicit solution of multi- contact mechanics: geometrically exact theory for arbi-
body, large deformation-frictional contact problems. trary shaped bodies. New York, NY, USA: Springer
Int J Num Meth Eng 1993; 36: 3451–3485. Science & Business Media, 2012.
152. Meguid SA and Czekanski A. Advances in computa- 156. Sathe S and Tezduyar TE. Modeling of fluid–structure
tional contact mechanics. Int J Num Meth Eng 2008; 4: interactions with the space–time finite elements: con-
419–443. tact problems. Comput Mech 2008; 43: 51–60.
153. Belytschko T and Neal MO. Contact-impact by the 157. Flores P and Ambrósio J. On the contact detection for
pinball algorithm with penalty and Lagrangian meth- contact-impact analysis in multibody systems.
ods. Int J Num Meth Eng 1991; 31: 547–572. Multibody Syst Dyn 2010; 24: 103–122.
154. Maruyama K. A procedure to determine intersections 158. Ganter MA and Isarankura BP. Dynamic collision
between polyhedral objects. Int J Comput Inform Sci detection using space partitioning. J Mech Des 1993;
1972; 1: 255–266. 115: 150–155.