Laboratory modelling of the wind-wave interaction with modified PIV-method

Laboratory experiments on studying the structure of the turbulent air boundary layer over waves were carried out at the Wind-Wave Flume of the Large Thermostratified Tank of the Institute of Applied Physics, Russian Academy of Sciences (IAP RAS), in conditions modeling the near water boundary layer of the atmosphere under strong and hurricane winds and the equivalent wind velocities from 10 to 48 m/s at the standard height of 10 m. A modified technique of Particle Image Velocimetry (PIV) was used to obtain turbulent pulsation averaged velocity fields of the air flow over the water surface curved by a wave and average profiles of the wind velocity. The main modifications are: 1) the use of high-speed video recording (1000-10000 frames/sec) with continuous laser illumination helps to obtain ensemble of the velocity fields in all phases of the wavy surface for subsequent statistical processing; 2) the development and application of special algorithms for obtaining form of the curvilinear wavy surface of the images for the conditions of parasitic images of the particles and the droplets in the air side close to the surface; 3) adaptive cross-correlation image processing to finding the velocity fields on a curved grid, caused by wave boarder; 4) using Hilbert transform to detect the phase of the wave in which the measured velocity field for subsequent appropriate binning within procedure obtaining the average characteristics.


Introduction
The interaction between wind flow and surface waves is one of central problems of the study and parameterization of exchange processes in boundary layers of the atmosphere and ocean [1]. Of special interest is the case of steep and breaking waves forming at a strong wind. But carrying out field measurements of air flows, wave parameters for such conditions is very difficult problem (see. [2]). That is why laboratory modeling is used [3][4][5] to simulate conditions up to severe weather. The main difficulties in the experimental study of a turbulent air flow over a wavy water surface in laboratory conditions are connected with measuring wind characteristics near the water surface, especially in wave troughs, where one can expect the appearance of the most interesting features of this flow, such as shielding and separation of the flow. The technique Particle Image Velocimetry PIV is best fit for measuring the air flow in wave troughs and water flows close to the surface. In [6][7][8], the experience of using the PIV technique for measuring the air flow velocity over a wavy surface was presented. In [8], the structure of average velocity fields in the air flow and their perturbations initiated by waves, as well as the structure of turbulent stresses, were successfully studied.
However, those measurements were carried out at low wind velocities.
The subject of this work is studying the characteristics of air flows, waves, under conditions of breaking of waves with the formation of sprays near a wavy surface, in particular, in wave troughs due to high wind speeds.
Earlier, measurements under conditions of strong winds with an equivalent wind velocity U 10 > 25 m/s in the laboratory modeling of extreme meteorological conditions were carried out only using contact gages (Pitot tubes and hotwires/hotfilms) at a consider able distance from surface wave crests [3][4][5]9]. The measurement region was positioned above the layer of constant fluxes, where the velocity profile is logarithmic and its parameters (the friction velocity and roughness height z 0 ) could be determined directly using the profile method. Therefore, and z 0 have to be determined using the self-similarity property of the velocity profile of the flow in channels, as it was done in [4]. However, to refine the self-similar form of the velocity profile, it is necessary to measure the flow velocity as close to the water surface as possible [10].
The measurement can be implemented only using noncontact methods basing on visualisation, e.g., PIV. This work presents measurement results for the air flow velocity over waves in laboratory conditions modeling wide range of the winds (up to strong and hurricane conditions) with use of a specially modified PIV technique.

