COMPUTER SIMULATION OF FLUTTER ELEMENTS OF A SYMMETRIC AIRFOIL USING THE VIRTUAL BOUNDARY METHOD FOR FLUID-STRUCTURE INTERACTION

Abstract. The objective of the present work is an numerical analysis of flutter elements of a typical airfoil using a classical method with direct eigenvalue method for the undamped case (forced oscillation). The dynamic coupling of the airfoil structure into two-dimensional aerodynamic flow was simulated with ultra-low Reynolds number so that the aeroelastic motion of the airfoil in the flow could have simulated in the time domain. The intention was to determine the flutter elements in the small wings (like a Micro air vehicle), using a twodimensional aerodynamic code based on Virtual Boundary Method, suitably coupled to the airfoil structural characteristics. The symmetric NACA 0012 airfoil was choosed to simulate the time history of flow parameters. The airfoil was considered as a rigid section, supported by translational and rotational springs, so that only heave and pitch degrees of freedom are permitted at the point of support. The pitching and heaving movements where simulated separately. The effects of forced oscillation were analyzed an attack angle of 0 f r heaving, and vertical oscillation of±0.1 of chord length; pitching of ±2 with frequency (sinusoidal), f s , of 1Hz, 2Hz, 5Hz and 10Hz. The effects of the reduced frequency (k), amplitude of forced oscillation (h) and the maximum non-dimensional flapping velocity (kh) on the thrust generation were analyzed. The pressure coefficient CP, lift L, lift coefficient CL, drag coefficient CD and pitching moment M about the support point were computed. The results obtained were compared and agree with the literature ones.


INTRODUCTION
The problem of calculating the response of a system subject to an external excitation (forced oscillation), since the mathematical model and of the structure properties (natural frequencies, damping factors and natural modes of vibration) from the input and output information, in general is a difficult problem to be solved and form part of the class of inverse problems.Among the methods that treat these problems, the analysis modal methods are characterized by direct calculation of the modal parameters of the structure and, computationally, these parameters can be determined directly from the simulation results, in the time domain or frequency domain.In the study of aerodynamics, aeroelastic effects are the interactions results of the distributed dynamic fluid loads and inertial, structural efforts and elastic reactions that can induce rupture behavior of the structure, including flutter.The phenomenon of flutter is usually observed in the wings and control surfaces, as these types of structures are subject to large distributed dynamic fluid loads produced from the dynamic deviations of a elastic structure in the undeformed state [1].The structural dynamic behavior of an aerolastic model can be described by their modal parameters considering three hypotheses: a) linearity of the dynamic behavior: the response of the structure to a combination of forces applied simultaneously is equivalent to the sum of responses of each force acting individually; b) Invariance in time: the physical parameters of the structure are constant; c) Observability: the relationship between input/output measurement contains sufficient information to determine the dynamic behavior of the model.
The dynamic system to be studied, shown in the model in Fig. 1, consists of a airfoil profile NACA Series 0012 subject to forced movement of vertical oscillation (plunging) and torsional (pitching).The study of this three-dimensional phenomenon reduced to a bidimensional aerodynamic flow is justified by considering a hypothetical wing of infinite extent such that the properties of the wing remain the same in any arbitrary section located along the length, where c is the the length of airfoil chord , x cp is the position of pressure center, x cm is the position of mass center, O is the position of the center of the twist axis (pivot) and h and α denote the plunge and pitch displacements, respectively.The integration of the distributed dynamic fluid loads are plunging, L (positive upward) and pitching moment M (positive clockwise -nose up) acting at the torsional point O.The moment of inertia of the section around the axis of torsion, pivot O, is given by I α .The rigidities of the lift and torsional springs are indicated by k h and k α and the structural damping coefficients are c h and c α , if the damping was considered.
Figure 1.Dynamic model with two degrees of freedom, [1]. the pitch angle is α and α e is the effective attack angle.Basically flutter is defined by a critical velocity V and a critical frequency ω f , where V is the velocity flow non-disturbed that focuses on the structure and ω f is the natural frequency of vibration (harmonic oscillations simple) of a given structure immersed in a flow under certain conditions of plunging.The solution of the vibration leads to a complex eigenvalue problem where two characteristic numbers determine the velocity and frequency.Few studies for this range of Reynolds number are found in the literature, for example, [2] analyzes oscillation with large amplitudes or, in case of flutter, [3] for Reynolds numbers above 10 3 .

FORMULATION AND NUMERICAL METHOD
The distributed dynamic fluid loads are distributed through the analysis of movements of pitch and elevation separately and, when coupled, these define the kinematic angle of attack or effective α e (t) = α α − α h .As [3], to a small Strouhal number (St ≤ 0.25), the contribution of the plunge motion to the kinematic attack angle is small and it follows the form of a sinusoidal function.The figure Fig. 2 show the plunge and pitch movements.We used a) b) Figure 2. Plunge, 0.1c, and pitch, ±2 o , movements.
the immersed boundaries methodology [4], particularly the virtual boundary method [5] for modeling fluid-structure interaction in a domain discretized by two meshes: one fixed mesh to represent the two-dimensional flow field (Eulerian) and one mesh of points to represent the immersed boundary (Lagrangian) that is independent and does not need to align with the mesh Eulerian, which enables simulation of flow over complex geometries [6] and/or local refinement of the flow.

Governing equations
The flow is modeled by Navier-Stokes equations with forcing term, mesh discretization in Cartesian and computational fluid-structure interaction in the model proposed by [4].Consider a flow of homogeneous and viscous fluid in a two-dimensional rectangular domain Ω = [0, L] × [0, L] with an immersed boundary represented by a closed curve Γ, described by X(s, t), with 0 ≤ s ≤ L b and with X(0, t) = X(L b , t), where L b is the length of the curve Γ.The governing equations are given by: In the equations Eq. (1) to Eq. ( 4), x = (x, y) is the position vector, u(x, t) = (u(x, t), v(x, t)) is the field of fluid velocity, p(x, t) is the pressure field and ρ(x, t) is the specific mass.The force acting on the flow (relative to dx = dxdy) is f(x, t) = (f 1 (x, t), f 2 (x, t)), while the force exerted on the immersed boundary (compared to ds) is F(s, t) = (F 1 (s, t), F 2 (s, t)).The equation Eq. ( 1) are the Navier-Stokes equations for incompressible flow and Eq. ( 2) is the continuity equation.The equations Eq. ( 3) and Eq. ( 4) represent the interaction between the flow and the immersed boundary.The Dirac delta function in the equations is a function composed of two other delta functions, δ 2 (x) = δ(x)δ(y).One "feedback" function ensures that the flow velocity is zero in the points defining the immersed boundary with the condition no-slip and can be expressed as where the constants negative α and β will be chosen with large enough magnitude to adjust the fluid velocity and the velocity near of the interface in order to obtain the expected physical behavior of the flow, δ 2 (x − x s ) is the Dirac delta function given by equation Eq. ( 8), X(s) are the lagrangian points arranged on the immersed boundary, F(X(s), t) is the Lagrangian force density (nonzero only at points of the immersed boundary) and f(x, t) is the Eulerian force.

Numerical Method
The fluid variables are defined in the Eulerian mesh N × N with x = (x i , y j ) = (ih, jh) for i, j = 0, 1, ..., N − 1, where h = ∆x = ∆y = L N is the length of each mesh range.It used the set o M Lagrangian points X = (X k , Y k ) with k = 0, 1, ..., M − 1 to discretize the immersed boundary, with initial spacing between points ∆s = L b M .The forcing term is defined at these points.We used an explicit scheme, where the force exerted on the border is calculated at the beginning of each step in time, with t n+1 = t n + ∆t and the numerical solution is given by: 1.The force field is calculated in points Lagrangian to the initial conditions.The force F n (s) is calculated using X n (s) in the immersed boundary and then the force F n (s) is used to determine f n (x) with the following equations: where the delta function is given by δ 2 h (x) = δ h (x)δ h (y) and discretized by 2. The update of the velocity field u n+1 (x) is performed by a Runge-Kutta time integration (4 a order): where V is a generic vector; the spatial variables (coupling pressure and velocity) are solved using the projection method and the convective derivatives are solved using the scheme of high order ("up-wind") VONOS, [7].
3. The velocity u n+1 is interpolated in the points of the immersed boundary and the new position is updated using X n+1 (s): where X n+1 (s) was calculated using the finite differences centered of the second order.

Modal analisys
For the airfoil section of mass m and moment of inertia about the center of mass, I cm , the kinetic energy, T , and the potential energy, P , of system are given by: The Lagrangian is given by The plunge, h, and pitch, α, movements are are obtained through Lagrange equation: which leads to the following equation of motion for the airfoil where M is the mass matrix, K is the stiffness matrix, U is the displacement vector and F is the vector of forces, then the equation of motion can be rewritten as: The stiffness matrix is diagonal because the foothold is at the center of flexion, where the plunge and pitch are statically decoupled.The lift and the aerodynamic moment in the pivot, per unit of length, are given by the expressions: where C L e C M are the lift and moment coefficients respectively.The equations of the motion (equations Eq.( 14) to Eq.( 16)) were resolved by classical method with direct approximation of eigenvalue to the case of undamped system, ie, this simplified approach can give a first estimate of the vibration limit where aerodynamic expressions are used to approach aerodynamic loads unstable, without incorporating damping conditions.The equation Eq.( 14), with forces given by equation Eq.( 17), reduces to: and the time-dependent displacement vector expressed by The nontrivial solution of the equation Eq.( 19) implies that Det K + A + λ 2 M = 0 which carries to three possibilities for the stability of the system: • for flow with subcritical velocity, with absence of damping, both eigenvalues (−λ 2 ) are real and positive and therefore the values for λ 2 are real and negative, ie the parameter λ = ±iω is purely imaginary.The two different corresponding values of λ purely imaginary return frequencies of the two circular modal branches in radians/second.The system vibrates constantly with harmonic motion simple in each of the two resulting branches from the vibration modes (vibration free) natural.
• after of critical velocity (u ∞ > V ), the parameters λ occur like complex conjugates λ = λ r ± iλ i in which one has positive real part λ r > 0 and the oscillations are unstable in the appropriate modal branch, characterized by increases in amplitude over time.
• the instability of the type of divergence is indicated by condition in which the imaginary part of λ disappears, λ i = ω = 0.
Accord [8], for small angles (x αm ≈ x cm ), then, the summarized form, the equations of motion for the typical section are given by: where f h is the natural frequency of translation uncoupled, f α is the natural frequency of torsion.The movements of translation h and the torsion α are taken as harmonics and there is a delay between them due to aerodynamic loads.For harmonic damped motion and using a linear approximation to the relationship between the aerodynamic forces and the variables of movement, the stability of the system is given by The values λ have the form where ω d is the frequency of damped oscillation and ζ is the damping ratio.So, for a given reduced frequency k, its possible to extract the eigenvalues, the damping ratio ζ and the velocity flow corresponding u ∞ .The reduced frequency is specified in aerodynamic simulation with forced oscillation for k = ω d u∞ , where c is the characteristic length (chord of section in the two-dimensional case).

RESULTS
For calculation of the eigenvalues λ were simulated cases of lifting and pitching, independently, with the parameters and characteristics shown in table Tab.1.The total force on the airfoil is a combination of viscous forces and of pressure.The pressure at each point of the airfoil is found from the C P and the viscous forces are given by the derivative of the tangential velocity with respect to normal, ∂u ∂n .Since the pressure and viscous forces are known in each point, they are numerically integrated to find the components of the total force and moment M of the airfoil.
We used a rectangular Cartesian grid with nonuniform grid 216 × 250 with ∆x min = 0.01, ∆Y min = 0.01, ∆x max = 0.35 and ∆Y max = 0.35 in a domain 6 × 3 dimensionless units, [9].The geometric center of the profile of NACA0012 airfoil was positioned at the point of coordinates 1.5 in direction x and 1.5 in direction y, Fig. 3 a) and b).The determination of the flutter depends on the structure elements of the system, such as mass and moment of inertia.Taking as reference the data presented for [10], not because it is a similar analysis, but for proximity between models treated, we will use the values shown in the table Tab.(1) and some parameters were chosen so that the frequency natural for plunge is 2Hz and 0.5Hz for pitching, such that f h and f α produce oscillations with the velocity reduced of the interest.Due to the limited availability of experimental results and simulation results in the literature for Re = 1000, was used, as reference, values taken from the work of [11] by to present frequency (7.16Hz) near to the studied in this work.
Besides the dimensioless parameters, took up the chord c as standard of scale, u ∞ as the characteristic velocity and ρu 2 ∞ 2 as pressure characteristic.The characteristic time is given by tu inf ty c , the characteristic velocity V * = u∞c π f and normalization of dimensionless equations provides the mass and moment of inertia ( M and Ī) and low frequencies ( fh and fα ).In practice, the values of f h and f α are chosen such that the airfoil undergo oscillations with velocity reduced of interest and M and I are tuned so that the loads will be higher (when reduced) or smaller (when increased) on the airfoil.The Table (2) shows the maximum values obtained for the lift and pitching moment and lift and plunging moment.For frequencies 0.5Hz, 0.75Hz, 1Hz and 2Hz was measured the reduced frequencies to be used in the calculation of eigenvalues.The results for the eigen- velocity is subcritical, ie, the eigenvalues have only imaginary parts.We also can see that there is a tendency to decrease the magnitude of the eigenvalues, that is, higher frequencies can lead to critical velocity of flutter and shown to be consistent as compared to those obtained by [11].

CONCLUSION
The present analysis with a simetric airfoil in a 2D flow (simulated by the Virtual Boundary Method) indicates that it is possible to simulate flutter parameters for ultralow Reynolds number, using classical eigenvalue approaches and direct integration methods in aerodynamic codes (solving Navier Stokes equations).The performance of an airfoil section NACA 0012 series subjected to forced oscillation (movement of lifting and pitching independently) in a flow with very low Reynolds number (Re = 1000) was studied using a twodimensional DNS simulation and showed that increasing the oscillation frequency increases the relationship lift-drag and has provided a standard for vibration prediction with frequency that dependents of parameters for preliminary determination of the flutter phenomenon.Also, modeling using modal analysis and immersed boundaries was effective for preliminary study of the phenomenon of flutter.

Figure 3 Table 1 .
Figure 3. a)Positioning of the profile in the mesh, b) Refinement of the mesh at the trailing edge.
(a), Fig.3(c), Fig.3(e) and Fig.3(g) show the values of C L for pitching and the development of sinusoid over time for the rotation α.Figures Fig. a

Figure 4 .
Figure 4. Results for determining the frequency reduced to pitching of ±2 o , plunging of 0.1c and frequencies of 0.5Hz, 0.75Hz, 1Hz and 2Hz.

Table 3 .
Results of simulations -eigenvalues