Unsteady Lifting Line Theory Using the Wagner Function for the Aerodynamic and Aeroelastic Modeling of 3D Wings

A method is presented to model the incompressible, attached, unsteady lift and pitching moment acting on a thin three-dimensional wing in the time domain. The model is based on the combination of Wagner theory and lifting line theory through the unsteady Kutta–Joukowski theorem. The results are a set of closed-form linear ordinary differential equations that can be solved analytically or using a Runge–Kutta–Fehlberg algorithm. The method is validated against numerical predictions from an unsteady vortex lattice method for rectangular and tapered wings undergoing step or oscillatory changes in plunge or pitch. Further validation is demonstrated on an aeroelastic test case of a rigid rectangular finite wing with pitch and plunge degrees of freedom.


Introduction
Closed-form solutions for the attached incompressible unsteady flow problem around a two-Dimensional (2D) airfoil exist in both the frequency domain [1] and in the time domain: • Wagner theory [2,3] • Finite state flow model [4] • Leishman unsteady state space representation [5] For three-Dimensional (3D) wings, there exists one closed-form solution for the unsteady aerodynamics of elliptical wings [2]. For general geometries, closed-form solutions are usually obtained either from strip theory (see for example Dowell [6]) or from panel methods, such as the Doublet Lattice Method (DLM) [7] or the Vortex Lattice Method (VLM) [8]. Strip theory is based on estimating the 3D unsteady loads by integrating 2D loads along the span. It therefore ignores the downwash induced by the trailing vortices and overestimates the lift; it is mostly used on lifting surfaces with very high aspect ratios, such as helicopter or wind turbine blades. The DLM can be used to estimate modal, frequency domain aerodynamic loads in the form of the generalized aerodynamic force matrix. This matrix is evaluated numerically at discrete reduced frequency values and is interpolated in order to estimate the aerodynamic loads at intermediate frequencies.
The generalized force matrix can be transformed to the time domain using the rational function, the Chebyshev polynomial or indicial function representations, again based on a set of discrete frequency estimations. Several efficient transformation methodologies have been developed, notably the minimum state approach [9], but they remain approximations. The Vortex Lattice Method can also be used to derive a generalized aerodynamic force matrix [10,11], which can then be transformed to the time domain.
Alternative closed-form solutions of the 3D attached flow problem include Reissner's method [12], which combined the Theodorsen and lifting line theories and was formulated in the frequency domain. Chopra [13] developed expressions for the lift, thrust and moment of lunate tails oscillating in pitch and plunge, based on the lifting line assumption that the flow is locally 2D around cross-sections of the wing's span, but that the local angle of attack is influenced by the vorticity in the wake. The work was limited to rectangular wings, and frequency domain solutions for the aerodynamic loads were obtained. Furthermore, Chopra and Kambe [14] formulated an unsteady lifting surface theory and applied it to wings with non-rectangular planforms, but the aerodynamic load calculations were still calculated in the frequency domain. James [15] developed an unsteady lifting line theory for wings of a large aspect ratio with smooth tip geometries (such as elliptical planforms). He used a matched asymptotic expansion approach to obtain solutions for impulsively-started motion, constant acceleration and sinusoidal oscillations in pitch, plunge or flap. Nevertheless, the method failed to yield total aerodynamic loads for wings with chords that jump abruptly to zero at the tip (such as rectangular or trapezoidal wings). Furthermore, Ahmadi and Widnall [16] argued that James's theory is only valid for low reduced frequencies and that its 3D results are incorrect. Van Holten [17] also used a matched asymptotic expansion to develop an unsteady lifting line theory for pitching wings and rotating blades, but Ahmadi and Widnall [16] again claimed that the work is only valid for low reduced frequencies. Phlips et al. [18] derived a time domain unsteady lifting line theory for flapping (but not pitching) wings.
Other frequency domain unsteady lifting line approaches were proposed by Dragos [19], Sclavounos [20] and, more recently, Drela [21]. State-space time domain models are usually quasi-steady, such as the models by Nabawy and Crowther [22,23]. The present paper details a robust, closed-form, time domain, 3D unsteady aerodynamic model that does not involve a transformation from the frequency domain. The approach is based on Wagner's 2D unsteady lift theory and Prandtl's lifting line theory and will be referred to as the Wagner Lifting Line (WLL) method. It was first proposed by Boutet and Dimitriadis [24], but is presented here in much more detail, including an aeroelastic extension. A similar approach was proposed slightly later by Izraelevitz et al. [25], but this technique was not extended to aeroelastic applications.

