**Fluid-Structure Interaction Techniques for Parachute**

Vinod Kumar1 and Victor Udoewa2 *1The University of Texas at El Paso 2American Association for the Advancement of Science Policy Fellow USA* 

#### **1. Introduction**

238 Fluid Dynamics, Computational Modeling and Applications

[19] Zhang Y. H.; Yang C.; Mao Z. S. Large eddy simulation of the gas-liquid flow in

[20] Jain M.; Paranandi M.; Roush D.; Göklen K.; Kelly W. J. Using CFD to understand

[21] Panneerselvam R.; Savithri S.; Surender G. D. CFD simulation of hydrodynamics of gas-liquid-solid fluidised bed reactor. *Chem. Eng. Sci.* 2009, 64, 1119

how flow patterns affect retention of cell-sized particles in a tubular bowl

sstirred tank. AIChE J. 2008, 54, 1963

centrifuge. *Ind. Eng. Chem. Res.* 2005, 44, 7876

Fluid-Structure Interaction (FSI) is an inescapable feature of the complicated flow physics where strong coupling between fluid dynamics and structural dynamics occurs. One such example is to understand the behavior of parachutes during canopy inflation and decent. Such simulations can substantially reduce the design costs of parachutes by reducing the number of rather expensive experiments/airdrop tests required. Additionally, FSI simulations can augment experimental approaches by providing detailed fluid flow and structural deformation characteristics of the parachute systems under various scenarios. There are however many computational challenges that one faces in performing the high fidelity FSI simulations. Some of these are turbulence modeling, fluid-structure interaction coupling, the convergence of a nonlinear iteration loop, capturing the physical discontinuity in the flow field, efficiently and accurately solving a large set of linear equations on parallel computers, and improving the parallel performance. Here, we present a high fidelity FSI technique that addresses the coupling issues on High Performance Computing (HPC) environments.

#### **2. Mathematical formulations**

FSI modeling of parachutes requires simultaneously solving the Navier-Stokes equations (a set of highly nonlinear Partial Differential Equations (PDEs)) for fluid dynamics, structural mechanics, and mesh motions. Parachutes are made of membrane type structure (usually nylon clothes) that goes through large deformation as a result of aerodynamic forces. This often causes the FSI simulations to break down because of the convergence issues of nonlinear and linear iterative solvers and large mesh stretching/distortions. It is, therefore, difficult to operate and requires in-depth knowledge of numerical techniques, fundamentals of fluid and structure dynamics, coupling behavior, programming languages and environments on High-Performance Computing (HPC) systems. Preparations required before FSI simulations for the parachute-systems are equally challenging and poses tremendous difficulties for novice users.

#### **2.1 Fluid dynamics (FD)**

The physics of fluid dynamics is mathematically represented by the Navier-Stokes equations for compressible and incompressible flows. These equations represent conservation of mass

$$\nabla.\mathfrak{u} = \mathbf{0} \text{ on } \Omega\_t \,\forall \, t \in (0, T)T \tag{1}$$

$$\rho \left( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u}. \nabla \mathbf{u} - \mathbf{f} \right) - \nabla \cdot \sigma = \mathbf{0} \text{ on } \Omega\_t \,\forall \, t \in \{0, T\} \tag{2}$$

$$\sigma(p,\mathfrak{u}) = -pI + \mathfrak{r} \text{ on } \Omega\_t \,\forall \, t \in (0, T) \tag{3}$$

$$
\mathfrak{r} = 2\mu \dot{\mathfrak{e}} \tag{4}
$$

$$\mathfrak{E}(\mathbf{u}) = \frac{1}{2} \left( (\nabla \mathbf{u}) + (\nabla \mathbf{u})^{\mathsf{T}} \right) \tag{5}$$

$$\frac{\partial \rho}{\partial t} + \nabla. (\rho \mathbf{u}) = \mathbf{0} \tag{6}$$

$$\frac{\partial(\rho \mathbf{u})}{\partial t} + \nabla.\left(\rho \mathbf{u} \mathbf{u}\right) + \nabla p - \nabla.\mathbf{T} = 0\tag{7}$$

$$\frac{\partial(\rho e)}{\partial t} + \nabla.(\rho e \mathbf{u}) + \nabla(p \mathbf{u}) - \nabla.\mathbf{T} \mathbf{u} + \nabla.\mathbf{q} = 0\tag{8}$$

$$\frac{\partial \mathbf{U}}{\partial t} + \frac{\partial \mathbf{F}\_l}{\partial x\_l} - \frac{\partial \mathbf{E}\_l}{\partial x\_l} = \mathbf{0} \tag{9}$$

$$\mathbf{U} = \begin{pmatrix} \rho, \rho u\_1, \dots, \rho u\_{n\_{\text{sd}}}, \rho e \end{pmatrix}^T \tag{10}$$

$$\mathbf{F}\_{l} = \left( u\_{l} \rho, \ u\_{l} \rho u\_{1} + \delta\_{l1} \ p, \dots, u\_{l} \rho u\_{n\_{\text{sd}}} + \delta\_{l n\_{\text{sd}}} \ p, u\_{l} (\rho e + p) \right)^{\text{tr}} \tag{11}$$

$$\mathbf{E}\_{l} = \begin{pmatrix} 0, \tau\_{l1}\rho, \ \ldots, \tau\_{l n\_{sd}}, -q\_{l} + \tau\_{lk}u\_{k} \end{pmatrix}^{\mathrm{T}} \tag{12}$$

$$\frac{\partial \mathbf{u}}{\partial t} + \mathbf{A}\_{\mathbf{l}} \frac{\partial \mathbf{u}}{\partial x\_{l}} - \frac{\partial}{\partial x\_{l}} \left( \mathbf{K}\_{\mathbf{l}\mathbf{j}} \frac{\partial \mathbf{u}}{\partial x\_{j}} \right) = \mathbf{0} \tag{13}$$

$$\mathbf{A}\_{\mathbf{l}} = \frac{\partial \mathbf{F}\_{\mathbf{l}}}{\partial \mathbf{u}} \tag{14}$$

$$\mathbf{K}\_{\rm li} \frac{\partial \mathbf{U}}{\partial x\_{\rm J}} = \mathbf{E}\_{\rm I} \tag{15}$$

$$
\rho^s \left( \frac{d^2 \mathbf{y}}{dt^2} + \eta \frac{d \mathbf{y}}{dt} - \mathbf{f}^s \right) - \nabla. \sigma^s = 0 \tag{16}
$$

Fluid-Structure Interaction Techniques for Parachute 243

Galerkin approximation methods are used for the finite element formulations. Galerkin methods are among the most commonly used weighted residual methods. In structural analysis, where often the minimization of energy is the underlying idea, the application of Galerkin methods leads to symmetric matrices and provides *optimal* results. By *optimal* we mean that the solution possesses the best approximation property. The difference between the approximate and the exact solutions is minimized with respect to a certain norm as

The situation, however, is very different in the presence of advective terms. The matrix associated with the advective term is non-symmetric and the best approximation property is lost (Brooks & Hughes, 1982). As a result Galerkin methods applied to these problems are far from optimal and show spurious node-to-node oscillations in the solutions, worsening with growing advection-domination. This not only leads to qualitatively incorrect results but also violates basic physical principles like the second law of thermodynamics (Hirsch, 1988). The pollution of the solution with oscillations is dependent on the domination of the advection terms over other terms of the differential equation. Domination of advection is determined by dimensionless numbers such as Reynolds (ratio of inertial to viscous terms) or Peclet (ratio of advection to diffusion terms). The larger these numbers are, the more dominant is the advection term and the stronger is the numerical node-to-node oscillations in the results. Brooks and Hughes (Brooks & Hughes, 1982) introduced the Streamline-Upwind/Petrov-Galerkin (SUPG) method and this method can be considered as the first successful stabilization technique to prevent oscillations in advection-dominated problems. The SUPG method introduces artificial diffusion in the streamline direction. The introduction of artificial diffusion is done in a consistent way. This can be interpreted as a modification of the test function in the advection direction. Therefore, the weak (variational) form still satisfies the exact

Another source of potential instabilities in standard Galerkin methods arises when velocity and pressure interpolation functions are not chosen from compatible spaces. Babuska and Brezzi (Babuska, 1973; Brezzi, 1974) showed that compatible spaces must satisfy the *inf –sup*  conditions. If they are not chosen from compatible spaces, then oscillatory behavior is observed, primarily in the pressure field. Unfortunately, much desirable function spaces are precluded due to Babuska-Brezzi conditions. Most computationally attractive combinations are the ones which employ equal order interpolations. Hughes et alproposed a consistent way to circumvent the *inf -sup* condition for the Stokes problem (Hughes, Franca, & Mallet, 1986). As a generalization, Tezduyar, et al., proposed a Pressure-Stabilized/Petrov-Galerkin (PSPG) formulation for finite Reynolds number flows (Tezduyar, Mittal, Ray, & Shih, 1992). The PSPG formulations reduces to the one proposed by Hughes et al. formulation in the

**3.1 Deforming spatial domain/stabilized space-time (DSD/SST) formulations for** 

In order to construct the \_nite element function spaces for the space-time method, the time interval ��� �) is partitioned into sub-intervals �� = ���� ����), where �� and ���� belong to an ordered sequence of time levels � = �� � �� ����� = �. Let Ω� = �� and Γ� = ��. The space-time slab �� is defined as the domain enclosed by the surfaces Ω�, ���, and P�, where P� is the surface described by the boundary Γ� as t traverses I�. As is the

shown by (Brooks & Hughes, 1982).

solution of the problem.

**incompressible flows** 

limiting case when the Reynolds number tends to zero.

