Fluid-Solid Coupling in Kinetic Two-Phase Flow Simulation

Real-life flows exhibit complex and visually appealing behaviors such as bubbling, splashing, glugging and wetting that simulation techniques in graphics have attempted to capture for years. While early approaches were not capable of reproducing multiphase flow phenomena due to their excessive numerical viscosity and low accuracy, kinetic solvers based on the lattice Boltzmann method have recently demonstrated the ability to simulate water-air interaction at high Reynolds numbers in a massively-parallel fashion. However, robust and accurate handling of fluid-solid coupling has remained elusive: be it for CG or CFD solvers, as soon as the motion of immersed objects is too fast or too sudden, pressures near boundaries and interfacial forces exhibit spurious oscillations leading to blowups. Built upon a phase-field and velocity-distribution based lattice-Boltzmann solver for multiphase flows, this paper spells out a series of numerical improvements in momentum exchange, interfacial forces, and two-way coupling to drastically reduce these typical artifacts, thus significantly expanding the types of fluid-solid coupling that we can efficiently simulate. We highlight the numerical benefits of our solver through various challenging simulation results, including comparisons to previous work and real footage.


INTRODUCTION
Fluid flows are all around us in our daily lives.Their complex dynamics and visual appeal are a source of awe for many, and an enduring challenge for computer graphics (CG) and computational fluid dynamics (CFD) researchers.While early fluid simulation offered only limited realism, the ever-more sophisticated techniques used in graphics for the animation of immiscible multiphase flows have successfully demonstrated the ability to produce typical complex behaviors such as bubbling or splashing realistically.To achieve this feat, various physically based simulation methods based on grids and/or particles have been developed to capture both the individual phases and their detailed interactions.More recently, the use of kinetic solvers based on lattice Boltzmann methods (LBM) have proved capable to handle large air-to-liquid density ratios and generate very turbulent flows, which had been out of reach for Navier-Stokes solvers in graphics.Yet, even if state-of-the-art fluid simulation of complex water-air interactions has reached cinematic realism at very reasonable computational cost, fluid-solid coupling is still limited to low-speed motions, preventing the simulation of common natural occurrences like the thrust of a propeller underwater or a stone skipping over a lake.
Reasons for the lack of robust and efficient coupling in multiphase flows are manifold.First, most graphics fluid simulators achieve efficiency through large time steps; although the resulting lack of accuracy is acceptable for visual purposes, large time steps are simply unable to capture the subtle interactions between two phases or the coupling between phases and obstacles, particularly for fast motions and turbulent flows.Second, enforcing a strong coupling (via a monolithic solver) is inefficient for turbulent multiphase flows, so partitioned coupling (for which a fluid solver and a deformable body solver are run separately, yet interfaced so as to exchange forces) is often favored in CG; but it often results in oscillating pressure near boundaries and inaccurate interfacial forces which eventually lead to blow-ups.Even the use of the immersed boundary method in CFD is known to fail for fairly turbulent flows and when dealing with thin shell models due to its inherent leakage through boundaries, making it poorly adapted to CG.
In this paper, we extend the phase-field and velocity-distribution based LBM approach of [Li et al. 2022] to offer two-way coupling and significantly improved robustness.We show that more accurate interface normal evaluations and phase field gradients dramatically ACM Trans. Graph.,Vol. 42,No. 4  enhance stability of fluid-solid coupling.By using a twice-finer phase-field discretization than the grid size of the fluid distribution function, we also significantly improve boundary handling and offer a simple strategy for dealing with dead and fresh cells.Our contributions bring robustness to multiphase flow simulation without impairing accuracy or efficiency, allowing us to simulate a far improved range of motions as shown for example in Figs. 2, 8, and 17.We also demonstrate that our solver reproduces well-known coupling behaviors, such as the fluttering or tumbling of coins in water (Fig. 16), or stone skipping over an air-water interface (Fig. 19) as it can handle the real-life density ratio between air and water as well as high Reynolds numbers flows very efficiently due to its massively-parallel implementation.

PREVIOUS WORK
Dealing with fluid-structure interactions has drawn significant attention in both CG and CFD due to its key role in capturing the fine interplay between obstacles and flows.We briefly review previous work on coupling of single-and two-phase fluid and solids.

