Eulerian Method for Multiphase Interactions of Soft Solid Bodies in Fluids

The challenges of simulating fluid–structure interaction (FSI) have been approached from many directions. One set of challenges stem from domain discretization: fluid problems on their own are amenable to solution on an Eulerian computational domain [1–5] and solid deformation is natural to compute within a Lagrangian framework [6–9]. To couple solid and fluid phases in one setting, several approaches have been proposed to mitigate this separation in methodology. One approach is to treat the solid with a standard Lagrangian finite-element framework and to apply an Arbitrary-Lagrange Eulerian method in the fluid [10–12], which remeshes the fluid domain in order to remedy issues of excessive mesh deformation. Other approaches include the family of immersed methods [10,11,13], which keep an ambient stationary Eulerian grid throughout, on which fluid flow is solved, as well as a moving collection of interacting material points representing the solid structure. Here, discretized delta functions are used to pass information between the grid and the nodes.

Some advantages of a fully Eulerian method—fluid and solid both computed on an Eulerian grid—can be directly seen. Since all phases are solved by sweeping through a single fixed mesh, there are computation time advantages. Topological changes are easy to manage on a fixed grid using level sets to track interfaces [4,14]. Furthermore, multiscale and multiphysics coupling can also have advantages when done on an Eulerian grid. One specific example, which will be discussed in more detail later, is the case of multiple solids making contact immersed in a fluid. Finding contact between two solid phases on an Eulerian grid can be achieved using gridwise distance functions or simply identifying grid points which become occupied by multiple solid phases during a trial step.

To achieve these goals, one must pose an Eulerian scheme capable of solving finite-deformation solid problems. The recently proposed reference map technique (RMT) is such an Eulerian framework, based on tracking the reference map field [15,16]. Other approaches for solid deformation on a fixed Eulerian grid include hypoelastic implementations [17,18], which may succeed for small elastic strains but lack a thermodynamically consistent form as needed for large deformations, and methods that directly evolve the deformation gradient tensor field as the primitive kinematic grid variable [19,20,21]. In a previous paper [16], the RMT demonstrated the capability of accurately solving hyperelastic solid deformation problems on a fixed mesh—including shock propagation problems and problems with varied boundary and initial conditions—up to second-order accuracy in space and time. It also provided the first demonstrations of using the method to solve fully coupled problems of FSI. There, the FSI method hinged on a sharp-interface representation, extending on that of the Ghost Fluid Method for fluid–fluid interaction [22]. Sharp methods make a distinct separation between each phase down to the subgrid level. In the current work, our efforts exploit a blurred interface, a simpler and computationally faster implementation, involving fewer numerical extrapolations. A blurred interface method uses a thin transition zone where one phase converts into the other. As the grid size decreases, the corresponding transition zone decreases, and results approach that of a sharp interface method. Here, we show that the blurred interface approach has advantages over the sharp, most important of which is the ability to represent key mechanical behaviors such as submerged solid–solid contact.

To satisfy subgrid jump conditions, sharp interface methods require a large number of gridwise extrapolations of kinematic and stress fields across the interface. These produce the “ghost values” for each phase, which represent an extension of each field into the region occupied by the other phase(s). The validity of a sharp interface method is limited by the quality of continued and progressive extrapolation. When used in the FSI method of the previous work [16], accrued extrapolation error can have a destabilizing effect—shots of pressure along the interface may erroneously appear when the interface crosses through grid cell boundaries, in response to the sampling of extrapolated values when new points enter the solid domain. Adding significant solid dissipation or surface tension can penalize these artifacts, but this can alter the physicality of the simulation. We have found that errors of this type often produce routine-ending numerical instabilities, posing a serious implementational issue. By contrast, in fluid–fluid interaction methods, this effect is less important due to the natural viscous damping within all phases.

