Velocity data-based determination of airfoil characteristics with circulation and fluid momentum change methods, including a control surface size independence test

An experimental method for determination of aerodynamic loads is presented. It is based on velocity vector field results obtained with Particle Image Velocimetry (PIV). As PIV is an optical measurement technique, the developed method for load determination can be defined as noninvasive. It is shown that the only information needed to estimate the lift and drag forces exerted on a body placed in the flow is a velocity distribution measured around the investigated object. Therefore, PIV results provide sufficient input experimental data to be used. Fundamental fluid mechanics theories were employed to develop algorithms for load estimation. Determination of the lift force is based on velocity circulation calculations. It is obtained by integrating the velocity field along a closed-loop encircling the body. An essential achievement made is the development of a procedure for finding an optimal size of the integration curve used for lift calculations. In the case of the drag force estimation, an analysis of fluid momentum changes has been used. The momentum deficit, within a given control volume containing the analysed aerofoil, is determined and related to the reaction drag force exerted on the body. Additionally, pressure field reconstruction based on velocity data, which enabled an application of small control surfaces and kept the drag estimation error at a satisfactorily low level, was introduced. The developed method was tested and verified with the reference computational fluid dynamics simulation results and applied further to the wind tunnel experimental data. A flow around the standard NACA0012 aerofoil at two flow regimes was investigated (Re=0.7×105\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.7\times 10^5$$\end{document} and Re=1.4×105\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1.4\times 10^5$$\end{document}). Lift and drag coefficient characteristics as a function of the angle of attack were obtained. An exceptional agreement between the experimental and reference numerical lift characteristics was attained (relative differences no larger than 5%). In the case of drag estimation, an acceptable level of similarity was observed (max. discrepancies below 20%).


Introduction
Slender, aerofoil-shaped bodies exposed to the flow have been an object of experimental studies since the beginning of aerodynamic investigations that started more than a hundred years ago. It is evidently necessary to carry out research on aerofoils as they are used extensively in many engineering applications such as aircraft wings, wind turbine blades, bridge decks, towers. In experimental fluid dynamics, the investigations of aerodynamic loads apart from measurements of flow field quantities such as velocity, pressure, and temperature are of special interest. A common experimental research practice is to use separate techniques for aerodynamic force and flow measurements (e.g. force balance measurements vs pressure probe measurements). However, it is obvious that there exists a physical link between the flow behaviour and the body forces generation mechanism. Therefore, the literature review shows numerous studies in which the aerodynamic loads are derived from the flow-field information.

Determination of loads based on the velocity field
A good example of a non-direct measurement approach is a well-established experimental technique for the determination of aerofoil profile drag from the momentum deficit in the wake described by Betz (1925). In this case, the method seems to be more accurate than direct force balance measurements as it does not include the induced drag always influencing the profile drag measurement. Another great advantage of this technique is a possibility to obtain a sectional drag and its variation with a wingspan. The above-mentioned method has been successfully developed and applied since the beginnings of wind tunnel measurements as it requires only simple pressure-based velocity measurement (e.g. with a rake probe). A more advanced method for load derivation based on flow variables has been proposed over the last two decades together with the development in high-resolution velocimetry techniques like Laser Doppler Anemometry (LDA) or PIV. A vast majority of researchers working in this field exploit the control-volume approach for the determination of both lift and drag. Aerodynamic forces can be obtained from an integration of flow variables over a control volume (and/or boundary surfaces) surrounding the object submerged in the fluid. Another method in use exploits the Kutta-Joukowski theorem which relates velocity circulation with aerodynamic lift.

Control volume approach
The flow parameters required for the control volume method are velocity, pressure, density, and viscous stress. The last one can be neglected as long as the close vicinity of the body is not considered. For the incompressible flow, the density is constant, thus the remaining components of the force equation require a distribution of velocity (and its derivatives) within the control volume and a pressure distribution over the surrounding surface. A detailed control volume theory with the comprehensive derivation of equations was presented in a series of works by Noca (1996), Noca et al. (1997), Noca et al. (1999). This method requires the knowledge of the pressure field p. For the sake of utilizing the method with velocity data only (obtained, for instance, by the PIV technique), the pressure term needs to be eliminated. An available solution is either an algebraic transformation of the momentum equation as shown by Noca et al. (1999) or an explicit evaluation of pressure with the momentum equation proposed by Unal et al. (1997).
A clarified summary of the above variants of the control volume approach was presented by van Oudheusden et al. (2006). The authors considered a simpler case, however, by reducing particular components of force equations. The time-resolved flow consideration was exchanged with a time-averaged flow approach, which is sufficient for many applications of technical interest. Such an approach simplifies computations as it is not necessary to evaluate flow variable integrals over the whole control volume. Surface integrals (or line integrals for the 2D case) along the closed contour encircling the body can be used instead. Apart from the contour method, drag was additionally obtained by an application of the classical wake approach as in the well-established method by Jones (1936). The authors reported significant improvement of drag estimation accuracy with this approach in comparison with the first method. For the contour method, unacceptably large errors were achieved due to a low relative value of the drag of the investigated aerofoil. A complex, multi-level comparison of different methods for calculation of time-averaged aerofoil loads was performed by Albrecht et al. (2012). The accuracy of various methods was evaluated on a common set of data. The comparison included methods based on the momentum equation requiring the knowledge about the pressure that was reconstructed on the basis of velocity data (from the differential momentum equation or the Poisson equation). The second class of methods investigated therein was the equations given by Noca et al., in which pressure was eliminated via integral identities. Both numerical and experimental data were used in this study. The authors evaluated an effect of the boundary position and size as well. An analysis of the presented results justifies the need to perform contour size independence tests (further details can be found in the original publication), especially for methods involving various approximation and pressure reconstruction algorithms.

Circulation approach
Several authors have recognized an alternative method to extract information about the aerodynamic lift from the velocity field around the aerofoil section. Instead of control volume with the momentum equation approach, the Kutta-Joukowski theorem (see Eq. (3)) can be used, taking advantage of the velocity circulation around the profile. In Gibbs et al. (2005), the velocity circulation was measured to illustrate a span-wise distribution of the lift force for a wing in formation flight. The method was used by Dobrev and Massouh (2011) to determine sectional loads of the rotating wind turbine blade. In Sharma and Deshpande (2012), the lift force generated by a two-dimensional thin flat plate at various angles of attack in the low Reynolds number flow was experimentally determined through the velocity circulation measured with the Laser Doppler Velocimetry (LDV) technique. As reported, the results proved the successful use of Kutta-Joukowski theorem in determining the lift even in the post-stall region of the lift curve without any need for data correction. The circulation approach can be effective in the case of viscous flows as shown in Li et al. (2015). A good example is a study by Wang et al. (2019), in which complex viscous unsteady wakes generated by flapping wings were investigated.

Study motivation
An analysis of the cited investigations proves that it can be beneficial to reduce the complexity of the measurement procedure by deriving some physical quantities from the others. A fewer devices are required; thus, the measurement procedure and data acquisition are simplified. In specific situations, parameters hard or impossible to be measured directly can be obtained (e.g. sectional loads on the wing). Certain measurement techniques like pressure taps or pressure probes provide information only at discrete points that cannot be changed during the test. Velocity-field-based methods, on the other hand, can be treated as a source of the data quasi-continuous in space if a fine enough spatial resolution is achievable. What is more, within such an approach, an easy synchronization between the measured and derived parameters is achieved. Finally, there is no need for intrusive instrumentation of the model and undesirable flow disturbances.
There are plenty of practical flow cases for which the method of load determination based on velocity fields can be more viable than standard methods (e.g. strain gauges-based measurement systems). These are freely moving objects, flexible bodies, flying and swimming animals, low Reynolds number flows. One should add here a set of experiments for which an application of standard force measurement systems can be problematic due to purely technical reasons, e.g. determination of loads on rotating blades - Berton et al. (2004), Dobrev and Massouh (2011).

