Numerical simulations of the flow with the prescribed displacement of the airfoil and comparison with experiment

The work is devoted to comparing measured data with the results of numerical simulations. As mathematical model was used mathematical model whitout turbulence for incompressible flow In the experiment was observed the behavior of designed NACA0015 airfoil in airflow. For the numerical solution was used OpenFOAM computational package, this is open-source software based on finite volume method. In the numerical solution is prescribed displacement of the airfoil, which corresponds to the experiment. The velocity at a point close to the airfoil surface is compared with the experimental data obtained from interferographic measurements of the velocity field. Numerical solution is computed on a 3D mesh composed of about 1 million ortogonal hexahedron elements. The time step is limited by the Courant number. Parallel computations are run on supercomputers of the CIV at Technical University in Prague (HAL and FOX) and on a computer cluster of the Faculty of Mechatronics of Liberec (HYDRA). Run time is fixed at five periods, the results from the fifth periods and average value for all periods are then be compared with experiment.


Introduction
The experimental methods in aeroelasticity have since long been complemented by attempts to design mathematical models and solve the numerical approximations using computers.The classical methods used in aerospace engineering include simplified (potential) flow models coupled with linear spring-mass systems [1][2][3].However, in more complex flow regimes, e. g. stall flutter with massive flow separation, these models cannot predict the aeroelastic response of the system and the stability boundaries accurately, and the airflow has to be modeled by Navier-Stokes equations.Turbulent flow has a three-dimensional character and so 3D simulations are preferable, but in a number of cases, two-dimensional models are still applied from practical reasons (mainly because the computational cost is drastically lower).
In turbulent flow simulations, three approaches are possible: direct numerical simulations (DNS), where the Navier-Stokes equations are discretized and solved directly and all turbulent scales are resolved.In many practical cases, e.g.[4], the Reynolds-averaged Navier-Stokes (RANS) approach is applied, where the Reynolds stresses can be modeled by a vast variety of turbulence models.
Because RANS models often give incorrect results (especially for unsteady and separated airflow), "laminar" simulations (actually DNS on an insufficiently fine mesh) are still widely used [5].A compromise between DNS and RANS are the Large Eddy Simulations (LES), where large coherent turbulent structures are resolved and the small-scale isotropic turbulence is modeled by means of a sub-grid scale model.
Among computational studies on flow-induced airfoil vibrations and flutter instabilities, a 2-DOF airfoil moving in both pitching and plunging mode was studied numerically for transonic flow by Geissler [6] using on a 2D Navier-Stokes equations solver and the Spalart-Allmaras turbulence model.A numerical investigation of the 2-DOF bending/torsion flutter characteristics of an airfoil in 2D transonic flow was carried out by Weber et al. [7] using RANS equations and Baldwin-Lomax turbulence model.The paper by Wang and Zha [8] investigates the NLR7301 airfoil limit cycle oscillation (LCO) in transonic flow caused by the flow nonlinearity using Detached Eddy Simulation (DES) of turbulence.The 2-DOF airfoil with freeplay nonlinearity in pitch was investigated numerically by Zhao et al. [9] for low, intermediate and high level of turbulence.Poirel et al. [10,11]  50,000 and 150,000.Both laminar and unsteady RANS calculations using the SST k-omega model with a low-Reynolds-number correction have been performed and found to produce reasonably accurate limit cycle pitching oscillations.
In case of this paper is carried out a parallel numerical solution of incompressible airflow on 3D meshes.This paper is focused on the comparison of experimental data obtained during measurements in aerodynamic tunnel in Nový Knín.In the experiments in the tunnel was measured the vibrations of airfoil NACA0015.Experimental data from the measurements were processed using program in Matlab for evaluation Mach-Zehnder interferograms.The result of the work is comparison of the experimental data with the results of numerical simulations and evaluates their agreement depending on the distribution of Mach number on the surface of the wing.