By switching the interfacial treatment to a blurred method, we present a compromise wherein we move away from the subgrid interfacial representation of a sharp technique in order to achieve a faster, simpler, and more stable method. Existing fluid–fluid blurred techniques [23,24] do not require any extrapolation. In the case of a fluid–solid problem as shown herein, using only one extrapolation—that of the reference map (described in more detail in Sec. 2.1)—we are able to simulate FSI with a hyperelastic solid body on a single fixed grid using a single velocity field for the whole computational domain. In addition, jump conditions do not need to be explicitly applied but are implicitly met inside the transition zone. Moreover, as discussed in Sec. 7, the fact that the method uses one universal velocity field valid for all phases permits a straightforward approach for simulating multibody contact of solids submerged in a fluid, a new capability for multiphase RMT simulation.

This paper proceeds as follows. We first review the fundamental relations of each phase. The RMT methodology is then described, and we present the selected stencils for our finite differences. Next, we define the transition zone and how the fields are appropriately mixed within its band of influence. We then discuss a number of FSI test results and the convergence of the method. An important ingredient in our scheme is a new spatial extrapolation procedure that enforces certain physical and smoothness requirements in the field extrapolation. This is presented and shown to remedy deleterious extrapolation artifacts near the interface that can occur otherwise. We then describe how the routine can be appended with a solid–solid contact subroutine by taking advantage of the signed distance-measures inherent within the existing level set fields used in distinguishing the interfaces. Two examples are provided of submerged solid–solid contact.

where $J=detF$ and a prime indicates the deviatoric part of a tensor. Equations (2), (5), and (7)–(9) are a closed Eulerian system for computing solid deformations. As described in Ref. [16], this system of equations can be recast in conservative form and implemented numerically in a discrete conservation scheme should nonsmooth solutions be expected.

where η is the viscosity and λ is the compressibility modulus.

The computations are carried out on a two-dimensional m × n grid of square cells with side length h, using plane strain conditions. The interface between the fluid and solid phases is described using the level set method [4,25], whereby an auxiliary Eulerian function $φ(x,t)$ is introduced such that the interface is given by the zero contour, $φ=0$⁠. At the start of each simulation step, the level set function is initialized to be the signed distance to the interface, using the convention that $φ<0$ inside the solid region and $φ>0$ in the fluid region. While there are a number of different numerical approaches for the level set method, we make use of the specific implementation described by Rycroft and Gibou [18]. The implementation contains a procedure for reinitializing the level set function so that it is a signed distance function to the zero contour, using a combination of the modified Newton–Rapshon iteration of Chopp [26] to update values of $φ$ adjacent to the interface, and a second-order fast marching method [4] to update values further away.

Because we develop a blurred interface technique, we introduce a transition zone corresponding to the region where $|φ|<wT$⁠, where w_T is a constant (Fig. 1). As described in detail later, the material response in the transition zone is modeled as a mixture of fluid and solid phases. The region $0<φ<wT$ is more fluid than solid, and the region $wT<φ<0$ is more solid than fluid. For a positive constant w_E such that w_E > w_T, we define the extended solid domain to be the region $φ<wE$⁠, which corresponds to a region with any fraction of solid response together with a thin band of points adjacent to the transition zone. For computational efficiency, our level set method implementation only stores and updates the values of $φ$ in a narrow band of grid points within a distance ±w_E of the interface.

Figure 1 shows the discretization of the simulation fields. Based on considerations of the finite-difference stencils presented later, a staggered field arrangement is used. Some fields are held at cell corners, which are indexed using i, j for i = 0,…, m and j = 0,…, n. Some fields are held at cell centers, which are indexed using [i, j] for i = 0,…, m − 1 and j = 0,…, n − 1. The level set function $φ$⁠, velocity field v, fluid density ρ^f, and combined density ρ are stored globally at cell corners in both the fluid and solid domains. The fluid stress is stored at globally at cell centers. In addition, in the extended solid domain, the reference map field $ξ$ is stored at cell corners, and the solid stress $σs$ and density ρ^s are stored at cell centers.

We now describe the progression of the blurred interface method. Starting at a time step n, the algorithm builds the kinematic fields for the next step n + 1. The stress fields are calculated at every time step as a midstep calculation. The algorithm below summarizes the major steps of the routine. Details are provided thereafter.