Gaps in the literature
Evident advantages and the application potential of velocity-based load determination techniques were not the only inspiration to continue the work on that subject. One of the major aims of this study was to create simple, effective, and accurate algorithms. For most of the investigations in this field, the aerodynamic force is based on more demanding procedures requiring either surface integral calculations or a pressure reconstruction. Unfortunately, this approach can result in lowering the accuracy of force determination due to an accumulation of numerical errors (when larger areas are used for integration), inaccurate pressure reconstruction algorithms, imperfections of the acquired PIV velocity fields (artefacts due to a lack of optical access or laser light reflections). For this study, therefore, it was of interest to investigate a simple and efficient circulation approach for lift calculations. To our knowledge, no study has yielded a detailed analysis of this technique. In Dobrev and Massouh (2011), where the circulation method was applied, the authors do not examine the size of the contour and no verification of the obtained lift results was proposed neither. In Gibbs et al. (2005), an attempt to expand the velocity circulation measurement procedure was proposed by using partial contours apart from the full ones. Large discrepancies in magnitudes of circulation were observed, however.
For the drag force determination, a wake analysis was used as it was confirmed to be more accurate than control volume methods -see van Oudheusden et al. (2006). We suggest even a further simplification in the method described in the next section. The most important advantage of this approach is that it can perform very well for incomplete velocity fields used as the input data.
The main practical problem we are confronted with is that the size and position of control volume/surface or 2D field boundaries (contours) used for velocity-derived data integration have been rarely studied in detail. An assumption of a negligible influence of size parameters is only true for the inviscid potential flow (as in the Kutta-Joukowski theorem, for instance) but for real flow conditions or experimental results with possible inaccuracies or artefacts, some additional actions should be undertaken.
In van Oudheusden et al. (2006), an attempt to include this influence was made by an estimation of the uncertainty in the load data from the contour integral while varying the distance of the contour to the aerofoil. For the wake method, the location of the drag line behind the trailing edge was changed. Such a procedure provided an estimate of the standard deviation of the obtained load results due to the contour modification. This approach, however, should be expanded with an alternative method of finding one optimal contour size to be demonstrated in further sections.

Airfoil characteristic reference data
For validation of in-house numerical simulation results, reliable reference data had to be found in the available literature. An interesting approach for the choice of proper experimental data containing aerodynamic characteristics of NACA0012 that can be used for verification purposes was proposed in McCroskey (1987). The author presents a critical assessment of the results from more than 40 2D wind tunnel experiments with a NACA0012 aerofoil. Experiments were carried out in the range of approx. Re = 4 × 10 5 to 4 × 10 7 both for subsonic and transonic flows. Surprisingly, it is reported that the scatter in the total ensemble of data is unacceptable and it is not clear which of the experimental data can be regarded as fully correct. Additionally, it is concluded that there exists no single experiment adequate for defining complete aerodynamic characteristics of the NACA0012 aerofoil or validating CFD codes. On the other hand, aggregating of collected data after a certain filtering process enabled ending up with a set of quantitative information. Based on the experiments conducted with the utmost care by a few research centres, plots of several parameters as a function of Reynolds and Mach numbers were drawn. This allowed one to obtain best-fitted curves and approximate them in an equation form. Two of them attracted a particular interest to be used herein for data validation. These are the lift-curve slope C L and the minimal drag coefficient C D 0 (at the angle of attack = 0 • ) as a function of the Reynolds number Re expressed by the following formulas The experimental data at Re = 1.7 × 10 5 - Jacobs and Sherman (1937) and Re = 3.6 × 10 5 - Sheldahl and Klimas (1981) were chosen for a comparison as shown in Fig. 4a. The first results were attained in the experimental study performed by Jacobs and Sherman in the NACA variabledensity wind tunnel at Langley Laboratory. The cited work was not included by McCroskey (1987); however, other experiments performed at Langley Laboratory, e.g. Abbott and Von Doenhoff (1959); Ladson (1988), were classified into a group of results characterized by the lowest scatter. (1) The results obtained by Jacobs and Sherman provide characteristics of a few basic NACA aerofoils over a wide range of the Reynolds number. The second mentioned experimental study of Sheldahl and Klimas was carried out in Sandia National Laboratories assigned by McCroskey to be the second-best in terms of the result reliability. Additionally, the synthetic drag data prediction for Re = 8 × 10 4 obtained by Sheldahl and Klimas with an airfoil section characteristics synthesizer computer code was included (see Fig. 4b).

Research methodology
The concept of methodology involved in the presented research is illustrated in the flowchart in Fig. 1. The process of developing a method for velocity-based aerodynamic loads determination consisted of two phases. Before applying a selected mathematical model to the experimental data, a training procedure was performed on synthetic CFD velocity fields (phase I). The aims of this stage were twofold: first, to verify the correctness of the governing equations presented in Subsect. 2.1 and second, to find optimal size parameters of contours used to compute integrals (see Sect. 3). From the moment the trained model was obtained, it could be applied to compute aerodynamic loads from the raw experimental PIV data (phase II). CFD was assumed to be a reliable source of the synthetic reference data mimicking PIV velocity fields. It provides complete information on the velocity at every node within the defined numerical domain. No discontinuity or missing data is expected within CFD results as opposed to the experimental data. PIV velocity fields can be contaminated with various imperfections such us unwanted artefacts influencing the magnitude of measured velocity or missing vectors.
Another reason to run numerical simulations is that additional tools to obtain the aerodynamic loads are provided. The lift and drag force values calculated with velocity-based algorithms can be compared to the forces from numerical simulation solutions. Additionally, the choice of CFD as a validation method for the 2D flow cases under investigations was considered as the most beneficial due to the weaknesses of available experimental methods. Firstly, in the case of aerodynamic balance measurements, it could be problematic to determine sectional loads easily. In such circumstances, total lift and drag forces need to be corrected due to a finite wing model span. 3D flow effects cause an increase in the total drag due to a significant contribution of the induced drag at the model tips. On the other hand, measurement of a static pressure distribution over the aerofoil surface involving pressure taps would allow for determination of sectional loads. In this case, however, viscous effects could not be observed and, thus, the skin-friction drag would be neglected. For flows around aerofoils in the case of no separation (at small angles of attack), a contribution of the skin-friction drag to the total drag value is predominant. Hence, such an experimental procedure needed to be rejected as well. A similar approach with an application of the CFD results for validation of velocity-based algorithms to evaluate loads was also successfully used by other researchers as indicated in Sect. 1.1, e.g. see Unal et al. (1997); Noca et al. (1999);van Oudheusden et al. (2006).

Mathematical model
The computational model for the determination of aerodynamic forces coefficients implemented in the developed method is based on the well-known concepts presented below.

