1 Introduction
Cables as subsystems appear in a wide range of engineering and scientific applications, such as suspension lines in the parachute dynamics fan2014simulation; tezduyar2010space; huang2018simulation; rabinovitch2019towards; sengupta2009supersonic; gao2016numerical; kim20093, offshore drilling and production risers jaiman2009fully; holmes2006simulation; herfjord1999assessment. and airborne refueling systems lofthouse2017cfd; ro2010modeling; zhu2007modeling; styuart2011numerical. When immersed in the flow, cable structures can be responsible for strong fluidstructure interactions that may significantly affect the performance of the system they are attached to. Their interactions with the fluid flow are significant, which cause drag reduction for the supersonic parachute system due to disturbing the front bow shock, and give rise to the vortexinduced vibration even fatigue damage of offshore systems. Therefore, to model it is important.
In computational fluid dynamics (CFD), the local fluid cable subsystem interactions consist of two parts, flowinduced forces on the cable and effect of the structural dynamic response of the cable on the nearby flow. Several efforts have been made to account it. Strip theory herfjord1999assessment; holmes2006simulation, in which computations are conducted on several cross sections, was applied to analyze the vortexinduced vibrations on deepwater risers. Oneway coupling is applied in tezduyar2010space, which constructed empirical flowinduced load formulas on the parachute suspension lines depending on the normal velocity and tabulated drag coefficients. The effect of suspension lines on the fluid flow was modeled by adding source terms based on the inertial and elastic forces of the structure in kim20093; kim20062. Besides, the ”brute force” approach jaiman2009fully was able to resolve complex phenomena due to the presence of a cable, however, the computational cost might be challenging, particularly when the size of its cross section is small compared to the characteristic length of the cable subsystem. Moreover, this issue is exacerbated by the fact that from a structural dynamics viewpoint, the geometry of a cable is typically modeled using line elements, which complicates the task of transferring information between its structural and fluid representations. For bodyfitted CFD meshes, a “dressing” approach based on phantom surface elements and massless rigid elements was proposed in huang2018simulation to address this issue. While it can also be used in non bodyfitted CFD frameworks, this approach has computational disadvantages that are identified and discussed in this work. For this reason, an alternative embedded boundary approach is proposed for computing cabledriven fluid structure interactions. This approach is more robust than the “dressing” approach and is characterized by a superior computational performance. It is based on: master/slave kinematics between the line representing the geometry of a cable typically found in finite element structural models, and a discretization of its true surface that is embedded in the CFD mesh; a highly accurate algorithm for computing the embedded surface displacement based on the beam displacement; and an energyconserving method for transferring distributed forces and moments acting on the nodes of the discrete surface to the nodes of the discrete line enclosed by that surface.
Another challenge for modeling the fluid cable subsystem interactions is the large deformation, caused by the large length diameter ratio or highly flexibility of these cables. To relieve this issue, embedded or immersed boundary methods (EBMs or IBMs) peskin1977numerical; mittal2005immersed; choi2007immersed; tseng2003ghost; wang2011algorithms are used in the present study, which are effective for highly nonlinear FluidStructure Interaction (FSI) problems. Both “dressing” approach and embedded boundary approach are incorporated into the inhouse Eulerian computational framework FIVER wang2011algorithms; lakshminarayan2014embedded; main2017enhanced; huang2018family— A Finite Volume method with Exact twomaterial Riemann Problems. To track the boundary layers around the cable subsystem, adaptive mesh refinement borker2019mesh is applied, which enables reasonable mesh resolution near the cable subsystem.
The remainder of the paper is organized as follows. First, the ”dressing” approach is introduced in section 2. Then, the embedded surface approach and its properties are discussed in section 3. Next, the embedded framework and governing equations are briefly overviewed in section 4. The performance of the discussed models are assessed in section 5 by comparing their results on a model airborne refueling system and a challenging supersonic parachute inflation problem. Finally, conclusions are offered in section 6.
2 Background: the “dressing” approach
The cable subsystems are generally modeled as onedimensional beam elements or cable elements, thanks to the large length diameter ratio, but the fluid domain is generally three dimensional. The detection of these onedimensional elements and defining proper structure normals are ambiguous in this scenario. Hence, it is challenging to capture the local fluidstructure interaction near the cable subsystem. Inspired by the fishbone aeroelastic models, the “dressing” approach considers a ”super” beam element (See fig. 1left), which consists of the onedimensional beam elements, massless rigid beam elements, and a massless discrete surface . The discrete surface represents the real geometry of the cable, where the fluidstructure interaction happens. Namely, the flowinduced forces are evaluated on this surface, and the effect of the cable on the flow is also enforced here through the transmission boundary conditions. The beam elements and the discrete surface are connected with massless rigid beam elements, which are all coupled and solved together by the structural solver. Besides modeling cable subsystems, thoughts of the equivalent beam element are widely used in simplified models, like veryflexible aircrafts patil2000nonlinear; palacios2010structural and turbine blades malcolm2003modeling; wang2014nonlinear.
The ”dressing” approach is able to handle the fluidcable subsystem interaction consistently, and the mesh generating procedure, although tedious, can be automated. However, the mass matrix entries corresponding to the nodes on the dressing surface are zero, due to the massless assumption. This singularity could be an issue for explicit structural integrators, which are needed for highly nonlinear fluidstructure interaction problems, like contact problems or damage mechanics.
3 Alternative embedded boundary approach
The zero mass singularity in the ”dressing” approach can be resolved through a matching procedure farhat1998load, which is common for handling nonconforming fluid/structure interfaces in aeroelastic computations. In the current scenario, the structure solver solves only the dynamics of the onedimensional beam elements and the fluid solver sees only the discrete surface , which allows evaluating the force load and moment precisely. To bridge the dynamics of the onedimensional beam elements and the discrete surface, a masterslave kinematics is applied. The slave surface follows the master beam, and meanwhile the load computed on the salve surface are transferred to the master beam. The detailed load and motion transfer algorithm is as follows:

