PaperlayerView exact revision on arXiv

Source-linked AI summary

Hybrid finite difference/finite element immersed boundary method

Boyce E. Griffith, Xiaoyu Luo

arXiv:1612.05916v2math.NAcs.CE

TL;DR

The paper addresses the immersed-boundary method’s reliance on Lagrangian meshes finer than the Eulerian grid. It introduces finite-element structural discretization with dynamically selected interaction points and retains a finite-difference fluid solver, achieving effective coarse structural meshes across benchmark fluid–structure problems while identifying a weak-form limitation.

  • Problem

    Conventional immersed boundary methods often require Lagrangian meshes finer than the Eulerian grid to avoid leaks, limiting independent structural and fluid discretizations.

  • Method

    The method couples finite-element structural mechanics to a Cartesian-grid finite-difference fluid solver through quadrature-defined interaction points inside structural elements.

  • Results

    The scheme enables Lagrangian meshes at least four times coarser than the Eulerian grid without leaks and can reduce errors substantially, including by an order of magnitude for rigid structures.

  • Takeaways & Limitations

    The coupling approach supports standard structural and incompressible-flow discretizations across elastic, rigid, actively contracting, and idealized cardiac mechanics problems.

  • Takeaways & Limitations

    The partitioned formulation lacks a discrete power identity implying energy conservation during Lagrangian–Eulerian interaction.

Abstract

from arXiv · show

The immersed boundary method is an approach to fluid-structure interaction that uses a Lagrangian description of the structural deformations, stresses, and forces along with an Eulerian description of the momentum, viscosity, and incompressibility of the fluid-structure system. The original immersed boundary methods described immersed elastic structures using systems of flexible fibers, and even now, most immersed boundary methods still require Lagrangian meshes that are finer than the Eulerian grid. This work introduces a coupling scheme for the immersed boundary method to link the Lagrangian and Eulerian variables that facilitates independent spatial discretizations for the structure and background grid. This approach employs a finite element discretization of the structure while retaining a finite difference scheme for the Eulerian variables. We apply this method to benchmark problems involving elastic, rigid, and actively contracting structures, including an idealized model of the left ventricle of the heart. Our tests include cases in which, for a fixed Eulerian grid spacing, coarser Lagrangian structural meshes yield discretization errors that are as much as several orders of magnitude smaller than errors obtained using finer structural meshes. The Lagrangian-Eulerian coupling approach developed in this work enables the effective use of these coarse structural meshes with the immersed boundary method. This work also contrasts two different weak forms of the equations, one of which is demonstrated to be more effective for the coarse structural discretizations facilitated by our coupling approach.

1. INTRODUCTION

The paper extends immersed boundary fluid–structure interaction with finite-element structural discretizations while retaining a finite-difference Eulerian solver. Its coupling scheme uses dynamically selected interaction points so structural meshes can be chosen primarily for structural accuracy rather than coupling requirements.

  • Structural modeling: The approach addresses general finite-deformation structural models beyond conventional fiber-based elasticity descriptions.Fiber models represent response only along a single material direction, whereas finite-element models support more general structural mechanics.
  • Method: The method combines a Cartesian-grid finite-difference scheme for incompressible flow with nodal finite elements for structural mechanics.It treats flexible hyperelastic structures and rigid structures with approximately imposed rigidity constraints.
  • Lagrangian–Eulerian coupling: Dynamically selected quadrature interaction points decouple direct fluid coupling from structural mesh nodes, overcoming the classical requirement for fine Lagrangian meshes.The finite-element nodes control interaction-point positions and represent the deformation of all material points.
  • Applications and evaluation: The method is evaluated on elastic, rigid, and actively contracting benchmark structures, including an idealized left-ventricle model.The elastic tests compare two weak formulations suitable for standard nodal finite-element structural mechanics.
  • Results: Coarse Lagrangian meshes can reduce discretization errors by an order of magnitude for fixed Eulerian spacing while improving volume conservation over the IFE method.The approach also requires sufficiently dense interaction points to prevent leaks and can accommodate large structural deformations through dynamic quadrature.

2. CONTINUOUS FORMULATIONS

