PaperlayerView exact revision on arXiv

Source-linked AI summary

An Immersed Boundary Method for Rigid Bodies

B. Kallemov, A. Pal Singh Bhalla, B. E. Griffith, A. Donev

arXiv:1505.07865v2math.NAcs.NA

TL;DR

Exact enforcement of rigid-body no-slip constraints requires solving a coupled fluid-and-force system, while fractional-step approaches impose them only approximately. This paper preconditions that constrained system using an analytical Schur-complement approximation; tests report robust solver behavior and good agreement for steady Stokes flow around cylinders.

  • Problem

    Fractional-step approaches impose the no-slip constraint only approximately, while exact enforcement requires solving a coupled system for fluid velocity, pressure, and rigidity forces.

    Source: §I · p. 4, §VIII · p. 50

  • Method

    The method solves the fully constrained fluid problem with a preconditioned Krylov solver and an analytical Schur-complement approximation based on near translational and rotational invariance.

    Source: §I · p. 4

  • Results

    The preconditioner leads to robust convergence when markers are sufficiently far apart for a well-conditioned mobility matrix, and cylinder-drag results show very good agreement with a boundary-integral method.

    Source: E. Steady Stokes flow around cylinders · p. 43, §VIII · p. 50

  • Takeaways & Limitations

    The method offers an alternative for steady or nearly steady Stokes suspensions and can handle nontrivial boundary conditions using standard fluid solvers.

    Source: §I · pp. 4–5

    Robust preconditioning is conditioned on keeping immersed-boundary markers sufficiently far apart to ensure a well-conditioned mobility matrix.

    Source: §VIII · p. 50

Abstract

from arXiv · show

We develop an immersed boundary (IB) method for modeling flows around fixed or moving rigid bodies that is suitable for a broad range of Reynolds numbers, including steady Stokes flow. The spatio-temporal discretization of the fluid equations is based on a standard staggered-grid approach. Fluid-body interaction is handled using Peskin's IB method; however, unlike existing IB approaches to such problems, we do not rely on penalty or fractional-step formulations. Instead, we use an unsplit scheme that ensures the no-slip constraint is enforced exactly in terms of the Lagrangian velocity field evaluated at the IB markers. Fractional-step approaches, by contrast, can impose such constraints only approximately. Imposing these constraints exactly requires the solution of a large linear system that includes the fluid velocity and pressure as well as Lagrange multiplier forces that impose the motion of the body. To solve this system efficiently, we develop a preconditioner for the constrained IB formulation that is based on an analytical approximation to the Schur complement. This approach is enabled by the near translational and rotational invariance of Peskin's IB method. We demonstrate that only a few cycles of a geometric multigrid method for the fluid equations are required in each application of the preconditioner, and we demonstrate robust convergence of the overall Krylov solver despite the approximations made in the preconditioner. We apply the method to a number of test problems at zero and finite Reynolds numbers, and we demonstrate first-order convergence of the method to several analytical solutions and benchmark computations.

I. INTRODUCTION

The paper replaces split or penalty-based rigid-body IB methods with a constrained formulation that enforces no-slip exactly and can include steady Stokes flow. It solves for fluid variables and body forces together, using a Schur-complement approximation to precondition the system; the semi-continuum formulation represents the body with discrete markers and the fluid as a continuum.

I. INTRODUCTION

  • Splitting methods impose no-slip only approximately and fail in the steady Stokes limit, motivating an unsplit solve with controlled accuracy.The paper also notes that splitting can permit flow to penetrate rigid obstacles.Source: §I · p. 3
  • The method directly solves a coupled velocity-pressure-force system and uses an analytical Schur-complement approximation as a scalable preconditioner.The preconditioner exploits Peskin IB's near translational and rotational invariance; the authors report computational complexity only a few times larger than the corresponding problem without rigid bodies.Source: §I · p. 4
  • The method is designed for broad Reynolds numbers and boundary conditions, while avoiding dense Schur-complement precomputation.The authors contrast their on-the-fly approximation with approaches that form the Schur complement densely and expensively.Source: §I · pp. 4–5

