Reference Map Technique for Incompressible Fluid–Structure Interaction

Reference Map Technique for Incompressible Fluid–Structure Interaction

Chris H. Rycroft1,2, Chen-Hung Wu1,3, Yue Yu4, Email address for correspondence:    Ken Kamrin3 Email address for correspondence:

We present a general simulation approach for fluid–solid interactions based on the fully-Eulerian Reference Map Technique (RMT). The approach permits the modeling of one or more finitely-deformable continuum solid bodies interacting with a fluid and with each other. A key advantage of this approach is its ease of use, as the solid and fluid are discretized on the same fixed grid, which greatly simplifies the coupling between the phases. We use the method to study a number of illustrative examples involving an incompressible Navier–Stokes fluid interacting with multiple neo-Hookean solids. Our method has several useful features including the ability to model solids with sharp corners and the ability to model actuated solids. The latter permits the simulation of active media such as swimmers, which we demonstrate. The method is validated favorably in the flag-flapping geometry, for which a number of experimental, numerical, and analytical studies have been performed. We extend the flapping analysis beyond the thin-flag limit, revealing an additional destabilization mechanism to induce flapping.

1 Introduction

Fluid–structure interaction (FSI) problems highlight a natural dichotomy in the simulation approaches for solids and fluids, where fluid problems tend to be solved using Eulerian-frame methods (Chorin 1967; Hirt et al. 1974; Versteeg & Malalasekera 1995; Tannehill et al. 1997) and solids with Lagrangian approaches (Zienkiewicz & Taylor 1967; Sulsky et al. 1994; Hoover 2006; Belytschko et al. 2013). An FSI simulation method must therefore bridge the gap between these two perspectives. For example, one set of FSI approaches treats both fluid and solid phases in a Lagrangian frame, with a finite-element representation in the solid and an adaptive Lagrangian mesh in the fluid (Rugonyi & Bathe 2001; Bathe 2007; Froehle & Persson 2015), or with both phases treated with a mesh-free approach (Rabczuk et al. 2010). An alternative methodology is to treat the fluid on a fixed Eulerian mesh and the solid with Lagrangian points, such as the family of immersed boundary methods (Peskin 2002; Griffith et al. 2009; Fai et al. 2013).

A fully Eulerian method whereby fluid and solid are both computed on a fixed grid has its advantages. Computation time benefits arise from both phases being treated on a single fixed background grid. The handling of multiple objects interacting or of topological changes to objects can be done with level set fields (Sethian 1999; Osher & Sethian 1988) rather than requiring complex on-the-fly Lagrangian remeshing. In addition, certain common conditions such as incompressibility are easier to implement in an Eulerian form. Lastly, fixed-grid approaches are well-suited to numerical analysis, such as a von Neumann stability analysis (LeVeque 2007).

The key challenge for a fully-Eulerian FSI method is to develop an Eulerian description of the solid. In a small strain limit, this can be achieved by writing the equations of linear elasticity in rate form, referred to as hypoelasticity (Truesdell 1955), which has formed the basis of several numerical techniques (Udaykumar et al. 2003; Rycroft & Gibou 2012; Rycroft et al. 2015). However, here our interest is in developing a large-deformation description of the solid, the more general approach in solid mechanics (Gurtin et al. 2010; Belytschko et al. 2013). In recent years, we have addressed the issue by developing an Eulerian-frame solid simulation approach called the Reference Map Technique (RMT) (Kamrin & Nave 2009; Kamrin et al. 2012; Valkov et al. 2015), which is based on tracking the reference map field—i.e. where material started from—on the Eulerian grid. The reference map field allows the finite deformation of the solid to be computed, from which the material stress is calculated according to a prescribed nonlinear constitutive law. This approach has shown the ability to simulate basic FSI and separately cover a span of desirable features. However, a single implementation covering all needed features for robust physical simulation—e.g. (i) numerical stability, (ii) second-order accuracy in space and time, and (iii) desirable physical traits such as the ability to model incompressible materials, objects with sharp corners, and activated media—has been lacking and non-trivial to produce. In this paper we present such a method and provide a variety of physical simulations using it, which extend our understanding of certain FSI problems.

To represent incompressible solids and fluids we have reformulated the numerical discretization using the projection method framework of Chorin (1967, 1968). In this method, to integrate the velocity field forward by a time step, an intermediate velocity field is computed where the incompressibility constraint is temporarily relaxed. After this, a Poisson problem is solved for the pressure, which is used to project the velocity to be divergence-free. The method has been extensively developed since Chorin’s original work (Brown et al. 2001). Here, we consider a modern second-order implementation described by Yu et al. (2003, 2007) in the context of inkjet printer nozzle simulation. This implementation incorporates a number of improvements, including the treatment of advective terms by Bell et al. (1989), and the approximate projection approach of Almgren et al. (1996) based on a finite-element discretization. We deliberately keep the fluid component of the simulation to match this existing implementation, to emphasize that the reference map technique does not require any special treatment of the fluid. However, we show that the discretization techniques can be generalized to simulate solids via the RMT, and we find that the advective discretization is especially well-suited to simulating the reference map update equations in a fashion more accurate than Valkov et al. (2015).

