An algorithm to estimate unsteady and quasi-steady pressure fields from velocity field measurements

A long-standing 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 time-dependent 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).

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⁵ 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 two-dimensional 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 two-dimensional, stationary bluff body and the flow created by a three-dimensional, self-propelled anguilliform swimmer. The first flow is used to characterize a quasi-steady 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, signal-to-noise ratio and out-of-plane 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.

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 quasi-steady 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 near-wake 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, nearest-neighbor 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 near-body pressure can be recovered by increasing the grid resolution so that the nearest-neighbor filter artifact on the body surface is limited to a smaller region very close to the body. Additional surface pressure calculations for the quasi-steady 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 quasi-steady 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.

Fig. 2 compares an instantaneous pressure field from the numerical simulation of a self-propelled 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 positive-y direction; the low–high pressure couple formed at the mid-body as it accelerates in the negative-y 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 quasi-steady 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 quasi-steady 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 out-of-plane effects for three-dimensional 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 three-dimensional 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 post-processing, 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 quasi-steady 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 quasi-steady calculation provides the most efficient tool for determination of the pressure field.

Although the present evaluation focused on two-dimensional 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.

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.

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).

The viscous term on the right-hand 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.

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⁵ 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 first-order spatial convergence of the method (see supplementary material Appendix S1).

An important limitation of the present algorithm that arises from the trade-off 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 non-zero value of pressure at the termination point. If the median of all paths terminating at that point on the domain boundary is non-zero, then the final pressure estimate at that point will also be non-zero. 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 trade-off in the algorithm design is a more than order-of-magnitude 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.

To validate the accuracy of the quasi-steady pressure estimates achieved using this algorithm, a numerical simulation of flow past a two-dimensional square cross-section cylinder at a Reynolds number of 100 was used. This numerical data set enabled quantification of the effects of spatial resolution, domain size and signal-to-noise 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 order-of-magnitude larger and smaller values of Δt. Where noted, nearest-neighbor 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 three-dimensional, self-propelled 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.

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).

The authors gratefully acknowledge S. Kern and P. Koumoutsakos for providing the three-dimensional numerical simulation of the self-propelled swimmer, and X. Liu for providing access to the pressure gradient integration code used in Liu and Katz (Liu and Katz, 2006).