Lift force
The phenomena of lift generation and its relation to circulation is a starting point for an explanation of the lift force determination procedure. A general remark is such that whenever the circulation around the body appears, the lift is produced as well. Joukowski proves that when a cylindrical body of an arbitrary cross section moves with the velocity U ∞ in the fluid of the density ∞ and the circulation around it has magnitude , then a lift force per unit span L ′ equal to the product presented in E. (3) is produced.
To define the term of circulation, a simple case of the circulatory flow can be considered. For such a flow, all fluid particles move around one centre point along the circular streamlines. The circulation can be calculated for an arbitrary streamline as a product of the velocity v and the circumference of the streamline 2 r . Generalizing this definition for any flow, the circulation is equal to a product of the mean velocity along an arbitrary closed curve and the length of that curve. Such a product has a constant value independent of the curve shape and size if the flow is vortex-free (potential). The above definition can be expressed as well in an integral form as shown in Eq. (4).
where u is a velocity vector and d is an infinitesimal length segment of the closed loop C.
The lift coefficient is a non-dimensional quantity obtained by dividing the lift force per unit length L ′ by the chord length c and the dynamic pressure 1 2 ∞ U ∞ 2 .
A simplified formula for the lift coefficient estimation based only on the circulation and the freestream velocity can be derived by substituting Eq. (3) into Eq. (5).