The projection method removes the Courant–Friedrichs–Lewy (CFL) condition (Courant et al. 1967) associated with pressure waves. This makes it possible to simulate a wide variety of problems in an intermediate Reynolds number regime (and potentially for high Reynolds problems should an adaptive background grid be used). As in Valkov et al. (2015), the level set field representing interface(s) is not explicitly updated, but is tied to where the boundary should be in the reference map field. However, here we switch to a regression-based extrapolation method, which is more stable, simpler, and allows shapes with corners to be considered. Some accuracy tests are provided, demonstrating second-order spatio-temporal accuracy. As a further test of this method for physical simulation, we consider the flag flapping stability problem, which has been studied extensively (Zhang et al. 2000; Watanabe et al. 2002; Zhu & Peskin 2002; Connell & Yue 2007). We are able to quantitatively reproduce the phase plot of stability for a thin flag (Connell & Yue 2007) with very good accuracy for Reynolds numbers below 1000, and reasonably good accuracy for Reynolds numbers above 1000. Our method also makes it possible to simulate flags with substantial thickness, which show a different instability mechanism due to vortex shedding from the tip. The transition between the thin and thick flag behaviors is captured and studied with our method. We also augment the approach to allow internal actuation of the solid bodies. With this addition, the method is well-suited to biolocomotion problems and we show an example of this tool by modeling a jellyfish-like swimmer. Another advantage of the method is the ability to perform many-body contact problems quickly but in a fashion that balances momentum carefully. We demonstrate this approach with an example of many objects of various sizes settling under gravity.

2 Theory

2.1 Overview of the reference map technique

We begin by considering the solid material, which we model using the large-deformation hyperelastic framework (Lubliner 2008; Gurtin et al. 2010). As shown in Fig. 1(a), we introduce an undeformed reference configuration for the solid at time with coordinate system . We then consider a time-dependent map from the undeformed configuration to the deformed state in the physical frame at time . The deformation gradient tensor is defined as


and represents how an infinitesimal element of the solid has been deformed and rotated. From here, a consititutive law


can be used to calculate the Cauchy stress in the physical frame, where represents any internal state variables such as plastic deformation. The material velocity then satisfies


where in this case, , and is the solid density in the undeformed configuration.

(a)(b)Blur zone, Solid, Fluid,
Figure 1: (a) Overview of the hyperelastic framework, whereby an initially undeformed solid with reference coordinate system undergoes a time-dependent mapping to its current configuration at time . (b) Overview of the reference map technique for simulating fluid–structure interaction on a fixed background grid. The sign of the level set function demarcates the solid and fluid phases. The blur zone, used to transition from solid to fluid stress, is defined as the region where .

The most commonly used approach to simulate hyperelastic solids is to introduce a deforming mesh on the solid, and then solve for the nodal displacements, from which (2.1) can be used to compute the stress (Belytschko et al. 2013). However, here we take an alternative approach of introducing the reference map field in the physical frame that represents the inverse mapping of . The field is initialized as , and satisfies the advection equation


The deformation gradient tensor is computed from the reference map field according to


from which the Cauchy stress is evaluated. Equations (2.1), (2.1), (2.1), & (2.1) then form a minimal system of equations for finite-strain hyperelasticity in an Eulerian frame. The reference map and velocity can be represented on a fixed grid. At each timestep equations (2.1) and (2.1) can be used to evaluate the Cauchy stress, after which equations (2.1) and (2.1) can be integrated forward in time. So far, this prescription is general, and could be solved using a variety of discretization approaches such as a finite difference method, finite volume method, or a discontinuous Galerkin method.

The reference map is a standard definition in solid mechanics (Gurtin et al. 2010), and it has been used in problems of inverse design (Govindjee & Mihalic 1996; Fachinotti et al. 2008), but it is not widely employed as a primary simulation variable in the physical frame. Fixed-grid approaches by Plohr & Sharp (1988); Trangenstein & Colella (1991); Liu & Walkington (2001) have been developed that use the deformation gradient tensor as a primary simulation variable.

2.2 Incompressible fluid–structure interaction

In this paper we employ the reference map technique to simulate incompressible fluid–structure interactions. We shall use the terms , , and to refer only to the deviatoric part of the stress, as the pressure field is now deformation independent and separately calculated. We make use of a globally defined velocity field that satisfies the incompressibility constraint


We consider a solid immersed within the fluid, and introduce a level set function (Sethian 1996; Osher & Fedkiw 2003) that is the signed distance to the solid–fluid interface with the convention that in the solid and in the fluid. The reference map is defined within the solid region only.

Let the fluid have density and dynamic viscosity . The fluid stress deviator at any gridpoint is given by


Kinematic viscosity is defined as . The deviatoric stress is then defined as a smooth transition between the fluid and solid stresses,




