Eulerian simulation of complex suspensions and biolocomotion in three dimensions.

We present a numerical method specifically designed for simulating three-dimensional fluid-structure interaction (FSI) problems based on the reference map technique (RMT). The RMT is a fully Eulerian FSI numerical method that allows fluids and large-deformation elastic solids to be represented on a single fixed computational grid. This eliminates the need for meshing complex geometries typical in other FSI approaches and greatly simplifies the coupling between fluid and solids. We develop a three-dimensional implementation of the RMT, parallelized using the distributed memory paradigm, to simulate incompressible FSI with neo-Hookean solids. As part of our method, we develop a field extrapolation scheme that works efficiently in parallel. Through representative examples, we demonstrate the method's suitability in investigating many-body and active systems, as well as its accuracy and convergence. The examples include settling of a mixture of heavy and buoyant soft ellipsoids, lid-driven cavity flow containing a soft sphere, and swimmers actuated via active stress.

Fluid–structure interactions (FSI) are at the heart of many physical and biological problems, including flexible structures in flow (1, 2), blood circulation in the heart (3, 4), animal locomotion (5, 6), and cilia motion (7, 8). The couplings between fluid and immersed solids give rise to complex nonlinear dynamics dependent on geometry and boundary conditions, material constitutive relations, and collective interactions among the solid objects. Analytical solutions are rare and limited to simplified settings in reduced dimensions; numerical methods for FSI have become indispensable for understanding these problems. Although methods in reduced dimensions have led to invaluable insights, more aspects of the rich FSI dynamics await to be revealed and understood by fully three-dimensional treatments.

In designing numerical methods for fluids and solids, Eulerian and Lagrangian perspectives are the more convenient choices, respectively, due to the differences in constitutive responses. Bridging between the two perspectives is a classic dilemma in developing numerical methods for FSI. Various frameworks have been proposed to resolve this dilemma. The immersed boundary method (IBM) (3, 9) and the family of immersed methods that it inspires (4, 10, 11) solve fluid on a fixed mesh, use Lagrangian points to represent the solid, and employ a coupling scheme between the two. Arbitrary Lagrangian–Eulerian methods use moving nodes for both phases and reposition the nodes to maintain mesh quality (12, 13).

There are also fully Eulerian methods, which typically recast the solid momentum balance equations in the rate form, matching the Navier–Stokes equation for the fluid, and solve the coupled equations in the Eulerian space. Formulations of fully Eulerian methods can be broadly divided into two categories. The hypoelasticity formulation (14–16) uses linear elasticity, whereas the hyperelasticity formulation employs a general large-deformation description in the solids. To compute solid stresses in hyperelasticity, various quantities have been used. Examples include level-set functions defining the fluid–solid interface (17), the deformation gradient tensor (18), the deformation tensor (19), and the solid displacements in the undeformed configuration (20–22).

The reference map technique (RMT) is a fully Eulerian, hyperelastic FSI approach that uses the reference coordinates as the primary simulation variables in the solids. In two dimensions (2D), the RMT has been demonstrated to simulate compressible fluid and solids (23, 24), handle contact between multiple solids (25), and simulate incompressible fluid and solids with complex geometries and actuation (26). The RMT also has many properties favorable for addressing computational challenges in three dimensions (3D). For instance, by using a single fixed regular grid for fluid and solid, the RMT is well suited to parallelization and allows for efficient memory usage, since no Lagrangian mesh topology needs be stored. In many-body problems, the regular grid structure makes contact among solid bodies easy to detect. Computational meshing of complex geometries is a difficult problem and an active area of research (27–29). By using level sets (30, 31) to represent the fluid–solid interfaces we circumvent the need for meshing, which becomes a compelling advantage in 3D. However, in all prior works on the RMT in the FSI context (24–26), the fast-marching method (FMM) (16, 31) has been used to reinitialize the level sets at every time step, as part of a necessary field extrapolation routine near the interfaces. The FMM is not well suited for parallel computing since it updates field values sequentially. Therefore, we developed a method for extrapolation that removes the need for level-set reinitialization and can be easily parallelized.

