The slip velocity of nearly neutrally buoyant tracers for large-scale PIV

The behaviour of nearly neutrally buoyant tracers is studied by means of experiments with helium-filled soap bubbles and numerical simulations. The current models used for estimating the slip velocity of heavy micro particles and neutrally buoyant particles are reviewed and extended to include the effect of unsteady forces and particle Reynolds number. The particle motion is analysed via numerical simulations of a rectilinear oscillatory flow and in the flow around an airfoil within a particle flow parameter space that is typical of large-scale PIV experiments. An empirical relation is obtained that estimates the particle slip velocity, depending on the particle-to-fluid density ratio, the particle Reynolds number and frequency of the local flow fluctuations. The model developed is applied to assess the slip velocity of helium-filled soap bubbles in a large-scale experiment conducted at the German–Dutch wind (DNW) tunnels in the flow around an airfoil, with chord Reynolds numbers up to three millions. Furthermore, a procedure is proposed that can be used to retrieve the bubbles mean density and dispersion from measurements of mean velocity and fluctuations, respectively.


Introduction
Particle image velocimetry (PIV) is generally deemed as a non-intrusive technique for measuring the instantaneous flow velocity, in contrast to probe-based techniques. However, the tracer particles must follow the flow faithfully for accurate measurements, a crucial requirement often taken for granted.

3
186 Page 2 of 24 Analysis of particle tracking accuracy originates from the pioneering work of Stokes (1851), in which he derived the equation for the drag force of a rigid sphere in a viscous flow at very small Reynolds number values, the so-called Stokes law. A few decades later, Boussinesq (1885a, b) and Basset (1888a, b) expanded the Stokes' analysis to the more general case of unsteady motion, with the particle velocity being an arbitrary function of time, however, still under creeping flow conditions. There have been several works following Stokes, Boussinesq and Basset, attempting to relax the assumption of creeping flow to achieve a more general dynamic equation that would describe the motion of particles with small, yet finite, particle Reynolds number. Maxey and Riley (1983) give a review of the most relevant developments in the century following Boussinesq's and Basset's papers and provide a more generic equation for a small rigid sphere, including the forces due to non-uniform flow. However, most useful equations up to date still rely on semi-empirical corrections (especially for the drag coefficient), resultant from experiments, simulations and analytical expansions of asymptotic solutions (Mei 1996;Magnaudet 1997;Michaelides 1997;Loth and Dorgan 2009).
Based on analytical models, Melling (1997) analysed the tracking characteristics of typical particles used for PIV, restricted to the case of particles much heavier than air, in which the Boussinesq-Basset's equation is greatly simplified. From these theoretical considerations, it was concluded that typical PIV particles with density in the range of 1000-4000 kg/m 3 should not exceed 1 µm in diameter for a frequency response of 10 kHz. Particle relaxation length and time were obtained across shock waves by Scarano and van Oudheusden (2003) in the supersonic regime and by Schrijer et al. (2006) up to Mach 7, returning a response time of submicron TiO2 particles in the range of 2 µs. A systematic study of micron-sized flow tracers conducted by Ragni et al. (2010), comprising several particles used for PIV, reported a time response of 2 µs for di-ethyl-hexylsebacat (DEHS) droplets of 1 µm median diameter. Based on the study of Samimy and Lele (1991), it can be inferred that velocity measurements with micrometre particles of 2 µs time response result in errors smaller than 2% for flow frequencies up to 100 kHz.
A known drawback of micron-sized tracers is the limited amount of light scattered from laser illumination, due to the small scattering cross-section (Adrian and Yao 1985), rendering them optically ineffective when the scale of the experiment is increased. For instance, three-dimensional measurements with tomographic PIV seldom exceed 100 cm 3 (Scarano 2013).
The amount of light scattering can be increased without significantly compromising the tracing fidelity by adopting larger particles with reduced weight, such that the particles become approximately neutrally buoyant. If the velocity gradients at the particle scale are negligible, spherical neutrally buoyant particles follow the flow perfectly (Mei 1996). In air, a particle of neutral density is obtained through the generation of soap bubbles filled with a light gas, usually helium (Hale et al. 1969;Bosbach et al. 2009) of about 0.5 mm diameter. Although particles cannot follow velocity fluctuations occurring at turbulent length scales smaller than their diameter, even if neutrally buoyant (Xu and Bodenschatz 2008), submillimetre helium-filled soap bubbles (HFSB) have been shown to accurately retrieve the mean velocity and turbulent fluctuations for wall distances larger than twice their diameter (Faleiros et al. 2018). The use of submillimetre HFSB has been reported to reflect 10 4 to 10 5 times more light than micrometre particles (Caridi 2018). In the past two decades, HFSB have been applied to several PIV experiments at metre scale: aircraft cabins (Bosbach et al. 2009), wind turbines , thermal plumes (Huhn et al. 2017), aeronautics , full-scale cyclists (Jux et al. 2018;Terra et al. 2019;Spoelstra et al. 2019), aeroacoustics (Lima Pereira et al. 2020) and bird flight (Usherwood et al. 2020).
HFSB can be produced as neutrally buoyant tracers by careful control of the supply rate of their constituents. However, small deviations from neutral buoyancy always occur in practice, leading to a slip velocity between the particle and the flow. Kerho and Bragg (1994) were the first to attempt to characterize the tracing fidelity of HFSB for aerodynamic experiments, measuring the slip velocity along the bubble trajectories in the flow around an airfoil. Errors of up to 10% of U ∞ were ascribed to a bubble generation process biased towards lighter-than-air particles. Based on dimensional considerations, Scarano et al. (2015) devised a practical approach for measuring the HFSB time response in the stagnation region of a cylinder. Although not explicitly mentioned, comparison to the equation of particle motion, as given by Mei (1996), shows that the approach assumes approximately equal acceleration between particle and fluid. Applying the procedure experimentally, a mean time response of approximately 10 µs was measured for submillimetre HFSB. The methodology was followed by many authors for the study of HFSB tracing fidelity (Morias et al. 2016;Gibeau and Ghaemi 2018;Faleiros et al. 2019;Gibeau et al. 2020). Morias et al. (2016) reported measurements of the standard deviation of HFSB time response of 40-50 µs, exceeding its mean. They also observed that the bubbles could be produced in two regimes, bubbling and jetting, with only the former resulting in bubbles of uniform size. Faleiros et al. (2019) performed a rigorous assessment of such regimes for several input values of helium, air and soap volume flow rates, resulting in general guidelines for the generation of HFSB in the bubbling regime, and the control of bubble size, density and production rate. A mean time response below 20 µs was found, corresponding to density deviations smaller than 10%. The time response dispersion of about 40 µs was found to be independent of the regime of production. Gibeau and Ghaemi (2018) and Gibeau et al. (2020) performed similar experiments and also obtained mean time responses below 10 µs for bubbles approaching neutral buoyancy; however, standard deviations as high as 180 µs were reported. The dispersion of the physical properties of HFSB tracers could therefore represent a more critical phenomenon than matching the ensemble average density of the tracers to the air density. Furthermore, the tracing fidelity of HFSB has also been verified in isotropic turbulence (Qureshi et al. 2008;Bourgoin et al. 2011) and wall-bounded turbulent flows (Faleiros et al. 2018), with the results showing accurate measurements of turbulent stresses, even for air-filled soap bubbles.
In most PIV experiments, it is desirable to know the expected particle slip velocity for estimating the measurement uncertainty. Although the method devised by Scarano et al. (2015) is adequate for time response estimation when the particle acceleration does not depart significantly from that of the surrounding fluid, it may not hold for many practical situations. For instance, in the work of Morias et al. (2016), when a 30% heavier-than-air HFSB approaches the cylinder stagnation point (within 10% of the cylinder diameter from the leading edge), the mean particle deceleration is up to twice as large as that of the air flow. The HFSB slip velocity has not been investigated when the particle acceleration deviates considerably from that of the flow. In fact, current research on HFSB focused mostly on the generation of neutrally buoyant bubbles rather than the estimation of the slip velocity of nearly neutrally buoyant tracers. The timeline in Fig. 1 summarizes the main investigations on HFSB tracing fidelity, focused on its quantification and control.
The present work surveys the available models for estimating the slip velocity of nearly neutrally buoyant particles and provides a generalized model within a particle flow parameter space that is typical of large-scale PIV experiments conducted with HFSB. The approach expands from simple predictions of HFSB flow tracing accuracy based on the Stokes regime, or other similar simplifications that neglect the unsteady forces, towards a model that incorporates the latter and allow for large acceleration differences.
For this purpose, the role of nonlinear terms in the equation of motion is examined. The slip velocity of HFSB cannot be simply estimated from a single time response parameter, as in the case of heavy small particles (Melling, 1997), because the unsteady forces also become important and bring additional complexity to the analysis. Although the slip velocity can be promptly calculated by realizing a numerical simulation (similar to that carried out in this paper), the intention is to provide means for judging the experimental velocity errors without the need of such timeconsuming computations. This is performed by analysing the motion of nearly neutrally buoyant particles through numerical simulations in a rectilinear oscillatory flow and in the flow around an airfoil with focus on the high-acceleration region in the vicinity of the leading edge. From these analyses, an extension of the current methods used to estimate the slip velocity is proposed. Fig. 1 Timeline of main contributions regarding control of HFSB generation (focused on the generation of neutrally buoyant bubbles) and quantification of their tracing fidelity. A complete overview of HFSB generation would also have to include the works of Hale et al. (1969), Okuno et al. (1993), Caridi et al. (2016), Gibeau et al. (2020), among others Additionally, a set of experiments in a large-scale aerodynamic wind tunnel is conducted to examine the tracer behaviour of HFSB at high-Reynolds-number flows, under similar conditions than the simulations. The potential of the proposed model to assess PIV measurement errors using HFSB is demonstrated via application to experimental data, from which the HFSB density is retrieved, helping to identify the source of velocity errors. In addition, the effects of the HFSB density dispersion, herein proposed as the main source of the previously measured time response dispersions (Morias et al. 2016;Faleiros et al. 2019), are also investigated.

