Source-linked AI summary
Accurate calculation of Stokes drag for point-particle tracking in two-way coupled flows
Jeremy Horwitz, Ali Mani
TL;DR
Two-way coupling disturbs the fluid velocity used to calculate Stokes drag, leaving the undisturbed value unavailable at the particle location. The paper estimates drag error and corrects the velocity estimate; in settling tests, the correction reduces time-averaged steady-state settling-velocity error by at least an order of magnitude versus standard interpolation schemes.
Problem
In two-way-coupled point-particle simulations, the disturbed fluid velocity is often substituted for the undisturbed velocity required by Stokes drag.
Source: §1 · pp. 2–3Method
The paper estimates drag error as particle size relative to grid spacing changes and corrects the undisturbed-velocity estimate using the computed disturbed field.
Source: §1 · p. 4, §4 · pp. 11–12Results
At least one order of magnitude smaller time-averaged steady-state settling-velocity error was achieved by the correction scheme than by trilinear and Lagrange-4 schemes.
Source: §6 · p. 32Takeaways & Limitations
Uncorrected schemes underpredict momentum exchange between particles and fluid because they use disturbed rather than undisturbed velocity to calculate Stokes drag.
Source: §6 · p. 32At early times, especially for St = 0.25, the correction can overestimate drag error and underestimate settling velocity because it is tuned to steady flow past a sphere.
Source: §5 · p. 27
Abstract
from arXiv · showhide
In this work, we propose and test a method for calculating Stokes drag applicable to particle-laden fluid flows where two-way momentum coupling is important. In the point-particle formulation, particle dynamics are coupled to fluid dynamics via a source term that appears in the respective momentum equations. When the particle Reynolds number is small and the particle diameter is smaller than the fluid scales, it is common to approximate the momentum coupling source term as the Stokes drag. The Stokes drag force depends on the difference between the undisturbed fluid velocity evaluated at the particle location, and the particle velocity. However, owing to two-way coupling, the fluid velocity is modified in the neighborhood of a particle, relative to its undisturbed value. This causes the computed Stokes drag force to be underestimated in two-way coupled point-particle simulations. We develop estimates for the drag force error as function of the particle size relative to the grid size. We then develop a correction method that estimates the undisturbed fluid velocity from the computed disturbed velocity field. The correction scheme is tested for a particle settling in an otherwise quiescent fluid and is found to reduce the error in computed settling velocity by an order of magnitude compared with common interpolation schemes.
1. Introduction
Point-particle methods model particle–fluid momentum exchange without resolving particle-scale fluid stresses, but two-way coupling disturbs the fluid velocity used to calculate Stokes drag. The paper develops drag-error estimates and a correction for the undisturbed velocity, then tests it against a settling-particle solution.
1. Introduction
- Point-particle simulations represent particles as momentum and energy sources in the Eulerian fluid, with equal and opposite forces coupling the phases.Fully resolving particle boundary layers is impractical for simulations containing many particles.Source: §1 · p. 2
- In two-way coupling, using the disturbed fluid velocity instead of the undisturbed velocity can significantly miscalculate Stokes drag.Stokes drag depends on the difference between particle velocity and undisturbed fluid velocity at the particle location.Source: §1 · p. 3
- The disturbed–undisturbed velocity difference is expected to scale as O(d_p/d_x) for small particles, and the paper develops drag-error estimates as this ratio increases.The proposed correction estimates undisturbed velocity from the computed disturbed field.Source: §1 · pp. 3–4
- The corrected scheme agrees very well with the analytical settling solution, while uncorrected schemes produce very large errors.Higher-order interpolation can perform worse because it estimates the disturbed velocity more accurately, which is a poorer estimate of the undisturbed velocity.Source: §1 · p. 4
- The correction is intended for regimes where Stokes drag, and possibly gravity, are physically appropriate; assumptions include low particle Reynolds number and small volume fraction.The authors also identify small particle size relative to the Kolmogorov scale, high particle-to-fluid density ratio, and low particle shear Reynolds number as relevant conditions.Source: §1 · pp. 4–5
2. Point-Particle Methodology
The methodology uses DNS with two-way momentum coupling and treats calculating undisturbed velocity as an estimation problem rather than ordinary interpolation. It also relates drag error to particle size relative to grid spacing and requires consistent interpolation and force projection for the discrete energy balance.
2. Point-Particle Methodology
- The study assumes the unladen fluid is well resolved and restricts its analysis to DNS, attributing drag error to estimating undisturbed velocity.The method uses a second-order finite-difference momentum solver on a staggered mesh and explicit fourth-order Runge–Kutta time stepping.Source: §2 · p. 5
- In two-way coupling, calculating undisturbed particle velocity is an estimation problem; without particle feedback, it is an interpolation problem.Particle and fluid data occupy Lagrangian and discrete Eulerian locations, respectively.Source: §2 · p. 6
- The error associated with using disturbed velocity may be small when particles are much smaller than grid spacing, and the ratio Λ = d_p/d_x is a key parameter.Grid refinement or nonuniform grids can increase Λ and make uncorrected Stokes-drag errors significant.Source: §2 · p. 7
- Interpolation and projection should use the same stencil weights for a consistent discrete energy equation; differing weights introduce an additional numerical dissipation error.Point-particle simulations do not resolve particle boundary layers, so their true fluid dissipation rate is not accurately calculated; matching total dissipation to fully resolved simulations remains an open question.Source: §2 · pp. 7–8
3. Error Estimation
The error-estimation tests show that using the disturbed velocity underestimates Stokes drag, with error increasing as particle size grows relative to grid spacing. The estimates use a stationary particle in uniform flow and distinguish the effects of projection from the velocity error itself.
3. Error Estimation
- The estimate uses uniform low-Reynolds-number flow past a stationary particle held at a cell face, eliminating interpolation error while its projected force slows nearby fluid.The tested particle has Rep = 0.1, with the drag evaluated from the disturbed velocity rather than the undisturbed inflow value.Source: §3 · pp. 8–9, Figure 2 · p. 8
- The symmetric projection places the entire force on one grid point, while Gaussian mollification spreads it across neighboring points to weaken the local disturbance.The Gaussian filter is normalized so its weights sum to one; its cutoff exceeds 1.5 grid spacings.Source: §3 · p. 9
- Drag-force error grows as O(Λ) and remains non-negligible at Λ = 0.1; Gaussian projection reduces it only modestly and does not ensure convergence to Stokes drag.Figure 3 reports steady errors for one-point and Gaussian projection schemes.Source: §3 · p. 9, Figure 3 · p. 10
- The reported values are steady, collocated-particle error estimates; moving-particle errors may differ with sampled positions, projection, and Stokes number.The simulations integrate for about 100 viscous relaxation times; spreading force over many grid points could reduce error but would be impractical for millions of particles.Source: Figure 3 · p. 10, §3 · p. 10
4. Proposed Correction Method
The proposed correction estimates undisturbed velocity by adding a calibrated correction for the particle’s own disturbance to the computed velocity. It represents that correction with a time-independent C field, then interpolates the field and velocity information for particles at arbitrary positions in a grid cell.
4. Proposed Correction Method
- The method estimates undisturbed velocity with a truncated expansion that corrects the computed disturbed velocity for the disturbance generated by the particle itself.This is distinct from a Faxén correction for disturbances created by other particles, and is expected to work when particle interactions are not leading-order effects.Source: §4 · p. 11
- Interpolating C is simpler than interpolating fluid velocity because C depends on particle size and position but is time-independent and uncontaminated by particle disturbance.The correction coefficient is calibrated by reverse engineering a fixed particle in uniform flow for each particle-to-grid size Λ.Source: §4 · p. 12
- For arbitrary particle positions, C values are determined within a grid cell; Stokes-flow symmetry reduces the eight vertices to six unique coefficients.The coefficients are tabulated for different Λ, and smooth cubic splines can estimate values for other particle sizes.Source: §4 · pp. 14–17, Table 2 · p. 16
- For Λ ≥ 0.25, the C field varies little, while its sensitivity to Λ is greater at smaller particle sizes.Errors in estimating C at small Λ are expected to matter less because those particles generate relatively small disturbance velocities.Source: §4 · p. 17
- The same six coefficients apply across grid cells and Cartesian velocity directions, with the C field transposed for the v and w components.The field is periodic and can be characterized by two layers.Source: §4 · p. 17, Figure 11 · p. 18
4.1 Algorithm Summary: to calculate undisturbed velocity for arbitrary particle location
- For arbitrary particle locations, the algorithm interpolates C and velocity information to estimate undisturbed velocity, then forms the particle drag force.The reported coefficients are consistent with trilinear interpolation for C and force projection; the velocity Laplacian is available at fluid cell faces for the momentum update.Source: §4.1 · pp. 18–19
5. Method Verification
The settling tests compare the correction scheme with trilinear and fourth-order Lagrange interpolation across particle sizes, Reynolds numbers, and Stokes numbers. The section also examines how the correction affects particle–fluid energy accounting and identifies boundaries for applying point-particle models.
5. Method Verification
- The verification uses a settling particle in an unbounded, quiescent fluid, with gravity tilted relative to the grid so its trajectory samples interpolation scenarios within mesh cells.The three-dimensional periodic box has 128^3 grid points; irrational ratios between gravity-vector components make the trajectory dense within each mesh element.Source: §5 · p. 20, Figure 12 · p. 20
- For Λ=1, the correction scheme’s maximum steady-state settling-velocity error was approximately 5%, compared with overshoots of approximately 75% for trilinear and 100% for Lagrange-4.The settling histories compare the schemes with the analytical solution across particle sizes.Source: §5 · pp. 22–23, Figure 13 · p. 22
- Higher-order interpolation can estimate the disturbed velocity more accurately yet estimate the undisturbed velocity less accurately; its projection weights can also intensify the nearby disturbance.The Lagrange-4 scheme assigns weights of 0.5625 to the two nearest grid points and -0.0625 to the next-nearest points.Source: §5 · pp. 23–24
- The correction performs robustly over the examined decade of Reynolds numbers, while low-Stokes-number cases show early settling-velocity undershoot, especially at St=0.25.For St=5 and 10, settling agrees excellently with the analytical solution; at low St, the steady-state-tuned correction can overestimate early drag, though settling approaches the analytical steady state.Source: §5 · pp. 24–27
- At least an order of magnitude less steady-state settling-velocity error was found with the correction scheme than with either interpolation scheme; the maximum was ~5.5%, versus ~76% for trilinear and ~102% for Lagrange-4.The reported maximum comparison is for Rep=0.05, St=10, and Λ=1.0; errors were ≤O(5%) across tested parameters.Source: §5 · p. 27
- For all tested parameters, normalized rms drift-velocity error was below 1% with the correction scheme; at St=0.25, Lagrange-4 had ~9.2% and trilinear almost 11% error.At St=10, drift was O(1%) for all schemes.Source: §5 · p. 28, Table 7 · p. 28, Table 8 · p. 28
5.1 Energetics
- In the energy accounting, a point-particle model’s artificial near-particle dissipation is offset at steady state by an opposite-sign term, while the first term represents analytical heat generation.The point-particle source term accounts for a portion of heat generation associated with unresolved particle boundary layers.Source: §5.1 · p. 30
- Point-force models necessarily produce artificial fluid structures near particles; extending similar corrections to other particle-force terms requires verification benchmarks for those more complex equations of motion.Added-mass and history terms also depend on undisturbed fluid information.Source: §5.1 · p. 31
6. Summary
The correction scheme estimates undisturbed fluid velocity from disturbed velocity data and improves settling-velocity accuracy over standard interpolation schemes. The summary also identifies the drag-error scaling, the scheme’s scope, and the implications of using disturbed velocity in two-way coupling.
6. Summary
- For uniform flow past a stationary point-particle, drag-force error is of the same order as particle size relative to grid spacing; wider force projection reduces but does not eliminate it.Only a projection kernel with infinite support could eliminate the error.Source: §6 · p. 31
- The correction scheme estimates undisturbed velocity from neighboring disturbed velocities using a non-convex interpolation with negative near-particle and positive farther-away weights.For arbitrary particle positions, six time-independent C constants characterize the scheme; the C field is not contaminated during simulation.Source: §6 · pp. 31–32
- Time-averaged steady-state settling-velocity error was at least one order of magnitude smaller with the correction scheme than with trilinear or Lagrange-4 interpolation.The comparison covered different particle sizes, Reynolds numbers, and Stokes numbers.Source: §6 · p. 32
- Uncorrected schemes under-predict particle-fluid momentum exchange because their disturbed-velocity slip is smaller than the slip based on undisturbed velocity.The authors suggest this may help elucidate disagreements between point-particle turbulence simulations and experiments on particle-driven turbulence modification.Source: §6 · p. 32
- The scheme is developed for momentum coupling in dilute particle-laden flows where Stokes drag is justified as the leading-order contribution.The presented implementation uses a uniform staggered grid; non-uniform grids require modification because C-field symmetry is partly broken.Source: §6 · p. 33
Appendix: Spline Interpolants for C-field
The appendix gives cubic-spline interpolants for the C field, with coefficients organized across six scalar values and five Λ intervals. It specifies the interval origins and ranges and identifies tables containing the C-field and spline coefficients.
Appendix: Spline Interpolants for C-field
- The C-field is smoothly interpolated with a cubic spline whose coefficients are indexed by six scalar values and five Λ intervals.The appendix identifies the interpolated field with Table 2 and lists coefficient tables for C1 through C6.Source: Appendix: Spline Interpolants for… · p. 37, Table A.1 · p. 37, Table A.2 · p. 38
- The five Λ intervals are [0.01,0.05], [0.05,0.1], [0.1,0.25], [0.25,0.5], and [0.5,1.0], with virtual origins [0.01, 0.05, 0.1, 0.25, 0.5].The appendix gives the interval-specific spline expression in equation (A.1).Source: Appendix: Spline Interpolants for… · p. 37
- Tables A.2–A.7 provide spline coefficients C_ijk for the six scalar values C1 through C6.Table A.1 lists C-field coefficients for different Λ values.Source: Appendix: Spline Interpolants for… · p. 37, Table A.1 · p. 37, Table A.2 · p. 38