Method
As mentioned in the Introduction, one of the main characteristics of lifting line theories is that the flow is two-dimensional around spanwise cross-sections, but the local angle of attack is affected by the downwash induced by the wake. In the classical version of the theory, the wing and wake are modeled using a superposition of horseshoe vortices, the strength of which is constant in time. In this way, the wake is straight and semi-infinite, and its strength varies in the spanwise, but not in the chordwise direction.
In Wagner's and Theodorsen's 2D unsteady aerodynamic theories, the wake is still straight and semi-infinite, but its strength varies in the chordwise direction since it is calculated from the unsteady Kutta condition. The wake is assumed to propagate with the free stream velocity, U, so that chordwise displacement and time are directly related by the equation x = Ut. A change in vorticity at the trailing edge that occurs at time t 0 will be reflected in the wake at a downstream distance U(t − t 0 ) at time t.
Lifting line and Wagner theories are not directly compatible because in the former, the chordwise strength of the trailing vortices is constant, while in the latter, it varies. In the present work, Wagner theory is applied to spanwise cross-sections, so that the strength of the wake varies in both the spanwise and chordwise direction. A quasi-steady version of lifting line theory is used in order to approximate the downwash velocity caused by the wake and to add it to the other sources of downwash used in Wagner theory.
There are three sources of downwash on finite wings in unsteady flow: • Geometric downwash due to camber and twist. • Downwash due to the motion of the wing (including the free stream and angle of attack).
• Downwash due to the three-dimensionality of the flow (wing tip vortex effect).
The 3D downwash will be calculated using Prandtl's lifting line theory, as detailed in Kuethe and Chow [26]. The other two downwash contributions will be modeled from Wagner's unsteady aerodynamic theory, as presented by Fung [3].

Lifting Line Theory
A truncated Fourier series is used to represent an arbitrary time-varying circulation distribution along the span of a flat plate wing: where a 0 is the lift curve slope at the wing's axis of symmetry, c 0 is the chord length at the wing's axis of symmetry, U is the free stream airspeed, a n (t) are the time varying Fourier coefficients, m is the number of terms kept in the series, θ comes from the substitution y = (s/2) cos(θ), y represents the location along the span and s represents the span of the wing. In the classical lifting line theory, the wing is split into m spanwise strips, and the order of the Fourier series is also m, so that the number of unknowns (Fourier coefficients a n ) is equal to the number of equations (spanwise strips). Izraelevitz et al. [25] proposed an alternative approach, whereby a horseshoe vortex representation was used to model the 3D wing circulation. Using Prandtl's lifting line theory, it is possible to compute the downwash w y caused by the 3D circulation distribution at a location y along the wing span: Glauert [27] evaluated the integral in this latest expression, so that the downwash w y can be computed as a function of the Fourier coefficients and the location y(θ) along the span: Note that this is a quasi-steady version of lifting line theory, since any instantaneous changes in vorticity over the wing affect the entire wake simultaneously.

Unsteady Kutta-Joukowski
It is possible to express the unsteady sectional lift coefficient as a function of a n (t) and location along the span y, using the unsteady Kutta-Joukowski theorem and considering a lumped spanwise vortex element, as explained by Katz and Plotkin [8] on page 439. The circulatory sectional lift coefficient becomes: where Γ is the vortex strength, c(y) is the chord at span section y and the unsteady termΓ comes from the unsteady Bernoulli equation. The vortex strength can be replaced by its Fourier series representation from Equation (1), to obtain: c 0 c a n + c 0 Uȧ n sin(nθ) Furthermore, the circulatory lift coefficient for the entire wing can be computed from: where S is the wing's surface area. The lift will cause a pitching moment around the pitch axis of each wing section, i.e., the axis around which the wing section can pitch. One can compute the circulatory moment distribution from the lift distribution, assuming that the sectional lift acts on the quarter chord. Refer to Figure 1, which shows a wing section with chord c and half-chord b at pitch angle α to a free stream U. The position of the pitch axis, x e , is defined with respect to the half-chord point. The circulatory sectional moment equation is simply: where x e and c are allowed to vary in the spanwise direction y. Note that the pitch axis is measured from the mid-chord point and is defined as positive if it lies downstream of that point as chosen by Theodorsen [1]. The total circulatory moment coefficient is:

Wagner's Sectional Circulatory Lift
The computation of the unsteady circulatory aerodynamic loads is based on the circulatory sectional lift c c l (t, y) response of an airfoil undergoing a step change in downwash ∆w(y) << U at span location y. The resulting step change in the lift coefficient can be expressed in terms of the Wagner function, Φ(t), as follows: where a 0 (y) is the lift curve slope of the local airfoil section, which is approximated by 2π for thin airfoils, and Φ(t) is Jones' [28] approximation of the Wagner function: with Ψ 1 = 0.165, Ψ 2 = 0.335, 1 = 0.0455 and 2 = 0.3.
Duhamel's principle can be applied to Equation (9) in order to express a continuous lift response as the time integral of infinitesimal step responses: The troublesome ∂w(τ,y) ∂τ term inside the integral can be removed by applying integration by parts, such that: The downwash w(t, y) must now be computed; it will depend on the kinematics of the wing. In this work, the wing is assumed to be rigid with pitch and plunge degrees of freedom, but flexible wings with bending and torsion modes can also be considered. Figure 1 defines the local plunge and pitch degrees of freedom, h(t, y) and α(t, y), respectively. A local downwash w y is added in order to represent the 3D downwash effects, so that w(t, y) becomes: where d is the non-dimensional distance between the mid-chord and the pitch axis, as defined by Theodorsen [1]. After combining Equations (11) and (12), The following changes of variables are performed in order to eliminate the integrals from Equation (13), where the variables z k (t) are called aerodynamics states.
Working through the integrals, we obtain: First order differential equations for the aerodynamic state variables can be derived using Leibniz's integral rule. As an example, the equation for z 1 (t, y) is: The equations for all the aerodynamic states are given by: Finally, the continuous unsteady circulatory lift coefficient can be written as follows: where:

Added Mass Effect
The non-circulatory sectional lift and moment coefficients, also known as the added mass effect, can be computed from the non-circulatory terms derived by Theodorsen [1]: The computation of the total non-circulatory lift C i l and moment C i m coefficients follows the same principle as Equations (6) and (8), and the final result is given in Equations (A26) and (A27) for a rectangular wing.

Assembling the Pieces Together
In order to compute the m Fourier coefficients a n , m spanwise wing strips are considered, as shown in Figure 2, which represents m strips along the wing span, their respective local aerodynamics variables z i , 3D downwash effects w y,i and chord section c i . For an arbitrary strip i, Equations (1), (18) and (19) can be combined to obtain the following system: The variables z i represent the local aerodynamic state variables for the i-th strip. Matrices C i , D i , E i and W i are the matrices computed in Section 2.3 where the chord c is replaced with its local value c i . Variable w y i represents the 3D downwash effect on the i-th strip; its value is given by Equation (3). Substituting for w y i and re-arranging, Equation (22) become: Applying Equation (23) to all m strips, a set of 7m ordinary linear differential equations for 7m unknowns (m for the Fourier coefficients a n and 6m for the aerodynamic states z 1 , . . . z m ) is obtained. Once this system of ODEs is solved and the coefficients a n (t) are evaluated, the lift distribution acting on the wing can be computed from Relation (5), the total lift from (6) and the total pitching moment from (8). As an example, the computation of the aerodynamic loads for a rectangular wing is fully detailed in the Appendix, Equations (A1)-(A29).
2.6. Asymptotic Behavior for a Rectangular Wing (c = c 0 ) By imposing stationarity, all derivatives in time vanish and the system of Equation (23) should reduce to the classical lifting line theory. After imposing stationarity on Equation (17), the aerodynamic states variables become: Injecting these expressions in Equation (23) and applying the necessary simplifications, the system becomes: which is the equation derived by Kuethe and Chow [26] for a rectangular wing, using Prandtl's lifting line theory for a constant pitch angle α(y) = α.
If the span of the wing is infinite, s = ∞, all 3D effects vanish and System (23) should reduce to the 2D Wagner solution. Looking at Equation (3), it is obvious that lim s→∞ w y = 0. It can then be shown from Equation (14) that z 5 = 0 and z 6 = 0. The system of Equation (23) is now reduced to: which is the classical 2D Wagner formulation, as given by Fung [3] for all arbitrary strips i.