II. SEMI-CONTINUUM FORMULATION

  • Unknown interaction forces enforce the prescribed marker motion by making the interpolated fluid velocity match each marker's known velocity.For the regularized marker constraint, flow may still leak between markers; a staggered-grid formulation can greatly reduce this leak.Source: §II · p. 8
  • The semi-continuum formulation represents the rigid body by discrete markers while retaining a continuum description of the fluid.The authors note that constraining a finite number of markers in an infinite-dimensional velocity field avoids the marker-count restriction of a fully discretized strong formulation.Source: §II · p. 9
  • The force interpretation is delicate for interior markers, and the numerical studies are restricted to stationary bodies.With finite inertia and a density mismatch, the inertial terms inside the body require modification; sufficiently many interior markers are also needed to prevent spurious angular momentum.Source: §II · pp. 9–10

III. DISCRETE FORMULATION

The formulation couples staggered-grid fluid equations to rigid-body markers through adjoint interpolation and spreading operators, then solves a constrained saddle-point system using its mobility Schur complement.

III. DISCRETE FORMULATION

  • The fluid uses a second-order staggered-grid discretization, with explicit advection, implicit viscosity, and viscous CFL number β.Source: §III · p. 10
  • Pressure and constraint force are solved as Lagrange multipliers enforcing incompressibility and rigidity within the coupled formulation.Source: §III · p. 10

A. Force spreading and velocity interpolation

B. Rigidly-constrained Stokes problem

C. Mobility matrix

  • The marker mobility matrix is the Schur complement, with each block mapping a marker force to the velocity induced at another marker.Source: C. Mobility matrix · p. 13
  • The method approximates the mobility matrix to build robust preconditioners for the mobility subproblem.Source: C. Mobility matrix · p. 14
  • The saddle-point system’s conditioning depends on the mobility spectrum; near-zero eigenvalues can represent physical null modes, such as sphere compression.Source: C. Mobility matrix · p. 14

D. Periodic steady Stokes flow

Periodic steady Stokes flow requires a solvability constraint and a force-balanced mean-velocity choice. The mobility preconditioner approximates pairwise interactions under simplifying boundary and grid-invariance assumptions, without guaranteeing positive definiteness for every marker configuration.

D. Periodic steady Stokes flow

IV. APPROXIMATING THE MOBILITY MATRIX

  • The exact mobility Schur complement is impractical to form, so the method uses a dense, low-rank approximation built from analytical pairwise mobility functions.Source: §IV · pp. 15–16
  • The approximation ignores boundary-condition specifics and assumes near translational and rotational invariance, leaving the Krylov solver to handle boundary and grid effects.Source: §IV · p. 15
  • The empirical mobility fit is not constrained to be positive definite for every marker configuration.Source: §IV · pp. 16–17

A. Physical Constraints

The paper relates pairwise marker mobility to fluid Green’s functions across steady, inviscid, and finite-viscosity regimes, then uses known asymptotic behavior to guide approximations.

A. Physical Constraints

  • The fluid Green’s function has the pairwise mobility form; steady Stokes flow gives the Oseen tensor, while finite β gives the Brinkmanlet.Source: A. Physical Constraints · p. 17
  • For nonzero Reynolds numbers, fβ(r) and gβ(r) have leading-order asymptotic r^-3 decay, corresponding to an electric dipole and originating in the incompressibility constraint.Source: A. Physical Constraints · p. 18
  • For steady Stokes flow, the RPY tensor approximates pairwise mobility and guarantees an SPD mobility matrix; for the six-point kernel its fitted hydrodynamic radius is 1.47 h.Source: A. Physical Constraints · p. 18

B. Empirical Fits

The empirical fits combine numerical mobility data with short- and long-distance constraints, accounting for periodic-domain corrections and discretization effects.

B. Empirical Fits

  • Periodic corrections can dominate the inviscid dipole contribution, while discrete solver and kernel effects matter at short distances, motivating empirical fits constrained by known limiting behavior.Source: B. Empirical Fits · pp. 19–20
  • The fitting procedure samples pairwise parallel and perpendicular mobilities by applying unit forces to randomly directed markers in a periodic domain.Source: B. Empirical Fits · p. 20
  • The normalized mobility data are well fit by smooth functions across distances, with separate fits for steady Stokes and finite β.Source: B. Empirical Fits · p. 20
  • For steady Stokes flow, the numerical mobility data are in reasonable agreement with the RPY tensor.Source: B. Empirical Fits · p. 21

V. LINEAR SOLVER

The constrained Stokes system is solved with an outer flexible Krylov method and a Schur-complement preconditioner that uses approximate fluid and mobility solves. The unconstrained fluid solver relies on multigrid-based Stokes preconditioning, which keeps its cost to a few times that of a scalar Poisson solve.

V. LINEAR SOLVER

  • Flexible GMRES solves the constrained Stokes problem because Krylov solves inside the preconditioner make it non-linear as an operator.The implementation uses IBAMR and PETSc, with uniform grids in this work.Source: §V · p. 21

A. Preconditioner for the constrained Stokes system

  • The constrained-system preconditioner approximately solves the fluid subproblem, computes marker slip, solves the mobility-based Schur complement, and optionally corrects the fluid solution.The mobility approximation is constructed as described in Section IV.Source: A. Preconditioner for the constrained… · p. 22
  • The full Schur preconditioner uses two inexact fluid solves per application, while the lower triangular version omits the correction solve and uses one.With Ns inner iterations, the two versions require 2Ns and Ns applications of the unconstrained Stokes preconditioner, respectively.Source: A. Preconditioner for the constrained… · pp. 22–23

B. Unconstrained Fluid Solver

  • The unconstrained Stokes solver with multigrid preconditioning costs only a few times as much as solving a scalar Poisson problem.A block lower triangular Stokes preconditioner can omit the pressure multigrid sweep and may be preferred at zero Reynolds number; including the Poisson velocity projection improves outer-solver performance.Source: B. Unconstrained Fluid Solver · pp. 23–24

C. Mobility Solver

Mobility solves use direct or factored approximations, while conditioning depends strongly on marker spacing and body representation. In convergence tests, the lower triangular preconditioner is more efficient overall, and the solver remains relatively robust under refinement, with a gradual increase in outer iterations.

C. Mobility Solver

  • For well-spaced markers with s/h ≳ 2, the mobility approximation is typically SPD and can be solved by Cholesky; an SVD pseudoinverse filters small or negative eigenvalues when needed.A factorization can be reused while the body remains fixed; effective mobility solves for multiple moving rigid bodies remain future work.Source: C. Mobility Solver · p. 24

VI. CONDITIONING OF THE MOBILITY MATRIX

  • For s/h ≳ 1, the nontrivial mobility eigenvalues remain bounded away from zero under refinement, while the condition number is O(N) for the spherical-shell example.The sphere also has a trivial zero eigenvalue corresponding to uniform compression.Source: §VI · p. 26
  • At s/h ≈ 1, the condition number reaches O(10^6–10^7), and at s ≈ h/2 the constrained Stokes problem cannot be solved in double precision; s/h ≈ 2 is a heuristic balance with small leak.The authors caution that this spacing recommendation is not a precise estimate.Source: §VI · p. 26
  • The discrete mobility and empirical approximation have similar spectra near s ≈ 2h, but both worsen in conditioning near s ≈ 1.5h; the RPY approximation is better conditioned at smaller spacings.These observations implicate both kernel regularization and discretization artifacts in ill-conditioning.Source: §VI · p. 28

A. Empirical convergence of GMRES

  • In steady Stokes flow, the optimal inner iteration count depends on conditioning: Ns = 2 is empirically optimal for the ill-conditioned filled sphere, while Ns has less effect for well-conditioned shells.Larger Ns can waste computational effort through oversolving; higher Reynolds number cases are easier to solve, especially for filled spheres.Source: A. Empirical convergence of GMRES · p. 30
  • The lower triangular preconditioner is more efficient overall under refinement; tighter solver tolerances require some fluid correction inside the preconditioner.With Ns = 4, FGMRES iterations rise slowly as grid and shell resolution increase.Source: A. Empirical convergence of GMRES · p. 31

B. Flow through a nozzle

In nozzle flow, the rigid IB method limits leakage across the body boundary much better than the splitting-based comparison, especially at low Reynolds number. The paper also checks steady Stokes flow between concentric shells, where hydrodynamic radii approach geometric radii as the grid is refined.

