On mathematical modelling of aeroelastic problems with finite element method

This paper is interested in solution of two-dimensional aeroelastic problems. Two mathematical models are compared for a benchmark problem. First, the classical approach of linearized aerodynamical forces is described to determine the aeroelastic instability and the aeroelastic response in terms of frequency and damping coefficient. This approach is compared to the coupled fluid-structure model solved with the aid of finite element method used for approximation of the incompressible Navier-Stokes equations. The finite element approximations are coupled to the non-linear motion equations of a flexibly supported airfoil. Both methods are first compared for the case of small displacement, where the linearized approach can be well adopted. The influence of nonlinearities for the case of post-critical regime is discussed.


Introduction
There are many technical or scientific applications in a wide range of technical disciplines, where the mathematical modelling of fluid-structure interaction (FSI) problems is important. Let us mention for instance a wide range of aeroelastic tools used in aerospace engineering to investigate the aircraft safety, see e.g. [1], aeroelastic stability of bridges in civil engineering, see e.g. [2], or scientific applications in biomechanics, e.g. [3] or [4]. The classical aeroelastic approach based on the asymptotic aeroelasticity is fast, efficient and reliable for thin airfoils, see e.g. [5], and that is why it is still very popular in the technical practice.
The aeroelasticians in practice use the linearized methods to find the solution of aeroelastic problems usually in terms of the prediction of the aeroelastic instability, see [6]. The use of linearized aerodynamics allows to determine the aeroelastic asymptotic stability in terms of the critical velocity, damping and frequency graphs, see [1], [7], or [8]. Although similar approaches are well studied it needs to be mentioned that the asymptotic stability determines the critical velocity, but does not guarantee the safety, see [9,10]. The external excitation can lead to transient growth and consequently to the structural failure, although the system is aeroelastically stable. Thus the use of an alternative accurate aeroelastic simulations is very attractive because it reduces the development risk, the number of experiments and the possibile future design modifications. which could provide additional information to that one provided from the linearized approach. ⋆ e-mail: petr.svacek@fs.cvut.cz However, the mathematical modelling of FSI problems in general is much more complicated, because one needs to consider the viscous possibly turbulent flow which is moreover in the interaction with the nonlinear behavior of the elastic structure. Moreover, the time changes of the flow domain due to the vibrating structure cannot be neglected and the methods applicable on moving domain or meshes must be employed, see e.g. [11]. Last, the coupled fluid-structure system for the fluid flow and for the oscillating structure needs to be solved simultaneously using a coupled strategy at any time instant, which possibly can create an artificial instability of such a scheme, see e.g. [12].
In applications the FSI problems are usually solved with the aid of the so-called partitioned or de-coupled approach, which approximates the fluid flow and structure motion by different solvers. The coupling is then enforced by using suitable conditions prescribed at the fluidstructure interface. These approaches can be characterized further either as strongly coupled scheme (if subiterations are used to enforce the interface conditions) or as loosely coupled one. Such a loosely coupled scheme was used e.g. in [13] for FSI problems tested on a proposed benchmark problem and the unconditional stability was derived based on an energy estimate. Nevertheless the loosely coupled algorithms were observed in many cases to be instable due to the poor fluid-solid coupling, see [12]. In order to get a stable scheme either the partitioned approach with relaxed strongly coupled algorithm or the monolithic approaches can be used, cf. [14,15] or [16]. Nevertheless, the high fidelity aeroelastic models are still very demanding to be solved and thus its direct use in industry is rather rare, see [17]. For the solution of each part in FSI problem a suitable numerical method needs to be applied. One of the most popular methods is the finite element method, but in its application for numerical simulation of incompressible flow problem one must overcome several sources of instabilities. Particularly, the possible instability caused by the incompatibility of the pressure and velocity pairs of finite elements needs to be taken into account, see e.g. [18], [19]. Furthermore, the other instability is caused by the dominating convection terms, see [20]. Let us mention that there is also another not so well-studied instability source in the Galerkin discretization method which is related to the possible poor resolution of pressure, see [21], [22], [23]. For the case of high Reynolds numbers the turbulence effects should be also taken into account using a wide variate of the mathematical models, see e.g. [24], [25], [26], [27] or [28] as few examples how to include a turbulence model in the numerical simulation.
The other thing is that for approximation of FSI problem the applied numerical method should be able to treat the moving domain/meshes. The most popular method is the arbitrary Lagrangian-Eulerian(ALE) method, see [29]. There are several possibilities how to apply the finite element method, see [30]. Usually in order to provide higher order accuracy or stability some additional assumptions on the ALE are required, see [31], [32]. Particularly the socalled geometrical conservation law is important, see [33].
In this paper the numerical analysis of an interaction of the incompressible viscous flow with a vibrating airfoil is considered, see e.g. [30], [34] or [35]. The attention is paid to an approximation on moving meshes with the aid of suitable arbitrary Lagrangian-Eulerian(ALE) method. For the case of high Reynolds numbers anisotropically refined mesh is used and the choice of the stabilizing parameters is discussed. Two finite element pairs are considered for the numerical approximation and the stabilization procedure for both of them is presented based on the so-called residual based stabilization, see e.g. [23]. A benchmark problem of mutual interaction of the airflow and a flexibly supported airfoil with two degrees of freedom (2DOF) (the bending and the torsion modes) is considered. The problem is solved using a partitioned scheme with strongly coupled/loosely-coupled scheme. The problem of the solution of the non-linear flow problem per-each time step is discussed. The paper is organized as follows: first the mathematical model is presented, then the weak form suitable for application of the finite element method is introduced, then the numerical approximation is presented and the numerical results are shown.