is a smoothed Heaviside function with a transition region of width . The detailed form of is not important, but the choice in (2.2) has been used elsewhere (Yu et al. 2003, 2007) and is twice differentiable. In order to calculate it is necessary to smoothly extend in the region , which is done using extrapolation methods that will be described in the following section. The density is also defined as a smooth transition between the solid and fluid, as


3 Numerical Method

The numerical procedure is based on the projection method of Chorin (1967, 1968) for solving the incompressible Navier–Stokes equations. Specifically, we consider a modern second-order method described by Yu et al. (2003, 2007) that is especially effective at dealing with advection, and incorporates a number of algorithmic advancements from Chorin’s original treatment.

The simulation domain is a rectangle of size by , and is divided into an grid of rectangular cells of size by . The velocity, the reference map, and the level set, are held at cell centers and are indexed as , , and , respectively, for and (Fig. 2(a)). The components of the velocity field are written as . Pressures are held at cell corners and are indexed as for and . In addition, the grid is padded by two layers of cells in each direction whose values are populated to enforce different boundary conditions.

Figure 2: (a) Arrangement of the fields within a simulation grid cell. The reference map , velocity , and level set field are held at the cell center, while the pressure is held at the cell corners. The level set field is bracketed to emphasize that it is not time-evolved independently, but is instead derived from the reference map. (b) Arrangement of the edge velocities and reference maps that are computed at the half-timestep to evaluate the advective terms.

Superscripts are used to denote timesteps. To advance the simulation forward from timestep to with interval , the following procedure is used. The reference map field is updated using


and an intermediate velocity is computed using


Here, the advective derivatives and are evaluated at the middle of the timestep using a second-order explicit Godunov scheme, described in Subsec. 3.1. Once the advective derivatives are evaluated, Eq. (3) allows to be computed. This allows the time-centered reference map to be defined as after which is computed using Eq. (3). From here, the Poisson problem for pressure is evaluated using


Following Almgren et al. (1996); Puckett et al. (1997), Eq. (3) is solved using a finite-element formulation, described in Subsec. 3.4. After this, the velocity is projected to be divergence-free using


where the gradient of is evaluated using a second-order centered difference formula.

3.1 Advective terms

To evaluate the advective terms appearing in Eqs. (3) and (3), a second-order explicit Godunov scheme is used. The same scheme is applied to both the velocity and reference map . Throughout this section, we denote to be a generic scalar component of either of these two fields. We also refer the reader to recent work by Jain & Mani (2017), which introduces an alternative numerical treatment for reference map advection.

3.1.1 Godunov upwinding

To begin, the gradients of the reference map and velocity at each cell center are evaluated using the fourth-order monotonicity-limited scheme of Colella (1985) described in Appendix A.1. Once the gradients are calculated at the center of each cell, edge-centered velocities and reference maps are created at using Taylor expansions to each of the four edges, which are indexed using half-integers as shown in Fig. 2(b). As an example, an extrapolation of the reference map to the right edge of the cell (with superscript ) is given by


where Eq. (2.1) has been substituted for the derivative. The extrapolation of the velocity to the right edge is




Equivalent procedures are used to compute extrapolations left, down, and up with superscripts , , and , respectively. To ensure tangential stability the terms with tildes in Eqs. (3.1.1) & (3.1.1) are computed differently using the procedure in Appendix A.2. After this procedure, each edge has velocities and reference maps from the two cells that adjoin it, and a Godunov upwinding procedure is used to select which values to use. On the vertical edge at ,


where is a generic component. Thus if the velocity field points rightward then the components are taken from the left cell, and if the velocity field points leftward then the components are taken from the right cell. The function is used when the two velocities are ambiguous. For the horizontal velocity (Case A), and for all other components (Case B). On an edge where a velocity boundary condition is applied (e.g. a no-slip condition) the corresponding edge velocity is set to exactly match the condition. In this paper we restrict to cases of localized solid objects that do not extend to the boundary and thus we do not apply special boundary condition treatment for edge reference map fields.

3.1.2 Marker-and-cell (MAC) projection

The edge velocities calculated in Subsubsec. 3.1.1 may not be precisely divergence free. We therefore apply an intermediate MAC projection step to ensure that the discrete flux entering any grid cell is exactly zero. Let be the edge velocities, and let be a cell-centered scalar field. We aim to make


divergence free. Taking the divergence of Eq. (3.0) yields


which is discretized as


Edge-based densities appearing in this equation are computed via linear interpolation from the two adjacent grid cells. At boundaries where a velocity boundary condition is applied, any derivative on the left hand side of Eq. (3.0) is omitted if it contains values that are out of range. If a pressure condition is applied, then a Dirichlet condition of is applied, where the factor of two arises because the edge velocities are time-centered.