Test Cases
The Wagner lifting line method is applied here to a rectangular and a tapered wing. The rectangular wing, shown in Figure 3, has constant chord c = 1 m, span s, surface S = cs and aspect ratio AR = s 2 /S. It has two degrees of freedom, a pure plunge h(t, y) = h(t) and a pure pitch angle α(t, y) = α(t) around its pitch axis. Two positions of the pitch axis are considered: one at the leading edge and one at the quarter chord.

Types of Motion
The degrees of freedom of the wings are subjected to two kinds of motion: step changes and sinusoidal oscillations. The step is expressed as the function f = A 1 − e −10t where t is the time and A is the position of the degree of freedom after the step. Steps ∆h = 0.1 m and ∆α = 5 • are separately applied to the plunge and pitch degrees of freedom.
Sinusoidal oscillations with ten distinct frequencies are tested separately for each degree of freedom, in order to assess the WLL model's frequency response. The oscillations are expressed as

Validation
The unsteady Vortex Lattice Method (VLM) is used as a reference solution to which the results obtained by the Wagner lifting line approach are to be compared. The particular implementation of the VLM used here is more thoroughly described by Dimitriadis et al. [29]. The difference between the solutions obtained from WLL and VLM is quantified using the Normalized Root-Mean-Square Deviation (NRMSD).
[%] (27) where N is the number of time instances, y v,r represents the lift or moment coefficient computed by the VLM at the r-th time instance and y w,r represents the lift or moment coefficient computed by the WLL approach at the r-th time instance.

Convergence
The differential Equation (23) are solved by means of a Runge-Kutta-Fehlberg fourth and fifth order numerical time integration technique. A time convergence analysis is therefore needed in order to minimize numerical integration errors without increasing the computation time too much. A convergence analysis is also performed for the vortex lattice method.
As shown in Equation (1), the WLL uses a truncated Fourier series with m coefficients, which correspond to the m spanwise strips shown in Figure 2. A spatial convergence study must also be carried out in order to determine the effect of the number of strips on the aerodynamic load predictions. Note that the number of states in the system is 7m; therefore, keeping m as low as possible is important. The NRMSD is used in order to determine if convergence has been achieved, such that: where y i,t represents the lift coefficient response at the t-th time instance for the i-th value of the convergence parameter. The latter can be either the time step tolerance of the Runge-Kutta-Fehlberg scheme or the number of strips; y re f represents the reference lift coefficient against which the convergence is assessed. This reference is computed for an appropriately high value of the convergence parameter.

Runge-Kutta Convergence
The time step tolerance controls the error of the Runge-Kutta time integration scheme. Given the solution arrays r 4 and r 5 of the respective fourth and fifth order Runge-Kutta estimates for a current time t i and time step ∆t, the Runge-Kutta-Fehlberg algorithm used in this work is given by: • if Tol < (r 5 − r 4 )(r 5 − r 4 ) T , reduce the time step to ∆t/2 • otherwise, go to the next time instance t i+1 = t i + ∆t and reset ∆t to its default value Figure 5 plots the variation of the error of Equation (29) with the tolerance Tol used in the Runge-Kutta-Fehlberg algorithm with respect to a reference Tol = 10 −9 . The solid line represents the convergence for a step case and the dashed line the convergence for a sinusoidal oscillation case with a reduced frequency k = 0.3. Assuming that e < 0.01% is good accuracy, the figure shows that a tolerance value of 10 −7 is sufficient to achieve convergence for both the step and sinusoidal motion. Step AR 6 Step AR 12 Step AR 18 Oscillation AR 6 Oscillation AR 12 Oscillation AR 18 Step AR 6 Step AR 12 Step AR 18 Oscillation AR 6 Oscillation AR 12 Oscillation AR 18 (b)   (29) with the number of strips for different aspect ratios and kinematics with respect to a reference number of strips m = 26. The time step tolerance is fixed to 10 −7 . It can be seen that, for all cases, a value of m = 20 is sufficient to achieve errors e < 0.01%. A still acceptable error of e < 0.1% can be achieved with m = 10 strips. The figure also shows that the higher the AR, the greater the number of strips necessary to reach the same level of convergence. As the aspect ratio increases, the spanwise lift distribution becomes flatter and decreases more quickly at the wingtips; the Fourier series of Equation (1) requires more terms in order to represent such lift distributions accurately. In the asymptotic case, a rectangular wing with infinite span can only be modeled if an infinite number of terms is used in the series; we have already shown in Equation (25) that the WLL model then reduces to the classical 2D Wagner formulation.

