Aeroelastic analysis of membrane airfoils and flexible-chord airfoils with trailing-edge flaps

This paper studies the static and dynamic aeroelastic characteristics of membrane airfoils and flexible-chord airfoils (deformable airfoils) with the emphasis on the effects of a trailing-edge (TE) flap which is a novel topic. Two modeling approaches are presented; the first method is the Rayleigh–Ritz method, and the second method is the finite element method which is an efficient method to study the TE flap effects. The two models are presented in the Laplace domain which enables the transient response analysis. The models adopt the potential flow aerodynamics based on the Prandtl–Glauert thin-airfoil theory and the Theodorsen’s unsteady theory. The airfoils are assumed to have small deformations, so linear structural models are used. The effect of the airfoils’ flexibilities on the static aeroelastic characteristics and the dynamic responses due to step and harmonic TE flap inputs is presented through a parametric study.


Introduction
The interest in the design of micro-air vehicles (MAVs) and the increasing flexibility of modern aircraft has led aerospace engineers and scientists to study the aerodynamic characteristics of flexible airfoils overs several decades. There are two types of flexible airfoils; the first one is the membrane airfoil which is a type of twodimensional structures with a very thin thickness such that its bending stiffness is negligible, so it supports the applied external transverse loads by the in-plane tension forces throughout the curvature of its deformed shape. In the field of aerospace engineering, parachutes, hang gliders, and some micro-air vehicles are examples of membrane structures and in nature, wings of insects and bats are perfect examples of membrane wings. The second type is the flexible-chord airfoil which is an airfoil with low bending stiffness such that it suffers from camber deformations due to the air flow. Such airfoils are used also for the design of MAVs with very thin plate-like wings or for larger morphing air vehicles where the airfoil internal structure flexibility is designed to be low to allow passive or active morphing capabilities [1].
The aeroelastic analysis of membranes has been an ongoing research since the beginning of the nineteenth century by studying parachutes and yacht sails performances [2][3][4][5]. Nowadays, the development of microair vehicles and the desire to imitate the performance of flying insects, birds, or animals has attracted the interests of many researchers in the past few years [6]. One of the earliest work in the field membrane airfoil aeroelasticity is the work of Thwaites [7] where an analytical solution was presented to get a relation between the static membrane deformation and the applied aerodynamic loading. The linear sail differential equation was used along with the potential flow theory. The work of Nielsen [8] also theoretically studied the linear static aeroelastic behavior of membrane airfoils and compared the results with experiments. The Fourier series expansion was used to solve the sail governing differential equation and the thin-airfoil theory was adopted to model the aerodynamic loading. The theoretical static deformation shape of the airfoil and the lift coefficient were in a good agreement with the experiments for the range of Reynold's number 5.6-7.8 (10) 5 and for maximum camber less than 15%. Jackson [9] assumed the deformed membrane to be a cubic function with two unknown coefficients. These coefficients were estimated by a method of aerodynamic forces balancing and geometric constraints. Newman and Low [10] used potential thin-airfoil theory as Nielson and Thwaites did, but they applied the Kutta condition at the leading edge instead of the trailing edge in an effort to simulate the effect of the boundary layer bubble at the trailing edge. The theoretical results were compared with experiments for four excess-length ratios, but not good agreement was found. Unlike the previously mentioned models where the potential thin-airfoil theory with continuous vortex sheet along the airfoil chord was used for the aerodynamic modeling, Greenhalgh et al. [11] used a discrete vortex sheet by dividing the membrane chord into number of segments with a concentrated vortex located at the quarter point of each segment. A good match with the experimental results was reported for small angle of attacks.
Newman and Paidousses [12] studied the dynamic instability (Luffing) of membrane sail under very low tension forces. An approximate analytical closed form solution using the potential flow theory was presented to specify the stability limits. Mavroyiakoumou and Alben [13,14] studied the stability of membrane airfoils with different types of boundary conditions: fixed-fixed, fixed-free, and free-free edges. Inviscid potential flow model including the wake effect of Saffman [15] was used. Starting by the membrane differential equation, a nonlinear eigen-value problem was derived that was solved iteratively to estimate the stability limits. Sygulski [16] presented an integral aeroelastic formulation that was solved by the boundary element method and a differential equation formulation that was solved by the finite difference method. Potential flow theory with wake modeling was used. Tiomkin and Raveh [17] presented a Laplace domain formulation based on the Fourier series solution adopted by Nielsen [8] and the Schwarz's unsteady aerodynamic model. The static and dynamic stability limits of membrane airfoils were studied for different values of the membrane mass density. Tsesana and Breuer [18] studied the steady and unsteady aeroelastic behaviors of flapping membrane airfoils at different frequencies. The membrane deformation was represented by a series solution, and the potential flow model of Wu [19] was chosen to estimate the unsteady aerodynamic loadings. The resulting aeroelastic ordinary differential equation was solved by the Chebyshev collocation method iteratively. Recently, Hussein [20] proposed a finite element formulation to study the aeroelasticity of membrane airfoils with elastic boundary conditions using the potential flow aerodynamics which is extended in this work to model membrane and deformable airfoils with TE flaps.
Another approach for the aerodynamic modeling of membrane airfoils is based on computational fluid dynamics models (CFD) to capture some effects that cannot be predicted by the inviscid potential flow method like the effect of the leading-and trailing-edge shape [21], the viscous effects [22][23][24], and large angle of attacks [25].
Moving to the analysis of deformable airfoils, Su [26] developed an aeroelastic model of flexible-chord airfoils using Legendre polynomial series representation of the airfoil camber deformation and the finite-states potential flow theory of Johnson and Peters [27] to capture the unsteady aerodynamic loadings on the airfoil. The Hamilton's principle was used to derive the system of differential equations describing the aeroelastic problem. Walker and Patil [28] studied the unsteady aerodynamics of deformable thin airfoils where the camber deformation is presumed and represented by Chebyshev polynomials series. An aerodynamic theory is derived to estimate the lift, moment, and thrust generated by a harmonically deformed airfoil. Murua et al. [29] analyzed the effect of the camber deformation on the dynamic instability (flutter) of airfoils. The airfoil camber deformation is predefined by a single symmetric quadratic function besides the two rigid-body degrees of freedom and the airfoil was supported by two translational springs at equal distances from the mid-chord point. Berci et al. [30] proposed a semi-analytical method to study the aeroelastic instability and response of flexible airfoils. The finite-states inflow theory of Peters was used to estimate the unsteady aerodynamic loading. The camber deformation is represented by Chebyshev polynomials, and the aeroelastic system of equation is derived using the Ritz method. Riso et al. [31] also presented a semi-analytical formulation for the unsteady aerodynamics of flexible flat plate airfoils. A potential flow conformal mapping technique is derived to estimate the aerodynamic coefficients for a cantilevered flexible airfoil with a predefined deformation. CFD approaches are also available in the literature to study the aerodynamics of deformable airfoils [32][33][34].
The present work studies the static and dynamic aeroelastic responses of membrane and flexible-chord airfoils with TE flaps which have not been covered before in literature. The aeroelastic formulation is derived and presented using the Rayleigh-Ritz approach and the finite element approach. The Theodorsen's unsteady aerodynamic theory and the Prandtl-Glauert thin-airfoil theory are used to present the aerodynamic loading in an efficient matrix form in terms of the airfoil degrees of freedom. The static, transient, and frequency responses of the flexible airfoils due to step and harmonic TE flap inputs are presented and discussed via a parametric analysis.

