Accelerated MD-Forcing Immersed Boundary Method

• The paper introduces an accelerated Richardson iteration approach that refines direct forcing to enforce no-slip conditions more accurately than one-shot methods.
• AMDF-IBM leverages Eulerian–Lagrangian coupling with iterative correction to bypass the costly full inversion of the moving-boundary coupling matrix.
• Stability is maintained by tuning a lumped body-motion parameter, offering robust guidelines for both stationary and moving rigid-body simulations.

Searching arXiv for the most relevant papers on accelerated multidirect-forcing immersed boundary methods and adjacent accelerated IBM formulations. The accelerated multidirect-forcing immersed boundary method (AMDF-IBM) is a diffuse-interface immersed boundary formulation for incompressible viscous flow with rigid moving boundaries or particles, in which the fluid is solved on a fixed Cartesian grid and the boundary is represented by Lagrangian surface points moving independently of that grid. Its defining feature is the use of an accelerated multi-direct-forcing iteration to enforce the no-slip condition more accurately than one-shot direct forcing, while retaining the Eulerian–Lagrangian interpolation and spreading structure of classical immersed boundary methods. In the 2025 analysis that formalized its accuracy and stability properties, AMDF-IBM is presented as a Richardson-iteration interpretation of multi-direct forcing, with accuracy controlled by the relaxation parameter and stability controlled by a single lumped body-motion parameter (Suzuki et al., 7 Jul 2025).

1. Definition and discrete setting

In the formulation studied for AMDF-IBM, the fluid satisfies the incompressible Navier–Stokes equations with an immersed-boundary volume force,

$\nabla\cdot \mathbf{u}=0,$

$$\rho_{\rm f}\left[\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{u}\right] =-\nabla p+\mu \nabla^2\mathbf{u}+\mathbf{g},$$

where $\mathbf{g}$ is a localized forcing term used to impose the body boundary condition. The Eulerian unknown is the fluid velocity $\mathbf{u}(\mathbf{x})$ on the Cartesian grid; the Lagrangian unknowns are the boundary-point positions $\mathbf{X}_k$, prescribed rigid-body velocities $\mathbf{U}_k$, and boundary forces $\mathbf{g}(\mathbf{X}_k)$ (Suzuki et al., 7 Jul 2025).

The Eulerian–Lagrangian coupling is defined by interpolation of velocity from grid points to boundary points and spreading of boundary forces back to the grid: $$\mathbf{u}^*(\mathbf{X}_k,t+\Delta t) = \sum_{\mathbf{x}} \mathbf{u}^*(\mathbf{x},t+\Delta t)\, W(\mathbf{x}-\mathbf{X}_k)(\Delta x)^d,$$

$$\mathbf{g}(\mathbf{x},t+\Delta t) = \sum_{k=1}^N \mathbf{g}(\mathbf{X}_k,t+\Delta t)\, W(\mathbf{x}-\mathbf{X}_k)\Delta V_k.$$

Here $d$ is the spatial dimension, $$\rho_{\rm f}\left[\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{u}\right] =-\nabla p+\mu \nabla^2\mathbf{u}+\mathbf{g},$$0 is the regularized delta or weighting function, and the boundary measure satisfies

$$\rho_{\rm f}\left[\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{u}\right] =-\nabla p+\mu \nabla^2\mathbf{u}+\mathbf{g},$$1

with $$\rho_{\rm f}\left[\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{u}\right] =-\nabla p+\mu \nabla^2\mathbf{u}+\mathbf{g},$$2 the body surface area, or the perimeter in two dimensions. The direct-forcing step enforces

$$\rho_{\rm f}\left[\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{u}\right] =-\nabla p+\mu \nabla^2\mathbf{u}+\mathbf{g},$$3

A central point in the AMDF-IBM literature is that one-shot direct forcing only approximately satisfies this condition, because the interpolation–spreading pair does not exactly invert itself. In matrix form, the exact discrete force-balance problem is

$$\rho_{\rm f}\left[\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{u}\right] =-\nabla p+\mu \nabla^2\mathbf{u}+\mathbf{g},$$4

with

$$\rho_{\rm f}\left[\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{u}\right] =-\nabla p+\mu \nabla^2\mathbf{u}+\mathbf{g},$$5

and

$$\rho_{\rm f}\left[\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{u}\right] =-\nabla p+\mu \nabla^2\mathbf{u}+\mathbf{g},$$6

Direct inversion would enforce the discrete no-slip constraint very accurately, but it is expensive for moving bodies because $$\rho_{\rm f}\left[\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{u}\right] =-\nabla p+\mu \nabla^2\mathbf{u}+\mathbf{g},$$7 changes with the boundary configuration (Suzuki et al., 7 Jul 2025).