Equation (3.0) results in a large linear system where is a sparse matrix, is the source term, and is a vector of the components . This is solved with a custom multigrid C++ library that employs multithreading using OpenMP. Since the field typically varies smoothly in time, the initial guess for the multigrid algorithm is computed as a linear interpolation from the previous two timesteps. Multigrid V-cycles are performed until the mean squared element in the residual vector reaches a required tolerance . We assume that velocities and densities are within several orders of magnitude of unity. Then an appropriate scale for an element of the residual is , and a tolerance of is used, where is the machine epsilon for double precision floating point arithmetic. Once the tolerance level is reached, one further V-cycle is performed to further improve accuracy. Typically 5–15 V-cycles are required.

3.1.3 Evaluation of the derivative

Once the MAC projection has been performed the time-centered advective term for the velocity and reference maps are evaluated as


where is a generic field component.

3.2 Level set update and reference map extrapolation

The simulation makes use of a cell-centered level set function for tracking the fluid–solid boundary. The level set routine is stored in a narrow band of points of width surrounding the interface (Sethian 1996; Rycroft & Gibou 2012) that is chosen to be large enough to contain the entire blur zone and perform finite difference calculations at all points in this region. The level set is used to extrapolate the reference map fields in the narrow band, and to calculate the mixing of stress and density according to Eqs. (2.2) & (2.2), respectively. Unlike typical applications of the level set method, the field is not explicitly updated, but is instead continually given by the reference map field using the procedure first described in Valkov et al. (2015).

3.2.1 Level set construction

For a given shape, define a continuous function of the reference map such that for reference map values in the solid, for reference map values outside the solid, and on the interface. During the timestep, the reference map field is computed inside the solid using Eq. (3), from which the half-timestep reference map is defined as . Both fields are extended into the narrow band fluid region using the extrapolation methods described in Subsubsec. 3.2.2.

To construct the new level set function , an auxiliary function is first computed in the narrow band such that . The zero contour of will lie at the fluid–solid interface, but this function itself may not satisfy the signed-distance property. To recover the signed-distance property, the level set is constructed from using the reinitialization procedure described in Rycroft & Gibou (2012). This procedure first considers points that straddle the interface, so that one of their orthogonal neighbors has a value of an opposite sign. Each straddling point is considered. The bicubic interpolant is computed, and the modified Newton–Raphson approach of Chopp (2001, 2009) is used to find the shortest distance vector from each straddling point to the interface , after which the level set function is initialized to . In extremely rare cases the root-finding method can fail, in which case the routine falls back on the first-order method described by Sethian (1996). For further details, see Rycroft & Gibou (2012).

With the straddling points of now initialized, the remaining points are filled in using the second-order fast marching method of Sethian (1996). In the fluid, the positive values are computed in order of increasing value, until reaching a cutoff that defines the width of the narrow band. The same procedure is used to fill in the negative values in the solid, until reaching a cutoff . After this procedure, the level set function is now a signed-distance function inside the narrow band. Note that these routines work reliably even if the function has a loss of regularity as some points: since the entire field is directly constructed, there is no possibility for instabilities to grow over time, as can happen in PDE-based update procedures. Identical methods are used to construct from .

3.2.2 Extrapolation

During the construction of the level set function, a list of non-straddling fluid points sorted in order of increasing value, is constructed, which is used for extrapolating the reference map from the solid into the fluid narrow band. Previous approaches to do this have employed PDE-based methods by defining a normal vector and extrapolating outwards from the object in the direction of (Aslam 2004; Rycroft & Gibou 2012). While these methods are well-suited to mathematical analysis, they require considerable bookkeeping for performing the finite difference calculations of and due to the fields only existing at certain grid locations. In previous work we have found this to be a source of difficulty (Valkov et al. 2015).

In the current work, we make use of the following alternative extrapolation procedure. Consider the points in increasing order of value. For a particular point at physical location :

  1. Set .

  2. Use least-squares regression to fit a linear map using all available reference map values at such that . Weight each value in the regression according to .

  3. If the linear map is degenerate then increment and return to Step 2. Otherwise, continue.

  4. Set .

This procedure is simpler than the PDE-based methods since it does not require extensive bookkeeping. Since the method uses all available values in a neighborhood, this repeated averaging results in substantial blurring if the extrapolation is continued far away from the interface. However, here, only values near the interface are required, and the averaging is beneficial, serving to damp out high-frequency modes that could be the source of instability. In Step 3, degeneracies occur only when the available points are colinear, in which case there is insufficient information to determine the linear map. In this case, Step 4 causes the algorithm to retry using more neighboring points.

3.3 Computation of stress

In order to evaluate the stress divergence terms that appear in Eq. (3) & (3.0), the stresses are first computed on the edges of each grid cell. The stress term in Eq. (3.0) is computed as


where and are the components acting on the vertical and horizontal edges, respectively.

3.3.1 Solid stress

To begin, the components of the Jacobian are computed using the second-order finite difference formulae


after which the deformation gradient is evaluated as


From here, any constitutive law could be used to evaluate the deviatoric stress, . In the current work, we employ the plane-strain incompressible neo-Hookean law,


where is the small-strain shear modulus.

3.3.2 Fluid stress

To evaluate the fluid stress, the gradients of the velocity on vertical edges are computed as