Flapless airfoils
The aeroelastic modeling of flexible airfoils without a trailing-edge flap can be done using the Rayleigh-Ritz approach which approximates the global deformation as a series solution with coefficients to be determined using the Euler-Lagrange equation which can be in terms of the elastic strain energy (U ), the Kinetic energy (K), and the external work done (W ) as follows: where η n are the coefficients of the series solution expansion which can be written for a membrane airfoil with (n s ) terms shown in Fig. 1 as follows: For a membrane airfoil with a chord (c) and mas per unit length ( ) that is subjected to a tension force (T ) and an aerodynamic pressure p(x, t), the elastic strain energy, the kinetic energy, and the external work can be written as follows: Substituting Eqs. (2) and (3) into Eq. (1) yields the governing equations of motion which can be written in the matrix form as follows: The mass matrix [M] and the stiffness matrix [K ] are given in Appendix. The pressure term in Eq. (4) has three components: the quasi-steady p QS (x, t), the added mass p M (x, t), and the wake component p W (x, t).
The quasi-steady aerodynamic pressure distribution is modeled using the Prandtl-Glauert thin-airfoil potential flow theory which estimates the pressure due to a sheet of bound vortices γ b (x, t) and a free stream dynamic pressure q ∞ 1 2 ρ ∞ V 2 ∞ using a finite series solution with n a +1 terms and a transformed coordinate θ a cos 1 − 2x c as follows [35]: The aerodynamic series coefficients can be estimated due to the membrane deformation, the airfoil angle of attack (α), and the free stream speed (V ∞ ) as follows: Substituting the series representation of the membrane deformation into Eq. (6) gives the following matrix form of the aerodynamic coefficients (see Appendix): Substituting Eq. (7) into Eq. (5) yields the following expression for the quasi-steady pressure contribution in Eq. (4): The Theodorsen's unsteady thin-airfoil theory is used to estimate the added mass aerodynamic pressure coefficient and the wake effect pressure due to a harmonic motion of the membrane surface at a certain reduced frequency k ωb/V ∞ which causes a downwash speed amplitude (x) as follows [36,37]: The Theodorsen's function C(k) and the inertia function (x * , ξ * ) are given in Appendix. The down wash speed amplitude ( ) is related to the membrane surface deformation amplitude (w) as follows: Substituting Eq. (9) into Eq. (10) will give the following expressions for the contribution of the added-mass pressure and the wake effect pressure in Eq. (4) as follows: Undeformed Airfoil For a flexible-chord airfoil shown in Fig. 2, the deformation can be represented using the free-free vibration mode shapes of an elastic beam as follows [38]:

Deformed Airfoil
The functions ϕ n (x) are as follows: where x s is the location of the point to which the translational and rotational springs are attached which is analogous to the elastic axis of a wing. The parameter λ n is obtained by solving the following solution: The elastic strain energy of the airfoil is written in terms of the bending stiffness (E I ), the translational and rotational springs stiffnesses (k w , k φ ) as follows: The expressions of the kinetic energy and the work done are the same as in Eq. (3), so following the exact procedure as in the case of a membrane airfoil, the aeroelastic equations of motion can be easily obtained.

Flapped airfoils
In order to add the effect of the TE flap shown in Fig. 3, an additional TE deflection function ϕ TE (x) is introduced which can be written as follows: For the flexible-chord airfoil, the assumed deflection function will be the same as defined in Eq. (14), while for the membrane airfoil, the functions have to be redefined to be as follows: Therefore, the final governing system of equations of a flapped airfoil can be written in the Laplace-domain matrix form through replacing iω by the Laplace variable (S) and using the Laplace form of the Theodorsen's function F(S) which yields the following form:

Finite element formulation
This section presents the finite element aeroelastic formulation for a membrane airfoil and a flexible-chord airfoil. The airfoil is discretized into (n s ) elements with length (l e ) as shown in Fig. 4. Starting with the modeling of a membrane airfoil, the vertical displacement is approximated over each element in terms of its nodal displacement vector − → q w e w 1e w 2e T through a shape function [N s ] as follows [20]: The principle of virtual work is used in this work to derive the elemental matrices. So, the virtual strain energy in a membrane element can be written in terms of the element stiffness matrix [K e ] and the element displacement vector − → q w e as follows: Similarly, the virtual work due the inertial loading (δW e m ) can be written in terms of the element mass [M e ] matrix as follows: And the external virtual work due to the aerodynamic pressure loading (δW p ) will be: As has been done in the Rayleigh-Ritz approach, the Prandtl-Glauert theory is used for the quasi-steady model. Therefore, the Prandtl's expansion coefficients can be written in terms of the membrane displacement vector − → q w as follows: The matrices [A s ] e and[A d ] e are of size 2×(n s + 1). The finite element formulation of the added-mass pressure coefficient can be written as follows: The membrane nondimensional coordinates (x * , ξ * ) are related to the elemental coordinates (s, η) as follows: So, the elemental external virtual work of the added-mass pressure δW M p will be: Similarly, the external virtual work of the wake term pressure δW W p will be: The aerodynamic modeling of a flexible-chord airfoil (without a flap) shown in Fig. 2 is exactly the same as a membrane airfoil except the definition of the shape function which will be than of a flexible beam. So, the nodal displacement vector will contain the nodal rotations besides the nodal translations as { q w } e w 1e θ 1e w 2e θ 2e T and the shape function of the vertical displacement will be as follows: The virtual elastic strain energy of an element can be written of its bending stiffness (E I ) follows: The effect of the translation spring (k w ) and the rotational spring (k φ ) is simply added to the entries of the stiffness matrix corresponding to the node where the springs are attached. Therefore, the final finite element aeroelastic formulation for a flapless airfoil will be as follows: The advantage of the finite element formulation over the Rayleigh-Ritz formulation is that it can be easily extended to model more complex geometries without any reformulations. So, an airfoil with a trailing-edge flap shown in Fig. 3 can be easily modeled by modifying Eq. (33) to be as follows: The aerodynamic force vectors due to the flap deflection f M β , f D β , and f S β can be easily computed from the aerodynamic matrices using the principle of superposition such that the total airfoil deformation can be decomposed into two components; the first component is the undeformed airfoil with a unit flap deflection angle { q w } (o) , and the second component is the deformed airfoil without a flap deflection { q w } (1) . Therefore, the flap aerodynamic force vectors are as follows:

Membrane airfoils
The aeroelastic formulations and codes are verified by considering a membrane airfoil with a chord (c) length of 1 m and a mass per unit ( ) of 0.5 kg/m. The membrane is subjected to a free stream with air density (ρ ∞ ) of 1.225 kg/m 3 and speed (V ∞ ) of 15 m/s at an angle of attack (α) of 3 degrees and an applied tension force coefficient (C T T /q ∞ c) of 3. For the Rayleigh-Ritz approach, 20 terms were considered for the membrane deformation, while for the FEM 50 elements were considered for the study and 40 aerodynamic terms were used for the Prandtl-Glauert series expansion. The static deformation shape of the membrane airfoil is shown in Fig. 5, and the values of the corresponding aerodynamic coefficients (the lift coefficient (C L ), the lift-alpha curve slope (C Lα ), the moment coefficient at the quarter chord (C M ), the leading-edge moment coefficient (C M L E ), the center of pressure (x cp ), the maximum camber (z cmax ), and the divergence speed (V D )) are presented in Table 1, and it can be seen that both the Rayleigh-Ritz and the FEM are in a good agreement with results of Nielsen [8].
Modal analysis was used for the dynamic analysis of the membrane airfoil. The first four mode shapes obtained from the FEM are shown in Fig. 6. Only heavy membranes exhibit flutter instability, so the membrane density ratio ( /ρ ∞ c) is set to 25 and the root loci of the dynamic system as the free stream speed increases from 0 to 18 m/s are computed as shown in Fig. 7. It can be seen that the system becomes dynamically unstable at speed of 14.5 m/s which is corresponding to a tension coefficient of 3.2 and it is the value reported by Tiomkin [17]. All modes exhibit reduction in their frequencies as the speed increases due to the loss of stiffness caused by the aerodynamic loading. Regarding the computational cost of each method, the two codes were written in Matlab and run on a Intel(R) Core(TM) i7-10870H CPU @ 2.20 GHz Machine. The time required to get the final dynamic aeroelastic system using the Rayleigh-Ritz (20 terms) method and the finite element method (50 elements) was 0.61 and 10.5 s, respectively, which is a significant difference, but as the geometrical complexity of the problem increases, the Rayleigh-Ritz will no longer be suitable. The running time of the finite element  Table 2 shows a parametric analysis for different flap lengths and combinations of the angle of attack and the flap deflection angle. Figure 8 shows the deformation shape of the membrane airfoil with/without flap deflection at a 3°angle of attack. The maximum camber shifts toward the LE for an upward flap deflection and vice versa for a downward flap deflection. Also, an inflection point appears near the TE for an upward flap deflection. Figures 9 and 10 show the variation of the lift and moment coefficient with the TE flap length for different tension coefficient values. Both the lift and moment coefficients exhibit nonlinear relations in terms of the flap length and the tension coefficient. The moment shows stronger nonlinearity compared to the lift and as the tension coefficient increases, this nonlinearity decreases.
The dynamic aeroelastic characteristics of the membrane airfoil are studied by considering the transient response due to a step and harmonic flap deflections using eight natural mode shapes. Figures 11 and 12 show the step response of the lift and moment coefficients, while Figs. 13, 14, 15, and 16 show the first four modes step responses. The higher modes have low aerodynamic damping leading to the superharmonics in the coefficients' step responses which can be eliminated by adding structural damping. The frequency response due to harmonic TE flap deflection is shown in Figs. 17, 18, and 19. Initially, as the reduced frequency increases, the lift coefficient amplitude decreases due to the influence of the wake; then, it increases as the reduced frequency approaches the aeroelastic natural frequency of the membrane. As the flap length and the tension coefficient increase, the peak lift amplitude increases relative to static values.