2. Accelerated multi-direct forcing as Richardson iteration

AMDF-IBM interprets multi-direct forcing as a Richardson iteration applied to the linear system $$\rho_{\rm f}\left[\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{u}\right] =-\nabla p+\mu \nabla^2\mathbf{u}+\mathbf{g},$$8: $$\rho_{\rm f}\left[\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{u}\right] =-\nabla p+\mu \nabla^2\mathbf{u}+\mathbf{g},$$9 where $\mathbf{g}$0 is the acceleration parameter and $\mathbf{g}$1 is the forcing-iteration count. Conventional multi-direct forcing is recovered at $\mathbf{g}$2; acceleration enters by multiplying the direct-forcing correction by $\mathbf{g}$3 (Suzuki et al., 7 Jul 2025).

The study considers three specific choices: $\mathbf{g}$4 For the Peskin kernels $\mathbf{g}$5 and $\mathbf{g}$6,

$\mathbf{g}$7

with

$\mathbf{g}$8

independent of $\mathbf{g}$9. Hence,

$\mathbf{u}(\mathbf{x})$0

This characterization makes the “accelerated” qualifier precise. AMDF-IBM is not accelerated merely because more forcing is applied; it is accelerated because the fixed-point iteration is relaxed with a parameter chosen to match the spectral scale of the discrete coupling operator. The same paper shows that for practical discretizations of Peskin-kernel IBM,

$\mathbf{u}(\mathbf{x})$1

so the Gsell correction and the relaxation based on $\mathbf{u}(\mathbf{x})$2 become effectively the same choice in practice (Suzuki et al., 7 Jul 2025).

This places AMDF-IBM within a broader family of direct-forcing enhancements. Unlike penalty formulations, it does not rely on stiffness to enforce boundary data. Unlike monolithic constraint formulations, it does not solve the coupled Eulerian–Lagrangian system exactly at each step. Its distinctive position is that it preserves the explicit direct-forcing workflow while accelerating convergence toward the exact discrete force solution.

3. Accuracy of no-slip enforcement

The accuracy analysis in AMDF-IBM is framed around the maximum no-slip mismatch at boundary points,

$\mathbf{u}(\mathbf{x})$3

Across all stationary and moving test cases examined, the no-slip error decreases monotonically with the number of iterations $\mathbf{u}(\mathbf{x})$4, is much smaller for $\mathbf{u}(\mathbf{x})$5 than for $\mathbf{u}(\mathbf{x})$6, and already improves by about one order of magnitude without iteration, $\mathbf{u}(\mathbf{x})$7, if $\mathbf{u}(\mathbf{x})$8 is optimally chosen (Suzuki et al., 7 Jul 2025).

For the two kernels explicitly examined, the near-optimal values are

$\mathbf{u}(\mathbf{x})$9

For a circular cylinder, the reported matrix characteristics are

$\mathbf{X}_k$0

$\mathbf{X}_k$1

which explains why $\mathbf{X}_k$2, $\mathbf{X}_k$3, and $\mathbf{X}_k$4 are numerically close.

The empirical scope of this result is unusually broad within the tested class. The study reports that the optimal acceleration parameter is essentially independent of boundary point spacing over

$\mathbf{X}_k$5

independent of whether points are evenly or unevenly distributed, and independent of body shape and spatial dimension for the examined circle, ellipse, and sphere problems. The stationary-boundary accuracy tests comprise a circular cylinder fixed in Poiseuille flow, an elliptical cylinder fixed in Poiseuille flow, and a sphere fixed in planar Poiseuille flow (Suzuki et al., 7 Jul 2025).

The same analysis also clarifies an often-misunderstood point about conditioning. For a circular cylinder with evenly spaced boundary points, $\mathbf{X}_k$6 decreases sharply when $\mathbf{X}_k$7, making $\mathbf{X}_k$8 ill-conditioned if one attempts direct inversion. The iterative method, however, remains usable. This distinction is important: AMDF-IBM improves practical no-slip enforcement without requiring explicit inversion of the full moving-boundary coupling matrix.

4. Stability of moving-body simulations

The stability contribution of AMDF-IBM concerns freely moving rigid bodies. Using the discrete translational body-motion equation with internal-fluid mass effect, the analysis identifies a single lumped parameter controlling body-motion stability: $\mathbf{X}_k$9 and, with forcing iteration, the controlling parameter becomes

$\mathbf{U}_k$0