The equation of particle motion
For a spherical heavy particle, the equation of motion, neglecting buoyancy force, reads as (Melling 1997): where ⃗ u slip = ⃗ u p − ⃗ u is the slip velocity, ⃗ u p and ⃗ u and are the particle and fluid velocity, respectively, p = p d 2 p ∕18 is the particle time response, p is the particle density, d p is the particle diameter and is the fluid dynamic viscosity. The term on the right-hand side is commonly referred to as Stokes' drag (Stokes 1851). Equation (1) ∕ is the particle Reynolds number, is the fluid density and is the fluid kinematic viscosity. Although the applicability of Eq. (1) seems limited, it describes accurately the motion of micrometre heavy particles in air, and it is commonly used for evaluating the tracing fidelity of PIV seeding (Melling 1997;Scarano and van Oudheusden 2003;Schrijer et al. 2006;Lazar et al. 2010;Ragni et al. 2010).
In the case of large nearly neutrally buoyant particles, the unsteady forces become relevant and the Boussinesq-Basset's equation with the corrections for finite Re p must be applied. The full unsteady equation of motion as proposed by Mei (1996) is: where ̂= p ∕ is the particle-to-fluid density ratio ( . is used throughout the paper for representing normalized variables), ⃗ g is the gravitational acceleration, Re p is an empirical relation to correct for deviations from the Stokes' drag law due to a finite Re p , K(t − ) is the history force kernel, d∕dt = ∕ t + ⃗ u p ⋅ ∇ is the time derivative on the particle trajectory and D∕Dt = ∕ t + ⃗ u ⋅ ∇ is the time derivative evaluated on the trajectory of fluid elements around the particle. Rigorously, the lower limit of integration of the history force is -∞. For simulation purposes, t − is the instant right before the particle is introduced into the flow. Furthermore, it is noted that the buoyancy force is usually negligible in external aerodynamic experiments, given that the acceleration due to the flow kinematics substantially exceeds the gravitational term, and is henceforth omitted.
In the study of aerodynamic flows, normalization of the equation of particle motion is usually based on a chosen characteristic velocity U 0 and length L 0 , with the reference timescale being derived from the latter. For instance, in airfoil flows the free stream velocity and the airfoil chord (or airfoil thickness) are typical choices. Such normalization is convenient in Eq. (1), as it results in a single non-dimensional parameter, the so-called Stokes number S k ≡ p U 0 ∕L 0 . Normalization of Eq. (2) based on U 0 and L 0 leads, however, to three non-dimensional parameters, S k , ̂ and Re p , increasing the complexity of the particle tracing fidelity analysis. Alternatively, Eq. (2) may be normalized based on a chosen reference velocity U 0 and the viscous timescale t 0 ≡ d 2 p ∕ : where the history force kernel may be calculated using the approximation (Mei 1994) with Re p evaluated at t −̂ . Notice that ̂ and Re p are the only non-dimensional parameters in Eqs. (3) and (4). Furthermore, although Re p represents the actual particle Reynolds number, often it is desirable to specify a reference Reynolds number that is independent of the slip velocity. For this purpose, a diameter-based Reynolds number