Flexible-chord airfoils
This section investigates the aeroelastic characteristics of deformable airfoils. The considered airfoil has a unit chord length and mass per unit length ( ) of 4 kg/m. The free stream flow has an air density (ρ ∞ ) of 1.225 kg/m 3 and speed (V ∞ ) of 15 m/s at an angle of attack (α) of 3 degrees. The airfoil flexibility is represented by the nondimensional bending stiffness (E I E I /q ∞ c 4 ) which is set to vary from 0.5 to infinity. The torsional and translational spring constants (k ∅ , k w ) are set to 800 N.m/rad. For the dynamic analysis, the translational spring constant is set to be equal to the torsional spring constant (k w k ∅ ) and for the static analysis k w 10k ∅ as it has no contribution to the results which will be more appropriate for deformation demonstration purposes.  the peak value at the first bending mode aeroelastic natural frequency. The peak amplitude increases as the springs position moves forward.

Conclusion
The aeroelastic performances of membrane airfoils and flexible-chord airfoils with TE flaps are studied in this paper. Two aeroelastic formulations are presented; the Rayleigh-Ritz formulation and the finite element formulation. The two presented formulations are very efficient as the fluid-structure interactions are presented in closed-form, easy to compute, aerodynamic matrices. The finite element formulation has the advantage of modeling complex geometries without any adjustments to its formulation, so it is more suitable to model flapped airfoils. The analysis of both airfoil types showed that the aerodynamic moment coefficient is more sensitive to the flexibility level and to the TE flap length. For flexible-chord airfoils, the position of the elastic axis (the rigid-body springs position) plays a major role in the effect of the flexibility on the aeroelastic characteristics. The flutter speed decreases with increased airfoil flexibility, except for certain elastic axis positions of flexiblechord airfoils where the flexibility increases the flutter speed. For harmonic TE flap input, the peak amplitude

Conflict of interest The author declares no conflict of interests of any kind.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.
Finite element formulation