Match each slave node on the embedded surface with a master point on some beam element using the closest point projection. The superscript of
indicates that this matching procedure might be a manytoone map. The distance vector between these two locations in the initial configuration is denoted by
, which is written as(1) 
At each time step, the displacement and velocity of the slave node are formulated as
(2) where and are displacement and rotation freedoms of the matched beam point, and are are its velocity and angular velocity, respectively. is the rotation matrix at the matched beam point, which depends on the rotation freedoms .

Meanwhile, the force and moment on the master point are induced from the nodal force of the slave node as follows,
(3) where is the number of slave nodes that are matched to .
The masterslave procedure effectively decouples the freedoms of the slave nodes from the structural solver, hence removes the zero mass singularities from the ”dressing” approach and meanwhile accelerates the structural solver. It is worth mentioning that, although decoupled, the aforementioned procedure guarantees conservativity. Consider a virtual displacement field of the flow on the fluidstructure boundary , the virtual work of the fluid traction acting on can be written as
(4) 
where and
are the pressure and stress tensor of the flow,
is the outward normal of , and the conservative load transfer algorithms farhat1998load is applied to compute the nodal force at each slave node. And is the number of master points. Substituting eqs. 3 and 2 into eq. 4 leads to(5) 
Finally, noting that eq. 5 is the virtual work of the onedimensional beam structure, therefore, the energy is conservative in the decoupled masterslave procedure. Moreover, the masterslave kinematics does not require that the matched master point coincides with the beam node, in such case the conservative load transferred algorithm in farhat1998load needs to be applied. Although we focus on cable subsystems with circular crosssections in the present work, this procedure can also be generalized for more complicated surfaces, like bridges kvamsdal1999two, turbine blades and flexible aircrafts.
4 Embedded computational framework
The aforementioned cablesubsystem modeling approaches are incorporated into Embedded or Immersed boundary methods (EBMs or IBMs) peskin1972flow; johansen1998cartesian; fadlun2000combined; kim2001immersed; lohner2004adaptive; tseng2003ghost; mittal2005immersed; choi2007immersed; wang2011algorithms; lakshminarayan2014embedded; uddin2014cartesian, which are also known under other names, including fictitious domain methods glowinski1995numerical and Cartesian methods aftosmis1998robust; johansen1998cartesian. Since they effectively handle FSI applications featured with large structural displacements, rotations, deformations, and/or topological changes farhat2010robust; wang2011algorithms; kramer2013fluid, and cable subsystems characterized by large length diameter ratio, generally undergo large deformations. Though the arbitrary EulerianLagrangian (ALE) method hirt1974arbitrary; donea1982arbitrary, incorporated with mesh motion algorithms or remeshing techniques kvamsdal1999two; herfjord1999assessment; jaiman2009fully, has also been applied to explore cable system dynamics. In the present work, we focus mainly on the Eulerian framework, but the cablesubsystem modeling approaches discussed above can straightforwardly be incorporated into the ALE framework, for problems with moderate displacements.
4.1 Governing equation
Let denote the fixed fluid domain in the Eulerian computational framework, the conservative form of the NavierStokes equations can be written with Einstein notation as
(6) 
here is the density, is the velocity vector, is the pressure, and is the total energy per unit volume, where is the specific internal energy. The viscous stress tensor and heat flux vector are defined as
where is the dynamic viscosity, is the Kronecker delta, is the thermal conductivity, and is the temperature of the gas. The system of equations 6 is closed by assuming that the gas is ideal and calorically perfect,
where and are the gas constant and the specific heat ratio.
The governing equations of the dynamic equilibrium of the structure, including the cablesubsystem, are written in a Lagrangian formulation with Einstein notation
(7) 
where denotes the initial configuration with the material coordinate . denotes the structural material density, denotes its displacement vector, denotes the deformation gradient tensor, denotes the second PiolaKirchhoff stress tensor, and is the vector of body forces acting on . Given a structural material of interest, the closure of eq. 7 is performed by specifying a constitutive law that relates the second PiolaKirchhoff stress tensor to the Green symmetric strain tensor
Dirichlet and/or Neumann boundary conditions are applied to the Dirichlet and Neumann part of the boundary of , as required by the problem of interest.
In addition to eq. 6, eq. 7 and their associated boundary conditions, the FSI problem resulting from the embedding of the structural system in is governed by the transmission conditions
(8) 
and
(9) 
where is the outward unit normal to the deformed configuration of the material interface , and . Since the fluid is assumed to be viscous, additional boundary conditions are specified on : these are the adiabatic or isothermal boundary conditions, and the appropriate boundary conditions for the turbulence model equations when these are presented.
4.2 Imposition of the transmission condition
Imposition of transmission conditions eqs. 9 and 8 in the Eulerian computational framework is the key factor, which distinguishes one embedded boundary method from another. This paper focus on the inhouse Eulerian computational framework— The Finite Volume method with Exact twomaterial Riemann Problems (FIVER) wang2011algorithms; farhat2008higher; main2017enhanced; huang2018family. This framework has been extensively verified and validated for complex multimaterial problems, including FSI problems with large deformations farhat2014ale, turbulent flows lakshminarayan2014embedded, or dynamic fractures wang2015computational.
4.2.1 Noslip boundary condition
The noslip boundary condition eq. 8 embodies the structure to the fluid coupling, i.e. the effect of the structure on the fluid. Its implementation in FIVER consists of two parts: the treatments for the inviscid flux and the viscous flux near the fluid structure interface. For the inviscid flux, a local fluidstructure Riemann problem is formulated and solved at the fluid structure interface wang2011algorithms, which guarantees the characteristic theory and captures well the shock wave or rarefaction wave near the interface. And for the viscous flux, the ghost node method fedkiw1999non; lakshminarayan2014embedded is applied.
4.2.2 Load computation
The force equilibrium condition eq. 9 embodies the fluid to the structure coupling, i.e. the effect of the flow on the structure. The flowinduced load on the embedded surface is generally evaluated first on each Gauss quadrature point of surface elements (See fig. 2). And then, these forces are assembled to the nodal forces. However, the Gauss quadrature point is generally in fluid primal cells that contain inactive nodes. Therefore, to obtain the flowinduced load requires extrapolation, which produces oscillations. Inspired by procedures in huang2018family, modifications are proposed to avoid oscillations in these quantities. The shifted Gaussian points are introduced by shifting each Gaussian point in the direction of the outward normal to the wall, to the point (See fig. 2) defined by
(10) 
where is the characteristic mesh size of the fluid primal cell containing the Gauss point . And the shifted Gaussian point is generally in a primal cell, consists of only active nodes.
The pressure force is computed based on the value on the shifted Gaussian point
(11) 
here the term is derived from the Taylor expansion of the pressure field at and the fact , obtained from the wall normal momentum conservation equation. The pressure at the shifted Gaussian point is obtained by interpolation in the fluid primal cell containing it. As for the shear force, the shear stress tensor depends on the velocity gradient, which is written in the local coordinate with corresponding to the wall normal direction and and corresponding to the other two tangential directions, as follows
(12) 
The components corresponding to the derivatives in the tangential directions are computed based on the finite element approximation
(13) 
where is the primal cell that contains the shifted Gauss point, and is the standard piecewise linear test function associated with each fluid node , The wall normal derivative, which is generally dominant, is approximated directly
(14) 
where is the structure velocity at . When it is necessary, the wall normal derivative can be replaced by the wall model results. The evaluation of the velocity gradients also requires only interpolations.
5 Applications
5.1 Airborne refueling system
The first problem considered here is a model problem from airborne refueling systems ro2010modeling; zhu2007modeling; styuart2011numerical, which is designed to compare the aforementioned ”dressing” approach and embedded boundary approach for modeling the cable system. A flexible hose that trails from the tanker aircraft is considered. Its setup is graphically depicted in fig. 3. The computational domain is a cube of size . The flexible hose is of length , and its length diameter ratio is . The top end is pinned, i.e. all six freedoms are fixed, and the other end is free. Therefore, the flexible hose is supposed to undergo large deformations, when interacting with high speed flows, which motivates the usage of Eulerian computational framework, instead of arbitrary Lagrangian Eulerian framework. The detailed parameters of the hose material are listed in table 1. The inflow conditions match the earth atmosphere of km height, which is also listed in table 1.
Parameter  Description  Value 