We employ a tight stencil to calculate the divergence of stress and the (constitutive input) gradients of $ξ$ and v as in Fig. 2, with more details on our stencil selection to be discussed later. The discretized gradient of the reference map is given by

and $J[i,j]n=(det∇ξ[i,j]n)-1$ would be applied to create $σ[i,j]sn$ and $ρ[i,j]sn$⁠. When put together, the divergence of the solid stress at i, j is calculated from nearby values of $ξ$ at i, j, at i − 1, j, at i + 1, j, at i, j − 1, at i, j + 1, at i − 1, j + 1, and at i + 1, j − 1.

where the discretization of ∇v uses the same stencil as in Eq. (15).

where the advection term uses the ENO discretization, and the discretization of the fluid divergence is chosen in a similar manner to Eq. (15).

respectively, where $φ[i,j]$ is calculated by averaging the four values on the grid cell corners.

where the above uses the average of four values to interpolate ρ^s onto the cell-cornered grid.

The above mixing procedure models the no slip condition between fluid and solid. Other interfacial conditions are possible within this framework, though we have yet to perform extensive tests on them. For example, perfect slip may be achievable by applying a shear traction elimination step after the above mixing is applied, which modifies the stress in the transition zone by smoothly decreasing the shear traction resolved onto the level set contours such that shear traction vanishes at $φ=0$⁠.

In order to define $ξ$ and consequently $σs$ in the portion of the extended solid domain where $φ>0$⁠, great care must be taken. This zone represents a region whose behavior is more fluid than solid, and, thus, is liable to undergo extensive shearing and nonlinear deformation. This is problematic because large nonlinear deformations reduce the accuracy of computation of $∇ξ$ necessary to calculate $σs$⁠, and, moreover, the solid stress can become unboundedly large as deformation builds up in this zone, which negates the purpose of the gradual switch-over to a fluid-dominated stress response in Eq. (25). For these reasons, we avoid advecting the reference map per Eq. (20) in this region and instead obtain values for it through extrapolation from the nearby grid points where $φ<0$⁠. This constitutes the only field extrapolation in our current method, and it is recalculated every time step. Note that in the velocity update step, only those values of stress within the transition zone (i.e., very close to the $φ=0$ contour) influence the outcome, which supports our usage of near-field extrapolation as an appropriate remedy. Our algorithm for consistent extrapolation is left for a detailed discussion in Sec. 6.

The above comprises a general fluid–solid interaction routine for simulating solids with finite-strain constitutive relations coupled to a fluid phase that obeys the compressible Navier–Stokes equations. We now verify the method and its accuracy using several test geometries, many of which (those of Sec. 4.1) under the previous sharp-interface approach [16] exhibited interfacial errors of the type previously described, even with the inclusion of surface tension and significant solid viscosity. Our current blurred-interface routine has also reduced the number of extrapolation steps significantly—20 scalar fields are extrapolated per time step in the sharp-interface approach, but now only the two components of $ξ$ are now extrapolated. Unlike the sharp approach, we do not require a specific routine to ensure interfacial jump conditions as they are implicitly satisfied in the transition zone. These comments are not to suggest that the numerical interfacial errors of the sharp approach are irreparable in that framework; it is possible they could be alleviated in part through integration of our improved extrapolation algorithm of Sec. 6 among other adjustments, but we reserve this topic for future work.

Throughout this paper, we make use of dimensionless simulation units. In all tests shown, unless otherwise stated, the following input material parameters are used. The initial fluid density is $ρ0f=1.0$⁠, the viscosity is η = 0.12, and the compressibility is λ = 60.0. The solid uses G = 10 and κ = 50. The artificial viscosity is η_a = 0.012. These parameters match our previous work [16], though it should be noted the previous work also included a significant separately added solid viscosity. The predominant Courant–Friedrichs–Lewy (CFL) criterion is due to the fluid viscosity, which is proportional to h²ρ^f/2η. The other criteria that need to be considered are the one due to compressibility of the solid $~hρs/κ$⁠, the one due to compressibility of the fluid $~hρf/λ$⁠, and the convective CFL $~h/max|v|$⁠. For all of our results, we use a time step of Δt = 0.05h²ρ^f/η, which satisfies the dominant CFL condition for the fluid viscosity and also satisfies the additional CFL conditions for all of the grids and parameters that are used. Throughout, we use a transition zone width of w_T = 3h and we use w_E = 6h. The simulations are written in C++ and use the OpenMP library to multithread many of the operations that scan over the entire grid of points.