The continuous formulation couples Eulerian fluid variables with Lagrangian structural mechanics, then develops strong and weak forms for elastic bodies and a penalty treatment for rigid structures.

  • 2.1. Immersed elastic bodies: The immersed elastic-body formulation represents fluid momentum, velocity, and incompressibility in Eulerian form while describing structural deformation and elasticity in Lagrangian coordinates.The structure occupies χ(U,t) within the physical domain, while the fluid occupies the complement.
  • 2.1. Immersed elastic bodies: The structural stress is expressed through the first Piola–Kirchhoff tensor, deformation gradient, Jacobian, and hyperelastic strain-energy models, while stress may also be specified directly.The methodology supports active tension models that use a stress response without requiring an energy functional.
  • 2.1.1. Strong formulation: The strong form combines incompressible fluid momentum with Eulerian elastic forcing generated by volumetric internal forces and interfacial transmission forces.Delta-function coupling spreads these Lagrangian force densities to the Eulerian domain, while velocity interpolation determines structural motion and enforces no-slip and no-penetration through motion.
  • 2.1.1. Strong formulation: Interface traction continuity requires elastic and fluid stress discontinuities to balance, making the transmission force singular where structural stress ends at the fluid–solid interface.The resulting elastic force is variationally equivalent to the divergence of the elastic stress.
  • 2.1.2. Weak formulations: Two weak formulations are introduced for standard nodal C0 finite elements because the Eulerian equations remain discretized by finite differences rather than weakly formulated.The partitioned formulation separates internal and transmission forces, whereas the unified formulation uses one volumetric body force.
  • 2.1.2. Weak formulations: The partitioned weak form separates smooth internal forcing from singular transmission forcing, while finite-element basis functions regularize the transmission force when projected volumetrically.The unified alternative combines both effects into a single total elastic force density and is related to earlier immersed finite element and variational formulations.
  • 2.2. Immersed rigid structures: For rigid structures, the fully constrained equations introduce a Lagrange multiplier for rigidity alongside pressure for incompressibility, forming an extended saddle-point problem.The paper instead uses a penalty formulation for a stationary structure, with stiffness and damping parameters approximating the constraint.
  • 2.2. Immersed rigid structures: As the stiffness penalty κ increases, the rigid structure approaches its initial configuration, while damping contributes through the penalty force.This replaces the specialized solution techniques required by the fully constrained formulation.

3. SPATIAL DISCRETIZATION

The method combines staggered-grid finite differences for Eulerian variables with nodal finite elements for structural variables, coupling them through quadrature-based interaction operators.

  • 3.1 Eulerian discretization: The Eulerian variables use a staggered-grid finite difference discretization, placing pressure at cell centers and velocity and force components on corresponding cell edges.Divergence, gradient, and Laplace operators are discretized with standard second-order finite differences, while nonlinear advection uses a staggered-grid PPM variant.
  • 3.2 Lagrangian discretization: The immersed structure is represented with a triangulated Lagrangian finite element mesh whose deformation, stresses, and force densities are computed from FE basis functions.The formulation uses standard Galerkin discretizations and Gaussian quadrature for the resulting integrals.
  • 3.3 Lagrangian-Eulerian interaction: Forces are evaluated at quadrature-defined interaction points and then spread to the Eulerian grid with a smoothed delta function, rather than being applied directly from structural nodes.This force-prolongation construction uses FE interpolation between nodes and interaction points before applying the standard IB spreading operation.
  • 3.3 Lagrangian-Eulerian interaction: The velocity-restriction operator is chosen so that nodal motion approximates the L2 projection of the Eulerian IB velocity field and is adjoint to force spreading.The adjoint relation yields the power identity associated with energy conservation during Lagrangian-Eulerian interaction.

4. IMPLEMENTATION

The implementation is provided in IBAMR, an open-source C++ framework for immersed-boundary fluid-structure-interaction models.

  • 4. IMPLEMENTATION: IBAMR supports distributed-memory parallelism and adaptive mesh refinement through integrations with SAMRAI, PETSc, hypre, and libMesh.These libraries provide much of the framework's parallel, solver, mesh-refinement, and finite-element functionality.

5.1. Thick elliptical shell

The thick-shell benchmarks test static and dynamic anisotropic and orthotropic materials, showing designed convergence and effective accuracy with coarse structural meshes. The partitioned formulation particularly improves pressure accuracy for coarse meshes in the orthotropic case.

  • 5.1.1. Anisotropic shell: The static anisotropic shell achieves second-order velocity convergence in all norms, with pressure convergence of second order in L1, approximately 1.5 in L2, and first order in L∞.The pressure rates reflect its C0 but not C1 regularity.
  • 5.1.1. Anisotropic shell: The dynamic anisotropic shell shows essentially second-order velocity convergence, pressure rates of 2, 1.5, and 1 in L1, L2, and L∞, and nearly second-order deformation convergence in L1 and L2.Deformation convergence is between first and second order in L∞, and higher-order structural elements can improve its robustness.
  • 5.1.1. Anisotropic shell: Across static and dynamic anisotropic tests, virtually identical errors for Mfac = 1, 2, and 4 indicate no appreciable accuracy loss with coarse structural meshes.The results also suggest that coarse Lagrangian meshes do not cause leaks at fluid-structure interfaces.
  • 5.1.2. Orthotropic shell: The static orthotropic shell has first-order velocity convergence, pressure convergence of first order in L1 and 0.5 in L2, and no L∞ pressure convergence because interface pressure is discontinuous.The orthotropic material model behaves as a fiber-reinforced solid with circumferential and radial fiber families.
  • 5.1.2. Orthotropic shell: In the dynamic orthotropic test, unified and partitioned formulations have similar accuracy for u and χ, while the partitioned formulation substantially improves pressure accuracy for relatively coarse Lagrangian meshes.The pressure improvement appears also to improve volume conservation.

5.2. Soft elastic disc in lid driven cavity

The soft-disc cavity benchmark evaluates volume conservation and formulation behavior during near contact with a moving boundary. The partitioned formulation generally conserves volume better, especially with coarse structural meshes, while volume errors converge at first order.

  • 5.2. Soft elastic disc in lid driven cavity: The benchmark uses a soft neo-Hookean disc in a lid-driven cavity, where the flow brings the structure nearly into contact with the moving upper boundary.The immersed boundary formulation handles this near contact automatically using a modified kernel.
  • 5.2. Soft elastic disc in lid driven cavity: Pressure and viscous-stress discontinuities at fluid-structure interfaces limit expected convergence rates to no better than first order in this cavity problem.The formulation also does not employ a specialized treatment for the cavity flow’s corner singularities.
  • 5.2. Soft elastic disc in lid driven cavity: The partitioned formulation generally provides superior volume conservation, especially for relatively coarse Lagrangian meshes, and its performance is relatively insensitive to mesh coarseness.Both formulations produce volume errors converging to zero at essentially first order.
  • 5.2. Soft elastic disc in lid driven cavity: The method’s volume-conservation results compare favorably with IFE results reporting losses up to 20%, or approximately 2.5% after an IFE modification.These comparisons concern the same test without, or with, a volume-conservation algorithm.

5.3. Flow past a cylinder

The cylinder benchmark shows that the method agrees quantitatively with prior results, while structural mesh refinement at fixed Eulerian resolution can reduce accuracy. Kernel choice changes the acceptable structural spacing, with three-point kernels remaining accurate across the tested spacings.

  • The benchmark uses a six-level adaptively refined Cartesian grid and quadratic structural elements with spacing approximately Mfac∆xfinest.The cylinder is modeled as a thin circular interface at Re = 200.
  • The present results show excellent quantitative agreement with earlier experimental and computational results for lift, drag, and Strouhal number at Re = 200.
  • For four- and six-point kernels, coarser structural meshes produce more accurate lift, drag, and vortex-shedding dynamics under fixed Eulerian resolution.Under simultaneous Lagrangian and Eulerian refinement, all relative structural spacings converge to the same dynamics.
  • For N = 32, the six-point kernel gives erratic results at Mfac = 1, whereas coarser structural discretizations avoid this behavior.
  • The three-point kernel produces comparable results for all tested Mfac values, unlike the four- and six-point kernels.Two-point piecewise-linear kernels reportedly behave similarly to the three-point kernel.

5.4. Idealized model of left ventricular mechanics

The method is applied to passive inflation and active contraction of an idealized left ventricle using finite-element structural models on Cartesian fluid grids. Both cases agree closely with consensus structural-code results, while active contraction shows approximately first-order grid convergence and transmural fiber rotation produces localized midwall stretch.

  • Model and discretization: The ventricular reference geometry is a truncated ellipsoid, discretized with trilinear Q1 hexahedral elements on a 64 × 64 × 64 Cartesian grid.
  • Active contraction: Active contraction produces torsion, and fiber strain is mostly compressive but includes slight midwall stretching because fibers rotate across the wall.
  • Passive inflation: The passively inflated ventricular model agrees within 1% with other structural mechanics codes for longitudinal, circumferential, and radial strains.
  • Active contraction: The active-contraction test uses a fiber-reinforced material with 180° transmural fiber rotation, active contractile stress, and a 15 kPa endocardial pressure load.
  • Active contraction: The actively contracting ventricular model agrees within 1% with grid-converged consensus results from other structural codes.Displacement differences are 2.5% between N = 48 and 64 and 1.8% between N = 64 and 96.

6. DISCUSSION AND CONCLUSIONS

The paper presents a finite-element immersed boundary method that permits structural meshes independent of the Cartesian fluid grid and supports elastic, rigid, and actively contracting structures. Benchmarks show that coarse structural meshes can improve accuracy, although the partitioned formulation lacks a discrete power identity.

  • The method couples structured or unstructured finite-element structural discretizations to a Cartesian-grid finite-difference scheme for Eulerian equations.It is demonstrated for hyperelastic, rigid, and actively contractile immersed structures.
  • Structural meshes at least four times coarser than the background Eulerian grid can avoid leaks, while fine meshes remain necessary for fine geometric features.The authors also report improved volume conservation over the immersed finite element method for some very coarse discretizations.
  • The partitioned formulation lacks a discrete power identity implying energy conservation during Lagrangian-Eulerian interaction.The unified formulation satisfies this identity and may be needed for an unconditionally stable implicit time-stepping scheme.
  • For circular-cylinder flow, fine structural discretizations can produce spurious drag, and acceptable structural spacing depends on the regularized kernel.Three-point kernels perform well at Mfac = 1, whereas four- and six-point kernels perform poorly at that spacing.
  • Projection onto finite-element shape functions filters velocity fluctuations below the structural mesh width, reducing lift and drag errors for relatively coarse structural meshes.The authors interpret these artifacts as consequences of overly dense structural meshes rather than intrinsic limitations of the immersed boundary methodology.

A.1. Basic time-stepping scheme

The time-stepping scheme advances the immersed structure, incompressible velocity, and pressure through a sequence of coupled updates. The resulting Crank-Nicolson-type fluid solve is handled with FGMRES and multigrid-based preconditioning.

  • The scheme advances χ, u, and p from time level n using a time increment ∆t, beginning with a preliminary deformed-structure configuration.
  • The nonlinear advection term is computed with a PPM-type approximation, and the unified formulation uses an analogous time-stepping scheme.
  • The update requires solving a Crank-Nicolson-type discretization of the time-dependent incompressible Stokes equations.
  • The coupled linear system is solved with flexible GMRES using a pressure-free projection method and inexact multigrid subdomain solvers as a preconditioner.

A.2. Initial time step

The initial time step uses a two-step predictor-corrector method because the standard time-stepping scheme requires lagged velocity and pressure values that are unavailable initially.

  • A.2. Initial time step: The initial time step uses a two-step predictor-corrector method instead of the lagged-value time-stepping scheme.The standard scheme cannot be applied because it requires time step-lagged values of u and p.
  • A.2. Initial time step: Because the initial pressure is unavailable, the method initializes the pressure guess with p = 0 for p^{n+1}.
  • A.2. Initial time step: The predictor-corrector procedure modifies the subsequent expression by using A^{n+1} in place of A^n.
Loading 1612.05916v2…