Experimental setup and measuring technique
Experiments were carried out on the Wind-Wave Flume of the Large Thermostratified Tank of IAP RAS (overview on Fig. 1 a). The length of the air flow channel with a cross section of 0.4 × 0.4 m is 10 m over the water surface. A detailed description of this setup and principles of design and control for the air flow in it were presented in [4,5]. The general scheme of the experiments is presented in Fig. 1 b. In addition to PIV methods, which were the main tool in the presented investigations, earlier approbated contact ways of measurements were used. In the working section of the channel, at a distance of 7 m from the entrance, average profiles of the air flow velocity were measured using a Pitot tube. The air flow in the channel was visualized using polyamide particles with an average diameter of 10 ȝm and density of 1.02 g/m3. The inertial time was 2 × 10 -4 s. The seeding device was similar to that used in [8] and positioned at the entrance of the channel at a distance of 6 m from the detection area. The test experiments demonstrated that the system brings no distortions to the wind flow in the region of measurements. To apply the PIV method (principal scheme shown on the Fig. 1 c), the motion of particles in the air in the region of measurements was illuminated by a laser sheet along the channel axis 8 m from the inlet. The laser sheet is formed by a cylindrical lens from a vertical laser beam (a continuous Nd YAG laser, 532 nm, 4 W). The width of the illuminated area was varied choosing the radius of the cylindrical lenses and their mutual arrangement. The motion of particles introduced into the air flow and water surface which were illuminated by the laser sheet was taken from the side by VideoSprint high speed camera positioned horizontally in a sealed box (see Fig. 1b). The camera was positioned horizontally and the level of its optic axis was above the level of the water surface by 8 cm. The focal plane was at a distance of 77 cm from the laser knife. The size of the taken area was from (66.4 -16 mm) × 256 mm (the horizontal size of a frame decrease with an increase in filming speed). A droplet removal system in the form of a metal tube blown with supplied compressed air was mounted from the internal side of the channel side wall through which the filming was performed. According to test experiments, the system introduced no distortions to the wind flow in the region of measurements.
The experiments were carried out at four values of air flow rates in the channel: 1.1, 1.6, 2.2, and 2.7 m 3 /s; as it is shown below, this corresponds to equivalent wind velocities at heights of 10 m, U 10 , 11, 20, 37, and 48 m/s, respectively. In two last cases, a strong breaking of waves with the formation of white caps and intense generation of sprays was observed (see photographs in Fig. 2). For each wind velocity, 30 implementations were taken, 3000-4500 frames in each implementation.
The filming speed was 1500, 3000, 5000, and 6000 fps; the exposure time was from 50 to 14 ȝs.

Fig. 1. a) Wind-Wave Flume of the Large Thermostratified
Tank of IAP RAS (air channel in the center) b) principal scheme of experimental setup in the working cross-section of the wind channel 1-laser, 2 -droplets blow-off system, 3high-speed camera, 4 -wind-wave channel. c) scheme of PIVmethod using.