where �� is the material density, **y** is the displacement vector, **f**s is the external force body forces, �� is the Cauchy stress tensor, and � is the mass proportional damping coefficient. The mass-proportional damping provides additional stability, but can significantly affect the dynamics of the structure. Here, we assume large displacements and rotations, but small strains for nonlinear analysis. The second Piola-Kirchhoff stress (force per unit area in the original configuration) tensor, **S**, and the Green-Lagrange strain (in the original configuration) tensor, **E**, are used to write the constitutive equations using the total Lagrangian formulation. Thus, stresses are expressed in terms of the 2nd Piola-Kirchoff stress tensor. The 2nd Piola-Kirchhoff stress tensor is related to Cauchy stress (force per unit area in deformed configuration) tensor, ��, by the kinematic transformation, � = �� �� ����� where ρ� is the density in the original configuration and **F** (=**gG**; **g**=covariant tensor in deformed configuration, **G**=contravariant tensor in original configuration). Firstly, we will assume linear stress-strain relations (Hookean materials) and plane stress conditions. The constitutive equations are given by

$$\mathcal{S}^{ij} = \{\bar{\lambda}G^{ij}G^{jk} + \mu\{G^{il}G^{jk} + G^{ik}G^{jl}\}\} E\_{kl} \tag{17}$$

where �̅= ��� ���� (� � � are the Lame constants). The Lame constants are related to the Young's modulus *Y* and the Poisson's � ratio by: �̅= �� (���)(����) and � = � �(���) . Dirichlet- and Neumann-type BCs are *y=g*s and **n**. ��=**t**s where *g*s is the specified displacements and **t**s is traction forces (shear stress from the fluid dynamics). The initial conditions are **<sup>y</sup>** = � ��**<sup>y</sup>** �� = �.

#### **2.3 Mesh deformation (MD)**

Fluid mesh is considered as elastic materials which deforms along with SD deformation. The governing equations mesh deformation is given by

$$\nabla.\sigma^{\mathsf{M}} = 0\tag{18}$$

#### **3. Finite element discretization**

The finite element method is a numerical tool for obtaining solutions to boundary value engineering problems governed by partial differential equations as described in the previous section. It is especially attractive for problems that involve complex geometries where other numerical methods, such as spectral or finite difference methods, are difficult to apply. The principle of the method is to replace an entire continuous domain by a number of sub-domains in which the unknown function is represented by simple interpolation functions with unknown coefficients. Thus, the original boundary-value problem with an infinite number of degrees of freedom is converted into a problem with a finite number of degrees of freedom or, in other words, the solution of the whole system is approximated by a finite number of unknown coefficients. Here, we construct a discretization of a weighted residual formulation in order to arrive at a linear matrix equation. The discretization can be applied in space and time. This formulation is called the space-time finite element method. Alternatively, one can carry out the discretization in space only. Such a method is called a semi-discrete formulation. Details of these methods can be found in a number of references (Tezduyar, Behr, & Liou, 1992; Shakib, 1988; Behr & Tezduyar, 1994).

where �� is the material density, **y** is the displacement vector, **f**s is the external force body forces, �� is the Cauchy stress tensor, and � is the mass proportional damping coefficient. The mass-proportional damping provides additional stability, but can significantly affect the dynamics of the structure. Here, we assume large displacements and rotations, but small strains for nonlinear analysis. The second Piola-Kirchhoff stress (force per unit area in the original configuration) tensor, **S**, and the Green-Lagrange strain (in the original configuration) tensor, **E**, are used to write the constitutive equations using the total Lagrangian formulation. Thus, stresses are expressed in terms of the 2nd Piola-Kirchoff stress tensor. The 2nd Piola-Kirchhoff stress tensor is related to Cauchy stress (force per unit

area in deformed configuration) tensor, ��, by the kinematic transformation, � = ��

constitutive equations are given by

**2.3 Mesh deformation (MD)** 

**3. Finite element discretization** 

where �̅= ���

��� = ��̅

modulus *Y* and the Poisson's � ratio by: �̅= ��

governing equations mesh deformation is given by

where ρ� is the density in the original configuration and **F** (=**gG**; **g**=covariant tensor in deformed configuration, **G**=contravariant tensor in original configuration). Firstly, we will assume linear stress-strain relations (Hookean materials) and plane stress conditions. The

Neumann-type BCs are *y=g*s and **n**. ��=**t**s where *g*s is the specified displacements and **t**s is traction forces (shear stress from the fluid dynamics). The initial conditions are **<sup>y</sup>** = � ��**<sup>y</sup>**

Fluid mesh is considered as elastic materials which deforms along with SD deformation. The

The finite element method is a numerical tool for obtaining solutions to boundary value engineering problems governed by partial differential equations as described in the previous section. It is especially attractive for problems that involve complex geometries where other numerical methods, such as spectral or finite difference methods, are difficult to apply. The principle of the method is to replace an entire continuous domain by a number of sub-domains in which the unknown function is represented by simple interpolation functions with unknown coefficients. Thus, the original boundary-value problem with an infinite number of degrees of freedom is converted into a problem with a finite number of degrees of freedom or, in other words, the solution of the whole system is approximated by a finite number of unknown coefficients. Here, we construct a discretization of a weighted residual formulation in order to arrive at a linear matrix equation. The discretization can be applied in space and time. This formulation is called the space-time finite element method. Alternatively, one can carry out the discretization in space only. Such a method is called a semi-discrete formulation. Details of these methods can be found in a number of references

(Tezduyar, Behr, & Liou, 1992; Shakib, 1988; Behr & Tezduyar, 1994).

���� (� � � are the Lame constants). The Lame constants are related to the Young's

������ ��������� � ����������� (17)

�� �� = 0 (18)

� �(���)

(���)(����) and � =

�� �����

. Dirichlet- and

�� = �. Galerkin approximation methods are used for the finite element formulations. Galerkin methods are among the most commonly used weighted residual methods. In structural analysis, where often the minimization of energy is the underlying idea, the application of Galerkin methods leads to symmetric matrices and provides *optimal* results. By *optimal* we mean that the solution possesses the best approximation property. The difference between the approximate and the exact solutions is minimized with respect to a certain norm as shown by (Brooks & Hughes, 1982).

The situation, however, is very different in the presence of advective terms. The matrix associated with the advective term is non-symmetric and the best approximation property is lost (Brooks & Hughes, 1982). As a result Galerkin methods applied to these problems are far from optimal and show spurious node-to-node oscillations in the solutions, worsening with growing advection-domination. This not only leads to qualitatively incorrect results but also violates basic physical principles like the second law of thermodynamics (Hirsch, 1988). The pollution of the solution with oscillations is dependent on the domination of the advection terms over other terms of the differential equation. Domination of advection is determined by dimensionless numbers such as Reynolds (ratio of inertial to viscous terms) or Peclet (ratio of advection to diffusion terms). The larger these numbers are, the more dominant is the advection term and the stronger is the numerical node-to-node oscillations in the results. Brooks and Hughes (Brooks & Hughes, 1982) introduced the Streamline-Upwind/Petrov-Galerkin (SUPG) method and this method can be considered as the first successful stabilization technique to prevent oscillations in advection-dominated problems. The SUPG method introduces artificial diffusion in the streamline direction. The introduction of artificial diffusion is done in a consistent way. This can be interpreted as a modification of the test function in the advection direction. Therefore, the weak (variational) form still satisfies the exact solution of the problem.

Another source of potential instabilities in standard Galerkin methods arises when velocity and pressure interpolation functions are not chosen from compatible spaces. Babuska and Brezzi (Babuska, 1973; Brezzi, 1974) showed that compatible spaces must satisfy the *inf –sup*  conditions. If they are not chosen from compatible spaces, then oscillatory behavior is observed, primarily in the pressure field. Unfortunately, much desirable function spaces are precluded due to Babuska-Brezzi conditions. Most computationally attractive combinations are the ones which employ equal order interpolations. Hughes et alproposed a consistent way to circumvent the *inf -sup* condition for the Stokes problem (Hughes, Franca, & Mallet, 1986). As a generalization, Tezduyar, et al., proposed a Pressure-Stabilized/Petrov-Galerkin (PSPG) formulation for finite Reynolds number flows (Tezduyar, Mittal, Ray, & Shih, 1992). The PSPG formulations reduces to the one proposed by Hughes et al. formulation in the limiting case when the Reynolds number tends to zero.

#### **3.1 Deforming spatial domain/stabilized space-time (DSD/SST) formulations for incompressible flows**

In order to construct the \_nite element function spaces for the space-time method, the time interval ��� �) is partitioned into sub-intervals �� = ���� ����), where �� and ���� belong to an ordered sequence of time levels � = �� � �� ����� = �. Let Ω� = �� and Γ� = ��. The space-time slab �� is defined as the domain enclosed by the surfaces Ω�, ���, and P�, where P� is the surface described by the boundary Γ� as t traverses I�. As is the

$$(\mathcal{S}\_u^\hbar)\_n = \left\{ u^\hbar \middle| u^\hbar \in \left[ H^{1\hbar} (Q\_n)^{n\_{\rm sd}} \right]\_\prime \ u^\hbar \equiv g\_u^\hbar \text{ on } (\mathcal{P}\_n)\_{\mathcal{g}^\prime} \ d = 1, \ldots, n\_{\rm sd} \right\} \tag{19}$$

$$(\mathcal{V}\_{\boldsymbol{u}}^{\hbar})\_{\boldsymbol{n}} = \left\{ \boldsymbol{w}^{\hbar} | \boldsymbol{w} \in \left[ H^{1\hbar} (\mathbb{Q}\_{\boldsymbol{n}})^{n\_{\mathrm{sd}}} \right], \; \boldsymbol{w}^{\hbar} \equiv \boldsymbol{0} \text{ on } \left( \mathbb{P}\_{\boldsymbol{n}} \right)\_{\mathcal{G}}, \; \boldsymbol{d} = 1, \ldots, n\_{\mathrm{sd}} \right\} \tag{20}$$