Equivalent stencils are used to compute velocity gradients on horizontal edges, after which the fluid stress is given by


where is the viscosity. Equation (3.0) is our standard approach for computing the fluid stress. However, we have also investigated a simplified calculation. Since , it follows that in the bulk of the fluid, the second term in Eq. (3.0) has zero contribution to . Hence an alternative formula is


This formula only requires evaluating the simpler stencil in Eq. (3.0). However, Eq. (3.0) is not strictly valid in the blur zone since taking the divergence Eq. (2.2) results in a non-zero contribution from the second term of Eq. (3.0).

Once all edge stresses are computed, the divergence is computed using


3.4 Finite-element projection

To solve the Poisson problem in Eq. (3), we make use of a finite-element formulation. The pressure is comprised of piecewise bilinear elements, and the velocity and density are piecewise constant on the grid cells. For a given pressure element the weak formulation of Eq. (3) is


where is the section of the boundary where inflow and outflow conditions are prescribed. Consider a particular bilinear element function located at a pressure point in the bulk of the domain. The first term of Eq. (3.0) is


and the second term is




This linear system is solved using the multithreaded custom C++ geometric multigrid library mentioned in Subsubsec. 3.1.2 using an error tolerance of .

3.5 Parameter choices and stability

Our test cases involve four physical parameters: solid shear modulus , solid density , fluid viscosity , and fluid density . In the solid, the shear wave speed is . The CFL condition requires that the simulation timestep be less than or equal to


In addition, performing a von Neumann stability analysis shows that the timestep must be less than or equal to


in order to resolve the viscous fluid stress. Inside the solid, we find that simulating stress using only Eq. (3.0) results in an instability—this should be expected since we are effectively solving a hyperbolic system using centered finite differences. To rectify this, we incorporate an extra small artificial viscous stress inside the solid. Based on dimensional considerations, the artificial viscosity should satisfy


where is a dimensionless constant. In addition, we also find that artificial viscosity is useful in the fluid–solid transition region. We therefore define the extra viscous stress as


where is a dimensionless constant. Based on a variety of tests, we use and throughout the paper. Since the purpose of this extra stress is to stabilize the numerical system, we employ the simpler form of fluid stress given in Eq. 3.0. Since scales linearly with grid spacing, and the simpler fluid stress only introduces a discrepancy in the blur zone, any errors that are introduced will reduce to zero as the grid is refined. The corresponding timestep restriction is


With these definitions in place, the simulation timestep is chosen to be smaller than the minimum of the three conditions in Eqs. (3.0), (3.0), & (3.0), so that


Here, and are padding factors that are smaller than one. For this paper we use and , so that the restrictions arising from the two physical stresses (I & II) are applied more stringently than the one for the artificial stress (III). Note that in the limit as , the artificial viscosity vanishes.

4 Results

Since our purpose is to demonstrate the numerical method as opposed to apply it to a specific problem, we make use of non-dimensionalized quantities for all of the results that we present. To connect the results to reality, the simulation parameters and output can be multiplied by appropriate length, time, and mass scales. Our results also focus on the case of equal grid spacing, .

4.1 A spinning flexible rotor

The first example demonstrates our method’s ability to handle sharp solid corners within a non-trivial FSI setting. It consists of a spinning flexible regular seven-pointed star that is centered on the origin and has a vertex at for , with density . The resolution is , the simulation domain is and periodic boundary conditions are used. The fluid and rotor are initially stationary. The region is used as a pivot. To excite the fluid, the pivot is rotated with an oscillatory motion with angle . This is done by applying an external tether force to the pivot region of


where and is a rotation matrix with angle . The spring constant is set to , which ensures that the natural frequency of the tether satisfies the timestep restriction imposed by the method.

The simulation was run from to using sixteen threads on a Linux server with two Intel Xeon E5-2683 2.10 GHz processors. For the given parameters, the timestep of was determined by the extra viscous stress in the solid. Simulation output was saved at regular intervals of . The total wall clock time for the simulation was 10.48 h. A total of 114,000 timesteps were performed, with each taking 0.33 s to compute. A substantial fraction of the computation time is spent performing the two linear solves. The MAC projections take on average 14.78 V-cycles and require 0.08 s per timestep. The finite element projections take on average 13.48 V-cycles and require 0.11 s per timestep.

Snapshots of the simulation are shown in Fig. 3. As the star begins to rotate, each point deforms clockwise, and vortices are shed from the points, which are visible at . By , the rotor is stationary, and the points are now deformed anti-clockwise due to the angular deceleration. As time progresses, the disturbance to the fluid becomes larger. By , the rotational symmetry of the fluid flow is lost, due to interactions across the periodic boundaries, which break the seven-fold symmetry. By , after two complete cycles of the oscillatory motion, there are vortices present throughout the fluid. Supplemental Movie 1 shows the complete simulation. To visualize the fluid motion, the movie also shows a number of tracers with trajectories . The tracers are initialized at random positions in the fluid and are updated using the ordinary differential equation


