Source-linked AI summary
A Stochastic Immersed Boundary Method for Fluid-Structure Dynamics at Microscopic Length Scales
P. J. Atzberger, P. R. Kramer, C. S. Peskin
TL;DR
The paper addresses how to incorporate thermal fluctuations into immersed-boundary simulations of flexible structures and fluids across widely separated time scales. It derives a stochastic immersed-boundary method that preserves fluid–particle statistics over long steps, including when fluid modes are underresolved. The method reproduces Boltzmann equilibrium statistics, correct three-dimensional diffusion scaling, and a long-time τ^-3/2 velocity-autocorrelation decay.
Problem
Microscopic fluid–structure systems combine thermal fluctuations with stiff, widely separated fluid and structural time scales, making conventional stochastic discretization difficult.
Method
The paper derives thermal fluctuations for the time-dependent Stokes immersed-boundary equations and introduces a stochastic discretization that preserves fluid–particle statistics and correlations over long time steps.
Results
The method reproduces correct Boltzmann equilibrium statistics, three-dimensional diffusion scaling, and a long-time τ^-3/2 algebraic decay in immersed-particle velocity autocorrelation.
Takeaways & Limitations
The physical checks indicate that the stochastic immersed boundary method can model coarse-scale thermal fluctuations in fluid–structure systems relevant to cellular and intracellular mechanics.
Takeaways & Limitations
The τ^-3/2 autocorrelation scaling is restricted to τ ≪ ρL^2/µ; at much longer times, the correlation decays exponentially because of the lowest fluid modes.
Abstract
from arXiv · showhide
In this work it is shown how the immersed boundary method of (Peskin2002) for modeling flexible structures immersed in a fluid can be extended to include thermal fluctuations. A stochastic numerical method is proposed which deals with stiffness in the system of equations by handling systematically the statistical contributions of the fastest dynamics of the fluid and immersed structures over long time steps. An important feature of the numerical method is that time steps can be taken in which the degrees of freedom of the fluid are completely underresolved, partially resolved, or fully resolved while retaining a good level of accuracy. Error estimates in each of these regimes are given for the method. A number of theoretical and numerical checks are furthermore performed to assess its physical fidelity. For a conservative force, the method is found to simulate particles with the correct Boltzmann equilibrium statistics. It is shown in three dimensions that the diffusion of immersed particles simulated with the method has the correct scaling in the physical parameters. The method is also shown to reproduce a well-known hydrodynamic effect of a Brownian particle in which the velocity autocorrelation function exhibits an algebraic tau^(-3/2) decay for long times. A few preliminary results are presented for more complex systems which demonstrate some potential application areas of the method.
1. Introduction.
The paper extends the immersed boundary method to thermally fluctuating microscopic fluid–structure systems, where broad spatial and temporal scales make detailed simulation difficult. Its stochastic method preserves relevant fluid contributions across resolved and underresolved regimes while reproducing key physical behavior.
- Background: The immersed boundary method provides a physically direct framework for coupling flexible structures with fluid dynamics at cellular and subcellular scales.Its applications include blood flow around heart valves, inner-ear wave propagation, and insect-flight lift generation.
- Motivation: Microscopic cellular systems require modeling thermal fluctuations and interactions among fluids, membranes, vesicles, polymers, and molecular motors.These phenomena span length scales from tens of microns to tens of nanometers, where thermal fluctuations can be significant.
- Motivation: Molecular-detail simulation is infeasible for complex cellular systems because resolving their broad active length and time scales is computationally costly.The paper therefore uses a coarse-grained approach intended to capture the most relevant dynamical features.
- Contribution: The method retains fluid dynamics, allowing thermally fluctuating dynamics to include subtle inertial effects and preserve the topology of flexible structures.The latter property supports simulations of polymers that do not cross themselves or each other, including links and knots.
- Contribution: The proposed stochastic discretization maintains good statistical accuracy when fluid modes are completely underresolved, partially resolved, or fully resolved.Only the immersed-structure degrees of freedom constrain the time step in the stated error regime.
- Validation: Theoretical and numerical checks indicate correct Boltzmann statistics, diffusion scaling, and long-time velocity autocorrelation behavior for thermally fluctuating immersed systems.The paper presents the stochastic immersed boundary method as a viable coarse-scale approach for cellular and subcellular biological phenomena.
2. Fluid-Particle Equations.
The fluid-particle equations model immersed structures as force-coupled components of an incompressible, time-dependent fluid. Thermal and structural forces enter the total fluid force density, with smoothed particle representations coupling force application and velocity interpolation.
- Fluid dynamics: At small Reynolds number, nonlinear advection is neglected while the fluid time derivative is retained to represent fast Brownian and structural timescales.The fluid is treated as Stokesian in its spatial regime, but not quasistatic in time.
- Force coupling: The total fluid force density combines forces from immersed structures with thermal fluctuations.Structural forces may arise from elastic deformation or external application and are transmitted to the fluid through the immersed structures.
- Immersed structures: Flexible structures and particles are represented as collections of M discrete elementary particles with positions and force laws determined by their structural properties.For the principal formulation, these forces derive from a conservative potential V({X}).
- Coupling operators: A smoothed delta function δa represents each elementary particle and maps particle forces to localized fluid force density.The same representation is used to interpolate nearby fluid velocity back to the elementary particle.
- Modeling assumption: The parameter a is treated as a physical elementary-particle size rather than a numerical parameter that vanishes with fluid-grid refinement.This distinguishes the smoothed-delta construction from its standard immersed-boundary use.
3.1. Summary of the Numerical Method.
The numerical method discretizes the coupled fluid–structure equations on a periodic grid and updates structural forces, fluid Fourier modes, and particle positions sequentially. Stochastic increments are generated to preserve the required thermal statistics and correlations.
- Discretization: The scheme uses finite differences on a periodic N-by-N-by-N fluid grid with spacing Δx = L/N.Fluid velocity and pressure are represented at grid points, while Fourier transforms connect physical and spectral representations.
- Time-step update: Each time step Δt first computes structural forces and their Fourier-transformed force density, then updates the fluid velocity in Fourier space.The fluid recurrence includes viscous dissipation, structural forcing, and thermal forcing.
- Thermal forcing: Thermal fluctuations are represented by Gaussian random variables whose variances are selected for the stochastic fluid-mode updates.The procedure also generates time-integrated fluid velocities for updating immersed-particle positions.
- Particle update: Particle positions are updated using time-integrated fluid velocities generated with correlations consistent with the previously computed fluid state.This construction links particle motion to the stochastic fluid update rather than treating their increments independently.
- Computational cost: The computational cost excluding application-specific forces is dominated by FFT and IFFT operations, requiring O(N^3 log(N)) arithmetic steps in three dimensions.The Fourier representation is therefore central to the method’s implementation.
3.2. Heuristic Discussion of the Numerical Method.
The method integrates fluid modes over each time step while preserving their effects on immersed structures, allowing fast fluid dynamics to be resolved or underresolved. Its stochastic increments maintain equilibrium variance and fluid–structure correlations across these regimes.
- Fluid update: The scheme’s Fourier-space fluid update integrates structural and thermal forces over a time step designed to tolerate partially resolved or underresolved fluid dynamics.Exponential factors encode the time-step dependence of viscous mode relaxation.
- Thermal increments: For a fluid mode with relaxation time 1/αk, the thermal velocity-increment variance approaches kBT/(ρL^3) for long time steps.For short time steps, the variance is proportional to Δt, giving increment magnitude proportional to √Δt.
- Incompressibility: The projection ℘⊥k enforces incompressibility in the structural and thermal contributions to the fluid update.The mode-dependent distinction in Dk is described as a technical consequence of the discrete Fourier transform.
- Particle motion: Particle positions are updated using time-integrated fluid velocities generated so that immersed structures retain the correct correlations with previously computed fluid velocities.The construction assumes structural forces remain constant over the time step.
- Derivation: The stochastic forcing is derived after spatial discretization so the semi-discrete system remains consistent with equilibrium statistical mechanics.The derivation separates spatial discretization from the continuous-time stochastic formulation before time discretization.
3.3. Derivation of the Numerical Method.
The method discretizes the incompressible fluid and its thermal forcing in Fourier space, enforcing real-valuedness and incompressibility while integrating fluid modes analytically over each time step. This yields a scheme whose time step is constrained by immersed-structure dynamics rather than fluid-mode time scales.
- The fluid equations are spatially discretized by finite differences and coupled to immersed particles through fluid-particle coupling equations.
- The particle velocity averages the fluid velocity over a finite-width region, keeping it finite in the continuum limit despite divergent pointwise thermal fluctuations.
- The same averaging weight is used to spread force to the fluid, ensuring energy conservation in the fluid-particle interaction.
- Fourier constraints: Fourier-mode projections enforce incompressibility by restricting nonzero modes to the subspace orthogonal to the wavevector, while special modes obey real-valuedness constraints.
- Thermal forcing: Thermal forcing is modeled as Gaussian white noise in Fourier space, with complex Brownian motions constrained to produce a real-valued velocity field.
- Thermal forcing: The fluctuation-dissipation relation sets thermal forcing strengths so constrained fluid modes reproduce Boltzmann statistics and equipartition.
- Thermal forcing: The zero Fourier mode is not thermally forced because it represents undamped whole-fluid translation and internal thermal fluctuations conserve total momentum.
- Time integration: Analytically integrating the fluid-mode expression permits fully resolved, partially resolved, or completely underresolved fluid modes, provided the time step is small relative to immersed-structure time scales.
4. Accuracy of the Method.
The method’s error is analyzed across fully resolved, fully underresolved, and partially resolved fluid dynamics. Formal estimates and simulations indicate accuracy for time steps below the immersed-structure motion timescale, even when fluid modes are underresolved.
- 4. Accuracy of the Method: The analysis considers fully resolved, fully underresolved, and partially resolved fluid-dynamics regimes.Each regime receives asymptotic error estimates for the particle, fluid-mode, and velocity-field discretizations.
- 4. Accuracy of the Method: The scheme preserves statistical contributions and correlations from underresolved fluid dynamics over time steps limited primarily by immersed-structure timescales.This treatment keeps local time-discretization error small even when the time step exceeds the fastest fluid-mode timescales.
- 4.1. Error Estimates for Time Steps which Fully Resolve the Fluid Dynamics: The method has strong first-order accuracy as the time step becomes small, although its particle-count error scaling is a pessimistic worst-case estimate.The worst case assumes all elementary particles are clustered; actual errors may be smaller depending on force interactions and clustering.
- 4.1. Error Estimates for Time Steps which Fully Resolve the Fluid Dynamics: In the force-free case, fluid modes are simulated exactly, while elementary-particle dynamics retain strong first-order temporal accuracy.The exact fluid-mode treatment avoids the finite-difference error contribution that grows when time steps exceed fluid-mode timescales.
- 4.2. Error Estimates for Time Steps which Underresolve All Fluid Modes: Underresolved-regime estimates can be smaller than fully resolved-regime extrapolations, so apparent lower powers of ∆t do not necessarily indicate worse accuracy.The relevant error-to-system-change ratios are much smaller than one in the underresolved asymptotic regime.
- 4.3. Error Estimates for Time Steps which Underresolve Only Some Fluid Modes: The method remains theoretically accurate for all ∆t ≪τmov(a), regardless of how completely the fluid dynamics are resolved.This includes transitions between resolved, underresolved, and partially resolved regimes.
5. Physical Behavior of the Method and Numerical Results.
The method is evaluated through diffusion, equilibrium, and autocorrelation analyses, showing physically correct diffusion scaling and accurate behavior even with time steps that underresolve fast fluid modes.
- 5.1. Diffusion of Immersed Particles.: The diffusion coefficient is derived from the fluid-velocity autocorrelation function and exhibits physically correct scaling with the relevant physical parameters.The estimate uses the Kubo formula and the Fourier representation of the velocity autocorrelation.
- 5.1. Diffusion of Immersed Particles.: Numerical simulations with particle sizes a = 1, 2, 3, 4, 5 compare diffusion estimates against the analytic prediction for long time steps that underresolve fast fluid modes.Each estimate used n = 10^4 sampled trajectories with ∆t = 10^3ns and t1 = 10^4ns.
- 5.1. Diffusion of Immersed Particles.: The simulated particle diffusion coefficient agrees well with the theoretical estimate 5.2.The comparison is made using numerical estimates and one-standard-deviation error bars.
- 5.1.1. Derivation of the Diffusion Coefficient.: The numerical method preserves the velocity-autocorrelation structure for finite time steps when ∆t ≪ τdiff(a), without requiring ∆t ≪ 1/αk.This permits accurate treatment of fluid dynamics over time steps longer than the fastest fluid-mode timescales.
5.2. Algebraic Decay of Velocity Autocorrelation Function.
The method reproduces the hydrodynamic algebraic decay of an immersed particle’s velocity autocorrelation over an intermediate time interval, while finite system size causes exponential decay at very long times.
- 5.2. Algebraic Decay of Velocity Autocorrelation Function.: The scaling arises from fluid-particle coupling, in which fluid momentum retains memory of the particle’s recent motion.The immersed boundary representation has slightly different constant prefactors because particles are represented differently from the physical model.
- 5.2. Algebraic Decay of Velocity Autocorrelation Function.: The algebraic regime is restricted to τ ≪ ρL^2/µ because of the finite size of the computational domain.For τ ≫ ρL^2/µ, the correlation decays exponentially according to the lowest-wavenumber fluid modes.
- 5.3. Equilibrium Statistics of Immersed Particles.: The equilibrium test compares the numerical radial distribution of a confined immersed particle with the corresponding Boltzmann distribution.The potential uses R1 = 125nm, R2 = 250nm, and c = 6kBT/(R2 − R1).
5.4. Osmotic Pressure of Confined Non-interacting Particles.
The section relates osmotic pressure to particle confinement, concentration, and thermal fluctuations. For non-interacting particles, the stochastic immersed boundary method recovers van’t Hoff’s law and supports pressure estimates from fluid or wall forces.
- Osmotic pressure arises when particles are confined by a boundary permeable to fluid but less permeable to particles.
- At equilibrium, van’t Hoff’s law relates osmotic pressure to the concentration of confined particles.
- The method predicts fluid pressure by averaging thermal fluctuations, particle forces, and pressure under Boltzmann statistics.
- The simulated fluid pressure obeys the local formulation of van’t Hoff’s law for particles in a conservative potential.
- Osmotic pressure can also be computed from the average force per unit area exerted by solute particles on the chamber wall.
- For hard-wall confinement, fluid-pressure and wall-force formulas agree up to a difference that vanishes as the boundary layer narrows; soft walls can produce different pressures.
5.5. Application: Simulation of Interacting Immersed Particles and Os-
The section applies the stochastic immersed boundary method to interacting particles, polymers, and motor-cargo systems. It shows how coupling, topology, and hydrodynamic loading alter osmotic pressure, concentration, and transport.
- The simulations vary osmotic pressure for particle dimers confined in a microscopic spherical chamber as binding strength changes.
- Strongly coupled dimers behave effectively as single particles, producing a pressure approximately half that of unbound monomers in accordance with van’t Hoff’s law.
- As dimer coupling increases, particle density near the confining region decreases and osmotic pressure drops nonlinearly.
- At intermediate coupling, dimer length scales comparable to the chamber and wall thickness produce pressures not well described by van’t Hoff’s law.
- The method preserves fluid-mediated velocity correlations between nearby immersed structures and supports simulations of polymer chains and knots.
- Increasing polymer knottedness restricts accessible configurations, reduces boundary-layer occupancy, and significantly lowers osmotic pressure.
- Large opposing fluid flows increase cargo drag and can nearly stall motor-protein transport.
6. Conclusion and Discussion.
The paper introduces a stochastic immersed boundary method that incorporates thermal fluctuations while permitting time steps with unresolved, partially resolved, or fully resolved fast dynamics. Theoretical and numerical checks support its physical fidelity and potential for biological fluid-structure simulations.
- The method systematically accounts for fast fluid-particle dynamics across underresolved, partially resolved, and fully resolved time-step regimes.
- Three-dimensional simulations reproduce the correct physical-parameter scaling of immersed-particle mean squared displacement.
- The method captures the algebraic τ^-3/2 decay of a Brownian particle’s velocity autocorrelation at long times.
- These checks indicate potential for modeling thermally fluctuating immersed structures in cellular and intracellular biological systems.
Appendix A. The Representation Function δa for Immersed Particles.
The appendix specifies the representation function δa used for immersed particles and describes its Fourier-space treatment for numerical analysis. Particle sizes are restricted to integer multiples of the lattice spacing to preserve numerical properties.
- The immersed boundary method represents elementary particles using a specified function δa.
- For three-dimensional systems, δa is defined componentwise for particles of size a.
- Particle sizes are restricted to a = n∆x, where n is a positive integer, to maintain good numerical properties.
- The appendix analyzes δa through its discrete Fourier coefficients, including their dependence on particle position relative to the lattice.
- The appendix derives the fluid-velocity autocorrelation by representing velocity modes in Fourier space and applying stochastic calculus.
Appendix D. Constants: Accuracy and Error Estimates.
The nondimensional Q factors in the Section 4 error estimates are approximately independent of physical parameters, but can be evaluated for specific simulation parameters.
- The nondimensional Q factors are approximately independent of the physical parameters.
- Specific physical parameters are used to compute estimated Q values for comparison between theoretical error estimates and numerical simulations.
- The parameter-specific Q evaluations use the system parameters listed in Table 4.2.
Tables.
The tables organize the method parameters, numerical-simulation values, and notation conventions used in the paper.
- Table 4.1 lists parameters of the method.
- Table 4.2 lists values used in numerical simulations.
- Table 4.3 defines notation conventions.