For all of the simulations presented here, domain boundary conditions are applied by fixing the values of v and ρ^f on the edge of the cell-cornered grid, where i = 0, j = 0, i = m, or j = n. We consider two types of boundary conditions: Dirichlet conditions where the boundary field values are fixed, and free boundary conditions where the boundary field values are linearly extrapolated from the two adjacent layers of interior points. For example, to implement a free boundary condition on an arbitrary field f_i,j at the top boundary, we set f_i,n = 2f_{i,n − 1} − f_{i,n − 2} for i = 0, 1,…, m. Note that a free condition on the fluid velocity component tangent to a boundary gives the perfect slip condition.

The first example we consider involves a bulky, highly deformable neo-Hookean rod. The rod is initially vertical and anchored at its top end. It is placed in a horizontal fluid flow, which causes it to deform significantly. The domain covers −2 ≤ x ≤ 2 and −2 ≤ y ≤ 2 using a square grid of size 240 × 240. The bar initially covers the rectangle −1.4 < x < 0.6 and $|y|<1.1$ and has semi-circular end caps. The anchored region is a circle with center (x, y) = (−1, 1.1) and radius r_a. During the simulation, the reference map in the anchored region is enforced to be constant and the velocity is enforced to be zero. We have simulated two different anchor sizes r_a = 0.25 (Fig. 3) and r_a = 0.15 (Fig. 4). The fluid inflow and outflow is controlled by applying a constant horizontal velocity of (u, v) = (0.24, 0) on the left and right sides of the domain window. Perfect slip boundary conditions are used on the top and bottom sides, where v = 0, and u and ρ^f are free. In each simulation, the fluid flow deforms the rod inducing large local stretches in the solid, with stretch ratios exceeding two in parts of the bar. In the case of the smaller anchor, the bar is less able to resist the incoming fluid flow and swivels out of the way to a greater extent as expected. The simulations each reach a steady state by t ≈ 20 where the rod remains in a static, bent configuration with steady fluid flow surrounding it. On an Apple MacBook Pro (Mid 2014) system with a quad-core 2.8 GHz Intel i7 processor, the simulation using r_a = 0.25 takes 997 s using 216,090 time steps using a single thread. If two, three, or four threads are used, the simulation takes 773 s, 705 s, and 635 s, respectively. Because this simulation only uses a small grid, the speedup from multithreading is only modest, since the overhead from creating threads is comparable the computational work done each time step. However, for some of the larger simulation grids considered later, multithreading becomes significantly more advantageous.

We also used this configuration with r_a = 0.25 to test the speed of the simulation method against the previous sharp interface method [16]. For a benchmark test simulating over the interval 0 ≤ t ≤ 2.5 with 3,200 time steps using a single thread, the sharp interface method took 94.01 s, while the new method took 53.97 s, corresponding to a speedup factor of 1.74. These two simulations were not perfectly comparable, since the memory organization and boundary implementation are slightly different. However, the results demonstrate that the present method is significantly faster, due to the simplification of the finite-difference stencils, a reduction in the number of simulation fields required, and fewer boundary extrapolations.