Here $\mathbf{U}_k$1 is the IBM acceleration parameter, $\mathbf{U}_k$2 is the solid-to-fluid density ratio, $\mathbf{U}_k$3 is the surface-to-volume ratio, and $\mathbf{U}_k$4 is the Eulerian grid spacing (Suzuki et al., 7 Jul 2025).

For canonical geometries this reduces to

$\mathbf{U}_k$5

$\mathbf{U}_k$6

The appendix derives an amplification factor $\mathbf{U}_k$7 for the total immersed-boundary force under $\mathbf{U}_k$8 forcing iterations, so that iteration effectively strengthens the body-force coupling. The practical criterion reported by the paper is

$\mathbf{U}_k$9

This result changes the interpretation of AMDF-IBM from a purely accuracy-oriented forcing method to a coupled fluid–body stability problem. Instability is favored by larger $\mathbf{g}(\mathbf{X}_k)$0, smaller $\mathbf{g}(\mathbf{X}_k)$1, coarser relative spatial resolution, and more iterations when these increase $\mathbf{g}(\mathbf{X}_k)$2. The moving-boundary validation set comprises a circular cylinder moving in Poiseuille flow, a sedimenting circular cylinder, and a sphere moving in planar Poiseuille flow. In these tests, stable and unstable cases collapse when plotted against $\mathbf{g}(\mathbf{X}_k)$3 or $\mathbf{g}(\mathbf{X}_k)$4, largely independently of $\mathbf{g}(\mathbf{X}_k)$5, $\mathbf{g}(\mathbf{X}_k)$6, $\mathbf{g}(\mathbf{X}_k)$7, $\mathbf{g}(\mathbf{X}_k)$8, $\mathbf{g}(\mathbf{X}_k)$9, $$\mathbf{u}^*(\mathbf{X}_k,t+\Delta t) = \sum_{\mathbf{x}} \mathbf{u}^*(\mathbf{x},t+\Delta t)\, W(\mathbf{x}-\mathbf{X}_k)(\Delta x)^d,$$0, boundary-point spacing $$\mathbf{u}^*(\mathbf{X}_k,t+\Delta t) = \sum_{\mathbf{x}} \mathbf{u}^*(\mathbf{x},t+\Delta t)\, W(\mathbf{x}-\mathbf{X}_k)(\Delta x)^d,$$1, gravitational acceleration $$\mathbf{u}^*(\mathbf{X}_k,t+\Delta t) = \sum_{\mathbf{x}} \mathbf{u}^*(\mathbf{x},t+\Delta t)\, W(\mathbf{x}-\mathbf{X}_k)(\Delta x)^d,$$2, and the kernel choice $$\mathbf{u}^*(\mathbf{X}_k,t+\Delta t) = \sum_{\mathbf{x}} \mathbf{u}^*(\mathbf{x},t+\Delta t)\, W(\mathbf{x}-\mathbf{X}_k)(\Delta x)^d,$$3 versus $$\mathbf{u}^*(\mathbf{X}_k,t+\Delta t) = \sum_{\mathbf{x}} \mathbf{u}^*(\mathbf{x},t+\Delta t)\, W(\mathbf{x}-\mathbf{X}_k)(\Delta x)^d,$$4 in the tested cases (Suzuki et al., 7 Jul 2025).

The resulting guideline is explicit. For best no-slip accuracy, choose

$$\mathbf{u}^*(\mathbf{X}_k,t+\Delta t) = \sum_{\mathbf{x}} \mathbf{u}^*(\mathbf{x},t+\Delta t)\, W(\mathbf{x}-\mathbf{X}_k)(\Delta x)^d,$$5

that is,

$$\mathbf{u}^*(\mathbf{X}_k,t+\Delta t) = \sum_{\mathbf{x}} \mathbf{u}^*(\mathbf{x},t+\Delta t)\, W(\mathbf{x}-\mathbf{X}_k)(\Delta x)^d,$$6