Vortex Lattice Convergence
The quality of the solutions obtained from the vortex lattice method is based on the the number of panels used to represent the wing in the spanwise and chordwise directions. The VLM is more sensitive to the chordwise than the spanwise number of panels [30], so a convergence for the number of chordwise panels is performed. Figure 6 plots the variation of the error of Equation (29) with the number of chordwise panels for a rectangular wing (AR = 6) oscillating around its leading edge with a reduced frequency k = 0.3. The reference number of chordwise panels is j = 100. The number of spanwise panels is set to 15. It can be seen that j = 75 is sufficient to achieve errors e < 10 −1 % for the rectangular wing case with a leading edge pitch axis. The figure also shows that the convergence is independent of the aspect ratio. This is logical since the convergence parameter is the number of chordwise panels.
All other test cases used in this work converge for j = 75, as long as the pitch axis does not lie at c 0 /4. A larger number of chordwise panels j = 100 is necessary for convergence when the pitch axis lies on the quarter chord because the moment loads need to converge to values very close to zero for the oscillating case with low reduced frequency. Convergence is therefore harder to achieve.

Lift and Moment Results for a Rectangular Wing
In this section, the WLL method is applied to the rectangular wing under different kinematic conditions. The lift and moment estimates are compared to predictions obtained by the VLM technique and by strip theory. In all cases, the VLM estimates were obtained from time-converged and spatially-converged simulations.
For the first comparison, the pitch axis is located at the leading edge and the wing undergoes a step change in pitch or plunge, as explained in Section 3.1. The resulting lift and moment responses are shown in Figure 7; the WLL estimates are in good agreement with the VLM results for both pitch and plunge motions. In contrast, the strip theory result is significantly overestimated in the pitch step case and underestimated in the plunge step case. Note that the agreement in pitching moment for the step pitch case is not as good as in the other results; small differences between the VLM and WLL predictions persist at steady-state conditions. This is due to modeling differences between the two approaches. WLL calculates the total stripwise lift and places it at the quarter chord. In contrast, the VLM calculates a lift force on the 3/4 chord point of each chordwise panel, and the total stripwise lift is the sum of all the lift forces on the same strip. The point of application of the VLM's total stripwise lift is exactly the quarter chord only for an infinite number of chordwise panels. Simulations where repeated for other positions of the pitch axis, but the results are not shown here for conciseness. All of these simulations resulted in WLL predictions that were in good agreement with the VLM estimates.
For the second comparison, the WLL, VLM and strip theory techniques are applied to a rectangular wing with an aspect ratio of six undergoing sinusoidal oscillations in pitch or plunge, as detailed in Section 3.1. Several reduced frequency values were tested, but only the results for k = 0.1 and k = 0.3 are presented here. The lift results are plotted in Figure 8 and the moment results in Figure 9. There is very good agreement between the WLL and VLM predictions for all frequencies, while the strip theory estimates are quite inaccurate at both values of k.  Figure 10a plots the NRMSD values calculated from Equation (27), for all the tested values of k. The difference between the WLL and VLM predictions increases with reduced frequency, but stays lower than 3% for all kinematics and aerodynamic loads. This level of difference is considered good for such frequencies. Repeating the simulations after moving the pitch axis to the quarter chord, as shown in Figure 10b, results in equally good agreement between the WLL and VLM aerodynamic load predictions.