The second test case we consider involves a rotor of neo-Hookean material anchored around a frictionless pivot at its center. It is deformed and ultimately set into steady rotation by incoming fluid. Each blade of the rotor comprised a rectangle of length 1.1 with a semi-circular end cap. The join between each pair of blades is smoothed with a quarter-circle with radius of 0.4. The fluid boundary condition consists of a prescribed fluid inflow of (u, v) = (0.2, 0) on the bottom half of the left edge of the computational domain and an outflow of (u, v) = (0, 0.2) in the top right half edge of the domain; at all other boundaries, perfect slip boundary conditions are used. The rotor and fluid have the same material parameters as Sec. 4. This particular geometry, when implemented under the previous sharp-interface method, displayed many erroneous pressure shots about the interface, as the particular shape involved has much curvature variation. As can be seen from Fig. 5, this artifact is nonexistent in the current scheme. By t = 14, the rotor motion and fluid flow approach a steady repeating cycle, with each cycle corresponding to one quarter-turn of the rotor. By the final snapshot in Fig. 5 at t = 250, the rotor has undergone approximately $5.5$ complete rotations without any interfacial perturbations. Figure 6 confirms that the rotor enters a steady repeating cycle, by showing the rotation angle θ_p and angular velocity ω_p of the pivot as a function of time. Figure 6(c) shows a plot of θ_p against ω_p, where after an initial transient period as the rotor begins to turn, a steady relationship between ω_p and θ_p develops, with the angular velocity being slightly higher during periods when a blade of the rotor is directly in front of the region of fluid inflow.

where $Tk=1/2$ if k = 0 or k = n and T_k = 1 otherwise, so that the sum calculates the trapezoidal rule. By interpreting the $|·|2$ operator appropriately, Eq. (27) can be applied to scalars, vectors, or tensors. Since the resolution of each smaller grid is an exact multiple of the reference grid, each grid point (i, j) on a smaller grid exactly coincides with a grid point $(i',j')$ on the 1800 × 1800 grid, and hence we take $fi,jexact=fi',j'$⁠. For an arbitrary cell-centered field f, the same definition is used but replacing i, j with [i, j]. In this case, the smaller grid points may not coincide with the 1800 × 1800 grid points, in which case we calculate f^exact using bilinear interpolation. Figure 7 displays the convergence properties of four key fields. The velocity convergence plot in Fig. 7(a) compares the velocity fields in the entire domain to that of the reference solution. For the convergence plots of $ξ$ in Fig. 7(b), the sum in Eq. (27) is only evaluated at grid points where $φ<0$ in both the smaller grid and the reference grid. The figure indicates that all fields are converging with at least a linear rate, as we expect for the chosen stencils.

Figures 7(c) and 7(d) show the convergence of the in-plane pressure p and the in-plane effective shear stress $q=|σ1-σ3|/2$ where σ₁ and σ₃ are, respectively, the maximum and minimum principal stresses. The norms are split into contributions from solid and fluid phases, determined by the sign of $φ$ on the coarse grid. In both phases, the stress converges at a similar rate to the velocity and reference map. For the pressure, the fluid phase has the larger errors, due to the high value of λ. However, for the effective shear stress, the solid phase has larger errors, suggesting that errors in solid shear stress are more significant than errors in fluid viscous stress.

It is interesting to zoom into a region containing the interface and compare simulations with a coarse grid to a more refined one. Figure 8 shows a comparison of the pressure fields at t = 4 in the upper right concave part of the rotor with the 200 × 200 grid and the 600 × 600 grid. As the grid size decreases, and likewise the transition zone width, the interfacial behavior approaches that of a sharp interface; to wit, the pressure field develops a rapid variation across the interface that approaches a discontinuity. Analytically, the stress component representing tension in the direction tangent to the interface need not be continuous across the interface, which is why the pressure field should adopt this feature in the sharp limit.

With basic demonstrations in hand, we now return to discuss the details of our reference map extrapolation step. Physically consistent extrapolation of $ξ$ is key to a proper interface representation, as we shall show. This section highlights two new subroutines we have developed for this purpose.

The extrapolation algorithm we use is based on that of Aslam [30], which takes advantage of the regularized level set $φ$ already stored on the grid. This routine has the benefit of providing linear extrapolation outside an arbitrary-shaped domain, by extrapolating fields in the outward-normal direction to the interface—the outward normal is obtained easily via $n∧=∇φ$⁠. An immediate requirement of such a routine is that the gridwise field being extrapolated be smooth near the edge of its known domain. Numerical oscillation error in the known domain, hence, invalidates the extrapolation procedure. This motivates our use of one-sided difference stencils in Sec. 3.1, which do not support odd–even oscillation error. Future work will explore higher-order, nonoscillatory stencils to achieve this same aim akin to the usage of staggered grids in computational fluid dynamics.

