Reference Map Technique for Incompressible Fluid–Structure Interaction
Abstract
We present a general simulation approach for fluid–solid interactions based on the fullyEulerian Reference Map Technique (RMT). The approach permits the modeling of one or more finitelydeformable 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 neoHookean 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 flagflapping geometry, for which a number of experimental, numerical, and analytical studies have been performed. We extend the flapping analysis beyond the thinflag 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 Eulerianframe 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 finiteelement 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 meshfree 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 onthefly Lagrangian remeshing. In addition, certain common conditions such as incompressibility are easier to implement in an Eulerian form. Lastly, fixedgrid approaches are wellsuited to numerical analysis, such as a von Neumann stability analysis (LeVeque 2007).
The key challenge for a fullyEulerian 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 largedeformation 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 Eulerianframe 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) secondorder 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 nontrivial 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 divergencefree. The method has been extensively developed since Chorin’s original work (Brown et al. 2001). Here, we consider a modern secondorder 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 finiteelement 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 wellsuited 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 regressionbased extrapolation method, which is more stable, simpler, and allows shapes with corners to be considered. Some accuracy tests are provided, demonstrating secondorder spatiotemporal 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 wellsuited to biolocomotion problems and we show an example of this tool by modeling a jellyfishlike swimmer. Another advantage of the method is the ability to perform manybody 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 largedeformation 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 timedependent map from the undeformed configuration to the deformed state in the physical frame at time . The deformation gradient tensor is defined as
(2.0) 
and represents how an infinitesimal element of the solid has been deformed and rotated. From here, a consititutive law
(2.0) 
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
(2.0) 
where in this case, , and is the solid density in the undeformed configuration.
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
(2.0) 
The deformation gradient tensor is computed from the reference map field according to
(2.0) 
from which the Cauchy stress is evaluated. Equations (2.1), (2.1), (2.1), & (2.1) then form a minimal system of equations for finitestrain 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. Fixedgrid 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
(2.0) 
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
(2.0) 
Kinematic viscosity is defined as . The deviatoric stress is then defined as a smooth transition between the fluid and solid stresses,
(2.0) 
where
(2.0) 
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
(2.0) 
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 secondorder 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.
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
(3.0) 
and an intermediate velocity is computed using
(3.0) 
Here, the advective derivatives and are evaluated at the middle of the timestep using a secondorder explicit Godunov scheme, described in Subsec. 3.1. Once the advective derivatives are evaluated, Eq. (3) allows to be computed. This allows the timecentered reference map to be defined as after which is computed using Eq. (3). From here, the Poisson problem for pressure is evaluated using
(3.0) 
Following Almgren et al. (1996); Puckett et al. (1997), Eq. (3) is solved using a finiteelement formulation, described in Subsec. 3.4. After this, the velocity is projected to be divergencefree using
(3.0) 
where the gradient of is evaluated using a secondorder centered difference formula.
3.1 Advective terms
To evaluate the advective terms appearing in Eqs. (3) and (3), a secondorder 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 fourthorder monotonicitylimited scheme of Colella (1985) described in Appendix A.1. Once the gradients are calculated at the center of each cell, edgecentered velocities and reference maps are created at using Taylor expansions to each of the four edges, which are indexed using halfintegers 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
(3.0) 
where Eq. (2.1) has been substituted for the derivative. The extrapolation of the velocity to the right edge is
(3.0) 
where
(3.0) 
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 ,
(3.0) 
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 noslip 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 Markerandcell (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 cellcentered scalar field. We aim to make
(3.0) 
divergence free. Taking the divergence of Eq. (3.0) yields
(3.0) 
which is discretized as
(3.0) 
Edgebased 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 timecentered.
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 Vcycles 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 Vcycle is performed to further improve accuracy. Typically 5–15 Vcycles are required.
3.1.3 Evaluation of the derivative
Once the MAC projection has been performed the timecentered advective term for the velocity and reference maps are evaluated as
(3.0) 
where is a generic field component.
3.2 Level set update and reference map extrapolation
The simulation makes use of a cellcentered 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 halftimestep 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 signeddistance property. To recover the signeddistance 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 rootfinding method can fail, in which case the routine falls back on the firstorder 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 secondorder 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 signeddistance 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 PDEbased update procedures. Identical methods are used to construct from .
3.2.2 Extrapolation
During the construction of the level set function, a list of nonstraddling 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 PDEbased 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 wellsuited 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 :

Set .

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

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

Set .
This procedure is simpler than the PDEbased 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 highfrequency 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
(3.0) 
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 secondorder finite difference formulae
(3.0) 
after which the deformation gradient is evaluated as
(3.0) 
From here, any constitutive law could be used to evaluate the deviatoric stress, . In the current work, we employ the planestrain incompressible neoHookean law,
(3.0) 
where is the smallstrain shear modulus.
3.3.2 Fluid stress
To evaluate the fluid stress, the gradients of the velocity on vertical edges are computed as
(3.0)  
(3.0) 
Equivalent stencils are used to compute velocity gradients on horizontal edges, after which the fluid stress is given by
(3.0) 
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
(3.0) 
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 nonzero contribution from the second term of Eq. (3.0).
Once all edge stresses are computed, the divergence is computed using
(3.0) 
3.4 Finiteelement projection
To solve the Poisson problem in Eq. (3), we make use of a finiteelement 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
(3.0) 
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
(3.0) 
and the second term is
(3.0) 
where
(3.0) 
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
(3.0) 
In addition, performing a von Neumann stability analysis shows that the timestep must be less than or equal to
(3.0) 
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
(3.0) 
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
(3.0) 
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
(3.0) 
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
(3.0) 
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 nondimensionalized 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 nontrivial FSI setting. It consists of a spinning flexible regular sevenpointed 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
(4.0) 
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 E52683 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 Vcycles and require 0.08 s per timestep. The finite element projections take on average 13.48 Vcycles 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 anticlockwise 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 sevenfold 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
(4.0) 
where is the bicubic interpolation of the velocity field, and the time integration is performed using the secondorder improved Euler method (Süli & Mayers 2003).
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
(4.0) 
where
(4.0) 
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 nearimpossible 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 
We calculated normalized error measures with respect to norms
(4.0) 
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 cellcornered, and hence each coarse gridpoint exactly coincides with a reference gridpoint. The velocity field is cellcentered, 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 fluidonly tests, A & B, are the most accurate, exhibiting clear secondorder convergence across all metrics. Results for the solidonly 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 secondorder 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 higherresolution tests. This resulted in solutions that were almost as accurate as the fluidonly tests, with clear secondorder 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 onedimensional 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 firstorder 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 onedimensional 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 rootfinding 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.
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 (noflapping) 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 stabletoflapping transition line:
(4.0) 
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 nonzero 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 nonzero may also play a role.
The behavioral switch from stable to flapping is also quite evident in the longtime 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.
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 thinfilament limit. We can explain this effect at least in part with bluffbody dynamics. When is nonnegligible, 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
The method also admits a simple approach for simulating actuated solids. This feature allows one to assign timedependent 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: