Unsteady lifting-line theory and the influence of wake vorticity on aerodynamic loads

Frequency-domain unsteady lifting-line theory (ULLT) provides a means by which the aerodynamics of oscillating wings may be studied at low computational cost without neglecting the interacting effects of aspect ratio and oscillation frequency. Renewed interest in the method has drawn attention to several uncertainties however. Firstly, to what extent is ULLT practically useful for rectangular wings, despite theoretical limitations? And secondly, to what extent is a complicated wake model needed in the outer solution for good accuracy? This paper aims to answer these questions by presenting a complete ULLT based on the work of Sclavounos, along with a novel ULLT that considers only the streamwise vorticity and a Prandtl-like pseudosteady ULLT. These are compared to Euler CFD for cases of rectangular wings at multiple aspect ratios and oscillation frequencies. The results of this work establish ULLT as a low computational cost model capable of accounting for interacting finite-wing and oscillation frequency effects and identify the aspect ratio and frequency regimes where the three ULLTs are most accurate. This research paves the way towards the construction of time-domain or numerical ULLTs which may be augmented to account for nonlinearities such as flow separation.


Nomenclature
The unsteady aerodynamics of wings is an increasingly important topic in modern aerospace research. On the smallest of scales, engineers are using flapping wing motions in the design of micro air vehicles (MAVs). At larger scales, oscillatory motion is seen as a means by which energy can be captured from water currents with reduced environmental damage. 56 And at larger scales still, high aspect ratio aircraft 44,68 and wind turbine designers 29 face challenges is understanding the consequences of wings flexing under load. The focus of this article is on finite-wing effects in unsteady aerodynamics. Early research in this topic was strongly motivated by dynamic stall in helicopters. 39 Studies in this topic, carried out using Computational Fluid Dynamics (CFD) and experiments are reviewed in Carr, 12 Carr et al. 13 and Ekaterinaris and Platzer. 21 More recent investigations (in the 21st century) into three-dimensional effects in dynamic stall using state-of-the-art experimental and numerical methods include Angulo et al. 2 and Visbal and Garmann. 64,65 All these studies concern the forces and flow around a finite aspect ratio wing undergoing pitch-plunge motion in a high Reynolds number, low reduced frequency regime. The modern aerospace applications listed earlier have inspired unsteady aerodynamics research in new regimes, at low/transitional Reynolds numbers and involving high reduced frequencies. These flows are typically massively separated and involve leading-edge vortex (LEV) formation and shedding. Finite-wing effects in these regimes, specifically concerning the interactions of tip vortices with leading-and trailing-edge vortices have been investigated using CFD and experiments, for translating wings, 20,38,43 pitching wings, 26,32,34,63,70,72 plunging wings, 8-10, 23, 66, 71 rotating wings 3,14,15,41,47,62 and wings subject to gusts. 4,17,48 These studies have shed much light on force production in unsteady finite-wings and on the contributions from circulatory, apparent-mass and vortical effects.
Theoretical or low-order methods in aerodynamics are required for rapid design analysis and optimization, and, later on in the engineering cycle, real-time simulation and use in control systems. The development of these is aided and inspired by experimental and computational studies, such as those referenced earlier. In 2D, Theodorsen's theory 60 provides the lift and moment for a thin, symmetric airfoil undergoing small harmonic pitch and plunge oscillations in potential flow. This theory has been shown to provide good results even in regimes outside the expected limit of validity such as at low Reynolds numbers, for large amplitudes of oscillation, and even when LEVs are present (so long as they are attached). 40,45 2D low-order models based on the discrete-vortex method have been used to successfully predict airfoils undergoing arbitrarily large kinematics and massively-separated flows. 52,53 In 3D, the simplest approach is strip theory. 37 Neglecting all three dimensional effects (for large aspect ratios) allows two dimensional methods to be applied to the wing in strips. In instances where this is a poor approximation, low-order 3D methods such as the unsteady boundary-element method, 42 unsteady vortex lattice method 30,36,44,58 or vortex particle methods 69 can be used, but at far greater computational cost due to the large numbers of interacting vortex elements required to form wake vortex structures. Eldredge and Darakananda 22 suggest forming these vortex structures with a very small number of time-varying elements instead, but this introduces difficulties. As with similar 2D models, the lack of natural vortex shedding mechanisms makes such an approach dynamically incomplete for long-term simulations, according to Darakananda and Eldredge. 18 There is however an intermediate to 2D and 3D methods. By augmenting 2D models with corrections gathered from a much simplified 3D model, lifting-line theory may be obtained. By avoiding full three dimensional interaction, the computational complexity of the problem is significantly reduced. Lifting-line theory assumes a high aspect ratio wing. The chord scale and span scale are so far removed that the problem can effectively be separated into detailed inner 2D solutions which interact via a simplified outer 3D solution.
The simplest form of lifting-line theory assumes straight wings and steady flow. Prandtl 50 is credited with the earliest lifting-line theory. This was formalized through matched asymptotic expansions by Van Dyke 61 and extended to curved and swept wings by Guermond. 27 Non-linear extensions to steady lifting-line theory have become an important design tool. 24 The desire for an unsteady lifting-line theory (ULLT) lead to work on the problem of oscillating wings. Such a problem is the three-dimensional analogue of the oscillating plate, analyzed by Theodorsen. 60 Assuming a high aspect ratio wing, Cheng 16 has identified five frequency ranges based on the wake wavelength λ , span (2s) and average chord (c): (i) very low frequency (2s << λ ), (ii) low frequency (2s = O(λ )), (iii) intermediate frequency (c << λ << 2s), (iv) high frequency (c = O(λ )), and very high frequencies (c >> λ ). Early authors in this topic such as James 33 and van Holten 31 assumed a uniform induced downwash over the chord (as in the steady case) for all frequencies. Ahmadi and Widnall 1 showed that this assumption was asymptotically valid only for very low frequencies, and further showed that in the low-frequency regime the induced downwash has a harmonic variation in the chordwise direction. Guermond and Sellier 28 showed that this conclusion is actually true for all frequencies, and also derived a theory applicable for wings with curvature and sweep and uniformly valid for all frequencies. Sclavounos 57 produced a ULLT that assumes a uniform induced downwash but is claimed to be valid at all but very high frequencies.
All these ULLTs assume harmonic kinematics, small amplitudes of motion, planar wakes, and use Theodorsen's theory to represent the 2D solution in chordwise strips along the span of the wing. However, unlike in 2D where the frequency-domain solution (Theodorsen) has a time-domain equivalent (Wagner's indicial response), 25 it has not yet been possible to derive time-domain equivalents for ULLTs.
Low-order models based on lifting-line theory in the time domain hence employ more assumptions and approximations in comparison with harmonic ULLTs. A common assumption made is to treat the 2D strip-wise solutions as unsteady (using nonlinear inner 2D methods which may apply to non-harmonic kinematics), while treating the trailing vortices behind the wing as steady (in essence, Prandtl's LLT applied to a time-varying wing circulation distribution). This assumption is referred to as the "pseudosteady assumption" in this article.
The first unsteady finite wing analysis in the time domain is due to Jones, 35 although it is only applicable to elliptic planforms. A time-domain ULLT is presented by Boutet and Dimitriatis 7 using an inner solution based on Wagner's method, 67 and a pseudosteady assumption. Ramesh et al. 54 introduce a lifting-line theory for a geometrically non-linear inner solution, again using a pseudosteady assumption. Devinant 19 presents a method by which the pseudosteady assumption can be numerically removed for wake wavelengths of the wing span scale or larger. Bird et al. 5 apply this method to a geometrically non-linear inner solution. Sugar-Gabor et al. 59 suggested a method where the unsteady wake was only accounted for in the outer 3D problem, allowing for roll and yaw kinematics. Lifting-line theory has also been used to model wings in large-scale simulations where resolving the chord-scale would be prohibitively expensive, for example in Caprace et al. 11 Strip theory (no 3D effects) and the pseudosteady assumption (steady 3D effects through Prandtl's LLT) are attractive owing to their mathematical simplicity and easy of implementation. However, the implications of the pseudosteady wake model -or any other simplified wake model -are not well understood. Whilst Guermond and Sellier's analysis 28 found that a pseudosteady wake model was only asymptotically valid in the very low frequency regime, the practical implications of this assumption in aerodynamic flows are less clear.
Additionally, only limited validation of unsteady lifting-line theory has been undertaken. Sclavounos 57 and Boutet and Dimitriatis 7 compared their methods to the unsteady vortex lattice method. This method shares some assumptions with unsteady lifting-line theory, and requires care in discretization to satisfy the Kutta condition. 55 Sclavounos compared results against the whole-wing lift coefficient of a rectangular and an elliptic wing undergoing oscillating kinematics. Boutet and Dimitriatis examined whole-wing lift and moments for indicial and oscillating kinematics. Guermond and Sellier 28 compared their work to a panel method for a pointy wing. Much of experimental and CFD work undertaken on oscillating wings (listed earlier) post-dates publication of unsteady lifting-line theories. These works show that the vortical topology of the wake is far more complex than might be assumed, even at relatively low amplitudes.
Therefore, this paper firstly investigates the accuracy of a ULLT based upon Sclavounos' interaction kernel in comparison to CFD for rectangular wings oscillating in heave and pitch. The use of Euler CFD avoids the explicit assumptions of unsteady vortex lattice methods used in previous research. Whilst lifting-line theory is not strictly valid in the context of rectangular wings, 61 they are of great practical relevance.
Secondly, this paper modifies this 'complete' ULLT in order to investigate the effects of simplified wake models and the accuracy of frequency domain unsteady lifting-line theory. The work is based upon a framework similar to that of Sclavounos because of its arrangement as the solution of an integro-differential equation (akin to Prandtl's LLT), unlike the work of Van Dyke 61 or Guermond and Sellier. 28 This also makes it more amenable to modification of the 3D interaction kernel. Understanding the trade-off between wake-model complexity and ULLT accuracy is useful for those building more complex numerical ULLTs.
The article is therefore laid out as follows: in section 3, ULLT is introduced, based upon the framework of Sclavounos with minor modifications. The 3D interaction kernel and alternative forms of the kernel representing different wake assumptions are discussed in section 3.4. The 'complete' ULLT is validated against CFD in heave and pitch for the Euler regime in section 4.1.1. Next, the impact of using simplified ULLTs is explored for whole wing loads in section 4.1.2, the wake vorticity distribution in section 4.1.3, and wing load distributions in section 4.1.4. A brief discussion of the results obtained using the simplified wakes is then given in section 4.1.5. Section 4.2 compares the ULLTs to experimental data, confirming that ULLTs are relevant even when the flow is not inviscid before the conclusions are finally presented in section 5.

Theoretical approach
Sclavounos' ULLT considers potential flow past a thin, unswept wing subject to a uniform freestream velocity (U) and undergoing small-amplitude, harmonic pitch (α)

The outer solution
The lifting line and its oscillating wake as seen in the outer domain are illustrated in figure 1.
The bound vorticity / circulation of the wing is described by the complex quantity, Kelvin's theorem dictates that the change in bound circulation is convected into the wake. Hence the spanwise (γ w y ) and streamwise (γ w x ) wake vorticity at the wing trailing-edge are Vorticity is convected only by the free stream, hence Now consider a two dimensional slice in the x-z plane, as shown in figure 1b. From this perspective, the lifting-line appears as a point vortex, with a wake on z = 0, It is in this plane that the Kelvin condition for spanwise circulation, given by equation 2, is satisfied. For a bound circulation of unit amplitude, the velocity potential is given by The first term represents the bound circulation at origin and the second term represents the 2D wake vortex sheet.
The streamwise vorticity, γ w x , due to equation 3 must still be included. The effects of this are given as the second term in the asymptotic expansion of the full outer velocity potential equation, derived by Sclavounos as where the kernel K(y) will be discussed in section 3.4. This second term represents an unsteady downwash accounting for the 3D interaction in the outer domain. Differentiating this term with respect to z gives the induced downwash due to finite wing effects. Since this second term is invariant with respect to x and linear in z, the downwash at any spanwise location is uniform in the x-z plane over the section chord. This is a simplification that restricts the asymptotic validity of the theory to span-scale wake wavelengths and longer (very-low and low frequencies) according to Guermond and Sellier. 28 The effect of this assumption in predicting loads on the wing are investigated and discussed in section 4 by comparing lift coefficients from ULLT against numerical simulations over a large frequency range.

The inner solution
In the inner domain, the detail of the flow around a chord section is considered. Since it is assumed that 2s ≫ c, and that changes in the flow happen on the span lengthscale, the problem can be considered as 2D. The problem consists of a symmetric airfoil at a certain spanwise section undergoing small-amplitude pitch and plunge oscillations, where h * 0 is plunge amplitude per unit chord, α 0 is pitch amplitude, ψ is the phase between plunge and pitch, and ω is the frequency of oscillation typically expressed as a chordwise or spanwise reduced frequency Figure 2 shows this problem, for which the solution is given by Theodorsen. 60 The results are expressed here in terms of a general unsteady thin-airfoil theory 36,52 which allows easy modification of the approach for various wing geometries and kinematics. The bound vorticity distribution over the airfoil section is taken as a Fourier series where θ is a transformation variable related to the chordwise coordinate as x i = 0.5c(1 − cos θ ), and A 0 , A 1 ...A N are time-dependent Fourier coefficients that vary over the span. The A 0 term represents the suction peak caused by the flow having turn around the airfoil leading edge, and the Kutta condition (zero vorticity at the trailingedge) is enforced implicitly through the form of the Fourier series. The Fourier coefficients are determined by enforcing the zero-normal-flow boundary condition on the aerofoil chord line. The bound airfoil circulation and the sectional lift and pitching moment coefficient from unsteady thin-airfoil theory are given by where x * m is the reference location per chord (ranging from 0 − 1) about which the pitching moment is calculated.
Ramesh 51 has derived the Fourier coefficients corresponding to Theodorsen's solution as where C(k) is Theodorsen's function, and S(k) is the Sears function, where K 0 (z) and K 1 (z) are modified Bessel functions of the second kind. 46 Q n (k) are wake coefficients, and W 3qc is the normal downwash at the airfoil three-quarter chord location, calculated from the kinematics as Using equations 7 and 17, the 2D normal downwash at three-quarter chord location at any spanwise location from plunge and pitch are The 2D bound circulation, lift coefficient and moment coefficient for pitch and plunge kinematics at any spanwise location are calculated from equations 10, 11, 12, 13, 18 and 19 as where H (2) 0 (z) and H 1 (z) are Hankel functions of the second kind, and x * p is the pivot location.

Matching of the inner and outer solutions
The asymptotic inner expansion of the outer solution velocity potential is given in equation 6. For the inner solution, the velocity potential is constructed as a combination of homogeneous and particular solutions. 57 The homogeneous solution (to model 3D induced effects) is determined from the interaction of the airfoil with a uniform downwash, and the particular solution is obtained from the 2D solutions for pitch and plunge kinematics derived above.
where the first and second terms are the particular and homogeneous solutions, respectively. The uniform downwash iωze iωt corresponds to unit plunge amplitude h 0 = h * 0 c = 1 for which the associated velocity potential is φ 2D hn with resulting bound circulation Γ 2D hn . The outer expansion of the inner velocity potential at large distances from the chord is given as The outer and inner expansions, equations 6 and 27, may now be matched. Equating the φ 2D terms allows an expression for the bound circulation distribution to be extracted as leaving The contributions to circulation from pitch and plunge may be added together since the theory is linear. The lifting-line integro-differential equation for the timevarying spanwise circulation is obtained from equations 28 and 29 as We note that all terms in the lifting-line equation contain the common factor e iωt . An approximate solution for the complex circulation amplitude Γ 0 (y) can be obtained by expressing it in a Fourier series where y = −s cos ζ . For problems where both the kinematics and planform of the wing are symmetric about y = 0, the even m terms can be neglected. For a rectangular wing, the problem is then reduced to numerically solving the equation, at collocation points distributed over the span. The 3D induced downwash is given by The Fourier solutions A 0,1...N (y,t) representing the inner solutions across the span may then be corrected for 3D effects by evaluating equation 13 with the total 3D downwash at the three-quarter chord location The allows the calculation of the stripwise circulation, lift coefficient and pitching moment (accounting for 3D effects) from equations 10, 11 and 12. The wing lift and pitching moment are where the 2D coefficients have been corrected to include the 3D correction as C l/m = C 2D l/m − FC l/m hn .

The kernel K(y)
The kernel K(y), first introduced in equation 6, represents the spanwise interaction of the inner solutions. It must account for the difference between the inner solution, that assumes that the flow is 2D, and the actual 3D nature of the problem via the outer solution. The kernels based on various underlying assumptions considered in this article are expressed below in terms of the non-dimensional spanwise distance y * = y/s. For strip theory, all three dimensional interaction is neglected. The spanwise wake vorticity γ w y is modeled in the inner solution, but the strip theory outer solution neglects both the shed streamwise vorticity γ w x , and the correction for the fact that the variation of γ w y over the span of the wing is not captured by the inner solution. Since all interaction between the inner domains is neglected, the strip theory kernel K 2D is A pseudosteady kernel accounts for the shed streamwise vorticity γ w x in the outer domain, as given by equation 3. However, it neglects the sinusoidal variation with respect to downstream coordinate given in equation 4. Again, the variation in γ w y with respect to span is not corrected for. The resultant pseudosteady kernel K P is equivalent to that of Prandtl: and the ULLT based on this kernel is abbreviated P-ULLT.
If the sinusoidal variation of γ w x with respect to x given by equation 4 is accounted for the streamwise vorticity kernel K x can be obtained. The Biot-Savart law can be applied to the streamwise vorticity field. The downwash on a point of the wing at y 0 due to a section of the wing d y is therefore γ wx : harmonic variation in x-direction C-ULLT K C γ wy : harmonic variation in x-direction (Complete) γ wx : harmonic variation in x-direction allowing K S to be obtained as where I n (x) and K n (x) are the modified Bessel functions of the first and second kind respectively, and L n (x) is the modified Struve function. 46 Once again this neglects variation in γ w y with repect to span in the outer solution. The ULLT based on this kernel is referred to the simplified ULLT (S-ULLT). Sclavounos obtained a kernel K C that accounts for both the shed streamwise vorticity γ w x and the 3D correction to the effects of the shed spanwise vorticity γ w y .
where E 1 (x) is the exponential integral 46 and The ULLT using Sclavounos' full kernel is denoted as the complete ULLT (C-ULLT) in this research. The different ULLT models described above may be summarized in terms of the way the wake is modeled in the outer solution as shown in table 2.
At low frequencies the unsteady solution approaches the pseudosteady solution. Accordingly lim As the oscillation frequency tends to infinity, 3D effects become negligible, and the kernels that include sinusoidal variation of the shed wake approach the strip theory solution: lim

Results
ULLT is derived from potential-flow theory based on incompressible flow with zero viscosity, and employs further simplifying assumptions of high aspect ratio and low reduced frequency. To determine the range of validity of ULLT across the relevant range of aspect ratio and reduced frequency, and to study the influence of the simplified wake models on the solution, the ULLT is first compared against numerical computations of the incompressible Euler equations. In section 4.2, the ULLT models are validated against previously published experimental data in realistic aerodynamic conditions. Comparative remarks about the ULLTs are then made in section 4.1.5.

Analysis with Euler Computational Fluid Dynamics
Numerical computations of the unsteady incompressible Euler equations are performed using the open-source CFD toolbox OpenFOAM. A body-fitted, structured computational mesh is moved according to prescribed plunge and pitch kinematics, and the time-dependent governing equations are solved using a finite volume method. A second-order backward implicit scheme is used to discretize the time derivatives, and second-order limited Gaussian integration schemes are used for the gradient, divergence and Laplacian terms. The pressure implicit with splitting of operators (PISO) algorithm implements pressure-velocity coupling. This in-house setup has previously been used with the incompressible Navier-Stokes governing equations to study leading-edge vortex shedding on finite wings, 5,6 and with the incompressible Euler equations to verify unsteady potential flow solutions for an airfoil. 51 In this section, we use Euler CFD solutions to validate the 3D aerodynamic loads obtained from ULLT, and to study the influence of ULLT wake models on the loads and load distributions. In the CFD setup, inviscid flow is considered with kinematic viscosity set to zero, and a slip boundary condition is employed for the moving wing surface. A NACA0004 section is chosen to best match the theoretical assumptions of thin section and Kutta condition at the trailing edge.
Three aspect ratios (8, 4, and 2) are considered. Cylindrical O-meshes for half the wings are constructed since the pitch and plunge kinematics considered are all symmetric about the wing root. The meshes have 160 cells around the wing section, with increased resolution near the leading and trailing edges. The wall-normal direction has 115 cells with the far-field extending 20 chord lengths in all directions around the section. In the spanwise direction, the aspect ratio 8, 4 and 2 wings have 218, 199 and 87 cells respectively over the wing, with increased resolution near the wingtip. For all three wings, the spanwise domain extends 5 chord lengths beyond the wingtip, with 100 cells in this region. A diagram of the mesh is shown in Fig. 3. Symmetry boundary condition is used for the circular domain at the wing root, and freestream (inlet/outlet) boundary conditions are used at the spanwise and wall-normal far-field domains. The freestream boundary condition behaves as a zero-gradient condition when fluid is flowing out of the boundary face, and as a fixed value condition (equal to freestream) when fluid is not flowing out. Fig. 3 Dimensions of the mesh used for CFD. The symmetry boundary condition is colored gray, and the freestream boundary condition is colored red. Harmonic pitch and heave kinematics with chord reduced frequencies k of 0.0, 0.125, 0.25, 0.5, 1.0 and 1.5 were simulated. Small amplitudes of oscillation, 0.01c for plunge and 1 • for pitch were used, again to best satisfy the theoretical assumptions and ensure attached flows. The CFD cases are shown in table 3.

Validation of lift and moment coefficients from C-ULLT
In this section the C L and C M (about mid-chord) obtained from Sclavounos' complete ULLT are validated against the results obtained using CFD for heave and pitch kinematics. For both heave and pitch, the loads are normalized by the oscillation amplitude (since ULLT is linear with respect to this parameter) and the wing has a zero mean pitch angle. For heave, the loads are further normalized by the chord reduced frequency (k) since as seen in equations 22 and 24, it is a common multiplier for C l h and C m h (unlike for pitch). This is done to better differentiate the various curves at low frequencies as seen below.

Heaving kinematics
The normalized amplitude and phase of the wing lift and moment coefficients for heave oscillations cases are shown in figure 4. The theoretical results and the CFD results are shown by lines and points respectively. Strip theory (2D Theodorsen solution) is also shown by the bold line as a reference; the differences between strip theory and the three ULLT curves show the influence of 3D unsteady induced downwash.
Examining the result from the CFD in figure 4(a), the normalized lift |C L |/kh * 0 is lower for lower aspect ratios across all frequencies. In the frequency range approximately under k = 0.5, |C L |/kh * 0 decreases for AE Rs 8 and 4, and remains roughly constant for AE R2. Above this frequency, the |C L |/kh * 0 curve increases for all AE Rs, with the gradients of the curves being slightly smaller for lower aspect ratios.
The C-ULLT (Sclavounos' complete ULLT) results in figure 4(a) broadly predict the trends found in CFD. In the regime where ULLT assumptions are best satisfied, at high aspect ratio 8 and for low frequencies, the comparison with CFD is excellent. As aspect ratio decreases and as frequency increases, the prediction worsens. At low frequencies, the errors in C-ULLT predictions for AE R4 are relatively small but larger for AE R2. As frequency increases, C-ULLT follows the trends of the curves from CFD, but overpredicts |C L |/kh * 0 . there is also an overestimation of the gradient, suggesting that added mass is overpredicted by C-ULLT. This error increases as aspect ratio becomes smaller. It is likely that this error is caused by the ULLT assumption of 2D flow near the wing tips. This assumption is broken for rectangular and elliptic wings. The phase of wing lift coefficient from CFD and C-ULLT are compared in figure 4(b) . In the limit of low frequency, the C-ULLT predicts that all aspect ratios would have the same C L phase lag of 90 • , but that the response of the phase with increasing k would be different according to aspect ratio. Higher aspect ratios increase the phase lag, whilst lower aspect ratios decrease it. These trends agree well with the CFD results. As frequency increases, the phase prediction obtained for higher aspect ratios remains good but an offset is seen for the AE R2 case.
The comparison of normalized pitching moment |C M |/kh * 0 between C-ULLT and CFD is shown in figure 4(c). Similar to lift, the predictions are best for the high aspect ratio 8 wing. As AE R decreases, the error increases, particularly at high frequencies. In both CFD and C-ULLT, the |C M |/kh * 0 curve approaches a limiting value with increasing k. However, contrary to the CFD results, the C-ULLT predicts that this value is independent of aspect ratio.
For the phase of C M shown in figure 4(d), the CFD shows that at all aspect ratios, the phase lead of C M initially decreases with k. This phase lead then starts increasing again in the region of k = 0.5. C-ULLT predicts the trends of the CFD correctly, particularly in the expected regime of validity (high aspect ratio, k ≈< 1).
Overall, the C-ULLT predicts the results and trends of the CFD reasonably well, albeit worse at low aspect ratio. Importantly however, even at such lower aspect ratios, it always provided a better prediction of the CFD result than could be obtained using strip theory.

Pitching kinematics
Next, the C-ULLT is validated against CFD for leading-edge pitching kinematics. The amplitude and phase of lift and moment coefficients for the pitch oscillations cases are shown in figure 5.
The CFD result for the normalized lift amplitude |C L |/α 0 in figure 5(a) initially decreases with respect to k at aspect ratio 8, stays constant at aspect ratio 4 and increases slightly at aspect ratio 2 (until about k ≈ 0.25). Lower aspect ratio leads to a lower |C L |/α 0 . As k increases further, |C L |/α 0 increases super-linearly. C-ULLT reflects the trends of the CFD, with the initial gradient and y-intercept being dependent upon aspect ratio. C-ULLT obtains a better prediction for all frequencies with increasing aspect ratio, with errors increasing with higher frequencies.
The CFD data in figure 5(b) shows that initially the C L starts in phase to the kinematics. As the frequency increases a phase lead develops. Lower aspect ratios result in a larger phase lead. The phase of C L is excellently predicted by C-ULLT at aspect ratio 8. The accuracy of the prediction worsens at lower aspect ratio, but it remains good, qualitatively reflecting the trends of the CFD.
CFD results for the normalized pitching moment amplitude in figure 5(c) show that |C M |/α 0 is lower for lower aspect ratios. C-ULLT correctly predicts the trends of the CFD. For high aspect ratio, it constantly over-predicts the CFD result by a small amount. At lower aspect ratios, the over-prediction is initially small before growing. In figure 5(d), CFD shows that C M is initially approximately in phase with the kinematics. As frequency increases, the C M leads the kinematics. At low frequencies, the higher aspect ratios have lower phase lead. As frequency increase the phases become more similar, approximately converging at k ≈ 1. As frequency increases further, the higher aspect ratios develop larger phase lead than the lower aspect ratios cases. C-ULLT provides a good prediction of the aspect ratio 8 and 4 arg(C M ) results at approximately k ≈≤ 0.5. As frequency increases, the prediction of C-ULLT worsens. Contrary to the CFD result, it predicts that the phase lead of C M always decreases with aspect ratio.
Overall, it was found that C-ULLT was capable of predicting the CFD results with good accuracy for rectangular wing pitching about the leading edge. The C-ULLT always provided a better prediction of the CFD results than strip theory, and was most accurate for AE Rs 4 and 8 with k ≈< 0.5.

Influence of ULLT kernel on lift coefficient
In section 4.1.1 it was found that the agreement between CFD and Sclavounos' complete ULLT was reasonably good for both lift and moment coefficients, for both pitch- ing and heaving kinematics. The agreement was excellent in the "ideal" regime (for ULLT) of high aspect ratio and low chord reduced frequency. In section 3.4 alternate ULLT models with approximated kernels representing the 3D unsteady induced downwash were introduced. These were the S-ULLT (simplified) and the P-ULLT (pseudosteady), whose assumptions are summarized in table 2. In this section, the implications of approximating the ULLT wake model/kernel are investigated by comparing the normalized lift coefficient amplitude from all three ULLTs against CFD for heaving kinematics. Figure 6 shows the |C L |/kh * 0 curves for the three aspect ratios of wing studied, across a range of chord reduced frequencies.
All of the ULLTs predict the same low frequency |C L |/kh * 0 limit -an expected property due to the interaction kernels of the C-ULLT and the S-ULLT becoming equivalent to the P-ULLT kernel in the low frequency limit. This limiting value reduces with aspect ratio of the wing. For all AE Rs in general, the prediction from C-ULLT is larger in value than that from S-ULLT, which in turn is larger than the prediction from P-ULLT. All methods show an initial negative slope in the normalized C L amplitude until k ≈ 0.5, and a positive slope thereafter.
At high frequencies k ≈> 1, |C L |/kh * 0 has a linear slope in all 3 ULLT. This is owing to the fact that the added mass contribution to |C L |/kh * 0 , which varies linearly with k (see equation 22), is dominant at high frequencies. The curves for the C-ULLT and the S-ULLT converge due to their identical high frequency limiting behavior. With sufficiently high k, they will converge with strip theory. The P-ULLT predicts a slightly lower gradient. The interaction kernel of the P-ULLT does not have a zero high frequency limit. Consequently, even at high frequency there is an erroneous 3D downwash correction, which luckily results in a better agreement with the results from CFD in this regime. This luck may be specific to the rectangular wing planforms studied in this paper. For wing planforms where the assumptions of ULLT are satisfied, this downwash would introduce error. Looking specifically at the AE R8 results presented in figure 6(a), initially the C-ULLT and the S-ULLT give a good prediction of the CFD result and the P-ULLT under-predicts the result. As frequency increases, the C-ULLT and the S-ULLT overpredict |C L |/kh * 0 with the over-prediction increasing with frequency. The P-ULLT provides the best match compared to the CFD result out of the ULLTs at k ≥ 1.
As aspect ratios decrease, as shown in figures 6(b) and 6(c) for AE Rs 4 and 2 respectively, there are two major trends of note.
Firstly, the initial negative slope of the |C L |/kh * 0 curves decreases, and there is a greater difference between the C-ULLT and S-ULLT curves. The difference between the C-ULLT and the S-ULLT curves represent the failure to correct for the changing wake γ w y distribution with respect to span in the outer domain of the S-ULLT. At lower aspect ratios, these corrections accounting for finite-wing effects become more important.
Secondly, the ULLTs appear to increasingly over-predict the |C L |/kh * 0 low frequency limit with decreasing aspect ratio. This effect is most prominent at aspect ratio 2, but explains why at aspect ratio 4 the S-ULLT unexpectedly appears superior to the C-ULLT. Lifting-line theory is based on the assumption of high aspect ratio, so increasing error with decreasing aspect ratio is expected.
Summarizing, all the ULLTs provide a better prediction of the CFD results than the strip theory (2D Theodorsen) result. In regimes where the basic ULLT assumptions of high AE R and low k are satisfied, the complete C-ULLT provides the best predictions. For lower AE Rs at low frequencies (k ≈< 0.5), the simplified S-ULLT provided better predictions than C-ULLT. At high frequencies k ≈> 1, for all aspect ratios in general, the pseudosteady P-ULLT provides the best results. This is because the 3D unsteady induced downwash in this method doesn't tend to zero at high frequencies (unlike in the other two ULLTs), which better reflects the results from Euler CFD.

Comparison of wake topologies in the ULLT kernels
To better illustrate the differences between the wake approximations in the three ULLTs, the vortex sheet strength in the wake from CFD was compared against the vortex sheet strength assumed for the outer domain of the ULLTs. This is shown in figure 7 for the AE R4 wing oscillating in heave with k = 0.5.
To obtain the plotted data from the CFD cases, the vorticity was integrated over a line in the z direction at each point on the x-y plane. For ULLTs, the vorticity as assumed in the outer domain is shown. The C-ULLT predicts the CFD results well. Some error is introduced, as would be expected, near to the wing tip. This error is most visible in the γ w x plots. In the C-ULLT, a singular distribution of γ w x is predicted due to the assumed remapped Fourier distribution of bound circulation. In the CFD, this streamwise wake vorticity remains finite, instead spreading out. In the plots of γ w y , the error near the wing tip in the C-ULLT prediction is most visible when comparing the difference in phase between the vorticity at the center of the span and at the tip. In the C-ULLT result, there is a larger phase difference than is observed in the CFD results.
The S-ULLT produces a similar γ w x field to the C-ULLT. There is a small phase difference between the C-ULLT result and the S-ULLT result. At the wing tip, where the streamwise vorticity has the greatest amplitude, the phase is very similar however. The S-ULLT does not model γ w y in the outer domain. Whilst γ w y is included in the inner domain, there is no correction for its variation over the span. This applies to the P-ULLT as well.
The P-ULLT also doesn't account for the sinusoidal variation of γ w x in the outer wake. However, it roughly matches the CFD γ w x distribution close to the wing (up to about 2 chord lengths from the trailing edge). Since the vorticity closer to the wing has a larger impact on the induced downwash, the P-ULLT obtains reasonable solutions despite this assumption.

Influence of ULLT kernel on spanwise lift distribution
An advantage of ULLT over strip theory is that it accounts for finite wing effects (with bound circulation going to zero at the wingtips) when computing force distributions over the wing. For aeroelastic analysis (for which strip theory is often used), these load distributions are of importance. This section compares the lift distributions obtained using the three ULLTs against that obtained from Euler CFD.
The distribution of lift and changes with respect to oscillation frequency. This occurs because of the growing importance of added mass-effects with increasing oscillation frequency, and also because of the dependence of the 3D interaction kernel K on span reduced frequency ν. For the C-ULLT and the S-ULLT the 3D unsteady induced downwash tends to zero as frequency increases.
To compare unsteady lifting-line theory to CFD, the lift distribution was extracted from the CFD data at 8 equispaced time instants over a single oscillation. A sine wave was then fitted to the data using a least squares method to determine the amplitude.
The lift distributions will be compared for rectangular wings of aspect ratios 8, 4 and 2 for heaving oscillations at three different chord reduced frequencies. First, a low frequency k = 0.125 where C-ULLT provides the best agreement with CFD (so long as aspect ratio is high) is studied. Next, a high frequency k = 1.5 where addedmass effects dominate and where P-ULLT was seen to make the best predictions is studied. Finally, an intermediate frequency k = 0.5 where the loads are influenced strongly by both circulatory and added-mass effects is examined.

Low frequency behavior
The lift distribution amplitudes from ULLTs are compared against Euler CFD and strip theory for the low frequency case (k = 0.125) in figure 8.
The CFD results show a mostly smooth distribution of lift coefficient, with the maximum |C l | being found at the wing root (y = 0). Near the tip, |C l | starts to decrease ever more rapidly, excepting a spike at the very wing tip. This is caused by separation at the sharp edges of the wing tip.
The difference between the |C l | curves predicted by the different ULLTs at k = 0.125 is small due the fact that they all tend to the P-ULLT kernel as frequency decreases. The results obtained by the C-ULLT and S-ULLT are in particular more similar. The C-ULLT predicts the highest |C l | followed by the S-ULLT and then the P-ULLT.
At aspect ratio 8 where ULLT is expected to be most valid, the C-ULLT matches the CFD result well in the center of the wing but has increasing errors as the wing tip is approached. As aspect ratio decreases, the C-ULLT over-predicts the |C l | at the center of the wing. Consequently at aspect ratio 4 either the S-ULLT or the P-ULLT provides the best prediction of the CFD result. At aspect ratio 2, all of the ULLTs over-predict |C l |. The P-ULLT, which always gives the lowest |C l |, is therefore is closest to the CFD result. In all cases, ULLT is superior to strip theory.
High frequency behavior Figure 9 shows a comparison of CFD and ULLT lift distribution predictions for the high frequency case k = 1.5.
In section 4.1.2, it has been noted that the pseudosteady ULLT gives the best results at high frequencies, as the incorrectly modeled downwash inadvertently reduced the lift overestimate due to the lifting-line assumption being poor at rectangular wing tips. This is further confirmed in figure 9; P-ULLT agrees best with CFD for all three aspect ratios while C-ULLT and S-ULLT over-predict the lift distribution. For AE R8, the P-ULLT prediction matches with CFD over most of the wing, until y * ≈ 0.7. As aspect ratio decreases, the distance from the wingtip where errors are present increases, as expected.

Intermediate frequency behavior
In the intermediate frequency range, the choice of best ULLT model is less clear. Figure 10 shows a comparison of CFD and ULLT lift distribution predictions at k = 0.5, where added-mass and circulatory effects are equally important.
The C-ULLT lift distribution curves have the same shape as those from CFD at all aspect ratios. However, there is an offset which increases as aspect ratio decreases. The simplified S-ULLT which in general predicts lower lift values than C-ULLT hence provides better predictions for AE Rs 8 and 4. P-ULLT which always predicts lower lift values than both C-ULLT and S-ULLT, provides the best prediction for AR2.

Choice of ULLT kernel
The comparison of ULLT models against CFD in section 4.1 provides guidance on choosing the most suitable ULLT depending on the problem parameters. In the regime where the assumptions used in ULLT perturbation analysis are best satisfied, i.e. at high aspect ratio and low chord reduced frequency, Sclavounos' C-ULLT works best. On the other hand, at high chord reduced frequencies, for all aspect ratios, the pseudosteady P-ULLT provides the best predictions. P-ULLT also provides better predictions than the other two ULLTs at low aspect ratio, across all frequencies. At higher aspect ratios and at intermediate frequencies, the S-ULLT prediction (which lies between the C-ULLT and P-ULLT predictions) provides the best agreement with CFD. The superiority of the P-ULLT and S-ULLT in comparison to the C-ULLT is likely due to the error from the rectangular wing tips. Strip theory was shown in all the results as a reference to illustrate the importance of 3D effects (unsteady induced downwash) in various regimes; ULLT universally provides better predictions than strip theory.
All the ULLTs are easy to implement and have very low computational cost in comparison with vortex lattice methods and CFD. C-ULLT and S-ULLT have more complex kernels than P-ULLT. For numerical equivalents to these theories, a method similar to the C-ULLT is most complex to implement, with Devinant 19 showing how it is necessary to cancel the component of spanwise wake vorticity present in both inner and outer domain. An method based on the S-ULLT alleviates this complexity, whilst remaining more accurate than a pseudosteady method.

Validation with experimental data
In section 4.1, CFD simulations of the incompressible Euler equations were used to validate unsteady lifting-line theory and to compare the three kernels representing the 3D unsteady induced downwash. Here, ULLT is validated against against experimental data (at finite Reynolds and Mach number) taken from NASA technical report 4632 49 to confirm that the assumption of inviscid flow does not invalidate its use for practical applications.
NASA TR-4632 49 contains data for an aspect ratio 10.1 rectangular wing with a NACA0015 section undergoing pitch oscillation about its quarter chord. The Reynolds number is 1.951 million and the Mach number is 0.288. A case with chord reduced frequency k = 0.133, average angle of attack 3.98 • , and pitch amplitude 4.35 • was selected. This data is presented in figure 90 of the report. In terms of physical quantities in the wind tunnel, the wing had a span of 60.62 inches, chord of 12 inches, oscillation frequency of 14.02Hz, and was subject to a free stream of 100.58 ms −1 .
The studies conducted in sections 4.1.2 and 4.1.4 indicate that for case described above which falls in the high-aspect ratio, low-frequency regime, (i) the C-ULLT provides the best prediction and (ii) the results from all three ULLTs are nearly the same (see figures 6 and 8). Hence, the experimental data is compared against only the C-ULLT prediction below.
The ULLTs assume a zero average pitch angle. To model a non-zero average, the results of the ULLT are summed with the steady results of Prandtl's lifting-line theory. 50 This is consistent because both theories are linear. Variation of sectional lift coefficient (C l ) with angle of attack, at different locations over the wing span, from C-ULLT (with steady solution added) and experiment are compared in figure 11.
The experimental data has both the largest C l amplitude and the highest average C l (compared to itself) nearest to the center of the wing at y * = 0.25 as shown in figure 11(a). There is a small phase difference between the angle of attack and lift coefficient, resulting in an elliptic curve. At y * = 0.475, shown in figure 11(b), the amplitude of C l remains similar, although the less elliptical curves suggests that the loads are now more closely in phase to the kinematics.
Closer to the tip at y * = 0.800, the amplitude of the C l obtained by the experiment is reduced. Figure 11(d) shows the result very to the wing tip at y * = 0.966. Both C l average and amplitude are further reduced in comparison to the inboard result. The curve obtained in noticeably non-smooth. These results agree with the distribution of |C l | obtained using CFD for lower aspect ratios. Namely, that the |C l | amplitude decreases near the wing tips.
Near the center of the wing at y * = 0.25 and y * = 0.475 (figures 11(a) and 11(b)), the match between theory (ULLT) and experiment is very good. Strip theory, shown for comparison, also provides very good agreement with CFD at these spanwise locations. This is in line with expectations for very high aspect ratio wings (10 in this case), where the behavior near the root would be nearly 2D.
Moving towards the tip, the improvement offered by ULLT over strip theory is clear. The ULLT/LLT slightly over-predicts both the mean and amplitude of C l at y * = 0.8, shown in figure 11(c), and also suggests a larger phase difference between the kinematics and loads than experiment. At y * = 0.986, the ULLT/LLT under-predicts C l . However, at both these spanwise locations near the wing tip, ULLT provides a much better prediction than strip theory.

Conclusions
In this paper, three unsteady lifting-line theories (ULLTs) have been obtained from a common framework: a 'complete' ULLT contains corrections for both components of vorticity in the outer domain, a novel simplified version that considers only the streamwise component, and a pseudosteady, Prandtl-like ULLT. These ULLTs were systematically compared both against each other and against Euler CFD results for various aspect ratio rectangular wings at various oscillation frequencies. The ULLTs were capable of capturing trends in wing lift coefficient with respect to both oscillation frequency and aspect ratio, always providing better results than 2D strip theory. They were also able to capture the changes in lift distribution across the wing with respect to frequency, despite being less accurate close to the wing tips. For aeroelastic analysis, this provides advantages over assumed liftdistribution models unable to account for changes in oscillation frequency.
Comparing between the three ULLTs, it was found that in regimes where ULLT assumptions are best satisfied (high aspect ratio and low reduced frequency), the complete ULLT is most accurate. At high frequencies of motion, the Prandtl-like pseudosteady ULLT provided the best prediction, although this may be luck specific to the rectangular wings studied.
By modifying the wake model, the computational cost could be reduced at the expense of reduced accuracy. The novel ULLT considering only the streamwise component of the vorticity in the outer domain was found to provide an excellent reduction in computational cost with minimal loss in accuracy. The pseudosteady ULLT reduced computational cost further still, but provided worse accuracy at low frequencies.
ULLT was also compared against experimental data, confirming that it is not only applicable to Euler/inviscid problems, but also to practical high Reynolds number flows.
The results obtained in this paper are important because they show how ULLT provides a low computational cost model capable of accounting for the interacting 3D aerodynamic effects of aspect ratio and oscillation frequency. The simplified wake models provide not only a means by which accuracy can be traded for reductions in computational cost, but also guidance on the construction of numerical and timedomain ULLTs. In particular, the streamwise vorticity ULLT demonstrates how an entire component of outer solution wake vorticity can be neglected (along with its self-canceling singularities), with only a small cost to the accuracy of the solution.