Structure vibrations
In this paper an airfoil section that can be vertically displaced by h (downwards positive) and rotated by angle α (clockwise positive) is considered, see Fig. 1. These two degrees of freedom h and α are also called the bending and the pitching, respectively. The nonlinear equations of motion then read (see [30]) where m is the mass of the airfoil, S α is the static moment around the elastic axis (EA), and I α is the inertia moment around EA. The parameters k h and k α denote the stiffness coefficients. On the right-hand side the aerodynamical lift force L(t) and aerodynamical torsional moment Under the assumption of small displacements (i.e. α ≈ 0, sin α ≈ 0 and cos α ≈ 1) the system of equations can be linearized as In order to solve system of ordinary differential equations (1) or (2) the aerodynamical lift force L(t) and the aerodynamical moment M(t) needs to be evaluated for any time instant t. These aerodynamical quantities can be expressed with the aid of flow quantitiesas pressure p, flow velocity u, fluid density ρ and fluid kinematic viscosity ν, or using a suitable aerodynamical model. The system of equations (2) is also completed by the prescription of suitable initial conditions at t 0 = 0 for α(t 0 ), h(t 0 ),α(t 0 ) andḣ(t 0 ).

Aerodynamical model
Two aerodynamical models are considered in this paper. First, the aerodynamical forces are determined with for the case of a thin oscillating airfoil with focus on the determination of the aeroelastic stability of a vibrating airfoil. Using the asymptotic stability approach only the critical velocity needs to be determined, which is characterized as the flow velocity for which the airfoil starts to harmonically oscillate. This means that the applied aerodynamical model does not have to describe the aerodynamical forces in general but only for the case of harmonically vibrating airfoil.
Second, the aerodynamical forces are evaluated using the surface integral of the density of the aerodynamical forces obtained from the flow velocity and pressure. The flow velocity and pressure are then assumed to be solution of the incompressible Navier-Stokes equations. The solution of the Navier-Stokes equations is approximated using a stabilized finite element method.