$$\left(\left(\mathfrak{S}\_p^h\right)\_n = \left(\mathcal{V}\_p^h\right)\_n = \left\{q^h \middle| q^h \in \left[H^{1h}(\mathbb{Q}\_n)^{n\_{\le d}}\right]\_\* \text{ w}^h \cong \equiv 0 \text{ on } (\mathbb{P}\_n)\_g\right\}\tag{21}$$

$$\begin{split} & \int\_{\Omega} \mathbf{w}^{h} \cdot \boldsymbol{\rho} \left( \frac{\widehat{\boldsymbol{\alpha}}^{h}}{\widehat{\boldsymbol{\alpha}}} + \mathbf{u}^{h} \cdot \nabla \mathbf{u}^{h} + \boldsymbol{f}^{h} \right) d\boldsymbol{Q} + \int\_{\Omega} \dot{\boldsymbol{\varepsilon}} \left( \mathbf{w}^{h} \right) \cdot \boldsymbol{\sigma} \left( \boldsymbol{p}^{h}, \boldsymbol{u}^{h} \right) d\boldsymbol{Q} + \int\_{\Omega\_{s}} \boldsymbol{q}^{h} \nabla \cdot \mathbf{u}^{h} d\boldsymbol{Q} + \\ & \quad + \sum\_{e=1}^{\left( \boldsymbol{u}\_{e} \right)\_{h}} \int \frac{\pi}{\rho} \left[ \boldsymbol{\rho} \left( \frac{\widehat{\boldsymbol{\alpha}} \mathbf{w}^{h}}{\widehat{\boldsymbol{\alpha}}} + \mathbf{u}^{h} \cdot \nabla \mathbf{w}^{h} + \boldsymbol{f}^{h} \right) - \nabla \boldsymbol{\sigma} \left( \boldsymbol{q}^{h}, \boldsymbol{w}^{h} \right) \right] \Bigg[ \int \rho \left( \frac{\widehat{\boldsymbol{\alpha}} \mathbf{u}^{h}}{\widehat{\boldsymbol{\alpha}}} + \mathbf{u}^{h} \cdot \nabla \mathbf{u}^{h} + \boldsymbol{f}^{h} \right) - \nabla \boldsymbol{\sigma} \left( \boldsymbol{p}^{h}, \boldsymbol{u}^{h} \right) \Bigg] dQ + \left. \tag{22} \\ & + \sum\_{e=1}^{\left( \boldsymbol{u}\_{e} \right)\_{e}} \int \boldsymbol{\delta\nabla} \mathbf{w}^{h} \boldsymbol{\rho} \boldsymbol{\sigma}^{h} \boldsymbol{d\$$

$$\tau = \left[ \left( \frac{2}{\Delta t} \right)^2 + \left( \frac{2 \left| \left| u^{\hbar} \right| \right|}{h\_e} \right)^2 + \left( \frac{4 \cdot \nu}{h\_e^3} \right)^2 \right]^{-\frac{1}{2}} \tag{23}$$

$$\mathcal{S} = \frac{\hbar\_{\varepsilon}}{2} \left| |u^{\hbar}| \right| z \tag{24}$$

$$\mathcal{Z} = \begin{cases} \left(\frac{Re\_u}{3}\right) \text{ if } Re\_u \le 3\\ 1 \text{ if } Re\_u > 3 \end{cases} \tag{25}$$

$$\mathbf{N}\_1(\mathbf{d}\_1, \mathbf{d}\_2, \mathbf{d}\_3) = 0 \tag{26}$$

$$\mathbf{N}\_2(\mathbf{d}\_1, \mathbf{d}\_2, \mathbf{d}\_3) = 0 \tag{27}$$

$$\mathbf{N}\_3(\mathbf{d}\_1, \mathbf{d}\_2, \mathbf{d}\_3) = 0 \tag{28}$$

$$\left[\frac{\partial \mathbf{N}\_1}{\partial \mathbf{d}\_1}\right]^l \Delta \mathbf{d}\_1^{l+1} + \left[\frac{\partial \mathbf{N}\_1}{\partial \mathbf{d}\_2}\right]^l \Delta \mathbf{d}\_2^{l+1} + \left[\frac{\partial \mathbf{N}\_1}{\partial \mathbf{d}\_3}\right]^l \Delta \mathbf{d}\_3^{l+1} = -[\mathbf{N}\_1(\mathbf{d}\_1, \mathbf{d}\_2, \mathbf{d}\_3)]^l \tag{29}$$

$$
\left[\frac{\partial \mathbf{N}\_2}{\partial \mathbf{d}\_1}\right]^l \Delta \mathbf{d}\_1^{l+1} + \left[\frac{\partial \mathbf{N}\_2}{\partial \mathbf{d}\_2}\right]^l \Delta \mathbf{d}\_2^{l+1} + \left[\frac{\partial \mathbf{N}\_2}{\partial \mathbf{d}\_3}\right]^l \Delta \mathbf{d}\_3^{l+1} = -[\mathbf{N}\_2(\mathbf{d}\_1, \mathbf{d}\_2, \mathbf{d}\_3)]^l \tag{30}
$$

$$
\left[\frac{\partial \mathbf{N}\_3}{\partial \mathbf{d}\_1}\right]^l \Delta \mathbf{d}\_1^{l+1} + \left[\frac{\partial \mathbf{N}\_3}{\partial \mathbf{d}\_2}\right]^l \Delta \mathbf{d}\_2^{l+1} + \left[\frac{\partial \mathbf{N}\_3}{\partial \mathbf{d}\_3}\right]^l \Delta \mathbf{d}\_3^{l+1} = -[\mathbf{N}\_3(\mathbf{d}\_1, \mathbf{d}\_2, \mathbf{d}\_3)]^l \tag{31}
$$

$$
\begin{bmatrix}
\mathbf{A}\_{11} & \mathbf{A}\_{12} & \mathbf{A}\_{13} \\
\mathbf{A}\_{21} & \mathbf{A}\_{22} & \mathbf{A}\_{23} \\
\mathbf{A}\_{31} & \mathbf{A}\_{32} & \mathbf{A}\_{33}
\end{bmatrix}^{l} \begin{bmatrix}
\Delta \mathbf{d}\_{1} \\
\Delta \mathbf{d}\_{2} \\
\Delta \mathbf{d}\_{3}
\end{bmatrix}^{l+1} = \begin{bmatrix}
\mathbf{b}\_{1} \\
\mathbf{b}\_{2} \\
\mathbf{b}\_{3}
\end{bmatrix}^{l} \tag{32}
$$

$$\mathbf{A}\_{\rm jk} = \frac{\partial \mathbf{N}\_{\parallel}}{\partial \mathbf{d}\_{\mathbf{k}}} \tag{33}$$

$$\mathbf{b}\_{\rangle} = -\left[\mathbf{N}\_{\rangle}(\mathbf{d}\_1, \mathbf{d}\_2, \mathbf{d}\_3)\right] \tag{34}$$

$$
\begin{bmatrix}
\mathbf{A}\_{11} & \mathbf{0} & \mathbf{0} \\
\mathbf{A}\_{21} & \mathbf{A}\_{22} & \mathbf{0} \\
\mathbf{A}\_{31} & \mathbf{A}\_{32} & \mathbf{A}\_{33}
\end{bmatrix}^{l}
\begin{bmatrix}
\Delta \mathbf{d}\_{1} \\
\Delta \mathbf{d}\_{2} \\
\Delta \mathbf{d}\_{3}
\end{bmatrix}^{l+1} = \begin{bmatrix}
\mathbf{b}\_{1} \\
\mathbf{b}\_{2} \\
\mathbf{b}\_{3}
\end{bmatrix}^{l} \tag{35}
$$

$$
\rho = \rho R \theta = (\chi - 1)\rho i = (\chi - 1) \left[ \rho e - \frac{\left\| \rho u \right\|^2}{2\rho} \right] \tag{36}
$$

$$[U\_1, U\_2, U\_3] = [\rho u\_1, \rho u\_2, \rho u\_3]\_{\text{interface}} \tag{37}$$

Fluid-Structure Interaction Techniques for Parachute 251

The parachute canopy is made of very flexible fabric materials. It responds quickly to any small changes in fluid behavior. Consequently, there is a strong interaction between the fluid and structure. Bigger the canopies, stronger will be the interactions. Stronger interactions make the problem very nonlinear in nature and hence increase the difficulty of

A single G-12 parachute is used to drop cargo weighing up to 2,200 lb. There is an interest in the Army for parachute systems that can drop cargos exceeding 2,200 lb. In order to achieve

The opening stage of a parachute to a fully inflated shape is a critical issue. Usually more than one stage is used to get the fully inflated shape of parachutes in a cluster configuration. One example of multi-stage opening is the use of drogue parachutes. First, a drogue parachute opens and then it pulls the other two G-12 parachutes. After all the parachutes are fully opened, the challenge is to safely drop the cargo into the hands of the right people. The strong parachute-to-parachute aerodynamic interaction can potentially destroy the efficiency of the system if it is not carefully designed. Various researchers have analyzed the cluster of parachute systems using both drop tests and computer simulations. Butler (Butler, 2001) presented the qualitative results from the airdrop tests using a drogue parachute to deploy the G-12 canopies. From the snapshots of a cluster of three G-12 canopy airdrop tests, one clearly observes strong interactions between fluid and canopies. Sahu, et al., (Sahu & Benney, 1997) carried out the numerical simulations using a quasi-static approach by imposing a symmetry condition for the three canopies in an attempt to predict the equilibrium configuration of a cluster of three half-scaled C-9 parachutes. Stein et al. presented results from the semi-discrete simulations for the aerodynamic interactions between the canopies of parachute clusters of varying numbers and arrangements. In these computations, the canopies were assumed to be rigid. In real life, however, the canopy of a parachute deforms in accordance with aerodynamic forces. The placement of parachutes in cluster configuration does not remain stationary. Therefore, FSI simulations are needed to replicate the real life scenario and to understand the dynamics of cluster of parachutes. Here, we present results from the FSI simulations of a cluster of two G-12 parachutes. The results from a cluster of three G-12 parachutes will be discussed in the next section. The DSD/SST formulations as discussed in the section 0, along with appropriate mesh-update strategies, allowed me to study the

The G-12 is a 64 ft diameter cargo parachute designed to deliver a payload of 2200 lb at a descent speed of approximately 28 ft/s. Clusters of G-12 parachutes are commonly used to deliver larger payloads. The G-12 is constructed with 64 suspension lines which extend to risers. For a single G-12 parachute, the confluence point of the risers is connected to a retraction cable which supports the payload. For a cluster of two G-12 parachutes, the retraction cable is connected to two cables, which connect to the confluence point of the

this objective, multiple parachutes in cluster configurations are often used.

**5. FSI Simulations of parachutes** 

**5.1 Cluster of two G-12 parachutes** 

interaction of canopies in a cluster.

risers for the two G-12 parachutes.

**Problem setup:** 

numerical modeling.

and fluid-structure interface meshes but the canopy mesh for the structural dynamics was very refined to resolve the spatial scales of the gores of the canopy. The same level of mesh refinement for the fluid-structure interface mesh (in case compatible case) was not suitable due to resulting large number elements of the fluid meshes and computationally expensive Navier-Stokes equation solver. Additionally, the length scale associated with fluid dynamics process is smaller than the scale required for structural dynamics to achieve the same level of fidelity. Level of refinement of the meshes were decided based on extensive parachute modeling experiences and should be decided on case-by-case basis. This resulted in incompatible meshes and consequently a mechanism to transfer the deformed coordinate and velocity information from SD to FD and pressure and shear stresses from FD to SD.

Whether one uses different types of elements or different levels of refinements, a projection mechanism is needed to transfer the information from fluid to structure and vice-versa. Errors (usually called projection errors) come along with the projection. Some of the projection techniques are the linear projection and the least-square projection. We used the least-square projection technique to achieve this goal. The least-square projection technique minimizes the sum of residuals in the entire domain. The least-square projection is achieved by solving the following minimization problem

$$\min \int\_{\Gamma^{\ell}} \left\| \mathbf{d}\_{S} - \mathbf{d}\_{F} \right\|^{\ast} d\Gamma \tag{38}$$

where �� (the fluid variables, e.g., pressure drop across the canopy) is projected to �� (the structure variables, e.g., pressure applied on the canopies). Integration is carried over the interface domain Γ� . Similar techniques are used to project the velocity and displacement from structure to the fluid-interface.

#### **4.5 Pressure ramping**

For some cases, the ramping of pressure was absolutely necessary to obtain converged solutions. It was desirable especially in those cases where the guess for the initial shape of the parachute canopies was not close to the real life parachute geometry. The convergence of nonlinear interactions at the start of the FSI simulations was very erratic in the absence of a ramping technique. Pressure ramping was successfully used in these cases. Here, the information exchange between structure and fluid is linearly ramped. The ramping is carried out such that the acting pressure on the canopy is equal to the guess pressure (used for inflating the parachute) at the start of the FSI simulations and equal to the fluid pressure at the end of the ramping. This is achieved by

$$
\Delta P\_S \leftarrow r(t)\Delta P\_G + \left(1 - r(t)\right)\Delta P\_F \tag{39}
$$

where ��� is the effective pressure applied on the structure, ���is the guess pressure used to get the inflated shape of the parachute, and ��� is the pressure drop across the canopies coming from fluid simulations. The ramping factor �(�) is assumed to be linear in time and varies from 0 to 1 in time ����� where ����� is the ramping time. Ramping time were choses case-by-case to achieve appropriate convergence of the Newton-Raphson iterations for the coupled non-linear FSI system.

#### **5. FSI Simulations of parachutes**

250 Fluid Dynamics, Computational Modeling and Applications

and fluid-structure interface meshes but the canopy mesh for the structural dynamics was very refined to resolve the spatial scales of the gores of the canopy. The same level of mesh refinement for the fluid-structure interface mesh (in case compatible case) was not suitable due to resulting large number elements of the fluid meshes and computationally expensive Navier-Stokes equation solver. Additionally, the length scale associated with fluid dynamics process is smaller than the scale required for structural dynamics to achieve the same level of fidelity. Level of refinement of the meshes were decided based on extensive parachute modeling experiences and should be decided on case-by-case basis. This resulted in incompatible meshes and consequently a mechanism to transfer the deformed coordinate and velocity information from SD to FD and pressure and shear

Whether one uses different types of elements or different levels of refinements, a projection mechanism is needed to transfer the information from fluid to structure and vice-versa. Errors (usually called projection errors) come along with the projection. Some of the projection techniques are the linear projection and the least-square projection. We used the least-square projection technique to achieve this goal. The least-square projection technique minimizes the sum of residuals in the entire domain. The least-square projection is achieved

��� � �|�� � ��|�

where �� (the fluid variables, e.g., pressure drop across the canopy) is projected to �� (the structure variables, e.g., pressure applied on the canopies). Integration is carried over the

For some cases, the ramping of pressure was absolutely necessary to obtain converged solutions. It was desirable especially in those cases where the guess for the initial shape of the parachute canopies was not close to the real life parachute geometry. The convergence of nonlinear interactions at the start of the FSI simulations was very erratic in the absence of a ramping technique. Pressure ramping was successfully used in these cases. Here, the information exchange between structure and fluid is linearly ramped. The ramping is carried out such that the acting pressure on the canopy is equal to the guess pressure (used for inflating the parachute) at the start of the FSI simulations and equal to the fluid pressure

where ��� is the effective pressure applied on the structure, ���is the guess pressure used to get the inflated shape of the parachute, and ��� is the pressure drop across the canopies coming from fluid simulations. The ramping factor �(�) is assumed to be linear in time and varies from 0 to 1 in time ����� where ����� is the ramping time. Ramping time were choses case-by-case to achieve appropriate convergence of the Newton-Raphson iterations for the

�

. Similar techniques are used to project the velocity and displacement

��� � �(�)��� � �1 � �(�)���� (39)

�Γ �� (38)

stresses from FD to SD.

interface domain Γ�

**4.5 Pressure ramping** 

from structure to the fluid-interface.

at the end of the ramping. This is achieved by

coupled non-linear FSI system.

by solving the following minimization problem

The parachute canopy is made of very flexible fabric materials. It responds quickly to any small changes in fluid behavior. Consequently, there is a strong interaction between the fluid and structure. Bigger the canopies, stronger will be the interactions. Stronger interactions make the problem very nonlinear in nature and hence increase the difficulty of numerical modeling.

#### **5.1 Cluster of two G-12 parachutes**

A single G-12 parachute is used to drop cargo weighing up to 2,200 lb. There is an interest in the Army for parachute systems that can drop cargos exceeding 2,200 lb. In order to achieve this objective, multiple parachutes in cluster configurations are often used.

The opening stage of a parachute to a fully inflated shape is a critical issue. Usually more than one stage is used to get the fully inflated shape of parachutes in a cluster configuration. One example of multi-stage opening is the use of drogue parachutes. First, a drogue parachute opens and then it pulls the other two G-12 parachutes. After all the parachutes are fully opened, the challenge is to safely drop the cargo into the hands of the right people. The strong parachute-to-parachute aerodynamic interaction can potentially destroy the efficiency of the system if it is not carefully designed. Various researchers have analyzed the cluster of parachute systems using both drop tests and computer simulations. Butler (Butler, 2001) presented the qualitative results from the airdrop tests using a drogue parachute to deploy the G-12 canopies. From the snapshots of a cluster of three G-12 canopy airdrop tests, one clearly observes strong interactions between fluid and canopies. Sahu, et al., (Sahu & Benney, 1997) carried out the numerical simulations using a quasi-static approach by imposing a symmetry condition for the three canopies in an attempt to predict the equilibrium configuration of a cluster of three half-scaled C-9 parachutes. Stein et al. presented results from the semi-discrete simulations for the aerodynamic interactions between the canopies of parachute clusters of varying numbers and arrangements. In these computations, the canopies were assumed to be rigid. In real life, however, the canopy of a parachute deforms in accordance with aerodynamic forces. The placement of parachutes in cluster configuration does not remain stationary. Therefore, FSI simulations are needed to replicate the real life scenario and to understand the dynamics of cluster of parachutes. Here, we present results from the FSI simulations of a cluster of two G-12 parachutes. The results from a cluster of three G-12 parachutes will be discussed in the next section. The DSD/SST formulations as discussed in the section 0, along with appropriate mesh-update strategies, allowed me to study the interaction of canopies in a cluster.

#### **Problem setup:**

The G-12 is a 64 ft diameter cargo parachute designed to deliver a payload of 2200 lb at a descent speed of approximately 28 ft/s. Clusters of G-12 parachutes are commonly used to deliver larger payloads. The G-12 is constructed with 64 suspension lines which extend to risers. For a single G-12 parachute, the confluence point of the risers is connected to a retraction cable which supports the payload. For a cluster of two G-12 parachutes, the retraction cable is connected to two cables, which connect to the confluence point of the risers for the two G-12 parachutes.

Fluid-Structure Interaction Techniques for Parachute 253

start the simulations with a first order-time accurate integration scheme. This helps to clear the start-up vortices quickly. Startup vortices are generated due to large differences in the initial conditions and the exact solutions. After the start-up vortices are cleared from the domain, the second order accurate time-integration scheme is applied. Now, the stage is set to start fluid-structure interaction simulations. Using the results from space-time simulations, a fluid-structure interaction computation is carried out. To remove the mismatch in the initial guessed prestressed configuration, we use the pressure ramping to

The pinned payload case corresponds to the wind tunnel testing where parachutes are fixed to a confluence point. The payload is pinned and is not allowed to move in any direction. The objective here is to understand the aerodynamics of two G-12 flexible canopies. Figure 4 presents the parachute cluster showing the deformed shapes and canopy pressure at different instances of time. Color coding ranging from blue to red represents low to high magnitudes of pressure, respectively. At the time t = 00:27s, two parachute canopies are close to each other. They move away from each other as time progresses. This may be because the initial configuration that I assumed to start the FSI simulations is not in equilibrium. The canopies become after (t = 03:36s) as a result of change in pressure distribution. One notices a strong interaction between the uid and the canopies. Pressure distribution keeps changing with time. Parachute canopies are made of very flexible fabric materials and modeled as a membrane structure in the FSI simulations. As a result, the canopies quickly respond to any change in pressure distribution by adjusting their shape as observed in Figure 4. As a result of dynamic behavior, two canopies come closer to each other and then move farther away. This is a time dependent phenomenon. At time t = 03:36s, these two canopies start going in conical motion in counterclockwise direction about their vertical axis. This motion can be clearly seen in Figs. 4.5 (left and right). They rotate by about 45o in 10:71s. Interestingly, this conical motion has also been observed in real life scenarios. I am not sure how the direction of rotation (clockwise or counter-clockwise) is chosen. I believe that a slight asymmetry in the mesh generated by the automatic mesh generator can give rise to counter-clockwise as being the preferred direction. In Figure 3: Cluster of two G-12 parachutes: Constructed (left) and prestressed (right) structural model. Figure 5 shows the pressure distribution on a plane cutting through the volume mesh and passing through the parachute canopy surfaces. As expected, there is a high pressure region inside of the canopy and lower pressure outside. This pressure gradient keeps the canopy inflated. In Figure 5(middle), a pair of vortices (the blue color spot signifies lower pressure at the core of a vortex) that are shed by the canopy can be seen in the downstream. These pairs of vortices shed from the left and right canopies are of different strengths implying that the aerodynamic responses of parachutes are not symmetric. Figure 5(right) shows a close view of the mesh viewed from the top and colored with pressure at t = 10:74s. Mesh is very \_ne close to the canopy surface. Whole simulations required 10 remeshes. Each remesh resulted in about 3 million elements and 0.5 million nodes. This implies that about 4 million unknown fluid variables were computed every nonlinear iteration in each time step. Total computational time required was about 200 hours using 32 processors on 16 nodes of a Linux cluster (2GB RAM/node,

soften the exchange of information.

1.7GHz P4 Xeon with 1.2GB/s Myrinet Switch).

**Pinned payload:** 

The structural model is composed of membranes, cables, and concentrated masses. The canopy is modeled with triangular membrane elements. Linear cable elements are used to model the suspension lines, radial reinforcements along the canopy, risers, and payload support cables. In each example, the payload is modeled with a single concentrated mass. Material properties are selected to be representative of the G-12. A model for the cluster of G-12 parachutes in constructed (unstressed) and inflated (prestressed) configurations is shown in Figure 3.

Several preparations are required for each fluid-structure interaction simulation. First, a stand-alone structural deformation simulation is carried out to determine the inflated (i.e. prestressed) shape of the G-12. The initial inflation pressure is assumed to be equal to stagnation pressure. From my experience, we observed that this pressure gives a better approximation for the initial shape of the parachute canopies. The prestressed configuration for the G-12 cluster is shown in Figure 3 (right). The unstressed configuration is shown in Figure 3 (left).

Fig. 3. Cluster of two G-12 parachutes: Constructed (left) and prestressed (right) structural model.

Using the parachute canopy from the prestressed configuration, a fluid mesh is generated and a stand-alone fluid simulation is carried out to obtain a developed flow. Fluid simulations are expensive and time consuming. Therefore, to arrive at developed flow quickly, at first semi-discrete (Johnson A. , 1995; Mittal, 1992) formulation is used to compute the flow field. Semi-discrete formulations are first order in time. Using this solution as the initial condition, space-time (Tezduyar T. E., Stabilized finite element formulations for incompressible flow computations, 1992) computations are carried out. We start the simulations with a first order-time accurate integration scheme. This helps to clear the start-up vortices quickly. Startup vortices are generated due to large differences in the initial conditions and the exact solutions. After the start-up vortices are cleared from the domain, the second order accurate time-integration scheme is applied. Now, the stage is set to start fluid-structure interaction simulations. Using the results from space-time simulations, a fluid-structure interaction computation is carried out. To remove the mismatch in the initial guessed prestressed configuration, we use the pressure ramping to soften the exchange of information.

#### **Pinned payload:**

252 Fluid Dynamics, Computational Modeling and Applications

The structural model is composed of membranes, cables, and concentrated masses. The canopy is modeled with triangular membrane elements. Linear cable elements are used to model the suspension lines, radial reinforcements along the canopy, risers, and payload support cables. In each example, the payload is modeled with a single concentrated mass. Material properties are selected to be representative of the G-12. A model for the cluster of G-12 parachutes in constructed (unstressed) and inflated (prestressed) configurations is

Several preparations are required for each fluid-structure interaction simulation. First, a stand-alone structural deformation simulation is carried out to determine the inflated (i.e. prestressed) shape of the G-12. The initial inflation pressure is assumed to be equal to stagnation pressure. From my experience, we observed that this pressure gives a better approximation for the initial shape of the parachute canopies. The prestressed configuration for the G-12 cluster is shown in Figure 3 (right). The unstressed configuration is shown in

Fig. 3. Cluster of two G-12 parachutes: Constructed (left) and prestressed (right) structural

Using the parachute canopy from the prestressed configuration, a fluid mesh is generated and a stand-alone fluid simulation is carried out to obtain a developed flow. Fluid simulations are expensive and time consuming. Therefore, to arrive at developed flow quickly, at first semi-discrete (Johnson A. , 1995; Mittal, 1992) formulation is used to compute the flow field. Semi-discrete formulations are first order in time. Using this solution as the initial condition, space-time (Tezduyar T. E., Stabilized finite element formulations for incompressible flow computations, 1992) computations are carried out. We

shown in Figure 3.

Figure 3 (left).

model.

The pinned payload case corresponds to the wind tunnel testing where parachutes are fixed to a confluence point. The payload is pinned and is not allowed to move in any direction. The objective here is to understand the aerodynamics of two G-12 flexible canopies. Figure 4 presents the parachute cluster showing the deformed shapes and canopy pressure at different instances of time. Color coding ranging from blue to red represents low to high magnitudes of pressure, respectively. At the time t = 00:27s, two parachute canopies are close to each other. They move away from each other as time progresses. This may be because the initial configuration that I assumed to start the FSI simulations is not in equilibrium. The canopies become after (t = 03:36s) as a result of change in pressure distribution. One notices a strong interaction between the uid and the canopies. Pressure distribution keeps changing with time. Parachute canopies are made of very flexible fabric materials and modeled as a membrane structure in the FSI simulations. As a result, the canopies quickly respond to any change in pressure distribution by adjusting their shape as observed in Figure 4. As a result of dynamic behavior, two canopies come closer to each other and then move farther away. This is a time dependent phenomenon. At time t = 03:36s, these two canopies start going in conical motion in counterclockwise direction about their vertical axis. This motion can be clearly seen in Figs. 4.5 (left and right). They rotate by about 45o in 10:71s. Interestingly, this conical motion has also been observed in real life scenarios. I am not sure how the direction of rotation (clockwise or counter-clockwise) is chosen. I believe that a slight asymmetry in the mesh generated by the automatic mesh generator can give rise to counter-clockwise as being the preferred direction. In Figure 3: Cluster of two G-12 parachutes: Constructed (left) and prestressed (right) structural model. Figure 5 shows the pressure distribution on a plane cutting through the volume mesh and passing through the parachute canopy surfaces. As expected, there is a high pressure region inside of the canopy and lower pressure outside. This pressure gradient keeps the canopy inflated. In Figure 5(middle), a pair of vortices (the blue color spot signifies lower pressure at the core of a vortex) that are shed by the canopy can be seen in the downstream. These pairs of vortices shed from the left and right canopies are of different strengths implying that the aerodynamic responses of parachutes are not symmetric. Figure 5(right) shows a close view of the mesh viewed from the top and colored with pressure at t = 10:74s. Mesh is very \_ne close to the canopy surface. Whole simulations required 10 remeshes. Each remesh resulted in about 3 million elements and 0.5 million nodes. This implies that about 4 million unknown fluid variables were computed every nonlinear iteration in each time step. Total computational time required was about 200 hours using 32 processors on 16 nodes of a Linux cluster (2GB RAM/node, 1.7GHz P4 Xeon with 1.2GB/s Myrinet Switch).

Fluid-Structure Interaction Techniques for Parachute 255

Figure 6 shows the shape of the parachutes and their fall in the vertical direction at three instants of time. In this figure, the fall position is determined with respect to an object falling with an assumed freestream velocity of 28 ft/s. The canopies are colored with pressure. With time they fall downward as expected implying that the terminal velocity in this case is higher than the assumed freestream velocity. Here, the cluster of parachutes with a payload of 4,400 lb is cruising at a terminal velocity of around 31 ft/s instead of 28 ft/s. Terminal velocity for a single parachute case is 28 ft/s. As expected, we lost the efficiency in cluster configuration. This implies that the impact velocity for the cluster at the time of landing would be higher than for a single parachute. This is not advisable for the safety of the payload. So, to find the terminal velocity of 28 ft/s in this case, we simulated a few more cases with lower payloads as shown in Figure 7. Following the trend of fall velocity for various payload weights from this figure, one notices that the lighter the payload, the lower the terminal velocity. One observes that the ideal weight for a payload would be about 3,600

t = 06:45 s t = 08:50 s t = 09:87 s

Fig. 6. Airdrop of payload using cluster of two G-12 parachutes: Parachute canopies colored

Fig. 7. Airdrop of a payload using a cluster of two G-12 parachutes: Fall velocity, Vz (ft/s)

v/s time, t (s) for different weight of payload.

lb. This weight will give a terminal velocity of 28 ft/s.

with pressure.

Fig. 4. G-12 parachute cluster with a fixed payload showing the deformed shapes and canopy pressure.

Fig. 5. G-12 parachute cluster with a fixed payload showing pressure (left and middle) and mesh colored with pressure (right) on a cross-section plane cutting through the parachutes.

#### **Airdrop with free payload:**

The pin was removed from the payload and the parachutes were released to fall with the payload weight. The total weight of the payload was 4,400 lb which is twice as heavy as that for a single parachute case. The flow results from a pinned case at t = 06:10s, when the FSI results were fully developed and free from startup conditions, were used as the initial conditions.

Fig. 4. G-12 parachute cluster with a fixed payload showing the deformed shapes and

t = 00:03 s, Top view t = 05:94 s, Side view t = 10:74 s, Top

Fig. 5. G-12 parachute cluster with a fixed payload showing pressure (left and middle) and mesh colored with pressure (right) on a cross-section plane cutting through the parachutes.

The pin was removed from the payload and the parachutes were released to fall with the payload weight. The total weight of the payload was 4,400 lb which is twice as heavy as that for a single parachute case. The flow results from a pinned case at t = 06:10s, when the FSI results were fully developed and free from startup conditions, were used as the

canopy pressure.

**Airdrop with free payload:** 

initial conditions.

Figure 6 shows the shape of the parachutes and their fall in the vertical direction at three instants of time. In this figure, the fall position is determined with respect to an object falling with an assumed freestream velocity of 28 ft/s. The canopies are colored with pressure. With time they fall downward as expected implying that the terminal velocity in this case is higher than the assumed freestream velocity. Here, the cluster of parachutes with a payload of 4,400 lb is cruising at a terminal velocity of around 31 ft/s instead of 28 ft/s. Terminal velocity for a single parachute case is 28 ft/s. As expected, we lost the efficiency in cluster configuration. This implies that the impact velocity for the cluster at the time of landing would be higher than for a single parachute. This is not advisable for the safety of the payload. So, to find the terminal velocity of 28 ft/s in this case, we simulated a few more cases with lower payloads as shown in Figure 7. Following the trend of fall velocity for various payload weights from this figure, one notices that the lighter the payload, the lower the terminal velocity. One observes that the ideal weight for a payload would be about 3,600 lb. This weight will give a terminal velocity of 28 ft/s.

Fig. 6. Airdrop of payload using cluster of two G-12 parachutes: Parachute canopies colored with pressure.

Fig. 7. Airdrop of a payload using a cluster of two G-12 parachutes: Fall velocity, Vz (ft/s) v/s time, t (s) for different weight of payload.

Fluid-Structure Interaction Techniques for Parachute 257

While soft landing systems have been shown to be effective for single parachutes, little is known about the retraction process for a cluster of parachutes. In the this simulation, the soft landing of a 4,400 lb cargo with a cluster of two G-12 parachutes is modeled to study the behavior during and after the retraction process. The soft landing computation is carried out after the parachute has reached a state of terminal descent. In this computation, a 26 ft retraction device is modeled with a cable that connects the confluence point of the two parachutes to the payload. Parachute-payload retraction is modeled by reducing the natural length of the retraction cable by 7.68 ft in 0.27 seconds. Finally, a series of computations are carried out with the retraction cable held at its reduced length after the retraction is completed to study the post retraction behavior of the G-12 cluster. In addition to the soft landing simulation, a second computation is carried out without soft landing retraction. This no-retraction case is a shorter computation and serves as a baseline for comparison with the soft landing case in Figure 8. The aerodynamic drag force on the G-12 cluster during and after soft landing retraction is shown in Fig. 4.8. The effect of the soft landing is apparent from the sharp increase in the drag force during retraction. It is important to note the dramatic drop in drag shortly after retraction ends. This drop is accompanied by canopy collapse and suggests that *harder* landings can be experienced if ground impact is delayed too long after retraction. The payload velocity during descent of the G-12 cluster with and without soft landing retraction is shown in Figure 8 (right). Decreased descent speed resulting from soft landing retraction is very evident, with the payload descent rate decreasing from 31.0 ft/s to 12.0 ft/s. A sequence of snapshots of the two G-12 parachutes during and after the soft landing, colored with the corresponding differential pressures on the canopy is shown in Figure 9. The first two snapshots (left and middle) correspond to times at the start of retraction and at the end of retraction. The final snapshot (right) corresponds to a time well after retraction has finished, with the canopies showing more severe parachute deformations. At the end, before the simulations are stopped, few gores collapsed resulting in gore-to-gore contact.

Further FSI simulations are not possible without a contact model.

Varying numbers of parachutes have been used in cluster configurations by the Army to test the efficiency of each configuration. The most commonly used configurations are clusters of two and three parachutes. A cluster of three parachutes is one of the simplest clusters of parachute configurations. One benefit of this configuration is that this configuration does not have the natural tendency of going in conical motion as observed for a two-parachute cluster. Stein et al. (Stein, et al., 2001; Stein K. , et al., 2003)] presented the results from the stand-alone fluid dynamics simulations of a cluster of three rigid parachutes. The three parachute case presented more challenges in starting the FSI simulations than the two parachute case. Pressure ramping along with smaller time step

A fluid mesh is generated by joining three pieces (domain box, refinement box and parachutes). Several preparation stages similar to the one discussed for a cluster of the two parachute case in the previous section are required. Geometry modeler software package called GAMBIT (Fluent, 2001), was used to model the computational domain and parachute

**5.2 Cluster of three G-12 parachutes** 

was used to overcome the startup problems.

**Problem setup:** 

#### **Soft landing:**

We carried out simulations of two G-12 cargo parachutes to study the soft landing behavior of such systems. *Soft landing* implies making the landing of paratroopers or cargos softer. Various techniques are used to achieve this. Pneumatic Muscle Actuator (PMA) is one technology to achieve this objective. In PMA, the risers are made of pneumatic muscles. When the cargo is very close to the ground the muscle is pressurized. The length of PMA starts contracting and in the process the canopies get pulled downward and the payload gets pulled upward reducing the impact velocity. Previous soft landing simulations focused on the soft landing for single T-10 personnel parachute and on comparisons with drop test data (Stein K. , Tezduyar, Sathe, Benney, & Charles, 2005). These simulations provide a level of confidence for simulations of soft landings with our computational methods. A follow-on simulation was carried out for a single G-12 cargo parachute with a 2,200 lb payload and a G-12 parachute weight of 130 lb (Stein K. , et al., 2003). In this example, soft landing is modeled by reducing the natural length of the retraction cable from 12.80 ft to 2.88 ft in 1 second.

Fig. 8. Drag force (left) and payload velocity with a cluster of two G-12 parachutes during and after soft landing retraction.

Fig. 9. G-12 parachute cluster at the start of pneumatic muscle contraction (left), at the end of contraction (middle) and after soft landing retraction showing the deformed shapes and canopy pressures.

We carried out simulations of two G-12 cargo parachutes to study the soft landing behavior of such systems. *Soft landing* implies making the landing of paratroopers or cargos softer. Various techniques are used to achieve this. Pneumatic Muscle Actuator (PMA) is one technology to achieve this objective. In PMA, the risers are made of pneumatic muscles. When the cargo is very close to the ground the muscle is pressurized. The length of PMA starts contracting and in the process the canopies get pulled downward and the payload gets pulled upward reducing the impact velocity. Previous soft landing simulations focused on the soft landing for single T-10 personnel parachute and on comparisons with drop test data (Stein K. , Tezduyar, Sathe, Benney, & Charles, 2005). These simulations provide a level of confidence for simulations of soft landings with our computational methods. A follow-on simulation was carried out for a single G-12 cargo parachute with a 2,200 lb payload and a G-12 parachute weight of 130 lb (Stein K. , et al., 2003). In this example, soft landing is modeled by reducing

the natural length of the retraction cable from 12.80 ft to 2.88 ft in 1 second.

Fig. 8. Drag force (left) and payload velocity with a cluster of two G-12 parachutes during

t = 11:11 s t = 11:38 s t = 13:10 s

Fig. 9. G-12 parachute cluster at the start of pneumatic muscle contraction (left), at the end of contraction (middle) and after soft landing retraction showing the deformed shapes and

**Soft landing:** 

and after soft landing retraction.

canopy pressures.

While soft landing systems have been shown to be effective for single parachutes, little is known about the retraction process for a cluster of parachutes. In the this simulation, the soft landing of a 4,400 lb cargo with a cluster of two G-12 parachutes is modeled to study the behavior during and after the retraction process. The soft landing computation is carried out after the parachute has reached a state of terminal descent. In this computation, a 26 ft retraction device is modeled with a cable that connects the confluence point of the two parachutes to the payload. Parachute-payload retraction is modeled by reducing the natural length of the retraction cable by 7.68 ft in 0.27 seconds. Finally, a series of computations are carried out with the retraction cable held at its reduced length after the retraction is completed to study the post retraction behavior of the G-12 cluster. In addition to the soft landing simulation, a second computation is carried out without soft landing retraction. This no-retraction case is a shorter computation and serves as a baseline for comparison with the soft landing case in Figure 8. The aerodynamic drag force on the G-12 cluster during and after soft landing retraction is shown in Fig. 4.8. The effect of the soft landing is apparent from the sharp increase in the drag force during retraction. It is important to note the dramatic drop in drag shortly after retraction ends. This drop is accompanied by canopy collapse and suggests that *harder* landings can be experienced if ground impact is delayed too long after retraction. The payload velocity during descent of the G-12 cluster with and without soft landing retraction is shown in Figure 8 (right). Decreased descent speed resulting from soft landing retraction is very evident, with the payload descent rate decreasing from 31.0 ft/s to 12.0 ft/s. A sequence of snapshots of the two G-12 parachutes during and after the soft landing, colored with the corresponding differential pressures on the canopy is shown in Figure 9. The first two snapshots (left and middle) correspond to times at the start of retraction and at the end of retraction. The final snapshot (right) corresponds to a time well after retraction has finished, with the canopies showing more severe parachute deformations. At the end, before the simulations are stopped, few gores collapsed resulting in gore-to-gore contact. Further FSI simulations are not possible without a contact model.

#### **5.2 Cluster of three G-12 parachutes**

Varying numbers of parachutes have been used in cluster configurations by the Army to test the efficiency of each configuration. The most commonly used configurations are clusters of two and three parachutes. A cluster of three parachutes is one of the simplest clusters of parachute configurations. One benefit of this configuration is that this configuration does not have the natural tendency of going in conical motion as observed for a two-parachute cluster. Stein et al. (Stein, et al., 2001; Stein K. , et al., 2003)] presented the results from the stand-alone fluid dynamics simulations of a cluster of three rigid parachutes. The three parachute case presented more challenges in starting the FSI simulations than the two parachute case. Pressure ramping along with smaller time step was used to overcome the startup problems.

