Effect of Manufacturing Inaccuracies on the Wake Past Asymmetric Airfoil by PIV

Abstract

The effect of manufacturing geometry deviations on the flow past a NACA 64(3)-618 asymmetric airfoil is studied. This airfoil is 3D printed according to the coordinates from a public database. An optical high-precision 3D scanner, GOM Atos, measures the difference from the idealized model. Based on this difference, another model is prepared with a physical output closer to the ideal model. The velocity in the near wake (0–0.4 chord) is measured by using the Particle Image Velocimetry (PIV) technique. This work compares the wakes past three airfoil realizations, which differ in their similarity to the original design (none of the realizations is identical to the original design). The chord-based Reynolds number ranges from $1.6\times {10}^4$ to $1.6\times {10}^5$. The ensemble average velocity is used for the determination of the wake width and for the rough estimation of the drag coefficient. The lift coefficient is measured directly by using force balance. We discuss the origin of turbulent kinetic energy in terms of anisotropy (at least in 2D) and the length-scales of fluctuations across the wake. The spatial power spectral density is shown. The autocorrelation function of the cross-stream velocity detects the regime of the von Karmán vortex street at lower velocities.

1. Introduction

Fluid flow is a highly non-linear problem that still lacks a reasonable solution. One of the general effects of non-linearity is the unpredictable response to even small perturbations or changes of the boundary conditions. The scale of possible responses to a small geometry perturbation ranges from almost zero effect to a linear response and up to a complex change of flow state. The famous butterfly effect represents this behavior with the example of a small butterfly that can alter the evolution path of a turbulent system [1]. As there is a large number of such micro-events, the evolution path of the entire system is unpredictable. However, the statistical properties can be predicted quite reasonably. This feature is used in modern computational fluid dynamics, which does not solve the non-linear Navier–Stokes equations on a fine mesh with the resolution of the Kolmogorov length-scale, but it solves only much larger cells with a direct link to the geometry of the boundary conditions under the assumption that the behavior at smaller scales follows some of the turbulence models.

The relatively low-cost computational fluid dynamics are used in thorough experiments in the area of industrial design and optimization. However, the computational methods need validation and verification [2], which are mainly based on the comparison with an appropriate experiment, which does not need to fit in all parameters, but at least the main geometrical and fluid properties might be met. This is where the problem with geometry becomes important—the object used in the experimental study is different from the desired design used in the computational approach. As the quality needs to increase, even smaller deviations between numerical and experimental results can be accepted. As the experiment has a semblance of truth, mathematicians have devoted a great deal of effort to developing even better and more accurate models in order to fit the experimental data measured with boundary conditions that are known only approximately.

In this study, we used a high-quality, contactless 3D optical scanner to measure the real geometry of boundary conditions (of course, with some uncertainty connected with each measurement technique) and compare the flow field past an airfoil produced at three levels of accuracy—the first variant represents the naive approach, and we sent to the 3D printer only the coordinates obtained from a public database [3]. If we had not used a 3D optical scanner, our study would have ended by declaring “we measured flow in the near wake past airfoil…”. However, the 3D scan discovered discrepancies from the ideal geometry, and therefore we attempted to prevent these by manipulating the geometry in such a way that the outcome could be brought closer to the ideal case. This step was repeated creating two more variants that were increasingly close to the ideal geometry. In any case, the best product was still different from the ideal geometry used in the processes of numerical optimization.

There are only very few other studies about manufacturing inaccuracies, and they are numerical only. Our study is experimental and compares models that have been actually created. A disadvantage of this work is that the ideal geometry does not exist in reality; hence, we can compare only different levels of manufacturing attention.

Some Recent Literature Concerning This Problem

Moreno-Oliva and coworkers [4] used the Laser Triangulation Technique to measure the real profile of a NACA0012 and FX 61-137 airfoil; one realization of the NACA0012 has been commercially made of metal, while the second has been 3D printed by the authors from PLA, and the profile FX 61-137 was measured at the real wind turbine blade at several heights. They used the measured profiles as an input for the Xfoil software [5] together with the theoretical profiles. They report a significant decrease of performance of the real wind turbine profiles when compared with the ideal geometry.

Ravikovich et al. [6] performed an optimization of fan blades with real geometry deviations from airfoil geometry. Their algorithm simultaneously solves mechanical integrity problems and the aerodynamic stability of the flow. The geometry data measured under cold conditions are transformed via numerical simulation into hot conditions, which allows the estimation of the influence of the measured geometry on the final fan parameters.

The effect of even small geometry modifications in turbine stator blades was studied by Klimko and Okresa [7]. Modifications in the stator wheel affect the reaction and efficiency of the entire stage. Data like those mentioned above are used in the numerical optimization [8] of the developed new stage.

Winstroth and Seume [9] studied numerically the effect of several artificial geometry deviations on four wind turbine airfoils. The geometry deviations cover mold tilt towards the leading or trailing edge, step changes at certain positions, sine waves placed at several positions on both pressure and suction sides or at the leading edge, and the thickening of the profile—in total, 40 modifications of four different airfoil profiles are used in the studied wind turbine. They found that the worst deviation (in terms of energy production) is the tilt towards the leading edge; they found also some cases with a positive impact.

