A fully-integrated lattice Boltzmann method for fluid-structure interaction

TL;DR

The paper addresses the challenge of simulating fluid–structure interactions involving deformable solids in complex suspensions using a fully Eulerian approach. It introduces the lattice Boltzmann reference map technique (LBRMT), which couples a finite‑strain solid description via the reference map technique to a lattice Boltzmann fluid on a single fixed grid, augmented by a smooth flux correction at the solid–fluid interface. Key contributions include a unified FSI formulation on an Eulerian grid, a novel LB boundary condition for moving interfaces with density differences, and a detailed numerical implementation with multi-body contact handling. The method is validated against a deformable solid in a lid‑driven cavity and demonstrated on rotating, settling, and mixing scenarios, highlighting its potential for studying collective behavior in soft matter and biofluid dynamics. The results indicate that LBRMT achieves accurate fluid and solid dynamics while offering parallel scalability and the ability to simulate large ensembles of deformable agents without remeshing, albeit under small Mach numbers and near‑unity relaxation time constraints.

Abstract

We present a fully-integrated lattice Boltzmann (LB) method for fluid--structure interaction (FSI) simulations that efficiently models deformable solids in complex suspensions and active systems. Our Eulerian method (LBRMT) couples finite-strain solids to the LB fluid on the same fixed computational grid with the reference map technique (RMT). An integral part of the LBRMT is a new LB boundary condition for moving deformable interfaces across different densities. With this fully Eulerian solid--fluid coupling, the LBRMT is well-suited for parallelization and simulating multi-body contact without remeshing or extra meshes. We validate its accuracy via a benchmark of a deformable solid in a lid-driven cavity, then showcase its versatility through examples of soft solids rotating and settling. With simulations of complex suspensions mixing, we highlight potentials of the LBRMT for studying collective behavior in soft matter and biofluid dynamics.

Figures (14)

• Figure 1: Illustration of fluid--structure interaction (FSI) methods.(A) Types of FSI methods based on the solid and fluid discretization frameworks. Mesh-free methods use particles to represent both phases; Lagrangian methods use unstructured adaptive meshes for both solids and fluids; Eulerian--Lagrangian methods use a fixed Eulerian mesh for the fluid, but moving Lagrangian markers for solids; and Eulerian methods only use one fixed computational grid for both phases. (B) Overview of the hyperelastic solid deformation and the lattice Boltzmann reference map technique for FSI simulations on a fixed computational grid. A time-dependent mapping $\bm{\chi}(\bm{X},t)$ is applied to an initially undeformed solid with a reference coordinate system $\bm{X}$, resulting in a deformed coordinate system at time $t$. The inverse mapping is the reference map $\bm{\xi}(\bm{x},t)$, which maps a deformed solid back to its initial configuration on the same fixed grid. A level set function $\phi(\bm{x},t)$ is employed to define solid geometries, whose signed value determines the solid and fluid phases. To transition between the two phases, a blur zone with half-width $\varepsilon$ is defined as $\lvert\phi\rvert<\varepsilon$ to smooth out the density, velocity, and stress field.
• Figure 2: Diagram of the $D_2Q_9$ lattice model and the smooth flux correction boundary condition.(A) The blue arrows represent the nine discrete velocities $\bm{c}_i$. Each $f_i$ is a probability distribution function of a particle velocity at $(i,j)$ in the direction of the blue arrow. The empty circles represent the neighboring nodes. In the streaming step, each post-collision $\widehat{f}_i$ moves from its original position in the direction of the arrow to its neighboring nodes. Its value then becomes the new $f_i$ of these nodes in the next timestep. (B) Illustration of the smooth flux correction (SFC) for solid--fluid interface with density difference. In order to remove the outgoing flux from the higher density region to the lower density region, we add a correction flux (green arrows) to the original outgoing populations (transparent red arrows) of a node to remove additional fluxes crossing the interface. We then add these additional fluxes back to $f_0$ to enforce mass conservation. The resultant outgoing populations are labeled with red arrows with green outlines. The amount of flux removed depends on the density differences between the two regions, which can be computed from the target density based on the blur zone (Eq. \ref{['2_3_drho']}).
• Figure 3: Illustration of one-dimensional smooth flux correction.(A) Without any constraints on the solid--fluid interface, density flux goes from the higher density region ($x_s$) to the lower density region ($x_f$). This flux blurs the density difference between the two phases in the streaming and equilibrium steps, eventually averaging out the density in the domain. (B) By adding a correction to the solid node closest to the interface, i.e. removing additional outgoing flux and putting it back to the green node, we can preserve the density difference. The outgoing flux, illustrated by the purple arrow, is the difference between the density flux from $x_s$ to $x_f$ (the red arrow) and the density flux from $x_f$ to $x_s$ (the blue arrow). The correction flux (the green arrow) is equal and opposite to the outgoing flux.
• Figure 4: Stencils for the reference map advection and solid stress computation. There are three key steps to compute the divergence of the solid stress at node $(i,j)$: (A) We first calculate the gradients of the reference map field $\partial\boldsymbol{\xi}/\partial\boldsymbol{x}$ from the reference map advection; (B) then we build the half-edge deformation gradient tensor $\boldsymbol{F}$ using the computed gradients; (C) after converting the half-edge $\boldsymbol{F}$ into half-edge solid stress $\boldsymbol{\tau}_s=\textbf{f}(\bm{F})$ with a constitutive relation f, we use the four half-edge $\boldsymbol{\tau}_s$ around node $(i,j)$ to construct $\nabla\cdot\boldsymbol{\tau}_s$. Each of these steps corresponds to a panel illustrating the stencils required for discretization, with the example of the left half-edge solid stress: (A) Three nodes are used for reference map advection in the $x$ direction when $u>0$, (B) six nodes are used for constructing the left half-edge deformation gradient tensor, and (C) nine nodes are involved for all four half-edge solid stresses.
• Figure 5: Illustration of the reference map extrapolation. The extrapolation procedure starts from the first layer in the extrapolation zone and then moves outward to the next layer after all nodes have been extrapolated. (A) The first layer $l=1$ is initialized by extending one node in four directions (up, down, left, right) of the exterior solid solids (red nodes). Layer 1 nodes have only been marked with their positions, meaning they are unset (empty orange nodes) with no extrapolated reference map values. After marking all nodes in Layer 1, we proceed to compute their extrapolated reference map values $\bm{\xi}_{\text{extrap}}$. (B) For the target node $(i,j)$ (empty orange square), we initialize a scan window centered at it with a half-width $r=2$. Within the scan window, we compute its extrapolated reference map values using the valid nodes (enlarged red nodes). (C) We perform the extrapolation layer by layer. In the scenario when the linear map is ill-defined or we find fewer than three valid nodes within the scan window, we gradually increase the scan window half-width by 1. For an example case of a target node in Layer 5 (empty gray square), its scan window has been increased to a half-width $r=5$ to include more valid nodes from previous layers (enlarged red, orange, yellow, light and dark green nodes) in the extrapolation procedure.

Theorems & Definitions (2)

• proof : Mass convervation
• proof : Momentum conservation