Beyond first-order accuracy in continuous-forcing immersed boundary methods, and their well-conditioned projection-based solution

1. Problem Statement

Immersed Boundary (IB) methods are widely used to solve partial differential equations (PDEs) on complex geometries without requiring body-fitted meshes. However, traditional continuous-forcing IB methods face two significant limitations:

1. First-Order Accuracy: Despite using high-order discretizations for the underlying PDEs, standard continuous-forcing methods typically achieve only first-order global accuracy. This degradation occurs because the regularization of singular source terms (using smoothed Dirac delta functions) and the interpolation of solution values to the interface introduce errors near the boundary. Previous attempts to improve accuracy often required heuristic corrections, non-physical interface thicknesses, or suppressing independent solution branches across the interface.
2. Ill-Conditioning: When enforcing boundary conditions (like no-slip) via projection methods (treating the IB force as a Lagrange multiplier), the resulting linear system for the surface force becomes severely ill-conditioned as the ratio of Lagrangian marker spacing to Eulerian grid spacing ($\Delta s / \Delta x$) decreases. This leads to spurious, high-frequency oscillations in the computed surface stresses and makes the system difficult to solve iteratively.

2. Methodology

The authors propose a refined IB methodology that simultaneously addresses accuracy and conditioning by reformulating the problem using composite solutions and Taylor series expansions.

A. Composite Solution Formulation

Instead of treating the IB forcing as a singular source term added to a single field, the authors define the solution as a composite field $u$ (or velocity $v$ and pressure $p$) constructed from distinct interior ($u^-$) and exterior ($u^+$) solutions separated by an interface $\Gamma$:
$$u = H^+ u^+ + H^- u^-$$
where $H^\pm$ are smoothed indicator (Heaviside) functions. By applying discrete product rules to the governing equations (e.g., Poisson or Navier-Stokes), they derive a governing equation for the composite field that explicitly includes terms representing the jumps in the solution and its derivatives across the interface.

B. Taylor Series Expansion and Higher-Order Terms

The core innovation lies in analyzing the behavior of the solution within the support of the regularized Dirac delta function (DDF).

• Governing Equation: The authors expand the interior and exterior solutions using a Taylor series about the discrete surface points within the DDF support. This reveals that standard continuous-forcing methods neglect specific terms involving the jump in the normal derivative of the solution ($[u_n]_\Gamma$).
• Constraint Equation: Similarly, the interpolation constraint (enforcing the boundary condition at the interface) is re-derived using Taylor series. This yields a corrected interpolation formula that includes a term dependent on the jump in the normal derivative and the indicator function gradient.
• Result: The new formulation introduces additional terms in both the governing equation and the constraint equation. These terms account for the non-smoothness of the solution at the interface and the smoothing effects of the DDF.

C. Projection-Based Solution with Well-Conditioning

The method employs a projection-based approach (similar to the Immersed Boundary Projection Method, IBPM) to solve for the unknown interface quantities (the jumps in normal derivatives and pressure).

• Schur Complement: By performing a block-LU factorization, the system is reduced to a Schur complement system for the unknown jumps.
• Regularization: Crucially, the inclusion of the Taylor series terms introduces a nonzero diagonal term in the Schur complement matrix.
  • In traditional methods, the Schur complement approximates an ill-posed Fredholm integral equation of the first kind, leading to ill-conditioning.
  • In the proposed method, the Schur complement approximates a well-posed Fredholm integral equation of the second kind.
• Outcome: This structural change eliminates the need for ad-hoc regularization parameters or post-processing smoothing, rendering the linear system well-conditioned even for small $\Delta s / \Delta x$ ratios.

3. Key Contributions

1. Beyond First-Order Accuracy: The paper demonstrates that smoothing and interface interpolation do not inherently limit IB methods to first-order accuracy. By systematically including neglected terms derived from Taylor series, the method achieves second-order convergence in canonical Poisson problems and near-second-order (slightly sub-second-order) in incompressible Navier-Stokes simulations.
2. Natural Conditioning: The method transforms the ill-posed force computation problem into a well-posed one through formal mathematical derivation rather than heuristic fixes. This results in smooth, physically accurate surface stresses that are robust to variations in the $\Delta s / \Delta x$ ratio.
3. Unified Framework: The approach retains the computational efficiency of standard continuous-forcing methods (using regularized DDFs on Cartesian grids) while removing the need for mesh regeneration or complex stencil modifications.
4. Generalizability: While demonstrated on Poisson and Navier-Stokes equations, the framework is applicable to any IB method utilizing regularized DDFs.

4. Results

The authors validated the method through numerical experiments:

• 1D and 2D Poisson Problems:
  • The proposed method achieved second-order convergence in the $L_\infty$ and $L_2$ norms when excluding points within the DDF support, and between first and second order when including them.
  • Traditional methods remained strictly first-order.
  • The condition number of the Schur complement for the proposed method remained low and stable as $\Delta s / \Delta x$ decreased, whereas the traditional method's condition number exploded, leading to unphysical oscillations in the forcing.
• 2D Circular Couette Flow (Navier-Stokes):
  • The method successfully simulated laminar flow between rotating cylinders.
  • Velocity errors showed convergence rates approaching second order for coarser grids.
  • The computed surface forces (shear stress) were smooth and accurate across a wide range of $\Delta s / \Delta x$ ratios (from 1.3 down to 0.7), whereas the traditional method failed to produce meaningful results at lower ratios due to numerical instability.

5. Significance

This work fundamentally reframes the limitations of continuous-forcing immersed boundary methods. It proves that:

• Accuracy is not inherent to the method type: The first-order limitation is a result of neglecting specific terms in the governing and constraint equations, not an unavoidable consequence of using smoothed delta functions.
• Conditioning is solvable: The ill-conditioning plaguing projection-based IB methods is a mathematical artifact of the formulation that can be resolved by correctly modeling the solution behavior near the interface.
• Practical Impact: The proposed algorithm offers a robust, high-accuracy alternative for simulating flows with moving or deformable bodies, particularly in regimes where high resolution of the boundary is required ($\Delta s / \Delta x < 1$), without the computational penalty of solving ill-conditioned systems or the complexity of body-fitted meshes.