then compute $$\mathbf{u}^*(\mathbf{X}_k,t+\Delta t) = \sum_{\mathbf{x}} \mathbf{u}^*(\mathbf{x},t+\Delta t)\, W(\mathbf{x}-\mathbf{X}_k)(\Delta x)^d,$$7 or $$\mathbf{u}^*(\mathbf{X}_k,t+\Delta t) = \sum_{\mathbf{x}} \mathbf{u}^*(\mathbf{x},t+\Delta t)\, W(\mathbf{x}-\mathbf{X}_k)(\Delta x)^d,$$8 and maintain it below about $$\mathbf{u}^*(\mathbf{X}_k,t+\Delta t) = \sum_{\mathbf{x}} \mathbf{u}^*(\mathbf{x},t+\Delta t)\, W(\mathbf{x}-\mathbf{X}_k)(\Delta x)^d,$$9. If it is too large, the paper recommends reducing $$\mathbf{g}(\mathbf{x},t+\Delta t) = \sum_{k=1}^N \mathbf{g}(\mathbf{X}_k,t+\Delta t)\, W(\mathbf{x}-\mathbf{X}_k)\Delta V_k.$$0, increasing $$\mathbf{g}(\mathbf{x},t+\Delta t) = \sum_{k=1}^N \mathbf{g}(\mathbf{X}_k,t+\Delta t)\, W(\mathbf{x}-\mathbf{X}_k)\Delta V_k.$$1, refining spatial resolution, or reducing effective iteration amplification. The analysis is, however, a reduced discrete-force scaling argument rather than a closed-form eigenanalysis of the full coupled fluid–body system, so its conclusions are presented as robust practical criteria rather than a complete modal stability theory (Suzuki et al., 7 Jul 2025).

5. Computational acceleration beyond the forcing iteration

AMDF-IBM addresses the forcing-side accuracy–stability tradeoff, but large-scale acceleration in immersed-boundary solvers has also been pursued through pressure-correction preconditioning, strong-coupling architecture, and GPU implementation. These developments are adjacent rather than identical to AMDF-IBM.

A direct-forcing moving-body solver based on a SIMPLE reformulation keeps pressure and immersed-boundary forces coupled, forms a pressure–forces correction system, and applies a block reduction to a primal Schur complement

$$\mathbf{g}(\mathbf{x},t+\Delta t) = \sum_{k=1}^N \mathbf{g}(\mathbf{X}_k,t+\Delta t)\, W(\mathbf{x}-\mathbf{X}_k)\Delta V_k.$$2

with a Laplacian-based preconditioner. Its central theoretical result is the spectral bound

$$\mathbf{g}(\mathbf{x},t+\Delta t) = \sum_{k=1}^N \mathbf{g}(\mathbf{X}_k,t+\Delta t)\, W(\mathbf{x}-\mathbf{X}_k)\Delta V_k.$$3

which explains why the iteration count for the preconditioned pressure–forces system is reported as nearly constant regardless of grid resolution, Reynolds number, and the number of spheres (Yovel et al., 25 Jan 2025). This is not classical multidirect forcing, but it targets the same no-slip enforcement bottleneck through a pressure-first coupled correction solve.

A different acceleration route is represented by GPU-native immersed-boundary solvers. The sharp-interface solver ViCar3D uses a ghost-cell method rather than direct forcing, but it demonstrates implementation strategies directly transferable to accelerated Eulerian–Lagrangian solvers: OpenACC, CUDA, CUDA-aware MPI, persistent device-resident data regions, batched cuBLAS inversion of many $$\mathbf{g}(\mathbf{x},t+\Delta t) = \sum_{k=1}^N \mathbf{g}(\mathbf{X}_k,t+\Delta t)\, W(\mathbf{x}-\mathbf{X}_k)\Delta V_k.$$4 matrices, and communication-aware halo exchange. The paper reports approximately $$\mathbf{g}(\mathbf{x},t+\Delta t) = \sum_{k=1}^N \mathbf{g}(\mathbf{X}_k,t+\Delta t)\, W(\mathbf{x}-\mathbf{X}_k)\Delta V_k.$$5 speedup in a node-to-node comparison, capability up to 200 million mesh points on one node with four GPUs, and maximum strong- and weak-scaling efficiencies of $$\mathbf{g}(\mathbf{x},t+\Delta t) = \sum_{k=1}^N \mathbf{g}(\mathbf{X}_k,t+\Delta t)\, W(\mathbf{x}-\mathbf{X}_k)\Delta V_k.$$6 and $$\mathbf{g}(\mathbf{x},t+\Delta t) = \sum_{k=1}^N \mathbf{g}(\mathbf{X}_k,t+\Delta t)\, W(\mathbf{x}-\mathbf{X}_k)\Delta V_k.$$7 (Kumar et al., 22 May 2025). Its relevance to AMDF-IBM is infrastructural rather than methodological.

For rigid-body fluid–structure interaction, a strongly coupled DF-IBM algorithm has been developed in which an inner iterative direct-forcing loop is embedded inside an outer implicit Newton–Euler coupling, while expensive pressure-correction steps are excluded from the inner FSI iterations. The method is stabilized by fixed relaxation of rigid-body kinematics and is explicitly motivated by low-density-ratio robustness (Farah et al., 27 Apr 2026). This suggests that, in moving-body settings, acceleration of the forcing stage alone may be insufficient unless the larger coupling architecture is also stabilized.