Drag correction
Several empirical relations have been proposed for the drag correction term Re p = C d ∕C d0 (Clift et al. 1978), where C d is the particle drag coefficient and C d0 = 24∕Re p is the drag coefficient of a rigid sphere given by the Stokes' law ( Re p ≪ 1 ). A commonly used expression given by Schiller and Naumann (1933) for which the drag coefficient is accurate within 5% for Re p ≤ 800: More accurate approximations were compiled by Clift et al. 1978, for different Re p ranges: where w = log Re p . There are different expressions for the drag coefficient in the case of a clean bubble (the flow slips along the surface) or due to particle deformation (Clift et al. 1978;Magnaudet and Eames 2000;Loth 2008). In this study, it is assumed that the particles are rigid non-deformable spheres, in which the no-slip condition at the surface applies. A soap bubble of 0.5 mm diameter in air maintains its sphericity (less than 10% difference in size between the two spheroid semi-axes) for slip velocities up to about 10 m/s (Faleiros 2021). The no-slip condition applies if the surface is sufficiently contaminated with surfactants. A soap bubble can be considered fully contaminated if the contaminants concentration is larger than 10 -2 g/l (Loth 2008). The amount of surfactants used to produce soap solutions for HFSB is far greater than this threshold (Hale et al. 1971).
Comparison of Eq. (5) and (6) (Fig. 2) shows that, up to Re p = 800 , the corrections differ less than 5% from each other. However, as Re p surpasses this limit, Eq. (5) underestimates drag significantly, with the estimate from Eq. (6) being 50% larger than that given by Eq. (5) at Re p = 5000 . In most PIV applications, Eq. (5) is sufficient, since Re p = 800 represents a slip velocity of 24 m/s for a half-millimetre particle in air. Furthermore, Fig. 2 shows that the assumption of Re p = 1 is limited to low-speed flows and small particles. The correction term is already twice this value ( Re p ≈ 2 ) for a slip velocity of 0.4 m/s under the same conditions.

Faxén terms and lift force
Strictly speaking, the flow velocity and acceleration should be calculated as averages along the particle surface ⃗ u S and over its volume D⃗ u∕Dt V , respectively. Faxén (1922), as cited in Michaelides (1997), derived approximations to these averages using a Taylor series expansion in the particle limit d p → 0 for ⃗ u ≈ ⃗ u p , resulting in extra terms that contain the Laplacian of the velocity (Mei 1996), the socalled Faxén terms. However, in the numerical simulations performed in this study, the flow field is derived from a potential flow solution. In an incompressible ( ∇ ⋅ ⃗ u = 0 ) and irrotational ( ∇ × ⃗ u = 0) flow, the Laplacian of the velocity is equal to zero. In fact, even in rotational flows the Faxén terms are usually negligible in comparison with the remaining terms (Mei 1996), only becoming relevant when the slip velocity approaches zero, and the effect of drag and history force have the least impact on the net force (Calzavarini et al. 2009). Therefore, they have been omitted from Eq. (2).
Additionally it is noted that the lift force has also been neglected. Shear lift or spin lift induced by shear is not relevant for the flow around an airfoil ( ⃗ = ⃗ 0) , only in Fig. 2 Comparison of drag correction given by Eqs. (5) and (6) the boundary layer. The influence of wall lift (Zeng et al. 2009), which exists even in a shear free flow, is only important at distances from the wall L W of O(d p ) and for low Re p . From the empirical relations given by Zeng et al. (2009), as reviewed by Shi and Rzehak (2020), it is estimated that for L W ∕d p > 2 the ratio of lift to drag coefficient C l ∕C d is less than 1% for Re p > 25 . In the simulations performed in Sect. 3.3, the particles approach the airfoil at 50 < Re p < 200 . Thus, wall lift can also be safely neglected here.

Simplifications to the equation of particle motion
In this section, Eq. (2) is studied under different assumptions in order to reproduce previous models available in the literature for the analysis of particle slip velocity (Sect. 2.3). Notice, however, that in numerical simulations presented in chapter 3, Eq.
(2) is used as presented.
If, for the sake of the discussion, the buoyancy and history forces are neglected, Eq. (2) becomes: In the study of micrometre heavy particles (Melling 1997) typically used for PIV, ̂≫ 1 and Eq. (7) is greatly simplified to: As Re p → 0 , Re p → 1 , and Eq. (8) returns to Eq. (1). This simplified form enables analysis of the slip velocity based only on particle acceleration and time response. Other simplifications to the equation of motion have attempted to enable similar analysis without the restriction of ̂≫ 1 . If the added mass force is also neglected, Eq. (7) becomes: In addition, for flow regions where the particle has approximately the same acceleration as the flow, the slip velocity may be written in a similar fashion to Eq. (1): Re p is the particle time response typically used in the study of HFSB tracing fidelity (Scarano et al. 2015;Morias et al. 2016;Faleiros et al. 2019;Gibeau et al. 2020). Similarly to Eq. (1), Eq. (10) provides an operational approach to the measurement of the particle time response through the simultaneous measurement of the slip velocity and the flow acceleration. Furthermore, if the particle size is also measured, e.g. by estimation of the bubble glare points distance (Morias et al. 2016;Faleiros et al. 2019), Eq. (11) yields an estimate of the particle density, a notably difficult parameter to measure.
In comparison with p , * p is proportional to the density difference p − , representing more appropriately the time response of nearly neutrally buoyant particles. However, it is emphasized that basing the analysis of HFSB tracing fidelity only on * p is an oversimplification that is only valid if the particle acceleration is approximately equal to the fluid acceleration.