#### **Problem setup:**

A fluid mesh is generated by joining three pieces (domain box, refinement box and parachutes). Several preparation stages similar to the one discussed for a cluster of the two parachute case in the previous section are required. Geometry modeler software package called GAMBIT (Fluent, 2001), was used to model the computational domain and parachute

Fluid-Structure Interaction Techniques for Parachute 259

downstream direction such that the equivalent freestream Mach number is 2.8. The computed values of Mach cone angles (�) with the analytical values are compared. As observed in the table from Figure 11, there is a good agreement between these two results.

> *Fixed wedge*

*Moving Wedge* 

54.790 52.420 Analytical (�) 54.520 52.000 Computed (�)

Fig. 11. Mesh with wedge (left), density at the beginning and end of the computations (center) and computed Mach cone angle (�) compared for fixed and moving wedge with

In this chapter, we reviewed some computational challenges in fluid structure techniques. FSI simulations consisted of the following primary components: a solution method for the fluid dynamics, a solution method for the structural dynamics, an intelligent way to efficiently move and update a finite element mesh of fluid, and strategies for the coupling of fluid, and structural dynamics along the fluid-structure interface. The coupling is achieved in staggered fashion, with the fluid and structure coupled iteratively within a nonlinear iteration loop, and with multiple nonlinear iterations improving the convergence of the coupled system. The mesh domain was treated as a pseudo-elastic solid. To handle large structural-deformations, fluid dynamics solver required a methodology which can handle time-dependent spatial domain deformations. We presented Deforming-Spatial Domain / Stabilized Space-Time (DSD/SST) method which have a built-in capability to handle deforming structures. Semi-discrete formulations were used to carry out structural dynamics simulations which are based on the principle of virtual work. The SD model usually consists of membranes along with cables and