Single-phase and free-surface solvers
Single-phase fluid-solid coupling.Early works handled fluid-solid coupling explicitly, semi-implicitly or fully implicitly [Yngve et al. 2000;Foster and Fedkiw 2001;Carlson et al. 2004;Takahashi et al. 2002;Génevaux et al. 2003;Klingner et al. 2006;Chentanez et al. 2006;Batty et al. 2007] by taking the Lagrangian solid velocity as a Neumann boundary condition for the fluid and integrating pressure forces at the solid surface.A fully Eulerian fluid-solid treatment was later introduced in [Teng et al. 2016].To deal with thin-shell solids, [Guendelman et al. 2005] proposed ray-casting to prevent leakage, while [Robinson-Mosher et al. 2008, 2009] improved on this approach using ghost cells and more accurate momentum transfer.Eulerian solvers with boundary-conforming meshes [Feldman et  2005a,b; Dai and Schmidt 2005;Klingner et al. 2006;Elcott et al. 2007] have also been suggested, followed by Lagrangian methods [Misztal et al. 2010;Clausen et al. 2013].More recently, the use of cut-cells using sub-cell information [Roble et al. 2005;Batty et al. 2007;Ng et al. 2009;Gibou and Min 2012;Weber et al. 2015;Edwards and Bridson 2014;Liu et al. 2015;Azevedo et al. 2016;Tao et al. 2022] has found wide adoption due to its reduced computational complexity and absence of fluid leakage for thin objects.
Free-surface fluid-solid coupling.A Lagrangian discretization of fluids through smoothed particle hydrodynamics (SPH) was first proposed by [Desbrun and Gascuel 1996;Müller et al. 2003], but its inherent visual blobbiness near the interface started a number of palliative measures [Desbrun and Cani-Gascuel 1998;Zhu and Bridson 2005].After the incompressibility of SPH was improved in [Solenthaler and Pajarola 2009;Bender and Koschier 2016], freesurface and solid boundary conditions were handled for particles more generally in [Schechter and Bridson 2012;Koschier and Bender 2017;Band et al. 2018], sometimes through penalty forces [Peer et al. 2015;Ihmsen et al. 2013].To generate uniform particle distributions and precise volume control, power particles [de Goes et al. 2015] and hybrid particle-grid transfers [Qu et al. 2022] were also introduced.Hybrid particle-grid methods have been favored in recent years for coupling [Fei et al. 2018;Hu et al. 2018;Fei et al. 2019;Fang et al. 2020;Fei et al. 2021], but they are limited to fairly viscous flows.A mix of FLIP and boundary element method (BEM) [Huang et al. 2021] was even formulated to simulate larger-scale coupling, while Shao et al. [2022] presented an unsmoothed aggregation algebraic multigrid (UAAMG) method to support large-scale simulation, but bubbles and two-way coupling were not handled.The use of levelset-based free-surface liquid and solid one-way coupling has also demonstrated success when used with adaptive octrees for high-resolution surface discretizations [Aanjaneya et al. 2017;Ando and Batty 2020] or cut-cells to handle static thin objects [Chen et al. 2020b].The volume-of-fluid (VOF) representation was also leveraged to design monolithic solvers [Takahashi andBatty 2020, 2022] for fluid-solid coupling, but such solvers have not been extended yet to multiphase flows and are currently rather inefficient.The lack of versatile and fast coupling techniques have also spawned a number of specialized solvers: [Ruan et al. 2021] proposed a three-way coupling method via a hybrid Eulerian-Lagrangian framework to simulate large surface tension flow in which liquid and solid are coupled by a Lagrangian membrane; [Xiong et al. 2022] proposed a Clebsch method to simulate free-surface vortical flow; [Liu et al. 2022] introduced a Lagrangian model for viscous solid-fluid interaction; and position-based dynamics (PBD) for surface tension [Xing et al. 2022] was even proposed, but it does not support turbulent multiphase flows.