Drag force
One can distinguish various contributions to the total drag force exerted on the body; however, the objective of this work is to focus on the profile drag (or the boundary-layer drag) which is a sum of the skin-friction and form drags. For a two-dimensional incompressible flow, it is a total force in the stream-wise direction exerted on the body. Simultaneously, a reaction force of the same magnitude and an opposite direction is applied to the fluid, giving rise to a reduction in its momentum. Thus, one can determine a value of the profile drag knowing the velocity distribution upstream and downstream of the aerofoil provided that the static pressure is constant. A value of force per unit span can be evaluated with the momentum-based method expressed by the following formula The above formula can be interpreted as a difference between the fluid momentum calculated at profiles I and II (refer to Fig. 7 (8) is obtained. It should be noted, however, that for the mass conservation to be valid, the height of profiles I and II need to be different (which may not be manifested enough on a simplified scheme in Fig. 7).
Finally, Eq. (9) can be used to rewrite Eq. 8 expressing the drag force in terms of a velocity deficit as below The integration of velocity is performed here only along the downstream profile. At vertical positions far enough above and below the aerofoil trailing edge level, the integrand achieves zero value. Due to this fact, from the practical point of view, it is enough to integrate only the velocity deficit region in the wake disregarding the velocity distribution in the remaining part of profile II. This simplified variant of the momentum-based formula will be referred to hereinafter as a wake method.
In the above considerations, the condition of a constant static pressure along upstream and downstream profiles needs to be satisfied, which is valid if measurements are performed at a relatively large distance from the body. Otherwise, one needs to take into account pressure forces in the fluid momentum analysis. Both velocity and pressure distributions along profile I and II should be assumed as non-uniform then. Taking this into consideration, the following general formula with an additional pressure integral term should be applied to determine the profile drag (provided no additional mass or momentum source enters the 2D control surface from top and bottom) Equations 7, 10 and 11 provide a set of tools for profile drag determination. The method of the drag determination based on PIV velocity fields does not include pressure measurements; thus, in the last case, an application of a pressure reconstruction with Bernoulli's formula for an incompressible flow is required (see Subsect. 3.3).
As a dimensionless measure of the drag, one can introduce the profile drag coefficient C D expressed by the following formula, analogically to the lift coefficient C L in Eq. (5)

CFD simulations
CFD simulations of the virtual wind tunnel were performed in the ANSYS 14.0 environment. The numerical study aimed to model the real experimental tests performed in the wind tunnel with a closed test section arrangement. The flow around the NACA0012 aerofoil located in the centre of the closed test section was investigated. Simulations were run for a wide range of the angles of attack ( = 0 • -17 • ) covering the linear part of lift characteristics, static stall, and a few angle settings with large flow separation from the aerofoil surface. Two freestream flow velocities were investigated-U ∞ = 11 m/s and U ∞ = 22 m/s , which correspond to the Reynolds numbers Re = 0.71 × 10 5 and Re = 1.42 × 10 5 , respectively (for the aerofoil chord c = 0.1 m).
A quasi-2D analysis, i.e. a flow around the infinitely long profile, was simulated. To do that, a thin slice of the geometry representing the cross section of the wind tunnel test section was prepared and meshed with one-element-layer across its depth. The size of the numerical domain corresponded to the dimensions of the wind tunnel test section, i.e. the height h = 1.36 m and the length l = 0.97 m. The test model-NACA0012 profile-was approximated with 140 points on which the spline curve was based.
A combination of the structured and unstructured hexahedral type of a grid was chosen for meshing the geometry. (11) The structured mesh was used for outer domains, whereas the unstructured one for the inner domain in the proximity of the NACA0012 profile and the wake region. A careful mesh refinement, as well as a mesh size independence test, were performed to guarantee a good quality of the numerical solution. The final mesh consisted of approx. 5.2 × 10 5 nodes and 2.6 × 10 5 elements. Transient analyses were run with a constant timestep -5 × 10 −4 s (the value chosen on the basis of the timestep length independence test). A uniformly distributed velocity profile with a medium turbulence intensity of 5% was set at the inlet boundary, and an average pressure condition was set at the outlet. All walls in the numerical model, i.e. an aerofoil surface and test section walls (top and bottom) were defined as smooth and non-slip. The symmetry boundary condition was applied to the sidewalls of the quasi-2D fluid, domain. A SST turbulence model with the standard definition of model coefficients was chosen Menter (1994). Additionally, the transitional turbulence-gamma-theta model was applied, which was motivated by the in-house experience gained as regards simulations of external flows at moderate Reynolds numbers ( Re = 1 × 10 5 ) Błaszczak and Sobczak (2010); Olasek and Karczewski (2012).
The residuals of, for instance, momentum and mass conservation equations were monitored during the simulation process in order to evaluate the convergence of the numerical solution. Stabilization of residuals below the desired level of 1 × 10 −7 was satisfied for all simulated cases varying in the aerofoil angle of attack. The chosen level can be regarded as an exceptional solution quality with a very tight convergence indicating low errors in the solution to the system of equations, which makes it suitable for scientific applications.
The procedure of the solution quality evaluation consisted in monitoring values of the non-dimensional wall distance y+ parameter at the aerofoil wall surface. The element size at the aerofoil surface was adjusted in a way allowing one to obtain its values close to or below 1, which guarantees that the flow is completely resolved near the wall without any approximation with wall function models.

Experimental stand & PIV measurements
The experimental measurements were performed in a wind tunnel facility at the Institute of Turbomachinery -[pl. Instytut Maszyn Przeplywowych] (IMP), Lodz University of Technology (TUL). It is an open-return, blow-down type wind tunnel. The closed test section dedicated to 2D airfoil flow measurements (height 1.36 m × width 0.29 m length × 0.97 m) was used in the presented investigations. For this configuration, the maximal achievable freestream flow speed was U max = 22.5 m/s. Upstream of the test section, the wind tunnel channel includes a converging nozzle (contraction rate ca. 3.4) as well as a stabilization section with a set of grids and a honeycomb. All these elements provide flow conditions of a good quality characterized by a uniform velocity distribution (a symmetrical 'top-hat' velocity profile with spatial variation less than 1.5% at U max for ca. 86% of the test section width) and the turbulence intensity below 1% measured in the middle of the test section at U max .
The model of the NACA0012 aerofoil was manufactured with the Electrical Discharge Machining (EDM) technique. The desired shape was cut out in aluminium, obtaining a model with a smooth surface of the measured roughness Ra = 0.84 m . The span of the profile was adjusted to the width of the test section l = 0.29 m (wall to wall arrangement) to avoid tip effects and guarantee 2D flow conditions.
An optical arrangement of the PIV apparatus was designed to determine the full velocity vector field around the aerofoil. This provided complete information on the velocity distribution along a closed-loop encircling the aerofoil to calculate the circulation and to estimate the lift. The velocity distribution at locations far enough upstream and downstream was captured as well, allowing for drag force determination.
The test section experimental setup is depicted schematically in Fig. (2). Wind tunnel tests were carried out for two free-stream velocities, namely: U ∞ = 11.6 m/s and U ∞ = 22.2 m/s (equivalent to the Reynolds numbers Re = 0.75 × 10 5 and Re = 1.43 × 10 5 , correspondingly). The NACA0012 profile was inclined at various angles of attack ∈ [−18 • , 18 • ] , changed every 1 • close to the stall angle, and with a larger step for the rest of its range. For each angle of attack, 100 double-frame images were recorded with two PIV cameras, which allowed one to determine 100 instantaneous velocity fields. These were further used to calculate the time-averaged resultant velocity field per each angle of attack.
For the measurement series, the delay time between two PIV frames was dt ≈ 75 s. It was adjusted to obtain an average displacement of tracer particles on raw images of around 10 pixels. The Davis 7.0 software was used only to calculate vectors. All manipulations on the obtained velocity data (such as contour maps generation, profiles plotting.) were performed in MATLAB, with the PIVMat Toolbox ver. 3.03. Crucial PIV pre-and postprocessing parameters are listed below -Image preprocessing: -subtracting the sliding background (with the length scale of around 10 pixels) improving the contrast of tracer particles on a raw image by removing the bright background. -Vector calculation parameters: -cross-correlation mode, -multi-pass algorithm with various interrogation window sizes -3 initial passes (64 × 64 pixels) and 2 final passes (32 × 32 pixels), -overlapping of windows − 25% for the initial pass and 50% for the final passes.
-Vector postprocessing: -filtering of poor quality vectors by deleting results with the peak ratio Q < 1.5 (i.e. vectors with the peak height lower than 150% of the noise value were rejected), -median filter, -3 × 3 smoothing filter (only for displaying and printing purposes) The laser sheet optics was illuminating the observed area from above, creating thus a light plane parallel to the flow (see Fig. 2a). Light-sheet thickness was adjusted to be equal to ca. 1 mm in the middle of the observed measurement plane. A mirror located at the bottom of the test section was used to reflect the laser sheet to illuminate a region at the bottom of the aerofoil. As presented in Fig. 2b, PIV cameras were placed in a vertical arrangement to observe the flow both above and below the aerofoil. Around 35% of cameras' fields of view overlapped (see the region marked in orange in Fig. 2b) to allow for easy merging of separate 2D PIV vector maps (for detailed image calibration procedure please refer to Appendix 1). Camera sensors were distanced around 0.8 m from the measurement plane. The camera lens aperture size f4 and the focal length of ca. 35 mm were used. Each camera observed a 0.35 × 0.35m field of view which provided an optical resolution of ca. 6.1 pix/mm (155ppi). For a specified interrogation window size and overlapping function, the final resolution of the velocity vector data was equal to 2.65 mm.

PIV error estimation
The accuracy of PIV velocity measurements depends on systematic and random errors. One can identify a few contributions to velocity uncertainty in the performed experiments that further propagate to the circulation and fluid momentum estimation.
Random components of uncertainties are mainly related to cross-correlation uncertainty and flow fluctuations. A typical fit of the correlation peak is associated with a 0.1 pixel error- Westerweel (1993). Taking into account the optical resolution of the PIV experiment of ca. 0.17 mm/pix, the number of taken PIV samples n=100 and the inter-frame delay time dt = 75 s, which translates to the velocity uncertainty of 0.022 m/s and the relative uncertainty of 0.1% (with respect to U ∞ = 22.2 m/s).
A measure of fluctuations of the flow can be assessed from vector fields with a standard deviation of velocity. It was shown that for a free-stream part of the flow u = 0.25 m/s, which corresponds to the uncertainty of 0.025 m/s and 0.11%, respectively, for n=100 samples. In the wake region, accurate averaging is harder due to higher flow fluctuations. This may significantly influence the uncertainty of the drag force determination, which relies mainly on the probing velocity from the wake. However, the present paper concludes in the next sections with specific conditions for which accuracy of the drag determination increases (i.e. the minimal distance of 2x chord lengths downstream of the trailing edge to probe the velocity data). At a certain distance downstream of the airfoil, a degree of fluctuations is lower, thus averaging is more accurate. For the downstream position > 2c , the standard deviation of velocity drops below 1.5 m/s at u = 25 m/s (mean velocity at this location) for the angle of attack = 5 • , 3.5 m/s at 22 m/s for = 10 • and 6.5 m/s at 21 m/s for = 15 • . This corresponds to the velocity uncertainty equal ca. 0.15 m/s, 0.35 m/s and 0.65 m/s, respectively, and the relative uncertainty of 0.68%, 1.6% and 2,9%, correspondingly.
The main systematic contribution to the velocity uncertainty identified was peak locking. While adjusting the optical parameters at the test stand, careful attention was paid to defocus slightly the image to obtain a seeding particle image of at least 2 pixels in size. Such a procedure helps to obtain good conditions under which peak locking should not exceed the 0.1-pixel accuracy - Westerweel (1993). This yields the velocity uncertainty of ca. 0.23 m/s that is 1% of the relative uncertainty with respect to U ∞ .

CFD reference results verification
The reliability of numerical simulations was confirmed by a comparison with the in-house PIV results as well as the literature references on airfoil characteristics. CFD and PIV contour maps of the velocity distribution are enclosed in Fig. 3. Three exemplary angles of attack are displayed to visualize the flow around the aerofoil at specific stages (attached flow, transitive stage close to static stall, and large separation). Additionally, plots of the lift and drag coefficient as a function of the angle of attack are presented in Figs. 4 and 5.
A qualitative analysis of numerical and experimental velocity fields indicates a very high level of agreement. Numerical predictions perfectly match the velocity distributions measured at the experimental stand. Patterns of the velocity increase on the suction side as well as a decrease on the pressure side are consistent. What is more, the separation region size is comparable in the CFD and the PIV for each angle of attack. A slight difference can be noticed only in the case of the wake height which is smaller on the CFD contour plots. This can be attributed to the turbulence model that tends to damp flow fluctuations more excessively than in the experimental case.
In Fig. 4a, apart from own numerical results and the literature data, a statistically derived best fit lift-curve slope C L (slopes for Re = 0.7 × 10 5 and Re = 1.4 × 10 5 are almost indistinguishable on the plot) according to Eq. (1) was displayed on the plots. In Fig. 4b, along with the drag characteristics, two points with the drag coefficient at the null angle of attack C D 0 calculated for Re = 0.7 × 10 5 and Re = 1.4 × 10 5 with Eq. (2) were added. For better recognition of the condensed data at low angles of attack, a detailed plot for the pre-stall condition (up to = 10 • ) is provided in Figure 5.
The evaluation of the data presented on C L and C D plots leads to several observations. Firstly, the stall angle was well predicted in the IMP CFD simulations. Similarly as in the experimental data, the critical angle can be identified in the However, the position of the maximal lift coefficient C L max remains in a very good agreement with the rest of the results.
The value of the maximal lift coefficient C L max achieved in the in-house numerical simulations for higher Reynolds numbers ( Re = 1.4 × 10 5 ) corresponds very well to the experimental data by Sheldahl ( Re = 3.6 × 10 5 ). The CFDbased value is underestimated with respect to the reference value by only 2.7% . In the case of the second reference characteristics by Jacobs ( Re = 1.7 × 10 5 ), CFD parameters are overestimated by ca. 13.2% and 2.6% for Re = 1.4 × 10 5 and Re = 0.7 × 10 5 , respectively. There are several reasons for the observed discrepancies. Firstly, it can be simply explained by a considerable difference in the Reynolds Numbers between particular data sets. Secondly, the conditions of reference experiments were not precisely defined or are unknown. For instance, the model roughness, which was not reported, exerts a significant influence on the quantitative lift results. Finally, the overestimation of CFD results with respect to the experimental data by Jacobs, obtained for the comparable Reynolds number, can be explained by 3D effects possibly occurring during the tests. The experimental maximal lift force coefficient could have been reduced due to the induced drag. Nevertheless, an expected tendency was observed for the lift coefficient value changes depending on the Reynolds number variation for the in-house CFD results. An evident fact is that an increase in the maximal lift with a Reynolds number raise occurs, which can be observed on collective plots in McCroskey (1987) or on the aerodynamic characteristics of specific aerofoils measured at various Reynolds numbers in Jacobs and Sherman (1937).
Considering the slope of the lift curves as a function of the angle of attack, a satisfactory agreement can be observed. The general trend was maintained, when the numerical and experimental results for the pre-stall part of the characteristics are compared. All the characteristics increase monotonically without any discontinuities. There are, however, certain discrepancies of various nature for low and high angles of attack. The CFD simulations significantly overestimate the lift data for a lower angle of attack referring to the results provided by Jacobs and Sherman and the empirical slope by McCroskey (similarly as in the case of C L max value, it can be caused by 3D effects like the induced drag occurring on the experimental stand). On the other hand, an exceptional agreement can be observed for > 6 • . Both the CFD and experimental curves start to diverge from the constant slope line up to the critical angle.
Furthermore, a difference in the shape of pre-stall portions of the lift curves can be observed. The CFD characteristics are fully concave for the whole range of up to the critical angle. The same concerns the experimental curve by Sheldahl and Klimas but is not so significant for the Jacobs and Sherman data. Such a variability in the lift curve slope is a characteristic feature for low and moderate Reynolds  number flows in contrary to the aerofoil experiments with Re > 10 6 , for which a typical linear lift characteristics can be observed before the critical angle is achieved. The aforementioned phenomena were reported, for instance, in Martínez-Aranda et al. (2016), where NACA0012 aerodynamic characteristics for low-to-moderate Reynolds numbers were investigated (it should be remarked, however, that the study concerned finite aerofoils with various aspect ratios instead of 2D aerofoils).
While analysing the drag coefficient characteristic plots, one can see some more significant discrepancies. However, a satisfactory agreement was still achieved between the IMP CFD data and the reference results in terms of curve trends. For all of the complete results (only numerical results cover the whole range of the angle of attack; no experimental data were available for moderate Re and post-stall angles), the curves can be split into three stages: I. slight, gradual drag increase up to the critical angle; II. rapid drag raise past the stall; III. further drag increment with a smaller slope.
A detailed observation of the drag curves before stall shown in Fig. 5 demonstrates a very good agreement between the IMP CFD results for Re = 1.4 × 10 5 and the reference data by Jacobs ( Re = 1.7 × 10 5 ). The value of the CFD-based minimal drag coefficient C D 0 differs only by 1.4% with respect to the given experimental data. When compared to the drag characteristics by Sheldahl for Re = 3.6 × 10 5 , the in-house numerical results are highly overestimated. An explanation for this phenomenon, similarly as in C L max case, is a relatively high difference in Reynolds numbers for the cases under investigation. As reported by McCroskey, C D 0 is highly dependent on Re changes and is becoming more sensitive under low Re conditions. A comparison with the synthetic drag results by Sheldahl for lower Re shows a better agreement; however, a significant shift by a regular offset can be observed. At the same time, a good agreement was achieved with the statistically derived parameter from Eq. (2) by McCroskey. The CFD prediction varies only by ca. 6% and 11% for the Re = 0.7 × 10 5 and Re = 1.4 × 10 5 simulation case, respectively.
Taking into account the observations given above, the in-house CFD results can be considered correct and reliable thanks to validation with the experimental reference data. A considerably high difference in the Reynolds number between the available experimental data and the performed numerical simulations was identified as the dominant source of the occurring discrepancies. This could not be avoided due to a lack of the proper reference data at the desired Re.

Lift calculation procedure & model training with the CFD data
The procedure used for the lift coefficient computation based on Eq. (6) was performed in a few steps, namely: (1) Definition of the contour encircling the aerofoil, (2) Mapping of the velocity from the CFD velocity field onto the contour (velocity data available only in the discrete positions corresponding to CFD mesh nodes); the final data resolution (velocity vectors spacing) used for all calculations was equal to 1mm (for detailed description of CFD data resolution please refer to Appendix 2), (3) Numerical integration of the velocity components tangent to the contour to determine the circulation value, (4) Determination of the lift coefficient for the given freestream air velocity and the airfoil chord length. According to the Joukowski theorem, the circulation has a constant value independent of the contour shape and size. However, a significant sensitivity of the lift result due to a variation in size parameters was observed. To find the optimal size of the integration path used for the calculations was a crux then. The basic, rectangular shape of the contour was tested in this procedure, and multiple calculations varying its size were performed. Rectangles were always centred at the point located halfway the airfoil chord line. Thus, changes in the contour height and width were made symmetrically. As a result of such batch calculations, a lift coefficient was obtained for the defined range of contour sizes. The results can be visualized in the form of 3D surface plots. Exemplary plots are shown in Fig. 6. Surface plots obtained for two angles of attack, = 7 • (in the linear part of the lift curve) and = 12 • (post-stall region) are presented.
An analogous pattern of the 3D surface plot is obtained for each angle of attack. A few characteristic plot elements were marked and described in Fig. 6b. An asymptotic character of the plotted surfaces can be observed. The lift results tend to achieve slightly lower values for larger sizes of the bounding contour, which can be explained by the dissipation of circulation far from the aerofoil. On the other hand, when contour boundaries approach the aerofoil, the computed lift value gradually rises. Its local maximum is reached for the height and width corresponding to the rectangle of a small size. The rise in the lift has naturally a certain limit. For most of the 3D plots under analysis, it can be observed that a sudden lift drop occurs below a certain height limit. It was recognized as the moment in which the contour crosses the boundary layer or a separation region for higher angles of attack. The lift decreases rapidly also at the relative contour width x∕c < 1 , which corresponds to the situation in which vertical contour segments do not cross the flow but a portion of the aerofoil instead. All surface regions with high gradients of the lift should be regarded as non-physical. The plateau in the plot indicates a location of the correct solution, and the only issue is to find the lift coefficient value closest to the expected real one.
The developed procedure required an adequate criterion to determine the correct value of the lift coefficient. The obtained results were tuned with the CFD reference data, for which the lift characteristics was precisely determined (see Fig. 4a). When the surface plot results are related to the desired coefficient values for the specific angle of attack, it could be seen that the best correlation was achieved when the maximal value from the surface plateau was assumed. However, special attention should be paid to avoid detecting the maximal lift values outside this region as shown in Fig. 6b. A conclusion can be drawn that the rectangle size should be chosen in such a way as to encompass the investigated aerofoil as closely as possible. On the other hand, the flow region with significant velocity disturbances or deficits cannot be crossed (separation, boundary layer, aerofoil itself, etc.). The same was observed true for the whole angle of attack range, which makes it a clear and unambiguous criterion for further determination of the lift coefficient from the experimental data.

Fluid momentum change method
To determine the drag force by an analysis of changes in momentum (Eq.7 & (11), a control surface concept needs to be introduced. In Fig. 7, a control surface with its specific elements is illustrated. For the drag model expressed by the given equations, top and bottom borders have to be defined along actual air streamlines. There is more freedom regarding a choice of the position and height of profiles I and II (upstream and downstream, respectively). However, as long as the position and height of profile I can be chosen arbitrarily, the determination of the profile II height is not so trivial. It depends strictly on endpoints of profile I, from which the position of endpoints of profile II is derived by tracking along the top and bottom bounding streamlines from the x I to x II streamwise location. The next step performed is a determination of the velocity distribution along profiles I and II. Only streamwise velocity components are taken into account as transverse components (tangent to the vertical borders of the control surface) do not contribute to the momentum transfer outside the control surface. An exemplary plot of the velocity distribution along profiles is shown in Fig. 8. The influence of the aerofoil is visible. Primarily, a significant fluid velocity deficit in the wake can be noticed. It is a direct evidence of the fluid momentum loss that is proportional to the drag force exerted on the aerofoil. Another observation that can be made is a non-uniform velocity distribution in the remaining part of the plot. This concerns the upstream distribution as well. Velocity values higher than the freestream velocity (for the presented CFD case, a constant value U ∞ = 11 m/s) were achieved in the upper portions of the plots. The opposite situation can be seen for the lower portions. The phenomena are even more evident for the oncoming flow, which indicates a significant upstream influence of the tested obstacle. The aerofoil inclined at a positive angle of attack generates an upward lift force due to the formation of lower and higher pressure regions on the top and bottom surface, respectively. At the same time, the opposite effect can be observed for the velocity distribution as the total pressure of the flow around the aerofoil is assumed to be constant (apart from the wake region). The non-uniform velocity distribution follows from the choice of profile positions in the proximity of the aerofoil, where the assumption of constant static pressure is invalid. To guarantee a uniform static pressure distribution, multiple PIV measurements (far upstream and downstream of the aerofoil) are needed or a relatively large velocity field around the aerofoil should be captured. Additional measurements are impractical and time-consuming, whereas capturing a larger field of view leads to a decrease in the measurement resolution and, thus, to lower accuracy. What is more, it can be impossible to obtain a control surface with the uniform static pressure along the border during an experiment in the closed wind tunnel test section configuration. Hence, a pressure reconstruction algorithm was implemented for drag computations. Static pressure at every point with a known velocity value can be evaluated with Bernoulli's formula for incompressible flows. For every investigated velocity field, a global parameter of total pressure is calculated with known values of global airflow parameters (ambient pressure, air density, and freestream velocity). In the case of CFD, all values are known beforehand, whereas for the experimental session, it should be measured in the laboratory additionally. Afterwards, a static pressure value is reconstructed at every sample point along profiles I and II, except for the wake. In the wake region, Bernoulli's equation is not valid as the total pressure changes. The pressure reconstruction procedure is replaced with an interpolation from the field outside the wake.

Wake method
In the simplified wake method, a pressure analysis is disregarded and only a downstream velocity distribution profile is taken into account. To compute the drag force, the velocity profile is integrated only in the wake region where the velocity deficit can be observed. According to Eq. (10), the integrand term achieves a null value if the velocity equals the freestream value at the discrete position under analysis. As a result, it is justified to perform integration only for the limited portion of the velocity distribution profile. The integration limits are set to the velocity deficit region only and the remaining part of data is rejected even though a minor deviation from the freestream velocity can be observed. It should not be treated as an oversimplification since velocity variation in the remaining portion of the profile can be considered as negligible (for most cases, the deviation is not higher than 0.5-1% compared with the velocity drop in the wake which is larger by an order of magnitude).
To perform the integration within the desired limits, an automatic velocity deficit extraction algorithm was developed. The procedure can be divided into a few stages, namely: (1) Definition of the velocity distribution u(y) data at a specified position downstream (identically as in the control surface approach); (2) Calculation of the local curve slope with a central difference approach; (3) Velocity deficit detection by filtering the velocity data based on the specified threshold of calculated slope (increment of the slope by an order of magnitude with respect to the mean value obtained for the velocity profile outside the wake); (4) Removing the extracted velocity deficit from the original set of velocity values; (5) Interpolation of the data into the obtained gap based on the remaining part of the undisturbed velocity profile; 6) Computing the area integral of the velocity deficit proportional to the drag force coefficient (interpreted graphically by finding the area between the velocity deficit and the interpolated velocity curves marked in green in Fig. 9).

Control surface size independence test
Similarly as in the case of the lift coefficient, an effect of the control surface size on the drag estimation can be observed. In order to quantify it, a series of tests were performed. For illustration of the control surface with characteristic dimensions, see Fig. 7. It concerned a variation of the following parameters: -control surface length l cs (with a constant control surface height h cs ∕c = 0.5), The tests included drag coefficient computations with three approaches introduced before, i.e. the momentum change method with a pressure analysis, the momentum change method on a constant pressure assumption and the wake velocity deficit method (notation used on the plots: C D mom (u, p) , C D mom (u) and C D wake for each case, respectively).
In Figures 10 and 11, results of the analysis are presented. One representative CFD velocity field obtained for the angle of attack = 10 • was used in the tests. All drag coefficients are presented in a relative form with reference to the drag coefficient value obtained in the CFD simulations. For the supplementary test with a square-shaped control surface, only the results of the momentum method with a pressure analysis are shown. The simplified momentum method with a uniform pressure distribution was not displayed as it did not provide reliable outcomes.
For the presented angle of attack results, all drag coefficient values are underestimated with respect to the reference CFD result (not a case for each AoA). For the majority of plots, an asymptotic character is observed with an increase in the particular size parameter. All plotted results in Figs. 10 and 11 were well fitted to an exponential asymptotic model of the form f (x) = A ⋅ e −Bx + C. The simplified momentum method C D mom (u) did not yield reliable results. The coefficient values are highly underestimated and significantly sensitive to the control surface size. A rapid drop in the drag coefficient with an increase in the control surface height h cs was observed (non-physical negative values of the coefficient obtained). The convergence of the drag coefficient was obtained only at the large control surface length l cs ∕c > 4 , however, at a highly underestimated level. No convergence was achieved both for the varying profiles I and II position up to the investigated range 3.5c.
The extended momentum method C D mom (u, p) is less sensitive to the variation in size parameters. Asymptotic trends were achieved for all plots within the displayed range. Drag coefficient fluctuations below a rigorous level of 1% were achieved for: -relative control surface length l cs ∕c ≥ 6, -relative control surface height h cs ∕c ≥ 6, -relative profile I position d I ∕c ≥ 2, -relative profile II position d II ∕c ≥ 2.5.
The wake method C D wake provided very similar results as in the case of the momentum method, including both the velocity and pressure analysis. The convergence of the drag coefficient plot ( C D variation less than 1% ) was attained for the downstream profile position d II ∕c ≥ 2 , which is a more favourable result.
The final test, in which the control surface length and height were altered simultaneously, provided very promising results in the case of the extended momentum method (refer to Fig. 11). A convergence was achieved for the square control surface of the side length as little as 2.5 c.
To conclude, the simplified momentum method provides relatively good results for control surfaces of large length and small height only. Such a shape of the control surface minimizes the error following from the fact of neglecting the pressure integration in the equation used. For this method, it is favourable to position profiles I and II far away from the aerofoil, where its influence is negligible and an assumption of constant static pressure is more realistic. Secondly, for a small profile height, the velocity distribution at the downstream position is limited mainly to the wake region. This means that a contribution of the velocity drop to the fluid momentum change is more significant than the effects due to pressure. With an increase in the profile height, the contribution of pressure becomes more meaningful. As this method seems to be valid only when the control surface is limited to the wake region, it is reasonable to use the wake method instead, which yields more accurate results with less effort and computation time. Despite the implemented simplifications, the wake method provides outstanding results. If the control surface method is to be applied, however, it is evident that a pressure analysis needs to be included. The results obtained with this approach are very consistent and discrepancies between the expected reference data are negligible even for small sizes of the control surface. Therefore, it can be successfully applied when the velocity field obtained experimentally is of a limited size with respect to the aerofoil chord length.

CFD results
C L characteristics After applying the trained models to all in-house CFD velocity fields, a lift coefficient characteristics was obtained. The circulation-based results compared to the CFD reference data are presented in Fig. 12. An exceptional agreement was achieved as regards both the qualitative and quantitative aspects. As shown in Fig. 12b, a relative difference is very inconsiderable. The maximal discrepancy of ca. 2% was achieved. A general dependency to underestimate the lift results with respect to the reference data can be observed for the pre-stall condition. A constant systematic error can be noticed. This can be explained by the inaccuracy of the batch calculation procedure. The value of the local maximum found within the surface plot depends on the resolution of the batch calculations. The observed trend changes around the critical angle ( = 11 • ÷ 12 • ), after which the circulation method begins to overestimate the lift values. This could be related to the effects of the turbulence model applied in CFD. After the static stall, a large separation occurs and more complex flow physics is harder to be modelled accurately.

C D characteristics
The performed size sensitivity tests yielded a set of clear guidelines as regards the minimal value of size parameters to be chosen for a particular method. Drag characteristics were computed for two Reynolds number cases analysed numerically. The extended momentum and the wake method results were used only with the most favourable size of the control surface and position of the downstream profile. The results are presented in Fig. 13. An influence of the Reynolds number can be observed. Slightly larger discrepancies occurred in the case of higher freestream airspeed with maximal differences of ca. 12% and 18% for the momentum and the wake method, respectively.
The plot with a relative difference shown in Fig. 13b reveals an interesting dependency. Drag results for the angle of attack in pre-stall (stall observed above = 12 • ) are underestimated with respect to the reference data. An absolute error value is low for small angles of attack and increases up to around = 7 • ) where drag values are underestimated most significantly. For larger angles of attack, the character of the drag estimation error changes. The drag characteristics in the post-stall region are overestimated with respect to the reference data. Such a change in the error trend is identical to the observations made in the case of lift results. Similarly, it needs to be related to the effects of the turbulence model applied in the numerical simulations and a more complicated physics of the flow when separation occurs.