While its speed and geometric generality are benefits of the Aslam method, we and others [18] have noticed that the hyperbolic nature of the extrapolation routine can cause extrapolated fields to develop striations—loss of smoothness in the direction parallel to the interface—even when the known domain data are smooth. If uncorrected, striations in the extrapolated values of $ξ$ cause oscillations in the corresponding $σs$⁠, which induce an artificial wrinkling motion that becomes apparent in the $φ=0$ level set (see Fig. 9). This phenomenon destabilizes the routine if the wrinkle curvature grows to a level that competes with the grid spacing.

We note that satisfaction of the above constraint is equivalent to ensuring that all material points initially within the solid remain within the solid.

Because the values of $ξ$ and $φ$ are obtained from different numerical routines, discretization error can cause the above constraint to be violated over time. Recall that $ξ$ is generated on points with $φ>0$ solely through linear extrapolation from the adjacent $φ<0$ domain, and hence the values of ψ near the zero contour of $φ$ are not more than first-order accurate. The comparison between them in Fig. 10 is indicative that a drift between ψ and $φ$ grows after the onset of the previously described wrinkling artifact. This in turn causes persistence of the wrinkled shape, because some material that began on the solid side of the $φ=0$ interface is falsely assigned as fluid, which permits the wrinkles to sustain themselves over time without the elastic response correcting their shape. This calls for a routine to ensure the physical requirement that the zero contours of $φ$ and ψ satisfy Eq. (28), i.e., remain “pinned,” to ensure that solid stays solid and fluid stays fluid. One could argue at this moment that $φ$ is a redundant field, in that ψ could operate just as well on its own at discerning phases. However, we recall that $φ$ also gives a measure of distance, as its gradient always has magnitude 1. This feature is exploited in the Aslam extrapolation scheme and is needed for the solid–solid contact algorithm of Sec. 6.

We now describe two midstep subroutines, which are applied after an initial Aslam extrapolation of $ξ$ and successfully resolve the interfacial wrinkling, phase exchange error, and instabilities relating to these two phenomena. The need to appropriately pin the reference map extrapolation field and the level set field was pointed out previously [16] but a general and accurate pinning technique was not employed in that work. Artificial surface tension, projection of the extrapolated reference map field, and solid damping were used to reduce the appearance of these problems in that work.

The first step in our procedure is to remove artificial striations in the extrapolated values, see Algorithm 2. It is a general algorithm that can be applied to reduce fluctuations in an arbitrary extrapolated field f. In the FSI routine, it is applied separately to ξ_x and ξ_y. The fields on the solid side of the interface are never affected; only the extrapolated values within the fluid domain are adjusted. The algorithm consecutively applies a diffusion stencil to improve the extrapolated values. The edge of the solid domain is sampled by the diffusion stencil and hence the resulting extrapolation is both smooth within the fluid domain and continuously extends data from the solid domain.

Upon completion of Algorithm 2, the second midstep routine starts by updating the level set values $φ(x,t)$ in the narrow band to be equal to $ψ(x,t)$⁠. This ensures that the zero contour of the level set is consistent with the reference map, but the new values of $φ(x,t)$ may not satisfy the signed distance function property, $|∇φ|=1$⁠. We therefore make use of the reinitialization routine described in Ref. [18], to rebuild the narrow-banded level set function and recover this property. Figure 11 displays the solid stress from the reference map and its extrapolation in the same test geometry but using our new two-part extrapolation-pinning routine. The formerly observed striated solid stress, interfacial wrinkles, and drift between fields have been eliminated.