The slip velocity in an oscillating flow field
It is instructive to first consider the case of a spherical particle immersed in a rectilinear oscillating flow field. The respective velocities are given by u p (t) =ũ p ( )e −i t and u(t) =ũ( )e −i t , where ũ p ( ) and ũ( ) are the particle and flow amplitudes of the oscillation, and = 2 f is the angular frequency.
Assuming negligible history force and added mass and Re p = 1 , which is acceptable for low Re p , substitution of the particle and flow oscillating velocity and acceleration in Eq. (9) gives: From Eq. (12), the particle tracing fidelity can be analysed from the ratio of the slip velocity to the flow velocity amplitudes: Following Mei (1996), u R , p ,̂ may be rewritten as a function of ̂ and only, where is a non-dimensional parameter with similar significance to the Stokes number: It can also be easily verified that if the added mass had been included, Eq. (7), then the amplitude ratio becomes: Furthermore, the analysis can be expanded by considering an approximation to the history force in the oscillating flow field for ≫ 1 (Mei 1996): Including added mass and history force, the amplitude ratio becomes: where the history force is represented by the term (1 − i).
For the purpose of illustration (Fig. 3), the tracing behaviour of lighter-( ̂ = 0.5) and heavier-than-air ( ̂ = 1.5) particles is modelled on the basis of Eq. (14), which neglects both added mass and history forces, and Eq. (18), in which both forces are included. The particle motion modelled by Eq. (14) (Fig. 3, left) exhibits a positive phase lag and some degree of amplitude modulation for the heavier-thanair conditions. In the lighter-than-air case, the situation is reversed: the amplitude of the velocity oscillations exceeds that of the flow and a negative phase lag is observed. These observations are in accordance with numerical simulations performed by Müller et al. (2001). The situation remains qualitatively unaltered when the particle motion is modelled , by Eq. (18), including the added mass and history force term ( Fig. 3, right). However, a notable contraction of the discrepancy between the particle and the fluid motion is observed. This is because the added mass and history forces are proportional to the slip acceleration and act as restoring forces, reducing the acceleration difference, and resulting in a higher tracing fidelity. Analysis of the spectral behaviour of the slip velocity amplitudes (Fig. 4) gives an overview of the role of the unsteady forces. In general, including the added mass force as well as of the history force results in reduced slip velocities for all ̂ and > 1. Furthermore, note that the slip velocity amplitude converges to a finite value for large . At large oscillation frequencies, the acceleration-driven terms are dominant, as |du∕dt|∕|u| = . Therefore, the Stokes' drag force, being slip velocity-driven, becomes negligible at high frequencies in comparison with the other forces. When added mass and history force are neglected, it is found from Eq. (14) (14), (16) and (18). The absolute value | | u R | | is multiplied by sign(1 −̂) for clarity of visualization. Lighter-than-air particles: positive side of the vertical axis. Heavier-than-air particles: negative side true for both Eq. (16) and (18), because the history force scales with -while the other forces scale with 2 -not contributing, therefore, to | | u R | | in the limit of → ∞ . The latter cannot be observed in Fig. 4 due to the short range of that is plotted.
This analysis shows the importance of including the unsteady forces in the study of nearly neutrally buoyant particle motion. Estimates of tracing fidelity based only on Stokes' drag become increasingly conservative as the frequency of the flow fluctuation increases, yielding overestimated predictions of the slip velocity. In addition, the estimates from Eq. (18) should also be considered conservative due to the assumption of Re p = 1 . The slip velocity is further reduced when drag correction is considered-see Eq. (8) and (10).

The slip velocity: from sinusoids to the flow around an airfoil
In the case of a particle moving around an object, the analysis of the slip velocity becomes geometry dependent and, therefore, case-specific. The problem can be generalized by making an analogy between the velocity modulation of a particle along its trajectory and that occurring in a sinusoid flow. Particle trajectories represent finite oscillations that can be approximated by a sinusoid of defined amplitude and frequency, within a finite time interval. This can be visualized in Fig. 5 for the flow around a cylinder. As the particle moves around the object, the streamwise component of velocity oscillates from a minimum value upstream, to a maximum on top of the cylinder, returning to the same minimum value downstream. In addition, the transverse component fluctuates in a similar fashion, but as if "90° out of phase", that is, the transverse component is maximum/minimum while the streamwise component is zero, and vice versa.
Since the two components of velocity are orthogonally out of phase, the dynamics experienced by a particle travelling around an object are analogous to that in a complex sinusoid. Thus, the streamwise and transverse components of velocity can be thought of as the "real" and "imaginary" parts of the flow. This analogy allows taking advantage of the property that the angular frequency may be obtained locally through the ratio between the absolute values of acceleration and velocity. Thus, a local frequency f and phase can be defined as: The free stream velocity is subtracted from the velocity components, to ensure that the particle velocity oscillates around zero. The flow frequency and phase defined in Eq. (19) are shown in Fig. 6. The approximation proposed is most adequate in the region − ∕2 < < ∕2 . Within this phase band, covering half cycle of a sinusoid, the frequency along a streamline remains approximately constant, similarly to a complex sinusoid.
The hypothesis set forth is that through the simple calculation of f and in the flow around an object, the slip velocity may be approximated through predictions obtained from a sinusoid flow at the same frequency and phase. The slip velocity components may then be obtained as  reducing the problem to finding the amplitude modulation ũ p ∕ũ and the phase shift s = p − for a given particle and flow frequency.

Numerical simulation set-up
The drag correction term Re p is not considered in Sect. 2.3, and the history force was obtained through an approximation only valid for ≫ 1 . In this section, the rectilinear oscillating flow field is simulated using the full equations of motion, i.e. Equations (3), (4) and (6), with the purpose of finding an empirical relation suitable for assessing the slip velocity of HFSB tracers, depending on ̂ , Re d and f = f d 2 p ∕ . The simulated input parameters (Table 1) include a wide range of applications. Considering slip velocities from 0.1 to 10% of the reference velocity (wave amplitude), Re p varies from the Stokes regime ( Re p = 0.01 ) to regions extending that expected for HSFB in subsonic aerodynamics ( Re p = 5000 ). The frequency range is selected such that the asymptotic convergence of the velocity amplitude becomes apparent (see Fig. 4). Note that this occurs at higher frequencies as Re d increases. The density range covers ± 30% density deviations from neutrally buoyant particles, including particles of small density differences, down to 1%, allowing the empirical fits to better capture the particle behaviour near neutral buoyancy.
Neglecting gravity force, and rearranging the terms, Eq. (3) may be rewritten for a 1D flow as: where û = e −ît . The terms C and â are introduced for integration purposes only. The reference velocity U 0 is the wave amplitude, and the reference time is t 0 = d 2 p ∕ . Multiplying Eq. (21) by exp t∕C and integrating it once yields û n+1 p , while integrating it twice yields the normalized particle position at x n+1 p : where the initial conditions of the integration are: The simulation is performed until the particle dynamics reach steady conditions. This is achieved by requiring that the slip velocity and the phase shift vary less than 0.1% for the duration of half a cycle.

Velocity amplitude modulation and phase shift
Curve fits are applied to the results of the simulation, aiming to provide a simple relation for the velocity amplitude modulation ũ p ∕ũ and phase shift s , and, consequently, the slip velocity from Eq. (20), as a function of ̂ , Re d and f . The amplitude modulation is found to be well described by the following relation: Although not valid in all limits, this expression does match a few important expectations: no amplitude modulation is observed for neutrally buoyant particles ( ̂= 1 ) or in the case of zero oscillation ( f = 0 ), i.e. ũ p ∕ũ = 1 in both cases. Furthermore, similarly to the slip velocity estimations given in Sect. 2.3, as f → ∞ , the amplitude modulation converges to a finite value: The results of the simulations for ũ p ∕ũ, including the curve fits given by Eq. (24) are shown in Fig. 7. The visual agreement of the curve fits and the simulation data suggest that the proposed equation is highly accurate within the specified conditions (Table 1). A quantitative measure of goodness of fit is provided by the statistical coefficient of determination R 2 : where the total ( SS tot ) and the residual ( SS res ) sum of squares are defined as: where z is the data being fitted and F is the result given by the curve fit. In general, the results demonstrate an accurate prediction of the amplitude modulation with R 2 > 0.98. The phase shift is more challenging to represent with a general empirical relation, and a less accurate fit is accepted to allow generalization (Fig. 8). The best empirical relation obtained reads as: The expression used to estimate the phase shift becomes zero for a neutrally buoyant bubble, and although it does not converge to proper limits for f → 0 or f → ∞ , it matches the simulated data with reasonable accuracy within the tested range ( R 2 > 0.8 for most cases).

The slip velocity around an airfoil leading edge
The hypothesis set in the beginning of this chapter is tested by performing a numerical simulation in the potential flow (obtained using XFOIL, Drela 1989) around the leading edge of an airfoil at incidence = 14°. The airfoil is a section of a Fokker 100 aircraft wing (model 5-6 with retracted flap in Reinders W 1994)   The data are gridded using square bins of 0.5% chord length. Within each bin, the scattered velocity is fitted using a quadratic function in both spatial coordinates, following the approach discussed in Agüera et al. (2016). The mean velocity is taken as the fitted value at the centre of the bin. For convenience (as it should become clear when other angles of attack are included in the experimental part), the coordinate system in the Eulerian frame of reference is switched to the airfoil coordinates, where x is the chordwise direction and y is the normal to chord direction.
The slip velocity obtained from the numerical simulation is shown in Fig. 9a. In comparison, the slip velocity obtained from Eq. (20), with ũ p and s obtained from Eqs. (24) and (28), is shown in Fig. 9b. The close agreement between the approximation and the simulated data supports the validity of the hypothesis that the flow around an object may be approximated by that of a complex sinusoid. The approximation given in Eq. (18) (Fig. 9c), based on Mei's work, overestimates the slip velocity at higher frequencies, as the particle approaches the model surface. Furthermore, estimating the slip velocity from Eq. (10), which requires negligible acceleration difference and neglects unsteady forces, leads to an overestimation of up to about 3000% of the simulated value (Fig. 9d).
A closer comparison of the empirical estimations herein developed with Mei's work is performed along the line y∕c = 0 in Fig. 10. The empirical expressions (24) and (28) give accurate estimations for f > 1 , being within 20% of the simulated value for 1 <f < 5 . The errors from the estimations given by Eq. (18)

Density estimation
Usually in a PIV experiment with HFSB, the bubble density is tuned using mass flow controllers set at appropriate flow rates that have been obtained through controlled experiments (Faleiros et al. 2019). However, it is good practice to verify this information during the measurements for an assessment of the errors. Hence, a procedure is defined to retrieve the bubbles' density from HFSB measurements for which a reference flow is available (e.g. DEHS measurements). The steps involved are outlined below by considering numerical simulations of HFSB tracers of given density in a reference potential flow solution around an airfoil. A least square optimization is used to retrieve the particle density from Eqs. (20), (24) and (28) presented, by estimating the slip velocity for several density values within ̂= [0.5, 1.5] and comparing with the numerical simulation for ̂= 0.8 . The best match with the simulated particle velocity is obtained by minimizing the sum of squares: Fig. 9 Comparison of slip velocity estimations with the numerical simulation results for a light particle ( ̂= 0.8 ) in the Eulerian frame of reference: a numerical simulation; b proposed approximation, Eqs. (20), (24) and (28); c approximation (adapted) from Mei (1996) in a rectilinear oscillating flow, Eq. (18); d no-slip acceleration assumption, Eq. (10). Notice that the colour scale is different in the latter case  Mei (1996). Right: Increase in the normalized frequency as the particle approaches the model where the subscript est stands for estimated value. This procedure is applied within fc∕U ∞ > 3 (Fig. 11, left) and ∕4 < < ∕2 (Fig. 11, right), which is deemed well modelled by the proposed approximation.
The SS | | | � ⃗ u slip | | | value has two local minimums (Fig. 12), one for ̂< 1 and another for ̂> 1 . The correct local minimum is selected by considering the minimum sum of squares of the streamwise slip velocity, SS û slip . Although the density estimate from the latter is less accurate than given by SS | | | � ⃗ u slip | | | , it evidently distinguishes between lighter-and heavier-than-air particles. From the minimum SS û slip , it is obtained that ̂< 1 . The local minimum in the lighter-than-air part yields an estimated density ratio of ̂e st = 0.81 , only 1.25% larger than the simulated value. This result indicates that this method can be used for relatively accurate estimation of particle density from measurements of the slip velocity.

Velocity fluctuations due to density dispersion
The consequences of HFSB slip velocity dispersion, which was experimentally quantified (Morias et al. 2016;Faleiros et al. 2019) through the standard deviation of * p ( ∼ 40 µs), are herein considered.

Sources of time response dispersion
The two main particle parameters affecting the slip velocity are p and d p , whose standard deviations and d , respectively, are the main sources of dispersion. The generation of HFSB in the bubbling regime is crucial to guarantee low-diameter dispersion. The coefficient of variation of the HFSB diameter CV d = d ∕d p is 3% in the bubbling regime, but as large as 13% in the jetting regime (Faleiros et al. 2019). However, no correlation between diameter and time response dispersions has been observed, with measurements of in both regimes being of the same order. This may be understood through Reynolds averaged decomposition of the time response as defined in Eq. (11), assuming Re p ∼ 1 for simplicity without losing generality. Notice that the measurements mentioned above have been performed in regions of low slip acceleration, validating this discussion. Neglecting second-order terms and assuming constant ̂ , it reads as: Subtracting * p from both sides, the time response dispersion is given as: Therefore, for a neutrally buoyant bubble ( ̂= 1 ) the diameter dispersion does not affect the time response dispersion. In fact, the diameter dispersion only affects if there is a substantial deviation from the neutral buoyancy condition. Even if the mean density of half-millimetre bubble deviates 10% from the fluid density, would still be 5 µs in the bubbling regime. In the worst case scenario, where jetting regime is present, then ̂= 1.1 results in ≈ 25μs . This supports the experimental observations that HFSB diameter dispersion is not the main drive causing time response dispersion.
Bubble density dispersion may occur independently of bubble size through variations of the soap film thickness. Direct measurements of HFSB density have not been taken Page 13 of 24 186 so far, only indirectly through measurement of bubble size and time response (Morias et al. 2016;Faleiros et al. 2019). Assuming exclusive density dispersion, the Reynolds average decomposition yields: where CV = ∕ p is the density coefficient of variation. Thus, for a distribution of bubbles with a mean density equal to that of the flow, constant diameter d p = 0.5 mm and time response dispersion of = 40 µs, the density coefficient of variation is CV = 0.043 (in air at NTP).
This effect may be translated to soap film thickness variation, through mass conservation analysis. Assuming the volume of the soap film to be much smaller than that of helium, such that the bubble volume equals the helium volume, it is possible to estimate the bubble thickness as: This means that for d p = 0.5 mm, p = air and soap = 1124 kg/m 3 , the soap film thickness is on average 77 nm. Additionally, if the changes in thickness are exclusively responsible for the density dispersion, then the film thickness standard deviation t is about 3.5 nm for CV = 0.043 . Therefore, a film thickness coefficient of variation ( CV t = t ∕t ) of only 5% is enough to result in 40 µs dispersion of the time response. Although this analysis remains to be verified, it is plausible to assume that the soap film thickness varies during the process of bubble formation with a coefficient of variation in the same order of magnitude as that measured for the bubble diameter.

Velocity fluctuation estimation
The time response dispersion results in velocity fluctuations that may be falsely interpreted as turbulence. In a two-component PIV measurement, the flow turbulence intensity, is inferred from the root mean square (RMS) of the particle velocity fluctuation, whose streamwise component (and similarly the transverse component) reads as where u � p = u p − ⟨u p ⟩ , E is the expected value and P is the probability of the outcome u p . In steady flows, u p,rms may be rewritten as: For an error assessment of this effect, the slip velocity may be estimated from Eq. (20), (24), (28), while the probability P may be obtained by assuming that the slip velocity dispersion is exclusively resultant from the particle density distribution. Assuming a Gaussian distribution N ∼ p , 2 , P p is calculated as where F p is the cumulative distribution function: The accuracy of this approximation for estimating u p,rms is demonstrated by repeating the simulation of Sect. 3.3 for a particle of normally distributed density with ̂= 0.8 and CV = 0.1 . As before, the data are gridded into bins by fitting the scattered data with a quadratic function. The velocity fluctuations are then obtained as the difference of the simulated particle velocity to the local fit value. This results in more accurate calculation of u p,rms in comparison with The normalized turbulence intensity ( Î = I∕U ∞ ) around the airfoil leading edge obtained from the numerical simulation is shown in Fig. 13 (left). For the conditions tested, I reaches about 3-4% of U ∞ around the leading edge. There is a good agreement between the simulated Î levels and that obtained from the proposed approximation (Fig. 13, middle), with the magnitude of the turbulence intensity error being well captured. Additionally, if a diameter coefficient of variation CV d of 10% d p -dispersion typical of jetting regimeis included in the simulations (Fig. 13, right), the results remain virtually unchanged. This supports the arguments given in the discussion of Sect. 4.1 that the influence of size dispersion on the slip velocity distribution is negligible in comparison with that of density dispersion.

Experimental procedure, image processing and uncertainty
The procedures developed in chapter 3 to estimate slip velocity and in chapter 4 to obtain velocity RMS from density dispersion are applied to evaluate the HFSB tracing fidelity in large-scale PIV measurements.

Set-up of experiments
The experiments are performed in the low-speed tunnel (LST) of the German-Dutch wind (DNW) tunnels, a closed-circuit tunnel with a closed test section of 3 m (height) × 2.25 m (width) cross-section, area contraction ratio of 9:1 and free stream turbulence level of approximately 0.03%. The 2D high-lift airfoil represents an outer wing section of the Fokker 100 aircraft (model 5-6, Reinders W 1994) of scale 1:4.96 and chord of 67.59 cm and was tested with retracted flap. The airfoil was installed vertically spanning the test section height. The measurements were taken at 15, 40 and 70 m/s free stream velocity and at three angles of attack α = {9°, 14°, 17°}. The planar two-component PIV system (Fig. 14) features two LaVision Imager sCMOS cameras (2560 × 2160 px 2 , 16 bit, 6.5 µm pixel pitch) equipped with 50 mm focal length objectives (lens aperture diameter of f/16 for HFSB and f/5.6 for DEHS). The cameras were installed on the top of the test section with their optical axis perpendicular to the laser sheet at a distance of about 1.5 m, yielding an optical magnification of 0.03, a digital imaging resolution of 0.2 mm/px and a combined FoV of 0.95 × 0.4 m 2 , covering the whole airfoil. A Quantel Evergreen 200 Nd:YAG laser (2 × 200 mJ/pulse at 15 Hz) was used for the particle illumination. The laser sheet thickness was 10 mm for HFSB and about 4 mm for DEHS. The laser power was set at 40% for the former and at 100% for the latter. As it can be observed, the imaging settings were at the limit in the case of DEHS, in terms of enhancing the optical signal, while the laser power had to be set to low power and the camera aperture fully closed for HFSB to avoid saturation. Thus, the volume achieved with HFSB could have been considerably larger. The acquisition and optical imaging conditions are summarized in Tables 2 and 3.

HFSB generation
The system of HFSB generation is composed of a fluid supply unit (FSU), fluid supply lines, flow resistors and bubble generators attached to the seeding rake. The in-house built FSU is composed of vessels and valves that can be operated remotely to pressurize and depressurize the fluid supply lines. Pressure flow controllers (coupled with mass flow meters) from Bronkhorst control the flow rates of helium, soap and air. Flow resistors guarantee equal mass flow distribution to the bubble generators. These are CNCmanufactured nozzles of 1 mm orifice diameter (Faleiros et al. 2019) designed by the Royal Netherlands Aerospace Centre (NLR). The bubble generator dimensions, working principle, regimes of generation and bubble properties have been recently studied by the authors (Faleiros et al. 2019). In the present experiments, the average volume flow rates per generator were 80 l/h of air, 9.5 l/h of helium and 9.5 ml/h of soap, yielding 30,000 bubbles/s (per generator) of nominal bubble density of 1.1 kg/m 3 ( ̂∼ 0.9 ) and mean diameter of 0.5 mm. The seeding rake, an array composed of six horizontal segments, where 42 bubble generators are distributed in intervals of 15 cm (Fig. 15), was installed in the settling chamber. Its influence on the flow was not quantified. However, the system was also in place during the DEHS measurements, with equivalent air flows through the bubble generator nozzles, to mitigate its influence on slip velocity measurements. The bubble system provides a seeded stream tube of about of 0.75 m (height) × 0.90 m (width) crosssection area, with an injection rate of approximately 1.3 million bubbles/s. The resultant stream tube of HFSB after the wind tunnel contraction is 0.25 m (height) × 0.30 m (width). The bubble concentration in the test section is 0.24 bubble/ cm 3 for U ∞ = 70 m/s, 0.42 bubble/cm 3 for U ∞ = 40 m/s and 1.1 bubbles/cm 3 for U ∞ = 15 m/s. This concentration level is not sufficient to perform spatial correlation analysis of the instantaneous images as it requires typically 5-10 particles per interrogation window, for which it would be necessary an interrogation volume of at least 5 cm 3 at the lowest speed, resulting in low spatial resolution.

Data quality
The raw images obtained with HFSB (Fig. 16, bottom) and DEHS (Fig. 16, top) elucidate the advantage of both particles, while DEHS measurements benefit from higher particle concentrations, HFSB are brighter, providing two orders of magnitude larger signal-to-noise ratios (SNR). The signal from DEHS particles is of the same order of magnitude of the noise level (300 counts), reaching an SNR of about 1.5. HFSB particles reach over 30,000 counts with SNR > 100. It is also noted a chordwise light intensity gradient (decreasing in intensity from left to right), observed in the DEHS raw images of both cameras. The light intensity scattered by particles upstream of the camera axis is higher than that scattered by particles downstream of it, due to a small (but significant) component of forward scattering (Raffel et al. 2018) from the former-there is a 20° angle between the light rays reaching the camera from the most upstream to the most downstream positions. The DEHS tracer signal is, therefore, lower near the downstream edges of the fields of view, compromising the DEHS data quality in these regions.

Image processing
The data obtained with DEHS were processed using the crosscorrelation algorithm from the LaVision software DaVis 10. The final interrogation window used is 48 × 48 pixels large (0.96 cm × 0.96 cm in physical space). With an overlap of 75% among adjacent interrogation windows, the vector spacing is 0.24 cm. Vectors whose absolute difference from the mean exceeded two standard deviations were excluded prior to obtaining statistics of the first and second moments of velocity. The data acquired with HFSB were processed using an inhouse algorithm developed with MATLAB. Particles were identified based on local maxima, their centre was obtained using a 3-point Gaussian fit (Raffel et al. 2018), which was used for pairing to the particle positions in the next frame according to the nearest neighbour criterion. The particle displacement (2 mm in the free stream) was sufficiently small compared to the average particle distance (about 20 mm at 40 m/s). The FoV was then gridded into square bins of 1.5 cm 2 for statistical analysis. The velocity moments were obtained in the same manner as described in chapters 3 and 4.

Uncertainty quantification
The uncertainty on the measured velocity with PIV is dominated by the uncertainty on the measured particle displacement, which is typically 0.1 px (Raffel et al. 2018). The particle displacements during the experiments were about 10 px, yielding a measurement error on the instantaneous velocity of 1%. Furthermore, assuming a Gaussian distribution for the random error, the expanded uncertainty estimate of the mean streamwise velocity u (analogously for the v component), respectively, with a 95% confidence level is given as (Benedict and Gould 1996;Sciacchitano and Wieneke 2016): where N(x, y) is the number of uncorrelated vectors at the grid ( x, y) location and u is the streamwise velocity standard deviation, which contains both the true velocity fluctuations and measurement errors. The uncertainty on the streamwise velocity component is almost everywhere about 0.1-0.2% of the local mean velocity (Fig. 17) with the exception of the separated region. Higher uncertainty at the leading edge region is observed for HFSB because of density dispersion. In addition, because the DEHS SNR is low, the uncertainty is also higher at the leading edge region, where the laser illumination is less intense, and a more restrictive mask was applied.

Uncertainty on the identification of HFSB particle centre
The HFSB glare point image is diffraction dominated (Appendix A1), which is well approximated by a Gaussian function (Adrian and Yao, 1985). The standard deviation of the Gaussian function that approximates the Airy disc is given as (Zhang et al. 2007): where is the laser light wavelength. This yields = 0.55 px for f # = 16 , = 532 nm and camera pixel pitch of 6.5 µm. The glare point size and the distance between the glare points are estimated to be 3.3 and 1.7 px, respectively (see Appendix A1). This means that the glare points partially merge, yielding a single particle image. The peak intensity of the glare point I peak reflected on the external surface of the bubble is on average 1.34 larger than that reflected on the internal part (Appendix 2) for HFSB, which leads to a local maxima at the centre of the brightest glare point. A sketch of the particle image considering these characteristics is shown in Fig. 18 (left). Measured particle images intensities within a window of 9 × 9 pixels around the local maxima (Fig. 18, right) support the proposed representation. Estimation of the particle centre at the local maxima leads to a position uncertainty of about d g ∕2 , where d g = d p ∕ √ 2 is the distance between the glare points. This is not concerning, as flow velocity changes should be negligible at the particle scale. The measured displacement is not affected by this effect, as long as the brightest glare point remains the same in both frames ( Fig. 19a and b). However, analysis of data from previous measurements (Faleiros et al. 2018), where the glare point images are distinguishable from each other, seems to indicate that this is only the case in 62% of the measured data (Appendix 2), while in 38% of the cases the uncertainty is of the instantaneous particle displacement is in the order of the glare point distance (1.7 px in this case). Taking this into account, the uncertainty on the mean particle displacement is estimated to be 0.7 px ( Δx = 0.62(0.1) + 0.38(1.7) ).
With the illumination at about 45 • angle with the flow direction (Fig. 14), it yields an uncertainty of about 0.5 px in the u and v velocity components. Fortunately, as shown in Fig. 19, cases c and d counteract each other, and the mean velocity is not affected significantly. This type of uncertainty mainly affects the velocity fluctuations and should be detected anywhere in the flow. As it is shown in Sect. 6.2, the turbulence intensity measured with HFSB in the free stream is comparable to that measured with DEHS. Therefore, this type of uncertainty has not affected the results significantly.

Slip velocity
The magnitude of the DEHS and HFSB velocities are shown in Fig. 20 for U ∞ = 40 m/s and α = 14°. The colour contours indicate a good agreement between both measurements. The overall agreement between the two measurements is verified by overlaying the HFSB velocity data with the DEHS velocity isolines (Fig. 21). Two measurement conditions are presented for the sake of conciseness: 40 m/s at α = 14° and 70 m/s at α = 17°. The results show a good agreement for the two velocity components in both flow conditions, although with larger deviations at 70 m/s and 17° incidence.  Overall, the slip velocity remains within 2% of U ∞ (Fig. 22), with the exception of the leading edge and separated regions. In most cases, the slip velocity in the leading edge region is essentially within 4 to 6%, reaching up to 12% of U ∞ near the surface. In the most extreme case (70 m/s and 17°), the slip velocity approaching the surface near the leading edge reaches up to 20% of U ∞ . It is noted that the slip velocity in the separated region reaches about 10% of U ∞ in some cases, in contradiction to previous results (Bourgoin et al. 2011;Faleiros et al. 2018). Close inspection of the data shows that the separation point shifted slightly in some cases between HFSB and DEHS, which could explain these differences. As the separation point has not been fixed through the use of a tripping device, a precise comparison of the slip velocity in the separated region is not possible.
The HFSB density is estimated through the least square optimization described in Sect. 3.4 (Table 4). A trend is observed which shows bubble density decreasing with wind tunnel speed. The bubbles are about 20-30% heavier than air at 15 m/s, within 10% difference from the neutral buoyancy condition at 40 m/s and about 20% lighter than air at 70 m/s. This systematic density variation with wind tunnel speed is not explained by the generation process. Changes in the tunnel total pressure do not affect the mass flow rate of helium, which is mass flow controlled. The pressure-controlled soap input of about 4 bars renders tunnel total pressure variations negligible (< 1% of the input pressure). Even though the viscosity of the soap film is sensitive to temperature changes, causing variations in the volume flow rate, those were monitored with a mass flow meter and counteracted by readjusting the helium mass flow for a constant helium-to-soap flow rate ratio ( Q He ∕Q soap ∼ 1000 ). An alternative explanation is that the HFSB density changed after generation, due to soap film evaporation and diffusion of helium and air through the soap film. Both the tunnel temperature (Table 4) and the bubble residence time (the time from generation until the measurement) influence these physical processes. Shrinking of HFSB, as a result of helium diffusion, and a colour shift from red to blue, attributed to soap film thinning, have both been observed by Huhn et al. (2017), while studying   (Fig. 23, bottom), based on the velocity fields measured with HFSB and Eqs. (20), (24) and (28). Those are compared to that of the measured slip velocity (Fig. 23, top). Despite measurement noise, there is sufficient agreement to support the validity of the density estimation procedure.