Navier-Stokes-based multiphase fluids
While ad hoc techniques proposed simply to add bubbles in the flow [Thuerey et al. 2007;Kim et al. 2010;Goldade et al. 2020], SPHbased methods [Solenthaler and Pajarola 2008;Ren et al. 2014;Yan et al. 2016;Yang et al. 2017] were shown capable of simulating multiple immiscible fluids of different densities, but only for viscous flows.Multiphase flows with a levelset [Kim 2010] or VOF [Cho and Ko 2013;Langlois et al. 2016] encoding of the interface have improved their treatment of discontinuous jumps in fluid density and pressure at the interface [Boyd and Bridson 2012;Kim et al. 2007;Losasso et al. 2006;Mihalef et al. 2006] through a variable density pressure projection [Kang et al. 2000] or the ghost fluid method [Hong and Kim 2005]; but again, most approaches along these lines focus on viscous fluids.There are also a few approaches relying on a diffuse-interface model [Song et al. 2005;Zheng et al. 2009] to avoid having to deal with a sharp interface.Hybrid grid-particle methods [Saye 2016[Saye , 2017;;Yang et al. 2021;Su et al. 2021] have also been used for multiphase flows.A Moving Least Square Reproducing Kernel Method (MLSRK) framework to simulate multiphase fluids and solids in a unified manner was also formulated [Chen et al. 2020a]; however, the need for small time steps renders the solver impractical for large scenes.A mesh-based Lagrangian approach to simulate multiphase flows was also proposed in [Misztal et al. 2013], but only examples of very low Reynolds numbers were demonstrated.

Kinetic fluid solvers
Recently, kinetic solvers based on statistical mechanics have been introduced in graphics [Guo et al. 2017;Li et al. 2019], building upon the significant improvements over the original lattice Boltzmann method (LBM) used in computational physics over the last ten years.LBM-based kinetic solvers have emerged as a great alternative to the traditional incompressible Navier-Stokes fluid integrators from graphics both for single- [Li et al. 2020;Lyu et al. 2021] and multiphase flows [Li et al. 2021[Li et al. , 2022] ] due, by and large, to their massively-parallel nature: their high computational efficiency when implemented on GPU [Chen et al. 2021] is extremely attractive, offering the computation of realistic flows in a matter of minutes on a single GPU; moreover, progress in LBM collision modeling has resulted in simulation accuracy superior to Navier-Stokes solvers.While coupling in single-phase incompressible flows has already reached a maturity which is arguably sufficient for many graphics applications [Li et al. 2020;Lyu et al. 2021], this is far from true in the case of multiphase flows.In fact, developing a numerical scheme able to simulate the type of high Reynolds numbers and large density ratios that water-air interactions require has always been a challenge Ð even in CFD where the use of the immersed boundary method [Xiao et al. 2022;De Rosis and Tafuni 2022] can only handle viscous multiphase flows.The first fluid solver to reach the desired generality for graphics was presented recently in [Li et al. 2022]; it also offered much improved stability in handling static boundaries by proposing the use of a velocity-based distribution function instead of a momentum-based one, to avoid jumps between the two fluid phases of different densities.However, the authors lamented that their contributions to łensure robustness for static obstacles do not suffice in the case of fast-evolving solid objectsž.This paper precisely addresses this limitation, in order to reach a far improved range of coupling in multiphase flow simulation.

OUR MULTIPHASE FLUID SOLVER
In this section, we recap the work of [Li et al. 2022] on which we base our solver, and detail the key changes we propose in order to significantly improve the numerical stability of fluid-solid coupling.

Fluid and phase-field equations
Model.The momentum and continuity equations for incompressible multiphase flows are often described macroscopically as: where  and u are respectively the spatially-varying fluid density and velocity of the two fluids,  is the hydrodynamic pressure enforcing incompressibility (Eq.( 1c)), Π is the viscous stress tensor, while F  and F  are body and surface tension forces.To replace the typically sharp encoding of the interface and associated boundary conditions at this interface, a diffuse-interface model, called phase field, tracks the interface between the two phases through a conservative, profile-preserving Allen-Cahn equation to replace Eq. ( 1a) as explained in [Li et al. 2022]: Due to the right-hand side term of this equation, the phase field maintains a predefined steep profile where the interface between the two fluids lies, from which the interfacial forces can be well approximated.For a mixture of two immiscible, incompressible fluids with respective density   and   and respective viscosity   and   (where the subscripts  and  are used to differentiate the high-density fluid from the low-density fluid), the density field  is expressed from the phase field via: (3) while the viscosity  is defined from phase viscosities   and   as: Lattice Boltzmann formulation.As described in [Li et al. 2022], the momentum and phase field equations for multiphase flows can be formulated through the time evolution of two distribution functions.The first one, a velocity-based distribution function g, encodes the flow in time through: where   represents g in the direction c  from the 327 lattice structure shown in Fig. 3(left), while Ω   and   account for the collision and forcing terms respectively in this direction.Solving Eq. ( 5) is done in three steps: streaming, collision and boundary handling.The collision term requires the evaluation of the equilibrium distribution g eq , which, for this velocity-based formulation, is: where    are the usual LBM lattice weights (see Fig. 3).The macroscopic hydrodynamic variables such as the normalized pressure  * Fig. 6.Hybrid moving bounce-back.In the streaming step, we accommodate fluid-structure coupling through bounce-back: the distribution function at a node x  whose link   ′ hits a solid boundary at an intersection point x  is made to bounce-back towards the opposite direction , along with a momentum exchange due to the solid's speed.
or the velocity u can be obtained respectively from the zeroth and first moments of the distribution function g.The second distribution function, h, discretizes the phase field equation as: where ℎ  is the distribution function in the direction d  from the 37 lattice structure (Fig. 3, right).Solving Eq. ( 7) is done through the same three-prong (streaming, collision and boundary handling) process, and the macroscopic phase field can be recovered at any time step  as the zeroth-order moment of h, i.e., While we adopt most of this framework laid out by [Li et al. 2022] for multiphase flows, we present a series of targeted improvements on boundary handling and a new, finer-resolution scheme for the phase field equation to significantly improve the treatment and robustness of fluid-solid coupling as we detail next.

Fluid-solid coupling
The pros and cons of the various approaches to handle boundaries are well documented by now [Lyu et al. 2021], with the simplicity of the moving bounce-back (MBB) scheme often winning over the  inability of the immersed boundary method (IBM) to handle thin objects in graphics applications.We also adopt an MBB-like approach in our work to allow for arbitrary objects; but before introducing what we call a hybrid moving bounce-back (HMBB) scheme, let's review the original moving bounce-back scheme for single-phase flows first.Since streaming simply updates the distribution function g as if the flow was locally unobstructed, one needs to alter this basic streaming for links intersecting a boundary.For a node at x  with a link   ′ crossing a static obstacle as shown in Fig. 6, we can simply redirect the streaming along the opposite link   , which reproduces a bounce of the fluid particles off the obstacle.More generally, if the obstacle is moving, the moving bounce-back scheme proposed by [Ladd 1994] is written for a momentum-based distribution as ) where   is the obstacle velocity at the link intersection   : it corresponds to a bounce of the distribution on the obstacle to which is added a momentum transfer due to the boundary motion.In our context of velocity-based LBM multiphase flow simulation, this approach is simply rewritten as: (10) This moving bounce-back approach no longer needs to consider the density discontinuity near the interface due to the use of velocitybased distributions Ð but it inherits the typical ring artifacts mentioned in [Lyu et al. 2021;Li et al. 2022] for high Reynolds numbers.To ensure numerical stability of our boundary treatment, we thus adopt the hybrid bounce-back scheme designed by [Li et al. 2022] for static obstacles, to which we add the momentum exchange when the obstacle is moving.Moreover, their approach relied on an equilibrium distribution based on the normalized pressure of a fluid node next to the obstacle, as this was deemed to be the best local approximation of pressure locally available.However, bounce-back often creates spurious pressure oscillations near solids, which renders the equilibrium distribution inaccurate.We thus employ a smoothed pressure formulated in [Sitompul and Aoki 2019] with a second order filter instead, to further improve the stability of our approach.In fine, our hybrid moving bounce-back leverages the equilibrium distribution function f (, p * ) coming from an obstacle boundary towards a fluid node along a boundary-crossing link based on the known solid velocity   and the most reliable normalized pressure available p * (x  ), resulting in a boundary update: where the filtered normalized pressure to enhance the stability is expressed as [Sitompul and Aoki 2019] Note that we set the blending coefficient in Eq. ( 11) to   = 0.1 in order to safely suppress pressure fluctuations near moving obstacles.This differs from [Li et al. 2022] which used a value of   based on pressure instead; our experience shows that such a strategy can in fact back-fire during impact due to pressure inaccuracy, so we opt for a safer pressure-independent value.

Resulting Forces on Solids
While our hybrid treatment of bounce-back accounts for the momentum exchange from solid to fluid, we also need to evaluate the exchanged momentum Δp from the fluid to the obstacle based on the distribution g.While many expressions have been proposed to evaluate this term with various degrees of accuracy, the kinetic-inspired approach of [Peng et al. 2016] which simply tracks momentum changes of fluid lattice particles with the fluidśsolid boundary during bounce-back is particularly simple to evaluate: when a link at x  is intersected by a boundary, the momentum transferred to the solid in the direction   ′ is  (x  , )   ′ (x  , ) +   (x  ,  + 1)   ′ .However, the simplicity of this evaluation is also its weakness: because it uses the local density  near boundaries, it is prone to inaccuracies: since LBM deals with weakly-compressible fluids, pressure waves near boundaries often lead to transient, yet exaggerated density changes.Our numerical tests indicate up to 20% density fluctuation on large impact of obstacles in water Ð e.g., in the cup drop example of Fig. 2. Therefore, we modify the expressions above by replacing the actual density  with a filtered density ρ computed as with 25   = 0 otherwise; i.e., we primarily take the average of the densities around x exactly as done for the pressure in Eq. ( 12), unless the density jump along a link is too significant (>25%, which could not be due to spurious oscillations alone) in which case we use the safer value at x  .Therefore, to obtain the force acting on the solid surface, one can sum up all these momentum changes to find the net loss of fluid momentum: where  x refers to the set of links (if any) starting at node x which intersect a boundary.Similarly, if x  denotes the barycenter of the obstacle, the total torque is written as: (15) Note that the force and torque expressions are in LBM space, so we need to map them to rigid-body physical space before sending them to the solid integrator, as detailed in [Li et al. 2020;Lyu et al. 2021].

Finer resolution for phase field equation
Existing works in LBM-based multiphase flows [Li et al. 2021[Li et al. , 2022] ] use the same LBM grid resolution for both flow and phase fields, albeit with a higher lattice count for the velocity.Since coupling between the incompressible flow and phase field equations happens through the macroscopic velocity u and the macroscopic phase field  and its gradient, our tests show that it is often a lack of accuracy in the interface normal that creates blowups.In this section, we thus detail how we employ two different grid resolutions (the phase field one being finer to improve interfacial accuracy) that we position in space in a staggered fashion for easy interfacing.
Grid ratio.Instead of implementing an adaptive grid refinement scheme to improve the resolution of the phase field, we simply decouple the flow and phase field resolutions: the distribution g is computed on a grid of size  3 , while the distribution h is computed on a ( ) 3 grid with a ratio  ∈ N.While any scaling  could be accommodated, we found that  = 2 is sufficient for graphics purposes and better in terms of cache efficiency: the 8-times finer spatial discretization for the phase field is enough to put velocity and phase field on equal footing in terms of accuracy.We will thus stick to this grid ratio in the remainder of this section.
Spatial evaluation on staggered grids.Since we now have two different grids on which to store and evolve the flow and phase fields, we need to determine how they must be positioned with respect to each other, and how to interface them so that unknown values of a field on a given grid can be quickly interpolated from the other grid.We opted for a staggered positioning to simplify interfacing: as Fig. 9 illustrates for the 2D case, every coarse node is centered on a fine grid cell.This makes interfacing the two grids quite straightforward and symmetric: Fig. 9. Staggered grids.Our LBM grids flow and phase fields are different; the phase field grid (with red nodes above) is twice as fine as the flow grid (with magenta nodes) for higher spatial accuracy of the fluid interface.Coupling between the two grids is particularly simple because of their staggered spatial locations: generally, if a phase value is needed at a flow node x  , an equal-weighted average of the corners of the corresponding fine grid is used (top), while if a flow value (like the macroscopic velocity) is needed at a phase field node x ℎ , a simple multilinear interpolation between the corners of the corresponding coarse flow cell is performed (right).
• if a coarse node x  storing the distribution g needs to evaluate any function of the phase field, it is first evaluated on the fine grid for all the corners of the fine grid cell corresponding to the coarse node x  and then averaged with equal weighting; • reversely, if a fine node x ℎ storing the distribution h needs to evaluate the velocity u locally, we use a trilinear interpolation from the corners of the coarse grid cell that x ℎ is in.
These two rules make the access of coupling values (macroscopic velocity and phase field) simple, even if the grids are not collocated.
Time integration.Having two different grid resolutions also requires a specific time integration process, as grid resolution and time-step sizes are coupled in LBM.Here again, the grid ratio being 2 makes for a simple approach: once distribution functions g  and h  are known at time   , the usual three-step (streaming, collision, boundary handling) LBM integration of the flow equation can be performed using the latest phase field evaluation derived from h  , leading to a new distribution function g +1 at time  +1 .We then proceed with two LBM updates of the phase field discretizations: we first compute h + 1 2 using the macroscopic velocity derived from g  , then compute h +1 using the averaged macroscopic velocity derived from 1 2 (g  + g +1 ).Fig. 10 summarizes this integration process.3) substeps (bottom), first using the macroscopic velocity at time   , then using the average of the velocity at   and the one at  +1 .Coupling between flow and phase field equations is indicated through hexagonal labels here.

Improving (un)normalized phase gradient
The evaluation of the normalized phase gradient ∇/|∇ | (also called the interface normal n) present in Eq. ( 2) is key to the correct time evolution of the phase field, and thus, it plays a major role in the accuracy of all the interfacial forces and coupling terms.Consequently, Li and colleagues [2021;2022] used a rotationally-symmetric evaluation ∇ [2nd]  offering an efficient second-order accurate approximation of the gradient of the phase field before normalizing it.However, even when computed on a finer grid, their approximation suffers from grid-induced oscillations: Fig. 11 (left) shows the type of ringing artifacts that even a nearly still fluid exhibits, made visually more obvious using a fast-changing colormap.While these patterned oscillations did not affect stability in their static obstacle context, their very existence is extremely likely to lead to instability when simulating dynamic solid coupling, see Fig. 12.
In [Geier et al. 2015], a quest for the evaluation of an interface normal n ≔ ∇/|∇ | which would not suffer from grid artifact led the authors to compute the normal not via differentiation, but via a direct central first-moment transform: where ∝ is used to indicate that a normalization is needed to return the unit normal n * .While this evaluation assumes a quasi-static phase field, it offers a very different and very localized way to approximate the normalized gradient.
Upon numerical testing, we found that these two different approaches to the normalized phase field gradient are in fact complementary.As Fig. 11 illustrates, n * is clearly far less accurate in the air portion of the fluid, but also far smoother (and ringing-free) in the bulk of the heavy fluid.Thus, we propose a new approximation scheme combining their strengths through: where the weighting based on the local phase field value ensures that the best estimate is used.
The evaluation at nodes of the grid storing the distribution g of the unnormalized gradient ∇ of the phase field is also crucial to prevent instabilities due to coupling.Near solid boundaries (i.e., on cut-cell nodes), we follow precisely the spatial evaluation strategy described in Sec.3.4; i.e., we evaluate ∇ [2nd]  at the corners of the fine phase grid cell before averaging their values on the coarse flow grid node: The images above visualize the magnitude of the gradient of the phase field for a nearly hydrostatic pool.While the finite-difference-based method of [Li et al. 2022] (left) gets a seemingly nice interface, the phase field gradient in fact contains ubiquitous grid-induced oscillations in the bulk of the heavy fluid as highlighted in a close-up view below.Reversely, the quasi-static approximation from [Geier et al. 2015] shows artifact-free phase gradients in the heavy fluid, but suffers from massive distortion near the light fluid (middle).By combining these two approximations (right), our weighted scheme for normal computation (Eq.( 17)) has no visible artifact anywhere.
it provides the most isotropic second-order accurate approximation in these regions of fluid-solid interactions, which we denote by ∇ [2nd]  [h] to express that the gradient is evaluated on the h grid before being transferred to the g grid.However, the same strategy in the bulk of the heavy fluid is ill-advised: as discussed earlier, small oscillations are present which can easily bring instability.We thus alter the gradient evaluation: we introduce an evaluation, denoted ∇ [2nd]  [g] which, instead, first compute the phase field  directly on nearby neighbors of the coarse flow grid, and only then use the second-order accurate evaluation of the gradient directly on the coarse flow grid.Note that this interpolation-then-differentiation approach corresponds to a larger spatial stencil than the previous differentiation-then-interpolation, and is thus less likely to suffer from the small oscillations mentioned above.Finally, we then use an evaluation which blends these two evaluations based on the value of the local phase field to determine whether we are in the bulk of the heavy fluid or not.That is, the gradient of the phase field is achieved through: Except for these new expressions, we strictly follow [Li et al. 2022] in their evaluation of interfacial forces, collisions, or wetting conditions Ð modulo the trilinear interpolations needed at times when interfacing the grids of g and h as described in Sec.3.4.

Cut-cell updates
Since we are focusing on two-way coupling in this paper, intersection between links and boundaries have to be computed on the fly as obstacles continuously moves in time.Moreover, since we adopted a bounce-back based approach to coupling as described in Sec.3.2, we have to track fresh vs. dead cells: traditionally, new distribution functions need to be assigned for fresh grid nodes (i.e., fluid nodes that were in the moving solid at the previous time step), while clean up of distribution functions need to be performed for dead cells (i.e.,  fluid nodes that become inside a moving solid); moreover, the choice of a łrefillingž interpolation scheme used to derive fresh values from nearby nodes (similar to the Gauss-Jacobi process for invalid nodes in the context of Navier-Stokes fluid solvers [Guendelman et al. 2005]) can greatly influence the simulation accuracy around solid boundaries.As we show next, our particular setup makes the update of cut-cells quite simple.
Cut-cell treatment for distribution g.Given our use of the hybrid moving bounce-back described in Sec.3.2, no specific treatment is required when updating the fluid flow: whether a grid node is covered by a solid object or not, the same exact process applies.A typical łfictitious flowž will happen inside objects, but with no negative impact on the physical behavior of the real flow.
Cut-cell treatment for distribution h.The treatment of the phase field is slightly more involved: since the phase equation and solid dynamics are only weakly coupled, even with a finer time step as we explained in Sec.3.4, leakage could happen as seen in Fig. 13.First, we need to rapidly detect grid nodes that are becoming łdeadž, i.e., nodes over which a solid object is suddenly stepping.This is achieved by keeping track, for each phase-field node, of the links that are cut by a solid object.Given our use of D3Q7, we only need 6 boolean values per node to store the presence of cut-links in an array cut-link [k].Then our detection of dead nodes becomes simple: if a current link ℓ  from a node x is intersecting an object, we check if in the previous time step, the link ℓ  at x was not cut but its opposite link ℓ  ′ was, i.e., (cut-link[k]==0)&&(cut-link[k']==1).Once this łdeadž test is done (Fig. 13), we update the cut-link array for the next time step, and we clean up the values of all dead nodes by setting all their ℎ  to zero.This corresponds to the equilibrium distribution of  = 0, meaning that we artificially set to light-air any node which becomes covered by a solid object to prevent artifacts.While other authors also suggest a specific treatment for fresh cells, our two half-steps for the integration of h (see Fig. 10) renders this step unnecessary: there is no need to specifically identify fresh cells since they will, like other regular cells, be automatically filled in via streaming from their neighbors .A final note: unlike for the distribution g, the bounce-back for h does not include a velocity, as the phase field is a scalar; cut-links are thus only performing a plain bounce-back: ℎ  (  ,  +1) =ℎ  ′ (  , ).
Solid velocity for phase advection.Finally, when a phase field node x ℎ computes the macroscopic flow velocity using trilinear interpolation from the corners of the bounding coarse flow cell, we replace the macroscopic velocity of every cut-cell corner by the solid's velocity evaluated at the link intersection point x  (see Fig. 9), to account for the presence of an obstacle during the collision step.

Implementation details
We give pseudocode of our simulation solver in Alg. 1 to summarize the order in which the various steps are performed, before going through the few remaining steps we did not discuss until now.
Memory usage.Because we use two different grid resolutions, comparing memory usage with [Li et al. 2022] is not entirely straightforward.However, for multiphase flows, one can argue that it is the resolution of the interface which matters visually.For a phase resolution of  × × , [Li et al. 2022] ends up using a total of 48 3 floats, while our method only uses 22.37 3 floats in practice: since we reduce the resolution of the D3Q27 lattice, the extra memory we save is vastly superior to the added cost of our cut-link[] array and the additional storage of the macroscopic velocity from the previous frame to provide a half-time averaged velocity for the phase field.Thus we save around 54% memory compared to [Li et al. Fig. 14.Onion Throwing.A 3D model of an onion is thrown into a water tank; as the initial splash affects its trajectory, bubbles are dragged along its path, and it comes to rest at the bottom of the tank.GPU-based implementation.Similar to [Li et al. 2022], our multiphase fluid-solid coupling solver was designed to maximize the locality of computations, thus benefitting from being implemented on GPU.The use of structure-of-arrays (SoA [Chen et al. 2021]) memory layout improves GPU memory access concurrency.For dynamic solids represented by a 3D mesh model, we need to perform link-mesh intersection at every iteration; to improve performance, we first rely on the bounding boxes of 3D models to determine cutcells, i.e., cells that straddle an object; second, we perform link-mesh intersections for all cut-cell nodes during streaming and wetting boundary handling, and update the grid cut-link array accordingly.Note finally that we also gain efficiency compared to [Li et al. 2022] by only using their second-order rotationally-symmetric ∇ evaluation, not their third-order one: the use of a finer phase grid makes the second-order approximation sufficient, and it reduces the size of the evaluation stencil significantly Ð refer to Tab. 1 for timings.Soft-start initialization.In our context of two-way coupling, we often encounter situations involving a fast-moving solid.Starting an LBM simulation in this case often creates transient compression waves which can rapidly generate numerical instabilities.We thus use the soft-start initialization strategy suggested in [Li et al. 2021], adapted to our multiphase and dynamic coupling context: for the first 20 frames, we use the full weighted rotationally-symmetric evaluation of ∇ and systematically threshold any large łimpulsež force as done in [Li et al. 2022].This soft-start initialization quickly dampens all transient compression waves.Afterwards, we return to the cheaper and more local second-order rotationally-symmetric evaluation of ∇ and never perform any form of thresholding.This is an important advantage over [Li et al. 2022]: our accuracy improvements of fluid-solid coupling remove the need for artificial thresholding as systematically performed in [Li et al. 2022], or for a more complex adaptive temporal approach [Li et al. 2020].

RESULTS AND LIMITATIONS
We now go through our tests and results (including one-way and two-way examples, with complex, thin and/or thick moving solids), before discussing the remaining limitations and possible future works.Note that our tests used a CPU-based realtime rigid body solver [Coumans and Bai 2021], while the two-phase flow simulations were run on a variety of graphics cards, including an NVIDIA A100 with 40Gb and an NVIDIA RTX A6000 with 48Gb.We also used these GPUs to render our results using a modified GPU renderer based on [Jan van Bergen 2021].

One-way coupling
One-way coupling is simpler to achieve than two-way coupling as only the effects on a fluid of the scripted motion of a solid object are simulated.Fig. 4, for instance, shows an airplane moving through a pool of water for different Reynolds numbers.The wake flow generated by the moving airplane generates bubbles and splashing commensurate with the scale implied by the Reynolds number.(In fact, in the sequence shown on top of this figure, water even splashes all the way to the ceiling of the water container, as can be seen in the video.)Similarly, Fig. 7 shows a car entering at 40 mph into a pool of water, with a clear air bubble forming in front of the car.Fig. 17 shows a spaceship emerging from underwater and generating splashes and bubbles in the process.In Fig. 20, a propeller placed on an air-water interface starts spinning and creates the typical łfrothingž action with bubbles, splashes and waves.Finally, Fig. 22 uses the same propeller, now placed higher up above the water.Its spinning creates an air flow which induces eddies and splashes.

Two-way coupling
Multiphase flow simulations involving two-way coupling are more challenging, and often lead to intricate motions.We now cover the various examples shown in this paper.
Coupling with solid objects.Fig. 14 shows an onion being thrown into a water container.When the onion breaks the water surface, a thin splash and bubbles are generated, and the motion of the onion is altered accordingly; bubbles are dragged along with the onion into the tank, before making their way back up to the surface.Fig. 21 demonstrates a similar scenario, but for a nozzle being dropped in a water tank this time.In Fig. 19, a coin-shaped stone is being skimmed on a still water surface, generating splashes and water waves.This two-way coupling example exhibits multiple bounces as expected from stone skipping.
Coupling with thin objects.Handling thin objects in the context of multiphase flows is challenging.We demonstrate a number of examples with thin shells to demonstrate how our simulator can seamlessly handle this case as well.For instance, Fig. 2 shows thin cups of various densities being dropped in a water tank: a light cup ends up floating, while heavier cups make their way down more or less fast, as the water rushes inside and pushes air away and makes the cup fall to the bottom of the tank.
Coupling with lighter objects.We also test the dropping of a very light paper boat onto water: here, surface tension keeps the paper boat afloat, either straight up, or upside down depending on how the boat was initially dropped.Surface-tension driven (i.e., capillary) water waves are also generated by the boat.Fig. 23 shows a light ACM Trans. Graph.,Vol. 42,No. 4, Article .Publication date: August 2023.
Ours [Li et al. 2022] Fig. 18.Comparison with [Li et al. 2022].We compare a one-way coupling example with two different propellers (left: off-centered axis of rotation to create more splashes; right: correct rotation axis) for the immersed boundary method of Li et al. [2022] (top) and our solver (bottom), which exhibits finer turbulent details (see close-ups) as it can handle lower air viscosity and more accurate pressure boundary conditions to create bubbles.disk (five times lighter than the stone skipping example from Fig. 19) which first bounces on a calm pool surface, before sliding for a while due to surface tension.
Multiple-object two-way coupling.While previous results only showed a single dynamic object interacting with the fluid, we can trivially handle arbitrary numbers of objects.Fig. 15 shows a water column being dropped onto the floor and blowing toy warriors away in the process; warriors flip over, and the water reacts to the warriors' movement as expected.Finally, Fig. 5 shows cows and bunnies of various densities being dropped into water.

Comparisons
Finally, we show comparisons between our solver and the previous work we built upon, as well as comparisons between real experiments and our attempts at reproducing them.
Comparison to [Li et al. 2022].Fig. 18 shows two different propellers rotating while being half-immersed in water.While one of the propellers is rotating exactly around its proper axis of rotation, the other one is made to wobble so as to generate many bubbles and splashing.To be able to compare with [Li et al. 2022] (which used the immersed boundary method on the few examples of moving objects that they provide) fairly, we match their phase field resolution.We observe that our solver (bottom) generates more bubbles and splashing, while their poor handling of the pressure boundary condition results in weaker coupling than expected (top).Note also that water splashes are very smoothed out because of the high air viscosity used in [Li et al. 2022] to ensure stability.Importantly, our solver can handle an eight-fold increase in propeller speed, while theirs fails as soon the rotation speed increases even slightly.Finally, our solver is also around 4.5 times faster on this example.Comparisons to real-life experiments.There are many exquisitelyfine phenomena induced by multiphase flows and solid coupling, and the ubiquity of high-speed cameras offers us a large playground to evaluate our numerical solver.First, we compare our solver using an example of a video sequence found online [DepositPhotos 2021] showing a key dropped into water: the key generates a veil of bubbles at an early stage; as the key drops down, the veil shrinks and turns into a complex shape, and our simulation captures the same qualitative behavior as seen in Fig 1 .Second, we also simulate a coin falling in water, a phenomenon that has well-documented types of trajectory depending on initial conditions, the dimensionless moment of inertia, and the Reynolds number [Heisinger et al. 2014].Fig. 16 shows a łflutteringž fall and a łtumblingž fall, matching the expected real behaviors as depicted in [Heisinger et al. 2014].

CONCLUSION
In this paper, we present a kinetic two-phase flow solver with twoway coupling which can simulate in a massively-parallel fashion the dynamic interaction between solid structures (be they thick or thin, light or heavy) and the two phases of a real-life fluid (e.g., with a density ratio of 800 between air and water) at low or high Reynolds numbers.We reach this milestone by altering key components of the recent work of [Li et al. 2022] in order to increase the resolution of the interface and the accuracy of crucial coupling forces based on the phase field gradient.As a result, we no longer need to use artificially high viscosity in air or force thresholding to maintain stability, even for large interaction speeds.This opens up a vast area of possible simulations of daily phenomena, of which we only scratched the surface here by showing the correct numerical capturing of stone skipping or coin falling to name a few.
Future works abound, however, if we wish to keep improving the efficiency and accuracy of fluid-structure coupling.For instance, our treatment of dead cells is particularly efficient, but does not guarantee exact mass preservation; the addition of Lagrangian particles could correct this limitation.Moreover, we have not tested our coupling with deformable solids, which may require additional efforts to prevent blowups.Finally, we found out that many realtime rigid body solvers have unacceptably poor accuracy in their collision treatment; their tendency to exhibit fast oscillations during impact for instance can be a source of instability when coupled with a two-phase flow solver at high Reynolds numbers.Further improving our coupling in order to accommodate these typical flaws of rigid body solvers without having recourse to a monolithic solver is one of the many interesting avenues of research.

Fig. 2 .
Fig. 2. Cups.Depending on its mass, a thin-shell cup dropped in a water tank will either float (top), or first float before water gushes in and sinks the cup (middle), or directly drop to the bottom (bottom).
Fig. 3. Lattice structures.Flow velocities are stored and updated in time on a D3Q27 lattice (left) for high-accuracy with lattice weights    , while a simpler (and more memory-efficient) D3Q7 lattice structure (right) with lattice weights    suffices to encode the phase field.

Fig. 4 .
Fig. 4. Airplane ditching.An airplane comes to a stop in a pool of water.Depending on the Reynolds number (top:  = 4, 000; bottom:  = 200, 000) chosen for the simulation, one can witness very different coupling behaviors, with bubbles and smooth waves (top) or much greater turbulence (bottom).

Fig. 5 .
Fig. 5. Raining Bunnies and Cows.A bunch of cows and bunnies of various densities are dropped in a water tank, exhibiting complex bubbles and swirl patterns (top).Changing the global scale of the scene increases the Reynolds number from 13,000 to 130,000 (bottom), thus showing more turbulence.

Fig. 7 .
Fig. 7. Car in a flood.We simulate a car driving around 40mph through high water before coming to a stop; the usual łair bubblež forms in front of the car, and complex water patterns appear in the wake flow region.

Fig. 8 .
Fig. 8. Floating paper boat.In this example, surface tension keeps a very light paper boat afloat.Depending on how the paper boat is initially dropped, it lands squarely (top), or shifts back into a proper floating position (middle), or simply flips around (bottom).Capillary waves are visible as well.

Fig. 10 .
Fig. 10.Time integration.A time update from   to  +1 is performed in three steps: (1) first, the flow equation is integrated using the phase field at time   (top); then the phase field is updated in two (2)+(3) substeps (bottom), first using the macroscopic velocity at time   , then using the average of the velocity at   and the one at  +1 .Coupling between flow and phase field equations is indicated through hexagonal labels here.

[Fig. 11 .
Fig. 11.Enhancing accuracy of interface normal.The images above visualize the magnitude of the gradient of the phase field for a nearly hydrostatic pool.While the finite-difference-based method of[Li et al. 2022] (left) gets a seemingly nice interface, the phase field gradient in fact contains ubiquitous grid-induced oscillations in the bulk of the heavy fluid as highlighted in a close-up view below.Reversely, the quasi-static approximation from[Geier et al. 2015] shows artifact-free phase gradients in the heavy fluid, but suffers from massive distortion near the light fluid (middle).By combining these two approximations (right), our weighted scheme for normal computation (Eq.(17)) has no visible artifact anywhere.

Fig. 12 .
Fig. 12. Ablation test for interface normal computation.On this example of a ball dropping in a water tank (left: visualizing the phase field; right: visualizing the velocity field), if the normalized phase gradient is not handled through Eq. (17), the simulation crashes when the ball hits the bottom (top); with our new interface normal evaluation (and no other change in the simulation), the dynamics is properly captured (bottom).

Fig. 13 .
Fig.13.Dead cut-cells.When phase grid node x ℎ gets covered by a solid (like the yellow disk here) between two time steps, it becomes a dead node, and its distributions ℎ  are then set to zero.

Fig. 15 .
Fig. 15.Toy Warriors.As water is dropped, toy Egyptian warriors are swept away.This two-way coupling example clearly shows how the obstacles affect the flow of water, and how the water displaces the obstacles.2022].We list memory consumption for all our examples in Tab. 1 for completeness.

ALGORITHM 1 :Fig. 16 .
Fig. 16.Coin toss in water.Based on the initial conditions described in[Heisinger et al. 2014], our solver exhibits the expected fluttering (left) or tumbling (right) behaviors of a coin falling into a water tank; we show the coin at regular time intervals to help compare with the illustrations courtesy of[Heisinger et al. 2014] given as insets.

Fig. 19 .
Fig. 19.Stone Skipping.A rotating disk labeled 'ACM' is thrown on a flat water surface with two slightly different initial velocities, exhibiting the usual stone-skipping behavior (left: multiple bounces, with the wake flow creating ripples; right: one bounce before sinking).

Fig. 20 .
Fig. 20.Water frothing.A propeller starts spinning near the water surface, quickly generating splashing, waves and bubbles.

Fig. 21 .
Fig. 21.Nozzle cleaning.A nozzle is dropped in water, generating bubbles and splashes.Note the water coming out of the outlet (middle frame).

Fig. 22 .
Fig. 22. Air-driven eddies.Same propeller as Fig. 20, but now placed above water; its spinning creates eddies and splashes on an initially still surface.

Fig. 23 .
Fig. 23.Disk Sliding.A light rotating disk is thrown over a calm pool, sliding and bouncing over the surface; compare to the stone skipping example of Fig. 19, using a heavier disk.

Table 1 .
Statistics.All examples run on a NVidia A6000 GPU, timings indicated for sequences at 60 frames-per-second.Fig.17.Spaceship surfacing from underwater.A spaceship quickly emerging from deep underwater generates complex air-water interactions, including bubbles, splashes, wetting, and backwash.