Force estimation from PIV fields

This page explains the two vorticity-based force estimation methods implemented in dpivsoft/Postprocessing.py, with their FEM support machinery in dpivsoft/meshTools.py. Both follow the formulation of Martín-Alcántara & Fernández-Feria (2019), written here in dimensional form:

  • the vortical impulse method (§4) — needs only the velocity/vorticity field;

  • the projection method (§5) — needs, in addition, an auxiliary potential solved once by FEM on a mesh around the body (§6).

1. Why compute forces from vorticity?

The direct route to the force on a body — the momentum balance over a control volume — requires the pressure on its boundary, and PIV measures velocity, not pressure. Vorticity-based formulations rewrite the force so that pressure is eliminated entirely: the force becomes moments of the vorticity field and quantities derivable from the measured velocity, all computable from a time-resolved planar PIV measurement.

The price is sensitivity: vorticity is a derivative of the measured field, so these methods amplify measurement noise and demand care near boundaries. That trade-off is inherent, not an implementation detail.

2. Notation and inputs

  • 2-D velocity field \(u(x, y, t)\), \(v(x, y, t)\) on a regular grid, as a time series (arrays of shape (ny, nx, nt)), with the reference frame centred on the object. A frame acceleration \(\mathbf{a}(t)\) can be supplied for non-inertial (moving-body) frames.

  • Vorticity \(\omega = \partial v/\partial x - \partial u/\partial y\) (out-of-plane component).

  • Fluid density \(\rho\) and dynamic viscosity \(\mu\), constant.

  • Forces are per unit span (2-D), returned as time series of \((F_x, F_y)\).

3. Computing the vorticity field

Postprocessing.vorticity(x, y, u, v, method) offers five stencils (following Raffel et al.), differing in noise behaviour and footprint:

method

stencil

output size

centered

2-point centred differences

\((N_y{-}2) \times (N_x{-}2)\)

leastsq

4-point least-squares derivative

\((N_y{-}4) \times (N_x{-}4)\)

richardson

4-point Richardson extrapolation

\((N_y{-}4) \times (N_x{-}4)\)

circulation

8-point circulation loop / area

\((N_y{-}2) \times (N_x{-}2)\)

curl

np.gradient-based curl

\((N_y{-}2) \times (N_x{-}2)\)

Wider stencils (leastsq, circulation) average over more points and are less noise-sensitive; centered has the smallest footprint. All five agree on smooth fields (e.g. solid-body rotation \(u = -y,\ v = x \Rightarrow \omega = 2\)).

Near solid walls the interior stencils are invalid (they would reach across the body). When a body outline is given, Object() builds the body mask and an ordered boundary path, and walls_vorticity() replaces the vorticity along the wall with one-sided estimates. Both force methods apply this automatically when solid_points is provided.

4. Method 1 — the vortical impulse method (ImpulseMethod)

4.1 Theory

For a body immersed in a finite control volume \(V\) with open outer boundary, the force can be written from the first moment of vorticity (the vortical impulse) plus surface corrections (Wu, Lu & Zhuang 2007; formulation as in Martín-Alcántara & Fernández-Feria 2019):

\[ \mathbf{F} \;=\; \underbrace{-\rho \frac{d}{dt} \int_V \mathbf{x} \times \boldsymbol{\omega}\, dV}_{F_i\;\text{(impulse)}} \;+\; \underbrace{\rho \int_V \mathbf{u} \times \boldsymbol{\omega}\, dV}_{F_v\;\text{(vortex force)}} \;+\; \underbrace{\rho V_s\, \mathbf{a}}_{F_{sol}\;\text{(solid volume)}} \;+\; F_{oe} \;+\; F_\mu \]

  • \(F_i\) — rate of change of the vortical impulse \(\rho \int \mathbf{x} \times \boldsymbol{\omega}\, dV\). In an unbounded domain this term alone gives the force; the rest corrects for the finite field of view.

  • \(F_v\) — the volume vortex force; in 2-D components \(\rho \iint (\omega v,\; -\omega u)\, dA\).

  • \(F_{sol}\) — contribution of the solid volume \(V_s\) when the object-centred frame accelerates with \(\mathbf{a}(t)\).

  • \(F_{oe}\) — the open-boundary flux: vorticity advected across the edges of the PIV window carries impulse out of the accounting; this term (moments of \(\omega\, \mathbf{u}\) over the four edges) restores it.

  • \(F_\mu\) — the viscous boundary term: moments of the vorticity gradient \(\mu\, \partial\omega/\partial n\) along the four edges (the diffusive counterpart of \(F_{oe}\); usually small at high Reynolds number).

4.2 Practical notes

  • Volume integrals use Simpson’s rule on the regular grid; boundary terms use one-sided vorticity gradients at the outermost grid lines.

  • The time derivative in \(F_i\) is the noise-critical step. Instead of raw finite differences, the impulse history is locally fitted with a quadratic over a 5-sample window and differentiated analytically (lower-order fallbacks at the first/last samples), filtering high-frequency noise out of the force signal.

  • The field of view should capture the vorticity connected to the force: strong vortices crossing the boundary are exactly what \(F_{oe}\) compensates, but the compensation is only as good as the edge data.

import dpivsoft.Postprocessing as post

X, Y, omega = post.vorticity(x, y, u, v, 'circulation')   # per snapshot
F_v, F_i, F_sol, F_oe, F_mu = post.ImpulseMethod(
    x, y, u, v, omega, rho, Vsol, mu, solid_points, t, accel)
F_total = F_v + F_i + F_sol + F_oe + F_mu

Examples/forces_tutorial.py runs this end to end.

5. Method 2 — the projection method (ProjectionMethod)

5.1 Theory

The projection method (Chang 1992) extracts the force by projecting the Navier–Stokes equation onto a harmonic potential tailored to the body. For each force direction \(i \in \{x, y\}\), define the auxiliary potential \(\phi_i\) by

\[ \nabla^2 \phi_i = 0, \qquad \left.\mathbf{n} \cdot \nabla \phi_i\right|_{\text{body}} = -\,n_i, \qquad \left.\phi_i\right|_{\text{far field}} = 0 \]

— the potential flow generated by unit motion of the body in direction \(i\). Multiplying the momentum equation by \(\nabla\phi_i\) and integrating by parts makes the pressure term vanish (harmonicity + the boundary conditions), leaving

\[ F_i \;=\; \underbrace{\rho \int_V (\boldsymbol{\omega} \times \mathbf{u}) \cdot \nabla\phi_i\, dV}_{F_v\;\text{(vortex contribution)}} \;+\; \underbrace{\mu \oint_{S} \omega\; \big[\mathbf{n} \times \nabla(\phi_i + x_i)\big]\, dS}_{F_\mu\;\text{(viscous surface)}} \;+\; \underbrace{-\sum_j M_{ij}\, a_j}_{F_{am}\;\text{(added mass)}} \]

  • \(F_v\) — a volume integral weighting the Lamb vector \(\boldsymbol{\omega} \times \mathbf{u}\) by the potential gradient; in 2-D components \(\rho \iint \omega\, (u\, \partial_y \phi_i - v\, \partial_x \phi_i)\, dA\). Grid points inside the body are excluded.

  • \(F_\mu\) — a line integral of the wall vorticity along the body surface (Surface_projection(): ordered boundary path

    • normals, Simpson’s rule along arc length).

  • \(F_{am}\) — the added-mass reaction in an accelerating frame, with \(M_{ij} = -\oint_S n_i\, \phi_j\, dS\) computed once from the potentials (compute_added_mass()).

Compared with the impulse method, the weight \(\nabla\phi_i\) decays away from the body, so distant vorticity — and distant measurement noise — contributes little: the integrand is naturally localised. The flip side is that the same weight concentrates the integrand near the wall, exactly where PIV resolution is poorest — the accuracy of the total force is limited by the near-body data, and \(\nabla\phi_i\) develops a very high localised spot at any sharp corner of the body, which a PIV grid cannot resolve. The method is therefore best suited to smooth bodies without corners and to measurements with fine near-wall resolution. There is also the cost that \(\phi_i\) must be computed for the actual body shape, which is what the FEM machinery provides.

5.2 What grad_phi is

ProjectionMethod() takes the two potential gradients stacked as an array of shape (2, 2, ny, nx) on the PIV grid, with grad_phi[i, j] \(= \partial \phi_i / \partial x_j\). It is produced by the FEM pipeline below and interpolated onto the PIV grid with projectionMesh2Grid().

6. The FEM machinery (meshTools.py)

The auxiliary problems of §5.1 are solved once per body geometry:

  1. mesh_generator() builds an unstructured triangular mesh around the body with gmsh: the body outline (circle, spline, or point list) becomes the "object" boundary, the outer rectangle the "outbound" boundary, with element size graded from fine at the wall (tmr) to coarse far away (tm).

  2. FEM_Solver() solves one Laplace problem with scikit-fem: P1 (linear triangle) elements; the body condition enters as a Neumann facet term, the outer boundary is clamped to \(\phi = 0\). Gradients are evaluated at element centres.

  3. projection_FEM_Solver() runs the solver for both unit directions, stacks \(\phi\) and \(\nabla\phi\), computes the added-mass tensor, optionally plots, and saves everything as an .npz.

  4. projectionMesh2Grid interpolates the element-centre gradients onto the (regular) PIV grid; Read_Mesh() reloads a saved .npz.

import dpivsoft.meshTools as mt

mt.mesh_generator(obj, dirSave, ...)                       # gmsh mesh
mesh, cells, elems, phi, grad_phi, added_m = \
    mt.projection_FEM_Solver("mesh.msh", dirSave)          # potentials + added mass
grad_phi_grid = mt.projectionMesh2Grid(elems, grad_phi, X, Y, points)

Ft, Fv, Fmu, Fam = post.ProjectionMethod(
    x, y, u, v, omega, grad_phi_grid, rho, mu, solid_points, added_m, accel)

Examples/mesh_tutorial.py walks through the mesh + FEM steps.

7. Choosing between the two methods

impulse (§4)

projection (§5)

extra inputs

none

FEM mesh + potentials for the body shape

pressure needed

no

no

far-field noise

enters the impulse moment directly

suppressed by the decaying \(\nabla\phi_i\) weight

FOV edges

explicit correction terms (\(F_{oe}\), \(F_\mu\))

integrand already small at the edges

time derivative

required (impulse history)

not required for \(F_v\)

near-body resolution

no special demand (global vorticity moments)

critical — the \(\nabla\phi_i\) weight concentrates the integrand at the wall

body shape

only through wall vorticity

built into \(\phi_i\); sharp corners create a localised \(\nabla\phi_i\) spike the PIV grid cannot resolve

Rule of thumb: for an accurate total force value the impulse method is usually the safer choice — it is mesh-free and places no special resolution demand near the body. The projection method concentrates everything where the measurement is weakest: it needs much finer PIV resolution near the wall, and it prefers smooth surfaces, since a sharp corner turns \(\nabla\phi_i\) into a very high localised spot. Its distinctive strength is interpretation rather than the number itself: the pointwise integrand \(\rho\, (\boldsymbol{\omega} \times \mathbf{u}) \cdot \nabla\phi_i\) is a map of where and how the vorticity contributes to the total force — a visualisation the impulse method cannot provide. Running both on the same data remains a useful consistency check.

8. Function reference

function

role

Postprocessing.vorticity

vorticity field, five stencils (§3)

Postprocessing.walls_vorticity

one-sided wall vorticity along the body path

Postprocessing.Object

body mask + ordered boundary path from a polygon

Postprocessing.ImpulseMethod

vortical impulse force (§4)

Postprocessing.ProjectionMethod

projection-method force (§5)

Postprocessing.Surface_projection

viscous surface integral of the projection method

Postprocessing.ControlVolume

classical control-volume momentum balance (for comparison)

meshTools.mesh_generator

gmsh mesh around the body

meshTools.FEM_Solver / projection_FEM_Solver()

auxiliary potentials \(\phi_i\), \(\nabla\phi_i\)

meshTools.compute_added_mass

added-mass tensor \(M_{ij} = -\oint n_i \phi_j\, dS\)

meshTools.projectionMesh2Grid

FEM gradients → PIV grid

meshTools.Read_Mesh

reload a saved FEM result (.npz)

References

  1. A. Martín-Alcántara, R. Fernández-Feria (2019), “Assessment of two vortex formulations for computing forces of a flapping foil at high Reynolds numbers,” Phys. Rev. Fluids 4, 024702. — The formulation both methods implement (here in dimensional form).

  2. J.-Z. Wu, X.-Y. Lu, L.-X. Zhuang (2007), “Integral force acting on a body due to local flow structures,” J. Fluid Mech. 576, 265–286. — The vortical impulse formulation.

  3. C.-C. Chang (1992), “Potential flow and forces for incompressible viscous flow,” Proc. R. Soc. Lond. A 437, 517–525. — The projection method and its auxiliary potentials.

  4. M. Raffel, C. E. Willert, S. T. Wereley, J. Kompenhans, Particle Image Velocimetry: A Practical Guide, 2nd ed., Springer (2007). — The vorticity stencils of §3.