Velocity fluctuations
The normalized free stream turbulence intensity Î ∞ measured with HFSB (Table 5) are comparable to DEHS values at 15 m/s, while exceed in about 0.5% at 40 and 70 m/s. This difference is ascribed to the two different processing techniques, namely cross-correlation analysis and particle tracking. When conducted over two frames, the latter suffers higher uncertainty (Raffel et al. 2018) and can be more significant for HFSB (Sect. 5.5.1).  In general, a good qualitative agreement between DEHS and HFSB velocity fluctuations is observed (Fig. 24). The separated regions have similar topology, with HFSB measurements of turbulence intensity peak levels exceeding the reference by about 2%, with the exception of the 40 m/s and 14° incidence case, where DEHS measurements overestimate the turbulence intensity levels close to the trailing edge. The light intensity gradient of DEHS images (Fig. 16) results in higher turbulence intensity in the downstream edges, where the SNR is lower, and is especially noticeable in the measurements at 40 m/s. The turbulence levels measured with HFSB close to the leading edge are between 2 and 4% of U ∞ in most cases, except for = 17°, where it reaches 6%. As discussed in chapter 4, velocity fluctuations in the leading edge region are attributed to bubble density dispersion. It is noted, however, that no significant changes in the velocity fluctuations have been observed between the simulations at 14° and 17° incidence (not shown), posing the question of whether the leading edge fluctuations measured at 17° are indeed exclusively due to bubble density variation or a combination of the latter with existent flow fluctuations that are known to occur at the leading edge in the onset of stall (e.g. Benton and Visbal, 2018).
The density coefficient of variation CV is estimated via a least square optimization analogous to the one performed for the mean density estimation. The normalized density dispersion is mostly within 10 and 15% of the mean density (Table 6), with the exception of the data for 17° incidence, where is about 20% of the mean density. The density dispersion at 17° might be overestimated, however, as mentioned above, due to possible extra flow fluctuations that occur in the onset of stall. Overall, the estimated values are at least twice as large as those estimated in Sect. 4.1. The larger values might be explained by the extra uncertainties resultant from the simultaneous operation of a multi-nozzle system.
The Î values measured with HFSB in the leading edge region (Fig. 25, top) are compared to Î est (Fig. 25, bottom), obtained using Eq. (34) to (38) and the estimated CV (Table 6). Within the limitations of measurement and processing errors, and of statistical convergence, the estimated values show a good agreement with the experimental data, supporting the validity of the method proposed for estimation of HFSB density dispersion. Furthermore, the extra turbulence intensity caused by HFSB density dispersion is shown to be a localized effect, only occurring within regions of high flow frequency.

Summary and conclusions
This study has contributed to a more systematic assessment and prediction of the slip velocity of nearly neutrally buoyant tracers in typical conditions expected for large-scale PIV measurements in wind tunnels. It was shown that the slip velocity is governed by three non-dimensional parameters (as opposed to a single Stokes number), namely the ratio of particle to fluid density ̂ , the diameter-based Reynolds number Re d and the local frequency f (normalized with the viscous timescale d 2 p ∕ ). An analogy has been made between the particle motion in a sinusoid flow and that of a particle travelling around an object, which provided generalization to the analysis of particle slip velocity in external aerodynamics. Empirical relations have been obtained through numerical simulations, allowing the estimation of the slip velocity without the need of time-consuming computations. A direct comparison to other available methods in the literature has shown significant improvement in the accuracy of the estimations, emphasizing the importance of considering unsteady forces and drag correction for the slip velocity estimation. The empirical relations obtained may be used to establish the limits of  Fig. 7 for the amplitude modulation). A summary of the main equations used for the slip velocity analysis is given in Table 7. The tracing fidelity of nominally 10% lighter-than-air HFSB is also verified through PIV experiments in the flow around an airfoil of 70 cm chord in typical conditions of subsonic wind tunnel measurements in aeronautics at 15, 40 and 70 m/s. The HFSB slip velocity was shown to be overall below 2% of U ∞ , with the exception of the high-acceleration region around the leading edge, where in most cases it is contained below 5% of U ∞ . Only in the most extreme measurement conditions, as in the onset of stall at 70 m/s, that the slip velocity close to the airfoil surface has reached values in the order of 10% of U ∞ . In general, the velocity errors were localized and had minor effects on the overall measurement quality.
The method has also been extended to allow evaluation of bubble density as a function of the slip velocity and of bubble density dispersion, through measurements of velocity fluctuations. Comparison of the experimental results with those from the numerical simulations indicate that the HFSB density has changed post-generation-remaining within 10-30% of the neutral buoyancy condition-most likely because of variations of wind tunnel temperatures and bubble residence time, which affect the density through evaporation of the soap film and helium diffusion. In addition, the HFSB density coefficient of variation with respect to the mean density was quantified to be approximately 10%, resulting in measurement errors of turbulence intensity around the leading edge of up to 5% of U ∞ for most test conditions.

Appendix 1 HFSB glare point size, distance and merging
When soap bubbles floating in the air are illuminated, the reflected light rays form two glare points on the image sensor of the recording device, resultant from reflections occurring externally and internally of the bubble (Fig. 26). The distance between them can be used to calculate the bubble diameter, while the midpoint yields the particle centre. Since the refractive index of helium and air are approximately the same, and the shift of the light ray direction within the soap film is negligible (the film thickness of a neutrally buoyant HFSB of 0.5 mm diameter is about 80 nm, Eq. (33)), refraction can be neglected. Thus, the spherical bubble diameter can be obtained from geometric considerations: where d g is the distance between the bubble glare points and is the angle between the incoming light and the imaging direction. If the camera is positioned perpendicularly to the laser light sheet ( = 90 • ) , as in a typical planar PIV set-up, then eq. (A1.1) reduces to: The size of a glare point g can be obtained geometrically as a function of the bubble size and imaging conditions. Consider the case of a bubble being illuminated at a straight (A1.1)