Extension of a sharp-interface immersed-boundary method for simulating parachute inflation

In this work, the sharp-interface immersed boundary (IB) method proposed by Mittal et al. (J Comput Phys 227(10):4825–4852, 2008) is extended to fluid-structure-inter-action (FSI) simulation of parachute inflation by utilizing several open-source tools. The method employs a Cartesian-grid ghost-cell methodology to accurately represent the immersed boundary, and it is suitable for solving moving-boundary flows with arbitrarily complex geometries. The finite-element code CalculiX is employed to solve the structural dynamics of the parachute system. The IB flow solver is coupled with CalculiX in a minimally-invasive manner using the multi-physics coupling library preCICE. The implicit fluid-structure coupling together with the Aitken adaptive under-relaxation scheme is considered to improve the numerical accuracy and stability. The developed approach is validated by a benchmark FSI case. Numerical experiments on the inflation process of several typical parachutes are further conducted. The breathing process, flow structure, canopy displacement and drag coefficient are analyzed to demonstrate the applicability of the present approach for simulating para-chute inflation.


Introduction
Parachute inflation plays a critical role in aerospace applications, particularly in the field of aeronautics and space exploration.Parachute inflation involves the complex fluidstructure interaction (FSI).High fidelity simulation of this interaction could be computationally expensive and challenging.
Numerical methods for simulating parachute inflation have undergone significant development over the years.Arbitrary Lagrangian-Eulerian (ALE) methods combine the advantages of both Eulerian and Lagrangian approaches, enabling the mesh to move and deform with the fluid flow while still retaining the benefits of a fixed-grid Eulerian approach, and therefore they have been widely applied for simulating large deformations of the parachute fabric [1][2][3][4].The deforming spatial domain/stabilized space-time (DSD/SST) method is a numerical technique developed by the T*AFSM team [5][6][7][8] for simulating parachute inflation.This method combines the advantages of the deforming spatial domain approach and the stabilized space-time formulation to accurately model and capture the dynamic behavior of the parachute during inflation.
The immersed boundary (IB) method [9] is also an attractive numerical method for FSI simulation of parachute inflation because it does not require dynamic generation of body-fitted meshes, making it suitable for tackling problems involving large structural deformations or displacements.Depending on the representation of the fluid-structure interface, the IB method can be classified into the diffuse-interface and sharp-interface categories.On one hand, the diffuse-interface IB method represents the fluid-structure interface as a diffuse or smeared region, which provides flexibility and computational convenience in handling complex geometries and deformations.Kim and Peskin [10,11] used the diffuse-interface IB method to study the semi-opened parachute in both two and three dimensions.Liu et al. [12] combined the diffuse-interface IB method with the lattice Boltzmann method to simulate the inflation process of several typical parachute systems.The authors developed a computational framework for parachute inflation based on the diffuse-interface IB method using the open-source IBAMR library in their previous work [13].On the other hand, the sharp-interface IB method represents the fluid-structure interface as a sharp boundary, which offers more accurate tracking and representation of well-defined interfaces.Boustani et al. [14,15] developed a high-fidelity FSI solver based on the sharp-interface IB method and the finite element method to simulate the supersonic parachute inflation phase of the ASPIRE SR01 flight test.Nevertheless, application of the IB method to parachute inflation is still rare compared to the ALE and DSD/SST methods.
Simulation of FSI problems can be performed with the monolithic or partitioned approaches.The monolithic approach [16] solves the fluid and structure equations as a single coupled system which allows for capturing the strong interactions between the fluid and structure.This approach offers inherent numerical accuracy and stability, but may be computationally expensive.The partitioned approach [17] treats the fluid and structure domains as distinct subproblems, and their respective equations are solved sequentially or in an alternating fashion.This approach offers computational efficiency as the fluid and structure solvers can be optimized independently and existing software can be leveraged for each subproblem through a coupling interface, but it may lead to stability issues when the added mass effect of the FSI system is strong.The partitioned approach is more popular than the monolithic approach because of its flexibility and computational efficiency.Several good coupling tools such as MpCCI [18], OpenPALM [19], and ADVENTURE_Coupler [20] have been developed and served as a coupling interface for partitioned FSI simulations.Recently, the open-source coupling tool pre-CICE has drawn a lot of attention because it can offer the complete coupling functionality required for fast development of a multi-physics environment using existing, single-physics, and black-box solvers with minimal effort [21].Thanks to these attractive features, it has been applied to various partitioned FSI simulations [22][23][24].
In this work, based on the idea of partitioned FSI, the sharp-interface IB method proposed by Mittal et al. [25] is extended to simulate parachute inflation by utilizing open-source tools.This IB method is based on a ghost-cell approach which can handle arbitrarily complex geometries.The cells whose nodes are inside the solid but have at least one neighbor in the fluid are referred to as ghost cells.The flow variables at these ghost cells are obtained by implicit satisfaction of the boundary conditions on the immersed boundary.The open-source finite-element code CalculiX (Version 2.19) [26] is employed to solve the nonlinear large deformation of the parachute system.The opensource multi-physics coupling library preCICE (Version 2.3.0)[27] is used to handle the interfacing between the in-house IB solver and CalculiX.To overcome the strong added mass effect induced by the thin and light canopy fabric, the implicit-coupling scheme which introduces sub-iterations within each time step is adopted.
The remainder of this paper is organized as follows.Section 2 reviews the sharp-interface IB method for solving the fluid flow.In Section 3, extension of the IB method for FSI simulations of parachute inflation is elaborated.Section 4 introduces the computational setup for the fluid and structure parts, respectively.Section 5 conducts a benchmark case to validate the present approach.In Section 6, the inflation process of several typical parachutes including the square parachute, cross parachute, double cross parachute, and round parachute are simulated to demonstrate the capability of the present approach.Section 7 concludes this study and shows some prospects for future work.