We presented results from the FSI simulations of a cluster of two G-12 parachutes under three for the parachute airdrop: fixed payload, airdrop, and soft-landing and cluster of three parachutes under fixed load conditions. Here, the parachute is represented as a structure composed of membranes, cables, and concentrated masses. The cables and membranes are assumed to have no flexural rigidity and experience large displacements and rotations. As a result, the interaction between a parachute system and the surrounding flow field is dominant in most of the parachute operations. Thus, the ability to predict parachute FSI is a challenge that must be faced in airdrop systems modeling. The G-12 parachute is a very large cargo parachute and its analysis presents convergence difficulties at the start of FSI simulations. A pressure ramping technique was proposed to deal with this start-up mismatch in the information exchange between the parachute's canopy and the fluid. It was found that a pressure ramping technique to smoothly transfer

analytical solution.

payloads.

**6. Concluding remarks** 

geometry and to generate the surface meshes. The volume mesh is generated by an unstructured automatic mesh generator (Johnson A. , 1995; Johnson & Tezduyar, Mesh update strategies in parallel finite element computations of flow problems with moving boundaries and interfaces, 1994; Johnson & Tezduyar, Advanced mesh generation and update methods for 3D flow simulations, 1999). A refinement boundary is used to get a satisfactory level of mesh refinements near the canopy's boundary to capture the flow physics accurately. The automatic mesh generator creates about 6 million tetrahedral elements and 1 million nodes of fluid mesh, resulting in approximately 8 million unknowns per nonlinear iteration at each time step with DSD/SST formulations. The size of the structure mesh is negligible compared to the fluid mesh. The Reynolds number, based on radius of a flat G-12 parachute, is 5.5 million. We used the parallel FSI solver on 128 Cray T3E processors to arrive at the results. The freestream velocity, which corresponds to terminal velocity, is assumed to be 20 ft/s.