We now consider a second solid phase and seek to simulate a fully coupled fluid–solid–solid interaction, solids interacting through contact while submerged in a fluid. We describe the procedure for two solid phases, but the procedure could be generalized to more solids. We use superscripts (1) and (2) to denote each solid phase. There are now two level sets, $φi,j(1)$ and $φi,j(2)$⁠, to describe each solid's boundary and measure distance to it. Each level set has a corresponding transition field, Ω⁽¹⁾ and Ω⁽²⁾, using the same w_T for both. We maintain two reference map fields, $ξi,j(1)$ and $ξi,j(2)$⁠, and corresponding solid stresses, $σ[i,j]s(1)$ and $σ[i,j]s(2)$⁠, within the extended domain of each associated solid. The velocity field, as before, is a single unique field for the whole computational domain, independent of phase. The main procedure, Algorithm 1, requires very little change to represent interacting submerged solids. All steps corresponding to the solid are now simultaneously performed on both sets of solid fields. The major changes we must prescribe are to adjust the mixing rule in step 8 to correctly cross-fade between the two solids and fluid and to add an extra routine to correctly set contact conditions. The latter consideration ultimately arises as an adjustment in Step 6.

On each cell-centered grid point, the quantity of each solid follows the previous behavior, and the remainder is filled with fluid. The coefficient of the third term remains non-negative so long as the solid phases do not penetrate. The same approach is applied in defining the density.

The contact algorithm focuses on the case of frictionless contact, but it is sufficiently general that other contact conditions could be implemented. In the blurred interface routine, we define contact between two objects whenever there is an overlap between the transition zone of one object and the interface of the other. We define penetration to occur if the zero contours of each level set pass through one another. Hence, a nonpenetrating contact routine permits the zero level sets of the two bodies to be separated by less than half the transition zone width—the “start” of blurred contact—and rapidly penalizes any closer approach of the two solids. When the solids are farther apart than half the transition zone width, the routine should have no effect. As the grid size shrinks, this description approaches the standard definition of nonpenetrating contact between sharp interfaces.

which is used within Eq. (21), where $n∧12i,j=±∇φ12i,j/|∇φ12i,j|$ is the vector normal to the contours of $φ12$⁠, pointing away from the midsurface. The prefactor k_rep must be chosen large enough to successfully repel a crossing of the interface-centers, but must be small enough that the stable time-increment of such an explicit contact routine also decreases with k_rep. In the following examples, we choose $krep=1/20$⁠.

Figure 13 shows six snapshots of pressure during the simulation. By t = 9, the pressure starts to build up between the two disks as they come into contact. By the time the moving disk reaches its leftmost position at t = 12.5 when x_d(t) = 0.1, the pressure has increased further. At t = 18, the moving disk has separated from the anchored disk, and a small region of negative pressure is created between the two disks as fluid. At t = 25, the moving disk comes to rest at its original position, and returns to its original undeformed circular shape.

Figure 12 displays the γ_i,j function at three time points and shows how the solids do not feel the interaction until the moment when the narrow δ_s function enters into the solids. Within the simulation, the function is only defined in the region where the two narrow bands of the level set overlap. The repulsion force succeeds in separating the two interface-centers but maintains contact in the sense defined above. Some small asymmetry is visible, which is expected due to the asymmetry in the stencils described in Fig. 2; this effect diminishes with grid spacing. As is evident in Fig. 13 using this technique, we are now capable of getting very large deformations in the solids through contact, without any problem of penetration or sticking of one solid onto the other.

The second example that we consider is a disk bouncing on an anchored rubber bar, all while submerged in fluid. The disk has radius of 0.4 and is initially centered at (x, y) = (1, 0). The bar comprises the rectangle −1 < x < 1, −0.9 < y < −0.1 along with semi-circular end caps and is anchored in a circle with radius of 0.2 at (x, y) = (−1, −0.5). For this example, the density of the disk is 20, the disk's initial downward velocity is −1, the density of the rod is 6, the gravity is 0.03, and the fluid viscosity is η = 0.04. All other parameters are unchanged from the previous results. A grid of 600 × 600 is used and the system is simulated over the range of 0 ≤ t ≤ 25.

Figure 14 shows nine snapshots of the in-plane pressure field during this simulation. The disk first reaches the rod at approximately t = 1.5, and a large region of positive pressure is visible between the two. By t = 2.5, the force exterted by the heavy disk deforms the bar into a U-shape, with a large tension visible on the bottom side of the bar. By t = 4.5, the right end of the bar has moved downward, and this motion pushes the disk upward so that it undergoes a bounce. By t = 7, the disk is fully separated from the bar. The disk then sinks and comes into contact with the bar at t = 15, before slowly sliding down the bar.