B. Flow through a nozzle

  • The nozzle test uses a 256 × 128 grid, markers about two grid spacings apart, and a pressure jump of 2 across the channel.The flow is evolved until essentially steady; Reynolds number is estimated from the maximum opening velocity and opening width.Source: B. Flow through a nozzle · p. 32
  • The flow ratio exceeds 0.99 with the rigid IB method at all Reynolds numbers, while the splitting method gives 0.935 at Re ≈19 and needs a 1250-fold smaller time step to reach 0.94 at Re ≈0.2.Figure 5 compares the rigid IB results in the left panels with the splitting method in the right panels; the latter shows stronger penetration at low Reynolds number.Source: B. Flow through a nozzle · pp. 33–34, Figure 5 · p. 33

C. Stokes flow between two concentric shells

  • The concentric-shell test prescribes a moving outer shell and stationary inner shell, with markers spaced about two grid cells apart.The inner and outer shells have a radius ratio of four and use the same recursive triangulation scheme.Source: C. Stokes flow between two concentric… · pp. 34–35
  • A marker-covered shell acts hydrodynamically like a sphere of radius Rh ≈ Rg + a, so comparisons with theory use hydrodynamic rather than geometric radii.For the six-point kernel, the marker hydrodynamic radius is approximately 1.47 h.Source: C. Stokes flow between two concentric… · p. 35
  • As the grid is refined, the effective hydrodynamic radii of both concentric shells converge toward their geometric radii at approximately a linear rate.The inner-shell radius is estimated from drag and the outer-shell radius from drag on the inner sphere using the previously estimated inner radius.Source: C. Stokes flow between two concentric… · p. 36, Table I · p. 36

1. Convergence of fluid flow (pressure and velocity)

For flow around a fixed sphere inside a moving spherical cavity, velocity and pressure errors converge at first order in L1 and L2 norms. Surface traction estimates converge more slowly, while slit-channel tests show consistent effective-radius behavior across resolutions.

1. Convergence of fluid flow (pressure and velocity)

  • Most flow-field error is localized near the inner shell and inside it, where a spurious pressure gradient produces a counter-rotating vortex.The computed error is compared with the theoretical solution; using geometric rather than hydrodynamic radii would make the errors much larger.Source: §1 · pp. 36–37, Figure 7 · p. 37
  • Velocity and pressure errors show asymptotically first-order convergence in the L1 and L2 norms.The pressure is not expected to converge in the L∞ norm because it jumps across the interface; the reported data are consistent with that expectation.Source: §1 · p. 38, Table II · p. 38, Table III · p. 38

2. Convergence of Lagrangian forces (surface stresses)

  • Pointwise tractions appear to approach theoretical values with grid refinement, but convergence is slow and the reported data do not establish its order.Some lower-resolution results show qualitatively incorrect behavior, and marker-grid nonuniformity produces scatter among pointwise values.Source: §2 · pp. 39–40, Figure 8 · p. 39
  • Low-order traction moments are more accurate than pointwise tractions, but wider marker spacing increases flow-field error by allowing penetration through the shell.The observed global traction-error reduction with wider spacing may be specific to steady Stokes flow.Source: §2 · p. 40

D. Steady Stokes flow around sphere in a slit channel

  • In the slit-channel tests, effective-radius data from different resolutions are consistent and approach the same limit as the periodic domain width increases.Even a model with 12 markers is reported to represent a sphere of effective radius Rh reasonably accurately across the tested boundary locations.Source: D. Steady Stokes flow around sphere in… · p. 41

E. Steady Stokes flow around cylinders

For steady Stokes flow around cylinders, the rigid IB method agrees well with boundary-integral results when effective hydrodynamic radii are used. The study also examines unsteady cylinder flow but restricts attention to steady states up to Re ≤100.

