Parallel numerical simulation of oscillating airfoil NACA 0015 in the channel due to flutter instability

The work is devoted to 3D and 2D parallel numerical computation of pressure and velocity fields around an elastically supported airfoil self-oscillating due to interaction with the airflow. Numerical solution is computed in the OpenFOAM package, an open-source software package based on finite volume method. Movement of airfoil is described by translation and rotation, identified from experimental data. A new boundary condition for the 2DOF motion of the airfoil was implemented. The results of numerical simulations (velocity) are compared with data measured in a wind tunnel, where a physical model of NACA0015 airfoil was mounted and tuned to exhibit the flutter instability. The experimental results were obtained previously in the Institute of Thermomechanics by interferographic measurements in a subsonic wind tunnel in Nový Knín.


Introduction
The interaction of elastic structure and fluids appears frequently in fields such as aerospace engineering, biomechanics or turbine design.In aerospace applications, the classical problem is the interaction of airflow with the elastically supported airflow.The airfoil can be approximated as a two-degree-of-freedom system, with vertical translation and rotation modes.During flutter instability, this movement is harmonic with constant (or increasing) amplitude, and the energy dissipated by the internal damping is compensated by energy transfer from the airflow [1][2][3].
The numerical solution of interaction of fluids and elastic structures are very often time consuming.Flow around an oscillating airfiol is complex and during large angles of -attack massive flow separation may occur.In aerospace engineering applications, flow is almost always turbulent.Two-dimensional models are still applied from practical reasons, mainly because the computational cost is drastically lower.
In the CFD for turbulent flow simulation is possible to use the following approaches [4].First is the Reynolds-Averaged Navier-Stokes (RANS) approach where the Reynolds stresses can be modelled by a vast variety of turbulent models.However, in many technical applications, especially when the airflow is separated, RANS models give incorrect results.The second possibility is the Direct Numerical Simulations (DNS), which solve directly the Navier-Stokes equations.Mesh must be fine enough, element size corresponds to the smallest scales of turbulence.DNS is very time consuming, and in most cases of practical interest, it is still not feasible.The number of elements scales with Re 9/4 and for explicit solvers, time steps are also very short (Courant number limitation).Third approach is the Large Eddy Simulation (LES), where the solutions of the Navier-Stokes equations are divided into two parts.Large coherent turbulent structures are solved directly and the small-scale isotropic turbulence is modelled 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 [5] 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. [10] using RANS equations and Baldwin-Lomax turbulence model.The paper by Wang and Zha [9] 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. [11] for low, intermediate and high level of turbulence.Poirel et al. [6] studied the low amplitude self-sustained pitch airfoil oscillations in incompressible flow by 2D numerical simulations for Reynolds numbers between 50,000 and 150,000.Both laminar and unsteady RANS calculations using the SST k-omega model with a low-Reynoldsnumber correction have been performed and found to produce reasonably accurate limit cycle pitching The current paper is focused on a parallel numerical simulation of oscillating airfoil in the channel.Results of numerical simulation are compared with experimental data.Specifically, the simulated distribution of velocity field on the surface of the wing vibrating in the channel due to flutter instability is compared to the experimental values.In the experiment, the velocity is evaluated from interferograms [8].

Geometry and mesh
The computational domain for numerical solution corresponds to the experimental setup.The airfoil is located in an experimental setup in the middle of the height of the channel 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).Size of the computational domain is 580x210x80 mm (figure 1).

Fig. 1. Geometry of computational domain.
Numerical solutions are computed on 3D and 2D mesh.Mesh in 3D case has 890 k elements and in 2D case has 139 k elements.During the simulation, the airfoil oscillates and the mesh changes shape, but the number of elements stays the same.Changing the position of the element nodes is solved of the Laplace equation for the mesh velocity w with spatially variable diffusivity γ: Different between 3D and 2D mesh is the arrangement of elements on the surface of the wing.In 3D case the elements follow the shape of wings (figure 2).

Mathematical model numeric solution and boundary conditions
Numerical simulation of flow around the airfoil can be described using the incompressible equation of continuity and incompressible N-S equations. (3 where u, p and ρ are fluid velocity, pressure, and density, respectively ߥ is kinematic viscosity.Movement of the airfoil is prescribed by boundary condition for the mesh position on boundary Γ wing .This boundary condition, used for prescription of movement with two degrees of freedom, was designed and implemented to OpenFOAM.The parameters are frequency, amplitude, and phase shift of translation and rotation modes.The kinematic data was identified from wind tunnel measurements, where for this particular case the frequency was 19.5 Hz, amplitude of the rotational movement was ± 17 °, and the amplitude of the plunging movement was ± 7 mm and phase shift between translation and rotation is 8 ° [7]. Following boundary conditions for velocity were prescribed: on input (Γ inlet ) is prescribed velocity u x = 147 m/s.On the walls of channel (Γ top , Γ bottom and Γ frontandback ) is prescribed velocity u = 0 m/s.On the surface of the wing (Γ wing ) is prescribed velocity which corresponds to local velocity of the airfoil surface.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.All walls (Γ inlet, Γ top , Γ bottom and Γ frontandback ) have the Neumann boundary condition ∂p / ∂n = 0.
Numerical simulation of the flow around the oscillating airfoil was solved in parallel in computational package OpenFOAM.The basis of this library is a finite volume method.As solver was used a combination of PISO and SIMPLE algorithms.

Results
The results of numerical simulations are averaged value for 3-5 periods of vibration (simulation time is 256.5 ms).The movement of wing is prescribed according to the measured data.Values of numerical solutions are shown in graphs, where on the x-axis is plotted position normalized to the chord length and on the y-axis the velocity.
The results of the velocity field distribution around the airfoil are not as good as in the case of pressures [12].Velocity in the numerical simulation is evaluated in the nodes of mesh, which are located two cells from the surface of the wing.The velocity field on the wing surface is significantly affected by the boundary conditions and the behaviour of the flow in the boundary layer.
In figures 3, 6, 8 and 10 are shown velocity field fields from 3D simulation in 2D sections in the middle of the computational domain.On graphs are shown result on the bottom surface of the wing and one case on the top surface (figure 5).On the top surface of the wing the numerical results match the experimental data much worse.The reason for this is: 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.
The velocity fields in figures 4, 6, 8 and 10 show that on the trailing edge is a big difference between results of 2D and 3D simulations.The first reason is the difference between the two meshes, where the 3D mesh includes a body-fitted layer better capturing the trailing edge.The second reason might be the difference between the physics of the real 3D flow and the 2D approximation, which disregards convection and diffusion in the zdirection.

Conclusion and discussion
The results of a 2D and 3D numerical model of flow past an oscillating airfoil have been compared with experimental data.The results show that the 2D approximation, which has clearly much lower computational demands, gives reasonable velocity distribution predictions around most of the airfoil surface with attached flow, except in the vicinity of the trailing edge.Thus it seems that a 2D model could be used as a computationally cheap alternative to the fully 3D model to get a fast first prediction of the flow field.
The numerical results in the case, where the flow is massively separated from the surface, depart from the experimental data significantly.However, it should be noted that in this case even the experimental values have to be taken with care: the airflow in the separation bubble is heavily turbulent, and the velocities are considerably lower than the free-stream velocity, which negatively influences the resolution of the interferometric method.
In future, it will be appropriate to use a compressible flow model for the cases when the free-stream velocity exceeds Mach numbers of about 0.3, and to incorporate a suitable turbulence model.First preliminary results are currently being finished in this aspect.

Fig. 4 .
Fig. 4. Velocity field in phase 002 on the bottom surface of the wing.

Fig. 5 .
Fig. 5. Velocity field in phase 002 on the top surface of the wing.

Fig. 7 .
Fig. 7. Velocity field in phase 004 on the bottom surface of the wing.

Fig. 9 .
Fig. 9. Velocity field in phase 006 on the bottom surface of the wing.