where is the bicubic interpolation of the velocity field, and the time integration is performed using the second-order improved Euler method (Süli & Mayers 2003).

Figure 3: Snapshots of vorticity in a simulation of flexible seven-pointed rotor being spun with an oscillatory motion a fluid. The thick black line marks the fluid–structure interface. The thin dashed lines are contours of the components of the reference map and indicate how the rotor has deformed.

4.2 Tests of convergence and accuracy

To study the accuracy of the numerical method, we performed a convergence test in the periodic domain using an initial incompressible velocity field of




This velocity field is designed to have features with a variety of length scales. We simulated up to , used a shear modulus of , a fluid density of , and employed the standard choices for extra viscosity and timestep selection. Using the same initial velocity field, we ran tests using (i) fluid only, (ii) solid only, (iii) a circle of radius 0.6 centered on , and (iv) a square of side length 1.2 centered on . We also examined the effect of viscosity, and a fluid/solid density ratio, and the scaling of the extra viscous term. The configurations of eight different tests are shown in Table 1.

Due to the complexity of the governing equations, it is near-impossible to write down an analytical solution to compare against for any test configuration. We therefore performed reference simulations using a grid. For each test, we then ran a suite of coarser simulations using grids where to compare against the reference results. Since each divides evenly into , the grid squares of these coarse simulations align with the reference simulations.

Test State CEV
A Fluid only No 2.28 2.29 2.28 2.32
B Fluid only No 2.32 2.32 2.32 2.33
C Solid only 1 No 1.57 1.57 1.56 1.58
C’ Solid only 1 Yes 2.35 2.40 1.13 2.30
D Square 1 No 1.81 1.69 1.11 1.26
E Circle 3 No 1.83 1.79 1.89 1.19
F Circle 1 No 1.80 1.80 1.97 1.21
F’ Circle 1 Yes 1.61 1.63 1.53 1.21
Table 1: Details of the eight convergence tests that were performed with model problem described in the text. Tests C’ and F’ were performed using constant extra viscosity (CEV) whereby the extra viscosity was held constant at the standard value for the lowest resolution grid, , as opposed to scaling linearly with the grid size. The last four columns give the exponents of convergence for velocity and pressure under different norms, based on a linear fit of the three finest-resolution data points in Fig. 4.
Figure 4: Plots showing the accuracy of the solutions for different grid sizes for the eight convergence tests given in Table 1. Accuracy is computed with respect to reference solutions on a grid. Four accuracy measures are shown: the velocity in the , , and norms, and the pressure in the norm.

We calculated normalized error measures with respect to norms


where is the area of the domain, and the ‘ref’ and ‘coa’ subscripts refer to the reference and coarse simulation fields, respectively. The integral is calculated using a direct sum over the field values in the coarser simulation grid. The pressure field is cell-cornered, and hence each coarse gridpoint exactly coincides with a reference gridpoint. The velocity field is cell-centered, so some coarse gridpoints may not align with a reference gridpoint, in which case the reference value is computed using bilinear interpolation. The errors associated with this interpolation are and are small compared to the errors to be measured.

Figure 4 shows convergence plots for the velocity in the , , and norms, plus the pressure in the norm; our discussion focuses on velocity, since the pressure can be viewed as a Lagrange multiplier enforcing the incompressibility constraint. For each set of data, Table 1 lists the corresponding rate of convergence, calculated using linear regression for the data from the three finest grids, . The fluid-only tests, A & B, are the most accurate, exhibiting clear second-order convergence across all metrics. Results for the solid-only test C are less accurate with error measures on the scale of . However, this test is substantially more challenging since the limit involves tracking elastic waves with zero dissipation. The velocity fields converge at order in the , , and norms.

Test C uses the usual procedure (Subsec. 3.5) for choosing extra viscosity, whereby it scales linearly with the grid size. This procedure is consistent with standard numerical schemes; for example, in the second-order Lax–Wendroff method (Lax & Wendroff 1960; LeVeque 2002) the stabilizing diffusive term scales linearly with the grid spacing. However, we also considered an alternative test C’ whereby the extra viscosity was chosen based on the grid and then held constant for the higher-resolution tests. This resulted in solutions that were almost as accurate as the fluid-only tests, with clear second-order convergence in the and norms. However, the convergence rate in the norm is reduced. Inspection of the results shows that that the maximum deviations are localized to a one-dimensional set of points where the velocity components are switching sign, thus resulting in a discrete switch in the timestep update and a lower convergence rate when measured in the norm.

The remaining tests, D, E, F, and F’ all involve fluid–structure interaction. For these tests, the rate of convergence is approximately in the and norms. Inspection of the results shows that the largest deviations occur at the fluid–structure interface. Since the blur zone is defined in terms of grid points, its width shrinks at higher resolution. This involves altering the underlying equations over a region of size , this results in a lower rate of convergence. However, since the fluid and solid discretizations are independently second order, is likely that an alternative boundary treatment—perhaps using a sharp interface approach (Gibou & Fedkiw 2005)—could yield improved results. Test E shows that a fluid–solid density ratio has little effect on the convergence rate. Test D shows that the square geometry does not affect the convergence rate in the and norms, but does result in first-order convergence in the norm due to localized effects at the corners.