Lift and Moment Results for a Tapered Wing
In this section, step and sinusoidal numerical tests are applied to the tapered wing with an aspect ratio of six described in Section 3. Initially, the pitch axis is located at the leading edge of the root chord, and the wing undergoes step changes in pitch or in plunge, as detailed in Section 3.1. The resulting lift and moment responses are shown in Figure 11. In all cases, the WLL and VLM predictions are in very good agreement, while the strip theory results are highly overestimated in the pitch step case. Moving the pitch axis to the quarter chord results in equally good agreement between the WLL and VLM aerodynamic load responses. Finally, the WLL and VLM are compared for the case of pitch or plunge sinusoidal oscillations at different reduced frequencies for the tapered wing. The pitch axis is located at the leading edge of the root chord. Figure 12a plots the variation of the NRMSD values between the two sets of predictions, for increasing k values. It can be seen that in all cases the maximum NRMSD stays below 5%. It is concluded that the WLL approach can predict accurately the aerodynamic load responses for oscillating tapered wings. Moving the pitch axis to the quarter chord results in equally good predictions, as seen in Figure 12b.

Computational Cost
The computation times of the WLL and VLM approaches are compared for a step and three sinusoidal oscillations in pitch, with reduced frequencies k = [0.1, 0.5, 1] and a total simulation time T f = 1.3 s. The following parameters have been used for each method : VLM: There are 15 spanwise and 20 chordwise panels, and the time step is ∆t = 10 −3 . The wake shape is prescribed in order to reduce computational cost. WLL: The computation parameters used are 15 strips, a tolerance Tol = 10 −8 and an initial time step ∆t = 10 −3 for the Runge-Kutta-Fehlberg 45 algorithm.
All the calculations were run on a MacBook Pro laptop computer with a quad-core Intel i7 processor rated at 2.5 GHz, running iOS Version 10.9.5. Table 1 shows that WLL calculation times are much lower than VLM numerical simulations to simulate 1.3 s of pitch step or sinusoidal oscillations with reduced frequency k = [0.1, 0.5, 1]. The benefit of an adaptive time step for WLL can also be seen as the computation time is lower for cases with low or high reduced frequencies. It should be stressed that the implementation of the two methods is completely different; the WLL is implemented purely in MATLAB, while the VLM is implemented as a combination of MATLAB and C code (mex functions). If the WLL were also implemented using compiled code, it would be even faster.

Aeroelastic Test Case
In order to further validate the Wagner lifting line approach, an aeroelastic test case is presented for a rigid rectangular wing with pitch and plunge degrees of freedom. The flutter speed will be computed as a function of the position of the pitching axis and the wing's aspect ratio using both the Wagner lifting line and the Vortex lattice method.