Experimental data processing
Determining the shape of wave surface for each frame is necessary for finding the velocity field by crosscorrelation processing with the adaptive PIV algorithm over a curvilinear grid in the immediate vicinity from the water surface. Earlier, to determine the surface shape by images from a highspeed camera, a stepwise algorithm [11] was developed based on the Canny method [12]. That method operated well under conditions of weakly breaking waves. However, with an increase in the air flow velocity, a transition to a rather intense breaking with the formation of sprays and whitecaps was observed. In this order, a combined method of measuring the elevation of the water surface was used. In this method, optical measurements were completed by contact measurements using a wire wave recorder mounted on the channel axis near the edge of the laser sheet. The records of the level elevation and high speed camera were synchronized. The final form is a combination of data obtained by contact and noncontact ways; with an increase in wind velocity, the role of contact measurements increased up to a full substitution of optical measurements for the case of flow rate of 2.7 m 3 /s (see Fig. 3). After finding the surface shape, velocity fields were calculated by the cross-correlation method over a curvilinear grid taking into account the current shape of the surface [8]. We used a modified PIV processing method based on adaptive finding for the shift of the cross-correlation function (from now on, briefly CCF) peak for rectangular elements of the image for two sequential frames. The space on the grid of the velocity field calculation was 3.2 mm. To increase the algorithm accuracy due to a decrease in the search window size, the processing was performed in two stages. At the first stage, using a window with dimensions of 128 × 64 px, the profile of the average horizontal shift of particles for each wind velocity was found.
At the second stage, the comparison window was shifted to the magnitude of the shift found at the first stage. The size of the search window could be less than the final shift of particles, which increased the spatial resolution of the method, accelerated the processing, and made it possible to obtain more data. The CCF maximum was searched by the adaptive method during two iterations by analogy with [8]. At the first passage, the shift over the window with a larger size (32 × 32 px or 6.4 × 6.4 mm) was approximately determined; then, an after search was performed with allowance for the calculated shift for a window with a lesser size (8 × 8 px or 1.6 × 1.6 mm). To increase the accuracy at the last step of cross correlation, an algorithm of subpixel approximation of the peak by a Gaussian twodimensional profile (see [8]) was applied; this allows one to take into account CCF values at points near the maximum. Data that do not satisfy the quality criteria were eliminated from further processing. In particular, windows with an insufficient quantity of particles were eliminated from the processing. For each image, points corresponding to high values of the brightness gradient were found using the Sobel operator. These points correspond to edges of particle images in the frames. For each window, the number of similar boundary pixels was calculated. The threshold value of the brightness gradient in the process of searching the edges was matched so that no boundary pixels fell in regions without particles. The cross correlation was performed only with windows in which at least one boundary pixel was present. In this process, if the number of boundary pixels considerably exceeded the average number, flares, spindrifts, or a foam line were most often present in part of the image corresponding to this region; as a consequence, the shift calculated for it was not taken into account in the process of further averaging. Regions for which the value of the CCF maximum considerably differs from the typical value in the implementation were also not taken into account, because such regions most often corresponded to the case of a particle escaping from the laser sheet during the interval between frames.
At strong winds, the inhomogeneity of seeding increased, which resulted in considerable discontinuities (the absence of data) in time dependences of velocities at fixed horizontal coordinates. Despite measures for the removal of droplets, optical distortions increased with an increase in the wind velocity. This resulted in decreasing of accuracy.
To correctly calculate average velocity fields and profiles, a method based on phasing the measured wind velocities at different levels from a wavy surface was used. The phase of the wave over which the point of velocity measurement is positioned was determined using the Hilbert transform over time realizations of water surface elevations presented in several equidistant vertical cross sections for each frame. There were from three to seven such cross sections in a frame (depending on the flow rate).The position of the surface for each vertical column of the velocity field was calculated using linear interpolation. Thus, each column has its own dependence of the surface position on time.
Wave number spectra of waving have a pronounced peak, and deviation of the water surface from the horizontal is determined mainly by the fundamental wave harmonic. For this reason, averaging over the phase of the fundamental wave harmonic was performed. Frequency filtering of time realizations obtained in the aforementioned way was performed over a band with a width of 2 Hz near the peak frequency determined for each wind velocity by time realizations from the wave recorder. Such filtered realizations with a removed constant component fit well for determining the wave phase by use of the Hilbert transform. Using the economical algorithm of the fast Fourier transform, the Hilbert transform permits one to determine the amplitude and phase of surface elevation as a function of time.
To obtain velocity fields averaged over turbulent fluctuations, conventional averaging at a fixed phase was performed. To reduce errors connected with the insufficient number of measurements, the data were binned in phase intervals of 18°, which yields 20 different binns of the phase.
The averaging was performed in two ways: (i) Averaging over fixed horizons (Fig. 4a). The velocity fields were divided into horizons with a step of 32 px (6.4 mm of the real height). For each horizon and each phase, velocities were accumulated from the whole realization and then averaged. (ii) Averaging over a fixed distance from the surface (Fig. 4b). For each frame in each column of the grid, a cell with the same number was chosen; in the curvilinear coordinate system used, this meant a fixed distance of the point from the surface. The data from corresponding similar cells over the height and phase were accumulated from the whole realization and then averaged.

Results
Pictures of turbulentͲpulsation averaged velocity fields of the air flow were obtained for both ways of averaging. In curvilinear coordinates, by analogy with [8], the vertical coordinate is counted from the position of the water surface at each instant of time. The position of the average surface level was taken as zero in rectangular coordinates. The horizontal coordinate is the wave phase for a given point ĳ recalculated using the value of the wavelength Ȝ determined by the dispersion relation for waves on deep water for a frequency corresponding to the peak in the center of the spectrum for each value of the wind velocity: 2 x