Mathematical model, geometry and boundary conditions
Flow around the airfoil can be described using the equation of continuity and N-S equations.In this case are selected equations for incompressible and unsteady flow.In a strong conservative form that is suitable for finite volume discretization where u, p and ρ are fluid velocity, pressure, and density, respectively is kinematic viscosity.The geometry is derived from the experimental setup, in which the measurements were carried out at the airfoil NACA0015 (figure 1).The airfoil is placed on elastic support allowing two degrees of freedom -pitch (rotation about the elastic axis, located in 1/3 of the chord length) and plunge (vertical motion).The model was placed in the wind tunnel of the Institute of Thermomechanics in Nový Knín.The measurement setup allowed to record vibration data, together with the flow fields measured by Mach-Zehnder interferometer and high-speed camera (more details can be found in [12,13].The measured value of the frequency of vibration was 19.5 Hz, the oscillation period is 51.3 ms.The amplitude of the rotational movement was ± 17°, and the amplitude of the plunging movement was ± 7 mm.Size of the computational domain is 580 210 80 mm. where t is the ratio of the greatest profile thickness to chord length, c is the chord length, x is the position of a point on the continuous profile along the chord from 0 to c, y is half the thickness profile for a given value of x.Following boundary conditions were prescribed: on input (Γ inlet ) is prescribed velocity u = 147 m / s.On the walls of channel (Γ top , Γ bottom and Γ frontandback ) and the surface of the wing (Γ wing ) is prescribed velocity u = 0 m / s.Due to large intensity of vorticity resulting from flow separation downstream of the wing surface a stabilized boundary condition is prescribed at outlet Γ out : / = 0, when velocity direction points outward of the domain, u = 0 m / s otherwise.
Boundary conditions for the pressure at the outlet (Γ out ) prescribe p = 0 Pa.All boundaries walls (Γ inlet ,Γ top , Γ bottom and Γ frontandback ) have the Neumann boundary condition ∂ p / ∂ n = 0.

Mesh and numeric solution
For numerical solution was used mesh with hexahedral elements.Number of elements in the mesh is 879124 hexahedral, 13140 polyhedral and 2688 prisms elements.The elements near the surface of the wing are refined to capture better the boundary layer.On the figure 2 is visible 2D slice of 3D mesh.During calculation, the mesh changes due to wing vibration.During this movement, the elements of the mesh deform.Changing the position of the node elements is the solution of the Laplace equation for the mesh velocity w with spatially variable diffusivity γ: For the calculation of the flow around the oscillating airfoil was used computational package OpenFOAM [14].It is an open-source library written in C++ language.The basis of this library is a finite volume method for structured and unstructured mesh.
The numerical scheme employs the 'pimple' algorithm, which is a combination of the standard SIMPLE and PISO algorithms.This algorithm is a combination of SIMPLE and PISO algorithm.Combines the advantages of these algorithms and allows thus greater time step.
The resulting linear system for the momentum was solved by bi-conjugate gradient method (PBiCG) with diagonal-based incomplete LU (DILU) preconditioning.For the pressure predictor and corrector steps, faster convergence was obtained using a geometric-algebraic multigrid (GAMG).

Results
The simulation time was 5 periods of vibration (256,5 ms), solutions was compared for fifth period and for average for all periods.The movement of the wing is the in numerical solution is prescribed according to the measured data.For comparison of the experimental data with the numerical solution it was necessary to choose the points at which the calculated data are compared with experimental data [14].
Points must be close on wing surface, because evaluation experiment is made using by interferographic method on the surface of the wing.For the numerical solution is prescribed zero flow velocity on the wing (Γ wing ).For this reason mesh nodes in the first and second layer were chosen (M1-M4, see figure 3).The experimental data is compared against the fifth period of vibration and against an average over all five periods.The results are shown in graphs, where on the x-axis is plotted position normalized to the chord length and on the y-axis the Mach number.
How you can see on graph (figures.5, 8 and 11), the numerical results match the experimental data reasonably on the bottom surface, where the flow is not separated.The numerical solution shows significantly lower Mach numbers in the stagnation zone, which is most likely caused by insufficient resolution of the interferometric images in this region.Better match is provided at the points in the second layer (points M2 and M4).On the top surface of wing (figures 6, 9, 12 and 15) the numerical results match the experimental data much worse.This can be caused by two effects: first, the flow patterns in the separated region are very complex and the interferographic images and the subsequent evaluation procedure might not capture the separated turbulent boundary layer and the real velocity field with sufficient resolution.Second, the current simulation is run without a turbulence model and with insufficiently fine mesh in the boundary layer.Thus, its predictions in the separation region have to be taken with caution.

Conclusion and discussion
A 3D finite volume model of airflow on a dynamic mesh is very complicate for numerical solution.The solution is time consuming.In the new time steps must by solution equation computing deformation of mesh.The solution 5 period oscillation was computed 156,5 hours on 26 processors Intel Xeon 2,66 GHz (supercomputer fox Technical University in Prague).In this case was used simple mathematical model without compressibility and turbulence.
The result shows that this model provides approximate match of the numerical solution and experimental data only in the regions, where there is no flow separation.In the separated regions, the results differ significantly.This is likely caused by the fact that for the purpose of this preliminary study, a simple incompressible model without any turbulence modelling was chosen.As the geometry is fully 3D, it is also highly problematic to have a sufficiently fine mesh when no massively parallel computational resources are available.

01080-p.5
The evaluation of the experimental flow velocity fields from the interferographic images, on the other hand, is also problematic, especially in the regions of high density gradients.
In the future, it will be appropriate to switch to the compressible flow model and possibly incorporate a suitable turbulence model.The validity of the experimental velocity values evaluated in the separated region should be assessed carefully, too.

Fig. 5 .
Fig. 5. Mach number in phase 002 on the bottom surface of the wing.

Fig. 6 .
Fig. 6.Mach number field in phase 002 on the top surface of the wing.

Fig. 11 .
Fig. 11.Mach number in phase 006 on the bottom surface of the wing.

Fig. 12 .Fig. 13 .
Fig. 12. Mach number in phase 006 on the top surface of the wing.

Fig. 14 .
Fig. 14.Mach number in phase 010 on the bottom surface of the wing.

Fig. 15 .
Fig. 15.Mach number in phase 010 on the top surface of the wing..