Direct Numerical Simulation of the Airfoil Segment’s Flutter and its Effect on the Aerodynamic Force

This article presents numerical simulation of planar potential flow around an airfoil with possibility of changing its shape. Two-dimensional unsteady flow model with scalar velocity potential, which allows us to calculate pressure distribution along an airfoil from Cauchy-Lagrange integral, is used. For this purpose, an airfoil contour is approximated by a complex cubic spline with possibility of displacement its vertices. This algorithm has been used in the context of fluid-structure interaction and has been applied successfully to determination of stability of an elastic airfoil segment interacting with a flow stream, so-called panel flutter problem. Calculation of external flow is carried out by vortex panel method with Kutta-Joukowski trailing edge condition, which makes mathematical solution unique. Using this method of approximation of an airfoil in combination with the method of discrete vortices provides a semi-analytical solution for complex potential for whole computational domain of air flow. This solution significantly accelerates process of numerical computation of time-averaged aerodynamic force as well as the dynamic stability problem for aeroelastic wing design and temporal evolution of its natural disturbances.


Introduction
Aeroelastic phenomena can occur to many engineering structures. In recent years, beneficial effects of aeroelasticity have received greater attention. For example, promise of new aerospace systems such as uninhabited air vehicles (UAVs) and morphing aircraft will undoubtedly be more fully realized by exploiting the benefits of aeroelasticity while mitigating risks.
Reliable aeroelastic analysis tools are needed to reduce expensive and time consuming experimental tests. In a later development an alternative approach has gained popularity which is known as "double-lattice"method [1,2] (DLM). The DLM is one of the most powerful tools for linear flutter analyses in subsonic regime. In this approach a lifting surface is divided into boxes (panels) and collocation is used for both downwash and pressure. For subsonic flows, numerical methods are in an advanced state of development. Paper by Rodden et al [3] contains an extensive discussion of the numerical techniques.
Roughen et al. [4] have compared results of several theoretical models with experimental data from a Benchmark Active Controls Technology (BACT) wing (see also [5]). The BACT wing is a rectangular planform with a NACA 0012 airfoil profile. The model has a trailing edge control surface extending from 45 to 75% span. Previously, Schuster et al. [6] had compared results from a Navier-Stokes CFD model (ENS3DAE) to these experimental data. Roughen et al. used an alternative Navier-Stokes CFD model (CFL3D) and a classical potential flow model (doublet-lattice). Correlations were made at several subsonic to transonic Mach numbers. As they note: . . . For the purely subsonic condition (M = 0.65) . . . . there is relatively good agreement between the doublet-lattice results, the Navier-Stokes results and the test data. This is not surprising because the flow is entirely subsonic and well behaved (there is no shock wave and no flow separation). . .
Another valuable correlation among several theoretical results and that of experiment is based upon the experimental work of Davis and Malcolm [7] (for more discussions, see Dowell [8]).
All above mentioned facts are reliable foundation for applying inviscid flow model to computation basic properties of fluid (air) flow.
Our fluid model is based upon Euler (momentum) equation of motion and does not cause significant difficulties for its numerical solution. Assuming flow is incompressible and irrotational and using the vortex panel method together with Cauchy-Lagrange integral for non-steady flow with proper boundary conditions allows us to relate dynamic pressure to scalar velocity potential.
The study of fluid-structure interaction, which known as panel flutter is a fundamental problem in modern aeroelastic-ity, first addressed by Fung [9] and Kobayashi [10]. This is a flutter resulting from fluid flow over an elastic plate. Dowell and coworkers [11][12][13][14] have done extensive numerical simulation of the panel flutter and mechanics of this problem. See also recent research results that have been obtained in the last decade in works of Vedeneev et al. [15][16][17][18][19].
Fluid-structure interaction solution requires that we solve fluid and structural models simultaneously for unknown pressure and body motion.
Airfoil profile shape is known to be an important parameter. This paper demonstrates a new airfoil approximation method with a sharp trailing edge by known control points (airfoil profile frame). This algorithm successfully applied in description of fluid-structure interaction of flow stream over an oscillating elastic segment (section) of an airfoil. Linear discretization of vorticity field and replacement of the airfoil by its approximate coordinates allows us to numerically calculate pressure distribution and scalar velocity potential both for steady and unsteady case. As a result of the linear discretization the problem is reduced to solution of a system of linear algebraic equations in unknown point vortices. To make the system complete an additional restriction is made, so-called Kutta-Joukowski trailing edge condition. This restriction states, in our case, that the point vortices on trailing edge must be the same magnitude and spinning in opposite directions to satisfy the Kutta condition.
Dividing entire flow field into two simply-connected domains enables to find a semi-analytical expressions for complex potential in each of these domains. These semianalytical expressions significantly accelerate computation of basic characteristics of fluid flow, as well as numerical calculation of time-averaged lift force of a wing.
The last chapter of the article is devoted to numerical study of stability of an airfoil which takes into account unsteady motion near oscillating segment. Dynamic system of fluidstructure interaction is described. Obtained results were studied from viewpoint of aerodynamic stability and bifurcation theory. The following regimes, which replace each other as free stream velocity increases, were obtained: stationary point in a proper phase space, limit cycle oscillations.