M O S
Examples of flow fields for both averaging are presented in Fig. 5.
The number of points satisfying the criteria of data quality is different for the same distance from the surface in different phases. This difference is especially seen near the surface, where the decrease in the quantity of data satisfying the quality criterion from the leeward side of the wave crest (negative values of the phase in Fig. 6) caused by decrease in the number of tracing particles in this region at the instant of shooting due to the shielding of the wind flow by the wave crest. Note that, for the cases of two high wind velocities (U 10 of 37 and 48 m/s), an inverse picture is observed: the number of points satisfying the criteria of data quality becomes larger from the leeward side of the crest. This can be connected with the appearance of sprays, which are intensely generated at velocities U 10 > 25 m/s. The droplets concentration from the windward side of the wave crest is considerably higher than from the leeward side, because the droplets are generated mainly near the top of the wave and whirled away by the wind. The crossͲ correlation algorithm for cells in this region shows a shift of droplets, not tracing particles, which are relatively not numerous there.  Since the droplet velocity is lower than the wind velocity, this leads to an erroneous underestimation of the air flow velocity near the surface. The different number of particles in different phases makes it necessary to apply conventional averaging over the phase to obtain correct results. Velocity values obtained using the Pitot tube always turned to be lower than those obtained in PIV measurements. Note that the PIV technique is a direct method of determining the air flow velocity; in contrast to it, the Pitot tube is an indirect method of estimating the flow velocity by pressure and the underestimation of the flow velocity by the Pitot tube probably indicates the presence of a systematic error.  Both the velocity profiles have a similar shape and close values; they agree with each other well in the region above wave crests. The main difference between them is that averaging over curvilinear coordinates can yield air flow velocities in wave troughs, while averaging over horizons can yield the vertical profile of the average wind velocity only above wave crests. In connection with this, the aerodynamic drag factor of the water surface was determined by dependences obtained by averaging over curvilinear coordinates.

Conclusion
Laboratory experiments on studying the structure of the turbulent air boundary layer over waves were carried out at the WindͲWave Channel, IAP RAS, under conditions modeling the near water boundary layer of the atmosphere under strong and hurricane winds. Air flow in the channel with a square cross section of 0.16 m2 was created by fan. The air flow rate took values of 1.1, 1.6, 2.2, and 2.7 m 3 /s, which corresponded to estimate values of equivalent wind velocities from 10 to 48 m/s at the standard height of 10 m.
A modified technique of Particle Image Velocimetry (PIV) was used to obtain turbulent pulsation averaged velocity fields of the air flow over the water surface curved by a wave. We emphasize that such remote methods permit one to obtain the air flow velocity field also below wave crests in wave troughs. Using wave phase averaging at fixed distances from the water surface, average profiles of the air flow velocity were obtained in curvilinear coordinates following the wave. Note that measurements could be reliably performed only using the PIV technique. Using contact measurements (Pitot tube), profiles of the air flow velocity over wave crests were measured in Cartesian coordinates. To compare measurement results obtained by two techniques, measurement data obtained by the PIV technique were also expressed in Cartesian coordinates. Both experimental methods yielded close results at distances of more than 10 mm from wave crests. But introducing the Pitot tube with a diameter of about 10 mm into the flow could resulted in considerable distortions of the flow close to the water surface. Parameters of the air flow (friction velocity and drag coefficient of the surface) could be further determined by extrapolating the logarithmic part of the velocity profile, and their dependencies on the wind speed could be obtained.
This work was supported by the Russian Foundation of Basic Research (No. 15-35-20953 17-05-00958, 14-05-91767 AF_ɚ, 16-55-52022), providing experiments was partial supported by Russian Science Foundation (Agreement No. 14-17-00667) and processing of the data was partially supported by Russian Science Foundation (Agreement No. 15-17-20009). Also work is supported by international project Air-Sea Interaction under Stormy and Hurricane Conditions: Physical Models and Applications to Remote Sensing (ASIST) ASIST within the FP7 program and of the Federal Target Program "Research and development in priority areas of Russian