Length  m  
Mass of unit length  
Outer diameter  m  
Inner diameter  m  
Young’s modulus  MPa  
Poisson’s ratio  
Altitude  km  
Density  
Pressure  Pa  
Temperature  K  
Mach 
The hose is discretized by 100 beam elements, and its circular cross section is represented by a hexagon shape, And the ”dressing” geometry and embedded boundary surface are discretized by triangular elements. The fluid domain is discretized by tetrahedra (initial mesh). To resolve the hose geometry, boundary layer and vortices with the large deformation, the fluid mesh is adaptively refined and coarsened by the newest vertex bisection algorithm stevenson2008completion; borker2019mesh. The distancebased and Hessianbased mesh adaptation criteria are applied for boundary layer refinement and vortex refinement separately. The characteristics mesh size near the hose is about m, roughly of the diameter. The vortices are tracked by the Hessian of the velocity magnitude with a minimum edge length m.
The viscosity of air is modeled using Sutherland’s viscosity law, which is with kg ms and K. The Reynolds number based on the hose diameter is about . Hence, the flow is assumed to have transitioned to the turbulent regime, which is modeled here using the SA turbulence model spalart1992one. Time discretization of the FSI problem is performed using the secondorder, timeaccurate, implicitimplicit fluidstructure staggered solution procedure developed in the work of farhat2010robust, and the fluidstructure coupling timestep is s. These FSI simulations start from the same quasisteady fluid state obtained from a simulation with a fixed rigid structure. The simulation time is . Figure 4 graphically depicts the time evolution of the adapted fluid mesh, the velocity magnitude of the fluid flow, and the hose displacement.
The hose is drifted by the highspeed flow, and meanwhile interacts with the vortices behind. The AMR well tracks the boundary layer and the flow features. And it is worth mentioning that both ”dressing” approach and embedded boundary approach deliver almost the same flow features and structure displacement fields. The time histories of the drag force on the hose due to the fluid load and the xdisplacements of the bottom node and the middle node obtained by both approaches are reported in fig. 5. Great agreements are achieved by both approaches.
5.2 Supersonic parachute inflation dynamics
Next, the inflation dynamics of a DiskGapBand (DGB) parachute in the lowdensity, lowpressure supersonic Mars atmosphere is simulated cruz2014reconstruction; derkevorkian2019effects; rabinovitch2019towards. Its main purpose is to understand the effects of the suspension line subsystem on the parachute performance during the deceleration process.
The DGB parachute system (See fig. 5(a)), which has successfully landed Curiosity Rover on the Mars, contains three main components: the canopy, which is made of F111 nylon, the suspension lines, which are made of Technora T221 braided cords, and the reentry vehicle cruz2014reconstruction. Material properties and geometric parameters of the parachute system are listed in table 2. The simulation starts from the line stretch stage with a folded parachute (See fig. 5(b)), an idealized analytical representation of the line stretch configuration and the corresponding prestress are employed for the structure initial state. The supersonic flow is incoming at , and Pa. Note that the specified freestream density and pressure conditions are similar to those observed in the standard Martian atmosphere.
Component  Properties 

