Flow pressure evaluation on generic surfaces by robotic volumetric PTV

An experimental approach for the measurement of the time-average fluid flow pressure over the surface of generic three-dimensional objects is presented. The method is based on robotic volumetric PTV measurements followed by the integration of the pressure gradient. The domain for pressure evaluation is subdivided in two parts: in the irrotational region the static pressure is obtained following Bernoulli relation; in the turbulent wake and close to the object the pressure gradient is integrated. An approach based on the total pressure distribution is proposed to estimate the boundary between these two regions. The method is first assessed with experiments around a sphere equipped with pressure taps. A criterion for minimum spatial resolution is formulated in terms of maximum ratio between bin size and local radius of curvature of the object. An experimental database from a three-dimensional problem of higher geometrical complexity is considered: the time-averaged flow field around a full-scale cyclist. The surface pressure distribution is discussed in connection to the topological features of near-surface streamlines and streamwise vortices.


Introduction
The measurement of surface pressure in aerodynamics is of paramount importance for many engineering applications, from ground transport to aviation and including wind energy. The interest stems from the fact that pressure is, next to skin friction, the only mechanism for aerodynamic forces being exerted to an object immersed in the fluid stream. For wind tunnel experiments, surface pressure taps are considered as the standard approach to instrument scaled models, providing accurate local values of the surface pressure. Installation of Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. pressure taps impacts the design and manufacturing of scaled models with multitudes of pressure orifices and highways of small tubes towards the transducer unit. In some cases the installation of pressure taps becomes challenging, in proximity of sharp edges, flexible objects or when studying the aerodynamics around the human body. Pressure sensitive paint (PSP) (Mclachlan and Bell 1995) offers the ability to map the distributed pressure field on the surface by non-invasive measurements. Most applications of PSP are reported in the compressible flow regime (Mclachlan and Bell 1995, Klein et al 2005, Gregory et al 2008, Bitter et al 2012, whereas at low speed the small pressure variations exhibited at the surface are not accurately described with PSP (Tagliabue et al 2017).
The principle of inferring the flow field pressure from velocity measurements (mostly PIV) has already been reported two decades ago by Gurka et al (1999), followed by a detailed review from van Oudheusden (2013). Many studies utilize planar PIV data both in the incompressible and compressible flow regimes (Gurka et al 1999, van Oudheusden et al 2006, Ragni et al 2009, 2011, Auteri et al 2015. In contrast with the advancement of PIV systems for measurements at industrialscale (Michaux et al 2018) and volumetric PIV techniques, the recent review work of Discetti and Coletti (2018), only reports few works capturing the full 3D3C velocity-and pressure field and for measurements often confined to relatively small volumes (Violato et al 2011, Pröbsting et al 2013. The combination of tracer particle technology for large-scale wind tunnel experiments (helium filled soap bubbles (HFSB) Scarano et al 2015) and volumetric velocimetry by coaxial arrangement  along with the advancements in 3D particle motion analysis (Shake-the-Box, Schanz et al 2016) has shown the potential to evaluate the velocity and pressure field in volumes of several litres. Schneiders et al (2016) reported the analysis of the pressure field downstream of a truncated cylinder and Huhn et al (2018) in the impinging jet, respectively.
The above works have dealt with the measurement over flat surfaces for the 3D studies, while the planar applications yielded measurements along a line over cylindrical objects and aerofoils.
PIV-based pressure measurements over generic threedimensional surfaces are challenged by the difficulty to simultaneously illuminate a complex curved surface and record images of the tracer particles travelling in the flow field above it. Coaxial volumetric velocimetry (CVV, Schneiders et al 2018) and its application in a robotic measurement approach (originally published as 'robotic volumetric PIV  herein referred to as 'robotic volumetric PTV') have been demonstrated for the analysis of the flow field and the near-surface flow pattern of a full-scale cyclist. Evaluating the surface pressure field in such conditions would be beneficial for the characterisation of the aerodynamic loads acting on the model and their optimisation (e.g. drag reduction). The insight into the local pressure distribution can thereby complement the analysis of 'global' loads as resulting from the wake survey approach (Dabiri 2005, Terra et al 2019 which, unlike the time-averaged robotic volumetric PTV measurements, can also provide insight into unsteady systems such as flying birds or spinning rotors (Ragni et al 2011, Stalnov et al 2015. An open question in the evaluation of the time-averaged surface pressure on generic geometries by means of robotic volumetric PTV is the treatment of solid boundaries. The work of Auteri et al (2015), approached a 2D problem uncoupling the velocity and pressure field in the Navier-Stokes equations and retrieving the pressure distribution up to an integration constant without the need of specifying Dirichlet boundary conditions. Lastly, after the seminal work of Liu and Katz (2006), and further developments (rotating parallel ray method, Liu et al 2016) a path integration strategy was demonstrated to mitigate effects of error accumulation by optimizing path selection. Most current studies indicate an agreement towards a global formulation of pressure integration solving the Poisson problem (van Oudheusden 2013, van Gent et al 2017. The effect of error propagation originating from spurious or under-resolved data at the surface, however, has not been fully understood. Moreover, it is recognized that the pressure gradient integration is unnecessary within regions where the flow regime is irrotational (e.g. free stream or outer stream), where simpler models (Bernoulli or isentropic flow relations) relating velocity and pressure can suffice. This principle has been demonstrated in the work from Ragni et al (2009), around a 2D airfoil in the transonic flow regime and in the study of Kurtulus et al (2007), around a square cylinder in the incompressible regime.
The present work describes an experimental approach to evaluate the time-averaged static surface pressure distribution on generic geometries by means of robotic volumetric PTV measurements. Two key elements are addressed in this study: the first is the method to partition the domain for pressure evaluation into a region where Bernoulli's equation is applied and a second one where spatial integration of the pressure gradient is necessary. This part of the work extends the study of Jeon et al (2018), who partition the PIV measurement domain for a sequential pressure integration based on the local measurement reliability, to generic 3D objects, yet for time-averaged analysis. Secondly, the problem of pressure discretisation on a surface that is not compliant with the Cartesian dataset of the velocity measurements is addressed.
The article opens with a brief summary of the working principles of robotic volumetric PTV. Section 2.2 deals with the evaluation of the 3D pressure field. The process of domain partitioning is addressed in detail in section 2.3, followed by the detailed pressure evaluation procedures in the established domains and on the model surface in section 2.4. Experiments are presented first in form of proof-of-principle with the measurement around a sphere. Subsequently, the database pertaining the flow around a full-scale cyclist is examined. Experimental uncertainties are quantified using conventional surface pressure taps over the sphere model.

Methodology and working principles
The pressure field and surface pressure evaluation process by robotic volumetric PTV may be split into two key components: the first step delivers the time-averaged, three-dimensional (3D), three-component (3C) velocity field and its statistical fluctuations, in terms of root-mean-square fluctuations and Reynolds stress components. In this respect, data acquisition and velocity field evaluation are summarized in section 2.1. The second step pertains the post-processing of the velocity data to yield the static pressure distribution in the measurement domain and on the solid object surface. The foundation for the pressure evaluation is the computation of the pressure gradient as discussed in section 2.2. The pressure field can subsequently be obtained by spatial integration of its gradient.
To avoid the gradient integration in regions where the flow regime is irrotational, the measurement domain is partitioned into a region where Bernoulli's equation (Anderson 2011) can provide the static pressure, and a domain where, instead, the pressure gradient needs to be integrated. Such domain partitioning was already proposed by Kurtulus et al (2007). The identification of the two respective regions, which are referred to as irrotational and rotational hereafter, is discussed in section 2.3. The above approach is intended to significantly reduce the length of the paths followed for the pressure gradient integration, in turn reducing the associated error accumulation. Following the discussion on domain partitioning, the pressure evaluation in the separate regions and on the object surface are discussed in section 2.4.

Velocity field measurements by robotic volumetric PTV
Robotic volumetric PTV builds upon the compactness of the CVV measurement probe manipulated by a robotic arm. The working principles (imaging and illumination) are akin to tomographic PIV (Elsinga et al 2006), whereas the data analysis is based on 3D particle tracking velocimetry ('Shakethe-Box' (STB), Schanz et al 2016). The compact hardware arrangement of imagers and the illumination source make robotic volumetric PTV suited for the rapid acquisition of the flow problem from several views.
The low tomographic aperture of the CVV arrangement yields inferior depth accuracy along the viewing direction, as discussed in Schneiders et al (2018). The latter is compensated by the increased velocity dynamic range of the STB analysis. Moreover, the seeding density in physical space (concentration, particles/cm 3 ) as well as in the image plane (N p expresses the number of particles per pixel or ppp) is more constrained, compared to large aperture tomographic systems (Lynch and Scarano 2014 reported experiments up to N p = 0.2 ppp). Early works employing CVV have been conducted with image particle density N p not exceeding 0.02 ppp , while only recently Gallar et al (2020), report CVV recordings at up to 0.09 ppp. As a result, the instantaneous concentration of tracer particles remains below 2 particle/cm 3 in a typical measurement with CVV.
Working principles of robotic volumetric PTV schematically illustrated in figure 1 are outlined in the following. The reader is referred to the works of Schneiders et al (2018), and Jux et al (2018), for a detailed discussions of the measurement system and calibration procedures.
Prior to data acquisition, the CVV probe and the robotic arm are calibrated, such that all measurement results can be presented in the wind-tunnel reference system. Note that calibration is only necessary once, and it can be maintained if the systems are not disassembled.
Position and orientation of the illuminating and 3Dimaging probe are controlled by robotic arm manipulation to allow a scanning of the domain of interest. Particle images are acquired in time-resolved mode, and particle tracking analysis is performed based on the Shake-the-Box algorithm (Schanz et al 2016). The scattered ensemble of particle trajectories is reduced to a Cartesian grid by averaging tracks belonging to a common cell (or bin) of linear size l v , as detailed hereafter. The result is a 3D-3C time-averaged velocity field. The Cartesian grid is latter referred to as the fluid mesh.
The spatio-temporal average of the scattered particle data evaluated in each bin is performed likewise to the method  proposed by Agüera et al (2016), whereby a linear polynomial function is used to fit the velocity distribution inside the bin, following the expression: The fit of a first-order distribution function carries two advantages over a simple arithmetic mean: firstly, velocity fluctuations (u') can be evaluated with respect to the local mean, rather than the cell-centered mean, resulting in an improved estimation of the fluctuations in the presence of a spatial velocity gradient inside a cell. Secondly, the polynomial coefficients a provide directly the local velocity gradient, which is required for the computation of the pressure gradient. Details of the method are discussed in the work of (Agüera et al 2016), including the effect of averaging schemes on the turbulent statistics.
The 3D-3C velocity field constitutes the input for the pressure evaluation procedure discussed in the subsequent sections. As such, the proposed method is not exclusively applicable for robotic volumetric PTV data only, but any other technique providing equivalent velocity field data.

Pressure field evaluation
The time-averaged pressure distribution is related to the incompressible, time-averaged 3D velocity field by the momentum equation. The evaluation of the time-average pressure gradient ∇p is based on the Reynolds decomposition of the velocity according to the following expression (van Oudheusden 2013): The viscous term is not included as it becomes negligible in the high Reynolds number regime (Re > 10 3 ) as discussed in previous studies (van Oudheusden et al 2007, Ghaemi et al 2012. The ensemble averaging of the particle track data provides u ′ , its spatial gradient, and u ′ u ′ , whereas the gradient of u ′ u ′ is estimated by central finitedifferences.

Boundary conditions.
Following the assessment of pressure evaluation methods for 2D PIV (Charonko 2010), the two most common approaches are spatial integration of the pressure gradient (Baur andKöngeter 2012, Liu andKatz 2006) and the solution of the Poisson equation for pressure resulting from the divergence of the momentum equation (de Kat and van Oudheusden 2012). The former approach is prone to accumulation of local measurement errors along the integration path and such errors can vary along different paths to the same point. On the other hand, the method can be tailored in a way that the regions of largest error (often the solid boundaries) are reached only at the end of the path. The Poisson approach provides a global solution and solves the previous problem of sensitivity to the integration path. The method, however, propagates errors arising from the boundaries in the entire domain, let them be Neumann (dp/dn, the pressure gradient in the direction of the surface normal n is specified) or Dirichlet (p is specified) condition. Figure 2 indicates the type of boundary conditions to be specified for the two general approaches.
An additional difficulty arises from the mismatch at the fluid-solid interface between the two systems of coordinates (fluid & surface mesh), which results into uncertainties for specifying the boundary conditions for the Poisson problem. This is illustrated in figure 3 for the case of a generic curved fluidsolid interface.
The illustration shows the effect of a coarse discretization introducing an uncertainty in the boundary condition specification and compromising the pressure field evaluation. Some effects of error propagation through incorrect boundary condition specifications in PIV-based pressure field calculations are discussed in the work of Pan et al (2016).

Dual-model for pressure evaluation.
The propagation of errors from the solid object boundaries is avoided when the values at the boundary do not take part in the pressure integration process. This can be realized by a unidirectional integration procedure (see e.g. Dabiri et al 2014). In this way, the pressure gradient can be integrated from the outer (fluiddomain) boundary towards the object of interest, as illustrated in figure 2 (right). The spatial integration of the pressure gradient over a long path, however is prone to the above discussed error accumulation.
To attenuate such effects, a dual-model approach for the evaluation of the static pressure field is envisaged: the measurement domain is partitioned into a region of irrotational  flow where the Bernoulli equation returns the pressure from the local velocity, and the remaining region (rotational flow) where the pressure gradient is integrated uni-directionally. The equations used in this approach read as The proposed concept is illustrated schematically in figure 4 and it is discussed in detail hereafter, starting with the criterion that identifies the irrotational flow domain.

Criteria to determine the irrotational domain
The goal of domain partitioning is the robust identification of the sub-region where Bernoulli's equation can be applied accurately. The problem translates in determining the interface between the rotational-and irrotational-domain.
Recalling Anderson (2011), a vector field u (x) is defined irrotational if ∇ × u = 0. For the latter, a (scalar) potential φ can be defined such that u = ∇φ. In the assumption of incompressible, inviscid and adiabatic flow, the total pressure remains constant, yielding Bernoulli's principle. The main question arising here is: how can the part of the flow domain with constant total pressure be identified?
The problem has been addressed in a previous publication of the authors, where the region of constant total pressure was approximated by geometry-and turbulence-based surrogates (Jux et al 2019). In the current work the principle of domain partitioning and subsequent pressure evaluation is examined on a dataset obtained for the flow around a sphere (see figure  5-left), which is presented in more detail in sections 3 and 4.
2.3.1. Approximated total pressure. A first approximation of the total pressure field can be obtained when the total pressure at the inflow boundary can be assumed as constant. The resulting distribution of total pressure coefficient, herein defined as  is computed, based on the measured, time-average velocity magnitude |u|, and the static pressure field (p − p ∞ ). The latter is obtained by spatial integration of the pressure gradient, following the pressure integration scheme detailed in section 2.4.2, starting from a constant total pressure boundary condition upstream of the sphere model. The resulting pressure distribution is illustrated in figure 5 (right), where C p,tot drops below unity in the wake region.

Definition of irrotational-rotational interface.
The firstlevel estimation of total pressure can serve to define the edge of the rotational region as to where C P,tot < 1. However, given the finite uncertainty in the determination of C P,tot a threshold needs to be selected for C P,tot below unity, that encompasses such uncertainty. A possible approach to determine the experiment specific threshold level is based on the analysis of the standard deviation of the total pressure coefficient in the irrotational domain σ CP,tot , for instance over a region upstream of the model where the condition of constant total pressure can be assumed valid with high confidence. For the present case figure 6 illustrates such choice, along with the resulting C p,tot contour when allowing for a reduction of 3 × σ CP,tot below unity.
For the case shown here, a threshold of ϵ CP,tot = 98.5% is determined based on a standard deviation of σ CP,tot = 0.5%. It should be retained in mind that defining the domain interface at any value of C P,tot below unity introduces a systematic underestimation of the total pressure at such interface (bias error of 1.5% in the present case).
Such bias error, can be mitigated or eliminated if the rotational domain is dilated, moving the interface towards the irrotational domain, with values of C P,tot closer to unity. The need and extent of such operation may be case-dependent and a detailed discussion is not given here for sake of conciseness. Nevertheless, a comparison of the surface pressure obtained with and without domain partitioning is presented in the results (section 4).

Pressure evaluation in the partitioned domain and on the model surface
The static pressure distribution in the measurement domain and on the object surface is evaluated making use of two models as outlined before.

Pressure evaluation by Bernoulli's equation.
In the steady irrotational fluid domain, the relation between velocity and pressure can be accurately approximated with Bernoulli's equation (i.e. uniform total pressure). The pressure relative to known free-stream conditions (u ∞ , p ∞ ) is then given bȳ where u is the experimentally measured time-averaged velocity. The direct relation between velocity and pressure yields a robust point-wise estimate of the static pressure as it does not depend on boundary conditions and local errors do not accumulate along a path. Figure 7 illustrates the share of the measurement domain that can be evaluated by this approach for the flow around a sphere.

Pressure integration.
For the rotational flow region the static pressure is obtained through spatial integration of the pressure gradient (∇p, see equation (2)), with the pressure at the edge of the irrotational region providing the necessary Dirichlet boundary condition.
The pressure gradient is integrated with a spatial marching process of the whole interface, intended to enforce consistency among adjacent paths. The boundary of the rotational part of the flow field is eroded, starting from an initial interface with the irrotational domain (figure 7). The pressure is assigned for points in the rotational domain of the fluid mesh with a minimum of 8 neighbour points (N min ) with known pressure considering a 3 × 3 × 3 kernel around the point of interest. Selecting a smaller value of N min results in a faster propagation of the integration front, whereas a larger value adds redundancy, as the local pressure estimate is built from a larger ensemble. Its effect on the pressure integration procedure is discussed for the 2D case in the work of (van Oudheusden 2008). For the selected points, the pressure gradient ∇p is integrated from the available neighbors x i to the evaluation point x e : With each integration step, the domain of known static pressure is dilated. Conversely, the part of the rotational domain where the pressure has not been integrated yet is eroded.
A 2D visualization of the spatial marching procedure is given in figure 8, accompanied by an animation in the online appendices of this article (video 1). The spatial marching algorithm advances its fronts until the static pressure is evaluated at all cells, yielding the pressure distribution in the flow field.

Surface pressure mapping.
For the evaluation of the pressure at the object surface, few markers are placed on the object that are detected by the measurement system and accurately position the model within the measurement domain. Once the model position is known, its surface is discretised at the desired spatial resolution, yielding a description of the surface by a point cloud x s . Here, a conventional surface mesh with triangular faces, stored in a stereolithography (.stl) file is considered (see e.g. Kai et al 1997), but any equivalent type of discrete surface description may be used.
The value of the pressure on each point of the discretized surface x s is obtained from its nearest neighbour in the flow field (Cartesian) mesh. The nearest neighbour is determined by considering the geometrical distance (Euclidean). Points within a 3 × 3 × 3 kernel centred around the nearest neighbour point are used to compute the static pressure at the surface node. The latter is determined following a similar approach as outlined in equation (6), namely by integration of the pressure gradient. The principle is schematically outlined in figure 9.
In this last part of the process, the pressure gradient from the nodes is used, rather than using the trapezoidal type integration used in the flow field. Therefore, ∆p i−e from equation (7) is changed into The pressure mapping is carried out for all surface nodes that have neighbouring fluid mesh cells with valid static pressure values, yielding the surface pressure distribution on the object of interest.

Wind tunnel experiments
The wind tunnel experiments presented in this work are conducted in the Open Jet Facility (OJF) at TU Delft. The closedcircuit wind tunnel features an octagonal 3:1 contraction, feeding a 2.85 × 2.85 m 2 test section. Operations are at atmospheric conditions, and the turbulence level in the nominal wind speed range of 4-35 m s −1 is within 1%.
Here the experimental apparatus and procedures for the sphere-flow case are described, whereas the system and setup for the cyclist experiment have been detailed in (17).
Tracer particles are generated by a 10-wing, 200-nozzle HFSB seeding rake installed in the settling chamber, downstream of the turbulence meshes. An individual nozzle in nominal working conditions produces 20 000 to 60 000 bubbles per second of 300 to 550 µm diameter (Scarano et al 2015, Faleiros et al 2018. Nozzles are spaced 5 cm apart, yielding a seeded streamtube of 0.5 × 1.0 m 2 cross section. A fluid supply unit (FSU) from LaVision GmbH controls the supply of air-, helium-and soap-to the generators.
The LaVision MiniShaker Aero coaxial volumetric velocimeter is installed on a Universal Robots UR5 collaborative robotic arm to form the robotic volumetric PTV system. The CVV probe comprises four CMOS imagers in a diamond arrangement, integrated with an optical fiber for light transmission in the center, see figure 10. A Quantronix Darwin-Duo Nd-YLF laser (2 × 25 mJ pulse energy at 1 kHz) serves as a light source for illumination (λ = 527 nm). The active sensor counts 640 × 476 pixels and it acquires image quadruplets at a rate of 821 Hz. The technical specifications of the imaging system are listed in table 1.   The above systems are used to observe the airflow around a 3D printed sphere of diameter d = 15 cm at 10 m s −1 free-stream velocity. The corresponding Reynolds number Re d = 9.9 × 10 4 , obeys the sub-critical flow regime (Achenbach 1972), where the boundary layer remains laminar until flow separation. The sphere is equipped with 15 static  pressure taps in the vertical symmetry plane, between −10 • and 130 • azimuth at an equidistant spacing of 10 • . The orifices feature 1 mm diameter. The sphere is supported from behind by a 30 mm diameter, circular, stainless steel tube. A vertical beam holds the support from 1.2 m downstream. The sphere sits 1 m above the ground. The experimental apparatus is shown in figure 11. The pressure tap measurements serve as a reference for the surface-pressure evaluations from the robotic volumetric PTV recordings.
During acquisitions, the CVV probe is oriented normal to the free-stream direction, at approximately 40 cm distance to the sphere center. Streamwise translations of the probe yield data acquisitions up-and downstream of the sphere. Each view contains 10 000 images (T aq = 12.1 s). The measurement volume per acquisition, whose shape is best described by a truncated pyramid, extends 40 cm in depth, spanning 12 × 16 cm 2 at its top end and up to 36 × 48 cm 2 at its base. The total measurement volume, combining multiple acquisitions, exceeds 50 l.
For the selected free-stream velocity and acquisition frequency, a particle tracer displaces by 12.2 mm in the free-stream, respectively 23 px on the camera sensors (at a nominal distance of z = 40 cm). Figure 12 shows a typical recording of particle images, with a particle image density of approximately N p = 0.01 ppp in the image center. Background illumination is eliminated by image pre-processing with a high-pass frequency filter (Sciacchitano and Scarano 2014) and a subsequent mean subtraction, which eliminates any residual background noise.
The particle images are analyzed with the procedure outlined in section 2, consisting of particle tracking in 3D and subsequent ensemble averaging.
For the given parameters, the characteristic flow frequency of 67 Hz is estimated from the sphere diameter D and the free-stream velocity u ∞ . An acquisition of T aq = 12.1 s therefore results in approximately N = 800 statistically uncorrelated samples. Further assuming a maximum fluctuation level in the sphere wake of the order of 30% with respect to the free-stream velocity (48), the 95% confidence level uncertainty in the ensemble averaged mean velocity is about 2% (0.2 m s −1 ).

Results
The sphere experiment is intended to assess the accuracy and resolution of the surface-pressure evaluation approach. The analysis of the cyclist data in section 4.2 shall be regarded as a demonstration of the applicability and scalability of the developed method to geometrically complex and large-scale experiments.

Flow field and pressure around a sphere
The time-averaged velocity field measurement is shown in figure 13 (left) in the vertical symmetry plane of the sphere. The corresponding pressure field as obtained by the outlined method is illustrated on the right side of figure 13.
Flow stagnation in figure 13 (left) is corresponded in figure 13 (right) with a relative pressure increase of 60 Pa, consistent with the flow dynamic pressure in the free-stream (q ∞ = 60 Pa). The flow accelerates around the object up to approximately θ = 80 • , where the onset of separation is observed to occur. In the rear region of the sphere, a circulatory motion is established with reverse flow conditions. The inspection of the 3D data returns a toroidal structure for the recirculation, consistent with the axial symmetry of the object geometry. The pressure distribution exhibits little variation past separation, with an annular region of mildly lower pressure (∆p = −30 Pa) in the separated recirculating flow (θ~150 • ).
The three-dimensional static surface pressure distribution is shown in figure 14, along with contours of the streamwise velocity field on orthogonal symmetry planes. From a careful analysis of the 3D surface pressure data, it is observed that the flow field is not fully axis-symmetric, but the flow acceleration in the horizontal plane appears to be more dominant as compared to the acceleration in the z-x plane. This translates into a further reduction in surface pressure of up to 7 Pa on the side of the sphere (negative y).
The level of symmetry of the surface pressure data can be quantified in figure 15, where the distribution of the pressure coefficient C P is plotted as function of azimuth angle θ and elevation level ϕ. The data shown in figure 15 exhibit some degree of variation along the elevation, which feature a delayed pressure decrease followed by an increased suction around ϕ = 90 • . The maximum disparity narrowly exceeds one contour level, which corresponds to 10% C P for the present data.
Because potential flow theory can be considered as ground truth only in a limited angular range (θ ⩽ 30 • ), the accuracy of the surface pressure distribution resulting from the 3D PTV data is compared to the direct measurements obtained by pointwise pressure taps. The latter are located in the z-x plane, as indicated in figure 13, corresponding to ϕ = 0 • elevation.
The pressure coefficient distribution along radial lines is shown for selected azimuth angles in figure 16, where the value at the sphere surface is compared to the recording from the pressure transducer. For 0 • and 30 • azimuth, the static pressure is increasing towards the sphere surface as expected from the flow deceleration. For larger azimuth angles, a reduction of C P towards the wall is observed instead. For the selected angles, the pressure obtained by robotic volumetric PTV correlates well with the transducer readings for θ ⩽ 60 • , whereas a difference of up to 0.17 C P is observed approaching flow separation.
Lastly, the effect of spatial resolution in the time-averaged velocity field measurement on the computed surface pressure distribution is assessed. The size of the averaging bin is varied in the range from l v = 10 to 50 mm. An overlap factor between adjacent bins of 75% is applied. The results are compared to the reference measurements by the wall orifices in figure 17. Two more solutions are included in the comparison: first, the surface pressure distribution obtained by full spatial integration, starting from a free-stream boundary condition far upstream of the sphere; and second, the value obtained by extrapolating the pressure from the closest point to the surface after solving the Poisson problem. The two latter are evaluated with l v = 10 mm.
For the forward part of the sphere (θ ⩽ 30 • ) potential flow theory is considered a reliable approximation of the flow, which is plotted (solid black) in figure 17 (left) for reference. In this region, the PIV pressure agrees well with both, potential flow theory and surface tap measurements. Differences here are smaller than 0.03 C P (1.8 Pa) for θ < 30 • , even at the coarsest resolution (l v = 50 mm). The Poisson solution instead underestimates the stagnation value by approximately 10%. The pressure obtained by full integration matches well at stagnation, and shows a gradually increasing error until θ = 40 • .
For larger angles, the wall taps provide the reference up to θ = 130 • . In the region of 30 • ⩽ θ ⩽ 60 • , the azimuthal pressure distribution reaches its maximum gradient. The pressure provided by the PIV measurements is lower and closer to the potential flow solution as compared to the tap readings. A maximum deviation of 0.15 C P is found at 40 • . Interestingly, there is no clear dependence of the results from the spatial resolution, except for the coarsest bin size (50 mm), with the overall effect of smoothing the surface pressure distribution along the azimuth.     A small region with maximum suction point is located at θ = 70 • according to the orifice measurements, with the PIV data approximately following it. The subsequent mild pressure recovery (adverse pressure gradient) downstream of this point is observed in the transducer data but not well captured by robotic volumetric PTV. The resulting difference is approximately 0.15 C P (9 Pa) in the separated wake region for the data obtained with the proposed method. The data obtained by full spatial integration shows an increasing mismatch in the wake region up to 0.26 C P .
Considering the above, one may conclude that for the measurement of the time-averaged surface pressure distribution, spatial resolution effects are not observed as long as the interrogation bin remains smaller than half the local radius of curvature of the test model. The latter can be regarded as a design criterion for pressure measurements by robotic volumetric PTV.

Surface flow and pressure over a full-scale cyclist
Aerodynamics in competitive cycling remains a vividly discussed topic, which is studied by means of experimental and computational fluid dynamics (CFD) nowadays. A summary of various works is presented in a review article by Crouch et al (2017). More recently, experimental studies on full-scale cyclists have benefited from the use of large-scale PIV to infer the flow field and even forces acting on the athletes (Spoelstra et al 2019, Terra et al 2019). It is noticed however, that experimental studies in cycling aerodynamics seldom report surface pressure measurements. An exception is the study of Crouch et al (2014), which highlights the challenges in measuring the surface pressure over the back of a cyclist mannequin by conventional wall orifices. Examples of CFD studies indicate interest in the pressure distribution over cycling athletes however, see, for example, the works of Blocken et al (2013), and Beaumont et al (2018), among others. In this context the PIV based pressure evaluation method can provide the suitable tool to measure surface pressure distributions on realistic wind tunnel models experimentally.  Regions of high pressure are clearly identified (figure 18 front-view) on the forward sections of the helmet, the rider's hands and fingers, the upper arms, as well as the extended leg and the right knee where the flow is stagnating. A slightly asymmetric pressure distribution is observed on the two arms, which is less easily explained than that occurring on the legs, where the stretched left leg faces the flow more bluntly as compared to the flexed right leg. Regions of low pressure are  Both experiments are carried out in the low-speed, incompressible flow regime, with the Reynolds number in the order of Re = 10 5 , depending on the length scale which is considered in case of the cyclist. The estimated uncertainty level for the mean velocity compares well. A clear difference is observed in the variability of the free-stream pressure, σ CP,tot , which is approximately three times higher for the cyclist case. This is attributed to the fact that the seeding generator in the cyclist case is placed downstream of the wind-tunnel contraction rather than inside the settling chamber as in the sphere case, justifying an increased turbulence level in the free-stream. This translates into a higher bias error resulting from the domain partitioning described in section 2.3, adding a bias error of 2.7% in total pressure coefficient, over the 1.5% for the case of the sphere.
A clear difference is recognized in the variety of length scales between the two cases: the sphere is uniquely characterized by its radius, whereas the athlete model features a range of scales as indicated in table 2. In the analysis of the sphere data (section 4.1) it was concluded that spatial resolution effects are not observed however, as long as linear bin-size l v remains smaller than half the local radius of curvature R. Assuming the same scaling holds for the cyclist, it is seen that the ratio of l v over R is below 0.5 for the majority of the athlete model (legs, shoulders, arms, helmet, torso) and only locally, e.g. at the ankles, the fingers, or the facial details, the ratio exceeds the critical value of 0.5. Therefore, it is assumed that, with exception of the aforementioned regions, the confidence in the surface pressure on the athlete model compares well to the sphere-case, where an uncertainty level of 0.15 C P (9 Pa) has been established. Accounting for the above-mentioned bias due to the pressure-variability in the free-stream it is expected that the surface pressure distribution on the athlete model is accurate to within 0.2 C P (approx. 24 Pa, based on the freestream dynamic pressure q ∞ = 0.5ρu 2 ∞ ). A further discussion on the uncertainty of the surface pressure measurements is presented in the appendix, where the pressure uncertainty is evaluated from the variability of the results obtained from different integration paths.

Local pressure-distribution.
Three distinct regions are addressed in more detail in the following with the purpose to understand the level of detail that can be represented with the present measurements: the athlete's face, his lower right leg, and the lower back. The face and head region presents a geometry with limited optical access. The right ankle-foot region features the source of a number of vortices, which dominate the local flow topology. Lastly, at the lower back a region of flow separation is expected. Figure 19 displays the pressure distribution in the athlete's head and face region. Areas of high pressure are prominent on the forward facing fingers, the centre of the helmet and the  biceps. Here, a stagnation line can be identified which divides the flow into portions passing on the inside and outside of the upper arm respectively. On the helmet instead, there is a stagnation point from which streamlines emanate radially. As the flow passes the athlete, the static pressure reduces towards the widest sections of hands, helmet and upper arm. On the side of the helmet, the pressure reduction is followed by a rapid pressure increase (adverse pressure gradient) due to the flow decelerating in front of the biceps and shoulder. Negative pressure coefficients, indicating local suction, are visible on the outside of the elbow, and past the shoulder onto the curved upper back.
The foot-ankle region is shown in figure 20. On the upstream section displayed in figure 20 (left), the highest pressure is observed at the knee and the forward facing instep of the shoe. A significant pressure reduction is seen on the outside of the shoe, while around the inclined shin, the pressure gradient is less steep. The strong pressure variation on the foot can be linked to the presence of a number of (counter-rotating) vortices originating from the shoe and the ankle region as indicated in figure 20 (centre). Viewing from downstream (figure 20-right) there is a minimum pressure region just above the ankle, separating two moderate pressure regions on the heel and the calf. The low pressure region corresponds to the location where the two shear layers from the lower leg roll-up into two streamwise, counter-rotating vortices, causing a local pressure reduction.
When discussing the surface flow topology on the athlete's extremities, a word of caution is appropriate: for the surface pressure evaluation we assume a known geometry of the model as stated in section 2.4. In the case of the cyclist this is the surface model employed to manufacture the mannequin by 3D printing (van Tubergen et al 2017). While the manufacturing process can be regarded accurate, with the layer-resolution in fused deposition modelling (FDM) typically being below 1 mm (Ultimaker BV 2019), model deviations can be introduced when installing the mannequin on the bike. In the specific case of the legs, a small degree of bending may occur when the feet are fixed to the pedals, which can result in a displacement of several millimetres compared to the reference geometry. An additional source of uncertainty is therefore introduced in the evaluation of surface flow properties, which has not been considered in the discussion in section 4.2.2. This type of uncertainty can be alleviated if the final model geometry is determined experimentally, e.g. by means of photogrammetry or 3D-scanning procedures.
Lastly, figure 21 features the static pressure distribution on the athlete's lower back. Furthermore, an iso-surface of u = 0.1 · u ∞ is plotted as an indicator for quasi stagnant flow. The near-surface streamlines indicate a complex flow pattern: the flow over the hips contains a clear upward component, resulting into an in-wash motion over the lower back. Both hip regions feature low static pressure, especially on the left-side of the stretched leg. The streamlines entering the lower back region past the hips are drawn into rotational flow patterns on either side of the back, see figure 21(a) which are linked to the so-called hip-vortices as first reported by Crouch et al (2014). The two hip-vortices can also be seen as a driver for the inwash over the athlete's back and the down-wash in his wake. Below the right hip-vortex core, indicated by point (c)

Conclusions
A novel pressure-evaluation approach is presented to examine static surface pressure distributions on generic three-dimensional geometries by means of ensemble averaged 3D PTV data as obtained by robotic volumetric PTV. The pressure field is retrieved in a dual-model approach, exploiting Bernoulli's equation in the region of constant total pressure, and integrating the mean pressure gradient elsewhere. The uni-directional integration of the pressure gradient avoids the specification of boundary conditions on the fluid-solid interface, which often deteriorates the result following the Poisson approach.
Experiments on a 15 cm diameter sphere model in a 10 m s −1 airflow return pressure distribution in good agreement with potential flow theory as well as pressure taps. Differences up to 0.15 C P are observed in the region of maximum pressure gradient and in the separated wake. It is observed that spatial resolution of the ensemble averaged velocity field appears not to be critical as long as the linear bin size is smaller than half the local radius of curvature.
Data reduction of the flow field data recorded around a fullscale cyclist model at 14 m s −1 yields the static pressure distribution on the athlete model. The pressure evaluation by robotic volumetric PTV is particularly suited to such applications, given the model complexity and the use of a mannequin wearing sports garment. The 3D surface pressure data has allowed to further understand the flow behavior previously inspected in terms of streamlines topology and vorticity structures.
Measurement accuracy relies on the constant total pressure field hypothesized at the inflow boundary. Secondly a known model geometry and position are also assumed. The former is common practice in wind tunnel experiments, whereas the second assumption may not hold true when using flexible models or bending under the action of the aerodynamic forces is observed.

Appendix A: Cyclist-pressure uncertainty analysis
In addition to the analysis presented in section 4.2.2, an a-posteriori uncertainty assessment is carried out, based on  variations in the integration path and the boundary condition at the rotational-irrotational domain interface.
First, the measurement domain, schematically illustrated in a top-view in figure A1(a), is partitioned following the procedure detailed in section 2.3, resulting in a two-part domain as shown in figure A1(b). Evaluation of the pressure based on the theory outlined in section 2.4 yields the results presented in section 4.2, which are herein used as the reference.
Secondly, the domain interface shown in figure A1(b) is split along the median-plane (z = 0 mm), yielding the partitions illustrated in figures A1(c) and (d), where the irrotational flow assumption is only applied on one side of the measurement domain, and the opposite side is instead spatially integrated. In this way, the extent of spatial integration is significantly increased, whereas the Dirichlet boundary condition for the integration process in cases (c) and (d) are independent of each other. In a third and final step, the surface pressure obtained on the athlete's left-hand side (LHS, z > 0 mm) for case (c) and his right-hand-side (RHS, z < 0 mm) for case (d) are compared to the reference result, indicating the variability of the resulting surface pressure with the applied change in boundary condition. The results are graphically compared for case (d) in figure A2.
Pressure differences are predominant only on the far-side with respect to the side where pressure integration is started (RHS for case (d); and LHS for case (c)), as the integration path to reach the far-side is significantly longer compared to the reference case. As such, the pressure differences on the extremities (legs and arms) are larger as compared to the athlete's back or helmet, which are closer to the median plane. Furthermore, it is seen that there is no clear bias, but both, higher pressure (red) is observed on the right leg and the downstream facing side of the upper arm, and lower (blue) surface pressures are seen on the right hip, the upstream facing upper arm and the right-hand-side of the face and helmet. Lastly, the statistical distributions from both cases (LHS case (c) and RHS case (d)) are combined in a histogram of the pressure variation ∆C P = C P,i -C P,ref in figure A3.
The data in figure A3 reveals a mean delta of 3% C P , with a root-mean-square value of 23%, which corresponds approximately with the hypothesized uncertainty of 0.2 C P following the scaling between sphere and cyclist-experiment as shown in section 4.2.2, thereby confirming the a-priori uncertainty estimate.