An Embedded Boundary Approach for Resolving the Contribution of Cable Subsystems to Fully Coupled Fluid-Structure Interaction
AAn Embedded Boundary Approach for Resolving the Contribution of CableSubsystems to Fully Coupled Fluid-Structure Interaction
Daniel Z. Huang a , Philip Avery c , Charbel Farhat a,b,c a Institute for Computational and Mathematical Engineering, Stanford University, Stanford, CA, 94305 b Mechanical Engineering, Stanford University, Stanford, CA, 94305 c Aeronautics and Astronautics, Stanford University, Stanford, CA, 94305
Abstract
Cable subsystems characterized by long, slender, and flexible structural elements are featured in numerousengineering systems. In each of them, interaction between an individual cable and the surrounding fluid isinevitable. Such a Fluid-Structure Interaction (FSI) has received little attention in the literature, possiblydue to the inherent complexity associated with fluid and structural semi-discretizations of disparate spatialdimensions. This paper proposes an embedded boundary approach for filling this gap, where the dynamicsof the cable are captured by a standard finite element representation C of its centerline, while its geometryis represented by a discrete surface Σ h that is embedded in the fluid mesh. The proposed approach is builton master-slave kinematics between C and Σ h , a simple algorithm for computing the motion/deformationof Σ h based on the dynamic state of C , and an energy-conserving method for transferring to C the loadscomputed on Σ h . Its effectiveness is demonstrated for two highly nonlinear applications featuring largedeformations and/or motions of a cable subsystem and turbulent flows: an aerial refueling model problem,and a challenging supersonic parachute inflation problem. The proposed approach is verified using numericaldata, and validated using real flight data. Keywords: cable dynamics, embedded boundary, fluid-structure interaction, immersed boundary,parachute inflation
1. Introduction
Cable subsystems appear in a wide range of engineering systems and scientific applications, such assuspension lines of parachutes and other atmospheric decelerators [1, 2, 3, 4, 5, 6], offshore drilling andproduction risers [7, 8, 9], and booms of airborne refueling systems [10, 11, 12, 13]. When immersed in aflow, cable structures can be responsible for strong Fluid-Structure Interactions (FSIs) that may significantlyaffect the performance of the system to which they are attached. For example, in the case of a dynamic,supersonic parachute inflation, FSIs involving the suspension lines give rise to a significant drag reductiondue to the disturbance of the front bow shock, while FSIs associated with marine drilling risers can lead tovortex-induced vibrations and fatigue damage of offshore systems. Consequently, the ability to accuratelyand affordably model such phenomena has the potential to become a valuable tool for the design andmaintenance of the aforementioned and related engineering systems.In the context of an individual cable, capturing a local, two-way, FSI interaction requires in general theaccurate prediction of: the flow-induced generalized forces and moments on the FE representation of thecable; and the effects of the structural dynamic response of the cable on the nearby fluid flow. Several effortshave been made to account for one or both of these two requirements. Strip theory, where computations areconducted on several cross sections, has been applied to the analysis of vortex-induced vibrations of deepwater
Email addresses: [email protected] (Daniel Z. Huang), [email protected] (Philip Avery), [email protected] (Charbel Farhat) a r X i v : . [ c s . C E ] N ov isers [8, 9]. Empirical flow-induced load formulae and tabulated drag coefficients were adopted in [2] forcomputing one-way FSIs associated with parachute suspension lines. In [6, 14], the effect of suspension lineson the fluid flow was modeled using source terms based on the inertial and elastic forces generated by thestructure. By contrast, a brute force approach based on a body-fitted mesh that resolves the boundary layerof a dynamic cable was used in [7] to successfully capture complex FSI phenomena driven by the cable.However: the computational cost of this nonadaptive approach can be overwhelming, particularly when thesize of the cross section of the cable is small compared to its characteristic length; and the robustness of itsunderlying mesh motion algorithm needed for maintaining at all times the Computational Fluid Dynamics(CFD) mesh body-fitted is usually limited to small-amplitude dynamics.In addition to the computational complexity issue mentioned above for the case of the brute force approachadopted in [7], the accurate prediction of a local, two-way coupled, FSI associated with a cable subsystempresents many other challenges. First and foremost, the task of transferring information between fluid andstructural meshes in order to discretize the governing fluid-structure transmission conditions is complicatedby the fact that in structural dynamics, the Finite Element (FE) representation of a cable is typicallycarried out using line elements, and therefore is topologically One-Dimensional (1D). To address this issue, a“dressing” approach based on “phantom” surface elements and massless rigid elements was developed by thethird author almost two decades ago [15]. This noteworthy approach is equally applicable in both contextsof body-fitted [16] and non body-fitted CFD meshes [3]: however, it has computational disadvantages thatare identified and discussed in this paper. For this reason, an alternative approach is proposed here forcomputing cable-driven FSI. This new approach is appropriate for the solution of both one-way and two-way FSI problems, is far more robust and computationally efficient than the brute force approach adoptedin [7], and is more comprehensive and user-friendly than the aforementioned dressing approach. It is basedon: master-slave kinematics between the discrete line representing the centerline of a cable typically foundin FE structural models, and a discretization of its true surface that is embedded (or immersed) in thecomputational fluid domain; an accurate algorithm for computing the displacement of the discrete surfaceof the cable based on the displacement and rotation of its centerline; and an energy-conserving method fordistributing the flow-induced forces and moments on the nodes of the FE representation of the cable, basedon the Degrees Of Freedom (DOFs) of this representation.Another challenge for cable-driven FSI is the treatment of the potentially large motions and/or defor-mations of the cable, due its potentially large length-to-diameter ratio and/or high intrinsic flexibility. Inthe general case, this challenge is best addressed by adopting the Eulerian computational framework forFSI which calls for an Embedded Boundary Method (EBM, e.g. [17, 18]) – also known as an ImmersedBoundary Method (IBM, e.g. [19, 20]), or a Ghost Fluid Method (GFM, e.g. [21, 22]) – for CFD and FSI.Here, the Finite Volume method with Exact two-material Riemann problems (FIVER) – which itself is anembedded boundary framework for multi-material computations rather than a mere EBM for CFD andFSI [18, 23, 24, 25, 26] – is the chosen EBM. Specifically, both the dressing approach and the alternativemaster-slave kinematic approach proposed in this paper are incorporated into the Eulerian computationalframework FIVER, and the Adaptive Mesh Refinement (AMR) method for EBMs described in [27] is appliedto track the boundary layers around the cable subsystem and keep them at all times well resolved.The remainder of this paper is organized as follows. First, the dressing approach for cable-driven FSIproblems is overviewed and discussed in Section 2. Then, the alternative master-slave kinematic approachis introduced and analyzed in Section 3. Next, the EBM for FSI in which both aforementioned approachesare incorporated is overviewed in Section 4, in order to keep this paper as self-contained as possible. Inparticular, small but effective novel contributions are made to this EBM and its AMR framework. Theperformance of the proposed approach for modeling cable-driven FSI is contrasted with that of the dressingapproach in Section 5, first for an idealized airborne refueling system and then for a challenging, dynamic,supersonic parachute inflation problem. Finally, conclusions are offered in Section 6.
2. Dressing approach for FSI models with disparate spatial dimensions
As already mentioned above, the structural dynamics of cable subsystems are generally modeled usingtopologically 1D beam elements, or special-purpose variants usually referred to as cable elements. This is2ue to their large length-to-diameter ratio. Semi-discretizations of cables based on beam elements havevastly superior computational efficiency than counterparts based on Three-Dimensional (3D) solid elements,provided that certain simplifying assumptions are appropriate – for example, the assumption that planarsections initially normal to the longitudinal axis of the cable remain planar and normal to this axis. On theother hand, a computational fluid domain is usually semi-discretized using 3D elements such as tetrahedra,hexahedra, and/or triangular prisms. Whether the FSI computations are performed using the ArbitraryLagrangian-Eulerian (ALE) framework and a dynamic, body-fitted CFD mesh [28, 29, 30], or the Eulerianframework, an embedded discrete surface and an EBM for CFD as well as FSI [18, 23, 24, 25, 26], somedifficulties if not challenges arise in the presence of this type of mixed semi-discretizations with disparatespatial dimensions. For instance, consider the case where the FSI computations are performed using the Eu-lerian framework, which is the more appropriate framework for highly nonlinear FSI problems characterizedby large structural motions and/or deformations as well as topological changes. Within this computationalframework, the detection of 1D structural elements embedded in the non body-fitted CFD mesh and the defi-nition as well as construction of associated geometric characteristics such as surface normals may be difficult,ambiguous, and/or computationally intensive. For this reason, capturing a cable-driven FSI using numericaltechniques developed for other types of FSI problems characterized by convenient, Two-Dimensional (2D)fluid/structure interfaces is neither necessarily robust, nor necessarily computationally efficient, or evenpractical.Developed initially for the so-called “fish-bone” aeroelastic computational models, the dressing approachfirst proposed in [15] for facilitating the enforcement of the fluid-structure transmission conditions at theinterface between fluid and structural semi-discretizations with disparate spatial dimensions can be describedas a modeling approach based on the concept of a dynamically equivalent superelement (see Figure 1-left).In this approach, a superelement consists of: • A 1D flexible beam or cable element located at the centerline of the cable. • A cylindrical assembly of: – A set of massless , zero-stiffness , surface elements of typical size h – referred to for this reason asphantom elements – and their nodes, located on the geometrically true surface of the cable, Σ. – A network of massless , rigid beam elements connecting each FE node on the centerline of thecable to the adjacent nodes of phantom elements.Let Σ h and N denote the global set of phantom elements and global network of rigid beam elements,respectively. From a structural dynamics viewpoint, each superelement defined above preserves the dynamicsof the corresponding beam or cable element: in this sense, it is dynamically equivalent to this element. Indeed,this superelement has the same mass and mass distribution as the cable element it is associated with, becauseall elements in Σ h and those of N are massless. It has also the same stiffness and stiffness distribution becauseall elements in Σ h have zero stiffness and those in N are rigid. Unlike a beam or cable element however,this superelement is equipped with a representation of the true geometry of the surface of the cable wherethe FSI of interest occurs, which facilitates the loads evaluation process. These loads are transferred to theFE representation of the (flexible) cable via the rigid (beam) elements in N .Specifically, Σ h can be used in both contexts of body-fitted and non body-fitted CFD meshes to directlycouple the computational fluid and structural domains at the fluid/structure interface using a matching (orpairing) procedure such as that described in [31], which is commonly used for enforcing the semi-discretetransmission conditions across a non-conforming fluid/structure interface. The flow-induced loads can beconveniently evaluated on Σ h and transferred to the FE structural representation of the cable, and thekinematics of the cable (including time-derivatives) can be equally conveniently transferred to the compu-tational fluid domain to update its wall boundary conditions. The massless rigid beam elements couplethe DOFs of the discrete surface Σ h with those of the FE semi-discretization of the cable using algebraicconstraint equations. These can be enforced using, for example, the Lagrange multiplier method, whichtransforms the overall governing semi-discrete equations of equilibrium of the structural subsystem into aglobal, Differential-Algebraic Equation (DAE). 3 igure 1: Schematics of the dressing approach based on dynamically equivalent superelements (left), and counterpart schematicsof the master-slave kinematic approach (right): the discrete surface Σ h representing the true geometry of the surface of thecable encloses the topologically 1D FE representation of the cable (red); its nodes are connected to the discretized cable bymassless rigid elements (green) in the dressing approach (left), and by kinematics constraints (yellow) in the master-slavekinematic approach (right); and { S ji } n i j =1 denotes the set of n i coplanar nodes in Σ h whose kinematics are slaved to those ofthe corresponding master point M i located at the intersection of the discrete cross section defined by the set of nodes { S ji } andthe FE representation of the cable. The dressing approach outlined above can be used to construct computational FSI models for a varietyof applications. These include, for example, the aeroelastic analysis of highly flexible aircraft, rotorcraft, andturbine blades. It is equally applicable to cable-driven FSI problems. Although tedious, the generation of theset of surface elements Σ h and the construction of the network of rigid elements N can be automated [32].Nevertheless, this approach has a few noteworthy shortcomings, particularly for cable applications: • Its practical implementation requires a certain discipline. For example, Σ h and N are best constructedso that each DOF attached to each node of a phantom element in Σ h is connected to a correspondingDOF of a rigid element in N , or otherwise to any other DOF of the FE structural model, in orderto avoid the generation of artificial mechanisms in the overall FE structural model – that is, theintroduction in this model of DOFs with zero associated stiffness. • It increases the number of structural DOFs by typically an order of magnitude. Given the computa-tional cost associated with modeling a turbulent flow problem, this issue may be relatively insignificant.However, it is by no means desirable. • It leads to a mass matrix of the overall FE structural model that is singular. Specifically, this matrixhas a zero entry at each DOF attached to a node of an element in Σ h , due to the massless nature ofsuch an element as well as the massless nature of any rigid beam element in N to which it is connected.Consequently, it unnecessarily complicates the implementation of various computational modules of anFSI solver: – For explicit time-integration schemes requiring a positive definite mass matrix – which are com-monly used for highly nonlinear FSI problems involving contact and/or crack propagation (e.g.,see [24]) – it necessitates the implementation of a cumbersome static condensation in the maintime-integration loop of the structural analyzer. – For implicit time-integration schemes, the constraint equations underlying the rigid elements in N introduce a destabilizing effect in the overall FE structural dynamics subsystem [33], due to its4AE aspect. The mitigation of this effect requires either the introduction in the time-integratorof artificial damping, which lowers its order of time-accuracy, or a special purpose formulation ofthe structural dynamics equations of equilibrium, which complicates its numerical implementation[33].By contrast, the alternative approach presented next for facilitating the enforcement of fluid-structuretransmission conditions at an interface between fluid and structural semi-discretizations with disparate spa-tial dimensions does not suffer from any such pitfalls.
3. Alternative master-slave kinematic approach
The issues associated with the dressing technique outlined above can be resolved through an alternativeapproach for achieving the same objective, where (see Figure 1-right): • A discrete representation Σ h of the surface of each cable is constructed as: a collection of sets ofcoplanar nodes { S ji } n i j =1 , where each set defines a discrete representation of an i -th cross section of thecable; and a set of surface elements that connect these nodes. • Σ h is embedded in the given CFD mesh associated with the overall computational fluid domain. • Three displacement DOFs are attached to each node S ji ∈ Σ h . • These DOFs are slaved to the translational as well as rotational motion of the cable at the point M i located at the intersection of the discrete cross section defined by { S ji } n i j =1 and the FE representation C of the cable ( M i is also the closest point in C to any node S ji ∈ Σ h ). • The translational as well as rotational motion of the cable at each point M i is computed by interpolationof the displacement and rotational DOFs attached to the nodes of the FE beam or cable element e i containing the point M i . Remark 1.
The notation Σ h adopted in the dressing approach for describing the discrete representation ofthe surface of a cable is reused here because the same discretization of the cable’s surface can be reused inthe alternative approach described herein. However, the dynamics of Σ h differ in both approaches. As in the dressing approach, the CFD mesh associated with the overall computational fluid domaincan be in this case body-fitted or non-body fitted. However, given the intrinsic advantages of the Eulerianframework for FSI problems characterized by large structural motions and/or deformations, attention isfocused here on the case of a non body-fitted mesh – and therefore on the case of a preferred EBM for CFDand FSI. Furthermore, for simplicity but without any loss of generality, the structural subsystem is assumedto consist of a single cable: this simplifies the description of the proposed method, but in no way limits itsscope of applications to a single cable or an assembly of cables.In the alternative approach described here, the chosen structural analyzer computes at each time-step thestate of the FE representation C of the cable subsystem. However, it is not given direct access to the preferredflow solver. Similarly, the flow solver is not given direct access to the computed state of the discretized cablesubsystem. Instead, the kinematics of Σ h – that is, the position and velocity of this embedded discretesurface – are slaved to the computed state of the FE representation C of the cable, and the flow-inducedloads are computed on Σ h before they are transferred to C . Specifically, these motion and loads transfersare performed as follows: • Preprocessing step: – Pair each slave node S ji ∈ Σ h , j = 1 , · · · , n i , with the closest element e i ∈ C containing themaster point M i . (Note that the superscript j highlights the surjective aspect of the function S ji −→ M i : specifically, n i ≥ S ji are typically paired with one point M i ). The initial5istance vector d ji between these two locations – that is, the distance vector at t = t – is definedby x S ji = x M i + d ji ⇔ d ji = x S ji − x M i where here and throughout the remainder of this paper, the bold font designates a vector quantityand x P denotes the position vector of a point P . Consistently with the modeling assumptions ofa FE beam element, the distance vector d ji is assumed to be time-independent. • At each time-instance t n : – Compute the displacement u nS ji and velocity ˙ u nS ji of the slave node S ji as follows u nS ji = u nM i + R ( θ nM i ) d ji − d ji and ˙ u nS ji = ˙ u nM i + ω nM i × R ( θ nM i ) d ji (1)where: the dot designates the time-derivative; u nM i and θ nM i denote the interpolated displacementand rotation vectors at the point M i , respectively; ˙ u nM i and ω nM i are the interpolated velocity andangular velocity vectors at the point M i , respectively; and R is the rotation matrix at M i anddepends on θ nM i . – Compute the force vector f nS ji at each slave node S ji as follows f nS ji = (cid:90) Σ nh (cid:16) − p n n n + τ n n n (cid:17) φ S ji d Σ h (2)where: p n and τ n denote the pressure and viscous stress tensor of the flow at time t n ; n n denotesthe outward normal to Σ nh at time t n ; and φ S ji denotes a local shape function associated withthe node S ji ∈ Σ h (for example, the characteristic function at node S ji in the case of a finitevolume semi-discretization of the flow equations). In practice, Equation (2) is evaluated using aquadrature rule in each surface element (see Section 4.2). – Compute the force vector f nM i and moment vector m nM i at the point M i as follows f nM i = n i (cid:88) j =1 f nS ji and m nM i = n i (cid:88) j =1 R ( θ nM i ) d ji × f nS ji (3) – Apply the load transfer method presented in [31] to deduce the following expressions of thegeneralized force and moment vectors acting on a FE node N j of C , f nN j and m nN j , respectively f nN j = (cid:88) i : N j ∈ E i f nM i φ dN j ( M i ) and m nN j = (cid:88) i : N j ∈ E i m nM i φ rN j ( M i ) (4)where φ dN j and φ rN j denote the standard shape functions associated with the displacement androtational DOFs at node N j ∈ C , respectively.The master-slave kinematic approach described above for enforcing the fluid-structure transmission con-ditions at an interface between fluid and structural semi-discretizations with disparate spatial dimensionseffectively decouples the DOFs of the slave nodes from the structural analyzer, thereby avoiding zero massand other destabilizing singularities that may arise in the dressing approach and reducing computationaloverhead. Yet, despite this decoupling, this alternative approach is globally conservative as shown below.Consider at time t n a virtual fluid displacement field δ u nF of the fluid/structure boundary Σ h . Thecorresponding virtual work performed by the fluid tractions acting on Σ h can be written as − δ W nF = (cid:90) Σ nh (cid:16) − p n n n + τ n n n (cid:17) · δ u nF d Σ h = n M (cid:88) i =1 n i (cid:88) j =1 f nS ji · δ u nS ji (5)6here: n M denotes the number of discrete cross sections of the cable represented in Σ h – or equivalently,the number of master points M i ; f S ji is given in Equation (2); and the no-slip displacement transmissioncondition (viscous flow) is assumed to hold ( δ u nF = δ u nS ).Let n N denote the number of FE nodes in C . From Equation (5), Equation (1), Equation (3), andEquation (4), it follows that − δ W nF = n M (cid:88) i =1 n i (cid:88) j =1 f nS ji · (cid:16) δ u nM i + δ θ nM i × R ( θ nM i ) d ji (cid:17) = n M (cid:88) i =1 n i (cid:88) j =1 f nS ji · δ u nM i + (cid:16) R ( θ nM i ) d ji × f nS ji (cid:17) · δ θ nM i = n M (cid:88) i =1 f nM i · δ u nM i + m nM i · δ θ nM i = n M (cid:88) i =1 f nM i · (cid:16) (cid:88) N j ∈ E i φ N j ( M i ) δ u nN j (cid:17) + m nM i · (cid:16) (cid:88) N j ∈ E i φ N j ( M i ) δ θ nN j (cid:17) = n N (cid:88) j =1 f nN j · δ u nN j + m nN j · δ θ nN j = δ W nS which shows that the master-slave kinematic approach proposed here for performing motion and load trans-fers in cable-driven FSI is globally conservative. Remark 2.
The master-slave kinematic approach described above is independent of the shape of the crosssection of the cable, albeit this shape affects the construction of the embedded discrete surface Σ h .
4. Eulerian framework for fluid-structure interaction
Again, both dressing and master-slave kinematic approaches for enabling the simulation of cable-drivenFSI problems characterized by disparate spatial dimensions of the computational fluid and structural domainscan work with body-fitted as well as non body-fitted CFD meshes. However, because long, slender cablestend to be highly flexible and therefore tend to undergo large motions and/or deformations, the focus ofthis paper is set on non body-fitted CFD meshes and therefore on the Eulerian framework for FSI equippedwith an EBM. Since all FSI computations reported in this paper are performed using the second-order EBMFIVER [17, 18, 23, 24, 25, 26], this Finite Volume (FV)-based EBM is overviewed below in order to keepthis paper as self-contained as possible. Specifically, attention is focused on the case where the embeddeddiscrete surface Σ h exclusively emanates from the representation of the surface of a cable in the master-slavekinematic approach described in Section 3. For the case of an arbitrary embedded discrete surface Σ h – andtherefore the case of more general fluid-structure transmission conditions than Equation (1) and Equation (4)– and/or for further details on the FIVER framework, the reader is referred to the aforementioned references. Consider a compressible, turbulent, viscous flow governed by the Navier-Stokes equations with turbulencemodeling. If these equations are written in conservation form, their semi-discretization by a second-order,vertex-based, FV method leads to the following system of ordinary differential equations V ˙ W + F ( W ) − G ( W , W t ) = (6)where W denotes the vector of semi-discrete, conservative, fluid state variables, V is the diagonal matrixstoring the dual cells (or control volumes), F denotes the vector of numerical convective fluxes, G denotesthe vector of numerical diffusive fluxes and source terms due to turbulence modeling, and W t denotes thesemi-discrete state vector of the turbulence model conservation variables.7 .1.1. Numerical convective fluxes Typically, the numerical convective fluxes are computed in a vertex-based FV method on an edge-by-edge basis. Let Ω F denote the embedding computational fluid domain, Ω Fh denote its discretization by astructured or unstructured, non body-fitted, CFD mesh with dual cells C i , and Ω Σ h h denote the restrictionof Ω Fh to the volume enclosed by the embedded discrete surface Σ h representing here the geometrically truesurface of the cable (see Figure 2). Let also V i denote a generic node of the embedding CFD mesh, K ( V i )denote the set of nodes V j of this mesh that are connected by an edge V i V j to V i , and ν ij denote the unitoutward normal to the control volume boundary facet ∂C ij connecting the centroids of the elements of theembedding CFD mesh sharing the nodes V i and V j (see Figure 2). Figure 2: Discretization Ω Fh of an embedding fluid domain, dual cell (control volume) C i , boundary facet ∂C ij , unit outwardnormal ν ij , and embedded discrete surface Σ h (two-dimensional case). At any node V i away from Σ h – that is, any node where ∀ V j ∈ K ( V i ) , V i V j (cid:84) Σ h = {∅} ⇔ V i V j ∈ Ω Fh \ Ω Σ h h – a second-order, vertex-based, FV method computes the vector of numerical fluxes F i as follows F i = (cid:88) V j ∈ K ( V i ) V i V j (cid:84) Σ h = {∅} Φ ij ( W i , W j , ν ij , EOS)where Φ ij denotes a numerical flux vector function associated, for example, with a second-order extensionof a first-order upwind scheme – such as Roe’s approximate Riemann solver [34] – based on the MonotonicUpwind Scheme Conservation Law (MUSCL) [35], W i and W j denote the semi-discrete fluid state vectorsat the nodes V i and V j , respectively, and EOS refers to the Equation of State characterizing the consideredcompressible fluid.In the vicinity of the embedded discrete surface Σ h defined here by V i V j (cid:84) Σ h (cid:54) = {∅} , a “mixed (dual) cell”problem arises. For example, if V i ∈ Ω Fh \ Ω Σ h h but V i V j (cid:84) Σ h (cid:54) = {∅} – or in other words, V i ∈ Ω Fh \ Ω Σ h h but V j ∈ Ω Σ h h – the computational dual cell C i attached to V i is in this case a mixed cell as it is partly occupiedby the fluid and partly occupied by the structure (see Figure 2 and Figure 3): such a situation complicatesthe semi-discretization process in C i , and in particular the evaluation of the numerical convective flux vectorfunction Φ ij across the boundary facet ∂C ij . The second-order EBM FIVER [25, 26] addresses this issue byreconstructing the semi-discrete fluid state vector W i at Σ h and solving there a 1D, exact, fluid-structure,half Riemann problem designed to correctly semi-discretize the convective flux in a mixed cell. At eachtime-step t n , this can be written as described below, with all superscripts referring to the time-instance t n omitted except in t n and t n +1 in order to keep the notation as simple as possible:1. Compute a linear reconstruction of the vector of semi-discrete, primitive, fluid state variables at the8ntersection point I ij of V i V j and Σ h w I ij = w i + ∇ w i · ( x I ij − x i ) (7)where w = ( ρ, v , p ), and ρ , v , and p denote the fluid density, velocity vector, and pressure, respectively.2. Construct and solve in the time-interval [ t n , t n +1 ], at the moving point I ij and along the axis (cid:126)ξ whoseorigin is the moving point I ij and direction is the normal n I ij to the moving material interface Σ h at I ij , the following local, 1D, exact, fluid-structure, half Riemann problem ∂ w ∂t + ∂∂ξ [ F ( w ) − ˙ u Iij w ] = 0 , ξ ∈ [0 , ∞ ) w ( ξ,
0) = w I ij , ξ ∈ [0 , ∞ ) (cid:0) v I ij := v Iij · n I ij (cid:1) ( t − t n ) = (cid:0) ˙ u Iij := ˙ u Iij · n I ij (cid:1) ( t − t n ) ∀ t ∈ [ t n , t n +1 ] (8)where: w ( ξ, t ) = (cid:0) ρ ( ξ, t ) , ( v ( ξ, t ) := v ( ξ, t ) · n I ij ) , p ( ξ, t ) (cid:1) ; ξ denotes the abscissa along the (cid:126)ξ axis de-fined above; w I ij = ( ρ I ij , v I ij , p I ij ) is the reconstructed vector of discrete, primitive, fluid-state variablesat the point I ij (see Equation (7)); F denotes the continuous convective flux vector; and ˙ u denotesthe velocity vector of Σ h and is given by the second of Equations (1). Assuming that ˙ u I ij is constantwithin the time-interval [ t n , t n +1 ], the solution of the above half Riemann problem can be derivedanalytically when the EOS of the fluid allows it, or otherwise efficiently computed numerically usingsparse grid tabulations [18]. In either case, this solution consists of two constant states separated inthe ( ξ, t ) plane by either a shock wave or a rarefaction fan (see Figure 3). In practice, however, onlythe solution at the point I ij , w R = (cid:16) ρ R I ij , ( v R I ij = ˙ u Iij ) , p R I ij (cid:17) , needs be computed and stored. It isworthwhile noting that for an inviscid flow, solving the above half Riemann problem enforces the sliptransmission condition at the point I ij through the third of Equations (8), and computes the valuesof the flow density and pressure at this point. For a viscous flow however, solving the half Riemannproblem (8) simply computes the fluid state vector at the intersection point I ij of the edge V i V j andΣ h in order to enable in Step 5 below the computation of the numerical flux vector function Φ ij acrossthe boundary facet ∂C ij of the mixed cell C i – and hence, the semi-discretization of the convectivefluxes in a mixed cell. In both cases, it enables the numerical approximation to capture any shock orrarefaction wave near the fluid/structure interface Σ h .3. Update the discrete fluid velocity vector at the point I ij by replacing its reconstructed normal compo-nent by v R I ij v (cid:63)I ij = v R I ij n I ij + (cid:0) v I ij − ( v I ij · n I ij ) n I ij (cid:1) and update the vector of semi-discrete, conservative, fluid state variables at the point I ij accordingly W (cid:63)I ij = W ( ρ R I ij , v (cid:63)I ij , p R I ij )4. Compute the vector of semi-discrete, conservative, fluid state variables at the intersection of the edge V i V j and the boundary facet ∂C ij – which is also the midpoint M ij along the edge V i V j (see Figure 3)– as follows W ij = α W i + (1 − α ) W (cid:63)I ij (9)where α = (cid:107) x I ij − x M ij (cid:107) (cid:107) x I ij − x i (cid:107) . Note that (9) implies an interpolation if I ij happens to be between M ij and V j , or an extrapolation otherwise. Also, see [26] for a stabilization of (9) when I ij is very close to V i .5. Finally, compute the numerical flux vector function Φ ij across the boundary facet ∂C ij of the mixedcell C i as follows Φ ij = Φ ij ( W i , W ij , ν ij , EOS) (10)9 igure 3: Construction and solution of a local, 1D, exact, fluid-structure half Riemann problem at the fluid/structure materialinterface (two-dimensional case).
In summary, the second-order FIVER method computes the vector of numerical convective fluxes F (see Equation (6)) on a non body-fitted mesh, away from an embedded discrete surface Σ h , exactly as astandard, second-order, vertex-based, FV method computes F on a body-fitted CFD mesh. In the vicinityof Σ h however, the second-order FIVER method computes the contributions to F of each mixed cell asdescribed in Equations (7)–(10). Remark 3.
From the above description of the semi-discretization of the convective fluxes, it follows thatfor inviscid FSI problems, the second-order EBM FIVER – and for this matter, any of the FIVER methodspublished so far – does not introduce, exploit, or rely on the concept of a ghost flow, ghost boundary, or aghost cell (for example, see [36]). This is in sharp contrast with most of the alternative EBMs, IBMs, andGFMs, including those cited in the introduction of this paper.4.1.2. Numerical diffusive fluxes
Like many FV-based methods, the second-order EBM FIVER – and for this matter, any FIVER methodpublished so far – approximates the diffusive fluxes using a FE-like approach – that is, on an element-by-element basis. For this reason, it computes the diffusive numerical fluxes on the primal cells of the embeddingCFD mesh – that is, on the “elements” of this mesh rather than its control volumes.In each primal cell of an embedding CFD mesh located away from the material interface Σ h , the second-order EBM FIVER computes the contributions to the vector of numerical diffusive fluxes G exactly as astandard, second-order, vertex-based, FV method computes these contributions on a body-fitted CFD mesh.In a mixed (primal) cell however, the second-order EBM FIVER distinguishes between ghost/real, (in-active/active, or occluded/unoccluded) fluid nodes. Essentially, it identifies first the ghost (a.k.a inactiveor occluded) fluid nodes and populates the velocity vector v and temperature T at these nodes of the em-bedding CFD mesh using the standard mirroring technique [36], or a combination of constant and linearextrapolations from the neighboring real (a.k.a. active or unoccluded) nodes [23]. For example, considerthe primal cell ( V m , V n , V k , V j ) of the scenario illustrated in Figure 4. In this case, the second-order EBMFIVER identifies V m as an occluded fluid node in this primal cell. Then, it populates v and T at this nodeusing a combination of constant and linear extrapolations of their counterpart values at the unoccludedneighboring nodes V i , V j , V k , V l , V n , V p and V q (see [23] for a discussion of the specific details, which dependon the prescribed temperature boundary condition (isothermal or adiabatic material interfaces) and thechosen turbulence model (RANS, LES)). Then, the second-order EBM FIVER computes the contributions10 igure 4: Global approach for populating variables of the fluid state vector at the ghost fluid node V m involving informationextracted from the real nodes V i , V j , V k , V l , V n , V p , and V q , vs. local approach featuring 4 different fluid state vectors at V m ,each computed using information extracted from 1 of the 4 elements connected to this node and to be used in the computationof the contributions of this element to the vector of numerical diffusive fluxes (two-dimensional case). of the mixed primal cell ( V m , V n , V k , V j ) to the vector of numerical diffusive fluxes G exactly as a standard,second-order, FV method computes the contributions to G of any primal cell of a body-fitted mesh. Remark 4.
Like most if not all other EBMs, IBMs, and GFMs, the second-order EBM FIVER uses theconcept of a ghost fluid node and populates variables of the fluid state vector at such a node. However,unlike the classical GFM [37] and many variants, the EBM FIVER populates ghost fluid nodes only in mixedprimal cells, and these arise only one graph distance away from the boundary of an embedded discrete surface.Hence, this EBM populates only a fraction of the ghost fluid nodes populated by GFM and related methods.
Here, a noteworthy update on how to populate variables of the fluid states at ghost fluid nodes of anembedding CFD mesh for the purpose of computing the numerical diffusive fluxes is provided. This updateis motivated by slender embedded discrete surfaces such as those associated with cables which typically leadto mixed primal cell configurations such as those illustrated in Figure 4. In this figure, the population of v and T at the ghost fluid node V m using information extracted from the fluid state vectors at the realneighboring nodes V i , V j , V k , V l , V n , V p , and V q is referred to here as a global approach for populatingvariables of a fluid state vector at a ghost node. This approach is adopted by many EMBs, IBMs, and GFMs(for example, see [21]). For CFD and FSI applications involving a slender embedded discrete surface, it is notappropriate however as it may lead to unphysical computations that promote numerical instability. Indeed,in the aforementioned example, the real fluid nodes { V i , V j , V k , V l } and the real fluid nodes { V n , V p , V q } lie ontwo different sides of Σ h where, even if the fluid medium is the same, the fluid flow may have totally differentstates: for example, it could be laminar on one side but turbulent on the other, or subsonic on one side buttransonic or supersonic on the other, etc. Hence, combining in this case components of the fluid state vectorsextracted from the union of these two sets of nodes, because they happen to be connected to the ghost fluidnode V m , in order to populate at this node counterpart variables of the fluid state vector is inappropriate.It explains the numerical instabilities that are observed when the global approach for populating v and T at a ghost fluid node is applied in CFD and FSI applications involving embedded slender structures.Here, the issue exposed above is addressed by proposing an alternative approach for populating variablesof the fluid state vector at a ghost fluid node, for cable-driven and other FSI problems featuring slenderstructural subsystems. This approach avoids the inappropriate mixing of widely different fluid state vectorswithout having to identify “sides” of an embedded discrete surface. It is labeled here as the local approach– in contrast with the global approach explained above – as it consists of populating at each fluid ghost11ode V m connected to one or multiple mixed primal cells one pair of velocity vector v and temperature T per mixed cell to which it is connected. For example in the case of the scenario illustrated in Figure 4, fourdifferent velocity-temperature pairs ( v , T ) are populated at V m as follows: • One pair ( v , T ) computed from information extracted from the fluid state vectors at the real nodes V i , V j , and V l , to be used for computing the contributions of the mixed primal cell ( V i , V j , V m , V l ) to thevector of numerical diffusive fluxes G . • One pair ( v , T ) computed from information extracted from the fluid state vectors at the real nodes V j , V k , and V n , to be used for computing the contributions of the mixed primal cell ( V j , V k , V n , V m ) to thevector of numerical diffusive fluxes G . • One pair ( v , T ) computed from information extracted from the fluid state vectors at the real nodes V l , and V p , to be used for computing the contributions of the mixed primal cell ( V l , V m , V p , V o ) to thevector of numerical diffusive fluxes G . • One pair ( v , T ) computed from information extracted from the fluid state vectors at the real nodes V n , V q , and V p , to be used for computing the contributions of the mixed primal cell ( V m , V n , V q , V p ) to thevector of numerical diffusive fluxes G . A critical component of the master-slave kinematic approach for cable-driven FSI described in Section 3is the computation of the force vector f nS ji at each slave node S ji as specified in (2), from which the force andmoments vectors at the master point M i , f nM i and m nM i (3), respectively, and the generalized loads f nN j and m nN j (4) can be deduced. This force vector is best evaluated by performing the integration in (2) directlyon the embedded, discrete surface Σ h [23] and approximating it using a Gauss quadrature rule. Such anapproximation can be written as f nS ji = (cid:90) Σ nh (cid:16) − p n n n + τ n n n (cid:17) φ S ji d Σ h ≈ n G (cid:88) k =1 (cid:36) k (cid:16) − p nG k n nG k + τ nG k n nG k (cid:17) [ φ S ji ] G k (11)where (cid:36) k denotes the k -th weight of the quadrature rule, and the subscript G k designates the evaluation ofa quantity at the k -th quadrature point G k associated with the weight (cid:36) k . Each Gauss quadrature point G k ∈ Σ h can be located in some primal cell C l of the embedding CFD mesh. Hence, p nG k and τ nG k are typicallycomputed by interpolation of fluid state variables attached to the nodes of C l . However, C l is typically amixed primal cell connected to ghost fluid nodes – that is, nodes of the embedding CFD mesh that areoccluded by Σ h (for example, see Figure 5) – and the populated the fluid state variables at these ghost nodesare typically tainted by nonphysical values due to the underlying extrapolations (see [26]). Consequently,the approximation (11) is typically fraught with spurious oscillations. For this reason, it was proposed in [26]to shift each Gauss point G k used in the approximation (11) in the direction of the outward normal n G k tothe material interface Σ h , up to the point G (cid:48) k defined by (see Figure 5) x G (cid:48) k = x G k + h n G k where h is the characteristic mesh size of the fluid primal cell C l .Indeed, the shifted Gauss point G (cid:48) k will be generally located in this case in a primal cell without anyghost fluid node, the interpolations of p nG (cid:48) k and τ nG (cid:48) k will not be tainted by any nonphysical data, and unlikethe approximation (11), the quadrature approximation f nS ji = (cid:90) Σ nh (cid:16) − p n n n + τ n n n (cid:17) φ S ji d Σ h ≈ n G (cid:88) k =1 (cid:36) k (cid:16) − p nG (cid:48) k n nG (cid:48) k + τ nG (cid:48) k n nG (cid:48) k (cid:17) [ φ S ji ] G (cid:48) k (12)will not exhibit any spurious oscillation. 12 igure 5: Shifting by a distance h of a Gauss point G k used for evaluating the flow-induced forces on the material interface Σ h (two-dimensional case). Furthermore, from ∂p G k ∂ n ≈
0, due to the conservation of momentum in the normal direction at thematerial interface, and the Taylor expansion of the pressure around the point G k , it follows that p G (cid:48) k = p G k + O ( h )The above result implies that the contribution of the pressure term to the approximation (12) is a third-orderapproximation of its counterpart for (11). However, a similar Taylor expansion shows that the contributionof the shear stress traction to the approximation (12) is a second-order approximation of its counterpartfor the approximation (11). Hence, the loads evaluated using the approximation (12) are a second-orderapproximation of the loads evaluated using (11). When the embedded discrete surface undergoes large motions and/or deformations, maintaining reason-able boundary layer resolution around Σ h requires special effort in the Eulerian computational framework forFSI. At high Reynolds numbers, it requires for this purpose AMR [38, 39, 40, 41, 42]. All FIVER methods,including the second-order FIVER method described above are equipped with a local edge refinement andcoarsening algorithm based on the newest vertex bisection (NVB) method [39, 40, 27], which enables theboundary layer and flow features to be efficiently tracked using a wall distance estimator and a Hessian-basederror indicator, respectively. For a cable subsystem however, fully resolving its boundary layer is generallyunaffordable, especially when the cable has a large length-to-diameter ratio.To address the issue raised above, a new criterion for marking edges for refinement is proposed here.This criterion constitutes a lightweight alternative to the wall distance criterion originally proposed in [27].Specifically, wherever an edge of the embedding fluid mesh Ω Fh is intersected twice by the embedded discretesurface Σ h – which indicates that the computational fluid domain is underresolved in this region – theedge is selected for refinement and subsequently bisected. This criterion leads to an affordable mesh. Mostimportantly, it enables capturing the cable-driven FSI effects. This new edge refinement criterion, which isreferred to in the remainder of this paper as the doubly-intersected edge criterion, is equally applicable andequally important for FSI problems where the structural system contains small-scale subsystems that aretypically unresolved by practical CFD meshes. It is illustrated in Section 5.2, where its effectiveness is alsodemonstrated. 13 . Applications The dressing approach and the proposed master-slave kinematic approach for cable-driven FSI are im-plemented in the nonlinear, fluid-structure simulation platform AERO Suite (for example, see [16]). Thissoftware suite includes, among other modules: the versatile, nonlinear, structural analyzer AERO-S; and thecomprehensive, compressible, flow solver AERO-F. It has been validated for numerous wind tunnel, flighttest, and other fluid-structure configurations pertaining to various underwater systems, high performancecars, and aircraft. It supports both the ALE and Eulerian computational frameworks for FSI.Here, the Eulerian computational framework of AERO Suite and its EBM FIVER are applied to thesimulation of two different FSI problems in order to illustrate, verify, and validate the proposed master-slavekinematic approach for cable-driven FSI computations: • An airborne refueling model problem, where the cable subsystem consists of a single flexible hose thatundergoes large motions, and FSI is captured using implicit-implicit fluid-structure computations. Inthis case, the dressing approach is used to verify the accuracy of the proposed master-slave kinematicapproach. • A dynamic, supersonic parachute inflation problem where the cable subsystem consists of 80 suspensionlines, the canopy and suspension lines undergo large motions and deformations, and FSI is capturedusing implicit-explicit fluid-structure computations due to the massive self-contact experienced by thecanopy during the inflation process. The dressing approach is not suitable for this problem for thereason discussed in Section 2. The proposed master-slave kinematic approach is validated in this caseusing real flight data.
First, an airborne refueling model problem [11, 12, 13] consisting of predicting the dynamic, aeroelasticresponse of a flexible hose that trails from an aircraft tanker to the unsteadiness of the flow surrounding it isconsidered. This problem, which is graphically depicted in Figure 6, is designed to compare the performanceof the proposed master-slave kinematic approach for modeling cable-driven FSI to that of the dressingapproach. The flexible hose has a length of 8 m, and a length-to-diameter ratio of
L/D = 119 .
4. Its linearelastic material and other geometrical properties are given in the top part of Table 1. The high-speed inflowconditions are given in the bottom part of this figure: they correspond to an altitude of 8 km above sea level.Air viscosity is modeled using Sutherland’s viscosity law with the constant µ = 1 . × − kg m − s − and reference temperature T = 110 . . × . Hence, the flow is assumed to have already transitioned to the turbulent regime, which is modeledhere using the Spalart-Allmaras turbulence model [43].Parameter Description Value L Length 8 m m L Mass per unit length 0.38 kg m − D Diameter 0.067 m E Young’s modulus 17 MPa ν Poisson’s ratio 0.42 H Altitude 8 km ρ ∞ Free-stream density 0.58 kg m − p ∞ Free-stream pressure 40000 Pa T ∞ Free-stream temperature 240 K M ∞ Free-stream Mach number 0.5
Table 1: Geometrical and material properties of the hose [11, 13] (top); and inflow conditions (bottom). igure 6: Airborne refueling model problem: hose and embedding computational fluid domain. The hose is discretized by 100 geometrically nonlinear beam elements, and its surface is represented bya uniform discretization characterized by 1,200 triangular elements and a uniform, hexagonal, discrete crosssection. It is pinned at one end – that is, all 3 translational DOFs of its FE model are fixed at one end –and unrestrained at the other. Hence, it is anticipated that this flexible hose will undergo large motions inthe presence of a high-speed airstream.The embedding computational fluid domain Ω F (see Figure 6) is chosen to be a 6 m × ×
10 m cube.It is initially discretized by 67,686 tetrahedra. During the unsteady FSI simulation, this embedding CFDmesh is adaptively refined or coarsened, as needed, using: a distance-based criterion to track and resolvethe geometry of the hose as well as the boundary layer; a velocity magnitude Hessian criterion to capturethe shedded vortices; and an algorithm based on the newest vertex bisection algorithm [40, 27] to adapt theedges of the CFD mesh. Consequently, the characteristic size of the embedding CFD mesh near the hose isabout 4 × − m – that is, roughly 1/20 th of its diameter. Away from the boundary layer, the minimumedge length of the embedding CFD mesh is set to 1.5 × − m.Time-discretization of the structural subproblem is performed using the implicit midpoint rule, that ofthe fluid subproblem is performed using the 3-point implicit Backward Difference Method (BDF), and bothdiscretizations are coupled using the second-order, time-accurate, implicit-implicit fluid-structure staggeredsolution procedure presented in [44]. Two FSI simulations are performed using the coupling time-step of∆ t F/S = 2 × − s in order to compute the aeroelastic response of the hose in the time-interval [0 , . s :one using the dressing approach; and another using the proposed master-slave kinematic approach.Figure 7 illustrates the time-evolution of the adapted embedding CFD mesh at the cross section of Ω F corresponding to the free end of the hose, and the velocity magnitude of the fluid flow and displacement ofthe hose at two different cross sections of the computational fluid domain: that corresponding to the middlesection of the hose; and that corresponding to the free end of the hose where the structural motion is thelargest. It reveals that in response to the high-speed flow, the hose drifts while interacting with the trailingvortices. This figure also shows that AMR effectively tracks the boundary layer and flow features.More importantly, Figure 7 shows that both the dressing and master-slave kinematic approaches deliveralmost the same fluid and structural responses. This conclusion is supported by Figure 8 which shows thatthe FSI simulation equipped with the master-slave kinematic approach produces the same time-histories ofthe lateral displacement of the hose and the drag force acting on it as its counterpart equipped with the15 igure 7: Structural motion and velocity magnitude field in cross sections of Ω F computed using: the dressing approach (top);and the master-slave kinematic approach (bottom) – Solution snapshots and adapted embedding CFD mesh at to t = 0 .
01 s, t = 0 .
02 s, t = 0 .
03 s, and t = 0 .
04 s (left to right): cross section of Ω F at the pinned end of the hose (top); cross section of Ω F at the middle section of the hose (middle); and cross section of Ω F at the free end of the hose (bottom). dressing approach. Figure 8: Time-histories of the drag (left) and lateral displacements of the hose at different sections (right) computed using:the dressing approach (blue); and the alternative master-slave kinematic approach (orange).
Next, the FSI simulation of the inflation dynamics of a Disk-Gap-Band (DGB) parachute system inthe low-density, low-pressure, supersonic Martian atmosphere is considered [45]. While such a simulationis crucial to the understanding of the effects of the FSI driven by the suspension line subsystem on the16erformance of the main parachute during the deceleration process, its main purpose here is to validate theproposed master-slave approach for cable-driven FSI using flight data from the landing on Mars of NASA’srover Curiosity.To this end, the DGB parachute system that successfully deployed in 2012 for the Mars landing ofCuriosity is considered (see Figure 9-left). This aerodynamic decelerator system consists of three maincomponents [45]: • The canopy, which is made of F-111 nylon. • The suspension lines, which are made of Technora T221 braided cords. • And the reentry vehicle.Its exact geometric and material properties are listed in Table 2.The simulation discussed herein starts from the line stretch stage where the suspension line subsystemis deployed but the canopy is folded (see Figure 9-right), and the entire system is prestressed by the foldingpattern. The incoming supersonic flow is at the state defined by M ∞ = 1 . ρ ∞ = 0 . − , and p ∞ = 260 Pa. Figure 9: Dynamic supersonic parachute inflation problem: system configuration (left); and embedding computational fluiddomain as well as embedded initial folded configuration (right).
Since the Martian atmosphere is mainly composed of carbon dioxide, the viscosity of this gas is modeledusing Sutherland’s viscosity law with the constant µ = 1 . × − kg m − s − and the reference temperature T = 240 K. The Reynolds number based on the canopy diameter is 4 . × . Hence, the flow is assumed tohave transitioned to the turbulent regime, which is modeled here using Vreman’s eddy viscosity subgrid-scalemodel for turbulent shear flow [47] with model constant C s = 0 . D Diameter 15.447 m t Thickness 7.607 × − m E Young’s modulus 0 .
945 GPa ν Poisson’s ratio 0.4 ρ C Density 1154.25 kg m − α Porosity 0.08Suspension lines L Length 36.56 m D Diameter 3.175 × − m E Young’s modulus 29.5 GPa ρ SL Density 1154.25 kg m − Table 2: True geometrical and material properties of a DGB parachute system [46, 45].
The aforementioned computational fluid domain is a box of size 200 m ×
160 m ×
160 m. It is initiallydiscretized by a mesh composed of Kuhn simplices [40, 27]. Specifically, this initial tetrahedral mesh contains2,778,867 nodes and 16,308,672 tetrahedra. During the FSI simulation reported below, AMR [49, 40, 27]is applied to track and resolve, as in the the case of the airborne refueling model problem, the boundarylayer and flow features. The doubly-intersected edge criterion mentioned in Section 4.3 is applied to eachsuspension line to resolve the flow around it. The specified characteristic mesh sizes near the reentry vehicle,suspension lines, and canopy are 2 . D = 3 . , .
8] s. The length of this time-interval is such that it covers the inflation process aswell as a few breathing cycles of the DGB parachute system. As stated above, the explicit central differ-ence time-integrator is applied to advancing in time the semi-discrete structural subsystem. On the otherhand, the implicit, 3-point BDF scheme is applied to time-integrate the semi-discrete fluid state. For thisreason, the fluid and structural discretizations are coupled for this simulation using the stability-preserving,second-order, time-accurate, implicit-explicit fluid-structure staggered solution procedure presented in [44].The fluid-structure coupling time-step is set to ∆ t F/S = 10 − s.Figure 10 graphically depicts the time-evolutions of the dynamic inflation of the DGB parachute and theflow Mach number around it. It shows that disturbed by the wake generated by the reentry vehicle, the bowshock in the front of the canopy vibrates along with the breathing cycles of the canopy. It also reveals thata jet-like flow is ejected through the canopy vent and gaps between the band and disk gores, and interactswith the turbulent wake behind the canopy. The parachute is fully inflated at approximately t = 0 .
23 s:after this time, it starts the breathing cycles expected from a violent, high-speed, dynamic, inflation process.Although the suspension lines are in this case very slender, they increase the geometric blockage in frontof the canopy, which slightly reduces the flow speed. Moreover, the contributions of these suspension lines toFSI can be observed in Figure 10-e to generate shock waves that alter the wake behind the reentry vehicle,18 igure 10: Time-evolutions of the deployment of the parachute DGB system and the associated flow Mach number field. and even disrupt the bow-shock. These observations are in perfect agreement with those documented in [4]for sub-scale parachute inflation experiments conducted in wind tunnels.A second FSI simulation is performed, but without the contribution of the suspension line subsystemto FSI. In other words, the dynamics of the suspension line subsystem are accounted for in the overall FEstructural model, but the effects of this subsystem on the fluid flow and associated FSI are ignored.Figure 11 reports the time-histories of the total drag force predicted by both FSI simulations describedabove. For reference, this figure also includes: the contribution of the suspension line subsystem to the totaldrag, for the first FSI simulation where this contribution is accounted for; and the time-history of the totaldrag generated by the parachute system of NASA’s rover Curiosity as measured during Mars landing [45].Regarding the first FSI simulation, the reader can observe that the drag force generated by the suspensionline subsystem is relatively negligible. From the total drag result of the second FSI simulation however,the reader can infer that even though the suspension line subsystem is directly responsible for a negligiblepart of the total drag, its effect on the total drag generated by the DGB parachute system is significant.Indeed, the disturbances created ahead of the canopy by the shocks induced by this subsystem disrupt thebow shock of the overall DGB parachute system, mix the high pressure and low pressure flows, and asshown in Figure 11, reduce the total drag force. Most importantly, the reader can also observe that thetime-evolution of the total drag predicted by the first FSI simulation, which accounts for the effect of thesuspension line subsystem on FSI, is in good agreement with the counterpart data measured during theMars landing of Curiosity (relative error of less than 10%): this provides one successful validation of themaster-slave kinematic approach proposed in this paper for modeling cable-driven FSI.19 igure 11: Time-histories of the total drag generated during the dynamic, supersonic parachute inflation process: NASA’s roverCuriosity data [45] (blue); second FSI simulation without the contribution to FSI of the suspension line subsystem (orange);first FSI simulation with the contribution to FSI of the suspension line subsystem (green); and contribution of the suspensionline subsystem to total drag extracted from the first FSI simulation (red).
6. Summary and conclusions
In this paper, an embedded boundary approach for computing the Fluid-Structure Interaction (FSI) ofcable systems or subsystems is presented. In this approach, the dynamics of a solid cable are captured by aFinite Element (FE) semi-discretization of the cable centerline based on conventional beam or cable elements;however, the exact geometry of the cable is also represented, via a discretization of the surface of the cablethat is embedded in the fluid mesh. The embedded discrete surface is endowed with kinematics, and theseare slaved to the motion and deformation of the master beam/cable elements. The flow-induced loads arecomputed on the embedded discrete surface and an energy-conserving method is applied to transfer theseloads to the FE representation of the centerline of the cable. The proposed approach for capturing cable-driven FSI is incorporated in an Eulerian framework for FSI equipped with adaptive mesh refinement andapplied to the solution of two highly nonlinear, turbulent FSI problems: an aerial refueling model problem;and a challenging, dynamic, supersonic parachute inflation problem associated with the Mars landing ofNASA’s rover Curiosity. For the first problem, the proposed approach is shown to reproduce the sameresults as an alternative approach known as the “dressing” method. For the second problem for which thedressing method is inapplicable, the proposed approach performs very well: it reproduces the total draggenerated by the canopy and suspension lines of the parachute system and measured in flight during theMars landing of Curiosity.
Acknowledgments
The authors acknowledge partial support by the Jet Propulsion Laboratory (JPL) under Contract JPL-RSA No. 1590208, and partial support by the National Aeronautics and Space Administration (NASA) underEarly Stage Innovations (ESI) Grant NASA-NNX17AD02G. Daniel Huang also thanks Trond Kvamsdal forenlightening discussions on this topic. 20 eferences [1] Fan Y, Xia J. Simulation of 3D parachute fluid–structure interaction based on nonlinear finite elementmethod and preconditioning finite volume method.
Chinese Journal of Aeronautics
International Journal for Numerical Methods in Fluids
Journal of Spacecraft and Rockets
Advances in Space Research
Computers & Fluids
Journal of Offshore Mechanics and Arctic Engineering
Journal ofGuidance, Control, and Dynamics
InternationalJournal of Solids and Structures
SIAM Journal onScientific Computing
AIAA Journal
International Journalfor Numerical Methods in Fluids
Journal of ComputationalPhysics
Journal of Computational Physics
Journal of Computational Physics
Journalof Computational Physics
Commun. Comput. Phys.
International Journal forNumerical Methods in Fluids
International Journal for Numerical Meth-ods in Engineering
Journal of Computational Physics
Journal of Computational Physics
Inter-national Journal for Numerical Methods in Fluids
Journal of Computational Physics
Computer Methods in Applied Mechanics and Engineering
Journal of Computational Physics
Computer Methods in Applied Mechanics and Engineering
User’s Reference Manual,
Computer Methods in Applied Mechanics and Engineering
Journal of Compu-tational Physics
Journal of Computational Physics
Computational gasdynamics . Cambridge university press . 1998.[37] Fedkiw RP, Aslam T, Merriman B, Osher S. A non-oscillatory Eulerian approach to interfaces in mul-timaterial flows (the ghost fluid method).
Journal of Computational Physics
Journal of Computa-tional Physics
SIAM Journalon Scientific Computing
Mathematicsof Computation
Journal ofFluids Engineering
Journal ofComputational Physics
International Journal for Numerical Methods in Engineering
Journal of Spacecraft and Rockets
Physics of Fluids
Computer Methodsin Applied Mechanics and Engineering
Unified multilevel adaptive finite element methods for elliptic problems . PhD thesis. Uni-versity of Illinois at Urbana-Champaign, ; 1988.[50] Huang DZ, Wong ML, Lele SK, Farhat C. A homogenized flux-body force approach for modeling porouswall boundary conditions in compressible viscous flows. arXiv preprint arXiv:1907.09632arXiv preprint arXiv:1907.09632