Log In Sign Up

An Embedded Boundary Approach for Resolving the Contribution of Cable Subsystems to Fully Coupled Fluid-Structure Interaction

Many engineering systems contain cables as subsystems including suspension lines for parachutes, cables in suspended bridges, risers in offshore platforms, airborne refueling systems, and so on. However, the interactions between fluid and cable subsystems receive little attention in the open literature. This work proposes an embedded surface approach, in which the dynamics of the cable is captured by beam elements typically found in finite element structural models, and the geometry of the cable is represented by an embedded surface. It is built on: master/slave kinematics between beam elements (master), and the embedded surface (slave); a highly accurate algorithm for computing the embedded surface displacement based on the beam displacement; and an energy-conserving method for transferring distributed forces and moments acting on the nodes of the discrete surface to beam elements. Hence, both flow-induced forces on the cable and effect of the structural dynamic response of the cable on the nearby flow are taken into account. Moreover, the proposed model can be easily incorporated in the Eulerian computation framework, which enables handling large deformations of the cable subsystem. Finally, the effectiveness of the proposed model is demonstrated using a model airborne refueling system and a challenging supersonic parachute inflation problem.


Fluid-beam interaction: Capturing the effect of embedded slender bodies on global fluid flow and vice versa

This work addresses research questions arising from the application of g...

An electromechanically coupled beam model for dielectric elastomer actuators

In this work, the Cosserat formulation of geometrically exact beam dynam...

Machine Learning Surrogates for Predicting Response of an Aero-Structural-Sloshing System