#### **Pinned payload:**

The computed results from the FSI simulations are shown in Figure 10. We notice that the symmetric distribution is lost as a result of the dynamic behavior of parachutes in this cluster configuration. The changes in pressure distribution on the canopies surfaces result in stronger parachute-to-parachute aerodynamic interactions. In fact, two of the parachutes came very close to each other at the end of simulations. The closer they got, the higher the frequency of remeshes were and the more difficult it became to perform FSI simulations. At the onset of contact, the simulations have to be terminated in the absence of a contact model.

Fig. 10. Cluster of three G-12 parachutes at three different instances of time. Canopies are colored with pressure.

#### **5.3 An example for compressible flow fluid-structure interaction simulation**

Our targeted goal is to study the behavior of parachute dynamic in compressible flow regime (i.e., flow with Mach number>0.3). Currently, the method has been tested for simple validation problems. Here, we present FSI results for a moving wedge problem in supersonic flow regime. The Mach number of the incoming fluid is 3.

At first, a solution for *fixed-wedge* case is arrived before switching moving algorithms. The *fixed-wedge* is used as the initial solution for the *moving-wedge*. The wedge is moving in the

geometry and to generate the surface meshes. The volume mesh is generated by an unstructured automatic mesh generator (Johnson A. , 1995; Johnson & Tezduyar, Mesh update strategies in parallel finite element computations of flow problems with moving boundaries and interfaces, 1994; Johnson & Tezduyar, Advanced mesh generation and update methods for 3D flow simulations, 1999). A refinement boundary is used to get a satisfactory level of mesh refinements near the canopy's boundary to capture the flow physics accurately. The automatic mesh generator creates about 6 million tetrahedral elements and 1 million nodes of fluid mesh, resulting in approximately 8 million unknowns per nonlinear iteration at each time step with DSD/SST formulations. The size of the structure mesh is negligible compared to the fluid mesh. The Reynolds number, based on radius of a flat G-12 parachute, is 5.5 million. We used the parallel FSI solver on 128 Cray T3E processors to arrive at the results. The freestream velocity, which corresponds to