The structure of our article is as follows: first, we describe the methods including 3D printing, 3D optical scanning, and Particle Image Velocimetry used to measure 2D velocities within a square area past the trailing edge. This paper focuses on data at zero angle of attack; some data at an angle of attack of $\alpha ={10}^{\circ }$ are in the Appendix A. The ensemble average stream-wise velocity, the path of the wake centerline, and the wake width based on the average stream-wise velocity are shown. The lift coefficient is impossible to determine on the basis of PIV data; therefore, a simple balance is introduced and the lift at zero angle of attack and at several Reynolds numbers is measured. The turbulent kinetic energy is explored to discuss if it is produced by fluctuations in the stream-wise or cross-stream velocity component. The length-scale of fluctuations is explored by a unique method developed at the University of West Bohemia in Pilsen. The autocorrelation function reveals the periodicity at low velocities and the decay of the integral length-scale at larger ones.

2. Materials and Methods

2.1. Three-Dimensional Printer

We used a commercially available 3D printer Prusa Mk 2.5. It uses the FDM (Fused Deposition Modeling) printing technology [10]. The used material is thermoplastic polymerized lactic acid (PLA) [11]. The nozzle temperature was $230{}^{\circ }\mathrm{C}$, and the height of a single layer was $0.15\mathrm{mm}$. The 3D printed airfoils are used mainly for testing and development purposes. At the current stage, 3D printing technology is not suitable for industrial manufacturing. However, it is important to know the impact of the typical 3D printing shortcomings.

2.2. Three-Dimensional Scanner

The optical 3D scanning method FPS (Fringe Projection Scanning [12]) is based on projecting light periodic stripes onto the surface [12,13], and the variations of their phase are used to estimate the variations of height across the measured surface. The commercial 3D optical scanner GOM Atos Core [14] is used in our study. Mendřický [15] measured the accuracy of a such scanner dependent on the calibration and the optical matteness of the measuring surface of standard spheres. He found that the accuracy is on the order of micrometers, and a matte coating adds around ${10}^{-5}\mathrm{m}$, while simply an old calibration can change the measured results systematically by around $2\times {10}^{-5}\mathrm{m}$. Li et al. [13] found that ambient light has only a small effect on the results, while Vagovský et al. [16] states that even the size of the scanned object influences the accuracy due to the focus blur. The measuring volume of our scanner is $300{\mathrm{cm}}^3$; thus, the airfoil of height $12.5\mathrm{cm}$ is small for such a volume, but still, our object size is much closer to the measured volume than in the case explored by Vagovský [16]. We did not use the matte coating because we used white PLA material; on the other hand, we had to paint the material black before the PIV measurement, which uses intense laser light, and we did not measure the thickness added by this painting.

2.3. Model

The coordinates of the NACA 64-618 airfoil were obtained from the public database Airfoiltools.com [3]. This airfoil is frequently used at the tips of the wind turbines [17] and has a reasonable lift even at high angles of attack (which is valid only for larger Reynolds numbers, as shown later). The trailing edge of this airfoil is sharp, and it points towards the high pressure side; see the black line in Figure 1. First, this airfoil was 3D printed only according to the point coordinates obtained from the public database [3]; this first naive realization is denoted A. Then, a high-precision 3D optical scanner GOM Atos [18] is used to measure the deviation from the ideal geometry; see Figure 1. There exists a huge systematic printing error—the missing trailing edge; see Figure 1c. This is caused by the finite thickness of the extrusion line, which does not fit into the very thin walls, and the program for generating the g-code does not allow the extruder to move into that location. The trailing edge is shorter by almost $3\mathrm{mm}$, which is $3.7\%$ of the chord length.

A next step was to make the trailing edge thicker in order to fill it with at least some material. This was performed by the Minkowski sum [19] of the airfoil profile with a circle of radius $0.2\mathrm{mm}$. This method led to an overly large thickening of the entire model, but the trailing edge was still shorter than it should be. This blind method is not reported here. The next trial, denoted B, was performed by artificial manipulation with the profile points at the trailing edge combined with the Minkowski sum with a circle of radius $0.1\mathrm{mm}$. Now, the trailing edge was shorter by around $0.2\mathrm{mm}$ only. On the other hand, the trailing edge displayed a large radius of around $0.3\mathrm{mm}$.

The best result was obtained by using the Minkowski sum with circle of $0.1\mathrm{mm}$ radius, adding new points into the trailing edge area and expanding it over the ideal edge. After printing, the model was polished manually in order to sharpen the trailing edge. As a side effect, the surface of this last realization, denoted C, was smoother than the previous ones. The problem with manual polishing is that (i) it is a non-repeatable procedure, and (ii) it is too expensive for possible industrial use.

2.4. Wind Tunnel

The University of West Bohemia in Pilsen has multiple wind tunnels [20,21]. However, due to the laser safety requirements, we are limited for the use of the PIV method to the “small” wind tunnel [22] with a transparent test-section of length $400\mathrm{mm}$ and cross-sectional size of $125\times 125\mathrm{mm}$. It is an open low-speed wind tunnel driven by a radial fan at the tunnel inlet. The smallest stable velocity used in this study is $U_{\mathrm{ref}}^{\min }=3.1\mathrm{m}/\mathrm{s}$; the largest used reference velocity is $30.6\mathrm{m}/\mathrm{s}$. The working fluid is atmospheric air with a temperature of around $20{}^{\circ }\mathrm{C}$, and the actual humidity was not measured.