PIV results
Algorithms for computations of aerodynamic forces, extensively tested and tuned with the reference CFD velocity fields, were applied to the PIV experimental data. In  PIV-based lift characteristics of the NACA0012 aerofoil are shown. The PIV results are compared to the CFD reference data. For easier evaluation of the relative difference between the experimental data and CFD, a margin of ±10% of CFDbased lift coefficient values was highlighted in green (a direct quantitative comparison was not always possible due to inconsistency in the investigated angles of attack). Such a level of agreement between the numerical and experimental results can be considered satisfactory.
The analysis of the obtained characteristics reveals an exceptional agreement. As shown in Fig. 14a, a very good coincidence with the reference curve was achieved for the experimental case at U ∞ = 11m/s. The majority of computed PIV lift results can be found within the defined margin of ±10% . Moreover, around half of the experimental points are characterized by a relative difference with respect to the reference CFD data no larger than ±5% . Slight discrepancies can be observed only for high angles of attack in the poststall regions. In the case of results for the higher flow speed U ∞ = 22 m/s shown in Fig. 14b, discrepancies are slightly larger, which can be especially observed for the post-stall part of the characteristics. The lift coefficient drop observed in the PIV results is less significant than in the CFD simulations. This indicates that during the experiment, separation (and, thus, a pressure drop and a lift decrease) was not formed as rapidly as in the CFD computations.
From the qualitative point of view, shapes of the lift characteristics are determined very accurately. The slope of the pre-stall portion of the curve is in a perfect agreement with the reference results. The concave character of the lift curve with an inflection point around the angle = 2 • and = −2 • predicted in the CFD simulations was confirmed by the PIV measurements. The position of the critical angle found in the PIV experiment fits the numerical findings with a reasonable tolerance of ±1 • for a lower Re case. It is harder to determine precisely the critical angle for a high Re experimental case at the positive AoA as the curve reaches a plateau. Nevertheless, the maximal value of the lift coefficient C L max was determined with very good accuracy. For lower Re relative differences with respect to the CFD, the reference value equalled ca. 4.5% and 3.5% for the positive and negative AoA, respectively. For the high Re case, C L max differed by only 1% and 2.5%.
The experimental results of the NACA0012 drag coefficient are shown in Fig. 15. The plotted PIV drag characteristics were obtained with the wake method only. The momentum method did not provide satisfactory results as the velocity reading from the upstream profile was not available at a distance large enough from the aerofoil. As an effect, the velocity distribution was highly affected by the presence of the model. A good agreement can be observed between the CFD and PIV results for moderate angles of attack in the pre-stall region. A more rapid drag increase at high angles can be observed for the experimental results. The symmetry of results is visible for both the WT flow speed cases in the pre-stall region. The expected influence of the Reynolds number on drag results is observed for the experimental data, i.e. a drag drop with an increase in the Reynolds number.
Detailed drag characteristic plots for a limited angle range are presented in Fig. 15b and c. The error margins of ±10% and ±20% with respect to the reference CFD data are highlighted. As shown, a satisfactory agreement of the results was achieved in the range of from −10 • to 10 • . Roughly, half of the results vary by no more than 10% from the reference data both for low and high Re cases.