Equipped with these enabling simplifications, we develop a three-dimensional implementation of the RMT parallelized via domain decomposition and then apply it to simulate several difficult FSI scenarios. Inspired by recent work on fluid-filled soft granular packings (32), we simulate a mixture of deformable particles, each of which can be either denser or lighter than the fluid. We also revisit the classic problem of lid-driven cavity flows (33–35), now containing an immersed deformable sphere, and provide detailed data to serve as benchmarks. There has been sustained interest in swimming and biolocomotion in recent decades, yet reduced-dimension models of the swimmer (36–38) are often the state of the art. We present a model of a swimming organism using full volumetric actuation, demonstrating the RMT as a simulation tool for exploring a broader range of swimming modalities.

Theory and Numerical Method

Hyperelastic Formulation in Solid Mechanics.

In the hyperelasticity framework of solid mechanics (39), a time-dependent mapping determines how the undeformed configuration, , is transformed to its current physical configuration, ; i.e., (Fig. 1 ). The deformation gradient is defined as . A constitutive relation defines the Cauchy stress response in the solid material. The solid momentum balance equation in rate form iswhere ρ is the solid density, is the solid velocity, and is a body force density. For an incompressible material and thus solid density is unaffected by the deformation. We assume that is sufficiently smooth so that its inverse, the reference map, can be defined; i.e., . The deformation gradient tensor becomes

Fig. 1.

Schematics of the reference map technique. (A) An initially undeformed solid with reference configuration undergoes a time-dependent mapping to its current configuration at time t. (B) A level-set function distinguishes the two phases that share a global velocity . A blur zone (yellow), defined by , is used to transition between phases. (C) The order in which is extrapolated is defined by layers. The first layer cells (green), e.g., cell (i, j, k), are noninterior neighbors to the interior cells (yellow), e.g., cell . Subsequent layers, e.g., second layers (red), are constructed in the same way, until the blur zone is filled or a physical boundary is reached. The 2D schematics are shown for clarity.

where is the gradient operator in physical space. The reference configuration is constant; therefore, ; i.e.,

We discuss coupling fluid and solid phases and imposing the incompressibility constraint next.

Blurred Interface Method and Monolithic Governing Equations.

Consider a domain Ω containing n immersed solid objects covering subdomains . Denote the fluid domain as Ω and the solid domain as , so that . The fluid–solid interface (hereinafter the interface) is denoted by . We introduce a global velocity in Ω, shared by both fluid and solid phases. This implies that the velocity is continuous and satisfies a no-slip boundary condition at . The incompressibility constraints in both phases imply

As a result, we can define a global pressure field p in Ω, and we need only consider the deviatoric part of the stress tensors, , in both phases. The momentum balance equation becomeswhere density ρ, stress , and body force density , depending on the positions in the domain, can belong to the fluid or the solid phase. In Ω, is the deviatoric viscous stress of a Newtonian fluid, , where μ is the fluid dynamic viscosity. In Ω, is the deviatoric solid stress , obtained from hyperelastic formulation mentioned above. Similarly, ρ and are defined via in Ω and in Ω. Eqs. – form the governing equations for the coupled FSI system, satisfying the flow equation and the elasticity equation, in the fluid and the solid phases, respectively. While Eqs. and are solved in Ω, to avoid excessive distortions, we restrict the solution of Eq. to only the solid domain Ω ().

Across the interface, quantities , and may be discontinuous and should be carefully treated with an interfacial coupling method. The basic RMT equations work with a variety of interfacial coupling procedures, including sharp and blurred interface methods (24–26). In this work we focus on a blurred interface method because it has several advantages over a sharp interface method; e.g., it is more stable to interfacial perturbations, more amenable to simulating immersed solid–solid contact, and easier to implement (25, 26). In the blurred interface method, we make use of a blur zone of width across the interface. The width of the blur zone scales with the grid spacing so that as grid spacing approaches zero, a sharp interface representation is recovered.