The Reynolds number used in this study was based on the above-mentioned reference velocity $U_{\mathrm{ref}}$, which fits the incoming velocity in the case of an empty test section of the wind tunnel. As the length-scale for the Reynolds number, the planned chord length of the airfoil was used. It was a constant: $c=80\mathrm{mm}$. Note, the real chord length slightly differed among realizations A, B, and C (see

Table 1. Dimensions of the studied airfoil realizations compared with the ideal case, denoted CAD in the first column. The blockage ratio of the wind tunnel is the ratio of the blocked projection area to the total cross-section area ($125\times 125\mathrm{mm}$); it is calculated with the caliper approach—i.e., the distance of the plane touching the suction side to the second plane touching the pressure side, when their normals are locked to the y-axis.

	Unit	CAD	A	B	C
Chord length c	mm	$80.00$	$77.01$	$79.83$	$80.08$
Profile thickness—maximum inscribed circle	mm	$14.36$	$14.34$	$14.55$	$14.47$
Profile thickness—caliper along y-axis	mm	$14.44$	$14.41$	$14.63$	$14.55$
Blockage ratio	%	$11.55$	$11.53$	$11.70$	$11.64$

for more details). The range of Reynolds numbers used in this study spanned one order of magnitude from $1.63\times {10}^4$ (corresponding to $U_{\mathrm{ref}}=3.1\mathrm{m}/\mathrm{s}$) up to $1.63\times {10}^5$, where $U_{\mathrm{ref}}=30.6\mathrm{m}/\mathrm{s}$.

2.5. Particle Image Velocimetry

The method of Particle Image Velocimetry [23] visualizes the wake structure in a small area just behind the airfoil; see Figure 2a. This method is based on the optical observation of the movement of small particles carried by the studied fluid. In our case, the particles are droplets of Safex with a diameter in the order of micrometers. When the droplets are small enough, they follow the fluid motion because the inertial forces scale with particle volume—i.e., as $\sim r^3$—while the viscous forces unifying the motion of particles and fluid scale with the particle surface—i.e., as $\sim r^2$. These particles are illuminated in a single plane by using a solid-state pulse laser, New Wave Solo, using the second harmonics of wavelength $\lambda =532\mathrm{nm}$ (corresponding to a green color). The single shot energy is $0.5\mathrm{J}$, and its duration is $5\mathrm{ns}$; thus, the peak power is $100\mathrm{MW}$. The laser produces a pair of pulses separated by tthe ime-interval ${10}^{-5}-{10}^{-4}\mathrm{s}$ dependent on the flow velocity. The laser beam is defocused by using cylindrical optics in order to illuminate the plane. The camera is oriented perpendicularly to this plane and focused into it. The Mk II Flow Sense CCD camera has resolution of $2048\times 2048$ pixels; this camera is able to resolve both pulses of the laser into a pair of expositions. The obtained pairs of photographs of particles are processed by using the Dantec Dynamic Studio software (details for readers familiar with this software: (1) first five frames are removed from each movie, (2) the minimum of image ensemble is calculated by using Image Min/Max function, (3) this minimum is subtracted from each image in the ensemble by using the Image Arithmetic function, (4) the velocity vectors are calculated by using the Adaptive PIV function with parameters: grid step size 32, minimum IA size 32, maximum IA size 32, validation based on universal outlier detection over neighborhood $5\times 5$, normalization $0.1$, acceptance limit $2.0$; the algorithm adapts to the particle density with a particle detection limit of 5 to the desired number of 10 particles per IA; it adapts the IA shape to velocity gradients with $\left|Ux\right|,\left|Vx\right|,\left|Uy\right|,\left|Vy\right|\le 0.1$ and $\sqrt{U{x}^2+V{x}^2+U{y}^2+V{y}^2}\le 0.2$, the convergence limit is $0.01$ pixel and the maximum number of iterations is 10. (6) The last step within the Dantec software is exporting via the Numeric Export function). The vector fields contain $64\times 64$ vectors, where a single grid point covers an area of $0.53\times 0.53{\mathrm{mm}}^2$ corresponding to $6.65\times {10}^{-3}$ times the chord length c. The analysis continues by using our custom-made software: one grid point is removed from each side, obtaining vectors with a resolution of $62\times 62$, because the vectors at the edge of the field of view are not trustworthy, as there is a relevant number of particles leaving or entering the field of view during the time between laser pulses. The second step is to filter out the snapshots, whose small-scale energy content is significantly greater than the average of the ensemble. This procedure serves to remove rough errors, whose smallest-scale energy content is typically large. This method is described in more detail in our previous article [24]. In fact, the current data have reasonable quality; therefore, this procedure removes at least one snapshot only for 40% of datasets. An example of the obtained instantaneous velocity field is displayed in Figure 3.

2.6. Lift Force Measurement