Linearized model.
The first model of the linearized aerodynamics is based on the assumption of thin airfoil, for which the aerodynamical forces acting on the harmonically oscillating airfoil can be evaluated. Using the potential function approach, conformal transformation and writting the solution in terms of known analytical solutions, the aerodynamical lift force and moment can be expressed, [7]. Moreover, these aerodynamical forces depend only on the reduced frequency k = bω/U ∞ and the far field velocity U ∞ (see [8]) as where x f denotes the distance of the elastic axis from the leading edge, b = c 2 denotes the midchord of the airfoil (see Fig. 2, d e = x f − c 4 is the distance of the elastic axis from the aerodynamical center, d m = x f − b is the distance of the elastic axis from the midchord position, d t = 3 4 c− x f is the distance of the elastic axis from the 75 % of the airfoil and C(k) is the Theodorsen function.

Incompressible Navier-Stokes equations
The flow velocity u and the pressure p needs to be determined for instance using a flow model. Here, the flow in the bounded time-dependent domain Ω t ⊂ R 2 is governed by the system of incompressible Navier-Stokes equations coupled with the continuity equation where u = u(x, t) is the flow velocity, p = p(x, t) is the pressure, ρ is the constant fluid density ρ ∞ , and µ is the fluid viscosity. In order to treat the flow in time dependent computational domain the arbitrary Lagrangian Eulerian method is used, for the sake of brevity the description of the ALE method is omitted here. The boundary of the computational domain ∂Ω t consists of mutually disjoint parts Γ I (inlet), Γ O (outlet), Γ Wt (the moving part, vibrating airfoil). The following boundary conditions are prescribed where the domain velocity w D is on Γ Wt equal to the local boundary velocity, u is the flow velocity, p is the kinematic pressure, n denotes the unit outward normal vector to Ω t and α − = min(α, 0) denotes the negative part of a real number α. The computations start from the initial condition u(x, 0) = u 0 (x) for x ∈ Ω 0 . The aeroelastic lift force and aerodynamical moment can be also determined using the flow velocity and pressure, i.e. the aerodynamical force read and the aerodynamical moment is given by where

Numerical simulations
Here, two different approaches were used for numerical simulation using the two aerodynamical models mentioned above. First, the model of linearized aerodynamics is used. In this case however two types of instability needs to be considered, the static instability called the divergence type of instability the dynamic type if instability known as flutter. For the finite element method application we focus here only on approximation of the fluid flow.

Static instability.
First, we consider the problem of static instability (called divergence type of instability). To this end we introduce the non-dimensional lift coefficient c L defined as where U ∞ denotes the far-field flow velocity, c denotes the airfoil chord, ρ denotes the constant fluid density and l denotes the depth of the considered airfoil section (the dimension of the airfoil section in the z-direction). Furthermore, for the thin airfoil the lift coefficient (approximately) depends on the angle α linearly as where α 0 denotes the angle of attack of zero lift force, which for the sake of simplicity is asssumed to equal zero in what follows, i.e. α 0 = 0. Similarly, the nondimensional aerodynamical moment around the aerodynamical center x AC (aerodynamical center is the point of the airfoil for which the aerodynamical moment is independent of the angle of attack α, for subsonic flow it is 1/4 of the airfoil) reads The aerodynamical lift force then reads L = πρU 2 ∞ clα and the aerodynamical moment with respect to the elastic axis can be expressed as Assuming the static case (i.e.α =ḣ = 0 similarly as the second derivatives) and using the expressions for L and M in (2) we get the homogenous system of linear equations whose solution depends on the value If this value is non-zero, there exists a unique solution. The interesting case is the zero value, which means that there is no (unique) solution of the system and consequently no static equilibrium exists. This means that the system is aeroelastically unstable with the divergence type of instability and the velocity U div for which it occurs is called the divergence type of instability. This velocity can be expressed as

Flutter determinant and pk-method
Using the linearization (3) in the linearized system of equations (1) leads to the system of two nonlinear equations with two parameters, the reduced frequency k and the far field velocity U ∞ . Based on this nonlinear system the critical flutter velocity can be determined and using the so-called p − k method the vibration frequency and damping for a given flow velocity U ∞ can be computed.

Finite element method for flow problem
In order to discretize the problem (4) in time, the equidistant partition 0 = t 0 < t 1 < · · · < T, t k = k∆t of the time interval [0, T ] is considered with the constant time step ∆t > 0. The time derivative in (4) is approximated by the second order backward difference formula The time discretized problem is than spatially discretized by the finite element method. To this end the approximations u k and p k are sought in the finite element spaces W h and Q h , respectively. In the practical computations these spaces are defined with the aid an admissible triangulation T h of Ω and the well known Taylor-Hood finite element is used. The discrete formulation contains the Galerkin terms defined for any U = (u, p) and V = (ϕ, q) by The stabilization terms containing the SUPG/PSPG stabilization are defined by and the stabilization term for div-div stabilization reads P(U, V) = K∈T h τ K (∇ · u, ∇ · ϕ) K with the choice of δ K and τ K according to [30]. This approach was previously numerically verified on a number of benchmark problems, see [30], and also it satisfy the optimal error estimates and thus it should be more precise compared e.g. to the classical finite volume methods. Then the non-linear stabilized discrete problem at a time instant t = t n+1 reads: Find p n+1 := p, u n+1 := u, U = (u, p) ∈ W h × Q h , such that u satisfies approximately the conditions (6 a-b) and holds for all V = (ϕ, q) ∈ X h ×Q h . The solution of the nonlinear problem is realized by the Oseen linearization. The nonlinear problem (12) is then coupled with the system (1) and solved simultaneously.

Numerical results
In this paper two aeroelastic problems were considered, see [30] and [7]. The first one is an example of the divergence type of instability and with the further increase of the far field velocity the flutter type of instability occurs. The numerical approximation of such an example can be a little bit difficult particularly for the post-critical velocities, where the behaviour of the aeroelastic is nonlinear and for instance its behaviour can be dependent on the initial displacement. The second example is the classical flutter type of instability, where the instability occurs due to the coupling between the two modes of bending and torsion.

Divergence type of instability case
A model problem originally studied in [30] is considered here. In this case the airfoil NACA 0012 was considered in interaction with incompressible fluid with the constant densityρ = 1.225kg/m 3 , the fluid kinematic viscosity ν = µ/ρ = 1.5 × 10 −5 m 2 s, the airfoil chord c = 0.3 m, the airfoil airfoil mass m = 1.7324 kg (related to the depth 1 m), the airfoil inertia moment I a = 0.0097458 kg m 2 , the airfoil static moment S a = −0.015592 kg m, the bending stiffness was k h = 2102.2 N / m and the torsional stiffness k a = 73.912 N m / rad. The elastic axis was located at 40 percent of the airfoil chord measured from the leading edge of the airfoil, x f = 0.4c = 0.12 m and the center of gravity of the airfoil was located at 37 percent of the airfoil, i.e. x cg = 0.37c. The flutter velocity was determined approximately as 42.5 m/s, whereas the system becomes unstable due to divergence type of instability at lower far field velocity U ∞ = 37.7 m/s. The linearized approach prediction (realized with the aid of MATLAB/Octave script) in terms of aerodynamical damping and vibrations frequency is shown in Fig. 3 in terms of frequency and in Fig. 4 in terms of damping. In the graphs the flutter velocity 42.5 can be identified from the graph as the lowest positive velocity with the zero damping coefficient. Further, the finite element method was applied on the same aeroelastic problem and the similar results in terms of frequency and damping was obtained for the sub-critical far field velocities in the range 5 -35 m/s, see Fig. 4. However, with the increase of the far field velocity a difference can be observed namely in the damping coefficient. This is particularly due to the fact, that the applied numerical approach is able to capture the divergence type of instability, which is not captured by the linearized approach. This also corresponds well with the results for the post-critical velocity U ∞ = 40 m/s and different initial displacement α(0) = 0.03, 0.3, 3 and 6 degrees. The aeroelastic response shown in Figs. 5, 6 have the typical divergence type of instability behaviour, which especially for the initial displacement of 6 degrees start to resemble rather the flutter type of instability. The flutter velocity determined by the linearized approach is however determined approximately as 42.5 m/s, see Fig. 3-4. This decrease of the flutter velocity due to the high initial displacement can not be captured by the linearized approach, but it well agrees with the real physical phenomena.

Flutter instability case
For the second considered case the airfoil NACA 0012 was considered in interaction with incompressible fluid with the constant densityρ = 1.225kg/m 3 , the fluid kinematic viscosity ν = µ/ρ = 1.

Conclusion
This paper described two possible strategies of solution of an classical aeroelastic problem. The first on based on the classical approach of linearized aerodynamics was susccessfully applied for solution of two aeroelastic problems from references. Its results were also compared to the other method described in this paper, which is the finite element approximations of the FSI problem. Both approaches were found to be very similar in terms of the predicted frequency and damping for sub-critical velocities. However for far field velocities close to the critical one it was found that the latter approach based on the finite element method can be influenced by the nonlinearities presented in the system as large displacements, high initial displacement etc.