Without loss of generality, we consider a single solid object whose boundary is defined by the zero contour of a signed-distance function in the undeformed configuration, . We define the time-dependent level-set function (Fig. 1). In the blur zone defined by , the material can be considered to be a mixture, and discontinuous quantities are blended to smoothly vary across the interface (40–42). Let denote a scalar or a component of a vector or tensor quantity in Ω. Define as the fluid value of and as the solid value of . Omitting for brevity, we blend Q and Q by where is a smoothed Heaviside function that has a continuous second derivative, defined as

The blending procedure is similar to interpolating between Eulerian and Lagrangian quantities in the IBM via integration using approximate delta kernels (3, 9). Due to the stress tensors being blended in the blur zone, the traction-matching condition is automatically satisfied at the interface.

To improve numerical stability, we define an artificial viscosity μ in the solid domain Ω. If μ scales with grid size, in the limit of very fine spatial resolution, we recover the undamped solid equation. If it is set to a grid size-independent constant, it is equivalent to simulating a Kelvin–Voigt viscoelastic solid. We also find that to stabilize the interface, it is necessary to apply an additional multiplicative factor γ to the artificial viscosity in the blur zone. These modifications give rise to an artificial viscous stress in the solid domain Ω and in the blur zone. Overall, , where ().

Extrapolation of .

Since is defined only inside the solid, it needs to be extrapolated to several grid cells outside of the interface to calculate derivatives in Eq. and near the interface. We describe our extrapolation method for a single solid occupying domain Ω, although it can be easily applied to any number of objects. First, we simplify the spatial order in which is extrapolated by making use of adjacency rules on a fixed grid. Consider a Cartesian mesh with grid cells indexed by ; we define neighbors of a central cell as those cells that share a common face with the central cell (Fig. 1). We first label cells with as the interior cells (or 0th layer cells). Then, for each interior cell, we label its unmarked neighbors as first layer cells, with index l = 1. We repeat this procedure to find subsequent layers , until we reach a physical boundary or the maximum number of layers, whichever occurs earlier. This procedure guarantees that cells in the current layer (e.g., 1st) are selected from unmarked neighbors of cells in the previous layer (e.g., 0th). The minimum number of layers should ensure all cells in the blur zone are extrapolated and necessary finite-difference operations can be performed near the interface, where the maximum can be arbitrarily chosen.

The extrapolation is then performed in ascending order of layers, but it can be computed independently for each cell within the same layer. To extrapolate at cell (i, j, k) in layer , we first use weighted least-squares regression to build a local linear model of the reference map as a function of physical coordinate . To find data points for the regression, we search within an r = 2 box centered at cell (i, j, k), where r is the search radius measured in units of grid cells. Reference map at cell (p, q, r), extrapolated or not, is used for the regression only if it is marked in a layer and . In the case of multiple solid objects, must emanate from the same object as . If the linear system is degenerate, we increase r by 1 and repeat the procedure, but this is rare in practice. Using appropriate weights is important to ensure the quality of the extrapolated values, especially when local deformations are large ().

In multibody simulations, the extrapolation procedure is applied to each object independently. We require that solid bodies do not coexist at a grid cell. However, the blur zone of an object is allowed to overlap with blur zones or interiors of other objects. Thus, at a single grid cell, there can be several reference maps, each belonging to a distinct object. Since extrapolated values are needed only in a small region near the interface, we design a custom data structure that is tailored to store extrapolated values efficiently in memory. Besides eliminating the need of reinitialization, the current extrapolation method has two additional advantages: 1) Layers can be defined given a definition of adjacency on the grid, and 2) the method is layer-wise and object-wise independent and thus easy to be parallelized.

Numerical Procedures and Implementation.