Wind-tunnel balance is designed to measure the lift force L on the airfoil. The basis of the balance is a resiliently suspended movable table. Elastic suspension is created with four parallel flexible elements. The elastic suspension has the lowest stiffness in the direction of the force L, and in all other directions, the stiffness is an order of magnitude higher. The displacement of the table is measured by an eddy-current sensor that responds to the distance to the steel target. Unwanted vibrations are dampened by eddy-currents in the copper sheet as it moves near to the neodymium super magnet. The output of the sensor is calibrated by a known force created by using a calibration weight in the lift direction.

3. Results and Discussion

3.1. Zero Angle of Attack—Average Velocities

The average stream-wise velocity component for four different Reynolds numbers is displayed in Figure 4. In that figure, the spatial distribution for each case (velocity and airfoil realization) is plotted as a colormap. The direct comparison of isotachs (contour lines of velocity) for each product is plotted in a single panel in order to see the differences between those products. Figure 5 shows the profile of the average stream-wise velocity at distance $10\mathrm{mm}=1/8c$ past the trailing edge. This plot contains less information than the previous one, but this information is organized in a different manner in order to show the differences better.

A first look at Figure 4 and Figure 5 reveals that the largest difference occurs at a Reynolds number of $4.08\times {10}^4$ (denoted 40.8 k, where k represents $\times {10}^3$ as a standard multiplicator of the SI system). This Reynolds number is close to the critical Reynolds number of transition between a wake dominated by the von-Kármán vortex street [25] and the wake dominated by continuous boundary layers developed along the airfoil, as will be apparent from the next figures analyzing turbulent kinetic energy (mainly the ratio $\sigma \left[v\right]/\sigma \left[u\right]$ in Figure 6). As was mentioned above, variant A is slightly smaller than B and C, and both are expanded by the Minkowski sum with a circle of radius $100\mathrm{mum}$. Therefore, it is natural to expect that the real Reynolds number for those variants is slightly larger than that for variant A; thus, the transition occurs at slightly smaller velocities.

The next largest discrepancy lies in the position of the wake past variant A. It is most apparent at a Reynolds number of $8.17\times {10}^4$, although it applies at others as well. This shift is caused by the missing trailing edge of the variant A; see Figure 1. As the trailing edge is slanted towards the pressure side (direction down in the figures), the entire airfoil ends at the top (in figures), and thus the entire wake is shifted up. It is even more apparent in the velocity profiles in Figure 5 and Figure 7. One may ask why this shift is less apparent at the lower Reynolds number—we think that this is due to the different mechanism of wake creation: the von Kármán vortex street is formed due to Kelvin–Helmholtz instability of the shear layer between the fluid, which passes far from the obstacle, and between the fluid, which is slowed by the presence of the obstacle. However, this difference starts at the leading edge, which is much more similar for different airfoil variants; see Figure 1b.

3.2. Balance Measurement of Lift Coefficient

According to the momentum balance principle, the lift force acting in one direction has to deflect the fluid motion in the opposite direction. However, the true integration of the deflected momentum requires the knowledge of the flow field around the entire airfoil, which is not available. In order to obtain some information about the lift, which is in any case the crucial feature of any airfoil, we measured this by using a force balance equipped with four flexible metal elements in order to allow deformation in a single direction only; see Figure 8. The resulting lift coefficient is shown in Figure 9 for an angle of attack equal to $0^{\circ }$ and for Reynolds numbers in the same range as used for the PIV measurement. This figure shows a reasonable match with the prediction by the Xfoil algorithm [5,26,27] at a Reynolds number of ${10}^5$, but the values are greatly different at lower and higher velocities. Differences among airfoil realizations show that the A variant exhibits the lowest lift, as expected. The geometrically best airfoil C has the best aerodynamic performance. Interestingly, the A airfoil has larger negative values of lift at low velocities. The large errors reported here are caused by the unwanted oscillations at higher velocities (despite the damping by using highly-conductive cooper and a strong magnet). At smaller values, a large error is caused simply by the uncertainty of the displacement sensor and of the analog–digital converter, because the forces are two orders of magnitude smaller than the range adapted to fit the expected forces at largest velocities.

3.3. Rough Estimation of Drag Coefficient

The time-average drag force $〈F_D〉$ can be roughly estimated according to the principle of momentum balance:

$$\frac{〈F_D〉}{\rho }={\int }_{\mathrm{wake}}\left(U_{\infty }-〈u〉\right)〈u〉\mathrm{d}S-{\int }_{\mathrm{wake}}{\sigma }^2\left[u\right]\mathrm{d}S+{\int }_{\mathrm{wake}}\left(p_{\infty }-〈p〉\right)\mathrm{d}S,$$

(1)

where $U_{\infty }$ is the velocity in the unperturbed state—see the work of Terra et al. [28,29] for more details; $\mathrm{d}S$ is the surface element across the wake, and ${\sigma }^2\left[u\right]$ is the variation of u. This fluctuation term originates from the averaging of the Reynolds-decomposed momentum balance equation, and it decreases the resulting force. The calculation of the last term with static pressure p is rather complicated [30], as all components of the velocity gradient are needed to be known in order to integrate the pressure Poisson equation [31]. Terra et al. [28] used tomographic PIV in order to obtain this term, and they found that, fortunately, the pressure term is important only up to the recirculation zone, and it disappears downstream; see Figure 13 in [28]. This observation allows us to ignore this term and estimate the drag force only by using the measured spatial velocity distribution as