This approach allows us to model frictionless nonsticking contact that avoids cross penetration of the solid phases. We could imagine adherence conditions, friction, or other contact laws through new definitions of a possibly evolutionary influence function centered at the midsurface. Other Eulerian contact approaches also exist; we can, for example, measure the amount of overlap a step would cause and then apply a correction based on the intersecting volume [31] and a minimization problem finding the optimal velocity field which solves the equilibrium condition while minimizing overlap.

This work has described a blurred-interface finite-difference method for FSI on a single Eulerian grid. The method is notable for its simplicity and speed. Our explicit algorithm invokes a computation of fluid and solid stress, which are then mixed in accordance with a transition function, as are the fluid and solid density. The solid stress and density are computed in Eulerian-frame with the aid of the reference map field, which is stored and updated throughout. An important pair of subroutines is also presented, to improve the representation of the reference map and level set fields near the interface. Our method is shown to converge and as grid-size decreases we recover results that properly display the signatures of a sharp interface. The framework we create extends to the case of solid–solid contact as well and we have given two examples of submerged contact for the case of nonpenetrating, frictionless, nonsticking behavior. Here, the key idea is to take advantage of the distance-function property of the level set fields about each solid to produce a short-range separation force that acts only within the solid interfacial regions when they are close enough to each other. We note that while our method has the key feature of enabling simulation of largely deforming soft solids, there is nothing preventing this approach in the case of stiff solids; one can choose the elastic constants at will, as long as the stable time step is selected correspondingly.

There are a number of important directions to consider from this point. From a numerics standpoint, we can go to higher convergence rates by moving to a higher-order set of finite-difference stencils, though we emphasize the importance of removing oscillatory numerical errors, as we discussed in Sec. 6, which can arise in higher-order stencils. It is also important to port our scheme into 3D, which should require very little algorithmic adjustment. We can also implement incompressibility constraints by adapting the projection method [32]. It would be worthwhile to parallelize the scheme beyond multithreading to take advantage of a distributed-memory architecture. While we have focused on Cartesian grids for ease, this is not a necessity, and it would be useful to consider general meshes that enable local refinement.

There are also several next steps to be taken from the applications standpoint. Because our solid formulation is rooted in finite-deformation mechanics, any thermodynamically consistent constitutive law should be within the range of simulation. For example, one could model solid behaviors that depend on evolving internal variables by expressing the evolution laws of said variables in Eulerian-frame and storing/updating those variables on the grid (and using our various extrapolation routines to maintain accuracy near solid boundaries). We are also actively pursuing submodeling capabilities within our approach, which model the domain using two different grids of different refinement levels, and enable us to refine a local zone without the need to decrease the time step in the coarse grid zone. This would enable local refinement near a rough fluid–solid interface and could be applied to multiscale problems such as nuclear fuel-rod fretting, which involve interactions between cooling fluid and roughened solid parts. To extend the solid contact routine to the case of many submerged solid phases could have considerable usage in modeling dense suspensions of soft particles. To be able to model contact-induced finite deformations of particles is crucial when the compliance of particles is too high to admit simple contact-force laws. Biological applications, particular at the cellular level, could offer another arena in which this approach could be of use, as the solid components are highly deformable and fluid permeates the system. The Eulerian form could be advantageous in the implementation of growth models [33], or simultaneous species diffusion, which is also amenable to Eulerian-frame.

B. V. and K. K. acknowledge support from the MIT Department of Mechanical Engineering. K. K. acknowledges support from the Consortium for Advanced Simulation of Lightwater Reactors (CASL), an Energy Innovation Hub for Modeling and Simulation of Nuclear Reactors under U.S. Department of Energy Contract No. DE-AC05-00OR22725. C. H. R. was supported by the Director, Office of Science, Computational and Technology Research, U.S. Department of Energy under Contract No. DE-AC02-05CH11231. We thank the anonymous reviewers for their detailed comments and feedback.