E. Steady Stokes flow around cylinders

  • More resolved cylinders have effective hydrodynamic radii closer to their geometric radii; filling the cylinder with markers enlarges the effective radius and worsens mobility-matrix conditioning.The authors therefore advise against interior markers for cylinders at zero Reynolds number.Source: E. Steady Stokes flow around cylinders · pp. 42–43, Table IV · p. 43
  • Reliable cylinder-array calculations reach relative gaps of about 10^-2 for the tested resolutions, while closer gaps cause overlapping kernels and ill-conditioning.The available gaps are also constrained by integer hydrodynamic-cell counts for a fixed cylinder model and marker spacing.Source: E. Steady Stokes flow around cylinders · p. 43
  • The normalized drag agrees very well with boundary-integral results across a broad range of volume fractions, and the close-packed results match the lubrication prediction over an increase exceeding six orders of magnitude.The IB drag does not truly diverge as the inter-cylinder gap tends to zero because of kernel regularization.Source: E. Steady Stokes flow around cylinders · pp. 43–44, Figure 10 · p. 43

F. Unsteady flow around cylinders

  • The unsteady-cylinder study focuses on steady flow at Reynolds numbers no greater than 100.The passage notes that increasing Reynolds number brings unsteadiness and vortex shedding, followed eventually by transition to three-dimensional flow.Source: F. Unsteady flow around cylinders · p. 44

1. Drag on periodic array of disks

The numerical tests examine drag and flow around periodic arrays of disks, cylinders, and spheres. They report agreement with theory or prior computations, while identifying resolution and body representation as relevant qualifications.

1. Drag on periodic array of disks

  • Using k0 = 49.2 and k2 = 0.24, the theoretical formula (27) agrees very well with our results and with k0 = 51.2 and k2 = 0.26 given in the caption of Fig. 1 in Ref..Source: §1 · p. 46

2. Flow past a periodic column of cylinders

  • Computed drag coefficients and wake lengths for cylinders at Re = 10 agree well with prior Lattice-Boltzmann and finite-difference results.The wake length is measured from the cylinder center to the stagnation point; the empty-shell model also shows an interior counter-rotating vortex pair that diminishes under refinement.Source: §2 · p. 47, Table V · p. 48

G. Flow past periodic arrays of spheres

  • Sphere-array drag shows the expected switch from quadratic to linear dependence on Re, with reasonable agreement to literature data that improves with resolution.At close packing, filled and empty-shell models differ substantially; discrete sphere asymmetry can produce a small O(Re) correction that may dominate the true quadratic correction.Source: G. Flow past periodic arrays of spheres · pp. 48–49, Figure 13 · p. 49

VIII. CONCLUSIONS

The conclusion presents the constrained IB method as applicable to rigid-body flows at zero and finite Reynolds numbers, while noting marker-spacing and scalability challenges. It identifies extensions to finite-element force representations and freely moving bodies as future directions.

VIII. CONCLUSIONS

  • The method enforces strict rigidity at zero and finite Reynolds numbers and shows good agreement with theory and literature values across standard two- and three-dimensional tests.The authors report that marker spacing should be approximately two fluid-grid cells for accurate and stable pointwise traction estimates.Source: §VIII · p. 50
  • Many-marker use remains computationally challenging: dense mobility solves have O(N^2) memory and O(N^3) time complexity, motivating fast-solver development.The inverse fast multipole method has preliminary results for systems with as many as 5 · 10^5 markers, but integration and the effects of approximate inner solves require further study.Source: §VIII · pp. 50–51
  • Extending the method to large numbers of freely moving rigid bodies requires a preconditioner for the enlarged saddle-point system and its additional Schur complement.The paper describes this as a first step toward that goal.Source: §VIII · p. 51
  • An immersed-boundary/finite-element force representation may reduce mobility conditioning sensitivity to marker spacing and yield smaller systems, but remains future work.The authors expect smoother force representations to improve accuracy.Source: §VIII · p. 52

Appendix

The appendix fits approximations to mobility functions in three dimensions for steady Stokes and finite-Reynolds-number cases. The fits show reasonable agreement with numerical data and relatively small scatter associated with translational and rotational invariance.

a. Steady Stokes flow

  • For steady Stokes flow, the appendix fits normalized mobility functions using numerical data and known small- and large-distance behavior.The fits are compared with the Rotne–Prager–Yamakawa tensor approximation; the empirical fits match the data better.Source: a. Steady Stokes flow · pp. 53–54