$$\frac{〈F_D\left(x\right)〉}{\rho H}\approx {\int }_{\mathrm{FoV}}\left(U_0-〈u\left(x\right)〉\right)〈u\left(x\right)〉-{\sigma }^2\left[u\left(x\right)\right]\mathrm{d}y,$$

(2)

where the integration is performed along the y-direction in the measured two-dimensional FoV at a fixed x-position, which is far enough to ignore the pressure term. In our case, $x=0.4c$ is chosen (i.e., the most downstream position within our FoV). H is the height of the prismatic airfoil. The background velocity $U_{\infty }$ is replaced by the background $U_0$ of the Gaussian fit of the velocity profile (Equation (5)).

There is another approach to estimate the drag coefficient by using the spatial distribution of the measured velocity across the wake. Antonia and Rajagopalan [32] offer an alternative procedure including the transverse velocity fluctuations as well:

$$\frac{〈F_D〉}{\rho H}\approx {\int }_{\mathrm{FoV}}\left(U_0-〈u〉\right)〈u〉+{\sigma }^2\left[u\right]-{\sigma }^2\left[v\right]\mathrm{d}y.$$

(3)

The fluctuations in the stream-wise direction are thought to increase the drag force, while in Terra’s definition (2), they decrease the force. At higher velocities, the estimations are very similar; both are dominated by the first term of the average momentum deficit, while at middle Reynolds numbers, the stream-wise fluctuations are more important, and Antonia’s approach gives larger values than Terra’s one. At low Reynolds numbers, the wake is dominated by transverse fluctuations associated with the vortex street decreasing the drag, which is revealed when using Antonia’s formula. The latter formula is widely used in the literature; e.g., Zhou et al. [33] compared the direct force measurement by using a load-cell with this formula applied to the transverse profiles measured by using the Laser Doppler Anemometry (LDA) technique. Mohebi et al [34] uses this formula to study the flow past a flat plate at high angles of attack.

3.4. The Wake Width and Centerline

By using the spatial distribution of the average stream-wise velocity, the wake centerline can be determined. It is a set of points at which the stream-wise velocity component is minimal; i.e., the velocity deficit $u_d=U_{\mathrm{ref}}-〈u〉$ is maximal. The maximum or minimum are found by zero derivation:

$$\frac{\mathrm{d}〈u〉}{\mathrm{d}y}=0$$

(4)

The respective isolines for different Reynolds numbers are shown in Figure 10. The main line along the axis shows the wake centerline, where the average stream-wise velocity is minimal. However, there are other structures visible as well: at low $\mathrm{Re}$, the maxima caused by the acceleration along the airfoil enforced by the mass balance are depicted. At the highest $\mathrm{Re}$, some noise appears in the bottom right corner.

The wake width ${\delta }_w$ can be calculated by using multiple approaches [35] depending on which physical property has to be explored [35]; this shows that for comparison purposes, the most suitable wake width is determined as twice the $\sigma$-parameter of a Gaussian function fitted to the ensemble-average stream-wise velocity profile:

$$u\left(y\right)\approx G\left(y\right)=U_0-u_d{e}^{-\frac{{\left(y-y_c\right)}^2}{2\sigma }},$$

(5)

where $U_0$ plays the role of a background, $u_d$ is the deficit velocity, $y_c$ is the y-coordinate of the centerline, and $\sigma$ can be interpreted as half of the wake width. The Gaussian function fits well, as suggested, e.g., by [36] (where the wake past the circular cylinder is explored by using Acoustic Doppler Velocimetry). The dependence of ${\delta }_w=2\sigma$ on the stream-wise distance from the airfoil trailing edge is displayed in Figure 11 together with the deficit velocity $u_d$ for each stream-wise distance x. We see that the wake width increases with distance; however, at small Reynolds numbers, there is a minimum of the wake width. This minimum would be more apparent when taking into account the fluctuations, as shown in [35]. The observed growth of wake width is approximately linear, but, according to Eames [36], the transition to $\sim x^{0.5}$ growth can be expected downstream. The power of the wake width growth rate is locked as the power of the velocity deficit decreases by the conservation of momentum deficit within the wake. However, this rule does not need to be valid in the near wake, where the pressure field plays its own role.

The evolution of ${\delta }_w$ with $\mathrm{Re}$ at a fixed distance of $x=0.4c$ is plotted in Figure 12, as well as the wake growth rate a and the velocity deficit at a fixed distance.

3.5. Turbulent Kinetic Energy

The turbulent kinetic energy (TKE) is calculated as

$$E_T\left(\overrightarrow{x}\right)=\frac{1}{2}〈u^2\left(\overrightarrow{x}\right)-{〈u\left(\overrightarrow{x}\right)〉}^2+v^2\left(\overrightarrow{x}\right)-{〈v\left(\overrightarrow{x}\right)〉}^2〉,$$

(6)