Discussion and conclusion
An effective method for lift and drag determination in the form of dimensionless force coefficients was developed. The method allows for computing sectional loads as it utilizes two-dimensional velocity fields. The validity of algorithms was confirmed for numerical results of the flow around the NACA0012 aerofoil in the subsonic regime (i.e. Re = 0.7 × 10 5 and Re = 1.4 × 10 5 ). The obtained level of agreement between the reference CFD loads and the velocity-based loads computed from the numerical fields is exceptionally high. The maximal discrepancy of ca. 2% in the lift coefficient curves was achieved. In the case of the drag estimation, the maximal discrepancies were at the level of around 15%.
For the lift computation, a velocity circulation analysis was used taking advantage of the Kutta-Joukowski theorem, which defines the dependency between the lift force L and the circulation . The method appeared to be very efficient and accurate as well as quick and easily applicable at the same time (a much less demanding technique than the other ones used in the majority of investigations in this field). It requires only information on a velocity distribution along the closed-loop encircling the investigated body. No pressure data are required, which makes it a perfect tool to apply to velocity field results.
Although according to the theoretical model, the circulation value is independent of the size of the taken contour, special attention was paid to optimization of this parameter. The analysis performed is highly valuable and unique as none of many researchers working in this field (refer to Sect. 1.1) has presented any results regarding this aspect. The results of the contour size influence are shown herein in a very clarified form. 3D surface plots (as shown in Fig. 6) were used. The analysis of such plots not only allows one to find the optimal lift coefficient value but enables an evaluation of the PIV/CFD velocity field quality. In a correctly plotted surface of lift coefficient values, a plateau should be easily identified as it has an asymptotic character. Every discontinuity, cavity, peak, irregularity, etc., in the given surface indicates an aberration (e.g. measurement noise, artefacts, null vectors) in the velocity distribution along a particular contour size. The results of the performed investigations revealed the best correlation with the reference data when the maximal value of the lift coefficient from the plateau portion was taken.
Two effective methods for drag determination were developed. The first one exploits the equation of momentum conservation and the control surface approach with two control profiles upstream and downstream of the body. The second is a simplified one, in which pressure terms are omitted and the upstream velocity distribution is assumed as uniform. Therefore, the velocity distribution in the aerofoil wake is analysed only. As indicated in the verification procedure with the CFD reference data, the momentum method appeared to be more accurate with a maximal error at the level of ca. 10% vs 15% achieved for the wake approach. Lower accuracy obtained with the wake method calculations follows from simplifications used in the final derivation formula of the theoretical model (see Sect. 2.1).
The effects of the control surface size on drag results were extensively tested. The results presented in Sect. 3.3 indicated a significant influence of size parameters. An analysis of the converging plots yielded a set of recommendations on the size parameters, which provided acceptable and optimal results within each drag estimation method. A satisfactory level of the results convergence is such at which drag coefficient variations are less than 1%. Crucial findings as regards the momentum method are such that the optimal shape of the control surface is a square with the minimal side length 2.5 times larger than the aerofoil chord length, whereas in the case of the wake method, the minimal downstream profile distance from the trailing edge should be larger than two chord lengths d II ∕c ≥ 2.
Drag computations with the momentum method require a static pressure integration over the control surface borders as shown in Eq. (11). As only the flow velocity is available, pressure had to be reconstructed on the basis of Bernoulli's equation. The approximation yielded satisfactory results.
An alternative version of the momentum method was tested. A fluid momentum change between the upstream and downstream profile was analysed without the pressure term (i.e. p = const-see Eq. (7)). In general, this approach provided unsatisfactory results. It appeared to be highly sensitive to the control surface size and gave acceptable results only for a specific condition (i.e. a large distance of profiles from the aerofoil and a low height of the control surface comparable with the wake height).
The final validation of the developed algorithms was carried out by an application of own experimental PIV velocity fields. An exceptional agreement between the velocity-based lift coefficient characteristics and the reference data was found. The shape of curves agrees perfectly (i.e. the slope of the linear part of the lift curve, the angle of static stall phenomena, the symmetry of the plots). From the quantitative point of view, the maximal discrepancies observed for all measurement points are no larger than ±10% , whereas the majority of results fit an error margin of ±5% . Similar accuracy can be observed for both the Reynolds numbers investigated. The error of the maximal lift coefficient value C L max estimation is highly satisfactory as well − 4.5% and 2.5% for low and high Re cases, respectively. A similar level of accuracy of the PIV-based lift characteristics was achieved in other investigations enclosed in Sect. 1.1, e.g. by van Oudheusden et al. (2006).
For drag estimation, a good level of coincidence with the reference data was obtained with differences no larger than ±20% . Although this can be concerned as a significant estimation error, the lift to drag ratio should be taken into account. The drag value for aerofoils is usually an order of magnitude lower than the lift force. As an effect, drag does not contribute significantly to the resultant aerodynamic force exerted on the aerofoil and its estimation error can be regarded as negligible.
An important issue is a reference data quality used for validating the developed methods of load determination. Due to a lack of the reliable experimental force measurement technique to be applied within this study, CFD simulations were used. The reason of the observed discrepancies between the attained CFD and PIV characteristics can be identified as uncertainty of the numerical results. The CFD results are expected to be less accurate in the post-stall condition, which is due to complex flow physics characterized by high-pressure gradients and large vortex structures changing in time. Additionally, a low Reynolds flow regime condition investigated imposes the need of simulating laminar-turbulence transition phenomena that can influence the solution accuracy and stability. What is more, some subtleties of an experimental model such as aerofoil surface roughness, flow instabilities have not been included in the numerical model. All of the above issues can be a potential source of the observed discrepancies between the numerical and experimental results at high angles of attack. A similar tendency was observed when the in-house numerical aerodynamic characteristics were compared with the literature experimental data. Taking this into consideration, CFD simulations are not a fully reliable source for validation of PIVbased lift results in the whole range of angles of attack under investigation. From the practical point of view, however, aerofoils are dedicated to operate under pre-stall conditions. As long as the linear portion of the lift curve is regarded, it is proven that the developed algorithms provide exceptional good performance. Th obtained results are characterized by very good accuracy, repeatability, and agreement with the reference CFD results.
Although satisfactory results were achieved, the fact of testing one aerofoil family only in a limited range of Reynolds numbers can be a questionable issue. To state the validity of algorithms under more general conditions, it could be vital to extend the experimental campaign on different aerofoil types (e.g. asymmetrical), tested preferably in higher Reynolds number flow regimes. The tests performed within the frame of these investigations, however, were adjusted to the IMP wind tunnel facility capability. The section size and the maximal achievable flow speed disable higher Reynolds number flows while keeping the aerofoil size at a reasonable level. The motivation behind choosing a simple NACA0012 aerofoil design was an easy access to a huge database of the experimental aerodynamic characteristics provided by various research centres and analysed comprehensively in McCroskey (1987).
Due to the development of a methodology providing algorithms for accurate aerodynamic coefficients determination, a concept of noninvasive forces measurements was proven to be valid. There are a few advantages of such a measurement technique in terms of wind tunnel testing applications, namely: -measurement of forces at the experimental stand, where the installation of standard measurement devices (balance system, pressure probes) is impossible due to difficult access, limited space, etc., -avoiding flow disturbances imposed by the measurement probes placed close to the tested object, which protects from influencing measured values, -measurement of aerodynamic forces exerted on rotating bodies, e.g. wind turbine blades, propellers, etc., -determination of sectional loads impossible in the case of balance measurements, which provide the total force applied to a wing, blade, etc.
distortions (due to lens or inclination of cameras) by mapping the final image onto the 2D plane (Fig. 16).