4.3 Flag flapping

Besides numerical convergence, as a test of the robustness of our approach and its accuracy across Reynolds numbers, we consider the example of flag flapping, a problem that has been studied extensively from experimental, numerical, and analytical perspectives (Zhang et al. 2000; Watanabe et al. 2002; Zhu & Peskin 2002; Connell & Yue 2007). Following the problem description and notation of Connell & Yue (2007),For consistency with Connell & Yue (2007) we fully adhere to their notation. However we draw attention to the reader that two symbols used in this section, (mass ratio) and (filament thickness), have different meanings than (dynamic viscosity) and (grid spacing) that are used in all other sections. we introduce a thin filament of length , thickness , density , and Young’s Modulus , clamped at its left endpoint and submerged in fluid of kinematic viscosity and density , flowing rightward with speed at infinity. Three dimensionless numbers can be introduced to study the dynamical behavior of the filament: the mass ratio , Reynolds number , and nondimensional bending rigidity . Unlike the previous numerical approaches that consider the filament to be a one-dimensional beam, our method uses a true continuum solid formulation so we can consider cases beyond the thin filament limit, such as a thick flag for which the parameter does not necessarily satisfy .

We first seek to determine if our method correctly captures the transition of the filament dynamics from stable to flapping in the limit of a thin filament. We consider a filament with , , , and . To set , we use the fact that in the linear elastic regime , and the moment of inertia is . We vary and in order to test a range of and Re. The simulation domain is set to be a rectangle with assigned rightward velocity of speed on the left boundary, vanishing pressure on the right boundary, and periodic boundary conditions on the top and bottom boundaries. We use a grid to represent the domain. The filament is modeled as rectangle with semicircular end caps. The filament is clamped at using the tethering methodology described in Sec. 4.1, with and in this case. We track the filament tip by introducing a tracer that starts from and is integrated according to Eq. (4.0). To prevent integration errors building up over time, the position of the tracer is periodically reset to satisfy using a Newton–Raphson root-finding method, where is the bicubic interpolant of the reference map field. The results of this section are based upon 556 simulations with different parameters that were run on a variety of Linux and Apple servers at Harvard University and the Lawrence Berkeley National Laboratory. Depending on parameters and computer speed, each simulation took been approximately 3–12 days using 4–6 threads. Simulations with smaller Re generally take longer, since resolving the fluid viscosity requires a smaller timestep.

Figure 5: Plot showing the steady-state oscillation amplitude of a thin flag with aspect ratio 20 and bending rigidity , as a function of the Reynolds number Re and mass ratio . The colors shown are based on a bilinear interpolation of a two-dimensional grid of simulations. The axis ticks show the sampled values of Re and , with more simulations being performed in parameter ranges of interest. The thin dotted lines are contours at spacings of for . The thick solid line is the stable-to-flapping transition formula (4.0) of Connell & Yue (2007).

To systematically evaluate the behavior of the filament, we take the Fourier transform of the perpendicular deflection of the filament tip over , and output the maximum Fourier amplitude, . If the filament is in the stable (no-flapping) regime and otherwise the filament is flapping, with serving as a scalar measure of the amplitude of the dominant flapping mode. Since our initial conditions are symmetric, the breakage of symmetry occurs due to numerical noise introduced by the multigrid solver, on the scale of the parameters and introduced previously. We also investigated explicitly breaking symmetry by applying an initial perturbation to the perpendicular velocity in the filament tip, but the calculations of were insensitive to this. Since the typical filament oscillation period is approximately 1.7, the simulations correspond to almost one hundred complete oscillations, which is sufficient time for the oscillation amplitude to reach steady state. Connell & Yue (2007) proposed an analytical formula for the stable-to-flapping transition line:


Connell & Yue (2007) validated this equation numerically using a direct fluid–filament coupling procedure, a procedure that itself was validated against experiments (Zhang et al. 2000; Watanabe et al. 2002). In Fig 5 we show the behavior of from our numerical simulations together with the analytical phase boundary above. For Reynolds numbers below 1000, there is very good agreement between the locus where goes non-zero and the analytical curve. When the transition predicted by the simulation happens at a slightly higher than predicted by Eq. (4.0). The most likely explanation for this is that numerical diffusion from the fluid advection effectively increases the fluid viscosity. However, other factors such as the finite domain size, the extensibility of the filament, and the non-zero may also play a role.

Figure 6: Simulations of a thin flexible flag anchored at in a fluid with mean velocity , at . The flag has an aspect ratio of 20. Three simulations with different parameters are shown: (a) stable with , (b) limit-cycle flapping with , and (c) chaotic flapping with . The thick black lines mark the fluid–structure interfaces. The thin dashed lines are contours of the components of the reference map and indicate how the flags deform. The colors show vorticity, using the same scale as Fig. 3.