where u and v are the stream-wise and cross-stream velocity components measured by using PIV, $\overrightarrow{x}$ is the position within the field of view, and the sharp brackets $〈·〉$ denote the averaging over the ensemble. Note that this definition of TKE contains the in-plane velocity components only; therefore, its value is roughly underestimated by the entire z-component term. The level of underestimation may vary with the regime of flow past the airfoil: at a smaller Reynolds number, the turbulent structures are large and oriented with the height of the airfoil, while at large Reynolds numbers, the turbulent structures might be small and isotropically oriented, and thus the planar definition of TKE would contain approximately two-thirds of the true TKE. Another important note is that the velocity $\overrightarrow{u}\left(\overrightarrow{x}\right)$ is already averaged over the interrogation area—it is not a true single-point quantity. Therefore, the entire amount of TKE of length-scales smaller than the interrogation area is missing.

Figure 13 shows the spatial distribution of TKE based on two measured velocity components. The first impression is the decrease of $\mathrm{TKE}/U_{\mathrm{ref}}^2$ with the Reynolds number (note the colorscale of the last row of Figure 13 is multiplied by 8). The worst product A creates a significantly smaller amount of TKE, mainly at the largest Re. A regime switch can be observed at a Reynolds number of $4.08\times {10}^4$, where the variant A lies in the previous stage with a single massive spot of TKE, while B and C wakes consist of a pair of maxima past the boundary layers.

This change of regime is even more apparent in the plot of $\sigma \left[v\right]/\sigma \left[u\right]$ in Figure 6. In respect to the original work of Romano [37], the logarithm is used in order to make the anisotropy symmetrical, because the logarithm maps the interval $(0,\infty )$ with a center in 1 to interval $(-\infty ,\infty )$ with a center in 0. Thus, the areas in which fluctuations in both investigated directions are statistically the same have a value around 0 plotted as white color; the areas in which the fluctuations in cross-stream directions are twice as strong as in the stream-wise direction have the value $+1$, and they are depicted by a blue color.

At a smaller Re, the strong dominance of fluctuations in the cross-stream direction is observed. This is caused by the von Kármán vortex street formed past the airfoil. The areas at the left-hand-side of the FoV, but not behind the airfoil, where the stream-wise fluctuations dominate, can be caused by the unsteadiness of the wind tunnel velocity, which is a known issue at smaller velocities. This fluctuations occur in the stream-wise direction only.

At $\mathrm{Re}=4.04\times {10}^4$, the regime changes in favor of the stream-wise fluctuations. However, there is a central strip of cross-stream fluctuation dominance, which is connected with the von Kármán vortex street, which is still present at those velocities. The central blue strip weakens with increasing Reynolds number, but it still exists at the largest explored velocity.

The dominance of stream-wise fluctuations is caused by the continuation of boundary layers formed along the airfoil. The length-scale of a fluctuation inside boundary layer is limited in the cross-stream direction by the thickness of the boundary layer, while the stream-wise fluctuations are not limited in their development. The level of anisotropy becomes stronger at the pressure side of the airfoil (bottom part of the figures) and at the largest explored Re.

The classical approach of investigating anisotropy [38] uses all three instantaneous velocity components measured at least in the single point by the three-wire Hot Wire Anemometer [39] or by using Stereo PIV [40,41,42]. In this experiment, we have only two velocity components, although Stereo-PIV measurement is planned in the future. In the case of 2D velocity components, the eigenvalues cannot be extracted; thus, we are limited to the ratio of fluctuations in two directions, as already used in the article by Romano [37]. Another approach is to use the so-called degree of anisotropy, where the difference plays a role instead of the ratio, as used, e.g., in the book [43]. The spatial distributions obtained by using the degree of anisotropy or by using the ratio of $\sigma \left[v\right]/\sigma \left[u\right]$ are very similar.

3.6. TKE by Length-Scale of Fluctuations

A better insight into the nature of turbulent kinetic energy can be achieved by determining which are the sizes of fluctuations producing the TKE [1]. This can be done by using our approach [24,44] of separating the fluctuations by length-scale inspired by the work of Agrawal and Prasad [45,46]. The spatial spectrum [24] is obtained without a need for temporal resolution [47]. The idea is quite simple: the spatially resolved velocity field $\overrightarrow{u}\left(\overrightarrow{x}\right)$ is convoluted with a band-pass-filter:

$${\overrightarrow{u}}_{lh}\left(\overrightarrow{x}\right)=\overrightarrow{u}\left(\overrightarrow{x}\right)\ast P_{lh}\left(\overrightarrow{x}\right),$$

(7)

where $P_{lh}\left(\overrightarrow{x}\right)$ is the band-pass filter, keeping structures of size larger than ${\sigma }_l$ and smaller than ${\sigma }_h$. It is obtained as a difference of two Gauss functions:

$$P_{lh}\left(\overrightarrow{x}\right)=\frac{1}{2\pi {\sigma }_l^2}{e}^{-\frac{x^2}{2{\sigma }_l^2}}-\frac{1}{2\pi {\sigma }_h^2}{e}^{-\frac{x^2}{2{\sigma }_h^2}}.$$

(8)

Then, the spatial distribution of turbulent kinetic energies of different bands ${\overrightarrow{u}}_{lh}\left(\overrightarrow{x}\right)$ is calculated. Figure 14 displays the combination of spatial distributions of TKE of three bands depicted by three basic colors: red for the TKE of the spatially smallest fluctuations of length interval 1.0–1.5 grid points, which corresponds to $7\times {10}^{-3}–1\times {10}^{-2}c$, c chord width; the green channel shows fluctuations of length interval 3–4 grid points, i.e., $(2.0–2.7)\times {10}^{-2}c$; and the blue colors represent the largest resolved fluctuations of sizes between 8–12 grid points—i.e., $(5.3–8.0)\times {10}^{-2}c$. The relative intensities $\widehat{I}\left(k\right)$ of each color channel have to be normalized according to the band interval $k/\Delta k^2$. Additionally, they are normalized by the famous Kolmogorov scaling [48] $k^{-5/3}$ in order to obtiain equal color intensities in the case when the power spectral density follows this scaling;

$$\widehat{I}\left(k\right)\sim E_T\left(k\right)·\frac{k}{\Delta k^2}·k^{\frac{5}{3}},$$

(9)

where the effective wave number of an explored band of fluctuations with a size between ${\sigma }_l$ and ${\sigma }_h$ is calculated as

$$k=\frac{2}{{\sigma }_l+{\sigma }_h}$$

(10)

and

$$\Delta k=\frac{1}{{\sigma }_l}-\frac{1}{{\sigma }_h}.$$

(11)

By using more bands within the spatial resolution of our data, the spatial power spectral density $E\left(k\right)/k$ can be reconstructed. It is important to note that, in respect to the classical spectra obtained from time-resolved point data (e.g., [49,50,51]), our method is limited to the range of field of view to the grid point. This interval covers only values from $0.53\mathrm{mm}$ to $32\mathrm{mm}$; i.e., one and a half orders of magnitude. The classical time-resolved point measurements typically cover tens of minutes by the resolution of tens of kilohertz covering five orders of magnitude in frequency [52]. Our approach shows the spectra calculated by using the entire field of view; each point is affected by all other points in the data ensemble. Moreover, the signal at larger scales feels larger with neighboring areas, while the smallest scale signal is almost local.

The spatial distribution of length-scale dependent TKE is shown in Figure 14. At lower Re, there are large areas displaying the relative large-scale source of TKE, which is displayed as a massive blue spot. This is caused by large almost laminar vortices in the forming von Kármán vortex street. At Re $4.08\times {10}^4$, the wake past the product A belongs rather to the previous regime, while the others belong to the next regime, characterized by continuing boundary layers containing small scale fluctuations. The wake is still affected by larger scale oscillations of the wake; thus, it is displayed in violet colors—it contains large scales and small scales, while the middle scales are depleted. At even higher Re, the importance of large scales decreases, and the asymmetry appears between the pressure and suction side of the wake—there are more middle-scale fluctuations (green) past the pressure side. At the highest explored Re $1.63\times {10}^5$, the large scales are weaker, and the wake consists of a pair of strips of middle and small-scale fluctuations (yellow) with a strip of small scales (red) in between. We mention again that the dominance of some color in Figure 14 signifies that there is slightly more energy in the corresponding length-scale-band than would be in the case of ideal turbulence following the five-thirds law. Therefore, the energy content of large-scale fluctuations at the highest Reynolds number is still larger than the two other bands combined in absolute numbers. This issue is more apparent in the plot of the power spectral density in Figure 16.

The power spectral density, shown in Figure 16, displays $k^{-3}$ scaling at the smallest Re. We mention that even steeper scaling $k^{-11/3}$ is typical for a two-dimensional arrangement of vortices [53], which is typically observed in the low Reynolds number wake past the circular cylinder. At higher Reynolds numbers, the spectra approach the $k^{-5/3}$ law. Note that the vertical axis is not normalized by $U_{\mathrm{ref}}^2$; it displays absolute numbers. Despite this, the wake at Re $1.63\times {10}^5$ displays less large-scale energy than the wake at Re $8.2\times {10}^4$, while at middle and small-scale wave numbers, the larger Re contains more energy.

3.7. Spatial Correlation

Another way of exploring the size of turbulent structures is the usage of the autocorrelation functions or the structure functions [54,55]. Here, the autocorrelation function of the cross-stream velocity component is used, because this component is active in the formation of the von Kármán vortex street observed at the lower Reynolds numbers. The autocorrelation of some quantity between two separate points $\overrightarrow{x}$ and $\overrightarrow{\xi }$ is defined as

$$R_{vv}\left(\overrightarrow{x};\overrightarrow{\xi }\right)=\frac{〈v^{\prime }\left(\overrightarrow{x}\right)·v^{\prime }\left(\overrightarrow{\xi }\right)〉}{\sigma \left[v\left(\overrightarrow{x}\right)\right]·\sigma \left[v\left(\overrightarrow{\xi }\right)\right]},$$

(12)

where $〈·〉$ signifies the ensemble averaging and $\sigma \left[·\right]$ is the standard deviation; $v^{\prime }$ is a fluctuation velocity, $v^{\prime }=v-〈v〉$. In Figure 17, the position vector $\overrightarrow{x}$ runs over the entire FoV, while the reference point $\overrightarrow{\xi }$ is fixed in the FoV center.