Details of the CFD data resolution
In the case of the CFD results, the density of data samples is determined in a process of mesh generation. The computational mesh used for the simulations was composed of a structured mesh with regular hexahedral elements and an unstructured mesh with tetrahedral and wedge elements of various size and shape. The size of elements was controlled by several sizing functions, thus the element edge length varied from 5mm in the far flow field through 1mm in the wake region downstream of the aerofoil up to 0.02mm in the boundary layer. In the script used for loads determination, non-uniformly distributed CFD node data were firstly mapped onto a regular grid similar to the one available in the PIV results. Such an approach allowed one to simplify algorithms and apply the identical methodology for data processing. A decision had to be made on a resolution of the new computational grid used in Matlab after resampling. An intuitive solution was to choose a resolution higher than that of the source data to avoid information loss. On the other hand, a higher resolution leads to a larger data file size consuming more disk space, RAM memory and computation time. In order to find the golden ratio, a grid size independence test was performed for different resolutions of the main grid varying from 0.1 to 10 mm. Certain calculations performed within the developed Matlab script (e.g. mapping of the velocity data onto the closed loop encircling the aerofoil used to determine circulation) require another interpolation. Thus, the independence test was performed for the step size of the numerical integration as well. The methodology of finding the optimal resolution consisted in performing lift coefficient calculations within the developed Matlab code for various resolutions of the main computational domain and the integration step size. Two tests were carried out independently. The main grid resolution independence test was performed for a constant value of the integration step size equal to 0.5mm. Analogically, various integration step sizes were investigated for a constant main grid resolution of 0.5 mm. In the figure below, a chart with results of this analysis is presented. The calculations were performed for an exemplary reference velocity vector field obtained from the CFD simulation of the NACA0012 aerofoil at the angle of attack = 7 • . The abscissa is the resolution in mm, the ordinate is the normalised lift coefficient. The lowest resolution (0.1 mm for the main grid and 0.05 mm for the integration step) was assumed as the reference value of the lift coefficient. The source of differences can be identified as a numerical error due to an approximation of data during the interpolation and the numerical integration method used (trapezoidal integration in this case). The convergence of normalised results is clearly visible for resolution values below 1 mm. Values of the lift coefficient obtained for the resolution equal to 1 mm differ from the reference values by ca. 0.009% and 0.1% for the computational grid