The RMT implementation in 3D (RMT3D) is developed in C++ and parallelized via domain decomposition using the message-passing interface (MPI) library. The current second-order accurate, explicit numerical schemes extend our previous work on the RMT implementation in 2D (26) and follow established discretizations for solving hyperbolic conservative laws (42, 43). In summary, we extend the variable arrangement on the grid, the finite-difference schemes to compute spatial derivatives, and the Godunov-type upwinding scheme to handle the advective parts of Eqs. and (44). An approximate projection method (43, 45, 46) and a marker-and-cell projection method are used to enforce the incompressibility constraint (Eq. on the velocity solution at each time step and on an intermediate velocity field between two time steps, respectively. Large linear systems resulting from the projections are solved using a custom geometric multigrid solver. See for details of the numerical procedures including the time-stepping algorithm, stability and time-step restrictions, and the Godunov-type upwinding calculations.

Results

In this section, we consider immersed viscoelastic neo-Hookean solids (constant μ) in various settings. We nondimensionalize the governing equations using appropriate length, time, and mass scales and simulate using isotropic grid spacing h in all test cases.

Parallel Computation Scaling.

We conduct strong and weak scaling tests to examine the efficiency of the parallel implementation. For a system of 2563 grid points, we find that the speedup in strong scaling tests saturates at a factor of as MPI ranks increase from 1 to 512, while the efficiency in weak scaling tests drops to as MPI ranks increase from 1 to 512. In general, the extrapolation routines take of the total compute time for each immersed object, supporting our claims of its efficiency and of the RMT’s advantage in simulating many-body FSIs. See for more details.

Convergence: A Prestretched Sphere.

We simulate an immersed sphere with radius 0.2, shear modulus G = 0.1 in a cubic domain with unit side length and no-slip boundaries. The sphere is prestrained with stretches in the x, y, z directions, respectively. At T = 0, the sphere starts to relax to equilibrium. After initial oscillations, it eventually comes to rest. We compute convergence rates to be 1.34, 1.44, and 1.04, for solid surface area, volume, and , respectively. , where J is the determinant of , measures the deviation from the incompressibility constraint in the solid. Although our numerical schemes are second-order accurate when applied to full fluid or full solid problems, the errors near the interface reduce the overall order of accuracy in coupled FSI problems (26). For more details see .

Settling and Floating.

The transport of rigid and deformable particles in fluid flow is central in many biological and physical systems. While analytical solutions are available for simple cases such as a rigid sphere in unbounded creeping flow, more complex settings, e.g., those with walls and multibody interactions, elude analytical approaches. As a validation, we apply the RMT to simulate a rigid sphere settling in a domain with square cross-section and find good agreement of the position and velocity of the settling sphere between simulations and experiments (47, 48). See for details. The RMT requires no special treatment to simulate (neutrally) buoyant solids, which is a common numerical difficulty suffered by partitioned FSI methods (49). This is an important advantage as many FSIs of interest involve such density ratios in the solid and fluid phases, e.g., problems in hemodynamics (50) and biomechanics (7). To demonstrate this, and the RMT’s ability to simulate complex suspensions (51), in Fig. 2 we show 150 ellipsoids settling and floating in a box, where the particle Reynolds number for the lighter particles, and for the heavier ones (Movie S1).

Fig. 2.

(A–D) Various stages of 150 ellipsoids settling/floating in a square cylinder (Movie S1). All ellipsoids have major axis R = 0.13 initially aligned with the z direction, aspect ratio 2:1:1. Half of the ellipsoids are buoyant, (yellow), and the others are denser than the fluid, (orange). Other parameters are . No-slip boundary conditions are applied on all the walls. (A) Ellipsoids initially at rest and randomly dispersed in the fluid. In B and C cross-sections at y = 0.5 are shown, and color corresponds to the y component of vorticity. Fluid–solid interfaces are plotted with a thick black line, and contours of reference map components ξ and ξ are plotted with black and blue dashed lines, respectively. The particle Reynolds number for lighter ellipsoids, and for the heavier ones. are, respectively, the time average velocity of lighter and heavier particles while they are in motion. They are computed only from particles that approximately traverse the entire domain vertically. (D) Ellipsoids are fully separated with the lighter ones on top and the heavier ones at the bottom. This simulation was run with 48 MPI ranks on Intel “Cascade Lake” nodes and took 19 h.

Lid-Driven Cubic Cavity with a Sphere.

In computational fluid dynamics the lid-driven cavity has long been an important benchmark problem (33, 52). Despite its simplicity, it exhibits rich flow dynamics due to varying geometries, boundary conditions, and Reynolds numbers. In stark contrast to the extensive studies on the fluid problem in both 2D and 3D, results on lid-driven cavities with deformable boundaries and immersed solids are much fewer (20, 34, 35, 53–55).

Here we investigate a neutrally buoyant deformable sphere in a lid-driven cavity. Shown in Fig. 3, the cavity has size and two span aspect ratios, and . The lid moves with velocity and no-slip boundary conditions are applied on the other walls. We rescale length and velocity by L and u, respectively. As a validation, we first simulate lid-driven cavity flows without a solid. Our results agree well with high-accuracy benchmarks (33) (). A deformable circle in a square lid-driven cavity in 2D has been investigated by Zhao et al. (34) and widely used as a validation case in later works (19, 54, 56, 57). In simulating a sphere in a cubic cavity, we first choose parameters similar to those in the 2D test case, i.e., , to highlight qualitative differences in 3D (Fig. 3 ). The middle cross-section of the deformed sphere (Fig. 3) is similar to the shape of the deformed circle at long time (54). The distinctions between the two cases are more apparent in their centroid trajectories. Although the centroid of a sphere in a cubic cavity (Fig. 3) also converges to a stationary point, spirals of its trajectory are much closer together than of the circle in 2D. There are several reasons for this, most notably the topological difference between an infinite cylinder and a sphere. More importantly, the reduced circulation due to lateral walls in the third dimension (52) allows the sphere to interact with the moving lid for longer before being advected back to the center of the cavity.

Fig. 3.

Simulations of a sphere in the lid-driven cavity. For all cases, parameters , and the flow Reynolds number . The sphere has radius 0.2 and is initially centered at . (A) Snapshots of a 3D simulation at various times, . (B) Contours of the reference map and the interface are plotted against a background of velocity magnitude (heatmap and contours). A cross-section at y = 0.5 and T = 50 is shown, . (C) The trajectory of the sphere centroid from T = 0 to T = 50 of the sphere in B. The color corresponds to the velocity magnitude of the centroid. (D and F) The same as B but for solid shear moduli G = 0.25 and 0.5, respectively. (E and G) The same as C but for solid shear moduli G = 0.25 and 0.5, respectively. (H–J) The 3D shape of a deformed G = 0.1 sphere at the closest approach to the top in cavity with , respectively. For additional data, see ; Movies S2–S5; and Dataset S1.

To show the effect of the lateral walls on the solid deformation, we simulate a sphere with G = 0.1 initially at in cavities with . Fig. 3 shows two views of the sphere at the closest approach to the top lid in each cavity. As the walls become farther apart, the sphere becomes less stretched lengthwise and less compressed vertically, again suggesting that a reduced circulation increases the strength of the interaction between the sphere and the moving boundary. Additionally, the centroid trajectories in longer cavities more closely resemble the trajectory in 2D (19) and have farther-spaced spirals ().

We also simulate spheres with varying shear moduli. Snapshots of a simulation with G = 0.03 are shown in Fig. 3. In Fig. 3 we show that the shear modulus significantly influences the centroid trajectory. Fig. 3 offers some intuition for the qualitative changes. As a stiffer sphere moves toward and along the top lid, it deforms less and is able to separate earlier from driving flow. Consequently, a stiffer sphere is advected more toward the center than toward the bottom and thus able to find the stationary position faster. We do not impose any repulsive forces on the sphere near the walls to avoid interfering with its dynamics. For coarse resolution and low G simulations, we find is necessary to stabilize the interface near the boundaries ().

Swimming.

Swimming has been an FSI problem of interest for decades (58–60). To address the difficulty of resolving the full FSI, especially motions of the phase boundary, modelers often apply simplifications such as asymptotic analysis and scaling arguments. A related, widely applied numerical approach in the low Reynolds number regime is to abstract the swimmer into a one-dimensional collection of regularized singularities, which yields good approximations the swimmer’s far field (6, 38, 61–64). However, this method requires specialization when the near field is of prime importance, such as swimming in crowded or confined environments (65, 66). We now use the RMT to simulate a finite-size swimmer in a small no-slip box, showing the method naturally resolves the near field.

Consider a neutrally buoyant swimmer with a cylindrical flagellum of length L and radius R and a spherical head of radius as shown in Fig. 4. We decompose the solid deviatoric stress into passive and active parts , where is the elastic stress tensor from previous sections. We seek to define a time-dependent active stress field as a function of body position to induce cyclic deformation in the form of a planar bending wave traveling along the swimmer’s body. First, we specify how body orientation is determined in the deformed frame. We assume the reference frame’s coordinate system is centered on the flagellum and aligned with the swimmer. Denoting the associated orthonormal basis as , we orient so points along the body toward the head and vertically up. The directions of these vectors in the deformed space, , where , are also body aligned as shown in Fig. 4.

Fig. 4.

Swimming of flagellar objects driven by active bending moments. (A) Schematic of a swimmer. Swimmers have a spherical head of radius R and a cylindrical body of length L and a body radius . The body is actuated by a traveling wave of bending moments B applied to vertical cross-sections (marked in red). (B) Expanded view showing the active in-plane stresses directed along the body’s long axis . The stress magnitude varies linearly with the height above the horizontal midplane (dashed blue line). (C and D) The swim speed U (C) and efficiency e (D) as a function of W, the bending wave amplitude, are presented for a variety of body shapes. Simulations in C and D, Top vary R with L = 1.5, and those in C and D, Bottom vary L with R = 0.15. (E–G) Body shapes, midstroke, are shown for several R, L combinations (Movies S6 and S7). Steady-state measurements are taken at time t = 10. For all cases, .

Next, we define the active stress in terms of the reference map components and , which denote the reference distance along the flagellum and above the midplane, respectively. Since bending moments are induced by axial stresses of opposite signs about the midplane, as shown in Fig. 4, we set . Setting the full deviatoric tensor, we let where for some wavenumber k and frequency ω. We introduce an amplitude parameter , where I is the cylinder’s cross-sectional area moment of inertia and B a bending-moment magnitude, to compare bending moments across a range of R. Finally, we let

Nondimensionalizing position, time, and stress as , and , we simulate swimmers with varying R, W, and L at constant ρ, ρ, μ, and R. We calculate the swim speed U (Fig. 4) and active power , where denotes time averaging over many cycles of oscillation. These give rise to an approximate Lighthill efficiency (Fig. 4), where C is a drag coefficient used to estimate the force required to tow the swimmer at velocity U (67). At large W, the efficiency scales as , consistent with . Midstroke body shapes for three swimmers are shown in Fig. 4 . We find that there is effectively no motion at L = 1, suggesting a minimum length is required for locomotion. Swimming speed U and efficiency e increase as the swimmer body size increases. This may be due to the scaling of Reynolds number describing the time-averaged motion, since it grows as over the range of simulated W values. The swimmer gait Reynolds number varies more slowly as . Both for the results in Fig. 4. For more details, see and Movies S6 and S7.

Discussion

The reference map technique is a flexible numerical method for FSI problems and is a useful complement to other numerical FSI approaches. The strengths of the method are in simulating volumetric elasticity, in handling complex geometries, and in many-body contact. It is also well suited to coupling to other physical processes that are naturally simulated on an Eulerian grid.

To understand the relative merits of the RMT, we compare it to the family of immersed boundary methods (3, 9) for simulating the lid-driven cavity problem with an elastic sphere on an grid. For the same resolution, the computational requirements for simulating the fluid and solid are likely similar, although the RMT would not require the solid mesh topology to be stored. The largest difference between the two methods is in interphase coupling. At each solid node in the IBM, it is necessary to transfer information to the fluid grid via a smoothed delta kernel, usually via a set of grid points, creating a amount of work. In the RMT it is necessary to perform extrapolation at the boundary of the solid, creating an amount of work. Thus, as the grid is refined, the overhead of the RMT will be lower, highlighting its advantages for efficiently solving problems involving volumetric elasticity.

By contrast, the RMT currently cannot handle reduced-dimensional elastic structures (10, 36) that are straightforward with the IBM. Furthermore, the implicit description of the interface makes it less suited to problems that have large gradients near the interface, such as flows at high Reynolds number. A possible future direction is to formulate the RMT using adaptive mesh refinement, which may mitigate this issue. Alternatively, it may be possible to formulate the RMT using a high-order sharp-interface approach (68), which would extend the range of applicability of the method.