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 Reference Chorin1967; Hirt, Amsden & Cook Reference Hirt, Amsden and Cook1974; Versteeg & Malalasekera Reference Versteeg and Malalasekera1995; Tannehill, Anderson & Pletcher Reference Tannehill, Anderson and Pletcher1997) and solids with Lagrangian approaches (Zienkiewicz & Taylor Reference Zienkiewicz and Taylor1967; Sulsky, Chen & Schreyer Reference Sulsky, Chen and Schreyer1994; Hoover Reference Hoover2006; Belytschko et al. Reference Belytschko, Liu, Moran and Elkhodary2013). 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 Reference Rugonyi and Bathe2001; Bathe Reference Bathe2007; Froehle & Persson Reference Froehle and Persson2015), or with both phases treated with a mesh-free approach (Rabczuk et al. Reference Rabczuk, Gracie, Song and Belytschko2010). 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 developed by Peskin (Reference Peskin1972a,Reference Peskinb, Reference Peskin1977, Reference Peskin2002), which have been extensively used to simulate membranes (Griffith et al. Reference Griffith, Luo, McQueen and Peskin2009; Fai et al. Reference Fai, Griffith, Mori and Peskin2013), and obstacles (Coquerelle & Cottet Reference Coquerelle and Cottet2008; Gazzola et al. Reference Gazzola, Chatelain, van Rees and Koumoutsakos2011; Engels et al. Reference Engels, Kolomenskiy, Schneider and Sesterhenn2015) via Brinkman volume penalization.
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 (Osher & Sethian Reference Osher and Sethian1988; Sethian Reference Sethian1999; Osher & Fedkiw Reference Osher and Fedkiw2003) 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 Reference LeVeque2007).
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 Reference Truesdell1955), which has formed the basis of several numerical techniques (Udaykumar et al. Reference Udaykumar, Tran, Belk and Vanden2003; Rycroft & Gibou Reference Rycroft and Gibou2012; Rycroft, Sui & Bouchbinder Reference Rycroft, Sui and Bouchbinder2015). However, here, our aim is to develop a large-deformation description of the solid, the more general approach in solid mechanics (Gurtin, Fried & Anand Reference Gurtin, Fried and Anand2010; Belytschko et al. Reference Belytschko, Liu, Moran and Elkhodary2013). This has attracted interest over the past three decades, with a variety of different approaches being explored in the literature, such as the conservative first-order method of Plohr & Sharp (Reference Plohr and Sharp1988), the deformation gradient-based method of Liu & Walkington (Reference Liu and Walkington2001) and the initial point set (IPS) method of Dunne (Reference Dunne2006). See § 2.2 for further references and a comprehensive discussion.
In recent years, we have contributed to this field by developing an Eulerian-frame solid simulation approach called the reference map technique (RMT) (Kamrin Reference Kamrin2008; Kamrin & Nave Reference Kamrin and Nave2009; Kamrin, Rycroft & Nave Reference Kamrin, Rycroft and Nave2012; Valkov, Rycroft & Kamrin Reference Valkov, Rycroft and Kamrin2015), which is based on tracking the reference map field – i.e. where material started from – as the primary simulation variable 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) good convergence properties for fluid and solid phases 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 (Reference Chorin1967, Reference Chorin1968). 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, Cortez & Minion Reference Brown, Cortez and Minion2001). Here, we consider a modern second-order implementation developed by Almgren, Bell and coworkers (Bell, Colella & Glaz Reference Bell, Colella and Glaz1989; Puckett et al. Reference Puckett, Almgren, Bell, Marcus and Rider1997; Almgren et al. Reference Almgren, Bell, Colella, Howell and Welcome1998). The implementation uses the discretization of Colella (Reference Colella1990) for accurately calculating advective terms. It employs the approximate projection approach of Almgren, Bell & Szymczak (Reference Almgren, Bell and Szymczak1996) based on a finite-element discretization. The implementation was subsequently extended to two-phase flows using a level set method (Sussman et al. Reference Sussman, Almgren, Bell, Colella, Howell and Welcome1999), which was used to simulate an inkjet printer nozzle (Yu, Sakai & Sethian Reference Yu, Sakai and Sethian2003, Reference Yu, Sakai and Sethian2007). In our approach 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 the approach of Valkov et al. (Reference Valkov, Rycroft and Kamrin2015).
The projection method removes the Courant–Friedrichs–Lewy (CFL) condition (Courant, Friedrichs & Lewy Reference Courant, Friedrichs and Lewy1967) 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). Following Valkov et al. (Reference Valkov, Rycroft and Kamrin2015), 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. Two physically motivated examples of a flexible rotor and a paddle are presented, and a variety of convergence and performance tests are provided in appendices B, C, D and E. Our numerical tests show that the method is first-order accurate for a typical fluid–solid simulation, with the interface being the largest source of error. With respect to the 
$L^{2}$ norm, second-order accuracy is achieved for a fluid-only discretization, and for a solid-only discretization with specific choices for numerical damping.
With the properties of the method established, we consider the flag-flapping stability problem, which has been studied extensively (Zhang et al. Reference Zhang, Childress, Libchaber and Shelley2000; Watanabe et al. Reference Watanabe, Suzuki, Sugihara and Sueoka2002; Zhu & Peskin Reference Zhu and Peskin2002; Connell & Yue Reference Connell and Yue2007). We can quantitatively reproduce the phase plot of stability for a thin flag (Connell & Yue Reference Connell and Yue2007) 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 behaviours 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 by modelling a jellyfish-like swimmer. Another advantage of the method is the ability to perform many-body contact problems quickly but in a fashion that accurately balances momentum. We demonstrate this approach with two examples of many objects of various sizes and densities settling under gravity. All computer code to generate the results in this paper is released as open source software (appendix F).
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 Reference Lubliner2008; Gurtin et al. Reference Gurtin, Fried and Anand2010). As shown in figure 1(a), we introduce an undeformed reference configuration for the solid at time 
$t=0$ with coordinate system 
$\boldsymbol{X}$. We then consider a time-dependent map 
$\unicode[STIX]{x1D74C}(\boldsymbol{X},t)$ from the undeformed configuration to the deformed state in the physical frame at time 
$t$. The deformation gradient tensor is defined as
$$\begin{eqnarray}\unicode[STIX]{x1D641}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D74C}}{\unicode[STIX]{x2202}\boldsymbol{X}},\end{eqnarray}$$and represents how an infinitesimal element of the solid has been deformed and rotated. From here, a constitutive law
$$\begin{eqnarray}\unicode[STIX]{x1D748}_{s}=\boldsymbol{f}(\unicode[STIX]{x1D641},\unicode[STIX]{x1D73B})\end{eqnarray}$$ can be used to calculate the Cauchy stress 
$\unicode[STIX]{x1D748}_{s}$ in the physical frame, where 
$\unicode[STIX]{x1D73B}$ represents any internal state variables such as plastic deformation. The material velocity 
$\boldsymbol{v}(\boldsymbol{x},t)$ then satisfies
$$\begin{eqnarray}\unicode[STIX]{x1D70C}\left(\frac{\unicode[STIX]{x2202}\boldsymbol{v}}{\unicode[STIX]{x2202}t}+(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{v}\right)=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748},\end{eqnarray}$$ where 
$\unicode[STIX]{x1D748}=\unicode[STIX]{x1D748}_{s}$ in this case, 
$\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{s}/(\det \unicode[STIX]{x1D641})$ and 
$\unicode[STIX]{x1D70C}_{s}$ is the solid density in the undeformed configuration.
Figure 1. (a) Overview of the hyperelastic framework, whereby an initially undeformed solid with reference coordinate system 
$\boldsymbol{X}$ undergoes a time-dependent mapping 
$\unicode[STIX]{x1D74C}(\boldsymbol{X},t)$ to its current configuration at time 
$t$. (b) Overview of the reference map technique for simulating fluid–structure interaction on a fixed background grid. The sign of the level set function 
$\unicode[STIX]{x1D719}(\boldsymbol{x},t)$ demarcates the solid and fluid phases. The blur zone, used to transition from solid to fluid stress, is defined as the region where 
$|\unicode[STIX]{x1D719}|<\unicode[STIX]{x1D700}$.
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. Reference Belytschko, Liu, Moran and Elkhodary2013). However, here we take an alternative approach of introducing the reference map field in the physical frame 
$\unicode[STIX]{x1D743}(\boldsymbol{x},t)$ that represents the inverse mapping of 
$\unicode[STIX]{x1D74C}$. The field is initialized as 
$\unicode[STIX]{x1D743}(\boldsymbol{x},0)=\boldsymbol{x}$, and satisfies the advection equation
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D743}}{\unicode[STIX]{x2202}t}+(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\unicode[STIX]{x1D743}=\mathbf{0}.\end{eqnarray}$$The deformation gradient tensor is computed from the reference map field according to
$$\begin{eqnarray}\unicode[STIX]{x1D641}=\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D743}}{\unicode[STIX]{x2202}\boldsymbol{x}}\right)^{-1},\end{eqnarray}$$ from which the Cauchy stress is evaluated. Equations (2.2), (2.3), (2.4) and (2.5) then form a minimal system of equations for finite-strain hyperelasticity in an Eulerian frame. The reference map 
$\unicode[STIX]{x1D743}(\boldsymbol{x},t)$ and velocity 
$\boldsymbol{v}(\boldsymbol{x},t)$ can be represented on a fixed grid. At each time step, equations (2.5) and (2.2) can be used to evaluate the Cauchy stress, after which (2.3) and (2.4) 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.
2.2 Other related approaches
The reference map is a standard definition in solid mechanics (Gurtin et al. Reference Gurtin, Fried and Anand2010), and it has been used in problems of inverse design (Govindjee & Mihalic Reference Govindjee and Mihalic1996; Fachinotti, Cardona & Jetteur Reference Fachinotti, Cardona and Jetteur2008), but it is not widely employed as a primary simulation variable in the physical frame. Fixed-grid approaches by Plohr & Sharp (Reference Plohr and Sharp1988), Trangenstein & Colella (Reference Trangenstein and Colella1991) and Liu & Walkington (Reference Liu and Walkington2001) have been developed that use the deformation gradient tensor 
$\unicode[STIX]{x1D641}$ as a primary simulation variable. Takagi and coworkers have developed a related approach based on using the left Cauchy–Green deformation tensor 
$\unicode[STIX]{x1D63D}=\unicode[STIX]{x1D641}\unicode[STIX]{x1D641}^{\text{T}}$ as the primary simulation variable since this quantity features in many constitutive laws (Sugiyama et al. Reference Sugiyama, Ii, Takeuchi, Takagi and Matsumoto2011; Takagi et al. Reference Takagi, Sugiyama, Ii and Matsumoto2011; Ii et al. Reference Ii, Sugiyama, Takagi and Matsumoto2012). Methods based on 
$\unicode[STIX]{x1D641}$ and 
$\unicode[STIX]{x1D63D}$ have the advantage of requiring a single spatial derivative to compute acceleration, whereas our approach requires two, in (2.5) and (2.3). However, the reference map contains additional information for locating the boundary and has fewer components to store. Using 
$\unicode[STIX]{x1D641}$ requires that a gauge constraint is imposed to ensure that it remains a valid gradient of a material mapping.
Cottet, Maitre & Milcent (Reference Cottet, Maitre and Milcent2008) and Maitre et al. (Reference Maitre, Milcent, Cottet, Raoult and Usson2009) independently developed simulation approaches using the reference map, and stated its potential to be used as part of a general fluid–structure simulation. Further work by Milcent & Maitre (Reference Milcent and Maitre2016) simulated an immersed interface with full membrane elasticity. Bellotti & Theillard (Reference Bellotti and Theillard2019) coupled the reference map to the level set method for improved tracking of interfaces in two-phase flow simulation.
The IPS method of Dunne (Reference Dunne2006) is closely related to our approach. The IPS method is based on using finite elements to track the solid displacements 
$\boldsymbol{u}$ as the primary simulation variable. However, the reference map field (referred to by Dunne as the set of initial positions) emerges as part of the computational procedure, via the relationship 
$\boldsymbol{u}=\boldsymbol{x}-\unicode[STIX]{x1D743}$. The IPS method has been developed further (Dunne, Rannacher & Richter Reference Dunne, Rannacher and Richter2010; Richter Reference Richter2013; Wick Reference Wick2013), and has the advantage that the Eulerian field 
$\boldsymbol{u}$ can be used to track the deformation of the fluid–structure interface, in a similar manner to the approach we discuss in § 3.2.
2.3 Incompressible fluid–structure interaction
In this paper we employ the reference map technique to simulate incompressible fluid–structure interactions. We shall use the terms 
$\unicode[STIX]{x1D749}$, 
$\unicode[STIX]{x1D749}_{s}$ and 
$\unicode[STIX]{x1D749}_{f}$ 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 
$\boldsymbol{v}(\boldsymbol{x},t)$ that satisfies the incompressibility constraint
$$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0.\end{eqnarray}$$ We consider a solid immersed within the fluid, and introduce a level set function 
$\unicode[STIX]{x1D719}(\boldsymbol{x},t)$ (Sethian Reference Sethian1999; Osher & Fedkiw Reference Osher and Fedkiw2003) that is the signed distance to the solid–fluid interface with the convention that 
$\unicode[STIX]{x1D719}<0$ in the solid and 
$\unicode[STIX]{x1D719}>0$ in the fluid. The reference map 
$\unicode[STIX]{x1D743}(\boldsymbol{x},t)$ is defined within the solid region only.
Let the fluid have density 
$\unicode[STIX]{x1D70C}_{f}$ and dynamic viscosity 
$\unicode[STIX]{x1D707}_{f}$. The fluid stress deviator at any gridpoint is given by
$$\begin{eqnarray}\unicode[STIX]{x1D749}_{f}=\unicode[STIX]{x1D707}_{f}(\unicode[STIX]{x1D735}\boldsymbol{v}+(\unicode[STIX]{x1D735}\boldsymbol{v})^{\text{T}}).\end{eqnarray}$$ Kinematic viscosity is defined as 
$\unicode[STIX]{x1D708}_{f}=\unicode[STIX]{x1D707}_{f}/\unicode[STIX]{x1D70C}_{f}$. The deviatoric stress is then defined as a smooth transition between the fluid and solid stresses,
$$\begin{eqnarray}\unicode[STIX]{x1D749}=\unicode[STIX]{x1D749}_{s}+H_{\unicode[STIX]{x1D700}}(\unicode[STIX]{x1D719})(\unicode[STIX]{x1D749}_{f}-\unicode[STIX]{x1D749}_{s}),\end{eqnarray}$$where
$$\begin{eqnarray}H_{\unicode[STIX]{x1D700}}(\unicode[STIX]{x1D719})=\left\{\begin{array}{@{}ll@{}}0 & \text{if }\unicode[STIX]{x1D719}\leqslant -\unicode[STIX]{x1D700},\\ \displaystyle \frac{1}{2}\left(1+\frac{\unicode[STIX]{x1D719}}{\unicode[STIX]{x1D700}}+\frac{1}{\unicode[STIX]{x03C0}}\sin \frac{\unicode[STIX]{x03C0}\unicode[STIX]{x1D719}}{\unicode[STIX]{x1D700}}\right) & \text{if }|\unicode[STIX]{x1D719}|<\unicode[STIX]{x1D700},\\ 1 & \text{if }\unicode[STIX]{x1D719}\geqslant \unicode[STIX]{x1D700},\end{array}\right.\end{eqnarray}$$ is a smoothed Heaviside function with a transition region of width 
$2\unicode[STIX]{x1D700}$. Here we use a twice-differentiable form for 
$H_{\unicode[STIX]{x1D700}}$ that has been employed in previous studies (Sussman, Smereka & Osher Reference Sussman, Smereka and Osher1994; Sussman et al. Reference Sussman, Almgren, Bell, Colella, Howell and Welcome1999; Yu et al. Reference Yu, Sakai and Sethian2003, Reference Yu, Sakai and Sethian2007). For more details on the choice of 
$\unicode[STIX]{x1D700}$ and the precise form of 
$H_{\unicode[STIX]{x1D700}}$, see appendix D. In order to calculate 
$\unicode[STIX]{x1D749}$ it is necessary to smoothly extend 
$\unicode[STIX]{x1D743}$ in the region 
$0<\unicode[STIX]{x1D719}<\unicode[STIX]{x1D700}$, 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
$$\begin{eqnarray}\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{s}+H_{\unicode[STIX]{x1D700}}(\unicode[STIX]{x1D719})(\unicode[STIX]{x1D70C}_{f}-\unicode[STIX]{x1D70C}_{s}).\end{eqnarray}$$ Formally, our limit of interest is when 
$\unicode[STIX]{x1D700}\rightarrow 0$, when there is a sharp interface between the fluid and solid. In our numerical method, we choose 
$\unicode[STIX]{x1D700}$ to scale proportionally with the grid spacing, and thus we approach this limit as the simulation resolution is increased.
3 Numerical method
The numerical procedure is based on the projection method of Chorin (Reference Chorin1967, Reference Chorin1968) for solving the incompressible Navier–Stokes equations. Specifically, we consider a modern second-order method described by Almgren et al. (Reference Almgren, Bell, Colella, Howell and Welcome1998) 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 
$W$ by 
$H$, and is divided into an 
$M\times N$ grid of rectangular cells of size 
$h_{x}$ by 
$h_{y}$. Following the work of Colella (Reference Colella1990), the velocity, the reference map and the level set are held at cell centres and are indexed as 
$\boldsymbol{v}_{i,j}$, 
$\unicode[STIX]{x1D743}_{i,j}$ and 
$\unicode[STIX]{x1D719}_{i,j}$, respectively, for 
$i=0,\ldots ,M-1$ and 
$j=0,\ldots ,N-1$ (figure 2a). The components of the velocity field are written as 
$\boldsymbol{v}=(u,v)$. Pressures are held at cell corners and are indexed as 
$p_{i,j}$ for 
$i=0,\ldots ,M$ and 
$j=0,\ldots ,N$. 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 
$\unicode[STIX]{x1D743}_{i,j}$, velocity 
$\boldsymbol{v}_{i,j}$ and level set field 
$\unicode[STIX]{x1D719}_{i,j}$ are held at the cell centre, 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-time step to evaluate the advective terms.
Superscripts are used to denote time steps. To advance the simulation forward from time step 
$n$ to 
$n+1$ with interval 
$\unicode[STIX]{x0394}t$, the following procedure is used. The reference map field is updated using
$$\begin{eqnarray}\frac{\unicode[STIX]{x1D743}^{n+1}-\unicode[STIX]{x1D743}^{n}}{\unicode[STIX]{x0394}t}=-[(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\unicode[STIX]{x1D743}]^{n+1/2}\end{eqnarray}$$ and an intermediate velocity 
$\boldsymbol{v}^{\ast }$ is computed using
$$\begin{eqnarray}\frac{\boldsymbol{v}^{\ast }-\boldsymbol{v}^{n}}{\unicode[STIX]{x0394}t}=-[(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{v}]^{n+1/2}+\frac{1}{\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D719}^{n+1/2})}\unicode[STIX]{x1D735}\boldsymbol{\cdot }[\unicode[STIX]{x1D749}(\unicode[STIX]{x1D743}^{n+1/2},\boldsymbol{v}^{n})].\end{eqnarray}$$ Here, the advective derivatives 
$[(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\unicode[STIX]{x1D743}]^{n+1/2}$ and 
$[(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{v}]^{n+1/2}$ are evaluated at the middle of the time step using a second-order explicit Godunov scheme, described in § 3.1. Once the advective derivatives are evaluated, equation (3.1) allows 
$\unicode[STIX]{x1D743}^{n+1}$ to be computed. This allows the time-centred reference map to be defined as 
$\unicode[STIX]{x1D743}^{n+1/2}=(\unicode[STIX]{x1D743}^{n}+\unicode[STIX]{x1D743}^{n+1})/2$ after which 
$\boldsymbol{v}^{\ast }$ is computed using (3.2). From here, the Poisson problem for pressure is evaluated using
$$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}^{\ast }=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\frac{\unicode[STIX]{x0394}t}{\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D719}^{n+1/2})}\unicode[STIX]{x1D735}p^{n+1}\right).\end{eqnarray}$$Following Almgren et al. (Reference Almgren, Bell and Szymczak1996) and Puckett et al. (Reference Puckett, Almgren, Bell, Marcus and Rider1997), equation (3.3) is solved using a finite-element formulation, described in § 3.4. After this, the velocity is projected to be divergence free using
$$\begin{eqnarray}\boldsymbol{v}^{n+1}=\boldsymbol{v}^{\ast }-\frac{\unicode[STIX]{x0394}t}{\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D719}^{n+1/2})}\unicode[STIX]{x1D735}p^{n+1},\end{eqnarray}$$ where the gradient of 
$p^{n+1}$ is evaluated using a second-order centred difference formula.
3.1 Advective terms
To evaluate the advective terms appearing in (3.1) and (3.2), a second-order explicit Godunov scheme is used. The same scheme is applied to both the velocity 
$\boldsymbol{v}$ and reference map 
$\unicode[STIX]{x1D743}$. Throughout this section, we denote 
$a$ to be a generic scalar component of either of these two fields. We also refer the reader to recent work by Jain & Mani (Reference Jain and Mani2017), 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 centre are evaluated using the fourth-order monotonicity-limited scheme of Colella (Reference Colella1985) described in § A.1. Once the gradients are calculated at the centre of each cell, edge-centred velocities and reference maps are created at 
$t+\unicode[STIX]{x0394}t/2$ using Taylor expansions to each of the four edges, which are indexed using half-integers as shown in figure 2(b). As an example, an extrapolation of the reference map to a vertical edge from the left (with superscript 
$L$) is given by
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D743}_{i+1/2,j}^{L,n+1/2} & = & \displaystyle \unicode[STIX]{x1D743}_{i,j}^{n}+\frac{\unicode[STIX]{x0394}t}{2}(\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D743})_{i,j}^{n}+\frac{h_{x}}{2}(\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D743})_{i,j}^{n}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D743}_{i,j}+\frac{1}{2}(h_{x}-u_{i,j}^{n}\unicode[STIX]{x0394}t)\unicode[STIX]{x1D743}_{x,i,j}^{n}-\frac{\unicode[STIX]{x0394}t}{2}(\widetilde{v\unicode[STIX]{x1D743}_{y}})_{i,j}^{n},\end{eqnarray}$$ where (2.4) has been substituted for the 
$\unicode[STIX]{x1D743}_{t}$ derivative. The extrapolation of the velocity from the left to this edge is
$$\begin{eqnarray}\displaystyle \boldsymbol{v}_{i+1/2,j}^{L,n+1/2} & = & \displaystyle \boldsymbol{v}_{i,j}^{n}+\frac{\unicode[STIX]{x0394}t}{2}\boldsymbol{v}_{t,i,j}^{n}+\frac{h_{x}}{2}\boldsymbol{v}_{x,i,j}^{n}\nonumber\\ \displaystyle & = & \displaystyle \boldsymbol{v}_{i,j}^{n}+\frac{1}{2}(h_{x}-u_{i,j}^{n}\unicode[STIX]{x0394}t)\boldsymbol{v}_{x,i,j}^{n}-\frac{\unicode[STIX]{x0394}t}{2}(\widetilde{v\boldsymbol{v}_{y}})_{i,j}^{n}-\frac{\unicode[STIX]{x0394}t}{2}\boldsymbol{a}_{i,j}^{n},\end{eqnarray}$$ where (2.3) has been substituted to replace the 
$\boldsymbol{v}_{t,i,j}^{n}$ term, and 
$\boldsymbol{a}_{i,j}^{n}$ is defined according to
$$\begin{eqnarray}\boldsymbol{a}_{i,j}^{n}=\left[-\frac{1}{\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D719})}\unicode[STIX]{x1D735}p+\frac{1}{\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D719})}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\right]_{i,j}^{n}.\end{eqnarray}$$ Equivalent procedures are used to compute extrapolations from the right, down and up with superscripts 
$R$, 
$D$ and 
$U$, respectively. To ensure tangential stability the terms with tildes in (3.5) and (3.6) are computed differently using the procedure in § 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 
$(i+1/2,j)$,
$$\begin{eqnarray}a_{i+1/2,j}^{n+1/2}=\left\{\begin{array}{@{}ll@{}}a_{i+1/2,j}^{L,n+1/2}\quad & \text{if }u_{i+1/2,j}^{L,n+1/2}>0\text{ and }u_{i+1/2,j}^{L,n+1/2}+u_{i+1/2,j}^{R,n+1/2}>0,\\ a_{i+1/2,j}^{R,n+1/2}\quad & \text{if }u_{i+1/2,j}^{R,n+1/2}<0\text{ and }u_{i+1/2,j}^{L,n+1/2}+u_{i+1/2,j}^{R,n+1/2}<0,\\ {\mathcal{F}}(a_{i+1/2,j}^{L,n+1/2},a_{i+1/2,j}^{R,n+1/2})\quad & \text{otherwise,}\end{array}\right.\end{eqnarray}$$ where 
$a$ 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 
${\mathcal{F}}$ is used when the two velocities are ambiguous. For the horizontal velocity 
${\mathcal{F}}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D6FE})=0$ (Case A), and for all other components 
${\mathcal{F}}(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D6FE})=(\unicode[STIX]{x1D6FD}+\unicode[STIX]{x1D6FE})/2$ (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 § 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 
$\boldsymbol{v}_{e}$ be the edge velocities, and let 
$q$ be a cell-centred scalar field. We aim to make
$$\begin{eqnarray}\boldsymbol{v}_{e}-\frac{1}{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D735}q\end{eqnarray}$$divergence free. Taking the divergence of (3.9) yields
$$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\frac{1}{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D735}q\right)=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}_{e},\end{eqnarray}$$which is discretized as
$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{1}{h_{x}^{2}}\left(\frac{q_{i+1,j}-q_{i,j}}{\unicode[STIX]{x1D70C}_{i+1/2,j}}-\frac{q_{i,j}-q_{i-1,j}}{\unicode[STIX]{x1D70C}_{i-1/2,j}}\right)+\frac{1}{h_{y}^{2}}\left(\frac{q_{i,j+1}-q_{i,j}}{\unicode[STIX]{x1D70C}_{i,j+1/2}}-\frac{q_{i,j}-q_{i,j-1}}{\unicode[STIX]{x1D70C}_{i,j-1/2}}\right)\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{u_{i+1/2,j}-u_{i-1/2,j}}{h_{x}}+\frac{v_{i,j+1/2}-v_{i,j-1/2}}{h_{y}}.\end{eqnarray}$$ 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 (3.11) is omitted if it contains 
$q$ values that are out of range. If a pressure condition is applied, then a Dirichlet condition of 
$q=\unicode[STIX]{x0394}t\,p/2$ is applied, where the factor of two arises because the edge velocities are time centred.
Equation (3.11) results in a large linear system 
$\unicode[STIX]{x1D63C}q=b$ where 
$\unicode[STIX]{x1D63C}$ is a sparse matrix, 
$b$ is the source term and 
$q$ is a vector of the components 
$q_{i,j}$, which we solve using the geometric multigrid method with a standard V-cycle iteration (Demmel Reference Demmel1997). Since the 
$q$ field typically varies smoothly in time, the initial guess for the multigrid method is computed as a linear interpolation from the previous two time steps. Multigrid V-cycles are performed until the root mean squared (r.m.s.) element in the residual vector 
$r=\unicode[STIX]{x1D63C}q-b$ reaches a required tolerance 
$T_{MAC}$. 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 
$r_{s}=4(h_{x}^{-2}+h_{y}^{-2})\unicode[STIX]{x0394}t$, and a tolerance of 
$T_{MAC}=10^{4}r_{s}\unicode[STIX]{x1D700}_{m}$ is used, where 
$\unicode[STIX]{x1D700}_{m}$ 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. See appendices C and F for more details on the multigrid code library and performance.
3.1.3 Evaluation of the derivative
Once the MAC projection has been performed the time-centred advective terms for the velocity and reference maps are evaluated as
$$\begin{eqnarray}\displaystyle [(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735})a]_{i,j}^{n+1/2} & = & \displaystyle \frac{u_{i+1/2,j}^{n+1/2}+u_{i-1/2,j}^{n+1/2}}{2}\frac{a_{i+1/2,j}^{n+1/2}-a_{i-1/2,j}^{n+1/2}}{h_{x}}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{v_{i,j+1/2}^{n+1/2}+v_{i,j-1/2}^{n+1/2}}{2}\frac{a_{i,j+1/2}^{n+1/2}-a_{i,j-1/2}^{n+1/2}}{h_{y}},\end{eqnarray}$$ where 
$a$ is a generic field component.
3.2 Level set update and reference map extrapolation
The simulation makes use of a cell-centred level set function 
$\unicode[STIX]{x1D719}_{i,j}$ for tracking the fluid–solid boundary. The level set routine is stored in a narrow band of points of width 
$2\unicode[STIX]{x1D719}_{W}$ surrounding the interface (Sethian Reference Sethian1999; Rycroft & Gibou Reference Rycroft and Gibou2012) 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 (2.8) and (2.10), respectively. Unlike typical applications of the level set method, the 
$\unicode[STIX]{x1D719}$ field is not explicitly updated, but is instead continually given by the reference map field using the procedure first described by Valkov et al. (Reference Valkov, Rycroft and Kamrin2015).
3.2.1 Level set construction
For a given shape, define a continuous function of the reference map 
$\unicode[STIX]{x1D719}_{0}(\unicode[STIX]{x1D743})$ such that 
$\unicode[STIX]{x1D719}_{0}<0$ for reference map values in the solid, 
$\unicode[STIX]{x1D719}_{0}>0$ for reference map values outside the solid and 
$\unicode[STIX]{x1D719}_{0}=0$ on the interface. During the time step, the reference map field 
$\unicode[STIX]{x1D743}^{n+1}$ is computed inside the solid using (3.1), from which the half-time-step reference map is defined as 
$\unicode[STIX]{x1D743}^{n+1/2}=(\unicode[STIX]{x1D743}^{n}+\unicode[STIX]{x1D743}^{n+1})/2$. Both fields are extended into the narrow band fluid region using the extrapolation methods described in § 3.2.2.
To construct the new level set function 
$\unicode[STIX]{x1D719}^{n+1/2}$, an auxiliary function 
$\unicode[STIX]{x1D713}^{n+1/2}$ is first computed in the narrow band such that 
$\unicode[STIX]{x1D713}^{n+1/2}=\unicode[STIX]{x1D719}_{0}(\unicode[STIX]{x1D743}^{n+1/2})$. The zero contour of 
$\unicode[STIX]{x1D713}^{n+1/2}$ 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 
$\unicode[STIX]{x1D719}^{n+1/2}$ is constructed from 
$\unicode[STIX]{x1D713}^{n+1/2}$ using the reinitialization procedure described by Rycroft & Gibou (Reference Rycroft and Gibou2012). This procedure first considers points 
$(i,j)$ that straddle the interface, so that one of their orthogonal neighbours has a 
$\unicode[STIX]{x1D713}^{n+1/2}$ value of an opposite sign. Each straddling point is considered. The bicubic interpolant 
$\unicode[STIX]{x1D713}_{c}^{n+1/2}$ is computed, and the modified Newton–Raphson approach of Chopp (Reference Chopp2001, Reference Chopp2009) is used to find the shortest distance vector 
$\unicode[STIX]{x1D71F}\boldsymbol{x}$ from each straddling point to the interface 
$\unicode[STIX]{x1D713}_{c}^{n+1/2}(\boldsymbol{x})=0$, after which the level set function is initialized to 
$\pm |\unicode[STIX]{x1D71F}\boldsymbol{x}|$. 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 (Reference Sethian1999). For further details, see the paper by Rycroft & Gibou (Reference Rycroft and Gibou2012).
With the straddling points of 
$\unicode[STIX]{x1D719}^{n+1/2}$ now initialized, the remaining points are filled in using the second-order fast marching method of Sethian (Reference Sethian1999). In the fluid, the positive 
$\unicode[STIX]{x1D719}^{n+1/2}$ values are computed in order of increasing value, until reaching a cutoff 
$\unicode[STIX]{x1D719}_{W}$ that defines the width of the narrow band. The same procedure is used to fill in the negative 
$\unicode[STIX]{x1D719}^{n+1/2}$ values in the solid, until reaching a cutoff 
$-\unicode[STIX]{x1D719}_{W}$. 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 
$\unicode[STIX]{x1D713}_{n+1/2}$ has a loss of regularity as some points: since the entire 
$\unicode[STIX]{x1D719}^{n+1/2}$ field is directly constructed, there is no possibility for instabilities to grow over time, as can happen in update procedures that use partial differential equations (PDEs). Identical methods are used to construct 
$\unicode[STIX]{x1D719}^{n+1}$ from 
$\unicode[STIX]{x1D743}^{n+1}$.
3.2.2 Extrapolation
During the construction of the level set function, a list of fluid grid points sorted in order of increasing value, 
$0<\unicode[STIX]{x1D719}_{1}<\unicode[STIX]{x1D719}_{2}<\cdots \,$ is constructed, which is used for extrapolating the reference map 
$\unicode[STIX]{x1D743}$ from the solid into the fluid narrow band. Previous approaches to do this have employed PDE-based methods by defining a normal vector 
$\boldsymbol{n}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}$ and extrapolating outwards from the object in the direction of 
$\boldsymbol{n}$ (Aslam Reference Aslam2004; Rycroft & Gibou Reference Rycroft and Gibou2012). While these methods are well suited to mathematical analysis, they require considerable bookkeeping for performing the finite-difference calculations of 
$\unicode[STIX]{x1D719}$ and 
$\unicode[STIX]{x1D743}$ 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. Reference Valkov, Rycroft and Kamrin2015).
In the current work, we make use of the following alternative extrapolation procedure. Consider the points in increasing order of 
$\unicode[STIX]{x1D719}$ value. For a particular point 
$(i,j)$ at physical location 
$\boldsymbol{x}_{i,j}$:
(i) Set
$r=3$.(ii) Use least-squares regression to fit a linear map
$\unicode[STIX]{x1D743}_{lm}(\boldsymbol{x})=\boldsymbol{A}x+\boldsymbol{B}y+\boldsymbol{C}$ using all available reference map values at $(i^{\prime },j^{\prime })$ such that $|i-i^{\prime }|\leqslant r,|j-j^{\prime }|\leqslant r$. Weight each value in the regression according to $\unicode[STIX]{x1D719}_{i,j}-\unicode[STIX]{x1D719}_{i^{\prime },j^{\prime }}$.(iii) If the linear map is degenerate then increment
$r$ and return to Step 2. Otherwise, continue.(iv) Set
$\unicode[STIX]{x1D743}_{i,j}=\unicode[STIX]{x1D743}_{lm}(\boldsymbol{x}_{i,j})$.
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 neighbourhood, 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 neighbouring points.
The approach described here makes it possible to simulate objects with sharp corners. The reference map is smoothly defined within the object, and using the above procedure allows it to be smoothly extended into the fluid domain from which solid stresses in the blur zone can be computed. The function 
$\unicode[STIX]{x1D719}_{0}(\unicode[STIX]{x1D743})$ that defines the object boundary need not be smooth itself, and can describe a shape with corners. (The IPS method of Dunne (Reference Dunne2006) has a similar capability to handle corners, although it is based on a harmonic continuation of the solid velocity, which is used to update the solid displacement field 
$\boldsymbol{u}$, from which the boundary can be located by examining 
$\boldsymbol{u}-\boldsymbol{x}$.)
3.3 Computation of stress
In order to evaluate the stress divergence terms that appear in (3.2) and (3.7), the stresses are first computed on the edges of each grid cell. The stress term in (3.7) is computed as
$$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }[\unicode[STIX]{x1D749}(\unicode[STIX]{x1D743}^{n})]=\frac{[\unicode[STIX]{x1D749}_{x}]_{i+1/2,j}^{n}-[\unicode[STIX]{x1D749}_{x}]_{i-1/2,j}^{n}}{h_{x}}+\frac{[\unicode[STIX]{x1D749}_{y}]_{i,j+1/2}^{n}-[\unicode[STIX]{x1D749}_{y}]_{i,j-1/2}^{n}}{h_{y}},\end{eqnarray}$$ where 
$\unicode[STIX]{x1D749}_{x}=(\unicode[STIX]{x1D70F}_{xx},\unicode[STIX]{x1D70F}_{xy})$ and 
$\unicode[STIX]{x1D749}_{y}=(\unicode[STIX]{x1D70F}_{xy},\unicode[STIX]{x1D70F}_{yy})$ 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
$$\begin{eqnarray}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D743}}{\unicode[STIX]{x2202}x}\right)_{i-1/2,j}=\frac{\unicode[STIX]{x1D743}_{i,j}-\unicode[STIX]{x1D743}_{i-1,j}}{h_{x}},\quad \left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D743}}{\unicode[STIX]{x2202}y}\right)_{i-1/2,j}=\frac{\unicode[STIX]{x1D743}_{i,j+1}+\unicode[STIX]{x1D743}_{i-1,j+1}-\unicode[STIX]{x1D743}_{i,j-1}-\unicode[STIX]{x1D743}_{i-1,j-1}}{4h_{y}}\end{eqnarray}$$after which the deformation gradient is evaluated as
$$\begin{eqnarray}\unicode[STIX]{x1D641}_{i-1/2,j}=\left(\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D743}}{\unicode[STIX]{x2202}\boldsymbol{x}}\right)_{i-1/2,j}\right)^{-1}.\end{eqnarray}$$ From here, any constitutive law 
$\unicode[STIX]{x1D749}_{s}=\boldsymbol{f}(\unicode[STIX]{x1D641})$ could be used to evaluate the deviatoric stress, 
$\unicode[STIX]{x1D749}_{s}$. In the current work, we employ the plane-strain incompressible neo-Hookean law,
$$\begin{eqnarray}\unicode[STIX]{x1D749}_{s}=\boldsymbol{f}(\unicode[STIX]{x1D641})=G(\unicode[STIX]{x1D641}\unicode[STIX]{x1D641}^{\text{T}}-{\textstyle \frac{1}{3}}\mathbf{1}(\operatorname{tr}\unicode[STIX]{x1D641}\unicode[STIX]{x1D641}^{\text{T}}+1)),\end{eqnarray}$$ where 
$G$ 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
$$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\unicode[STIX]{x2202}\boldsymbol{v}}{\unicode[STIX]{x2202}x}\right)_{i-1/2,j}=\frac{\boldsymbol{v}_{i,j}-\boldsymbol{v}_{i-1,j}}{h_{x}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\unicode[STIX]{x2202}\boldsymbol{v}}{\unicode[STIX]{x2202}y}\right)_{i-1/2,j}=\frac{\boldsymbol{v}_{i,j+1}+\boldsymbol{v}_{i-1,j+1}-\boldsymbol{v}_{i,j-1}-\boldsymbol{v}_{i-1,j-1}}{4h_{y}}. & \displaystyle\end{eqnarray}$$Equivalent stencils are used to compute velocity gradients on horizontal edges, after which the fluid stress is given by
$$\begin{eqnarray}\unicode[STIX]{x1D749}_{f}=\unicode[STIX]{x1D707}_{f}(\unicode[STIX]{x1D735}\boldsymbol{v}+(\unicode[STIX]{x1D735}\boldsymbol{v})^{\text{T}}),\end{eqnarray}$$ where 
$\unicode[STIX]{x1D707}_{f}$ is the viscosity. Equation (3.19) is our standard approach for computing the fluid stress. However, we have also investigated a simplified calculation. Since
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0$, it follows that in the bulk of the fluid, the second term in (3.19) has zero contribution to 
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}_{f}$. Hence an alternative formula is
$$\begin{eqnarray}\unicode[STIX]{x1D749}_{f}^{(simp)}=\unicode[STIX]{x1D707}_{f}\unicode[STIX]{x1D735}\boldsymbol{v}.\end{eqnarray}$$This formula only requires evaluating the simpler stencil in (3.17). However, equation (3.20) is not strictly valid in the blur zone since taking the divergence in (2.8) results in a non-zero contribution from the second term of (3.19).
Once all edge stresses are computed, the divergence is computed using
$$\begin{eqnarray}[\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}]_{i,j}=\frac{\unicode[STIX]{x1D749}_{i+1/2,j}-\unicode[STIX]{x1D749}_{i-1/2,j}}{h_{x}}+\frac{\unicode[STIX]{x1D749}_{i,j+1/2}-\unicode[STIX]{x1D749}_{i,j-1/2}}{h_{y}}.\end{eqnarray}$$3.4 Finite-element projection
To solve the Poisson problem in (3.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 
$\unicode[STIX]{x1D713}$ the weak formulation of (3.3) is
$$\begin{eqnarray}-\int _{\unicode[STIX]{x1D6FA}}\boldsymbol{v}^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\,\text{d}x\,\text{d}y+\int _{\unicode[STIX]{x1D6FA}}\frac{\unicode[STIX]{x0394}t}{\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D719}^{n+1/2})}\unicode[STIX]{x1D735}p^{n+1}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\,\text{d}x\,\text{d}y=-\int _{\unicode[STIX]{x1D6E4}_{1}}\unicode[STIX]{x1D713}\boldsymbol{v}^{BC}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}S,\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6E4}_{1}$ is the section of the boundary where inflow and outflow conditions are prescribed. Consider a particular bilinear element function 
$\unicode[STIX]{x1D713}$ located at a pressure point 
$p_{i,j}$ in the bulk of the domain. The first term of (3.22) is
$$\begin{eqnarray}\frac{h_{y}(u_{i,j}^{\ast }+u_{i,j-1}^{\ast }-u_{i-1,j}^{\ast }-u_{i-1,j-1}^{\ast })}{2}+\frac{h_{x}(v_{i,j}^{\ast }-v_{i,j-1}^{\ast }+v_{i-1,j}^{\ast }-v_{i-1,j-1}^{\ast })}{2}.\end{eqnarray}$$ For the constant-density case when 
$\unicode[STIX]{x0394}t/\unicode[STIX]{x1D70C}$ can be taken out as a prefactor, the second term of (3.22) is
$$\begin{eqnarray}\unicode[STIX]{x1D706}_{a}p_{i,j}+\unicode[STIX]{x1D706}_{b}(p_{i-1,j}+p_{i+1,j})+\unicode[STIX]{x1D706}_{c}(p_{i,j-1}+p_{i,j+1})+\unicode[STIX]{x1D706}_{d}\mathop{\sum }_{\substack{ k=\pm 1 \\ l=\pm 1}}p_{i+k,j+l},\end{eqnarray}$$where
$$\begin{eqnarray}\unicode[STIX]{x1D706}_{a}=\frac{4(h_{x}^{2}+h_{y}^{2})}{3h_{x}h_{y}},\quad \unicode[STIX]{x1D706}_{b}=\frac{-2h_{y}^{2}+h_{x}^{2}}{3h_{x}h_{y}},\quad \unicode[STIX]{x1D706}_{c}=\frac{-2h_{x}^{2}+h_{y}^{2}}{3h_{x}h_{y}},\quad \unicode[STIX]{x1D706}_{d}=\frac{-h_{x}^{2}-h_{y}^{2}}{6h_{x}h_{y}}.\end{eqnarray}$$ The expression in (3.24) can be generalized for the case of non-constant density. Dirichlet conditions on pressure may also be imposed as essential boundary conditions (Johnson Reference Johnson2009). The resultant linear system is solved using the same multigrid library introduced in § 3.1.2 using an error tolerance of 
$T_{FEM}=\unicode[STIX]{x1D706}_{a}10^{4}\unicode[STIX]{x1D700}_{m}$. In cases where no Dirichlet conditions are used, the pressure field is projected at each step so that it has zero mean.
3.5 Parameter choices and stability
Our test cases involve four physical parameters: solid shear modulus 
$G$, solid density 
$\unicode[STIX]{x1D70C}_{s}$, fluid viscosity 
$\unicode[STIX]{x1D707}_{f}$ and fluid density 
$\unicode[STIX]{x1D70C}_{f}$. In the solid, the shear wave speed is
$c_{s}=\sqrt{G/\unicode[STIX]{x1D70C}_{s}}$. The CFL condition requires that the simulation time step be less than or equal to
$$\begin{eqnarray}\unicode[STIX]{x0394}t_{\text{I}}=c_{s}\min \{h_{x},h_{y}\}=\sqrt{\frac{G}{\unicode[STIX]{x1D70C}_{s}}}\min \{h_{x},h_{y}\}.\end{eqnarray}$$In addition, performing a von Neumann stability analysis shows that the time step must be less than or equal to
$$\begin{eqnarray}\unicode[STIX]{x0394}t_{\text{II}}=\frac{\unicode[STIX]{x1D70C}_{f}}{2\unicode[STIX]{x1D707}_{f}(h_{x}^{-2}+h_{y}^{-2})}\end{eqnarray}$$in order to resolve the viscous fluid stress. Inside the solid, we find that simulating stress using only equation (3.16) results in an instability – this should be expected since we are effectively solving a hyperbolic system using centred 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
$$\begin{eqnarray}\unicode[STIX]{x1D707}_{e}=\unicode[STIX]{x1D705}_{e}\unicode[STIX]{x1D70C}_{s}c_{s}\max \{h_{x},h_{y}\},\end{eqnarray}$$ where 
$\unicode[STIX]{x1D705}_{e}$ 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
$$\begin{eqnarray}\unicode[STIX]{x1D749}_{e}(\boldsymbol{x})=\unicode[STIX]{x1D707}_{e}[(1-H(\unicode[STIX]{x1D719}(\boldsymbol{x})))+q(1-\unicode[STIX]{x1D700}H^{\prime }(\unicode[STIX]{x1D719}(\boldsymbol{x})))]\unicode[STIX]{x1D735}\boldsymbol{v},\end{eqnarray}$$ where 
$q$ is a dimensionless constant. Based on a variety of tests, we use 
$q=1$ and 
$\unicode[STIX]{x1D705}_{e}=0.4$ 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 (3.20). Since 
$\unicode[STIX]{x1D707}_{e}$ 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 time-step restriction is
$$\begin{eqnarray}\unicode[STIX]{x0394}t_{\text{III}}=\frac{\unicode[STIX]{x1D70C}_{s}}{2\unicode[STIX]{x1D707}_{e}(1+q)(h_{x}^{-2}+h_{y}^{-2})}.\end{eqnarray}$$ With these definitions in place, the simulation time step 
$\unicode[STIX]{x0394}t$ is chosen to be smaller than the minimum of the three conditions in (3.26), (3.27) and (3.30), so that
$$\begin{eqnarray}\unicode[STIX]{x0394}t=\min \{\unicode[STIX]{x1D6FC}_{pad}\unicode[STIX]{x0394}t_{\text{I}},\unicode[STIX]{x1D6FC}_{pad}\unicode[STIX]{x0394}t_{\text{II}},\unicode[STIX]{x1D6FD}_{pad}\unicode[STIX]{x0394}t_{\text{III}}\}.\end{eqnarray}$$ Here, 
$\unicode[STIX]{x1D6FC}_{pad}$ and 
$\unicode[STIX]{x1D6FD}_{pad}$ are padding factors that are smaller than one. For this paper we use 
$\unicode[STIX]{x1D6FC}_{pad}=0.4$ and 
$\unicode[STIX]{x1D6FD}_{pad}=0.8$, so that the restrictions arising from the two physical stresses (I and II) are applied more stringently than the one for the artificial stress (III). Note that in the limit as 
$h_{x},h_{y}\rightarrow 0$, 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, 
$h_{x}=h_{y}=h$.
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 centred on the origin and has vertices at 
$(L\cos (2\unicode[STIX]{x03C0}k/7),L\sin (2\unicode[STIX]{x03C0}k/7))$ for 
$k\in \mathbb{Z}$, with length scale 
$L=0.62$, density 
$\unicode[STIX]{x1D70C}_{s}=3$ and shear modulus 
$G=24$. The resolution is 
$800\times 800$, the simulation domain is 
$[\!-1,1\!)^{2}$ and periodic boundary conditions are used. The fluid and rotor are initially stationary. The region 
$r=|\boldsymbol{x}|<0.16$ is used as a pivot. To excite the fluid, the pivot is rotated with an oscillatory motion with angle 
$\unicode[STIX]{x1D703}(t)=\unicode[STIX]{x03C0}(1-\cos t)$. This is done by applying an external tether force to the pivot region of
$$\begin{eqnarray}\boldsymbol{f}_{teth}(\boldsymbol{x})=K_{teth}H_{\unicode[STIX]{x1D700}}(r-r_{teth})(\boldsymbol{x}-\unicode[STIX]{x1D64D}_{\unicode[STIX]{x1D703}(t)}\unicode[STIX]{x1D743}(\boldsymbol{x})),\end{eqnarray}$$ where 
$r_{teth}=0.16$ and 
$\unicode[STIX]{x1D64D}_{\unicode[STIX]{x1D703}(t)}$ is a rotation matrix with angle 
$\unicode[STIX]{x1D703}(t)$. The spring constant is set to 
$K_{teth}=10^{-2}\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x0394}t^{-2}$, which ensures that the natural frequency of the tether satisfies the time-step restriction imposed by the method. The fluid has density 
$\unicode[STIX]{x1D70C}_{f}=1$ and viscosity 
$\unicode[STIX]{x1D707}_{f}=10^{-3}$. The r.m.s. angular velocity of the rotor is 
$\unicode[STIX]{x1D714}_{rms}=\unicode[STIX]{x03C0}/2$, and hence a characteristic tip velocity is 
$\unicode[STIX]{x1D714}_{rms}L$. Hence, we define the Reynolds number for this simulation as
$$\begin{eqnarray}Re=\frac{\unicode[STIX]{x1D70C}_{f}L(\unicode[STIX]{x1D714}_{rms}L)}{\unicode[STIX]{x1D707}_{f}}\approx 600.\end{eqnarray}$$ The simulation was run from 
$t=0$ to 
$t=4\unicode[STIX]{x03C0}$ using 16 threads on a Linux computer with dual 10-core 2.2 GHz Intel Xeon E5-2630 processors. For the given parameters, the time step of 
$\unicode[STIX]{x0394}t=1.105\times 10^{-4}$ was determined by the extra viscous stress in the solid. Simulation output was saved at regular intervals of 
$\unicode[STIX]{x03C0}/150$. The total wall clock time for the simulation was 
$6.53$ h. A total of 114 000 time steps were performed, with each taking 206 ms to compute. A substantial fraction of the computation time is spent performing the two linear solves. The MAC projections take on average 13.75 V-cycles and require 43.1 ms per time step. The finite-element projections take on average 11.91 V-cycles and require 48.1 ms per time step.
Figure 3. Snapshots of vorticity 
$\unicode[STIX]{x1D714}$ in a simulation of a flexible seven-pointed rotor being spun with an oscillatory motion in 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. The dark blue dotted circle shows the pivot region. Simulation parameters are 
$(\unicode[STIX]{x1D70C}_{f},\unicode[STIX]{x1D707}_{f},\unicode[STIX]{x1D70C}_{s},G)=(1,10^{-3},3,24)$.
Snapshots of vorticity 
$\unicode[STIX]{x1D714}=\frac{1}{2}(\unicode[STIX]{x2202}_{x}v-\unicode[STIX]{x2202}_{y}u)$ in the simulation are shown in figure 3. The vorticity is computed on each grid cell corner, using centred finite differences of the velocities in the four adjoining grid cells. As the star begins to rotate, each point deforms clockwise, and vortices are shed from the points, which are visible at 
$t=4\unicode[STIX]{x03C0}/15$. By 
$t=\unicode[STIX]{x03C0}$, 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 
$t=2\unicode[STIX]{x03C0}$, the rotational symmetry of the fluid flow is lost, due to interactions across the periodic boundaries, which break the sevenfold symmetry. By 
$t=4\unicode[STIX]{x03C0}$, after two complete cycles of the oscillatory motion, there are vortices present throughout the fluid. Supplemental Movie 1 available at https://doi.org/10.1017/jfm.2020.353 shows the complete simulation. To visualize the fluid motion, the movie also shows a number of tracers with trajectories 
$\boldsymbol{x}(t)$. The tracers are initialized at random positions in the fluid and are updated using the ordinary differential equation
$$\begin{eqnarray}\frac{\text{d}\boldsymbol{x}}{\text{d}t}=\boldsymbol{v}_{bic}(\boldsymbol{x}(t)),\end{eqnarray}$$ where 
$\boldsymbol{v}_{bic}$ is the bicubic interpolation of the velocity field, and the time integration is performed using the second-order improved Euler method (Süli & Mayers Reference Süli and Mayers2003).
4.2 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. Reference Zhang, Childress, Libchaber and Shelley2000; Watanabe et al. Reference Watanabe, Suzuki, Sugihara and Sueoka2002; Zhu & Peskin Reference Zhu and Peskin2002; Connell & Yue Reference Connell and Yue2007). Following the problem description and notation of Connell & Yue (Reference Connell and Yue2007), we introduce a thin filament of length 
$L$, thickness 
$h\ll L$, density 
$\unicode[STIX]{x1D70C}_{s}$ and Young’s Modulus 
$E$, clamped at its left endpoint and submerged in fluid of kinematic viscosity 
$\unicode[STIX]{x1D708}_{f}$ and density 
$\unicode[STIX]{x1D70C}_{f}$, flowing rightward with speed 
$V$ at infinity. (For consistency with Connell & Yue (Reference Connell and Yue2007) we fully adhere to their notation. However, we draw attention to the reader that 
$h$ (filament thickness) has a different meaning than 
$h$ (grid spacing) used in all other sections. Furthermore, 
$\unicode[STIX]{x1D707}$ (mass ratio) is distinct from 
$\unicode[STIX]{x1D707}_{f}$ (dynamic viscosity).) Three dimensionless numbers can be introduced to study the dynamical behaviour of the filament: the mass ratio 
$\unicode[STIX]{x1D707}=\unicode[STIX]{x1D70C}_{s}h/\unicode[STIX]{x1D70C}_{f}L$, Reynolds number 
$Re=VL/\unicode[STIX]{x1D708}_{f}$ and non-dimensional bending rigidity 
$K_{B}=EI/(\unicode[STIX]{x1D70C}_{f}V^{2}L^{3})$. 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 
$h$ does not necessarily satisfy 
$h\ll L$.
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 
$L=1$, 
$h=0.05$, 
$K_{B}=0.001$, in a fluid of density 
$\unicode[STIX]{x1D70C}_{f}=1$. To set 
$K_{B}$, we use the fact that in the linear elastic regime 
$E=3G$, and the moment of inertia is 
$I=h^{3}/12$. We vary 
$\unicode[STIX]{x1D70C}_{s}$ and 
$\unicode[STIX]{x1D708}_{f}$ in order to test a range of 
$\unicode[STIX]{x1D707}$ and 
$Re$. The simulation domain is set to be a 
$[-1,5]\times [\!-1,1\!)$ rectangle with assigned rightward velocity of speed 
$V=1$ on the left boundary, vanishing pressure on the right boundary and periodic boundary conditions on the top and bottom boundaries. We use a 
$1824\times 608$ grid to represent the domain. The filament is modelled as a rectangle 
$0<x<1$, 
$-h/2<y<h/2$ with semicircular end caps. The filament is anchored at 
$(0,0)$ using the tethering methodology described in § 4.1, with 
$\unicode[STIX]{x1D703}(t)=0$ and 
$r_{teth}=h/4$ in this case. We track the filament tip by introducing a tracer 
$\boldsymbol{x}(t)$ that starts from 
$(1,0)$ and is integrated according to (4.3). To prevent integration errors building up over time, the position of the tracer is periodically reset to satisfy 
$\unicode[STIX]{x1D743}_{bic}(\boldsymbol{x}(t))=(1,0)$ using a Newton–Raphson root-finding method, where 
$\unicode[STIX]{x1D743}_{bic}$ 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 approximately 3–12 days using 4–6 threads. Simulations with smaller 
$Re$ generally take longer, since resolving the fluid viscosity requires a smaller time step.
To systematically evaluate the behaviour of the filament, we store the perpendicular tip deflection 
$y(t)$ over the interval 
$t\in [120,160]$. Since the typical filament oscillation period is approximately 1.7, the simulations correspond to almost one hundred complete oscillations, and hence the interval 
$[120,160]$ is sufficient for the oscillation amplitude to reach a steady state. The normalized Fourier transform is given by
$$\begin{eqnarray}{\tilde{y}}(k)=\frac{1}{40}\int _{120}^{160}\text{e}^{\text{i}kt}y(t)\,\text{d}t.\end{eqnarray}$$The maximum Fourier amplitude is given by
$$\begin{eqnarray}A=\max _{k\in [0.1,50]}|{\tilde{y}}(k)|.\end{eqnarray}$$ Typically the oscillation frequency is 
$2\unicode[STIX]{x03C0}/1.7\approx 4$, and the range of 
$k$ in (4.5) is chosen to broadly cover the possible values. If 
$A\approx 0$ the filament is in the stable (no-flapping) regime and otherwise the filament is flapping, with 
$A$ 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 
$T_{MAC}$ and 
$T_{FEM}$ 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 
$A$ were insensitive to this. Connell & Yue (Reference Connell and Yue2007) proposed an analytical formula for the stable-to-flapping transition line
$$\begin{eqnarray}\unicode[STIX]{x1D707}=\frac{1.3Re^{-1/2}+K_{B}4\unicode[STIX]{x03C0}^{2}}{1-0.65Re^{-1/2}2\unicode[STIX]{x03C0}-0.5K_{B}8\unicode[STIX]{x03C0}^{3}}.\end{eqnarray}$$ Connell & Yue (Reference Connell and Yue2007) validated this equation numerically using a direct fluid–filament coupling procedure, a procedure that itself was validated against experiments (Zhang et al. Reference Zhang, Childress, Libchaber and Shelley2000; Watanabe et al. Reference Watanabe, Suzuki, Sugihara and Sueoka2002). The above formula is obtained without consideration of certain effects, such as possible variations of tension along the flag and the presence of global lift forces attempting to realign the flag with the flow, although methods to include these phenomena exist in certain limits (Argentina & Mahadevan Reference Argentina and Mahadevan2005). In figure 4 we show the behaviour of 
$A$ 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 
$A$ goes non-zero and the analytical curve. When 
$Re\geqslant 1000$ the transition predicted by the simulation happens at a slightly higher 
$\unicode[STIX]{x1D707}$ than predicted by (4.6). 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 
$h$ may also play a role.
Figure 4. Plot showing the steady-state oscillation amplitude 
$A$ of a thin flag with aspect ratio 20 and bending rigidity 
$K_{B}=0.001$, as a function of the Reynolds number 
$Re$ and mass ratio 
$\unicode[STIX]{x1D707}$. The colours shown are based on a bilinear interpolation of a two-dimensional grid of simulations. The axis ticks show the sampled values of 
$Re$ and 
$\unicode[STIX]{x1D707}$, with more simulations being performed in parameter ranges of interest. The thin dotted lines are contours at spacings of 
$(n/50)^{2}$ for 
$n\in \mathbb{N}$. The thick solid line is the stable-to-flapping transition formula, equation (4.6) of Connell & Yue (Reference Connell and Yue2007).
The behavioural switch from stable to flapping is also quite evident in the long-time flow fields, shown in figure 5. Small values of 
$\unicode[STIX]{x1D707}$ and 
$K_{B}$ result in stable behaviour, characterized by a straight filament and fluid flow that is symmetric about the filament axis (figure 5a). Upon crossing the transition, periodic undulatory filament motions develop with a fluid vortex street shed from the filament (figure 5b). Increasing 
$Re$ and 
$\unicode[STIX]{x1D707}$ even further reveals a chaotic filamentary motion, which was also observed by Connell & Yue (Reference Connell and Yue2007) (figure 5c, Supplemental Movie 2). The chaotic regime coincides with a drop in 
$A$ shown in the top right of figure 4 because the tip deflection no longer has a clear single dominant oscillatory mode. Because the filament is modelled 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.
Figure 5. Simulations of a thin flexible flag anchored at 
$(0,0)$ in a fluid with mean velocity 
$\boldsymbol{v}=(1,0)$, at 
$t=160$. The flag has an aspect ratio of 20. Three simulations with different parameters are shown: (a) stable with 
$(\unicode[STIX]{x1D707},K_{B},Re)=(0.04,0.001,400)$, (b) limit-cycle flapping with 
$(\unicode[STIX]{x1D707},K_{B},Re)=(0.16,0.001,1400)$ and (c) chaotic flapping with 
$(\unicode[STIX]{x1D707},K_{B},Re)=(0.32,0.001,3000)$. 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 small dark blue circles show the anchored regions. The colours show vorticity, using the same scale as figure 3.
Figure 6. Plot showing the steady-state oscillation amplitude 
$A$ as a function of the Reynolds number 
$Re$ and mass ratio 
$\unicode[STIX]{x1D707}$, for (a) flags with 
$K_{B}=0.002$ and aspect ratio 10, and (b) flags with 
$K_{B}=0.004$ and aspect ratio 5. The colours shown are based on a bilinear interpolation of a two-dimensional grid of simulations, using the same scale as figure 4. The axis ticks show the sampled values of 
$Re$ and 
$\unicode[STIX]{x1D707}$, with more simulations being performed in parameter ranges of interest. The thin dotted lines are contour at spacings of 
$(n/50)^{2}$ for 
$n\in \mathbb{N}$. The thick solid line in (a) is the stable-to-flapping transition formula for thin flags, equation (4.6) of Connell & Yue (Reference Connell and Yue2007). For (b), the formula is out of range and the entire parameter space is in the stable region.
Figure 7. Simulations of a thick flexible flag anchored at 
$(0,0)$ in a fluid with mean velocity 
$\boldsymbol{v}=(1,0)$, at 
$t=160$. The flag has an aspect ratio of 5. Two simulations with different parameters are shown: (a) vortex shedding with 
$(\unicode[STIX]{x1D707},K_{B},Re)=(0.04,0.004,750)$ and (b) limit-cycle flapping with 
$(\unicode[STIX]{x1D707},K_{B},Re)=(0.28,0.004,750)$. 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 dark blue dotted circles show the anchored regions. The colours show vorticity, using the same scale as figure 3.
We explore the importance of aspect ratio by introducing 
$R=h/L$ as an independent dimensionless group. We observe that as one departs from the 
$R\ll 1$ regime, adherence to (4.6) is diminished. In figure 6 we show results for 
$R^{-1}=10$ and 5. 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 
$R$ 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 
$D$ undergoes a transition from a laminar flow to a periodic vortex street as 
$DV/\unicode[STIX]{x1D708}_{f}$ grows beyond 
${\sim}50$ (Lienhard Reference Lienhard1966). In our case, the flag thickness acts like 
$D$, 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 
$Vh/\unicode[STIX]{x1D708}_{f}=Re\times R>50$ we see only flapping states for any choice of 
$\unicode[STIX]{x1D707}$ or 
$Re$. Figures 7(a) and 7(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.3 Tests for a range of elastic modulus
To demonstrate that the method works across a wide range of shear moduli, we consider a piston-like geometry where a flexible paddle is pushed through a fluid-filled cavity. The domain is 
$-1\leqslant x\leqslant 1$ and 
$0\leqslant y\leqslant 5$ and no-slip boundary conditions are used on all sides. The grid size is 
$160\times 400$, the fluid density is 
$\unicode[STIX]{x1D70C}_{f}=1$, the fluid viscosity is 
$\unicode[STIX]{x1D707}_{f}=10^{-3}$ and the solid density is 
$\unicode[STIX]{x1D70C}_{s}=2$.
Figure 8. Snapshots of pressure 
$p$ for a simulation where a flexible paddle is pushed through a fluid-filled cavity. The paddle has shear modulus 
$G=100$ and is anchored on its right end. 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 paddle deforms. The dark blue dotted circle shows the region of prescribed displacement. Simulation parameters are 
$(\unicode[STIX]{x1D70C}_{f},\unicode[STIX]{x1D707}_{f},\unicode[STIX]{x1D70C}_{s})=(1,10^{-3},2)$.
Figure 9. Snapshots of pressure 
$p$ at time 
$t=6.6$, in a sequence of simulations where flexible paddles with different shear moduli 
$G$ are pushed through a fluid-filled cavity. The paddles are anchored on its right end. 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 paddles deform. The dark blue dotted circles show the regions of prescribed displacement. The colours show pressure using the same scale as figure 8. Simulation parameters are 
$(\unicode[STIX]{x1D70C}_{f},\unicode[STIX]{x1D707}_{f},\unicode[STIX]{x1D70C}_{s})=(1,10^{-3},2)$.
A rectangular flexible paddle of width 1.6 and height 0.4 is initially centred at 
$(0,0.4)$. The displacement of the paddle is prescribed in a circular region of radius 0.15 centred on 
$(0.6,y_{p}(t))$, using the same tethering procedure as in (4.1). The simulation is run for a duration of 
$T=20$, and the circular region moves vertically according to 
$y_{p}(t)=0.4+4.2Y(t/T)$ where
$$\begin{eqnarray}Y(\unicode[STIX]{x1D70F})=\left\{\begin{array}{@{}ll@{}}4\unicode[STIX]{x1D70F}^{2}(3-4\unicode[STIX]{x1D70F})\quad & \text{for }t<{\textstyle \frac{1}{2}},\\ 1\quad & \text{for }t\geqslant {\textstyle \frac{1}{2}}.\end{array}\right.\end{eqnarray}$$ Hence, for 
$0\leqslant t<10$, the paddle is dragged through the fluid, and for 
$10\leqslant t\leqslant 20$ the circular region is stationary and the paddle equilibrates.
Simulations are run using shear moduli from 
$G=1$ to 
$G=10^{7}$, following the standard choices for time step and extra solid viscosity described in the main text. Figure 8 shows a sequence of snapshots of pressure in the simulation with 
$G=100$. As the paddle is pushed through the cavity, it is bent downward due to the pressure of the fluid. Vortices are shed from the paddle tip, creating regions of low pressure visible at 
$t=10$ and 
$t=20$. Figure 9 shows a sequence of snapshots at 
$t=6.6$ for the full range of shear moduli. For 
$G=1$ and 
$G=10$ the paddle deforms so strongly that there is little pressure build-up in the upper part of the domain. However, for 
$G\geqslant 10^{3}$ the pressure build-up is large, and fluid must push through the thin gaps on either side of the paddle. For 
$G=10^{7}$, the paddle becomes near rigid, so that the fluid flow becomes almost symmetric even though the paddle’s motion is only prescribed on the right side. In these simulations the time step is set by the limit from the extra solid viscous term (3.30). The total number of time steps scales according to 
$\sqrt{G}$ and thus the simulation for 
$G=10^{7}$ takes approximately 3100 times more computational resources than that for 
$G=1$.
4.4 Solid actuation
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 modelling active media such as swimmers. Unlike the tethering approach used in § 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, 
$B_{a}$, one writes the desired actuated deformation gradient 
$\unicode[STIX]{x1D641}_{a}(\boldsymbol{X}\in B_{a},t)$, which can then be equivalently expressed in Eulerian frame as 
$\unicode[STIX]{x1D641}_{a}(\boldsymbol{X}=\unicode[STIX]{x1D743}_{a}(\boldsymbol{x}\in b_{a},t),t)$ for 
$b_{a}$ the image of 
$B_{a}$ in the Eulerian frame. At any point 
$\boldsymbol{x}\in b_{a}$, the constitutive relation is adjusted by replacing all references to 
$\unicode[STIX]{x1D641}(\boldsymbol{x},t)$ with 
$\unicode[STIX]{x1D641}(\boldsymbol{x},t)\unicode[STIX]{x1D641}_{a}(\boldsymbol{x},t)^{-1}$. In an isotropic hyperelastic system, for example, this effectively distorts the region’s rest configuration to the distortional state given by 
$\unicode[STIX]{x1D641}_{a}$. If at any moment in time a configuration of the actuated domain differs from the intended actuated configuration, a stress given by 
$\boldsymbol{f}(\unicode[STIX]{x1D641}(\boldsymbol{x},t)\unicode[STIX]{x1D641}_{a}(\boldsymbol{x},t)^{-1})$ emerges that moves the system toward the actuated deformation (where 
$\boldsymbol{f}$ is defined from (3.16)). 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 Kröner–Lee decomposition used in plasticity (Kröner Reference Kröner1960; Lee Reference Lee1969), where a tensorial state variable 
$\unicode[STIX]{x1D641}_{p}$ is introduced and the elastic deformation gradient, which produces the stress, is given by 
$\unicode[STIX]{x1D641}\unicode[STIX]{x1D641}_{p}^{-1}$. But unlike 
$\unicode[STIX]{x1D641}_{p}$, which evolves under a constitutive flow rule, here we assign 
$\unicode[STIX]{x1D641}_{a}(\boldsymbol{x},t)$ directly.
Figure 10. Six successive snapshots of the flapping swimmer (
$Re\approx 200$), with colours showing vorticity 
$\unicode[STIX]{x1D714}$. 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 line marks the fluid–structure interface. The thin dashed lines are contours of the components of the reference map and indicate how the swimmer deforms. Simulation parameters are 
$(\unicode[STIX]{x1D70C}_{f},\unicode[STIX]{x1D707}_{f},\unicode[STIX]{x1D70C}_{s},G)=(1,5\times 10^{-4},4,10)$.
A similar, albeit reduced-dimensional approach has been used to model muscle contraction in swimming lampreys and other narrow-body swimmers. Tytell et al. (Reference Tytell, Hsu, Williams, Cohen and Fauci2010), simulated a lamprey swimming in a two-dimensional (2-D) geometry. The lamprey is modelled with three connected parallel filaments, the outer two of which obey a one-dimensional viscoelastic model with an additional user-defined contractile stress (McMillen, Williams & Holmes Reference McMillen, Williams and Holmes2008). Actuated bending of the lamprey occurs through asymmetric contraction of the filaments. An alternative approach, as used by Gazzola, Argentina & Mahadevan (Reference Gazzola, Argentina and Mahadevan2015) and Patel, Bhalla & Patankar (Reference Patel, Bhalla and Patankar2018), involves direct assignation of an external bending moment on the swimmer cross-section, which in a stiff limit equates to directly assigning the swimmer shape through time. The approach we describe above could be seen as a 2-D (or potentially 3-D) generalization of approaches like these, permitting possibly more through-thickness spatial variation in actuation. The implementation shown here could be made more realistic by including dissipation within the swimmer, neuro-muscular signalling, and contractile-only forcing as was done by Tytell et al. (Reference Tytell, Hsu, Williams, Cohen and Fauci2010), McMillen et al. (Reference McMillen, Williams and Holmes2008) and Patel et al. (Reference Patel, Bhalla and Patankar2018).
As an example, we consider a flapping swimmer (figure 10, Supplemental Movie 5). The swimmer is a rectangle of width 
$W=0.5$ and height 
$H=0.052$ with circular end caps, initially centred on 
$(0,-0.8)$, which we choose to be the location of the origin. We choose the actuated domain, 
$B_{a}$, to be a centred subregion within the swimmer, comprising a rectangle of width 0.28 and height 0.042 with circular end caps. The following actuation is applied:
$$\begin{eqnarray}\unicode[STIX]{x1D641}_{a}(\boldsymbol{X},t)=\left(\begin{array}{@{}cc@{}}\text{e}^{-\unicode[STIX]{x1D6FC}(\boldsymbol{X},t)} & 0\\ 0 & \text{e}^{\unicode[STIX]{x1D6FC}(\boldsymbol{X},t)}\end{array}\right),\end{eqnarray}$$where
$$\begin{eqnarray}\unicode[STIX]{x1D6FC}(\boldsymbol{X},t)=-\unicode[STIX]{x1D706}X_{y}H_{\unicode[STIX]{x1D700}}(-d)\sin ^{8}\unicode[STIX]{x1D714}t=-\unicode[STIX]{x1D706}\unicode[STIX]{x1D709}_{y}(\boldsymbol{x},t)H_{\unicode[STIX]{x1D700}}(-d)\sin ^{8}\unicode[STIX]{x1D714}t\end{eqnarray}$$ and 
$d$ is the signed distance from the Eulerian boundary of 
$b_{a}$. By blurring the boundary of the actuated domain under 
$H_{\unicode[STIX]{x1D700}}(-d)$, it should be noted material positioned up to 
$\unicode[STIX]{x1D700}$ away from the true boundary of 
$b_{a}$ will receive some actuation stress. The parameters used in the simulation are 
$\unicode[STIX]{x1D714}=2\unicode[STIX]{x03C0}/8$, 
$\unicode[STIX]{x1D700}=2.5h_{x}$, and 
$\unicode[STIX]{x1D706}=(\log 2.2)/0.021$. Thus the maximum stretch on the top boundary is 2.2. The simulation uses a 
$1200\times 1200$ grid in 
$[\!-1.5,1.5\!)^{2}$ with periodic boundary conditions. The densities, dynamic viscosity, and solid shear modulus are 
$(\unicode[STIX]{x1D70C}_{f},\unicode[STIX]{x1D707}_{f},\unicode[STIX]{x1D70C}_{s},G)=(1,5\times 10^{-4},4,10)$.
By actuating the flapper in this fashion, the Lagrangian domain 
$B_{a}$, which comprises roughly half the area of the body, is forced to bend periodically in time. The unactuated portion of the swimmer remains passive and flaps as an elastic body in response to be being conjoined to the actuated region. The swimming flapper achieves a Reynolds numbers of 
$Re=V_{solid}^{max}~W/\unicode[STIX]{x1D708}_{f}\sim 200$. Its ability to translate its centre of mass by swimming evidences that this example is not near zero Reynolds number; vortex shedding can be seen for each flap.
4.5 Multi-body contact
Since the reference map technique does not employ moving meshes, it is particularly well suited to problems involving many objects coming into contact. This capability would be useful for a variety of problems, such as studying colloidal mixtures with soft, deformable particles.
To generalize the method to 
$N$ objects, we introduce independent reference maps 
$\unicode[STIX]{x1D743}^{(1)},\unicode[STIX]{x1D743}^{(2)},\ldots ,\unicode[STIX]{x1D743}^{(N)}$ with the ‘
$(j)$’ suffix being used to denote any quantity associated with object 
$j$. For the purposes of exposition, we assume each field is defined as a separate globally defined function that is extrapolated separately, although in reality each reference map only need be defined in a local neighbourhood of each object. Each reference map is updated using (3.1). For a given 
$\unicode[STIX]{x1D743}^{(j)}$, the solid stress 
$\unicode[STIX]{x1D749}_{s}^{(j)}$ is computed using the methods of § 3.3.
When two or more objects come together, their blur zones may overlap, and thus it is necessary to generalize the definition of global stress that was given in (2.8). At a given point, define 
$\unicode[STIX]{x1D706}^{(j)}=1-H_{\unicode[STIX]{x1D700}}(\unicode[STIX]{x1D719}^{(j)})$ to be the solid fraction of object 
$j$. Then the stress is given by
$$\begin{eqnarray}\unicode[STIX]{x1D749}=\left\{\begin{array}{@{}ll@{}}\displaystyle \unicode[STIX]{x1D749}_{f}+\mathop{\sum }_{i}\unicode[STIX]{x1D706}^{(i)}(\unicode[STIX]{x1D749}_{s}^{(i)}-\unicode[STIX]{x1D749}_{f})\quad & \displaystyle \text{if }\mathop{\sum }_{i}\unicode[STIX]{x1D706}^{(i)}\leqslant 1,\\ \displaystyle {\displaystyle \frac{\displaystyle \mathop{\sum }_{i}\unicode[STIX]{x1D706}^{(i)}\unicode[STIX]{x1D749}_{s}^{(i)}}{\displaystyle \mathop{\sum }_{i}\unicode[STIX]{x1D706}^{(i)}}}\quad & \displaystyle \text{if }\mathop{\sum }_{i}\unicode[STIX]{x1D706}^{(i)}>1.\end{array}\right.\end{eqnarray}$$ If only one object is present, this definition exactly matches (2.8). If several objects are present, then they each contribute to the global stress, with the fluid stress filling any unassigned fraction. In rare situations (e.g. three objects meeting at a point) the solid fractions may total more than one. In this case, 
$\unicode[STIX]{x1D749}$ is taken as a weighted average of the solid stresses, and the fluid stress does not contribute at all. The global density field is defined using the same mixing procedure as in (4.10).
In our tests, we have found that independently updating 
$N$ reference maps and computing a global stress according to (4.10) is sufficient to perform multi-body simulations. However, since the simulation employs a single globally defined velocity field, it becomes problematic when shapes become very close together, since it is hard for them to separate as they move according to the same underlying velocity. Similar behaviour has been noted in the literature on the immersed boundary method (Lim & Peskin Reference Lim and Peskin2012; Krishnan, Shaqfeh & Iaccarino Reference Krishnan, Shaqfeh and Iaccarino2017), which also employs a single global velocity field for the movement of structures. To rectify this, we introduce a small contact stress (in addition to the stress of (4.10)) when the blur zones of two objects overlap, which penalizes the interfaces from becoming too close together. We first define a contact force function of
$$\begin{eqnarray}f(x)=\left\{\begin{array}{@{}ll@{}}\displaystyle \frac{1}{2}\left(1-\frac{x}{\unicode[STIX]{x1D700}}\right)\quad & \text{if }x<\unicode[STIX]{x1D700},\\ 0\quad & \text{if }x\geqslant \unicode[STIX]{x1D700}.\end{array}\right.\end{eqnarray}$$ Now, consider the stress calculation at an edge that is within the blur zones of two or more solids. Consider a pair of the solids 
$(i)$ and 
$(j)$. Using finite differences, compute a unit normal vector
$$\begin{eqnarray}\boldsymbol{n}=\frac{\unicode[STIX]{x1D735}(\unicode[STIX]{x1D719}^{(i)}-\unicode[STIX]{x1D719}^{(j)})}{\Vert \unicode[STIX]{x1D735}(\unicode[STIX]{x1D719}^{(i)}-\unicode[STIX]{x1D719}^{(j)})\Vert _{2}},\end{eqnarray}$$ where 
$\Vert \cdot \Vert _{2}$ denotes the Euclidean norm. The contact stress is defined as
$$\begin{eqnarray}\unicode[STIX]{x1D749}_{col}=-\unicode[STIX]{x1D702}\min \{\,f(\unicode[STIX]{x1D719}^{(i)}),f(\unicode[STIX]{x1D719}^{(j)})\}(G^{(i)}+G^{(j)})(\boldsymbol{n}\otimes \boldsymbol{n}-{\textstyle \frac{1}{2}}\mathbf{1}),\end{eqnarray}$$ where 
$\unicode[STIX]{x1D702}=4$ is a dimensionless constant, the 
$G^{(i)}$ are object-dependent shear moduli, and the 
$\mathbf{1}$ term is included to make the stress trace free. In the rare case where the edge is within three or more solid blur zones, the calculation is repeated for each pair, and each contribution 
$\unicode[STIX]{x1D749}_{col}$ is added to the global stress.
These collision stress terms induce forces that push apart objects when they become close. Formulating the collision interaction as an additional stress is advantageous since it immediately ensures that total momentum of the entire simulation is numerically conserved. The method is not sensitive to the exact functional form of 
$f$ in (4.11). An alternative formulation is to directly use the transition function, 
$f(\unicode[STIX]{x1D6FC})=1-H_{\unicode[STIX]{x1D700}}(\unicode[STIX]{x1D6FC})$, but we find that the faster growth of the function in (4.11) when 
$\unicode[STIX]{x1D6FC}$ becomes smaller than 
$\unicode[STIX]{x1D700}$ yields better results in our test cases.
Figure 11. Snapshots of vorticity 
$\unicode[STIX]{x1D714}$ in a simulation of 42 squares sedimenting in a fluid-filled box. The thick black lines mark the fluid–structure interfaces. The thin dashed lines are contours of the components of the reference map defined in each object and indicate how the squares deform. Simulation parameters are 
$(\unicode[STIX]{x1D70C}_{f},\unicode[STIX]{x1D707}_{f},\unicode[STIX]{x1D70C}_{s},G,g)=(1,10^{-3},3,2,1.5)$.
Figure 11 shows snapshots from a multi-body simulation in a non-periodic box 
$[-1,1]^{2}$ using a resolution of 
$1280\times 1280$ with fluid density 
$\unicode[STIX]{x1D70C}_{f}=1$ and dynamic viscosity 
$\unicode[STIX]{x1D707}_{f}=10^{-3}$. Forty-two squares with shear modulus 
$G=2$ and density 
$\unicode[STIX]{x1D70C}_{s}=3$ are inserted at random positions in the box, with side lengths chosen uniformly over the range 
$[0.1,0.4]$. Any squares that lie within a distance of 
$0.1$ of another square are rejected, and are chosen again. At 
$t=0$, each square is set to initially spin with angular velocity chosen uniformly from the range 
$[-5,5]$. A gravitational acceleration of 
$g=1.5$ in the negative 
$y$ direction is applied, so that the squares sediment at the bottom of the box. The full simulation is shown in Supplemental Movie 6.
Figure 12. Snapshots of (a) density 
$\unicode[STIX]{x1D70C}$, (b) pressure deviation 
$p^{\prime }=p+\unicode[STIX]{x1D70C}_{f}gy$ and (c) vorticity 
$\unicode[STIX]{x1D714}$ in a simulation of 75 flexible rods of variable density 
$\unicode[STIX]{x1D70C}_{s}\in [0.4,1.6]$ rearranging in a fluid-filled box. The thick black lines mark the fluid–structure interfaces. Simulation parameters are 
$(\unicode[STIX]{x1D70C}_{f},\unicode[STIX]{x1D707}_{f},G,g)=(1,2\times 10^{-3},1.3,1)$. The plot of 
$p^{\prime }$ at 
$t=0$ is based upon taking a small time step of 
$\unicode[STIX]{x0394}t=10^{-6}$ and computing pressure based on the finite-element projection.
A benefit of the reference map technique is that it can handle both neutrally buoyant solids, and solids that are lighter than the surrounding fluid, without any modification. To demonstrate this, we consider a second multi-body example with 75 solids made of rectangles of length 0.44, thickness 0.044 and rounded end caps, all of which are initially vertically aligned. The solid densities are randomly chosen uniformly over the range 
$[0.4,1.6]$ and a gravitational acceleration of 
$g=1$ in the negative 
$y$ direction is applied. A simulation grid of 
$640\times 1280$ is used on the domain 
$x\in [-0.75,0.75],y\in [-1.5,1.5]$, the dynamic viscosity is 
$\unicode[STIX]{x1D707}_{f}=2\times 10^{-3}$, the shear modulus is 
$G=1.3$, and the fluid density is 
$\unicode[STIX]{x1D70C}_{f}=1$. Figure 12 shows snapshots of the density, the deviation of pressure from the background gradient due to gravity on the fluid 
$p^{\prime }=p+\unicode[STIX]{x1D70C}_{f}gy$, and the vorticity at five time points. In the particular random sample chosen, the average rod density is 
$\bar{\unicode[STIX]{x1D70C}}_{s}\approx 0.904$. Thus the average density in the simulation is slightly lower than 
$\unicode[STIX]{x1D70C}_{f}$ and hence there is a small positive gradient in 
$p^{\prime }$ in the 
$y$ direction at 
$t=0$. The solids separate into two families, with solids with 
$\unicode[STIX]{x1D70C}_{s}<\unicode[STIX]{x1D70C}_{f}$ rising to the top of the domain, and solids with 
$\unicode[STIX]{x1D70C}_{s}>\unicode[STIX]{x1D70C}_{f}$ sinking to the bottom of the domain. While most rods have separated out by 
$t=40$, it takes a long time for the separation process to fully complete, since several rods are close to neutrally buoyant and the reduced gravity they experience is small. By 
$t=120$, all rods have completely separated into two families although there is still some residual movement visible. The density field, pressure deviation field, and vorticity field of the full simulation are shown in Supplemental Movie 7, Supplemental Movie 8, and Supplemental Movie 9, respectively.
5 Conclusion
Herein, we have presented a robustly accurate, yet straightforward to implement, reference map technique, which has allowed us to study a variety of FSI problems using a single background grid. It augments the multi-phase fluid framework of Almgren, Bell and coworkers (Bell et al. Reference Bell, Colella and Glaz1989; Puckett et al. Reference Puckett, Almgren, Bell, Marcus and Rider1997; Almgren et al. Reference Almgren, Bell, Colella, Howell and Welcome1998) by allowing general finite-deformation solid models to be coupled directly to a fluid. In doing so, it maintains a number of the advantages of working on a fixed Eulerian grid that are enjoyed in fluid simulation methods. The practicality and usefulness of this approach is demonstrated in various tests. It is shown to capture the flapping phase diagram for thin flags and the transition from thin- to thick-flag behaviours, which highlights the role of new mechanisms to initiate flapping. Additional physics, such as actuation of solids, is straightforward to implement with a user-described actuated deformation gradient. This capability is used to model a swimming object with realistic internal driving. The ability to model objects with sharp corners is typically a challenge in Eulerian approaches, but here it can be done by exploiting the reference map field near the edge of the object. We also present an improved contact algorithm, which we use to simulate situations with many soft interacting objects submerged in a fluid.
There are a number of future directions. One of clearest applications is in biomechanics, with the simulation of systems of many interacting, actuated cells. We also foresee modelling solids beyond hyperelasticity, such as plasticity, thermal material models and growth. These modifications can be done through the inclusion of new state variables in the solid and/or the addition of a heat diffusion equation; there are clear advantages to implementing thermal diffusion in the Eulerian frame. Beyond extensions to three dimensions, there are opportunities to use the approach for dimensionally reduced models such as membranes and shells by restricting the reference map to a lower dimensional set. Regarding contact modelling, the reference map field could be used to instruct formulations for more advanced contact problems, including friction and self-contact. Lastly, it is a major goal to extend the approach to allow for non-persistent material boundary sets, as occurs in fracture. It may be possible to represent crack surfaces through intersecting level set fields and to couple this capability with physical traction–separation relations to generate new surface material as cracks advance.
Acknowledgements
We thank the three anonymous reviewers for constructive feedback on this work. C.H.R. was partially supported by the Applied Mathematics Program of the U.S. DOE Office of Science Advanced Scientific Computing Research under contract number DE-AC02-05CH11231. Y.Y. acknowledges support from the National Science Foundation under Award No. DMS-1620434 and the Class of 1968 Junior Faculty Fellowship from Lehigh University. K.K. acknowledges support from Army Research Office grants W911NF-16-1-0440, W911NF-18-1-0118, and W911NF-19-1-0431.
Declaration of interests
The authors report no conflict of interest.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2020.353.
Appendix A. Additional numerical details
A.1 Monotonicity-limited derivative
The gradients of the reference map and velocity appearing in (3.5) are computed using the fourth-order monotonicity-limited scheme of Colella (Reference Colella1985). For the derivative of a generic component 
$a_{i,j}$ in the 
$x$ direction, finite differences
$$\begin{eqnarray}D^{c}(a)_{i,j}=(a_{i+1,j}-a_{i-1,j})/2,\quad D^{+}(a)_{i,j}=a_{i+1,j}-a_{i,j},\quad D^{-}(a)_{i,j}=a_{i,j}-a_{i-1,j}\end{eqnarray}$$are introduced, from which the limiting slope is defined as
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}_{lim}(a)_{i,j}=\left\{\begin{array}{@{}ll@{}}\displaystyle 2\times \min (|D^{-}(a)_{i,j}|,|D^{+}(a)_{i,j}|)\quad & \text{if }D^{-}(a)_{i,j}D^{+}(a)_{i,j}>0,\\ 0\quad & \text{otherwise.}\end{array}\right.\end{eqnarray}$$The second-order limited slope is then
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}_{f}(a)_{i,j}=\min (|D^{c}(a)_{i,j}|,\unicode[STIX]{x1D6FF}_{lim}(a)_{i,j})\times \operatorname{sign}(D^{c}(a)_{i,j})\end{eqnarray}$$from which the fourth-order monotonicity limited derivative is defined as
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}^{4}(a)_{i,j}=\min \left(\frac{|8D^{c}(a)_{i,j}-\unicode[STIX]{x1D6FF}_{f}(a)_{i+1,j}-\unicode[STIX]{x1D6FF}_{f}(a)_{i-1,j}|}{6},\unicode[STIX]{x1D6FF}_{lim}(a)_{i,j}\right)\times \frac{\operatorname{sign}(D^{c}(a)_{i,j})}{h_{x}}.\end{eqnarray}$$ The 
$y$-derivative is evaluated similarly.
A.2 Tangential derivatives
To ensure stability, the tangential derivatives appearing in (3.5) and (3.6) are computed using
$$\begin{eqnarray}\displaystyle & \displaystyle (\widetilde{v\unicode[STIX]{x1D743}_{y}})_{i,j}^{n}=\frac{\tilde{v}_{i,j-1/2}^{adv}+\tilde{v}_{i,j+1/2}^{adv}}{2}\frac{\tilde{\unicode[STIX]{x1D743}}_{i,j+1/2}-\tilde{\unicode[STIX]{x1D743}}_{i,j-1/2}}{h_{y}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle (\widetilde{v\boldsymbol{v}_{y}})_{i,j}^{n}=\frac{\tilde{v}_{i,j-1/2}^{adv}+\tilde{v}_{i,j+1/2}^{adv}}{2}\frac{\tilde{\boldsymbol{v}}_{i,j+1/2}-\tilde{\boldsymbol{v}}_{i,j-1/2}}{h_{y}}, & \displaystyle\end{eqnarray}$$where the terms with tildes are computed using a preliminary Godunov upwinding step where stress, pressure and tangential derivatives are neglected (Yu et al. Reference Yu, Sakai and Sethian2003). Extrapolations to a vertical edge from the left are given by
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D743}}_{i+1/2,j}^{L,n+1/2}=\unicode[STIX]{x1D743}_{i,j}+{\textstyle \frac{1}{2}}(h_{x}-u_{i,j}^{n}\unicode[STIX]{x0394}t)\unicode[STIX]{x1D743}_{x,i,j}^{n}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\boldsymbol{v}}_{i+1/2,j}^{L,n+1/2}=\boldsymbol{v}_{i,j}^{n}+{\textstyle \frac{1}{2}}(h_{x}-u_{i,j}^{n}\unicode[STIX]{x0394}t)\boldsymbol{v}_{x,i,j}^{n}, & \displaystyle\end{eqnarray}$$ and with extrapolations to the other edges given similarly. On each edge, with the selection procedure of (3.8) is used, with Case A used for 
$\tilde{\boldsymbol{v}}^{adv}=(\tilde{u} ^{adv},\tilde{v}^{adv})$ and Case B used for 
$\tilde{\unicode[STIX]{x1D743}}$ and 
$\tilde{\boldsymbol{v}}$.
Appendix B. Tests of convergence and accuracy
B.1 Overview of the test configurations
To study the accuracy of the numerical method, we performed a convergence test in the periodic domain 
$[\!-1,1\!)^{2}$ using an initial incompressible velocity field of
$$\begin{eqnarray}\boldsymbol{v}(\boldsymbol{x},0)=\mathop{\sum }_{k=0}^{5}(-1)^{k}\boldsymbol{v}_{vor}\left(x-\frac{-5+2k}{6},y-\frac{-5+2k}{6},2(k+1)\right),\end{eqnarray}$$where
$$\begin{eqnarray}\boldsymbol{v}_{vor}(\boldsymbol{x},\unicode[STIX]{x1D706})=(-\sin \unicode[STIX]{x03C0}y,\cos \unicode[STIX]{x03C0}x)\text{e}^{-\unicode[STIX]{x1D706}(2-\cos \unicode[STIX]{x03C0}x-\cos \unicode[STIX]{x03C0}y)}.\end{eqnarray}$$ This velocity field is designed to have features with a variety of length scales. We simulated up to 
$t=0.5$, used a shear modulus of 
$G=1$, a fluid density of 
$\unicode[STIX]{x1D70C}_{f}=1$, and employed the standard choices for extra viscosity and time-step selection. Using the same initial velocity field, we ran tests using (i) fluid only, (ii) solid only, (iii) a circle of radius 0.6 centred on 
$(-0.1,0)$ and (iv) a square of side length 1.2 centred on 
$(-0.1,0)$. We examined the effect of viscosity and the fluid/solid density ratio.
The configurations of eight different tests are shown in table 1. In our tests, we also considered two different models for the scaling of the extra viscosity. Our primary tests A–F follow § 3.5 and choose it to scale linearly with the grid size. This procedure is consistent with standard numerical schemes; for example, in the second-order Lax–Wendroff method (Lax & Wendroff Reference Lax and Wendroff1960; LeVeque Reference LeVeque2002) the stabilizing diffusive term scales linearly with the grid spacing. However, we also considered alternative tests C’ and F’, whereby the extra viscosity is viewed as a physical dissipation within the solid and is therefore held constant rather than scaling with the grid spacing. It is chosen based on the 
$360\times 360$ grid and then held constant for the higher-resolution grids.
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 
$5040\times 5040$ grid. For each test, we then ran a suite of coarser simulations using 
$N\times N$ grids where 
$N\in \{2520,1680,1260,1008,840,720,630,560,504,420,360\}$ to compare against the reference results. Since each 
$N$ divides evenly into 
$5040$, the grid squares of these coarse simulations align with the reference simulations.
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, 
$360\times 360$, as opposed to scaling linearly with the grid spacing. The last four columns give the exponents of convergence for velocity 
$\boldsymbol{v}$ and pressure 
$p$ under different 
$L^{q}$ norms, based on a linear fit of the three-parameter error model of (B 9) that incorporates a Richardson extrapolation correction. The proportion of Richardson extrapolation correction is shown in italics in brackets.
Figure 13. Differences between the velocity fields in the reference simulation (using a 
$5040\times 5040$ grid) and the coarsest simulation (using a 
$360\times 360$ grid). Plots are shown at 
$t=0.5$ for six of the convergence tests. The colours in each panel are normalized differently by a maximum value 
$V$. The thick black lines mark the fluid–structure interfaces. The thin dashed lines are contours of the components of the reference map. Simulation parameters are 
$(\unicode[STIX]{x1D70C}_{f},\unicode[STIX]{x1D707}_{f},\unicode[STIX]{x1D70C}_{s},G)=(1,10^{-3},1,1)$.
We calculated normalized error measures with respect to 
$L^{q}$ norms
$$\begin{eqnarray}E_{q}^{p}=\left(\frac{1}{A}\int _{\unicode[STIX]{x1D6FA}}|p_{ref}(\boldsymbol{x})-p_{coa}(\boldsymbol{x})|^{q}\,\text{d}\boldsymbol{x}\right)^{1/q},\quad E_{q}^{\boldsymbol{v}}=\left(\frac{1}{A}\int _{\unicode[STIX]{x1D6FA}}\Vert \boldsymbol{v}_{ref}(\boldsymbol{x})-\boldsymbol{v}_{coa}(\boldsymbol{x})\Vert _{2}^{q}\,\text{d}\boldsymbol{x}\right)^{1/q},\end{eqnarray}$$ where 
$A=4$ 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 centred, 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 
$O(h^{3})$ and are small compared to the errors to be measured.
Figure 13 shows plots of the difference in velocity fields between the reference simulation and the simulation on the coarsest grid, for six of the convergence tests at 
$t=0.5$. The colours in the panels are normalized differently, with differences for tests A and C’ being much smaller than those for the other tests that are shown. In the fluid-only test A, the largest errors are present on the diagonal line 
$x=y$, where the initial vortices are located. Some higher errors are visible on thin curved lines, which is a consequence of the switching between cases in the advection discretization. Test C shows the errors for the solid-only simulation, which are about two orders of magnitude larger than test A. In test C’ where the extra viscosity is held fixed, the additional dissipation allows a closer match to be achieved. Tests D, F and F’ show that the largest errors are all near the fluid–solid interface. Unlike models C and C’, there is limited difference between models F and F’, since the errors near the interface dominate.
Figure 14. Plots showing the accuracy of the solutions for different grid sizes 
$h$ for the eight convergence tests given in table 1. Accuracy is computed with respect to reference solutions on a 
$5040\times 5040$ grid. Four accuracy measures are shown: the velocity in the 
$L^{2}$, 
$L^{1}$ and 
$L^{\infty }$ norms, and the pressure in the 
$L^{2}$ norm. The symbols correspond to the computed errors with (B 3) and the lines correspond to the fitted three-parameter error model using (B 9).
B.2 Calculating convergence rates
Figure 14 shows convergence plots for the velocity in the 
$L^{2}$, 
$L^{1}$ and 
$L^{\infty }$ norms, plus the pressure in the 
$L^{2}$ norm; our discussion focuses on velocity, since the pressure can be viewed as a Lagrange multiplier enforcing the incompressibility constraint. If the method is of order 
$s$, and the reference solution is treated as exact, then the errors scale according to
$$\begin{eqnarray}E_{q}^{a}(h)\approx Bh^{s},\end{eqnarray}$$ for some constant 
$B$ where 
$a$ is either 
$p$ or 
$\boldsymbol{v}$. However, in reality the reference solution will not be exact. In particular, one could apply Richardson extrapolation (Richardson Reference Richardson1911; Hairer, Nørsett & Wanner Reference Hairer, Nørsett and Wanner1993; Heath Reference Heath2002) and propose that the numerical solution has a Taylor series in 
$h$, so the leading-order error term is of the form 
$h^{s}$. Specifically, let 
$f(h)\in C^{1}(\unicode[STIX]{x1D6FA})$ be a representation of a component (
$u$, 
$v$ or 
$p$) of the numerical solution computed with grid spacing 
$h$ so that it agrees with the numerical values at the grid points, and smoothly interpolates between them. Then
$$\begin{eqnarray}f(h)=f_{0}+f_{1}h^{s}+O(h^{s+1}),\end{eqnarray}$$ for some smooth functions 
$f_{0},f_{1}\in C^{1}(\unicode[STIX]{x1D6FA})$. Under this assumption the error scales according to
$$\begin{eqnarray}E_{q}^{a}(h)=B(h^{s}-h_{ref}^{s})+O(h^{s+1}).\end{eqnarray}$$ However, for the current method, equation (B 5) is not precisely true, since there are several steps in the numerical method are not Taylor expandable to higher order. The advective terms in the discretization rely on discrete switching between different cases, which manifests as the lines of higher error in tests A and C’ in figure 13. When a grid point leaves the solid, the 
$\unicode[STIX]{x1D743}$ switches from a time-integrated value to an extrapolated value, causing a small, discrete jump in the field, potentially contributing to errors near the boundary.
From figure 13 it is apparent that most of the errors in the domain are smooth, and regular. We therefore propose a model whereby the simulation domain is split into 
$\unicode[STIX]{x1D6FA}_{T}$ where the Richardson error model, equation (B 6), is applied, and 
$\unicode[STIX]{x1D6FA}_{S}$ where the original error model, equation (B 4), is applied. This leads to a three-parameter error model
$$\begin{eqnarray}E_{q}^{a}(h)=B(h^{s}-\unicode[STIX]{x1D6FC}h_{ref}^{s})+O(h^{s+1}),\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6FC}\in [0,1]$ is the proportion of the Richardson error correction. Taking the logarithm of (B 7) yields
$$\begin{eqnarray}\log E_{q}^{a}=C+s\log h\log \left(1-\unicode[STIX]{x1D6FC}\frac{h_{ref}^{s}}{h^{s}}\right),\end{eqnarray}$$ where 
$C=\log B\in \mathbb{R}$. Define 
$X=\log h_{ref}$, and let 
$(x_{k},y_{k})$ be pairs 
$(\log h,\log E_{q}^{a}(h))$ for each lower resolution that was considered. Then 
$s$, 
$C$ and 
$\unicode[STIX]{x1D6FC}$ may be determined by solving the constrained nonlinear least squares problem
$$\begin{eqnarray}\unicode[STIX]{x1D719}(C,s,\unicode[STIX]{x1D6FC})=\mathop{\sum }_{k}(C+sx_{k}\log (1-\unicode[STIX]{x1D6FC}\text{e}^{s(X-x_{k})})-y_{k})^{2}.\end{eqnarray}$$ Equation (B 9) fits all 32 data sets in figure 14 accurately, and the Richardson term correctly captures how most data sets curve downwards. The convergence rates, and proportion of Richardson correction are shown in table 1. The fluid-only tests, A and 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 
$10^{-3}$, and a convergence rate of one is seen across all four metrics in table 1. This is due to the linear scaling of the extra viscosity, which effectively results in changing the physical parameters as the grid spacing changes, approaching the limit of a non-dissipative process as 
$h\rightarrow 0$. In test C’ where the extra viscosity is held fixed, and the physical parameters remain the same, a convergence rate of two is achieved in the 
$L^{1}$ and 
$L^{2}$ norms. Thus second-order accuracy of the solid discretization is achieved, but only for the case where damping is a fixed physical parameter. It remains an open question to find a second-order discretization for a perfectly non-dissipative solid.
Figure 15. Snapshots of vorticity 
$\unicode[STIX]{x1D714}$ in a simulation of three-pronged rotor being spun with an oscillatory motion in 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. The dark blue dotted circle shows the pivot region. Simulation parameters are 
$(\unicode[STIX]{x1D70C}_{f},\unicode[STIX]{x1D707}_{f},\unicode[STIX]{x1D70C}_{s},G)=(1,10^{-2},3,48)$.
The remaining tests, D, E, F and F’ all involve fluid–structure interaction. In the 
$L^{1}$ and 
$L^{2}$ norms, the convergence rate is approximately 
$1.3$ for tests D, E and F, and 
$1.0$ for test F’. As seen in figure 13 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 
$O(h)$, and 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 Reference Gibou and Fedkiw2005; Francois et al. Reference Francois, Cummins, Dendy, Kothe, Sicilian and Williams2006) – 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 
$L^{2}$ and 
$L^{1}$ norms, but results in a lower convergence rate in the 
$L^{\infty }$ norm due to localized effects at the corners.
Appendix C. Performance tests over a range of resolutions
The simulations that are shown in the main text make use of high resolution to ensure accurate results. Here, we demonstrate that the method can work over a wide range of resolutions, and we examine its performance. We consider a three-pronged object whose boundary is described in polar coordinates by
$$\begin{eqnarray}r(\unicode[STIX]{x1D703})=R\frac{1+\unicode[STIX]{x1D6FC}_{3}\cos 3\unicode[STIX]{x1D703}+\unicode[STIX]{x1D6FC}_{6}\cos 6\unicode[STIX]{x1D703}}{1+\unicode[STIX]{x1D6FC}_{3}+\unicode[STIX]{x1D6FC}_{6}}.\end{eqnarray}$$ The simulation domain is 
$|x|\leqslant 1,|y|\leqslant 1$ and no-slip boundary conditions are used on all sides. An 
$N\times N$ grid is used. Parameters of 
$\unicode[STIX]{x1D70C}_{f}=1$, 
$\unicode[STIX]{x1D70C}_{s}=3$, 
$\unicode[STIX]{x1D707}_{f}=10^{-2}$, 
$G=48$ are used, and the simulation duration is 
$T=4\unicode[STIX]{x03C0}$. The shape is parameterized with 
$(R,\unicode[STIX]{x1D6FC}_{3},\unicode[STIX]{x1D6FC}_{6})=(0.8,0.5,0.125)$. The shape is rotated via a pivot centred at the origin of radius 0.2, whose angle is set to
$$\begin{eqnarray}\unicode[STIX]{x1D703}(t)=\left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x03C0}(1-\cos t)\quad & \text{if }t<2\unicode[STIX]{x03C0},\\ 0\quad & \text{if }t\geqslant 2\unicode[STIX]{x03C0},\end{array}\right.\end{eqnarray}$$ following the same pinning method as in (4.1) with a stiffness constant of 
$K_{teth}=10^{-1}\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x0394}t^{-2}$. The r.m.s. angular velocity for 
$T\in [0,2\unicode[STIX]{x03C0}]$ is 
$\unicode[STIX]{x1D714}_{rms}=\unicode[STIX]{x03C0}/2$, Hence, we define the Reynolds number to be
$$\begin{eqnarray}Re=\frac{\unicode[STIX]{x1D70C}_{f}R(\unicode[STIX]{x1D714}_{rms}R)}{\unicode[STIX]{x1D707}_{f}}\approx 100.\end{eqnarray}$$ Figure 15 shows a snapshot of vorticity for six different time points for an intermediate resolution of 
$N=240$. Since the viscosity is higher by a factor of ten from the example in § 4.1, fewer vortices are visible. Figure 16 shows snapshots at 
$T=4\unicode[STIX]{x03C0}/3$ for a range of resolutions from 
$N=20$ to 
$N=1280$. At lower resolution, the accuracy of the solid deformation and the fluid flow is reduced, but the flow is qualitatively similar, and the simulation runs successfully.
Figure 16. Snapshots of vorticity 
$\unicode[STIX]{x1D714}$ in the three-pronged rotor simulation at time 
$t=4\unicode[STIX]{x03C0}/3$ using an 
$N\times N$ computational grid for six different values of 
$N$. 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. The colours use the same scale as in figure 15. The dark blue dotted circle shows the pivot region. Simulation parameters are 
$(\unicode[STIX]{x1D70C}_{f},\unicode[STIX]{x1D707}_{f},\unicode[STIX]{x1D70C}_{s},G)=(1,10^{-2},3,48)$.
Figure 17(a) shows a comparison of wall clock times for a Linux computer with dual 10-core 2.2 GHz Intel Xeon E5-2630 processors. The total simulation time varies over a large range, from 8.1 s for 
$N=20$ to 67.7 h for 
$N=800$ when using a single thread. Using multiple threads creates a substantial speedup, reducing the time for 
$N=800$ to 20.9 h for four threads and 9.4 h for 16 threads; these times correspond to parallel efficiencies of 80.1 % and 44.9 %, respectively. Some parts of the simulation (e.g. the extrapolation routine) are not threaded, and contribute to this loss of efficiency.
Figure 17(b) shows the time-step restrictions due to the liquid viscosity, the shear wave CFL condition and the extra solid viscosity. For small grid sizes, the time step is set by the extra solid viscosity, but for 
$n>640$ the fluid viscosity provides the strongest restriction. Because of this, the total work scales like 
$N^{3}$ for most of the data in figure 17(a), but will eventually transition to 
$N^{4}$ once the fluid viscosity becomes important.
Figure 17. Performance of the simulation code for the three-pronged rotor in appendix C when run on simulations with different 
$N\times N$ grids. (a) Total wall-clock (WC) time as a function of 
$N$ for different numbers of threads. (b) The three time-step restrictions based on the shear wave speed, extra viscosity and fluid viscosity as a function of 
$N$. (c) The average number of V-cycles required to solve the marker-and-cell (MAC) linear system and the finite-element method (FEM) linear system as a function of 
$N$. (d) The wall-clock time for performing a single V-cycle on different grid sizes, using one and sixteen threads, as a function of 
$N$.
A large fraction of the computation time is spent solving the two linear systems for the marker-and-cell projection (§ 3.1.2) and the finite-element projection (§ 3.4). These are solved using a C++/OpenMP multigrid library (appendix F). The library works on any grid size, and uses a hierarchy of grids where the grid dimensions are halved at each successive level. On very small grids, it is inefficient to use threads due to the performance overhead of thread initialization. Because of this, the library self-tunes prior to use to determine the grid level at which to enable threading. Figure 17(c) shows the average number of V-cycles that are performed to reach the required error tolerances, showing that there is limited variation as 
$N$ changes. Figure 17(d) shows the time required to perform a V-cycle as a function of 
$N$. As expected the single-threaded times scale like 
$N^{2}$, proportional to the degrees of freedom. For 
$N\leqslant 100$ the time for 16 threads is similar to the single-threaded performance, since the grids are not large enough for threading to be enabled. However, for large grids, threading results in more than a 
$10\times$ speedup. Across all simulation sizes and thread numbers, the MAC and finite-element linear systems each take between 20 %–30 % of the total computation time.
Appendix D. Effect of the blur zone
As described in § 2, the numerical method uses a transition region of width 
$2\unicode[STIX]{x1D700}$ to blur between the fluid and solid stresses at an interface. The form of the blurring function is given by (2.9) based on previous studies (Sussman et al. Reference Sussman, Smereka and Osher1994, Reference Sussman, Almgren, Bell, Colella, Howell and Welcome1999; Yu et al. Reference Yu, Sakai and Sethian2003, Reference Yu, Sakai and Sethian2007), and throughout the paper we use 
$\unicode[STIX]{x1D700}=2.5h$. Here, we explore the effect of the blur zone width, by running the three-pronged rotor of appendix C using different values of 
$\unicode[STIX]{x1D700}$.
Figure 18. Close-ups of vorticity 
$\unicode[STIX]{x1D714}$ the three-pronged rotor simulation at 
$t=2\unicode[STIX]{x03C0}$ using a 
$240\times 240$ grid for varying values of the blur zone half-width, 
$\unicode[STIX]{x1D700}$. 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. The colours use the same scale as in figure 15. The dark blue dotted quarter-circle shows the pivot regions. Simulation parameters are 
$(\unicode[STIX]{x1D70C}_{f},\unicode[STIX]{x1D707}_{f},\unicode[STIX]{x1D70C}_{s},G)=(1,10^{-2},3,48)$.
Figure 18 shows close-ups of the vorticity field for four different values of 
$\unicode[STIX]{x1D700}$ at 
$t=2\unicode[STIX]{x03C0}$. At this time point, the clockwise motion of the rotor is arrested, which causes momentarily large shear stresses. If 
$\unicode[STIX]{x1D700}=0$, so that there is a hard transition between fluid and solid stress, then this creates numerical noise in the vorticity field near the interfaces. Since vorticity is based on first derivatives, the variations in the underlying velocity field are smaller. In these simulations the tethering force (4.1) is also smoothed using the same value of 
$\unicode[STIX]{x1D700}$, and hence numerical noise is also visible at the edge of the pivot region. Despite the noise, the simulation with 
$\unicode[STIX]{x1D700}=0$ runs adequately. As 
$\unicode[STIX]{x1D700}$ increases the noise at the interface is progressively blurred out. Figure 19(a) shows a plot of the maximum vorticity field in the simulation over time for differing values of blur zone width, highlighting the trade-off between additional noise for small 
$\unicode[STIX]{x1D700}$, and excessive blurring for large 
$\unicode[STIX]{x1D700}$. Our default value of 
$2.5h$ is an acceptable compromise between the two limits.
The precise form of smoothed delta and Heaviside functions has been studied in the literature, and we also considered the delta function studied by Tornberg & Engquist (Reference Tornberg and Engquist2003, Reference Tornberg and Engquist2004) where
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D700}}^{alt}(\unicode[STIX]{x1D719})=\max \left\{0,\frac{1}{\unicode[STIX]{x1D700}}-\frac{|\unicode[STIX]{x1D719}|}{\unicode[STIX]{x1D700}^{2}}\right\}.\end{eqnarray}$$ Equation (D 1) satisfies several discrete moment conditions (Beyer & LeVeque Reference Beyer and LeVeque1992) and has better convergence properties for some problems. We calculated an associated Heaviside function using integration, and ran a variant of convergence test F from appendix B. The measured convergence exponents are within 
$\pm 0.03$ of the original test. Hence, for the current method, errors at the boundary are dominated by other processes, such as extrapolation. Designing better smoothing and extrapolation is an interesting avenue for further study.
Figure 19. (a) Plots of the maximum value of vorticity field as a function of time, in the three-pronged rotor simulations using six different values of the blur zone half-width 
$\unicode[STIX]{x1D700}$. (b) Plots of the difference in the Jacobian 
$J$ from unity in the solid for the three-pronged rotor simulations for several 
$N\times N$ grids. A nonlinear square root scale is used on the vertical axis to account for the substantially larger differences in 
$J$ for the simulation without the MAC projection.
Appendix E. Volume conservation
Since the reference map field represents an incompressible solid, the determinant of the deformation gradient 
$J=\det \unicode[STIX]{x1D641}$ should remain equal to unity throughout the simulation, but in the numerical scheme this property will only be maintained approximately. This issue is encountered in other incompressible fluid–structure interaction approaches, such as the immersed boundary method where the Lagrangian description of the solid is updated from the fluid velocity and may experience volumetric changes over time (Wang, Zhang & Liu Reference Wang, Zhang and Liu2009; Griffith Reference Griffith2012; Vadala-Roth et al. Reference Vadala-Roth, Acharya, Patankar, Rossi and Griffith2020).
To investigate the magnitude of volumetric deviations, we simulate the three-pronged rotor of appendix C and compute the field 
$J-1$. The field is evaluated on all cell corners in the interior of the solid, where the bilinear interpolation of 
$\unicode[STIX]{x1D719}$ is negative. Calculating 
$J$ requires the gradient of the reference map, which is computed at each cell corner using centred finite differences of the four adjoining cell-centred reference map fields, some of which may be based on extrapolation.
Figure 20. Snapshots of volumetric change in the solid, 
$J-1$, in several simulations of the three-pronged rotor on 
$N\times N$ grids. 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. The dark red dotted circle shows the pivot region. The bottom right panel is from a simulation where the MAC projection is switched off. Simulation parameters are 
$(\unicode[STIX]{x1D70C}_{f},\unicode[STIX]{x1D707}_{f},\unicode[STIX]{x1D70C}_{s},G)=(1,10^{-2},3,48)$.
Figure 20 shows snapshots of 
$J-1$ for several simulation configurations. The top row shows three snapshots on a 
$240\times 240$ grid. Small volumetric deviations of the order of 
$10^{-2}$ are visible and are concentrated near the fluid–solid interface, and at the edge of the pivot region. Figure 19(b) shows a plot of 
$\Vert J-1\Vert _{2}$ as a function of time. The volumetric deviations grow rapidly up to 
$t=\unicode[STIX]{x03C0}/4$ but then remain relatively stable; in particular, for 
$t>2\unicode[STIX]{x03C0}$ when the pivot stops rotating, 
$\Vert J-1\Vert _{2}$ tends to a constant value. The volumetric deviations have the expected dependence on resolution: as shown in figures 19 and 20 the errors are larger for 
$N=160$ and smaller for 
$N=640$. The values of 
$\Vert J-1\Vert _{2}$ at 
$t=4\unicode[STIX]{x03C0}$ scale approximately like 
$h^{1.1}$, consistent with the convergence rates in table 1. We find that the MAC projection substantially improves the volume conservation of the reference map: switching off the projection results in larger deviations overall, plus a tendency for them to grow once the pivot stops rotating (figure 19). This can cause considerable errors near the pivot, as shown in the bottom right panel of figure 20.
Appendix F. Computer code
All results in this paper were created using a custom C++ code that uses OpenMP for multithreading. The code is released as an open source software package, IncRMT, via a GitHub repository at https://github.com/chr1shr/incrmt. As discussed in appendix C, a large fraction of the computation time is spent solving the MAC projection and finite-element projection to enforce incompressibility. This is done using the Templated Geometric Multigrid library, which is available on GitHub at https://github.com/chr1shr/tgmg.