Appendix A: Approximating the mobility in three dimensions

  • The empirical mobility fits show small scatter around the numerical data for steady Stokes and finite-β cases, indicating approximate translational and rotational invariance.For finite β, the authors report good agreement between the numerical data and fits.Source: a. Steady Stokes flow · p. 54, b. Nonzero Reynolds numbers · p. 55

b. Nonzero Reynolds numbers

  • For finite β, the fits treat self-mobility separately from off-diagonal mobility and use an empirical trial-and-error formula that is not constrained to guarantee an SPD approximation.The coefficients are fitted at selected β values, interpolated for intermediate values, and β > 1000 uses the steady Stokes fit.Source: b. Nonzero Reynolds numbers · pp. 54–55

Appendix B: Approximating the mobility in two dimensions

The two-dimensional mobility approximations are built around physical asymptotics and empirical fits. Periodic boundaries complicate steady-Stokes behavior, and the numerical data show boundary-sensitive scatter at large distances.

2. Empirical fits

  • Empirical fitting formulas use known short- and long-distance asymptotics, with coefficients fitted at β values from 0 to 10 and linear interpolation between them.Larger β values use the steady-Stokes fit.Source: Appendix B · p. 56, §2 · pp. 57, 59

1. Physical constraints

  • In two dimensions, steady-Stokes flow has a logarithmically growing Green’s function, so the divergent part must be subtracted and the large-distance behavior cannot be specified for an infinite system.The viscous corrections are expected to decay exponentially with the physical scale h√β.Source: §1 · pp. 56–57
  • For inviscid flow, the physical constraint uses the field of a dipole in two dimensions.The fitting approach also accounts for the logarithmic large-distance behavior of the two-dimensional Stokes Green’s function.Source: §1 · p. 56

2. Empirical fits

  • Steady-Stokes mobility in a periodic two-dimensional domain is boundary-sensitive, and its diagonal value f∞(0) diverges logarithmically with system size, so it is treated as an input parameter.The effective hydrodynamic radius is defined through this diagonal value.Source: §2 · p. 57
  • The two-dimensional numerical data follow the theoretical asymptotic behavior, but show substantial boundary-driven scatter at larger distances.Figure 15 compares the empirical fits and numerical data for the 6-point kernel.Source: §2 · p. 57, Figure 15 · p. 58

Appendix C: Stokes flow between two concentric spheres

This appendix considers steady Stokes flow between concentric spherical shells. It describes the shell-motion setup and the resulting velocity and traction fields.

Appendix C: Stokes flow between two concentric spheres

  • The test configuration places a stationary inner sphere of radius a inside a concentric outer shell of radius b = a/λ moving with velocity V.The viscosity is set to unity.Source: Appendix C · p. 59
  • The velocity between the shells is obtained from Brenner’s expressions for this geometry.The pressure’s mean is set to zero to remove its null mode.Source: Appendix C · p. 60
  • The inner-sphere traction is evaluated as the stress jump across the inner shell, using spherical coordinates aligned with the flow.The coordinate directions are specified by the radial and polar unit vectors.Source: Appendix C · p. 60

Appendix D: Imposing Physical Boundary Conditions

The boundary treatment extends the domain with ghost cells so standard immersed-boundary interpolation and spreading can be used near physical boundaries. Its spreading operator is chosen to preserve energy conservation for homogeneous boundary conditions.

Appendix D: Imposing Physical Boundary Conditions

  • Near physical boundaries, the method modifies interpolation and spreading to ensure interpolated velocity vanishes at a no-slip boundary and marker mobility decreases monotonically toward it.The chosen approach is Yeo and Maxey’s, selected for simplicity and kernel independence.Source: Appendix D · p. 61
  • The method fills ghost cells, interpolates as usual in the extended domain, and defines spreading as the adjoint operation that accumulates ghost-cell contributions back into the interior.The adjoint choice ensures energy conservation for homogeneous boundary conditions.Source: Appendix D · p. 61
  • The same extension operator handles boundary conditions in both the fluid solver and near-boundary interpolation and spreading, simplifying implementation without lowering the fluid solver’s second-order accuracy.For homogeneous Dirichlet no-slip boundaries, ghost values use mirror inversion.Source: Appendix D · p. 62
Loading 1505.07865v2…