Aeroelastic Equations of Motion
The structure simply consists of a rigid wing with two degrees of freedom in pitch and plunge. The structural equations of motion are given by: where q = [h α] T , h is the plunge displacement, α the pitch displacement, m w the mass, S w the static imbalance around the pitching axis, I w the moment of inertia around the pitching axis, k h and k α the stiffnesses of the springs providing restoring loads in the plunge and pitch degrees of freedom, respectively, and L(t) and M(t) are the lift and moment around the pitching axis computed using the VLM or WLL approaches. The wing is chosen as an aluminum rectangular flat plate with chord c 0 = 0.1 m, thickness h = 0.005 m and span s; the distance between the pitch axis and the mid-chord is x e . The mass matrix components can then be computed as: where ρ al = 2300 kg/m 3 is the density of aluminum. The spring stiffnesses for the two degrees of freedom are chosen such that the uncoupled, wind-off natural frequencies of the system are f h = 1 Hz and f α = 5 Hz. The stiffnesses are then given by: The normal force L(t) and moment around the pitching axis M(t) are computed from Equations (6) and (8), together with the added mass effects described in Expressions (20) and (21). The integrals are approximated using the trapezoidal rule.
Finally, the complete linear aeroelastic system composed of Equations (30) and (22) can be written in first-order form as:ẋ The aeroelastic system matrix H is derived in Appendix A, while x represents the system states and is defined as: where a i is the i-th Fourier coefficient, z i are the local aerodynamic state variables for the i-th strip: and m is the number of strips. Consequently, the total number of states is 7m + 4.
Finally, the WLL flutter solution is obtained by computing the eigenvalues of matrix H(U, x e ) as a function of the airspeed U. An indirect search procedure is employed to pinpoint the airspeed at which one pair of complex conjugate eigenvalues becomes purely imaginary, which is the definition of the flutter speed. The number of strips used to estimate matrix H(U, x e ) is m = 20, and the flutter solution is obtained for several values of the position of the pitch axis x e .
The VLM flutter solution is obtained using the modal frequency domain version of the method, as detailed in [11]. Rigid body modes are chosen, one for the plunge and one for the pitch. The mode shapes are given by: where w h (x, y) is the plunge mode shape and w α (x, y) the pitch mode shape. The elements of the mass matrix are then obtained from: The resulting flutter problem is of the form: where A s is the structural mass matrix, E s is the structural stiffness matrix and Q(k) is the frequency-dependent generalized aerodynamic force matrix generated by the VLM approach. The flutter problem is solved using the p − k method.

Results
The resulting flutter speed and frequency values are plotted against the position of the pitch axis in Figure 13 for two aspect ratios: AR = 4 and AR = 10. Figure 13a plots the flutter airspeeds and shows that the VLM and WLL predictions are in good agreement with each other for both aspect ratios. Figure 13b plots the flutter frequency predictions; the agreement is still very good as the highest frequency discrepancy is of the order of 5%. It can be concluded that the Wagner lifting line method can predict accurately the flutter of a wing with a finite span.

Conclusions
The WLL method results in a closed-form, state-space representation of the unsteady aerodynamic loads acting on finite rectangular and tapered wings of different aspect ratios, under attached incompressible flow conditions. The technique combines Wagner's 2D unsteady lift theory, Prandtl's lifting line theory, the unsteady Kutta-Joukowski theorem and the added mass terms from Theodorsen's analysis. Sample simulations on wings with and without taper have shown very good agreement between the WLL predictions and VLM simulation results. The method can also be readily applied to wings with twist and camber. Sweep is more problematic, since lifting line theory has to be modified in order to work in the presence of sweep. This modification will be addressed in future work.
The VLM approach is still more general than the WLL technique, as it can easily represent sweep. The advantage of WLL is the fact that the resulting aerodynamic loads are written in state space form, as functions of the structural and aerodynamic states. They can therefore be easily included in aeroelastic and flight dynamic calculations. In contrast, the VLM or DLM techniques result in time-marching simulations or, if using a modal frequency domain technique, in Equation (33), which is a hybrid time-frequency domain equation that must be transformed to the time domain in order to carry out aeroservoelastic calculations. In the present examples, the wings were rigid with discrete degrees of freedom, but flexible wings with generalized modes can also be treated. Finally, the WLL calculations are significantly faster than time domain VLM numerical simulations.
Using Expression (A2), one can write: L c (t) = TA y i a n + TD y iȧ n (A24) As the chord and the location of the pitch axis are constant with span for a rectangular wing, the circulatory moment in pitch M c is simply:

. Non-Circulatory Loads
It is assumed that the wing is rectangular; therefore, the non-circulatory load coefficients are identical to the sectional lift coefficients described in Section 2.4. They are written as: Appendix A.6. Total Aerodynamic Loads The total lift applied on the wing can be written as: = TA y i a n + TD y iȧ n + A lq + B lq (A28) Similarly, the total pitching moment is: = (c/4 + x e )TA y i a n + (c/4 + x e )TD y iȧ n + A mq + B mq (A29)

Appendix A.7. Structural Equations
From Expressions (A28) and (A29), Equations (30) can be written as: m w S w + A l q + TD y iȧ n = −TA y i a n − B lq − k h 0 q (A30) S w I w + A m q + (c/4 + x e )TD y iȧ n = −(c/4 + x e )TA y i a n − B mq − 0 k α q (A31)