The behavioral switch from stable to flapping is also quite evident in the long-time flow fields, shown in Fig. 6. Small values of and result in stable behavior, characterized by a straight filament and fluid flow that is symmetric about the filament axis (Fig. 6(a)). Upon crossing the transition, periodic undulatory filament motions develop with a fluid vortex street shed from the filament (Fig. 6(b)). Increasing Re and even further reveals a chaotic filamentary motion, which was also observed in Connell & Yue (2007) (Fig. 6(c), Supplemental Movie 2). The chaotic regime coincides with a drop in shown in the top right of Fig. 5 because the tip deflection no longer has a clear single dominant oscillatory mode. Because the filament is modeled as a thin continuum body of isotropic elastic media as opposed to an inextensible beam, we observe filament extension in the initial moments of the simulation as the imposed fluid flow applies a net rightward traction.

Reynolds number Re(a)Reynolds number Re(b)
Figure 7: Plot showing the steady-state oscillation amplitude as a function of the Reynolds number Re and mass ratio , for (a) flags with and aspect ratio 10, and (b) flags with and aspect ratio 5. The colors shown are based on a bilinear interpolation of a two-dimensional grid of simulations, using the same scale as Fig. 5. The axis ticks show the sampled values of Re and , with more simulations being performed in parameter ranges of interest. The thin dotted lines are contour at spacings of for . The thick solid line in (a) is the stable-to-flapping transition formula (4.0) for thin flags of Connell & Yue (2007). For (b), the formula is out of range and the entire parameter space is in the stable region.
Figure 8: Simulations of a thick flexible flag anchored at in a fluid with mean velocity , at . The flag has an aspect ratio of 5. Two simulations with different parameters are shown: (a) vortex-shedding with and (b) limit-cycle flapping with . The thick black lines mark the fluid–structure interfaces. The thin dashed lines are contours of the components of the reference map and indicate how the flags deform. The colors show vorticity, using the same scale as Fig. 3.

We explore the importance of aspect ratio by introducing as an independent dimensionless group. We observe that as one departs from the regime, adherence to (4.0) is diminished. In Fig. 7 we show results for and . In general, thick flags have a smaller stable domain than would be predicted by the thin-filament limit. We can explain this effect at least in part with bluff-body dynamics. When is non-negligible, the thickness of the flag allows the solid geometry to act as a bluff body over which the fluid is driven to flow. Flow over a fixed cylinder of diameter undergoes a transition from a laminar flow to a periodic vortex street as grows beyond (Lienhard 1966). In our case, the flag thickness acts like , and once a vortex street is induced off the bluff back end of the flag, the oscillatory force it induces necessitates flapping. We reiterate that this physical source of oscillatory forcing emerges only when flags are thick enough to act as a bluff body. Consistent with this expectation, when we see only flapping states for any choice of or Re. Figures 8(a) and 8(b) show simulation snapshots of bulky flags with low and high mass ratios, respectively. Simulations of these two cases are shown in Supplemental Movie 3 and Supplemental Movie 4, respectively.

4.4 Solid actuation

Figure 9: Six successive snapshots of the flapping swimmer (), with colors showing vorticity . A subregion within the solid body is actuated to bend periodically and the remaining solid is passive. The motion induces the flapping body to swim. The thick black lines mark the fluid–structure interface. The thin dashed lines are contours of the components of the reference map and indicate how the swimmer deform.

The method also admits a simple approach for simulating actuated solids. This feature allows one to assign time-dependent internal deformations to subregions of a solid, which is useful for modeling active media such as swimmers. Unlike the tethering approach used in Section 4.1, which assigns the full motion of a region by adding an external body force in that region, here what is done is to add extra internal stresses to achieve a desired shape change in a subdomain, without adding net external force. To actuate a particular (Lagrangian) solid region, , one writes the actuated deformation gradient , which can then be equivalently expressed in Eulerian frame as for the image of in the Eulerian frame. At any point , the constitutive relation is adjusted by replacing all references to with . In an isotropic hyperelastic system, for example, this effectively distorts the region’s rest configuration to the distortional state given by . If at any moment in time a configuration of the actuated domain differs from the intended actuated configuration, a stress given by emerges that moves the system toward the actuated deformation. One could in principle assign a stiffer response in the actuated domain if a faster conformation is desired, but we have found it to be sufficient to use the same underlying hyperelastic constitutive model in the actuated and passive subregions of the solid. This approach is similar to the multiplicative Kroner–Lee decomposition used in plasticity (Kröner 1960; Lee 1969), where a tensorial state variable is introduced and the elastic deformation gradient is given by . But unlike , which evolves under a constitutive flow rule, here we assign directly.

As an example, we consider a flapping swimmer (Fig. 9, Supplemental Movie 5). The swimmer is a rectangle of width and height with circular end caps, initially centered on , which we choose to be the location of the origin. We choose the actuated domain, , to be a centered subregion within the swimmer, comprising a rectangle of width 0.28 and height 0.042 with circular end caps. The following actuation is applied: