# 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](#forces-sec-4)) — needs only the velocity/vorticity field; - the **projection method** ([§5](#forces-sec-5)) — needs, in addition, an auxiliary potential solved once by FEM on a mesh around the body ([§6](#forces-sec-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)$. (forces-sec-3)= ## 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, {func}`~dpivsoft.Postprocessing.Object` builds the body mask and an ordered boundary path, and {func}`~dpivsoft.Postprocessing.walls_vorticity` replaces the vorticity along the wall with one-sided estimates. Both force methods apply this automatically when `solid_points` is provided. --- (forces-sec-4)= ## 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. ```python 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. --- (forces-sec-5)= ## 5. Method 2 — the projection method (`ProjectionMethod`) (forces-sec-5-1)= ### 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 ({func}`~dpivsoft.Postprocessing.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 ({func}`~dpivsoft.meshTools.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 {func}`~dpivsoft.Postprocessing.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 {func}`~dpivsoft.meshTools.projectionMesh2Grid`. --- (forces-sec-6)= ## 6. The FEM machinery (`meshTools.py`) The auxiliary problems of [§5.1](#forces-sec-5-1) are solved once per body geometry: 1. **{func}`~dpivsoft.meshTools.mesh_generator`** builds an unstructured triangular mesh around the body with [gmsh](https://gmsh.info): 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. **{func}`~dpivsoft.meshTools.FEM_Solver`** solves one Laplace problem with [scikit-fem](https://github.com/kinnala/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. **{func}`~dpivsoft.meshTools.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; {func}`~dpivsoft.meshTools.Read_Mesh` reloads a saved `.npz`. ```python 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](#forces-sec-4)) | projection ([§5](#forces-sec-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 | |----------|------| | {func}`Postprocessing.vorticity ` | vorticity field, five stencils ([§3](#forces-sec-3)) | | {func}`Postprocessing.walls_vorticity ` | one-sided wall vorticity along the body path | | {func}`Postprocessing.Object ` | body mask + ordered boundary path from a polygon | | {func}`Postprocessing.ImpulseMethod ` | vortical impulse force ([§4](#forces-sec-4)) | | {func}`Postprocessing.ProjectionMethod ` | projection-method force ([§5](#forces-sec-5)) | | {func}`Postprocessing.Surface_projection ` | viscous surface integral of the projection method | | {func}`Postprocessing.ControlVolume ` | classical control-volume momentum balance (for comparison) | | {func}`meshTools.mesh_generator ` | gmsh mesh around the body | | {func}`meshTools.FEM_Solver ` / {func}`~dpivsoft.meshTools.projection_FEM_Solver` | auxiliary potentials $\phi_i$, $\nabla\phi_i$ | | {func}`meshTools.compute_added_mass ` | added-mass tensor $M_{ij} = -\oint n_i \phi_j\, dS$ | | {func}`meshTools.projectionMesh2Grid ` | FEM gradients → PIV grid | | {func}`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](https://doi.org/10.1103/PhysRevFluids.4.024702)," *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](https://doi.org/10.1017/S0022112006004551)," *J. Fluid Mech.* **576**, 265–286. — The vortical impulse formulation. 3. C.-C. Chang (1992), "[Potential flow and forces for incompressible viscous flow](https://doi.org/10.1098/rspa.1992.0077)," *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](https://doi.org/10.1007/978-3-540-72308-0)*, 2nd ed., Springer (2007). — The vorticity stencils of [§3](#forces-sec-3).