Abstract
We describe and characterize a method for estimating the pressure field corresponding to velocity field measurements such as those obtained by using particle image velocimetry. The pressure gradient is estimated from a time series of velocity fields for unsteady calculations or from a single velocity field for quasisteady calculations. The corresponding pressure field is determined based on median polling of several integration paths through the pressure gradient field in order to reduce the effect of measurement errors that accumulate along individual integration paths. Integration paths are restricted to the nodes of the measured velocity field, thereby eliminating the need for measurement interpolation during this step and significantly reducing the computational cost of the algorithm relative to previous approaches. The method is validated by using numerically simulated flow past a stationary, twodimensional bluff body and a computational model of a threedimensional, selfpropelled anguilliform swimmer to study the effects of spatial and temporal resolution, domain size, signaltonoise ratio and outofplane effects. Particle image velocimetry measurements of a freely swimming jellyfish medusa and a freely swimming lamprey are analyzed using the method to demonstrate the efficacy of the approach when applied to empirical data.
INTRODUCTION
A longstanding challenge for empirical observations of fluid flow is the inability to directly access the instantaneous pressure field using techniques analogous to those established to measure the velocity field. Recent approaches have made significant progress, especially in the measurement of pressure associated with unsteady fluid–structure interactions (e.g. Hong and Altman, 2008; Jardin et al., 2009a; Jardin et al., 2009b; David et al., 2009; Rival et al., 2010; Rival et al., 2011; David et al., 2012; Tronchin et al., 2012; van Oudheusden, 2013; Liu and Katz, 2013). However, prior efforts have not achieved explicit pressure estimation for moving bodies with timedependent shape, such as those characteristic of animal locomotion and feeding. The pressure field of swimming animals is complicated by the interaction between pressure associated with vortices in the flow and the irrotational pressure field due to acceleration of the body, often referred to as the acceleration reaction or added mass (Daniel, 1984).
Existing methods for empirical pressure estimation often require relatively complex measurement techniques such as multicamera or timestaggered, multiexposure particle image velocimetry (PIV) (Jensen and Pedersen, 2004; Liu and Katz, 2006). In addition, significant computational costs can be associated with the postprocessing required to derive the pressure field from measurements of the velocity or acceleration fields. These postprocessing approaches generally fall into one of two categories. In the first case, the pressure field is computed as a solution to a Poisson equation, e.g. in an inviscid flow: (1) where p is the pressure, u is the velocity vector, ρ is the fluid density, and D/Dt is the material derivative, i.e. the time rate of change of an idealized infinitesimal fluid particle in the flow. Solution of Eqn 1 poses challenges in practice because measurement errors accumulate due to the required temporal and spatial derivatives of u, the condition number (i.e. sensitivity) of the Laplacian operator (Golub and Van Loan, 1996), and measurement uncertainty in the boundary conditions, especially at fluid–solid interfaces (Gurka et al., 2000). For attached flows at high Reynolds numbers, the Neumann boundary condition specifying the pressure gradient at fluid–solid interfaces is given by the boundary layer approximation as ∂p/∂n≈0, where n is the direction of the local normal surface vector (Rosenhead, 1963). However, for separated flows at moderate or low Reynolds numbers, such as those commonly found in animal locomotion, a priori determination of the appropriate fluid–solid boundary conditions for solution of Eqn 1 can be intractable.
A second category of approaches for pressure field estimation is those based on direct integration of the pressure gradient term in the Navier–Stokes equation, e.g. for incompressible flow: (2) where ν is the kinematic viscosity of the fluid. The pressure difference between two points in the domain is determined by integration of Eqn 2 between the two points. For example, the difference in pressure between two points x_{1} and x_{2} is given by: (3) Because measurement errors accumulate along the path of integration from x_{1} to x_{2} in Eqn 3, various techniques have been employed to make this approach less sensitive to measurement uncertainty. A common strategy is to take advantage of the scalar property of the pressure field, such that its local value is independent of integration path. Therefore, each independent integration path that arrives at a point in the flow is in principle an independent estimate of the pressure at that point, provided that measurement errors are uncorrelated. By polling a large number of integration paths, an estimate of the local pressure can be achieved. For example, one successful method (Liu and Katz, 2006) uses an iterative scheme that averages 2m(n+m)+2n(2m+n) integration paths on an m×n grid in order to estimate the instantaneous pressure field.
The aforementioned iterative scheme, while effective in limiting the influence of measurement errors, still incurs a relatively high computational cost. For example, for a 128×128 grid of velocity vectors that is commonly acquired using PIV, the method requires 1.6×10^{5} integration paths per iteration of velocity field integration; and several iterations can be required for convergence of the method (Liu and Katz, 2006). Furthermore, if each integration path is taken as a straight line through the domain, then the method requires interpolation of the estimated pressure gradient field in order to evaluate integration path points that do not coincide with the original data grid. While these requirements are not necessarily prohibitive for twodimensional calculations, they are time consuming and are indeed a showstopper for extension of the method to three dimensions.
We present a simple yet demonstrably effective approach for pressure estimation that is in the spirit of the second category of pressure estimation methods. The method is validated by using two numerically simulated flows: flow past a twodimensional, stationary bluff body and the flow created by a threedimensional, selfpropelled anguilliform swimmer. The first flow is used to characterize a quasisteady implementation of the algorithm, in which the pressure field is estimated from a single velocity field measurement. The second flow demonstrates the ability of the method to accurately estimate the pressure on unsteady, deformable bodies such as those of relevance in animal locomotion. Both flows are used to characterize the method, including its numerical convergence properties and sensitivity to domain size, signaltonoise ratio and outofplane effects. Furthermore, we apply the method to PIV measurements of a freely swimming jellyfish medusa and a freely swimming lamprey, and show that this tool can be applied to the type of measurement data commonly acquired in research.
The reader is strongly encouraged to proceed to the Materials and methods section of the paper before continuing to the Results and Discussion.
RESULTS AND DISCUSSION
Quasisteady pressure estimation
Fig. 1 compares an instantaneous pressure field from the numerical simulation of flow past a stationary bluff body with the pressure field estimated from the corresponding velocity field using the quasisteady form of the present algorithm. A vector field spatial resolution of D/16 (where D is the side length of the bluff body) is used in the horizontal and vertical directions to mimic typical PIV measurements. The salient features of the flow, especially the high pressure on the upstream face of the bluff body and the low pressure in the shear layers and nearwake vortices, are well captured by the algorithm (see supplementary material Appendix S1 for discussion of discrepancies in the far wake). Furthermore, the maximum and minimum pressures in the field are in quantitative agreement (see supplementary material Fig. S1). To be sure, nearestneighbor Gaussian smoothing creates a spurious thin layer of undefined pressure at the fluid–solid interface and moves the pressure peak on the upstream face of the body away from the interface. However, the correct nearbody pressure can be recovered by increasing the grid resolution so that the nearestneighbor filter artifact on the body surface is limited to a smaller region very close to the body. Additional surface pressure calculations for the quasisteady case (see supplementary material Appendix S1) are based on a velocity vector spacing of D/64. Note that a similar increase in resolution using a PIV camera would require a concomitant reduction in the measurement window size due to limits on camera pixel resolution.
Additional characterization of the quasisteady algorithm is detailed in supplementary material Appendix S1, including: analysis of spatial convergence; the relative contribution of each integration path to the median pressure field; robustness to measurement noise; and the effects of domain size, fluid viscosity and fluid–solid interfaces.
Unsteady pressure estimation
Fig. 2 compares an instantaneous pressure field from the numerical simulation of a selfpropelled anguilliform swimmer with the pressure field estimated from the corresponding velocity field using the unsteady form of the present algorithm. A vector field spatial resolution of L/42 (where L is the length of the swimmer) is used in the horizontal and vertical directions. No smoothing is applied to this data set in order to contrast the results with those in the previous section and to limit the spatial extent of the region of undefined pressure near the fluid–solid interface. The algorithm is effective in capturing: the high–low pressure couples formed on the sides of the swimmer head and tail as they accelerate in the positivey direction; the low–high pressure couple formed at the midbody as it accelerates in the negativey direction; and the pressure in the wake vortices.
The importance of the unsteady term in Eqn 9 (see Materials and methods) is illustrated by comparison with the pressure field estimated using the quasisteady approximation, shown in Fig. 2C. Low pressure in the wake vortices is captured, but the high–low pressure couples on the body surface due to the body added mass are missing entirely, as is the high pressure in the wake due to vortex added mass (Dabiri, 2006). The comparison is further quantified in Fig. 3, which plots the pressure on a contour surrounding the swimmer and immediately adjacent to the region of undefined pressure. At each of the four phases of the swimming cycle shown, good agreement is achieved between the pressure computed in the numerical simulation and the pressure estimated from the velocity field using the unsteady algorithm. By contrast, the pressure estimated by the quasisteady algorithm is erroneous everywhere except near the forming wake vortex at the tail.
Additional characterization of the unsteady algorithm is provided in supplementary material Appendix S2, including analysis of temporal convergence and outofplane effects for threedimensional flows.
To demonstrate the efficacy of the present method for analyzing empirically measured velocity fields, Fig. 4 shows measured velocity and vorticity fields for the freely swimming jellyfish and lamprey (Fig. 4A,C) along with the corresponding pressure fields estimated using the unsteady algorithm (Fig. 4B,D). The full measurement domain is shown in both cases; the velocity vector field is plotted at half of the full resolution. Only the left half of the jellyfish body is visible in the measurement domain; its exumbrellar surface is indicated by a black curve in Fig. 4A,B. The full lamprey body is visible in Fig. 4C,D.
In both cases, the pressure field derived from the velocity field measurements captures key features near the body surface and in the wake. In particular, the jellyfish data set indicates low pressure in the forming starting vortex and high pressure where the bell margin is accelerating inward and pushing the adjacent fluid. The results are consistent with the measured vorticity field (Fig. 4C), with the region of low pressure corresponding to the core of the starting vortex. The presence of low and high pressure regions near the bell margin is also in agreement with previous numerical simulations of a swimming jellyfish with similar body shape and kinematics (Sahin et al., 2009).
The lamprey data set shares similarities with the threedimensional numerical model shown previously. The vorticity and pressure fields are less smooth and show finer structure in the empirical measurements, which is attributable in part to the Reynolds number being approximately four times higher than that of the numerical simulation.
The ease of implementation of this algorithm, in terms of both data acquisition and velocity field postprocessing, and its relatively low computation cost (see supplementary material Appendix S2) gives it the potential to find use in a broad range of problems of interest in biological fluid mechanics. Because the temporal filter implemented in the unsteady algorithm does add considerable time to the pressure calculation (cf. supplementary material Fig. S10), in practice one should first evaluate the results of both the quasisteady and the fully unsteady implementations of the algorithm on a sample of the data of interest to determine whether unsteady effects are important. If they are not, then the quasisteady calculation provides the most efficient tool for determination of the pressure field.
Although the present evaluation focused on twodimensional velocity fields, it is straightforward to extend the algorithm to three dimensions by the addition of a limited number of new integration paths consistent with the geometry in Fig. 5. In that case, even greater reductions in computation expense can be achieved relative to existing methods as a result of the relatively small total number of required integration paths and the elimination of associated velocity field interpolation during integration of the pressure gradients.
A free MATLAB implementation of this algorithm is available at: http://dabiri.caltech.edu/software.html.
MATERIALS AND METHODS
Material acceleration estimation
The instantaneous fluid particle acceleration Du/Dt required for calculation of the pressure gradient in Eqn 2 is estimated by advecting idealized infinitesimal fluid particles in the measured velocity fields. For quasisteady estimation, the material acceleration is derived from a single velocity field as: (4) where i=1, 2… m×n (i.e. the dimensions of the velocity grid), x_{i} are the positions of fluid particles coincident with the grid points in the PIV velocity field, and are the positions of those fluid particles after being advected by the instantaneous velocity field for a period Δt: (5)
In order for Eqns 4 and 5 to remain valid, Δt is limited to values much smaller than the characteristic time scale of the flow, yet sufficiently large that there is a measurable change in the fluid particle velocity.
For many flows, especially those involving accelerating or deforming bodies, the aforementioned constraint on Δt cannot be satisfied. For these inherently unsteady fluid–structure interactions, we derive the material acceleration from two sequential velocity fields as: (6) where (7) Eqn 7 is akin to a Crank–Nicolson (i.e. trapezoidal) scheme for the particle positions, in contrast to the forward Euler scheme in Eqn 5. Hence, the convergence of the method with time step is second order (Crank and Nicolson, 1947).
The primary source of measurement error in this type of unsteady estimate of the material acceleration Du/Dt arises from temporal noise in the measured velocity components at each node in the velocity field. We address this by applying a temporal filter to the time series of velocity fields, which results in a smoothing spline approximation u* to the velocity u at each node in the velocity field. The spline approximations are defined such that they minimize, for each component of u, the parameter: (8) where t=1…N is the temporal sequence of velocity fields to be analyzed, u_{τ} is a velocity vector corresponding to velocity field t in the sequence, u_{τ}* is the splineapproximated value of the same velocity vector for the same velocity field in the sequence, t_{min} and t_{max} are the temporal bounds on the sequence of velocity fields, and ϕ is a weight between the first and second terms and has a value between 0 and 1. In effect, the parameter S_{u} quantifies both the deviation of the spline approximation from the original data (i.e. the first term) and the total curvature magnitude of the spline approximation (i.e. the second term). For ϕ=0, only the second term is minimized, resulting in a leastsquares fit with zero curvature, i.e. a linear fit to the data. For ϕ=1, only the first term is minimized, giving a cubic spline fit that passes through each original data point. In all that follows, we set ϕ=0.05, a value we have identified as enabling effective temporal noise filtering without discarding true temporal trends in the measurement data.
Further characterization of the temporal filter is provided in supplementary material Appendix S2. In particular, it is shown that the use of the temporal filter increases the order of temporal convergence above second order, as anticipated by theory (Atkinson, 1968).
It is worth noting that the distinction between the quasisteady and unsteady approaches can be made explicit by decomposing the material acceleration into its Eulerian components: (9) The quasisteady approximation in Eqns 4 and 5 implicitly neglects the first term on the righthand side of Eqn 9, whereas the unsteady calculation retains it.
The viscous term on the righthand side of Eqn 2 is computed using centered finite differences between adjacent nodes in the velocity field. The effect of the viscous term is evaluated in the context of a numerical simulation described in the validation section.
Pressure gradient integration
Whereas previous methods that integrate the pressure gradient via many integration paths assign to each grid point the arithmetic mean of the many integrations, in the present approach the paths are polled by taking the median. The median is less sensitive to grossly erroneous values that may arise on a few of the integration paths due to localized measurement errors or to localized errors created by the aforementioned material acceleration approximations in Eqns 4, 5, 6 and 7. Hence, this approach enables a significant reduction in the total number of integration paths per frame that are required to achieve accurate pressure estimates. Fig. 5 illustrates the paths used presently. Eight families of integration paths are used, with each family originating at the domain boundary and propagating toward each grid point from the left (L), upper left (UL), top (T), upper right (UR), right (R), lower right (LR), bottom (B) and lower left (LL), respectively.
Only eight integration paths (one per family) per grid point are used, for a total of 8m×n paths per velocity field. For the aforementioned example grid of 128×128 velocity vectors, 1.3×10^{5} integration paths are required, a 20% reduction from existing optimal methods (Liu and Katz, 2006). More importantly, the integration paths are constrained to include only grid points coincident with the original velocity field. For example, the UL integration path is composed of alternating integration steps in the −y and +x directions, originating at the domain boundary and terminating at each grid point. Hence, no interpolation is required in order to integrate the pressure gradient field. Furthermore, portions of many of the paths are redundant, facilitating fast calculation using simple matrix manipulations. A forward Euler spatial integration scheme is used throughout, resulting in firstorder spatial convergence of the method (see supplementary material Appendix S1).
An important limitation of the present algorithm that arises from the tradeoff between speed and accuracy is that it assumes the pressure is zero at the point on the outer domain boundary where each integration path is initiated. This does not imply, however, that the final pressure estimate is constrained to be zero at the boundaries. Integration paths that originate from the other domain boundaries and terminate at a given boundary may estimate a nonzero value of pressure at the termination point. If the median of all paths terminating at that point on the domain boundary is nonzero, then the final pressure estimate at that point will also be nonzero. Note that for all points in the domain, the final pressure estimate is relative to a zero reference pressure, as that is the pressure at the origin of each integration path. The impact of these assumptions on the robustness of the technique is quantified in supplementary material Appendix S1, and it is shown to be modest for the external flows tested. At the same time, the net result of this tradeoff in the algorithm design is a more than orderofmagnitude reduction in computational time compared with previous methods (see supplementary material Appendix S2).
A common source of localized error that can affect pressure estimates is the presence of solid objects in the flow. Typical PIV measurements are often unreliable in the region close to solid objects, which compromises pressure integration paths that cross the fluid–solid interface, especially in previous methods that average the erroneous data instead of discarding it via median polling (or in Poisson solvers that rely on the pressure gradient at the fluid–solid interface as a boundary condition). In the present algorithm, integration paths that cross a fluid–solid interface in the flow can be nullified by assigning the nodes nearest to the interface an undefined pressure gradient. Hence, when that value is integrated along any integration path, the pressure value for that path also becomes undefined and therefore does not contribute to the median calculation.
Validation data sets
To validate the accuracy of the quasisteady pressure estimates achieved using this algorithm, a numerical simulation of flow past a twodimensional square crosssection cylinder at a Reynolds number of 100 was used. This numerical data set enabled quantification of the effects of spatial resolution, domain size and signaltonoise ratio, while providing a known pressure field standard for comparison (see supplementary material Appendix S1). The numerical simulation was executed using a solver that computes on arbitrary polyhedra (Ham and Iaccarino, 2004). In the present case, a regular Cartesian mesh was utilized and subsequently interpolated onto coarser grids of varying sizes typical of PIV data. The viscous term in Eqn 2 was retained in all of the calculations to demonstrate the robustness of the median polling approach to errors normally associated with application of the Laplacian operator. For all calculations of Eqns 4 and 5 in this validation, we set Δt=0.01h/U_{max}, where h is the mean grid spacing and U_{max} is the maximum flow speed in the measurement domain. The results described below were insensitive to orderofmagnitude larger and smaller values of Δt. Where noted, nearestneighbor Gaussian smoothing was applied both to the pressure gradient before integration and to the resulting pressure field.
The accuracy of the fully unsteady pressure estimates was validated by using a published numerical simulation of a threedimensional, selfpropelled anguilliform swimmer (Kern and Koumoutsakos, 2006). The Reynolds number based on swimmer length and speed was ~2400. Time steps between sequential velocity fields from 0.02T to 0.08T (where T is the swimming stroke duration) were studied to quantify the temporal convergence of the method. The validation results described in the Results section are based on calculations of Eqns 6 and 7 using velocity fields separated by 0.02T.
Empirical data sets
The present method was also applied to PIV measurements of a freely swimming Aurelia aurita Linnaeus 1758 jellyfish medusa and a freely swimming Anguilla rostrata Lesueur 1817 lamprey to demonstrate the performance of the algorithm with empirical data inputs and, in the case of the jellyfish, without treatment of fluid–solid interfaces. The swimming Reynolds numbers of the jellyfish and lamprey were ~1000 and 10,000, respectively, and the time between sequential velocity fields was 5 ms (t/T=0.013) and 4 ms (t/T=0.015), respectively. Details of the PIV implementation can be found in published literature (Colin et al., 2012).
ACKNOWLEDGEMENTS
The authors gratefully acknowledge S. Kern and P. Koumoutsakos for providing the threedimensional numerical simulation of the selfpropelled swimmer, and X. Liu for providing access to the pressure gradient integration code used in Liu and Katz (Liu and Katz, 2006).
FOOTNOTES

Author contributions
J.O.D. conceived of the algorithm; J.O.D., S.B., B.J.M., S.P. and S.C. designed and executed the experiments; J.O.D. drafted the paper.

Competing interests
The authors declare no competing financial interests.

Funding
This research was supported by Office of Naval Research awards N000140810918 and N000141010137 to J.O.D. Deposited in PMC for immediate release.

Supplementary material
Supplementary material available online at http://jeb.biologists.org/lookup/suppl/doi:10.1242/jeb.092767//DC1
 © 2014. Published by The Company of Biologists Ltd
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.