Airfoil approximation 2.1 Problem formulation
Consider widely used and computationally simplest interpolation spline of the third-degree. Also, we shall proceed from the assumption that knots of the spline p k are a = p 0 , p 1 , p 2 , . . . , p m−1 , p m = b, [a, b] ∈ R such that ξ k = ξ(p k ) ∈ C, k = 0, 1, . . . , m and besides ξ 0 = ξ m is tangency condition (Figure 1).
Let us define a complex cubic spline as a piecewisepolynomial complex-valued function of a real variable: So, the complex cubic spline function is defined as where obeying continuity conditions for the function itself and its first two derivatives: II (continuity conditions)   Substituting functions and their derivatives into the continuity conditions and assuming for simplicity An obvious next step is to evaluate unknown coefficients a k , b k , c k , d k from this system. We can eliminate the a k , b k , d k and reduce the problem in the unknowns c k .
We write expression for finding the coefficients c k , but reader is referred to the literature [20] for these more detailed formulations. Eliminating a k , b k , d k from this system gives A few general comments should be made about solving the system (2). For k = 2 coefficient c 0 is involved in this equation, which we put equal to zero (so-called, dummy coefficient). In order to apply standard tridiagonal matrix algorithm, to solve this tridiagonal system of equations, coefficient c m also required is equal to zero. So, we first consider c m = 0, but, after the system (2) is solved, we express c m from the contact condition: c m = 3d 1 h 1 −c 1 . We do not take into account the continuity condition for the 1 st derivative in the point , (derivative at point of contact (sharp edge) does not exist). This omission does not affect subsequent results. In this case, the complex cubic spline (1) just loses property to be a "natural".

Results of the approximation
Let us put the h k = 1 for simplicity. In this case, system (2) takes the form: Right sides of equations (3) are complex, while coefficient matrix is real. One may solve two systems for real and imaginary parts of the right sides of equations separately (due to linearity) and combine these results to obtain solution for the Advantage of this method is possibility of varying the control points to get a curve of arbitrary shape.
3 Vortex panel aerodynamic model

Mathematical formulation
This paper presents vibrations (flutter) of a particular segment of an airfoil, and this requires a certain method of approximation by the complex cubic spline. Airfoil contour is replaced by its approximated coordinates from the first chapter: where n = m × (number of points per interval) and the vortex panel method was used to construct complex potential.
Introducing complex potential as superposition of uniform flow at an angle α and potential due to effects of vorticity field around the airfoil we obtain where Here, i is imaginary unit and complex variable notation z = x + iy is used. Unless otherwise stated, α = 0 and . . . ds means integrating over the airfoil contour.
Complex conjugate velocity V is defined by In order to prevent confusion, let us also define complex velocity U = V . So, we have where u, v are fluid velocity components.
Replacement of continuous airfoil by its approximate coordinates z k from the first section gives polygon with n sides, called "panels". Let us number the panels clockwise with the 1 st panel starting on lower side, at trailing edge and the n th panel ending on trailing edge from upper side (Figure 4). Assume that each panel represents a planar vortex sheet with linearly varying strength and such that end strength of each panel is the same as starting strength of the next panel: The only exception is γ 1 = γ n+1 . Therefore, integral in right side of (5) may be computed analytically: Evaluating integral (6) and substituting above expression into (5) gives The only unknowns of this problem are strengths γ k , k = 1, . . . , n + 1. They can be found by solving n + 1 equations consisting of • n equations, expressing condition, that velocity at n mid-points of the panels (called "collocation points") is tangential to the panel; • one equation expressing the Kutta-Joukowski condition, which means, physically, that flow leaves sharp edge smoothly with no singularity. But, mathematically, it makes solution of Neumann boundary problem in outer region unique.
Collinearity condition of two complex numbers z 1 and z 2 states where Re(•) is real part.
Apply this condition at the n collocation points (midpanels) ( Figure 5) After various manipulations we obtain system of linear, algebraic equations where l = 1, . . . , n. Note, that for k = l in (8) we must take It means that angle between z 0 l − z l+1 and z 0 l − z l equals π ( Figure 5). Figure 5. Collocation geometry of two-dimensional aerodynamic model.
To find solution of system (8) with satisfactory accuracy and adequate level of smoothness, we must take about one hundred points of approximation.
Let us consider an explicit Runge-Kutta method with second-order accuracy. Advancement from time instant t 0 to the time instant t 0 + ∆t in this method is made according to the following expression: wherein as an initial condition one can take any point z 0 outside airfoil contour, which streamline should be started with. Figure 7 demonstrates streamlines relative to airfoil. In previous chapter we obtained formula for complex conjugate velocity (7). Let us integrate it with respect to z: where For simplicity, in what follows, an integration constant in (10) will be omitted without essential loss of generality. Principal difficulty of evaluating w(z) by applying (10) is presence of logarithmic term ln(z − z k+1 ).

Defining of single-valued branches
Complex logarithmic function is defined by Since infinitely different values of Arg(z) are obtained by successively encircling origin z = 0, complex logarithmic function is infinitely multi-valued.
In any simply-connected domain there exist infinitely many single-valued branches of multi-valued complex logarithm and its values of Arg(z) differ by integer multiples of 2π. And each such branch is a bijective function.
To keep the function single-valued, an artificial barrier that cannot be crossed must be introduced. This barrier is called branch cut, direction of which is arbitrary and largely a matter of convenience.
Let us divide entire flow domain, which is 2 − connected, into two simply-connected domains D + and D − , boundary of which consists of separatrix and the airfoil contour itself (Figure 8). Figure 8. Two simply-connected domains -D + and D − Therefore, in each of these domains we can define singlevalued branch of Ln(z − z k+1 ): with following equality: It was found, experimentally, that branch cuts of ln + and ln − are rays from origin to infinity at angles − π 2 and π 2 , respectively ( Figure 9).  Direct Numerical Simulation of the Airfoil Segment's Flutter and its Effect on the Aerodynamic Force So, among infinitely many branches of multi-valued function (10) two single-valued branches were determined: and

Comparison of methods
Once complex velocity U is known, we can find scalar velocity potential ϕ from condition, which states that U = ∇ϕ.
Let us identify the velocity potential ϕ as a function of z(t) , where z(t) = (x(t), y(t)) is coordinates of fluid element along an arbitrary curve in whole flow domain at the time instant t.
Taking into account that ∂ϕ ∂t ≡ 0 (steady flow), we obtain Let us apply the same explicit Runge-Kutta scheme to integrate the scalar velocity potential: However, it should be noted that above formula can be used if only coordinates of the curve, along which the integration is performed, are known a priori. Or, for example, these two recurrence relations (9) and (15) must be solved simultaneously for z m and ϕ m with given ∆t. Figure 10 demonstrates result of numerical integration scheme (15) with initial conditions ϕ 0 ± = Re[w ± (s 0 ± )], (Fig. 10B), and two branches of scalar velocity potential, which were calculated by using (12) and (13), (Fig. 10C), along upper S + and lower S − streamlines, respectively, where S ± ∈ D ± . Agreement between these methods was found to be reasonably good. Difference between these two approaches is negligible and is of the order of 1 × 10 −7 . Figure 10. Comparison of numerical Runge-Kutta scheme (B) with our analytical solution ϕ ± (z) = Re[w ± (z)] (C).

Circulation
We know that circulation Γ may be calculated as where z t.e. is trailing edge. Consider the following expression or, given the fact that w ≡ ϕ + iψ, where ϕ, ψ ∈ R, it becomes in addition, ψ + (z t.e. ) − ψ − (z t.e. ) ≡ 0, because z t.e. is coordinate of streamline (airfoil contour). On the other hand, considering (12) and (13), we have or, taking into account (11) where, obviously, n k=1 γ k +γ k+1 2 ∆s k ∈ R. Thus, we have obtained formula for the circulation Γ: We may obtain absolutely the same result by integrating vorticity field γ(s) over the airfoil contour: Physical interpretation of (16): γ k +γ k+1 2 is average value of complex velocity on the k th panel.

2D aeroelastic design
For further description of airflow-structure interaction an idealized physical model of elastic airfoil design was developed. Let us assume that all points z 1 , z 2 , . . . , z n+1 of the airfoil are interconnected by elastic weightless rod to provide elastic deformation. The main purpose of this paper is numerical study of average lift force change caused by elastic vibrations of a segment of the airfoil (Figure 11).
To discuss these phenomena we must first develop dynamic theoretical model. Let us fix static (stationary) airfoil with its coordinates In what follows, superscript * is related to this stationary airfoil, which remains motionless in any time.
Without loss of generality consider one of control points of the stationary airfoil ξ * δ , by which airfoil contour was approximated, and two adjacent points z * δ−i and z * δ+i , which lie on both sides of the ξ * δ (Figure 12). Consider a small deviationξ δ of the point ξ * δ from its initial position, after which the point moves to ξ * δ +ξ δ . At the same time we believe, that points z * δ−i and z * δ+i stay in their same locations. Describe motion of the disturbed point ξ δ by the following equations: with initial conditions which are natural from viewpoint of the experiment: where F δ is resultant force acting on the ξ δ ; coefficient m represents measure of inertia of an oscillating segment. It should be noted that this deviationξ δ must be the same order as the arc length of airfoil: Next step is to find all external applied forces and internal elastic forces acting on the point ξ δ . Neglecting gravity, it is clear that major contribution to the total motion is made by aerodynamic pressure F dyn (force per unit arc length) and tension force T res (Figure 12).

Resultant vector T res
We put parameter i = 1. Now we must express the resultant force T res in terms of the z * δ−1 , ξ * δ , z * δ+1 and ξ * δ +ξ δ . All explanations will be carried out for T 1 (Figure 12), for T 2 analogous reasoning leads to similar result.
Hooke' s Law reads tion) of the arc length between ξ * δ and z * δ+1 , E is modulus of elasticity.
We obtain

Dynamic pressure. The Cauchy-Lagrange integral
This perturbationξ δ triggers an aerodynamic flow so that now we have to take into account unsteady motion near the oscillating segment.
Cauchy-Lagrange integral may serve for the same purpose as Bernoulli's equation. If potential of velocity ϕ is known then pressure may be found with the aids of this integral. The practical value of Cauchy-Lagrange integral is that it allows one to relate p to ϕ. Using it we obtain We may evaluate f (t) by examining fluid at some point where we know its state. For example, far away from body. Then f (t) becomes U 2 ∞ 2 + p∞ ρ ≡ f ∞ and we find that f ∞ is a constant independent of space and time. Hence finally Let us express pressure p: However, to describe motion of the ξ δ we need total derivative of ϕ with respect to t: Expressing from (24) partial derivative: and substituting it into (23) gives where c 3 is length of the oscillating segment, which due to very small oscillations (amplitude of the vibrations is much smaller than characteristic length of airfoil) is assumed constant.

Simulation Results
In this section, results of simulation will be analysed using theory of aerodynamic stability and nonlinear dynamical systems. Unless otherwise stated, m = 1, ρ = 1 and as investigated control point ξ * δ we take one of control points of lower surface of airfoil ( Figure 13).

2D solution and its stability
Simulations showed ( Figure 14) that for relatively small values of free stream velocity U ∞ stable solutions of the system (17)- (22) are realized.
For the control parameter U ∞ = U * ∞ and fixed E, solution of (17)-(22) suffers a bifurcation (Fig. 15B). This system has a phase space, and it was crucial for its understanding to choose proper projection of this space. So, in order to clarify the behaviour of the system, it was instructive to consider trajectories of perturbed point ξ δ .
Numerical simulations of (17)- (18) with initial conditions (19)- (22) confirmed birth of a stable cycle. This is illustrated in Fig. 15B (left), where it is shown that the solution acquired periodic oscillations in time. It is interesting to note, that trajectory of perturbed point, in this case, takes form of a straight segment (Fig. 15B, right).
This attracting-repelling scenario ran over many times. However, it is not our purpose to probe deeply into this problem here; for a thorough treatment reader is referred to Dowell [11] and Bolotin [21]. The reader may wish to consider another specific issues of the theory of aeroelastic instability (see, for example, Vedeneev, Shishaeva et al. [22][23][24][25][26][27][28][29][30][31]) Fig. 16A demonstrates a family of plots for various values E where oscillations of scalar potential are presented. One can see that mean values of oscillations, as modulus of elasticity increases, tend to their limiting value of potential of steady-state airfoil in point ξ * δ (Fig. 16B).

Lift
It is interesting to examine total force (lift) L on this aeroelastic effect. Pressure distribution (23) may be integrated along airfoil contour to obtain the total force (lift) on a wing: where n is normal unit vector pointing into the wing and e j is vertical unit vector, normal to freestream direction.
In our case (α = 0), we have e j ≡ i and (25) becomes so that now n y is y-component of outward-pointing unit normal vector.

Conclusions
Results of numerical simulation of aerodynamic instability, within context of airflow-structure interaction, were presented and understood from viewpoint of bifurcation and stability theory. In context of nonlinear system theory, a limit cycle oscillation LCO is one of the simplest dynamic bifurcations, "a first stop on the road to chaos", so to speak.
Despite the fact that we considered a highly simplified model (somewhat analogous to elastic design of an airfoil), it was possible to bring out some of essential features of this type of problem, thus the model serves a useful purpose in introducing this type of flutter.
In this paper, oscillations of airfoil segment have a natural (uncontrolled) character, so that they provides an insignificant augmentation of lift force. Up to a certain magnitude of deviation of the perturbed point oscillations result in a transfer of energy from flow to body. This transfer may be viewed as equivalent to a flow-induced negative aerodynamic damping. For larger deviations, however, there occurs a transfer of energy from the body to the airstream. Control of this process with the aim of more effective increasing of total lift force is the main purpose of next studies.
Optimal mechanism of the interaction for energy transfer from airstream to oscillating airfoil should include at least two oscillating segments on lower and upper surfaces of a wing. This mechanism should be able to control both dissi-pation and forced oscillations of segments themselves.
This control implies a specific connection between the friction forces and flow-induced forces, acting on oscillating sections of airfoil, to pump energy from airstream into vibration and thus sustain the presumed motion.