Flow visualisation in swimming practice using small air bubbles

A non-harmful system to visualize the flow around an entire swimmer in a regular swimming pool is developed. Small air bubble tracers are injected through the bottom of the swimming pool in a prescribed measurement area. The motion of these bubbles, which will be largely induced by the swimmer’s motion, is captured by a camera array. The two-dimensional velocity field of the water at arbitrary planes of interest can be resolved using a refocussing method in combination with an optical visualisation method, based on particle image velocimetry, which is commonly used in fluid dynamics research. Using this technique, it is possible to visualize coherent flow structures produced during swimming; it is demonstrated here for the dolphin kick.


Introduction
New insights in propulsion mechanisms and resistance can be gained by studying the swimmer from a hydrodynamical point of view. The motion of the swimmer drives fluid motion, which is reflected in the wake of the swimmer. Studying the flow structures emerging in this wake can give insight into the forces both experienced and applied by the swimmer, and thus may be used as a diagnostic for the efficiency of the swimming technique. The strength, size and direction of motion of coherent vortex structures and the interaction of the swimmer with these structures are all expected to be related to the swimming efficiency.
Within the area of experimental fluid dynamics several techniques exist to visualize and quantify the velocity field of the flow. A pervasive technique is particle image velocimetry (PIV) [2,14], in which the flow is seeded with neutrally buoyant tracer particles. The motion of these particles is captured with one or more cameras. The velocity field can be evaluated with spatial correlations between subsequent frames resulting in a coarse-grained velocity field.
In the last decade it has become more common to apply flow visualisation techniques in advanced studies of swimming propulsion [8,9,12,[17][18][19], optionally in combination with force measurements, to observe the flow phenomena in the neighbourhood of the swimmer and to analyse the propulsive mechanisms. For example, the unsteady flow field around the hand of a swimmer during front crawl swimming has been evaluated in more detail by Matsuuchi [12]. The vortices generated during underwater undulatory swimming have been visualized by Hochstein et al. [8,9]. They suggested that vortex re-capturing may be used to enhance propulsion in human undulatory swimming. Takagi et al. [17][18][19] have clarified the propulsive mechanisms during front crawl swimming and the sculling technique by a combination of pressure measurements at the hand and visualisation of the unsteady flow field around (robotic) hands. They concluded that the unsteady flow phenomena are essential in generating high propulsive forces. For instance, high propulsive forces during sculling motion are possible due to vortex re-capturing.
So far, all studies concerning flow visualisation in human swimming have been performed in a laboratory setting. In this study, a non-harmful system to visualize the flow around a swimmer in a regular swimming pool is developed, rendering these techniques available for practical use. The application in practice is accompanied with several restrictions and requirements, which makes the development challenging.
Our aim is to measure the velocity field in planes that can be selected during processing of the acquired images. For this we use synthetic aperture PIV (SA-PIV) [5], which has been proven to be suitable for analysing unsteady 3D biological flows [11,13]. Instead, we only use the 3D information to select planes of interest. This reduces the processing time because one does not have to process the complete 3D field for PIV. Moreover, if just planes are considered instead of volumes a very tight calibration is not required in the SA-PIV technique. This puts less pressure on the accuracy of the 3D calibration [5], which cannot be guaranteed in the swimming pool.
Usually in PIV experiments, small neutrally buoyant particles (which are typically illuminated by a laser light sheet) are used to act as tracer particles for the flow. However, this is not allowed in the swimming pool. Alternatively, small air bubbles have been chosen as tracer particles.
In Sect. 2 of this paper the experimental setup and the methods to reveal the flow velocity field are presented. In Sect. 3 the proof of principle of the application of this technique in swimming practice is shown by means of visualizing the flow in the wake of a swimmer performing different styles of the dolphin kick.

Methods
We first discuss the ideas behind synthetic aperture refocussing in Sect. 2.1. Subsequently, we present the experimental setup in the swimming pool in Sect. 2.2, followed by an outline of the post-processing in Sect. 2.3. For more details about the methods and setup we refer to the dissertation corresponding to this work [20].

Synthetic aperture refocussing
The synthetic aperture refocussing technique relies on the fact that the apparent position of a certain object relative to another object is shifted when viewed from different positions. With this shift, known as parallax, depth information can be obtained. In our application this depth information is used to select planes of interest in the velocity field. A schematic representation of the working principle of the synthetic aperture refocussing technique is given in Fig. 1.
In synthetic aperture techniques, multiple cameras are used, all of them recording the same measurement volume. Since all cameras view the measurement volume from a slightly different position, the images of particles on different cameras are mutually shifted, as illustrated by the blue dot on plane Z 2 and the red dot in Z 1 in Fig. 1a. With the help of a calibration [20], the shifts ( d A and d B in Fig. 1b) related to each different camera position and each plane of interest throughout the measurement volume can be determined.
During preprocessing this information can be used to align the images of different cameras into one refocused image with a narrow depth of field on a chosen focal plane, as illustrated in Fig. 1c. When refocussing happens on plane Z 1 , the red dot will appear in focus with a high intensity, while particles elsewhere in the measurement volume, like the blue dot, are out of focus due to the parallax of the cameras and appear repeatedly with low intensity. This offers possibilities to filter out the plane of interest. In summary, while all bubbles are in focus on all cameras, "refocussing" on the chosen plane is done by shifting and stacking images.

Setup
The experimental setup consists of a bubble system to generate the bubble tracers and a camera system to capture the motion of the bubbles. The bubbles are illuminated by the ambient light (daylight and regular pool lighting). A black canvas on the opposite wall to the cameras enhances the appearance of the bubbles. A schematic plan view of the bubble and camera system is given in Fig. 2.  The position of the measurement volume was prescribed (the third lane relative to the side wall at Z = 0 ). The visualisation of the flow around an entire swimmer requires a large field of view (FOV) and thus a large measurement volume. The choice of cameras and lenses was attuned to the principles of the synthetic aperture refocussing technique [5]. Besides that the design rules for PIV were also considered in the design of the system.

Camera system
Six monochrome cameras (Sony: XCG-CG240) with a full frame resolution of 1920 × 1200 pixels at a frame rate of 41 Hz were placed in a hexagonal frame with 0.3 m spacing in the double side wall of the swimming pool (see Fig. 2). To obtain higher frame rates of 50 Hz, the height of the active region on the light sensor chip was cropped to 956 pixels. All cameras were connected with ethernet cables (bandwidth: 100 gigabit/s) to the computer and trigger box to synchronize. Each camera (with lens) was placed in an underwater casing and aligned to a central point in the measurement volume [ (X, Y, |Z|) = (7.5, 1.05, 5.95) m ] to have as much overlap in the FOV of the cameras as possible. All cameras have a 16 mm lens (KOWA: LM16HC F 1.4 f 16 mm). This camera-lens combination results in a depth of field of ∼ 2 m and a FOV of ∼ 3.1 × 1.5 m 2 in the measurement volume (at 5.95 m from the cameras). The minimum focal plane spacing Z using this configuration is 24 mm [5], which is sufficient when considering multiple thick slices within the measurement volume to explore with 2D-PIV.
In this setup one pixel corresponds to ∼ 1.6 mm . In a typical PIV experiment, interrogation windows of 32 × 32 pixels are used with a 50 % overlap region [2,14]. This would result in a spatial resolution of ∼ 25 mm , which should be sufficient to observe coherent vortex structures with sizes of ∼ 125 mm . Considering the PIV design rule for in-plane motion (displacements of 1 4 window size) velocities up to 0.65 m/s can be resolved with ease [2].

Bubble system
The air bubbles are generated by the bubble system, which was fully integrated in the bottom of the swimming pool (see Fig. 2). The bubble system covers an area between 5 and 10 m from the starting platform, so that the FOV of the cameras is in the center of the generated bubble curtains. The bubble system consists of five parallel PVC tubes with a length of 4.5 m, placed 0.25 m apart. Each tube has a series of small conically shaped holes with a smallest diameter of 0.2 mm and a separation distance of 0.02 m along its length. The airflow is supplied by a compressor (Leonardo 201) placed in the basement of the swimming pool. The airflow through the bubble system is regulated per PVC tube with a flowmeter (Kytola EK4A) to control the number of bubbles in the measurement volume. Five homogeneous bubble curtains are produced with an air flow rate of ∼ 0.01 m 3 /min per tube. The air bubbles have a diameter of approximately 4 mm (particle image diameter of 2.5 pixels). The bubbles have an average vertical spacing of 12 mm, assuming an estimated rising speed of ∼ 0.25 m/s [10]. Therefore, the seeding density in a thick slice above each tube corresponds to 10.6 bubbles per interrogation window, which is in accordance with the PIV guidelines [14].

Calibration
For the calibration, a black canvas of 3.0 × 1.5 m 2 with a pattern of white dots was traversed in the water, through 3.1 m the measurement volume above the bubble system, at specified distances ( |Z k | ) from the camera array ranging from 5450-6450 mm with increments of 125 mm. At each distance a record was made with the six cameras. The shifts in the camera views were established with the calibration for each camera at every |Z k | [5,20]. Mapping functions were obtained by second order polynomial fits through the set of pixel coordinates of the grid points and corresponding world coordinates. The polynomial coefficients were found to depend linearly on Z, despite the coarseness of the calibration with manual grid positioning. Given this linear nature, interpolation of the polynomial coefficients for different Z was straightforward, and the mapping functions of any desired XY-plane in the measurement volume could be determined. More technical details and quantitative information about the calibration, its performance and the application can be found in the thesis of van Houwelingen [20].

Post-processing
The quality of this refocussing technique highly depends on the intensity distributions throughout the recordings of different cameras. Processing the recording is therefore an essential step to effectuate the refocussing. Based on the ideas of previous studies on SA-PIV [5,11,13], several preprocessing operations have been applied on the recordings of single cameras and the refocussed image. All processing was performed in Matlab R2015b.
1. The background was removed by subtracting an average background recording of about 50 frames for each camera; 2. The intensity distribution was normalized to create uniform brightness across the images. This so-called flat field image was created by dividing by an average of bubble recordings made up of about 500 frames; 3. Optionally, when the illumination of the measurement volume turns out to be too low, an asymptotic weighting function can be applied to improve the contrast of the bubbles; 4. With the imwarp function in Matlab, the recordings were transformed to world coordinates (mm) using the mapping functions following from the calibration. Simultaneously, the recordings were cropped towards regions where the images of the different cameras overlap; 5. To achieve equal pixel intensity distributions across the images of different cameras, the histograms of the images from different cameras were equalized using the standard Matlab function histeq; 6. The contrast of the image was improved to enhance the visibility of the bubbles by the imsharpen function in Matlab, which uses a Gaussian low-pass filter [15]; 7. The refocused image was generated by summing the images of different cameras [5]; 8. To reveal the bubbles in the plane of interest, the refocused image was thresholded on pixel intensities corresponding to the mean +3 , with the standard deviation of the pixel intensity across the image. Usually, in thresholding, all pixels with an intensity lower than the threshold are replaced with zero [5]. Here the eight neighbouring pixels of the pixel exceeding the threshold were kept. Therefore, the contrast distribution of the bubbles is partially preserved.

PIV
The PIV analysis in this study was applied with PIVtec software (PIVview version 3.6.2). The (final) window size was chosen 32 × 32 pixels with a typical overlap of 50%. A fast Fourier transform (FFT) correlation was applied between the interrogation windows, where the highest frequency components of the resulting spectrum were removed with a Nyquist frequency filter. To reduce the noise, the correlation was repeated twice with slightly shifted correlation planes. Through multiplication, the different correlations were combined into a single signal and the random noise peaks are reduced [7]. A multi-grid interrogation method, with an initial window size of 128 × 128 pixels, was applied to improve the spatial resolution of the velocity field. Image shifting was allowed to shift (deform) the windows in between the different interrogation passes according to the displacement data of the previous pass. Sub-pixel shifts were obtained using a Gaussian pixel interpolation scheme. To obtain the velocity field with sub-pixel accuracy, the correlation peak was detected using a least-square Gaussian fit [1]. Within the PIVtec software, the raw velocity data from the PIV analysis are subjected to an outlier detection method based on a normalized median filter to find spurious vectors [23]. A bi-linear interpolation scheme is used to calculate the replacement vectors at the location of the outliers. A Gaussian-weighted interpolation is used when several neighbouring vectors are also outliers [3]. Due to the use of an interpolation scheme, the reliability of the quantitative data is not optimal, but at this stage for this purpose it was convenient for the proof of principle of the visualisation method.
Further post-processing and analysis were performed in Matlab. The rising speed of the bubbles can be assumed equal, because they are of similar size. Based on the PIV results, the mean rising speed of the bubbles was found to be 0.33 ± 0.03 m/s , which was determined by averaging all velocity vectors resulting from the PIV analysis on that plane. This mean rising speed was subtracted to obtain a better view of the flow velocity induced by the swimmer. Vorticity is a key quantity to understand propulsion. Vorticity can be computed from a measured velocity field through differentation, which can produce large errors. The presence of vortical structures can be clearly visualized by means of the vorticity component perpendicular to the plane of visualisation. This vorticity component z is defined as: with the velocity vector with v x and v y the velocity components in the x and y direction. Rather than calculating the vorticity by means of differentiating the velocity, the vorticity was obtained by integrating the velocity around an enclosed area following Stokes' theorem [14]: with X Y the area of a grid cell and i,j the circulation in point (i, j) estimated by: with U, V the velocity components coinciding with the contour around point (i, j). This method has the advantage of being less sensitive to errors compared to differentiation, because it uses velocity information from eight vectors [14]. The uncertainty of eight velocities contributes to the uncertainty in the vorticity, with the gain in accuracy depending on the correlation properties of the velocity field. Below, we will estimate this uncertainty and argue that the measured vorticity structures are well above the noise.

Results
An attempt has been made to visualize and analyse the flow around a swimmer performing kicks. One experienced swimmer volunteered to participate in this study and has given written informed consent. The ethical officer of the Eindhoven University of Technology approved the design of this study. The swimmer was instructed to swim different styles of kicks through the center of the measurement volume ( |Z| = 5950 mm ) at a depth of approximately 1 m, see video in the online supplementary material. The vortices produced during the kicks could be visualized using PIV. We expect to observe cross-sections of vortex rings (a patch of positive and negative vorticity) shed at the end of each kick [6,8,9,21,22]. Typical vorticity plots of these trials (slow, fast, low frequency + large amplitude and focus on powerful up-kick) are shown in Fig. 3a-d. Repetitive vortical structures are created after each downkick and were present in the vorticity plots. Due to selfinduced motion, the vortices move downward and appear under an incline.
In Fig. 4, we show on the velocity field and vorticity component perpendicular to this plane for a representative vortical structure. In particular, the jet produced by the vortex ring is clearly visible in the velocity vector field. The radius R of the vortex rings produced during the kicks is approximately 0.1 m (see Figs. 3 and 4). Because of the finite size, the vortex ring in Fig. 4 is not dominantly present in the other planes of interest ( |Z| = 5450, 5700, 6200, 6450 mm).
Since vorticity was determined by velocity differences over small distances, a point of concern is its uncertainty . This uncertainty depends on the random fluctuations u and v of the measured velocities, and their correlation properties in a way that has been documented by [16]. Their formula can be generalized readily to our averaged vorticity (Eq. 2), while u , v and their correlation can be measured from the background velocity field. This field was not constant, but fluctuates due to the fluctuations of the rising bubble velocity. With values u ≅ v ≅ 0.1 m/s in the experiments we find ≅ 1.3 s −1 . This is a crude estimate as u and v are not the same everywhere. Even with = 3 s −1 the vortical structures which we associate with kicks are way above the background.
An approximation of the impulse of the vortex ring can be estimated by: where is the circulation of the vortex ring (based on one of the patches in a 2D cross-section), R is the vortex ring radius measured from center to core and is the unit vector of the impulse in the axial direction normal to the plane of the vortex [13]. The impulse of the vortex ring in Fig. 4 is approximately 10.1 kg m/s ( ∼ 0.32 m 2 ∕s , R ∼ 0.1 m ). The accuracy of these values depends on the interpolation scheme used to obtain the velocity field. Nevertheless, it is a reasonable estimate of the generated impulse.

Discussion
The appearance of the vortical structures in Fig. 3 is confirmed by the fact that they exist for a long period of time (visible in a series of frames) and they appear on positions expected from the raw records (see supplementary video material for the raw video records). Moreover, the typical vorticity signature, consisting of a patch of clockwise and a patch of counter-clockwise rotation, resembles an intersection of a vortex ring. This is in agreement with previous observations in the literature [6,8,9,21,22]. Unfortunately, the applicability is not valid for all strokes yet. The vortices in the 'high frequency + small amplitude' and the 'focus on  powerful down-kick' trial did not appear in the PIV analysis, most probably because the vorticity is weaker, smaller or less coherent in those cases. It is peculiar that in most cases only the vortices after the down (extension) kick could be made visible. It seems that the vortices produced after the up-kick are weaker, smaller or less coherent, and are not captured by the PIV analysis. From the literature it is known that significantly more thrust is produced during the extension kick. During the up (flexion) kick, smaller and less coherent vortices are produced [6,22], which is in agreement with the results in this study and moreover explains the poorer visibility. This difference occurs due to a joint asymmetry within the human body, and a larger projected frontal area during the flexion kick [6,22]. Another reason can be found in the combination of a weak secondary circulation induced by the bubble curtains directed upward ( ∼ 0.04 to 0.10 m/s), and the self-induced motion of a vortex ring ( | | ∼ 0.18 m/s for R ∼ 0.1 m , ∼ 0.32 m 2 ∕s and = 0.05 m [4]). The selfinduced motion is also directed nearly straight upward for a possible vortex ring after the up-kick. Hence, these structures disappear from the field of view more rapidly.
Based on the results of this study, it can be concluded that the SA-PIV with bubbles in a regular swimming pool works properly to visualize four flow structures within a limited range of size and vorticity. The measurement environment unfortunately does not yet allow accurate, repeatable validation of experiments to define these limits precisely. However, some additional quantitative and qualitative experiments have been conducted with the system to test the performance of the system [20].
Regarding the accuracy in these experiments, most difficulties arise at areas where noise originates from an inhomogeneous bubble distribution, and not directly from additional bubble dynamics. The bubble curtains show some largescale collective motion, leading to empty areas in planes of interest. To improve the system, it is suggested to supply the bubbles by means of a large perforated plate (2D) rather than perforated tubes (1D), to achieve a more equal bubble distribution throughout the complete volume. The art of the refocussing technique lies in finding a good balance between filtering out-of-plane bubbles and retaining in-plane bubbles for the PIV analysis. The addition of homogeneous illumination of the measurement volume in such a way that the bubbles appear with a high intensity in all cameras can improve the results. Moreover, it would be interesting to test the implementation of pattern recognition for optimizing the algorithm, since the hexagon formation of the cameras is apparent in the refocussed images of the bubbles.
At this moment, just the flow in the wake behind the swimmer can be visualized, because the swimmer distorts the PIV analysis too much. In the future, it will be valuable to develop a dynamic mask suitable for these experiments to distinguish the flow features induced by the swimmer and the distortions originating from the swimmer's body (which currently cannot be left out of the PIV analysis) [13]. With that, the flow closer to the body of the swimmer can be captured as well, and the ability of resolving smaller coherent flow structures produced by the swimmer can be tested.

Conclusion
The objective of this study is to design, build and test a prototype system to visualize the flow around a swimmer in practice with the help of SA-PIV. The proof of principle is shown for a swimmer performing different styles of kicks. The expected vortex rings after each down-kick can be visualized in most cases. However, the quantitative results must be carefully interpreted at this stage, because inhomogeneity of the bubble distribution compromised particle image velocimetry. In future research, the system must be further optimized to increase the range of flow structures that can be visualized.