This study demonstrates the feasibility of developing machine learning (...

An Image-Based Fluid Surface Pattern Model

This work aims at generating a model of the ocean surface and its dynami...

A consistent co-rotational formulation for aerodynamic nonlinear analysis of flexible frame structures

The design of structures submitted to aerodynamic loads usually requires...

Immersed boundary simulations of fluid shear-induced deformation of a cantilever beam

We derive a mathematical model and the corresponding computational schem...

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 fluid-structure 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 vortex-induced 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, flow-induced 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 vortex-induced vibrations on deepwater risers. One-way coupling is applied in tezduyar2010space, which constructed empirical flow-induced 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 body-fitted 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 body-fitted 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 cable-driven 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 energy-conserving 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 Fluid-Structure Interaction (FSI) problems. Both “dressing” approach and embedded boundary approach are incorporated into the in-house Eulerian computational framework FIVER wang2011algorithms; lakshminarayan2014embedded; main2017enhanced; huang2018family—- A Finite Volume method with Exact two-material 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 one-dimensional beam elements or cable elements, thanks to the large length diameter ratio, but the fluid domain is generally three dimensional. The detection of these one-dimensional elements and defining proper structure normals are ambiguous in this scenario. Hence, it is challenging to capture the local fluid-structure interaction near the cable subsystem. Inspired by the fish-bone aeroelastic models, the “dressing” approach considers a ”super” beam element (See fig. 1-left), which consists of the one-dimensional beam elements, massless rigid beam elements, and a massless discrete surface . The discrete surface represents the real geometry of the cable, where the fluid-structure interaction happens. Namely, the flow-induced 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 very-flexible aircrafts patil2000nonlinear; palacios2010structural and turbine blades malcolm2003modeling; wang2014nonlinear.

The ”dressing” approach is able to handle the fluid-cable 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 fluid-structure interaction problems, like contact problems or damage mechanics.

Figure 1: Schematics of a ”super” beam element from the ”dressing” approach (left) and the ”matched” beam from the embedded boundary approach (right). The discrete surface , representing the real cable geometry, encloses the one-dimensional beam elements (red). They are connected by the massless rigid beam elements (green) in the ”dressing” approach (left).

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 non-conforming fluid/structure interfaces in aeroelastic computations. In the current scenario, the structure solver solves only the dynamics of the one-dimensional 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 one-dimensional beam elements and the discrete surface, a master-slave 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 many-to-one map. The distance vector between these two locations in the initial configuration is denoted by

    , which is written as

  • At each time step, the displacement and velocity of the slave node are formulated as


    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,


    where is the number of slave nodes that are matched to .

The master-slave 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 fluid-structure boundary , the virtual work of the fluid traction acting on can be written as


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


Finally, noting that eq. 5 is the virtual work of the one-dimensional beam structure, therefore, the energy is conservative in the decoupled master-slave procedure. Moreover, the master-slave 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 cross-sections 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 cable-subsystem 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 Eulerian-Lagrangian (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 cable-subsystem 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 Navier-Stokes equations can be written with Einstein notation as


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 cable-subsystem, are written in a Lagrangian formulation with Einstein notation


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 Piola-Kirchhoff 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 Piola-Kirchhoff 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




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 in-house Eulerian computational framework— The Finite Volume method with Exact two-material Riemann Problems (FIVER) wang2011algorithms; farhat2008higher; main2017enhanced; huang2018family. This framework has been extensively verified and validated for complex multi-material problems, including FSI problems with large deformations farhat2014ale, turbulent flows lakshminarayan2014embedded, or dynamic fractures wang2015computational.

4.2.1 No-slip boundary condition

The no-slip 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 fluid-structure 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 flow-induced 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 flow-induced 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


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.

Figure 2: Evaluation of force loads on . and are the Gauss point and the shifted Gauss point associated to a structure boundary or fluid structure interface element.

The pressure force is computed based on the value on the shifted Gaussian point


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


The components corresponding to the derivatives in the tangential directions are computed based on the finite element approximation


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


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.

Figure 3: Airbone refueling system: problem setup.
Parameter Description Value
Length  m
Mass of unit length
Outer diameter  m
Inner diameter  m
Young’s modulus  MPa
Poisson’s ratio
Altitude  km
Pressure  Pa
Temperature  K
Table 1: Airborne refueling system: material properties of the hose ro2010modeling; styuart2011numerical and the inflow conditions.

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 distance-based and Hessian-based 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.

Figure 4: Airbone refueling system: snapshots of the velocity magnitude field computed using the ”dressing” approach (top) and the embedded boundary approach (bottom): , , , and (left to right).

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 second-order, time-accurate, implicit-implicit fluid-structure staggered solution procedure developed in the work of farhat2010robust, and the fluid-structure coupling time-step is  s. These FSI simulations start from the same quasi-steady 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 high-speed 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 x-displacements of the bottom node and the middle node obtained by both approaches are reported in fig. 5. Great agreements are achieved by both approaches.

Figure 5: Drag histories (left) and x-displacements (right) on the bottom nodes and middle nodes of the airborne refueling hose predicted by the ”dressing” approach (blue) and the embedded boundary approach (yellow).

5.2 Supersonic parachute inflation dynamics

Next, the inflation dynamics of a Disk-Gap-Band (DGB) parachute in the low-density, low-pressure 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 F-111 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 free-stream density and pressure conditions are similar to those observed in the standard Martian atmosphere.

(a) Parachute geometry.
(b) Embedded CFD mesh and initial (deflated) configuration of the parachute.
Figure 6: Supersonic parachute inflation dynamics: problem setup.
Component Properties
Canopy Diameter =
Thickness =
Porosity =
Suspension lines Length =
Diameter =
Density =
Table 2: Geometrical and material properties of the parachute lin2010flexible; cruz2014reconstruction.

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 thin-shell 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 F-111 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 cross-section 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 steady-state solution at is computed (See fig. 7-a). Next, the fluid state is initialized using this quasi steady-state 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 second-order, time-accurate, implicit-explicit fluid-structure staggered solution procedure developed in the work of farhat2010robust and the fluid-structure coupling time-step 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. Jet-like 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. 7-e) that alter the wake behind the reentry vehicle, and even disrupt the bow-shock, which agrees with the observations in sengupta2009supersonic from sub-scale parachute inflation experiments.

Figure 7: Time-evolution of the parachute and the flow Mach number.

The time-histories of the drag force computed by post-processing 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 bow-shock, 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 .

Figure 8: Drag histories of the supersonic parachute inflation dynamics problem: curiosity rover data cruz2014reconstruction (blue bold line), simulation with suspension line fluid interaction (orange), simulation with suspension line fluid interaction (green), suspension lines (red)

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 under-resolved slender structures, like turbine blades and fuselage.


Daniel Z. Huang, Philip Avery, and Charbel Farhat 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) under Early Stage Innovations (ESI) Grant NASA-NNX17AD02G. Daniel Z. Huang thanks Kvamsdal Trond for helpful discussions.