One of the natural properties of the correlation function is that it inherits some aspects of the basic function; e.g., its periodicity. Therefore, the periodic von Kármán vortex street is reflected as a periodic function of autocorrelation with minima in the distance of half of the spatial period. Note that the distance of the minimum is shorter in the upstream direction and longer in the downstream direction (e.g., at Re $1.63\times {10}^4$, Figure 18a, the minimum in the upstream direction occurs at a. distance of $-0.15c$ for the case A and at $-0.13c$ for B and C, while in the upstream direction, the first minimum is $0.18c$ for A and $0.17$ for B and C). With increasing Re, the dominant periodic behavior disappears, and the $R_{vv}$ starts to represent the spatial decay of coherence.

The level of spatial coherence is characterized by the integral length-scale, which is defined as the integral of $R_{vv}$ up to some distance M, where $R_{vv}$ reaches zero, $R_{vv}\left(M\right)=0$

$$L_{vv}={\int }_0^M{R}_{vv}(x-\xi )\mathrm{d}x.$$

(13)

However, reaching zero is not straightforward, as all experimental data contain noise. Thus, the autocorrelation oscillates around zero; see Figure 18, later panels. Azevedo et al. [56] explore several possibilities of the limits of this integral. Here, the simplest key is chosen: the integration is stopped at the distance where $R_{vv}$ reaches the value $1/e\approx 0.367$. This choice underestimates the integral length-scale, as shown in [56], but it does this in a most systematic way.

The values of integral length-scales integrated up to the distance where $R_{vv}\left(x\right)=1/e$ are plotted in Figure 19. Panel (a) of Figure 19 shows $L_{vv}$ along the stream-wise direction in millimeters. $L_{vv}$ decays systematically and continuously with increasing Re. Note that the interpretation of $L_{vv}$ depends on the regime: at the periodic regime, where $R_{vv}$ is a periodic function with the period of Strouhal shedding, it might represent only a constant fraction of the spatial Strouhal period (suppose that $R_{vv}=cosax$, where a is $2\pi$ times the spatial Strouhal frequency, then $L_{vv}={\left[a^{-1}sinx\right]}_0^M$, where $M,R_{vv}\left(M\right)=e^{-1}$, thus $M=a^{-1}arccos{e}^{-1}$. Rewrite $sinx=\sqrt{1-{cos}^{2}x}$; then, $L_{v}v=a^{-1}sinarccos{e}^{-1}=a^{-1}\sqrt{1-e^{-2}}=0.923a^{-1}$. ). At the regime dominated by continuing boundary layers, $L_{vv}$ might represent the statistical size of coherent structures, which is the original motivation for this quantity. However, a systematic shortening of $L_{vv}$ is observed without some evidence of regime change. However, if $L_{vv}$ is calculated along the cross-stream direction (panel (c) in Figure 19), the switch of regimes is evident. At lower Reynolds numbers, the observed $L_{vv}$ in the cross-stream direction is as large as the measured area. When the large von Kárm’an street disappears, the value of $L_{vv}$ in the cross-stream direction converges to the range observed in the stream-wise direction (panel (a) of Figure 19). The loss of the large-scale vertical correlation appears at the already discussed regime change at $\mathrm{Re}=4.04\times {10}^4$.

The autocorrelation function $R_{vv}$ is not symmetric for upstream and downstream directions, as shown in panel (b) of Figure 19. The asymmetry reaches the value of $15\%$ in favor of the downstream direction. This effect is caused by the increasing average velocity along the stream as the wake widens.

4. Conclusions

We aimed to discover how manufacturing inaccuracies affect the flow topology. The wake behind a single prismatic airfoil NACA 64(3)-618 has been studied experimentally by using the Particle Image Velocimetry (PIV) method. The airfoil has been realized in three copies with a different similarity to the ideal model. The real geometry has been measured by using the 3D optical scanner GOM Atos. We focused on the zero angle of attack at chord-based Reynolds numbers ranging from $1.6\times {10}^4$ to $1.6\times {10}^5$.

Generally speaking, the effect of model quality is weaker than the effect of the Reynolds number—the average flow pattern differs more among different velocities than among different models. The deviations in model shape shift the transition between turbulent regimes, which is most observable at a Reynolds number of $4.1\times {10}^4$, where the variant denoted A displays the same regime as for smaller velocities, while the other two realizations B and C produce wakes similar to those observable at higher velocities. This regime change is apparent in the average flow topology (Figure 5), in the decrease of average velocity deficit (Figure 12), in the spatial distribution of the turbulent kinetic energy (Figure 13), in the analysis of isotropy (Figure 6), or in the length-scale of structures producing TKE (Figure 14). This issue can be important if the flow machine (e.g., wind turbine) was designed to operate at the edge of some regime.

The worst variant A differs from the ideal geometry mainly in the quality of trailing edge: it is shorter by $2.8\mathrm{mm}$ ($3.5\%$ of the chord length), and it lacks the entire small trail pointing to the pressure side of the airfoil; see Figure 1. Thus, it produces the smallest lift (Figure 9), but also the smallest drag (Figure 20), as the other variants are slightly expanded by using the Minkowski sum. Despite this huge discrepancy, the flow characteristics of this realization are comparable to the better models. The wake centerline orientation is similar to others (Figure 10), the wake thickness reaches similar values (Figure 12), and even the length-scales of fluctuations in the wake display similar patterns (Figure 14) and spectra (Figure 16) as the better versions.

In future, we plan to compare these results with a numerical simulation covering not only the scanned geometry but also the ideal design, which was not accessible by this experiment.