Sharp-interface IB method
The goal of this work is to extend the sharp-interface IB method [25] for simulating parachute inflation, and therefore a review of the detailed numerical method is given in this section.
The governing equations for the air flow are the 3D unsteady incompressible Navier-Stokes equations, where u i , ρ f , p f , ν f represent the velocity, density, pressure, and kinematic viscosity, respectively.The equations are discretized in space using a cell-centered collocated arrangement of the primitive variables.The three-substep fractional step method [28] is used to integrate the equations in time.In the first substep, a modified momentum equation is solved and an intermediate velocity is obtained.A second-order Adams-Bashforth scheme is employed for the convective terms while the implicit Crank-Nicholson scheme is used to discretize the diffusion terms to eliminate the viscous stability constraint.The second substep requires the solution of the pressure correction for the Poisson equation, which is solved with a highly efficient geometric multigrid method.Once the pressure correction is obtained, the pressure and velocity are updated in the third substep.
In the sharp-interface IB method, construction of the interpolation stencil near the boundary requires identifying whether a mesh node belongs to the fluid or solid phase.As shown in Fig. 1, the fluid and solid cells are first identified, and then the ghost cells whose centers are inside the solid but have at least one neighbor in the fluid are identified.The flow variables at the ghost cells are not known a priori, and (1a) therefore an image point in the fluid phase is constructed, which is the mirror of the corresponding ghost-cell center along the surface normal.The flow variables at the image point are evaluated by a trilinear interpolation between the surrounding fluid nodes.The midway point between the ghost-cell center and the image point is referred to as the boundary intercept point, where the boundary conditions can be precisely satisfied.After that, the flow variables at the ghost-cell can be obtained by extrapolation between the corresponding image point and the body intercept point.Note that the thickness of the canopy membrane is usually much smaller than its characteristic length, making it difficult to realize the above identification procedure.Therefore, similar to the treatment in [29], an artificial thickness is introduced along the canopy membrane to expand the 2D membrane into a 3D body.Figure 2 shows the 2D schematic of the boundary treatment near the canopy membrane.By specifying an artificial thickness h which is around three cells wide, three different types of cells can be identified according to the position of the expanded geometry, i.e., fluid cells, solid cells, and ghost cells which are one layer inside of the expanded geometry.The boundary intercept point on the Fig. 1 2D schematic of the interpolation stencil of the sharp-interface immersed boundary method Fig. 2 2D schematic of the boundary treatment near the canopy membrane expanded geometry is found by extending the fluid node or ghost node in a direction normal to the physical surface of the geometry.The boundary conditions at the boundary intercept points on the expanded geometry are assumed to be the same as those at the physical geometry.