Canopy  Diameter = 
Thickness =  
=  
=  
=  
Porosity =  
Suspension lines  Length = 
Diameter =  
=  
Density = 
Since the Martian atmosphere is mainly composed of carbon dioxide, the viscosity of this gas is modeled using Sutherland’s viscosity law – that is, – with kg ms and K. The Reynolds number based on the canopy diameter is . Hence, the flow is assumed to have transitioned to the turbulent regime, which is modeled here using the Vreman turbulence model vreman2004eddy, with model constant . The embedding computational fluid domain of size is initially discretized by a mesh composed of Kuhn simplices stevenson2008completion; borker2019mesh. Specifically, this initial tetrahedral mesh contains vertices and tetrahedra. Adaptive mesh refinement (AMR) mitchell1988unified; stevenson2008completion; borker2019mesh is applied to track the boundary layer and the flow feature, the characteristic mesh sizes near the reentry vehicle, suspension lines and the canopy are cm, mm, and cm, respectively, and the size in the wake regime or near the shock is cm. The canopy of the DGB parachute consists of band gores and disk gores. These are discretized here by geometrically nonlinear thinshell ANDES elements militello1991first (although the membrane stiffness of the fabric is significantly larger than its bending stiffness, both stiffnesses are considered here). The canopy is made of F111 nylon with void fraction, and the permeability is modeled by the homogenized porous wall model proposed in huang2019homogenized. There are suspension lines, and each of them is discretized by geometrically nonlinear beam elements. And the reentry vehicle is treated as a fixed rigid body, which is unseen from the structure model. Due to the massive contact of the parachute during its inflation, the explicit central difference scheme is applied for the time discretization of the structure solver, therefore only the embedded surface approach in section 3 is capable of modeling the fluid suspension line subsystem interactions. The crosssection geometry of suspension lines is assumed to be circle with radius mm and represented as the hexagon.
First, the parachute configured in its initial position shown in fig. 5(b) is assumed to be rigid and the quasi steadystate solution at is computed (See fig. 7a). Next, the fluid state is initialized using this quasi steadystate flow solution, the deformable structure is initialized using its initial state described above, and the FSI problem associated with the inflation of the DGB parachute is simulated in the time interval s, during which the inflation process is expected to have completed and several breath cycles are captured. Time discretization of the FSI problem is performed using the secondorder, timeaccurate, implicitexplicit fluidstructure staggered solution procedure developed in the work of farhat2010robust and the fluidstructure coupling timestep is roughly s.
Figures 7 graphically depict the time evolution of the parachute and the flow Mach number field. The bow shock in the front of the canopy, disturbed by the wake generated by the reentry vehicle, vibrates along with the breath cycles of the canopy. Jetlike flow, ejected through the canopy vent and gaps between the band/disk gores, interact with the turbulent wakes behind the canopy. The parachute is fully inflated at s, and it starts breathing after that time.
Although the suspension lines are very thin, they increase geometric blockage in the front of the canopy region, which slightly reduces the flow speed. Moreover, in the interactions with flows, these suspension lines can generate shocks (See fig. 7e) that alter the wake behind the reentry vehicle, and even disrupt the bowshock, which agrees with the observations in sengupta2009supersonic from subscale parachute inflation experiments.
The timehistories of the drag force computed by postprocessing the FSI simulation results is plotted in fig. 8. Moreover, drag performance of a simulation without considering the effect of the suspension lines is also presented. The drag force on the suspension lines is negligible, roughly , however, the shocks from the suspension lines create disturbances ahead of the canopy which in turn disrupt the parachute bowshock, mix the high pressure and low pressure flows, and reduce the drag force. The drag force data cruz2014reconstruction collected from the Curiosity Rover during its Mars landing is also presented in fig. 8. The peak drag prediction, which is related to the material failure, is in good agreement with the landing data, and the relative error is less than .
6 Conclusions
In this paper, we have proposed a simple and practical embedded surface approach to model the fluid cable subsystem interactions in engineering systems. This approach is based on master/slave kinematics with the master beam elements and a slave embedded surface, which represents its real geometry. The fluid solver sees only the slave embedded surface which is used to evaluate the load, while the structure solver sees only the master beam elements. Moreover, it is suitable for Eulerian computational framework. Hence, combining with adaptive mesh refinement techniques, the approach effectively handles large deformations of the cable subsystem. Its application to the simulation of a model airborne refueling system and a challenging supersonic parachute inflation problem has clearly demonstrated its practicality and reasonable computational efficiency. An interesting area for future work would be to extend current model for underresolved slender structures, like turbine blades and fuselage.
Acknowledgments
Daniel Z. Huang, Philip Avery, and Charbel Farhat acknowledge partial support by the Jet Propulsion Laboratory (JPL) under Contract JPLRSA No. 1590208, and partial support by the National Aeronautics and Space Administration (NASA) under Early Stage Innovations (ESI) Grant NASANNX17AD02G. Daniel Z. Huang thanks Kvamsdal Trond for helpful discussions.