Proposal for a sub-grid lubrication correction between resolved immersed boundaries, to sit above the existing soft-sphere contact model in m_collisions.fpp. Most of the machinery this needs already exists, so the change is fairly contained. Sketching the design here before writing code, so the modelling choices can be argued with first.
Why
When two resolved particles approach, the fluid film between them thins below the grid spacing. The grid then cannot generate the lubrication resistance that should decelerate the approach, so particles come into contact too easily and too fast. The consequences are wrong collision statistics, wrong effective restitution, and wrong suspension rheology — and the error grows with volume fraction, precisely where resolved simulation is most worth doing.
For scale, at $\phi = 0.3$ the mean surface-to-surface gap between neighbours is about $0.29D$: 5.7 cells at $D = 20$ cells per diameter. An activation gap of $2\Delta x$ therefore covers roughly the bottom third of a typical approach, so this is a genuine correction rather than a wholesale replacement of the resolved physics. At $\phi = 0.5$ the mean gap falls to 1.7 cells and essentially every pair interaction would be model-driven; that is an argument for not running there and calling it resolved, not an argument for a better model.
This is also the cheapest route to a viscous validation case, which the IBM currently lacks.
Model
Three regimes, keyed on the surface gap $h$:
| gap |
active |
rationale |
| $h > h_c$ |
nothing |
grid resolves the film |
| $0 < h \le h_c$ |
lubrication correction |
grid under-resolves; supply the deficit |
| $h \le 0$ |
existing soft-sphere contact |
physical contact |
Normal component, following Nguyen & Ladd (2002) and the IBM formulation of Costa, Boersma, Westerweel & Breugem (PRE 92, 053012, 2015):
$$\mathbf{F}^{\mathrm{lub}}_{ij} = -6\pi\mu,R_{\mathrm{eff}}^{2},(\mathbf{u}_{ij}\cdot\mathbf{n})\left[\frac{1}{\max(h,,h_{\min})} - \frac{1}{h_c}\right]\mathbf{n},
\qquad R_{\mathrm{eff}} = \frac{a_i a_j}{a_i + a_j}$$
The $-1/h_c$ term is what makes this a correction: the force goes to zero continuously at the activation gap, so what the grid already produces is not double counted. $h_{\min}$ stands in for surface roughness and bounds the $1/h$ divergence; $h_{\min} \sim 10^{-3}D$ is conventional.
In MFC's non-dimensionalisation $\mu \to 1/\mathrm{Re}$, and m_ibm.fpp already reaches this as fluid_pp(fluid_idx)%Re(1) under the viscous flag.
Tangential (shearing) lubrication is a weaker, logarithmic divergence. Proposing to leave it out of the first pass and revisit if rheology work needs it; it would pair naturally with the Coulomb friction already in s_apply_ib_collision_forces_soft_sphere.
Why this is cheap here
Everything expensive already exists:
- pair detection and the
collision_lookup table
f_local_rank_owns_location(centroid_1), which already makes exactly one rank responsible for each pair- the force/torque accumulation path and its GPU atomics
patch_ib(i)%vel and %angular_vel for the relative surface velocity
Implementation sketch
1. Parameters (m_global_parameters.fpp, plus toolchain/mfc/params/definitions.py and the case schema)
integer :: lubrication_model ! 0 = off, 1 = normal-only correction
real(wp) :: lub_activation_gap ! h_c, in cells; default ~2
real(wp) :: lub_min_gap ! h_min, in particle diameters; default 1e-3
2. Widen pair detection. Both detectors currently record a pair only once the surfaces already overlap:
! m_collisions.fpp, s_detect_ib_collisions_n2
if (norm2(distance_vec) < patch_ib(pid1)%radius + patch_ib(pid2)%radius) then
This needs to become ... + h_c so that approaching-but-not-touching pairs enter the table. In the ghost-point detector s_detect_ib_collisions the equivalent change is the search stencil, currently a 3x3x3 cell cube around each ghost point, which catches roughly one cell of separation; $h_c = 2\Delta x$ wants 5x5x5. Note this only matters once that detector is re-enabled (see #1703).
3. Force routine. Either extend s_apply_ib_collision_forces_soft_sphere or add a sibling called from the same dispatcher. The branch structure falls out of the existing overlap_distance:
overlap_distance = patch_ib(pid1)%radius + patch_ib(pid2)%radius - norm2(normal_vector)
if (overlap_distance > 0._wp) then
! existing soft-sphere contact, unchanged
else if (lubrication_model > 0 .and. viscous) then
gap = -overlap_distance
if (gap < h_c) then
R_eff = a_i*a_j/(a_i + a_j)
u_n = dot_product(vel1 - vel2, normal_vector) ! includes omega x r, as the contact path already does
mu = 1._wp/fluid_pp(fluid_idx)%Re(1)
F_lub = -6._wp*pi*mu*R_eff**2*u_n &
*(1._wp/max(gap, h_min) - 1._wp/h_c)
! accumulate on pid1, equal and opposite on pid2, same atomics as contact
end if
end ifs_apply_collision_forces does not currently receive fluid_pp, so it either gains an argument or the module caches 1/Re at initialisation. The latter is probably tidier and avoids touching the call site.
4. Sign convention. Worth stating explicitly in the code, because the normal in s_apply_ib_collision_forces_soft_sphere is built as centroid_2 - centroid_1, i.e. running from $i$ to $j$. The lubrication force must resist approach, so it is repulsive when $u_n < 0$ and attractive when the pair separates.
Timestep and stiffness
$F \sim 1/h$ is stiff as $h \to h_{\min}$, and an explicit update will either need a reduced dt or a semi-implicit treatment of the normal component. Since the normal lubrication force is linear in $u_n$, it can be folded into the particle velocity update implicitly at low cost, which is what most implementations do. This interacts with the sub-stepping that the contact model already wants at high $\phi$, so the two are probably best designed together rather than separately.
Validation
- Gondret, Lance & Petit (Phys. Fluids 14, 2002) — a single sphere bouncing off a wall in a viscous fluid, wet restitution versus Stokes number. One particle, cheap, and the accepted benchmark for this correction. It also gives the IBM a viscous validation case, which it does not currently have.
- ten Cate et al. (Phys. Fluids 14, 2002) — sphere settling in a container with PIV data, for the approach dynamics.
Expected behaviour: below $\mathrm{St} \approx 10$ the film absorbs the impact and there is no rebound; above $\mathrm{St} \approx 100$ the wet restitution converges on the dry value and the correction matters only on approach. For a density ratio of 50 at $\mathrm{Re}_p \approx 44$, $\mathrm{St} \approx 244$, so a dense-suspension case sits in the regime where the contact model does not need recalibrating for film dissipation. That is convenient, but it also means Gondret should be run across a range of $\mathrm{St}$, including the low end, or the model is only exercised where it matters least.
Open questions
- Should $h_c$ be specified in cells or in particle diameters? Cells is the physically meaningful choice, since it tracks what the grid can resolve, but diameters makes cases portable across resolutions.
- Sphere-wall lubrication as well as sphere-sphere?
s_detect_wall_collisions and s_apply_wall_collision_forces_soft_sphere already exist, and Gondret is a wall case, so the wall variant is arguably needed first.
- Is a normal-only correction acceptable for the intended rheology work, or is tangential needed from the start?
- Polydisperse and non-spherical bodies: $R_{\mathrm{eff}}$ as written assumes spheres. The soft-sphere model already assumes spheres, so this is consistent, but it should be documented as a restriction rather than inherited silently.
Related: #1703 (IBM force exchange and the collision detector currently in use).
Proposal for a sub-grid lubrication correction between resolved immersed boundaries, to sit above the existing soft-sphere contact model in
m_collisions.fpp. Most of the machinery this needs already exists, so the change is fairly contained. Sketching the design here before writing code, so the modelling choices can be argued with first.Why
When two resolved particles approach, the fluid film between them thins below the grid spacing. The grid then cannot generate the lubrication resistance that should decelerate the approach, so particles come into contact too easily and too fast. The consequences are wrong collision statistics, wrong effective restitution, and wrong suspension rheology — and the error grows with volume fraction, precisely where resolved simulation is most worth doing.
For scale, at$\phi = 0.3$ the mean surface-to-surface gap between neighbours is about $0.29D$ : 5.7 cells at $D = 20$ cells per diameter. An activation gap of $2\Delta x$ therefore covers roughly the bottom third of a typical approach, so this is a genuine correction rather than a wholesale replacement of the resolved physics. At $\phi = 0.5$ the mean gap falls to 1.7 cells and essentially every pair interaction would be model-driven; that is an argument for not running there and calling it resolved, not an argument for a better model.
This is also the cheapest route to a viscous validation case, which the IBM currently lacks.
Model
Three regimes, keyed on the surface gap$h$ :
Normal component, following Nguyen & Ladd (2002) and the IBM formulation of Costa, Boersma, Westerweel & Breugem (PRE 92, 053012, 2015):
The$-1/h_c$ term is what makes this a correction: the force goes to zero continuously at the activation gap, so what the grid already produces is not double counted. $h_{\min}$ stands in for surface roughness and bounds the $1/h$ divergence; $h_{\min} \sim 10^{-3}D$ is conventional.
In MFC's non-dimensionalisation$\mu \to 1/\mathrm{Re}$ , and
m_ibm.fppalready reaches this asfluid_pp(fluid_idx)%Re(1)under theviscousflag.Tangential (shearing) lubrication is a weaker, logarithmic divergence. Proposing to leave it out of the first pass and revisit if rheology work needs it; it would pair naturally with the Coulomb friction already in
s_apply_ib_collision_forces_soft_sphere.Why this is cheap here
Everything expensive already exists:
collision_lookuptablef_local_rank_owns_location(centroid_1), which already makes exactly one rank responsible for each pairpatch_ib(i)%veland%angular_velfor the relative surface velocityImplementation sketch
1. Parameters (
m_global_parameters.fpp, plustoolchain/mfc/params/definitions.pyand the case schema)2. Widen pair detection. Both detectors currently record a pair only once the surfaces already overlap:
This needs to become$h_c = 2\Delta x$ wants 5x5x5. Note this only matters once that detector is re-enabled (see #1703).
... + h_cso that approaching-but-not-touching pairs enter the table. In the ghost-point detectors_detect_ib_collisionsthe equivalent change is the search stencil, currently a 3x3x3 cell cube around each ghost point, which catches roughly one cell of separation;3. Force routine. Either extend
s_apply_ib_collision_forces_soft_sphereor add a sibling called from the same dispatcher. The branch structure falls out of the existingoverlap_distance:s_apply_collision_forcesdoes not currently receivefluid_pp, so it either gains an argument or the module caches1/Reat initialisation. The latter is probably tidier and avoids touching the call site.4. Sign convention. Worth stating explicitly in the code, because the normal in$i$ to $j$ . The lubrication force must resist approach, so it is repulsive when $u_n < 0$ and attractive when the pair separates.
s_apply_ib_collision_forces_soft_sphereis built ascentroid_2 - centroid_1, i.e. running fromTimestep and stiffness
dtor a semi-implicit treatment of the normal component. Since the normal lubrication force is linear inValidation
Expected behaviour: below$\mathrm{St} \approx 10$ the film absorbs the impact and there is no rebound; above $\mathrm{St} \approx 100$ the wet restitution converges on the dry value and the correction matters only on approach. For a density ratio of 50 at $\mathrm{Re}_p \approx 44$ , $\mathrm{St} \approx 244$ , so a dense-suspension case sits in the regime where the contact model does not need recalibrating for film dissipation. That is convenient, but it also means Gondret should be run across a range of $\mathrm{St}$ , including the low end, or the model is only exercised where it matters least.
Open questions
s_detect_wall_collisionsands_apply_wall_collision_forces_soft_spherealready exist, and Gondret is a wall case, so the wall variant is arguably needed first.Related: #1703 (IBM force exchange and the collision detector currently in use).