The computed results from the FSI simulations are shown in Figure 10. We notice that the symmetric distribution is lost as a result of the dynamic behavior of parachutes in this cluster configuration. The changes in pressure distribution on the canopies surfaces result in stronger parachute-to-parachute aerodynamic interactions. In fact, two of the parachutes came very close to each other at the end of simulations. The closer they got, the higher the frequency of remeshes were and the more difficult it became to perform FSI simulations. At the onset of contact, the simulations have to be terminated in the absence

Fig. 10. Cluster of three G-12 parachutes at three different instances of time. Canopies are

Our targeted goal is to study the behavior of parachute dynamic in compressible flow regime (i.e., flow with Mach number>0.3). Currently, the method has been tested for simple validation problems. Here, we present FSI results for a moving wedge problem in

At first, a solution for *fixed-wedge* case is arrived before switching moving algorithms. The *fixed-wedge* is used as the initial solution for the *moving-wedge*. The wedge is moving in the

**5.3 An example for compressible flow fluid-structure interaction simulation** 

supersonic flow regime. The Mach number of the incoming fluid is 3.

terminal velocity, is assumed to be 20 ft/s.

**Pinned payload:** 

of a contact model.

colored with pressure.

downstream direction such that the equivalent freestream Mach number is 2.8. The computed values of Mach cone angles (�) with the analytical values are compared. As observed in the table from Figure 11, there is a good agreement between these two results.

Fig. 11. Mesh with wedge (left), density at the beginning and end of the computations (center) and computed Mach cone angle (�) compared for fixed and moving wedge with analytical solution.

### **6. Concluding remarks**

In this chapter, we reviewed some computational challenges in fluid structure techniques. FSI simulations consisted of the following primary components: a solution method for the fluid dynamics, a solution method for the structural dynamics, an intelligent way to efficiently move and update a finite element mesh of fluid, and strategies for the coupling of fluid, and structural dynamics along the fluid-structure interface. The coupling is achieved in staggered fashion, with the fluid and structure coupled iteratively within a nonlinear iteration loop, and with multiple nonlinear iterations improving the convergence of the coupled system. The mesh domain was treated as a pseudo-elastic solid. To handle large structural-deformations, fluid dynamics solver required a methodology which can handle time-dependent spatial domain deformations. We presented Deforming-Spatial Domain / Stabilized Space-Time (DSD/SST) method which have a built-in capability to handle deforming structures. Semi-discrete formulations were used to carry out structural dynamics simulations which are based on the principle of virtual work. The SD model usually consists of membranes along with cables and payloads.