FSI implementation
In this work, the in-house sharp-interface IB solver is coupled with an open-source finite-element solver CalculiX [26] for FSI simulation of the parachute inflation.An open-source coupling library preCICE [27] is employed for interfacing between the IB solver and CalculiX.Details about the implementation of the FSI approach are presented in this section.A parachute system typically consists of a canopy membrane and several suspension lines.Since the deformation of the suspension lines has little effect on the surrounding fluid, only the canopy membrane is involved in the two-way FSI coupling.

Structure solver
The structural dynamics of the parachute is governed by the Navier equation with a linear stress-strain relationship, where ρ s is the structure density, d the displacement, σ the stress tensor, and b 0 the body force.
The stress tensor is expressed by where E = 1 2 (∇ T d + ∇d) is the Green-Lagrange strain tensor, ( , µ) the Lamé con- stants, and I the second rank identity tensor.
The structure equation is discretized with the finite element method in CalculiX.Direct integration dynamic analysis is selected as the analysis type, and the equations of motion are numerically integrated using the α-method [26].

Coupling using preCICE library
The preCICE library offers parallel communication means, data mapping schemes, and methods for transient equation coupling (explicit or implicit) [27].With these features, preCICE can couple existing open-source codes or in-house codes in a minimally-invasive manner by using adapters.An adapter refers to a software component or module that enables the coupling of different solvers.The purpose of an adapter is to provide an interface between the preCICE library and a specific simulation code, allowing the exchange of data and the synchronization of the simulation processes.A schematic coupling of the IB solver with CalculiX using preCICE is shown in Fig. 3.Note that the Cal-culiX adapter is officially provided by preCICE, while the IB adapter is built by following (2) the instructions on the website of preCICE (https:// preci ce.org/).Detailed implementation of the coupling using the preCICE library is discussed as follows.

Data mapping
There are several interpolation-based data mapping methods provided in preCICE for information exchange between the fluid and solid subdomains.Since the canopy membrane can be treated as a 2D manifold embedded in the 3D space, an accurate and straightforward data mapping method is developed in this work.Taking the canopy membrane of a round parachute as an example, a nodal mapping method between the fluid and solid is shown in Fig. 4. First, on the solid phase, the canopy surface s is discretized with the quadratic triangular shell element in CalculiX, and the node numbering (1-2-3-4-5-6 in the counter-clockwise direction) for each element is also shown in a zoomed-in view.Then, the same node is used and each shell element is further split into four triangular elements (1-4-6, 2-5-4, 3-6-5, 4-5-6 in the counter-clockwise direction) to discretize the canopy surface f on the fluid phase.In this way, the solid and fluid meshes at the interface can be exactly mapped without involving any interpolation error.Although this data mapping method requires a relatively dense mesh for the solid subdomain, the total computational time would not increase significantly because the main contribution of the computational time comes from the fluid subdomain.

Coupling scheme
For a typical canopy membrane with characteristic length L and thickness h, the mass ratio between the solid and fluid is m * = ρ s h/ρ f L ≪ 1 , which will cause significant added-mass effect [36], and therefore the following implicit-coupling scheme with under-relaxation iteration is adopted, such that the numerical accuracy and stability of FSI coupling at each time step can be improved.
where φ represents the pressure, velocity or displacement at the fluid-structure interface, φ is the value predicted by the structure solver, k is the iteration number, and α is the relaxation factor between 0 and 1.To improve the convergence rate of implicit coupling, the Aitken acceleration technique [36] is introduced to dynamically determine the value of α.

Intercommunication
Throughout the simulation, intercommunication between the fluid and solid subdomains, such as parallel data transfer and synchronization, is efficiently handled by the preCICE library.
The procedure of implicit coupling between the IB solver and CalculiX with the help of the preCICE library is summarized in Fig. 5.A step-by-step explanation is as follows: (4) Fig. 5 Procedure of implicit FSI coupling with the preCICE library (1) At time step n, the IB solver computes forces on the fluid subdomain.
(2) Fluid forces are mapped onto the solid surface mesh.
(3) preCICE communicates the force data to CalculiX.(4) CalculilX computes the displacement with the mapped boundary conditions on the solid subdomain.(5) preCICE communicates the displacements back to the IB solver.
Steps (1)-( 5) are repeated with the Aitken adaptive under-relaxation until convergence is reached.After that, the IB solver maps back the received displacement field and computes new forces according to the updated fluid-solid interface at time step n + 1 .The FSI simulation stops when the specified time-step number N max is reached.

Fluid
As illustrated in Fig. 6, the size of the fluid domain is 10L × 5L × 5L in the x, y, and z directions, where L is the characteristic length of the canopy and x is the direction of the incoming flow.The computational domain is respectively discretized with 128, 64, and 64 nodes in the x, y, and z directions.The canopy membrane is placed in a box with dimensions 1.5L × 1.25L × 1.25L in the x, y, and z directions, respectively.The box is discretized with a uniform mesh of 0.02L, and gradually stretched to the far field in the x, y, and z directions.Previous studies have shown that this mesh resolution is adequate for FSI simulation of a membrane [31,37].The left side of the domain is the inlet with a uniform velocity U ∞ .The right side of the domain is the outlet with the Neumann velocity boundary condition.No-slip velocity boundary conditions are specified at other sides of the domain.A zero reference pressure is specified at the outlet, and all other boundaries are specified with the Neumann pressure boundary conditions.
The Reynolds number is defined as Re = U ∞ L/ν f .Due to the challenge of IB meth- ods for simulating turbulent flows [38], the Reynolds number in this study is reduced to Re = 1 × 10 3 without considering the turbulence effect.Similar treatment has been adopted in previous studies [11,39,40].The time-step size is chosen as 2 × 10 −3 L/U ∞ , and the cor- responding maximum Courant-Friedrichs-Lewy (CFL) number is approximately 0.3 − 0.4.

Structure
In this work, we consider simple parachute systems which consist of the canopy membrane, suspension lines, and reinforcement belts (when necessary).As illustrated in Fig. 7, we use Gmsh (Version 4.1) [41] and several Python scripts to develop an automatic procedure for preparing the ingredients required by CalculiX, namely, the mesh, boundary condition, and material property for each component of the parachute.Note that this procedure is versatile for any type of parachute system.
Each component of the parachute is discretized with an appropriate element type, i.e., shell element for the canopy membrane, and beam element for suspension lines and reinforcement belts.CalculiX automatically expands those 2D elements into 3D during computation, and therefore the thickness for each component should be predefined.Moreover, material properties such as density and elastic modulus should also be specified for each component of the parachute.
There are two types of displacement boundary conditions specified for the parachute system.One is the fixed constraint imposed on the six degrees of freedom of the node where the suspension lines intersect.The fixed constraint equations are expressed as, (5a) The other is the tied constraint imposed on the six degrees of freedom of the nodes between the canopy membrane, suspension lines and reinforcement belts.The tied constraint equations are expressed as, where the superscripts M and S represent the master and slave nodes, respectively.A typical uninflated round parachute model is shown in Fig. 8 to further explain the procedure developed in Fig. 7.The canopy membrane is discretized with the quadratic shell element, and the suspension lines and reinforcement belts are discretized with the quadratic beam element.Each radial line of the canopy is discretized by a set of uniformly-spaced control points {M 0 , M 1 , M 2 , M 3 , M 4 } , and these points are treated as master nodes.The intersection node of the suspension lines is imposed with the fixed constraint expressed by Eq. ( 5).The other end of each suspension line is tied with the canopy membrane by imposing the tied constraint expressed by Eq. ( 6) between the master node M 0 and the slave node SL 0 .Reinforcement belts usu- ally coincide with radial lines of the canopy to maintain the inflated shape.Each reinforcement belt is discretized by a set of control points with the same spacing as that for the radial line, i.e., {SR 0 , SR 1 , SR 2 , SR 3 , SR 4 } in the figure, and the reinforcement belt is tied with the canopy membrane by imposing the tied constraint expressed by Eq. ( 6) between the node pairs {M i , SR i } with i = 0, 1, 2, 3, 4 .By specifying the above constraints, the canopy membrane, suspension lines and reinforcement belts can be assembled as a simple round parachute system as shown in Fig. 8.

Validation
In this section, the benchmark case of a flapping flag in a uniform flow [12,37,42] is considered to validate the present FSI approach.A schematic of the flapping flag in a uniform flow is shown in Fig. 9.The initial shape of the square flag is a flat plate with length L, inclined of α = 0.1π with respect to the inflow direction.The leading edge of the flag is pinned, and the other edges are free to move.
The Reynolds number is Re = U ∞ L/ν = 200 , where U ∞ is the uniform flow veloc- ity, and ν the kinematic viscosity.The size of the computational domain is 8L × 8L × 4L in the x, y and z directions, respectively.The computational domain is discretized with a non-uniform mesh, where local refinement is introduced around the flag and in the wake.The minimum mesh spacing is x = y = z = 0.02L around the flag.The non- dimensional time step used is �tU ∞ /L = 2 × 10 −3 .
The flag is discretized with 1250 quadratic triangular shell elements.The mass ratio between the flag and the surrounding fluid is ρ s h/ρ f L = 1.0 , where h is the flag thick- ness.The non-dimensional bending rigidity is , where the Poisson's ratio is ν s = 0.4.
Figure 10 plots the time history results of the flapping flag in a uniform flow and compares them with those obtained by Tian et al. [37].These results include the transverse displacement of the middle point at the trailing edge of the flag (point A in Fig. 9), and the lift and drag coefficients which are normalized by 1  2 ρ f U 2 ∞ L 2 .Good agreement can be observed from these figures.Moreover, the peak-to-peak excursion amplitude A/L, and the Strouhal number St = fL/U ∞ measured at point A with f the flapping frequency, are compared with the reference data [37,42,43] in Table 1.Good agreement can also be observed from this table.
Furthermore, the deformation of the flag at different time instants within one flapping period is illustrated in Fig. 11, and the instantaneous vortical structure of the flapping flag is shown in Fig. 12, which are qualitatively consistent with those reported in the literature [37,42].

Results and discussion
In this section, the inflation process of several different parachutes is simulated by the present approach to demonstrate its capability.For each parachute system, the thickness of the canopy is h = 5 × 10 −3 m , and the cross section of each suspension line or rein- forcement belt is approximated as a circle with a radius of r = 1 × 10 −3 m .Note that for each case in this study, the uniform velocity is downscaled to U ∞ = 1 m/s to match the Reynolds number Re = 1 × 10 3 , and therefore the elastic modulus is scaled accordingly such that the non-dimensional results can be compared with the reference data.Moreover, each case is simulated for a sufficiently long time to make sure that the results are  converged, and the stable values of the force coefficients and displacements are obtained by averaging over the last five seconds with small variations on the time trace curve.

Inflation of a square parachute
In this subsection, the inflation process of a square parachute in a uniform flow is simulated.Figure 13a shows the inflation model of the square parachute.The length of the canopy is L = 0.8128 m .The canopy membrane is discretized with 1250 triangular shell  2.
The drag coefficient of the parachute is defined as where F D is the total drag force on the canopy surface, and A is the nominal area of the canopy.Time history of the displacement on the canopy center (highlighted red in Fig. 13a) and the drag coefficient is plotted in Fig. 14.The displacement is measured in the x direction.From the figure, we can observe several breathing-in-and-out cycles during the inflation process.The time phase between the drag coefficient at the trough and the center displacement at the crest agrees very well with each other, whereas there is a slight phase difference between the drag coefficient at the crest and the center displacement at the trough.The maximum and minimum displacements are respectively about 0.4 m and 0.12 m.As there is no benchmark data to verify the fidelity of the solution for parachute studies, only a range of reference values is provided.For a square parachute, the drag coefficient typically falls within the range of 1.0 to 2.5 [44,45].The stable value (7)   of C D is approximately 1.26, which is close to the values reported in [12,13] with the same geometric properties.The distinct deformation of the square parachute at several different time instants is shown in Fig. 15.The breathing phenomenon of the parachute can be clearly observed due to the complex interaction between the canopy and the airflow.

Inflation of a cross parachute
In this subsection, the inflation process of a cross parachute is simulated.Figure 13b shows the inflation model of the cross parachute.The cross canopy is constructed by The cross canopy is discretized with 1280 triangular shell elements.There are twenty suspension lines tied with the canopy by the tied constraint, whereas the other end of each suspension line is imposed with the fixed constraint.For each suspension line, the initial length projected in the x direction is the same as that for the square parachute.The material properties of the cross canopy and suspension lines are the same as those for the square parachute listed in Table 2.
Time history of the displacement of the center point on the canopy (highlighted red in Fig. 13b) and the drag coefficient is plotted in Fig. 16.The breathing phenomenon and time phase change for the cross parachute are similar to those for the square parachute.As time evolves, the breathing effect tends to decay, and the amplitude of dX gradually decreases to a stable value of around 0.32 m.A cross parachute typically has a drag coefficient in the range of 0.6 − 0.85 [46,47].The stable value of C D is about 0.79, which is close to the values reported in [12] (approximately 0.7) and [13] (approximately 0.75) with the same geometric properties.
To further observe the breathing phenomenon, the distinct deformation of the cross parachute at several different time instants is shown in Fig. 17. Figure 18 shows the vortex iso-surface of flow around the cross canopy at t = 20 s , where the vortex shedding can be observed behind the canopy.

Inflation of a double cross parachute
To demonstrate the adaptability of the present approach for parachute clusters, the inflation process of a double cross parachute is simulated in this subsection.
Figure 13c shows the deployment of the double cross parachute.The layout of the double cross parachute is obtained by rotating the single cross parachute about the z axis with an angle of 2π/15 and −2π/15 , respectively.Note that contact between the canopies is not considered in this case.The material properties and boundary conditions imposed on the canopy and suspension lines are the same as the single parachute.Figure 20 shows the vortex iso-surface of flow around the cross canopies at t = 20 s .Compared with the single cross parachute, the flow structure is more complex for this  The time history of the center point displacement for each parachute is plotted in Fig. 21a, where the index numbers 1 and 2 are denoted as in Fig. 13c, and the center point for each parachute is denoted as A and B in Fig. 13c, respectively.The displacement of point A in the x direction is almost the same as that for point B, while the displacement of point A in the y direction is almost the opposite of that for point B. The maximum and stable values of the center point displacement in the x direction for each parachute are about 0.25 m and 0.12 m, respectively.The absolute maximum and stable values of the center point displacement in the y direction for each parachute are about 0.14 m and 0.098 m, respectively.Although the magnitude of the y component displacement is generally smaller than that for the x component displacement, the magnitude of stable values for the x and y component displacements are close.Note that the maximum value of the x component displacement occurs at around 2.8 s, while the absolute maximum value of the y component displacement occurs at about 21.2 s.The time history In addition, the x component displacement and drag coefficient for the single cross parachute in Section 6.2 are added in Fig. 21 for comparison.Since the initial deployment between the single and double cross parachute is different, phase differences can be clearly observed between them.At the beginning phase, the variation of dX for the single cross parachute is higher than that for the double case, and this is because the breathing phenomenon is more pronounced for the single cross parachute.As time evolves, the variation for both the single and double cases gradually decreases, and the stable values of dx between them are close to each other.A similar phenomenon can be observed for the drag coefficient curves of the single and double cases.Compared with the single-cross-parachute deployment, the double-cross-parachute deployment can provide larger drag forces.

Inflation of a round parachute
In this subsection, the inflation process of a round parachute in a uniform flow is simulated.The geometric model and boundary conditions for the round parachute are elaborated in Fig. 8.There are 28 suspension lines and reinforcement belts tied with the canopy.A vent hole is designed on the apex of the canopy to allow for airflow ventilation and pressure equilibrium.The radii of the vent hole and canopy are r = 0.04 m and  3.
Time history of the displacement of the center point on the canopy and the drag coefficient is plotted in Fig. 22.Compared to the square parachute or the cross parachute with the same characteristic size, we can observe more breathing-in-and-out cycles during the inflation process.The time phases between the center point displacement and the drag coefficient agree well with each other.The maximum and stable displacements are about 0.24 m and 0.15 m, respectively.The maximum and stable values of C D are around 1.45 and 0.91, respectively.
The distinct deformation of the round parachute at several different time instants is shown in Fig. 23.Except for the typical breathing phenomenon of the parachute, we can also observe that the inflated shape of the canopy can be well preserved by the reinforcement belts.The fundamental frequencies for the square parachute, the cross parachute, and the double cross parachute are all around 0.15 Hz.The fundamental frequency for the round parachute is 0.3 Hz, which is approximately twice as that for the other parachutes.In addition, Fig. 25c compares the Von Mises stress at the fixed end of the cables for each parachute.The magnitude of the stress for the square parachute is much larger than that for the other parachutes.This is because the deformation of the square parachute is the largest among all parachutes, and therefore a larger stress is required at the fixed end of cables to balance the load equilibrium.For each parachute configuration, the phase change of the stress is consistent with that for the x component displacement of the center point.

Conclusions
The IB method is suitable for modeling moving-boundary flows involving complex geometries.Based on the idea of partitioned FSI, this work extends the sharp-interface IB method proposed by Mittal et al. [25] to simulate parachute inflation by utilizing several open-source tools.
The open-source finite-element code CalculiX is employed for its capability in modeling the nonlinear large deformation.The canopy and suspension lines are discretized with quadratic shell and beam elements, respectively.In addition, an automatic and versatile procedure is developed for preparing ingredients of the parachute required by Cal-culiX.The procedure is implemented by Gmsh and Python scripts, and the ingredients include the mesh, boundary conditions, and material properties for each component of the parachute.The open-source multi-physics coupling library preCICE is adopted to handle interfacing between the IB solver and CalculiX.A and accurate data mapping approach is developed to associate the fluid mesh with the solid mesh at the interface.To overcome the strong added mass effect induced by the canopy fabric in a uniform flow, the implicit-coupling scheme with the Aitken acceleration technique is considered to improve the numerical accuracy and stability.Intercommunication such as parallel data transfer and synchronization between the fluid and structure interfaces is also efficiently handled by preCICE.Therefore, the IB solver and CalculiX can be coupled in a minimally-invasive manner by using the preCICE library.The developed approach is validated by a benchmark FSI case of a flapping flag in a uniform flow.
The inflation process of several typical parachutes such as the square parachute, the cross parachute, the double cross parachute and the round parachute is simulated by the present approach.The breathing process, flow structure, canopy displacement and drag coefficient for each case are analyzed to demonstrate that the present approach is applicable for simulating parachute inflation.The research work presented in this paper enriches the numerical methods for simulation of parachute inflation.
Note that the Reynolds number for a real parachute is typically around 10 5 to 10 6 [11], whereas the Reynolds number in this study is reduced to 10 3 and the turbulence effect is not considered.Future work will focus on extending the present FSI approach to turbulent flow regimes by the large eddy simulation and wall modeling technique [48,49].

Fig. 3
Fig. 3 Schematic coupling of the IB solver with CalculiX using preCICE

Fig. 6
Fig. 6 Computational mesh of the fluid domain

Fig. 7
Fig. 7 Automatic and versatile procedure for preparing ingredients of the parachute required by CalculiX

Fig. 9
Fig. 9 Schematic of a flapping flag in a uniform flow

Fig. 10
Fig. 10 Comparison of the time history results of the flapping flag in a uniform flow.a The transverse displacement of point A at the after edge of the flag, b the lift coefficient, c the drag coefficient

Fig. 11 Fig. 12
Fig. 11 Deformation of the flag in one flapping period

Fig. 13
Fig. 13 Geometric model and mesh for different parachutes

Fig. 14
Fig. 14 Time history results of the square parachute

Fig. 16
Fig. 16 Time history results of the cross parachute

Fig. 17
Fig. 17 Deformation of the cross parachute

Fig. 19
Fig. 19 Deformation of the double cross parachute

Fig. 21
Fig. 21 Time history results for inflation of the double cross parachute.a The center point displacements, b the force coefficients

Figure 24 illustrates
Figure 24 illustrates the vortex iso-surface of flow around the round canopy at t = 20 s , where vortex shedding can be well captured behind the canopy.Finally, comparison of several quantities of interest among the square parachute, the cross parachute, the double cross parachute and the round parachute is shown in Fig. 25. Figure 25a compares the time history of the x component displacement of the center point for each parachute.The magnitudes of dX for the cross parachute, the double cross parachute, and the round parachute are close to each other, and the maximum and stable values of dX are around 0.24 m and 0.15 m, respectively.The magnitude of dX for the square parachute is almost twice as that for the other parachutes, where the maximum and stable values of dX

Fig. 22 Fig. 23
Fig. 22 Time history results of the round parachute

Fig. 25
Fig. 25 Comparison of time history results among all parachute configurations.a The x component displacement of the center point, b the fundamental frequency, c the Von Mises stress at the fixed end of the cables

Table 1
Comparison of the excursion amplitude A/L and the Strouhal number St for the flapping flag

Table 2
Material properties of the canopy and suspension lines for the square parachute

Table 3
Material properties of the canopy, suspension lines and reinforcement belts for the round parachute