A further adjacent development is a fast sharp-interface IBM for moving and deformable bodies that preserves the discrete Laplacian matrix and changes only the right-hand side of the projection Poisson problem through a consistent mass correction. It does not use multidirect forcing, but it shows that projection-induced boundary error can be reduced without sacrificing direct Poisson solvers (Vagnoli et al., 8 Jun 2026). A plausible implication is that accelerated multidirect forcing and projection-consistency corrections address complementary sources of no-slip degradation.

6. Related formulations, alternatives, and scope

AMDF-IBM sits within a broader spectrum of immersed-boundary boundary-enforcement strategies. The following comparison is methodological rather than hierarchical.

Method family	Boundary enforcement	Relation to AMDF
AMDF-IBM	Accelerated Richardson iteration for direct forcing	Canonical accelerated multidirect forcing (Suzuki et al., 7 Jul 2025)
BTDF	One-shot direct forcing with thickened boundary shell	MDF-like accuracy at nearly DF cost (Jiang et al., 2018)
Accelerated SIMPLE direct forcing	Coupled pressure–force correction with Laplacian preconditioning	Pressure-side acceleration, not MDF iteration (Yovel et al., 25 Jan 2025)
Sharp-interface ghost-cell IBM	Reconstruction in ghost and fresh cells	Not direct forcing; relevant for GPU/HPC implementation (Kumar et al., 22 May 2025)
IBDL / constrained-force IB	Well-conditioned or monolithic constraint solve	Operator-level alternatives to repeated forcing (Leathers, 2022)

One important misconception is that “accelerated immersed boundary method” is synonymous with AMDF-IBM. It is not. The boundary thickening-based direct-forcing method (BTDF) was proposed as a low-cost surrogate for MDF, IVC, and RKPM. In the reported cylinder benchmark, BTDF with no iteration outperforms MDF with four iterations, while MDF converges to the same solution when enough iterations are taken (Jiang et al., 2018). This is an acceleration of boundary enforcement by operator modification rather than by faster convergence of a Richardson iteration.

A second misconception is that AMDF-IBM is interchangeable with sharp-interface ghost-cell methods. It is not. Ghost-cell solvers reconstruct state variables in ghost and fresh cells and embed those reconstructions into the discretization, whereas AMDF-IBM uses diffuse-interface interpolation and spreading of direct-forcing corrections. The methodological distinction is substantial even when both are accelerated on modern hardware (Kumar et al., 22 May 2025).

A third misconception is that repeated explicit forcing is the only route to stronger no-slip enforcement. The Immersed Boundary Double Layer method replaces the regularized single-layer, first-kind operator of classical constraint formulations with a regularized second-kind operator, yielding small and essentially mesh-independent GMRES iteration counts in the reported Helmholtz, Brinkman, Poisson, and Stokes tests (Leathers, 2022). Likewise, unsplit constrained rigid-body immersed-boundary formulations enforce no-slip exactly at Lagrangian markers through a monolithic saddle-point solve with Schur-complement preconditioning rather than through repeated forcing corrections (Kallemov et al., 2015, Liska et al., 2016). These are not AMDF methods, but they solve the same boundary-enforcement problem at a different algebraic level.

The scope of the current AMDF-IBM stability theory is also delimited in the source study. Its results are best interpreted as applying to rigid-body motion, weakly coupled IBM–LBM schemes, the tested Reynolds-number ranges, the tested geometries of cylinder, ellipse, and sphere, and the Peskin kernels $$\mathbf{g}(\mathbf{x},t+\Delta t) = \sum_{k=1}^N \mathbf{g}(\mathbf{X}_k,t+\Delta t)\, W(\mathbf{x}-\mathbf{X}_k)\Delta V_k.$$8 and $$\mathbf{g}(\mathbf{x},t+\Delta t) = \sum_{k=1}^N \mathbf{g}(\mathbf{X}_k,t+\Delta t)\, W(\mathbf{x}-\mathbf{X}_k)\Delta V_k.$$9. The paper does not establish universality for arbitrary kernels, highly deforming bodies, strongly nonlinear transient collision-rich particulate flows, or strongly coupled implicit FSI schemes (Suzuki et al., 7 Jul 2025). Within that domain, however, AMDF-IBM provides a compact synthesis: the optimal acceleration parameter is set by the kernel-dependent spectral scale of the discrete coupling operator, while stability is governed by the lumped parameter $d$0, or $d$1 when multiple forcing iterations are used.