We presented results from the FSI simulations of a cluster of two G-12 parachutes under three for the parachute airdrop: fixed payload, airdrop, and soft-landing and cluster of three parachutes under fixed load conditions. Here, the parachute is represented as a structure composed of membranes, cables, and concentrated masses. The cables and membranes are assumed to have no flexural rigidity and experience large displacements and rotations. As a result, the interaction between a parachute system and the surrounding flow field is dominant in most of the parachute operations. Thus, the ability to predict parachute FSI is a challenge that must be faced in airdrop systems modeling. The G-12 parachute is a very large cargo parachute and its analysis presents convergence difficulties at the start of FSI simulations. A pressure ramping technique was proposed to deal with this start-up mismatch in the information exchange between the parachute's canopy and the fluid. It was found that a pressure ramping technique to smoothly transfer

Fluid-Structure Interaction Techniques for Parachute 261

Hirsch, C. (1988). *Numerical Computation of Internal and External Flows: Fundamentals of* 

Hughes, T. J., Franca, L. P., & Mallet, M. (1986). A new finite element formulation for

Johnson, A. (1995). Mesh Generation and Update Strategies for Parallel Computation of

Johnson, A. A., & Tezduyar, T. E. (1994). Mesh update strategies in parallel finite element

Johnson, A. A., & Tezduyar, T. E. (1999). Advanced mesh generation and update methods

Kumar, V. (2005). Advanced Computational Techniques for Incompressible / Compressible

Le Beau, G., & Tezduyar, T. (1991). Finite element computation of compressible flows with

Mittal, S. (1992). *Stabilized Space-Time Finite Element Formulations for Unsteady Incompressible* 

Sahu, J., & Benney, R. (1997). *Prediction of terminal decent characteristics of parachute clusters* 

Shakib. (1988). *Finite Element Analysis of the Compressible Euler and Navier-Stokes Equations.*

Stein, K. (1999). Simulation and Modeling Techniques for Parachute Fluid-Structure Interactions. *Aerospace Engineering, University of Minnesota, Ph.D. thesis*. Stein, K., Benney, R., Tezduyar, T., Kumar, V., Thornburg, E., Kyle, C., et al. (2001).

Stein, K., Tezduyar, S., Sathe, S., Senga, M., Ozcan, C., Soltys, T., et al. (2003). Simulation of

Stein, K., Tezduyar, T. E., Sathe, S., Benney, R., & Charles, R. (2005). Fluid-structure

Stein, K., Tezduyar, T., Kumar, V., Sathe, S., Benney, R., Thornburg, E., et al. (2003).

Tezduyar, T. E. (1992). Stabilized finite element formulations for incompressible flow

Tezduyar, T. E., Behr, M., & Liou, J. (1992). A new strategy for finitie element computations

Aerodynamic interaction between multiple parachute canopies. *Proceedings of the* 

parachute dynamics during control line input operations. *Proceedings of the 17th AIAA Aerodynamic Decelerator Systems Technology Conference.* Monterrey, California.

interaction modelling of parachute soft-landing dynamics. *International Journal for* 

Aerodynamic interactions between parachute canopies. *Journal of Applied* 

involving moving boundaries and interfaces - the deforming-spatial-domain/space time procedure: I. the concept and the preliminary tests. *Computer Methods in* 

computational fluid dynamics: I. symmetric forms of compressible Euler and Navier-Stokes equations and the second law of thermodynamics. *Computer Methods* 

Flow Problems with Moving Boundaries and Interaces. *Aerospace Engineering,* 

computations of flow problems with moving boundaries and interfaces. *Computer* 

the SUPG formulation. In *Advances in Finite Element Analysis in Fluid Dynamics*

*Numerical Discretization* (Vol. 1). Chichester: John Wiley and Sons.

Fluent, I. (2001). *GAMBIT 2, CFD Preprocessor.* Lebanon, NH: Fluent, Inc.

*in Applied Mechanics and Engineering, 54*, 223-234.

*Methods in Applied Mechanics and Engineering, 119*, 73-94.

(Vols. FED-Vol-123, pp. 21-27). New York: ASME.

*using CFD.* Technical report, AIAA Paper No. 97-1453.

Ph.D. thesis, Stanford University, Mechanical Engineering.

*First MIT Conference on Computational Fluid and Solid.* Boston.

*Numerical Methods in Fluids, 47*, 619-631.

computations. *Advances in Applied Mechanics*, 1-44.

*Applied Mechanics and Engineering*, 339-351.

*Mechanics, 70*, 50-57.

for 3D flow simulations. *Computational Mechanics, 23*, 130-143.

Fluid-Structure Interactions. Rice University: Rice University.

*Flows Involving Fluid-Body Interactions.* University of Minnesota.

*University of Minnesota, Ph.D. thesis*.

the information between the fluid and structure solves the convergence problem encountered at the start of FSI simulations.

It was shown that one can ideally achieve an impact velocity as small as 12 ft/s using a softlanding technique. However, for the 2-parachute cluster system, this minimum impact velocity is only achieved at a certain time after contraction. If one waits too long after contraction, the descent speed begins to increase, again, due to canopy collapse creating a harder landing. The amount of contraction is limited by the canopy strength because the loading on the canopy increases throughout the contraction process.

In a 3-parachute cluster system, we were able to predict that two parachutes collide as we observe in drop tests. We believe that this is caused by unsteady pressure distribution on the canopy surface as depicted in Figure 10. One can devise a control mechanism to prevent the collision of two canopies. One idea is to experiment with different control mechanisms such as heat-induced porosity-altering techniques for selected gores to introduce a sideways force.

The results from the FSI simulations of a cluster of three G-12 parachutes were also presented. The FSI simulations of both clusters of two and three parachute systems were found to suffer from contact issues. Two of the parachutes in this case came very close to contact. The FSI simulations faced a numerical barrier on the onset of contact. Similar contact issues were observed in the case of two parachute simulations when one or more gores collapsed. Further FSI simulations were not possible in the absence of a contact model.

In future research, we need to improve the structural modeling capabilities for the cables, risers, and membranes by instantaneously addressing the non-isotropic and non-linear deformation of the flexible structures caused by changes in the fluid dynamics forces to imitate the physical reality. Additionally, the addition of a contact model will enable the formulations to continue after contact is made to simulate what happens physically when contact is made while the parachutes are still falling through the air. All such additions to the code will greatly increase the cost, so increased accuracy in the model simulations depends on computational resources and continued innovation in clock speeds and floating point operations per second (FLOPS). Another route of research to address the contact is to experiment with different control mechanisms to prevent contact.

#### **7. References**

Babuska, I. (1973). The finite element method with Lagrange multipliers. *Numerische Mathematik, 20*, 179-192.


the information between the fluid and structure solves the convergence problem

It was shown that one can ideally achieve an impact velocity as small as 12 ft/s using a softlanding technique. However, for the 2-parachute cluster system, this minimum impact velocity is only achieved at a certain time after contraction. If one waits too long after contraction, the descent speed begins to increase, again, due to canopy collapse creating a harder landing. The amount of contraction is limited by the canopy strength because the

In a 3-parachute cluster system, we were able to predict that two parachutes collide as we observe in drop tests. We believe that this is caused by unsteady pressure distribution on the canopy surface as depicted in Figure 10. One can devise a control mechanism to prevent the collision of two canopies. One idea is to experiment with different control mechanisms such as heat-induced porosity-altering techniques for selected gores to

The results from the FSI simulations of a cluster of three G-12 parachutes were also presented. The FSI simulations of both clusters of two and three parachute systems were found to suffer from contact issues. Two of the parachutes in this case came very close to contact. The FSI simulations faced a numerical barrier on the onset of contact. Similar contact issues were observed in the case of two parachute simulations when one or more gores collapsed. Further FSI simulations were not possible in the absence of a contact

In future research, we need to improve the structural modeling capabilities for the cables, risers, and membranes by instantaneously addressing the non-isotropic and non-linear deformation of the flexible structures caused by changes in the fluid dynamics forces to imitate the physical reality. Additionally, the addition of a contact model will enable the formulations to continue after contact is made to simulate what happens physically when contact is made while the parachutes are still falling through the air. All such additions to the code will greatly increase the cost, so increased accuracy in the model simulations depends on computational resources and continued innovation in clock speeds and floating point operations per second (FLOPS). Another route of research to address the contact is to

Babuska, I. (1973). The finite element method with Lagrange multipliers. *Numerische* 

Behr, M., & Tezduyar, T. (1994). Finite element solution strategies for large-scale flow simulations. *Computer Methods in Applied Mechanics and Engineering, 112*, 3-24. Brezzi, F. (1974). On the existence, uniqueness and approximation of saddle-point problems arising from Lagrangian multipliers. *Modelisation Math. Anal. Number, 8*, 129-151. Brooks, A., & Hughes, T. (1982). Streamline upwind/Petrov-Galerkin formulations for

Butler, M. C. (2001). Additional applications of bat sombrero slider technology. *Proceedings* 

*of AIAA 16th Aerodynamic Decelerator Systems Technology Conference.*

convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations. *Computer Methods in Applied Mechanics and Engineering, 32*,

loading on the canopy increases throughout the contraction process.

experiment with different control mechanisms to prevent contact.

encountered at the start of FSI simulations.

introduce a sideways force.

model.

**7. References** 

199-259.

*Mathematik, 20*, 179-192.

Fluent, I. (2001). *GAMBIT 2, CFD Preprocessor.* Lebanon, NH: Fluent, Inc.


**Part 3** 

**Heat Transfer, Combustion, and Energy** 

