Frequency-Domain Lifting-Line Aerodynamic Modelling for Wing Aeroelasticity

Abstract

A frequency-domain lifting-line solution algorithm for the prediction of the unsteady aerodynamics of wings is presented. The Biot–Savart law is applied to determine the normalwash generated by the wake vorticity distribution, whereas steady and unsteady airfoil theories (Glauert’s and Theodorsen’s, respectively) are used to evaluate the sectional aerodynamic loads, namely the lift and pitching moment. The wake vorticity released at the trailing edge derives from the bound circulation through the Kutta condition and is convected downstream with the velocity of the undisturbed flow. The local bound circulation is obtained by the application of the Kutta–Joukowski theorem, extended to unsteady flows. Assuming a bending and torsion wing, this paper provides the aerodynamic matrix of the transfer functions, relating the generalised aerodynamic loads to the Lagrangian coordinates of the elastic deformation. Its rational approximation yields a reduced-order state-space aerodynamic model suitable for an aeroelastic stability analysis and control purposes. The numerical investigation examines the influence of both the wake shed/trailed vorticity modelling and different approximations of the Kutta–Joukowski theorem for unsteady flows on the aerodynamic transfer functions given by the developed frequency-domain lifting-line solver. The accuracy of the solver is assessed by comparison with the predictions obtained by a three-dimensional boundary-element-method solver for potential flows. It is shown that, at least for the frequency range considered, regardless of the approximation of the Kutta–Joukowski theorem applied, the formulation based on the Theodorsen theory provides predictions that are in very good agreement with the results from the boundary element method for a slender wing. This agreement worsens as the wing aspect ratio decreases. A lower level of accuracy is obtained by the application of the sectional loads given by the Glauert theory. In this case, the predictions are more sensitive to the approximation used to express the Kutta–Joukowski theorem for unsteady flows.

1. Introduction

One of the main topics of aircraft design has always been the prediction of aerodynamic loads. The first attempts in this field were made at the beginning of the last century by Glauert [1], Prandtl [2], Theodorsen [3], Wagner [4] and several other researchers, who developed theoretical methods to calculate analytically the lift, drag and pitching moment of airfoils and finite wings, both for steady and unsteady flows. Later, technological progress led to the definition of aerodynamic computational solvers based on numerical methods, which permit to study complex bodies in a wide range of flight conditions. Among these methods, the lifting-line theory (LLT) has proven to be a very effective tool for the preliminary design of aircraft aerodynamics that, although less accurate than CFD solvers, provides a reliable aerodynamic loads prediction at a significantly lower computational cost. Initially developed by Prandtl for the study of steady flows [2], it consists of considering the wing represented by a single lumped vortex representing the bound circulation placed at the quarter-chord line, combined with a system of trailed vorticity that takes into account the spanwise variation in the wing circulation. These trailed vortices induce a distributed normalwash over the wing that alters the local angle of attack. The application of the widely known Kutta–Joukowski theorem for airfoils in steady flows relates the circulation to the local wing lift.

Limited to a narrow range of wing geometries, the original Prandtl model was later improved by several researchers. For instance, Weissinger adapted the LLT to swept wings [5], and his method was simplified by Blackwell [6]. Prössdorf and Tordella extended the LLT to generic curved wings [7], inspiring the further work of Wickenheiser and Garcia on morphing wings of arbitrary curvature, chord distribution and airfoil cross-sections [8]. The model proposed by Weissinger was later modified by Owens to include nonlinear aerodynamics [9], and this type of application has become an object of interest for many research activities (see, for example, [10,11,12].

Next, the LLT was modified in order to include the analysis of unsteady flows, with the introduction of the wake shed vorticity effects. Sclavounos [13] developed one of the first unsteady lifting-line theories (ULLTs), based on the separation of the solution in an outer one, confined in the far-field wake and taking into account the three-dimensional effects (a Prandtl-like theory), and an inner one, due to the portion of the wake near to the wing, where the flow can be assumed to be two dimensional (a Theodorsen-like theory). This model was valid only for rectangular and elliptic wings, while Chiocchia et al. [14] extended the method presented by Prössdorf and Tordella [7] to curved wings in oscillatory motion.

One of the most interesting works dealing with the application of a ULLT approach is that presented by Drela [15], where an approximate extension of the Kutta–Joukowski theorem to unsteady flows was introduced. Drela’s model provides an unsteady-flow version of the Kutta–Joukowski theorem derived from the integration of the Bernoulli theorem under the hypotheses of small, low-frequency flow perturbations (a similar extension to unsteady flows of the Kutta–Joukowski theorem is also presented in [16] for application to the Vortex Lattice Method, VLM).

Born for the study of simple wings, the LLT was applied to several different lifting aerodynamics configurations. For example, Dumitrescu and Cardoş [17] developed a ULLT approach suitable for horizontal-axis wind turbines (HAWTs), including the dynamic stall model of Beddoes–Leishman, while Sebastian and Lackner [18] derived a free-wake model for the study of offshore floating wind turbines. Hindman et al. [19], instead, presented an LLT approach for the ice formation on wings, whereas Liorbano et al. [20] modified the Prandtl LLT using a finite element method approach to study nonplanar wings.

The works of Sclavounos and Drela on ULLTs inspired many authors. Among them, Sugar-Gabor developed a method capable of analysing complex wing geometries with the inclusion of winglets, based on an expanded vector form of the Kutta–Joukowski theorem proposed by Drela [21]. Furthermore, Boutet and Dimitriadis [22] developed a ULLT defined in the time domain through the combination of the Prandtl LLT with the Wagner unsteady airfoil theory applied at each wing section, coupled with Drela’s version of the Kutta–Joukowski theorem. Almost in the same period, Izraelevitz et al. [23] presented a similar method consisting of a state-space adaptation of the ULLT, where both the Wagner indicial response function and the circulation are expressed through a small number of additional states. Recently, Bird and Ramesh [24] presented a development of the ULLT introduced by Sclavounos that provides a less complex model, maintaining a good level of accuracy even at high frequencies.

Following the research line of these last works, the aim of this paper is to present frequency-domain LLT-like formulations based on distributed loads given by (steady or unsteady) sectional theories, combined with the normalwash generated by the wake vorticity derived either from the Kutta–Joukowski theorem or its exact extension to linear unsteady aerodynamics [25]. They are applied to evaluate the aerodynamic matrix collecting the transfer functions that relate the Lagrangian coordinates of a bending–torsion wing to the corresponding generalised aerodynamic forces. Furthermore, it is shown that, sampling this aerodynamic matrix in a given range of frequency, it is possible to determine its rational matrix approximation and hence identify a state-space aerodynamic reduced-order model [26] that, coupled with the wing structural dynamics equations, yields a state-space formulation for wing aeroelasticity that may conveniently be applied for a stability analysis and control design purposes.

In the following, the detailed derivation of the proposed frequency-domain LLT aerodynamic modelling, the definition of the corresponding aerodynamic matrix and its rational matrix approximation procedure are described. Then, the numerical investigation examines the influence on the evaluated aerodynamic transfer functions of the different sectional load theories considered, of the wake vorticity modelling and of different approximations of the Kutta–Joukowski theorem for unsteady flows. The results, shown in terms of the magnitude and phase of the components of the aerodynamic matrix, are compared with the predictions determined by a computational tool based on a boundary element method (BEM) for three-dimensional potential flows extensively validated in the past [27,28,29,30,31,32,33,34].

2. Frequency-Domain LLT-like Modelling

The proposed frequency-domain LLT-like solver is defined as follows. The spanwise distribution of loads is predicted either by the Theodorsen sectional theory for unsteady flows or by the Glauert theory for steady flows. In the former case, the wake is represented only by trailed vortices related to the spanwise bound vorticity distribution (three-dimensional effect), similarly to the Prandtl LLT approach [2] (shed vorticity effects are already included in the Theodorsen formulation). Instead, when the sectional loads are predicted by the steady-flow formulation, the wake is represented by trailed and shed vortices (related, respectively, to spanwise distribution and time variation in the bound circulation).

Note that, when the Glauert theory is applied, the proposed formulation can be considered as the extension of the Prandtl LLT to unsteady flows (indeed, the wake is enriched by the inclusion of the shed vorticity, but the sectional loads are still modelled by a single lifting line). Instead, when the Theodorsen theory is considered, the applied approach can be more appropriately considered as a sort of strip theory for which the aerodynamics of each section of the three-dimensional configuration are described through a two-dimensional model corrected by three-dimensional wake effects and is closely related to the time-domain model presented by Boutet and Dimitriadis [22].

The wake normalwash evaluated at the control points located on the $3/4$-chord line is combined with the normalwash related to the wing motion to define the kinematic input to the sectional aerodynamic theory used to determine the load distribution.

2.1. Aerodynamic Loads

Dividing the wing into a prescribed number, $N_B$, of finite sections, discretising the wake into $N_W$ streamwise elements, adopting a frequency-domain representation for the wake normalwash induced by the wake vortices and recalling that the wake vorticity is convected downstream following the material points, it is possible to derive an expression of the following type relating the bound circulation distribution to the wake normalwash at the $N_B$ sectional control points

$${\tilde{\mathbf{w}}}^i=\mathbf{D}\left(k\right)\tilde{\mathbf{\Gamma }}$$

(1)

where ${\tilde{\mathbf{w}}}^i$ denotes the vector collecting the wake-induced normalwash at the $N_B$ control points and $\tilde{\mathbf{\Gamma }}$ collects the $N_B$ values of the discretised spanwise bound circulation distribution. The $[N_B\times N_B]$ matrix D, which is transcendental in the reduced frequency $k=\omega b/V$ (with $\omega ,V$ and b denoting, respectively, the angular frequency, the undisturbed air velocity and the semichord length), is obtained by application of the Biot–Savart law, with elements defined as

$$D_{lm}\left(k\right)=\sum _{n=1}^{N_W}{e}^{-ik{\tau }_n}{r}_{ln}{S}_{nm}$$

(2)

where ${\tau }_n$ denotes the “age” of the n−th wake vortex (time delay), $r_{ln}$ is the influence coefficient between the n−th wake vortex and the l−th collocation point and $S_{nm}$ are the components of a topological matrix selecting the bound circulation associated to the considered wake strip.

Next, let us consider the Theodorsen formulation for the aerodynamic loads of an airfoil [3]. For the lift force, it yields

$$\tilde{L}={\tilde{L}}_c+{\tilde{L}}_{nc}=2\pi \rho bVC\left(k\right){\tilde{v}}_n^{3/4}+ik\pi \rho b^2{\tilde{v}}_n^{1/2}$$

(3)

where $\rho$ is the air density, ${\tilde{L}}_c$ and ${\tilde{L}}_{nc}$ are the circulatory and noncirculatory parts of the lift, respectively, $C\left(k\right)$ is the lift-deficiency function and ${\tilde{v}}_n^{1/2}$ and ${\tilde{v}}_n^{3/4}$ are the normalwash at the mid-chord and $3/4$-chord points, respectively.

A similar expression is provided for the pitching moment in terms of the section kinematics [3]

$$\tilde{M}=-ik\pi \rho b^2{\tilde{v}}_m+b\left(1/2+a\right){\tilde{L}}_c$$

(4)

with a representing the dimensionless distance, with respect to the semichord b, of the wing elastic axis from the mid-chord point, (e.g., if the airfoil shear centre is at the quarter-chord $a=-1/2$, the pitching moment does not depend on the circulatory lift). See Appendix A for details about ${\tilde{v}}_m$.

Expressing the section kinematics in terms of the wing elastic displacements, the circulatory lift at a given j-th wing section (and corresponding contribution to the pitching moment) can be expressed in the following way to take into account the three-dimensional wake vorticity effects

$${\tilde{L}}_{c_j}=2\pi \rho bVC\left(k\right)\left[{\tilde{w}}_j^i+{\tilde{v}}_{n_j}^{3/4}\left({\tilde{w}}_j,{\tilde{\theta }}_j\right)\right]$$

(5)

where ${\tilde{w}}_j$ and ${\tilde{\theta }}_j$ denote, respectively, plunge displacement at the collocation point of the j−th section and pitch rotation of the same section while ${\tilde{w}}_j^i$ is the wake-induced velocity at the same point. Collecting in the vector L_c the values of the circulatory lift at the $N_B$ wing sections, Equation (5) is recast in the following matrix form

$${\tilde{\mathbf{L}}}_c=\mathbf{C}\left(k\right)\left[{\tilde{\mathbf{w}}}^i+{\tilde{\mathbf{v}}}_n^{3/4}\right]=\mathbf{C}\left(k\right)\left[{\tilde{\mathbf{w}}}^i+{\mathbf{W}}_n^1\left(k\right)\tilde{\mathbf{g}}\right]$$

(6)

where ${\mathbf{W}}_n^1$ is the $[N_B\times 2N_B]$ matrix relating the normalwash at every wing section to the corresponding values of sectional plunge and pitch that are collected in the following $[2N_B\times 1]$ vector (see Appendix A for details)

$$\tilde{\mathbf{g}}={\left\{\begin{array}{ccccc}{\tilde{w}}_1& {\tilde{\theta }}_1& \cdots & {\tilde{w}}_{N_B}& {\tilde{\theta }}_{N_B}\end{array}\right\}}^T$$

(7)

In addition, C is a diagonal transcendental $[N_B\times N_B]$ matrix that, for a rectangular wing, can be simply expressed as

$$\mathbf{C}=2\pi \rho bVC\left(k\right)\mathbf{I}$$

(8)

with I denoting the unit matrix (different values of the diagonal terms appear in the case of a non-rectangular wing).

Note that two different vorticity models can be used to define the sectional lift:

• Model 1 uses the complete Theodorsen’s lift-deficiency function (that implicitly takes into account shed vortices effect) and, therefore, the wake-induced velocity is evaluated including trailed vortices only in Equation (1);
• Model 2 uses the quasi-steady approximation $C\left(k\right)=1$ (so that the Theodorsen theory reduces to the Glauert theory for airfoils in steady flows [1]) and, therefore, the wake-induced velocity is evaluated considering both shed and trailed vortices.

The choice of the model affects both matrix C, through the function $C\left(k\right)$, and the matrix D (in Equation (1)), through the contribution of the influence coefficients, $r_{ln}$.

Similarly, noncirculatory lift and pitching moment can be expressed in the following matrix forms

$${\tilde{\mathrm{L}}}_{nc}={\mathrm{W}}_n^2\left(k\right)\tilde{\mathrm{g}}$$

(9)

$$\tilde{\mathrm{M}}={\mathrm{W}}_n^3\left(k\right)\tilde{\mathrm{g}}+\mathrm{B}{\tilde{\mathrm{L}}}_c$$

(10)

where ${\mathrm{W}}_n^2$ and ${\mathrm{W}}_n^3$ are $[N_B\times 2N_B]$ matrices that relate, respectively, the noncirculatory lift and the contribution of the pitching moment independent from the circulatory lift to the sectional bending and torsion degrees of freedom (see Appendix A for details), while B is a diagonal $[N_B\times N_B]$ constant-coefficient matrix providing the effect of the arm of the circulatory lift with respect to the elastic axis position, namely for rectangular wings

$$\mathrm{B}=b\left(1/2+a\right)\mathrm{I}$$

(11)

2.2. Relation between Bound Circulation and Circulatory Lift

As stated in Equation (1), the definition of wake vorticity requires the knowledge of the bound circulation spanwise distribution that, in lifting-line theories, has to be related to the spanwise distribution of the circulatory lift. Here, this is accomplished considering both the widely used Kutta–Joukowski theorem for steady flows and through its extension to unsteady linear aerodynamics recently proposed by some of the authors [25], whose mathematical formulation moves from the work of Theodorsen for the solution of the velocity potential for circulatory flows around thin rectilinear airfoils in harmonic motion [3]. Specifically, the extended Kutta–Joukowski theorem is expressed in terms of a frequency-response function between bound circulation and circulatory lift, $G\left(k\right)$, which is a transcendental function of the reduced frequency through a combination of Bessel functions [25], and reads as

$${\tilde{L}}_c\left(k\right)=\rho VG\left(k\right)\tilde{\mathrm{\Gamma }}\left(k\right)$$

(12)

Therefore, in this work, two models for the definition of the circulatory lift bound circulation relation are considered:

• The Kutta–Joukowski theorem for steady flows with

  $$G\left(k\right)=1$$

  (13)
• The extension of the Kutta–Joukowski theorem to unsteady flows with

  $$G\left(k\right)=ik{e}^{ik}{K}_1\left(ik\right)$$

  (14)

where $K_1$ denotes the first-order modified Bessel function of the second kind.

From Equation (12), it is possible to relate the bound circulation at each wing section to its sectional circulatory lift and obtain the following compact matrix form

$$\tilde{\mathbf{\Gamma }}=\mathbf{H}\left(k\right){\tilde{\mathbf{L}}}_c$$

(15)

where H is a diagonal transcendental $[N_B\times N_B]$ matrix, whose elements depend on the Kutta–Joukowski model and, for a rectangular wing, is given by

$$\mathbf{H}=G^{-1}\left(k\right)/\rho V\mathbf{I}$$

(16)

Equation (6), combined with Equations (1) and (15), yields the following matrix expression relating the distribution of the wing sectional degrees of freedom to the circulatory lift distribution

$${\tilde{\mathbf{L}}}_c={\left[\mathbf{I}-\mathbf{C}\left(k\right)\mathbf{D}\left(k\right)\mathbf{H}\left(k\right)\right]}^{-1}\mathbf{C}\left(k\right){\mathrm{W}}_n^1\left(k\right)\tilde{\mathbf{g}}$$

(17)

Furthermore, combining Equation (17) with Equations (9) and (10), provides the following $[2N_B\times 1]$ vector, collecting all the sectional aerodynamic loads along the wing

$${\tilde{\mathbf{f}}}_s=\left\{\begin{array}{c}\tilde{\mathbf{L}}\\ \tilde{\mathbf{M}}\end{array}\right\}={\mathbf{E}}_g\left(k\right)\tilde{\mathbf{g}}$$

(18)

where the $[2N_B\times 2N_B]$ matrix, relating the sectional aerodynamic loads to the sectional plunge and pitch motions, is defined as

$${\mathbf{E}}_g=\left[\begin{array}{c}{\left[\mathbf{I}-\mathbf{C}\left(k\right)\mathbf{D}\left(k\right)\mathbf{H}\left(k\right)\right]}^{-1}\mathbf{C}\left(k\right){\mathbf{W}}_n^1\left(k\right)+{\mathbf{W}}_n^2\left(k\right)\\ \mathbf{B}{\left[\mathbf{I}-\mathbf{C}\left(k\right)\mathbf{D}\left(k\right)\mathbf{H}\left(k\right)\right]}^{-1}\mathbf{C}\left(k\right){\mathbf{W}}_n^1\left(k\right)+{\mathbf{W}}_n^3\left(k\right)\end{array}\right]$$

(19)

2.3. Aerodynamic Matrix and Reduced-Order Modelling for Aeroelastic Purposes

For aircraft in symmetric steady flight conditions, the aeroelastic response of the wing is typically studied by considering half-wing modelled as an equivalent bending–torsion beam, clamped at the root (midspan section) and free at the tip.

Thus, to define the aerodynamic matrix of the transfer functions relating the generalised aerodynamic loads to the Lagrangian coordinates of the elastic deformation of the wing, the bending and torsion displacements are first expressed as linear combinations of suitable spatial shape functions, and then the sectional aerodynamic loads are projected onto these functions. For the sake of simplicity, and without loss of generality, bending, $w(x,t)$, and torsion, $\theta (x,t)$, are expressed in terms of one single-shape function, $\psi \left(x\right)$ and $\gamma \left(x\right)$, respectively, with x denoting the spanwise coordinate. Therefore, plunge and pitch of the wing sections appearing in Equation (7) are given by

$$\left\{\begin{array}{c}{\tilde{w}}_j=\tilde{w}(x_j,k)=\psi \left(x_j\right){\tilde{q}}_1\left(k\right)\hfill \\ {\tilde{\theta }}_j=\tilde{\theta }(x_j,k)=\gamma \left(x_j\right){\tilde{q}}_2\left(k\right)\hfill \end{array}\right.$$

(20)

where $x_j$ are the spanwise coordinates of the collocation points, and $q_1$ and $q_2$ represent the Lagrangian coordinates of the wing elastic displacement.

Projecting the sectional forces onto the shape functions yields the vector, f, of the generalised forces as

$$\tilde{\mathbf{f}}={\mathbf{X}}_L{\tilde{\mathbf{f}}}_s$$

(21)

where X_L is the following $[2\times 2N_B]$ constant-coefficient matrix

$${\mathbf{X}}_L=\Delta x\left[\begin{array}{ccccccccc}\psi \left({\xi }_1\right)& & \cdots & & \psi \left({\xi }_{N_B}\right)0& & \cdots & & 0\\ 0& & \cdots & & 0\gamma \left({\xi }_1\right)& & \cdots & & \gamma \left({\xi }_{N_B}\right)\end{array}\right]$$

(22)

with ${\xi }_j=x_j/l$ denoting the dimensionless coordinates of the collocation points (l represents the wing half-span), and $\Delta x$ denoting the spanwise length of each wing section. In addition, the application of the plunge and pitch distributions in Equation (20) provides the following expression for the vector of the sectional displacements in Equation (7)

$$\tilde{\mathbf{g}}={\mathbf{X}}_R\tilde{\mathbf{q}}$$

(23)

where $\tilde{\mathbf{q}}$ is the vector of the Lagrangian coordinates, and ${\mathbf{X}}_R$ is the following $[2N_B\times 2]$ constant-coefficient matrix

$$\begin{array}{c}{\mathbf{X}}_R=\left[\begin{array}{ccc}\psi \left({\xi }_1\right)& & 0\\ 0& & \gamma \left({\xi }_1\right)\\ ⋮& & ⋮\\ \psi \left({\xi }_{N_B}\right)& & 0\\ 0& & \gamma \left({\xi }_{N_B}\right)\end{array}\right]\end{array}$$

(24)

Then, combining Equations (18), (21) and (23) gives the $[2\times 2]$ aerodynamic matrix

$$\mathbf{E}\left(k\right)={\mathbf{X}}_L{\tilde{\mathbf{E}}}_g\left(k\right){\mathbf{X}}_R$$

(25)

such that

$$\tilde{\mathbf{f}}=\mathbf{E}\left(k\right)\tilde{\mathbf{q}}$$

(26)

The matrix E is a frequency-domain transcendental aerodynamic operator. The presence of the transcendental contributions is due to three sources: $\left(i\right)$ the Biot–Savart law used to evaluate the normalwash generated by the wake vorticity (see Equation (1)), $\left(ii\right)$ the Theodorsen theory used to evaluate the sectional aerodynamic loads (see Equation (8)) and $\left(iii\right)$ the application of the unsteady Kutta–Joukowski condition (see Equation (15)). This implies that, in this form, the aerodynamic matrix is unsuitable for aeroelastic stability analysis and control applications.

In order to overcome this difficulty, it is convenient that the following rational matrix approximation form of the aerodynamic matrix along the imaginary axis of the complex plane is

$$\mathbf{E}\left(k\right)\cong {\left(ik\right)}^2{\mathbf{A}}_2+ik{\mathbf{A}}_1+{\mathbf{A}}_0+\mathbf{Q}{(ik\mathbf{I}-\mathbf{A})}^{-1}\mathbf{R}$$

(27)

where A₂, A₁ and A₀ are $[2\times 2]$ matrices, A is a diagonal $[N\times N]$ matrix, where N is the number of aerodynamic poles that can be arbitrarily increased to better approximate the aerodynamic transfer functions, and Q is a $[2\times N]$ matrix and R is a $[N\times 2]$ matrix. All these are constant-coefficient matrices which can be determined by a least-square approximation technique [26].

Indeed, the combination of Equation (27) with Equation (26) followed by the transformation into time domain yields, for $\tau =tV/b$,

$$\left\{\begin{array}{c}\mathbf{f}\left(\tau \right)={\mathbf{A}}_2\ddot{\mathbf{q}}+{\mathbf{A}}_1\dot{\mathbf{q}}+{\mathbf{A}}_0\mathbf{q}+\mathbf{Q}\mathbf{r}\hfill \\ \dot{\mathbf{r}}=\mathbf{A}\mathbf{r}+\mathbf{R}\mathbf{q}\hfill \end{array}\right.$$

(28)

which represents the state-space model of the wing generalised forces derived from the frequency-domain LLT solver presented in this work, perfectly suitable for aeroelastic stability analysis and control purposes. In Equation (28), $\mathbf{r}$ denotes the vector of N additional aerodynamic states introduced by the rational contribution in Equation (27).

3. Numerical Investigation

The main purpose of this section is twofold: $\left(i\right)$ to examine the effects of both the sectional load modelling and steady–unsteady versions of the Kutta–Joukowski theorem on the evaluation of the aerodynamic matrix through the presented LLT-like solver (see Section 2.1 and Section 2.2) and $\left(ii\right)$ to assess the accuracy of this aerodynamic matrix by comparison with the predictions provided by a BEM computational tool for three-dimensional lifting bodies in arbitrary motion. Note that this computational tool was extensively validated in the past against both the experimental and numerical data available in the literature, concerning both fixed and rotary wing configurations [27,28,29,30,31,32,33,34]. Finally, in order to obtain from the developed LLT-like aerodynamic operator a reduced-order model of particular interest for aeroelastic applications, the rational matrix approximation of the aerodynamic matrix is carried out and validated.

The numerical results present the evaluation of the magnitude and phase of the aerodynamic matrix components within the reduced frequency range $k\in [0,2]$. They concern an untwisted rectangular wing with a chord length $c=1$ m and three aspect ratios corresponding to half-span lengths $l=5$ m, $l=2.5$ m and $l=1.25$ m. The undisturbed-flow dynamic pressure is assumed to be equal to unity, whereas the pitching and torsion motions are assumed to occur around the local aerodynamic centre line.

Both the results presented obtained by the LLT-like tool and those provided by the BEM solver are numerically converged ones, in terms of the body and wake discretisation elements (with the latter suitably related to the specific k value of the analysis).

Two different cases of increasing complexity are considered for the definition of the bending and torsion shape functions, namely $\psi$ and $\gamma$:

• Spanwise constant shape functions

  $$\left\{\begin{array}{c}\psi \left({\xi }_j\right)=1\hfill \\ \gamma \left({\xi }_j\right)=1\hfill \end{array}\right.$$

  (29)

  for which the generalised forces correspond to the total lift and pitching moment around the aerodynamic centre line;
• Shape functions expressed as the following approximations of the first bending and torsion natural modes of vibration of a clamped-free beam (named the NMV shape functions, in the following)

  $$\left\{\begin{array}{c}\psi \left({\xi }_j\right)={\xi }_j^2\hfill \\ \gamma \left({\xi }_j\right)={\xi }_j\hfill \end{array}\right.$$

  (30)

3.1. Validation of the Proposed LLT-like Models

For the constant shape functions and $l=5$ m, Figure 1 and Figure 2 present the transfer functions $E_{11}$ and $E_{12}$ (namely the lift due to the heave motion and the lift due to the pitching motion), with S-T and U-T representing, respectively, the results obtained by the sectional loads described by the Theodorsen theory combined with the steady and unsteady versions of the Kutta–Joukowski theorem, and with S-G and U-G representing, respectively, the results obtained by the sectional loads described by the Glauert theory combined with the steady and unsteady versions of the Kutta–Joukowski theorem.

For both transfer functions, the S-T and U-T predictions are practically coincident and in very good agreement with the three-dimensional BEM solution. Instead, the differences between the S-G and U-G results arise both in terms of the magnitude and phase, with the S-G predictions better predicting the transfer function magnitudes and U-G better predicting the phases, as related to the BEM outcomes. These results can be explained by the small three-dimensional effect contribution to the wake normalwash by the trailed vortices, for the considered value of the wing aspect ratio. Indeed, for the S-T and U-T solutions, the different versions of the Kutta–Joukowski theorem affect only the intensity of the trailed vortices and are practically coincident and in agreement with the three-dimensional BEM results. On the other hand, the S-G and U-G results are more evidently affected by the versions of the Kutta–Joukowski theorem applied because these contribute to the definition of the shed vortices that, in these cases, are present in the wake to take into account the wake effects related to the flow unsteadiness. These results demonstrate that, for this wing configuration, the analytical shed vorticity contribution included in the Theodorsen formulation is capable of capturing unsteady-flow effects much better than the simulations combining the steady-flow load description with the numerically evaluated shed vorticity effects.

Next, Figure 3 shows the magnitude of the transfer functions $E_{21}$ and $E_{22}$ relating the pitching moment to the heave and pitching motion, respectively. In these cases, all the solvers considered provide almost identical results (note that the coincident null-phase predictions are omitted for the sake of conciseness). This is an expected result because the pitching motion is assumed to occur around the $1/4$-chord point, thus deleting the contribution of the circulatory lift (see Equation (4)) which, in turn, is the term affected by the application of the Kutta–Joukowski theorem (through the matrix H; see Equation (15)), and by the wake vorticity (through matrices C and D; see Equations (8) and (1)). The predictions from the LLT-like solvers and from the BEM tool are in very good agreement, with some discrepancies appearing at high values of the reduced frequency.

Then, Figure 4 and Figure 5 show the transfer functions related to the NMV shape functions. In this case, the S-T and U-T solutions are not coincident and present small discrepancies, particularly in terms of the magnitude, with respect to the BEM outcomes. These differences with respect to the results in Figure 1 and Figure 2 are due to the fact that the NMV shape functions are such to increase the three-dimensional wing-tip effects, thus increasing, at the same time, the influence of both the trailed vorticity and the version of the Kutta–Joukowski theorem applied.

Nevertheless, the S-T and U-T results are still very accurate. On the other hand, the S-G and U-G solutions present a behaviour similar to that shown in Figure 1 and Figure 2.

Finally, concerning the prediction of the transfer functions $E_{21}$ and $E_{22}$, Figure 6 proves that, as expected, the adoption of the NMV shape functions does not alter the full agreement between all the proposed solutions (the corresponding results perfectly overlap).

3.2. Wing Aspect Ratio Effect

Now, the capability of the proposed LLT-like solvers is assessed for lower aspect ratio wings, still considering the NMV shape functions for the definition of the aerodynamic matrix. For $l=2.5\mathrm{m}$, the LLT-like predictions are compared with the BEM solutions in Figure 7, Figure 8 and Figure 9, while Figure 10, Figure 11 and Figure 12 present the results concerning the wing with a half-span of $l=1.25\mathrm{m}$.

The correlations among the different LLT-like solutions confirm the same features shown for the higher aspect ratio wing (although smaller differences appear between the S-G and U-G results), with those based on the sectional loads provided by the Theodorsen theory in better agreement with the BEM outcomes. However, these figures clearly highlight that the accuracy of the LLT-like simulations tends to decrease as the wing aspect ratio decreases and, hence, the three-dimensional effects correspondingly increase. This is particularly evident in terms of the magnitude of the aerodynamic transfer functions. In addition, it is interesting to note that a significant difference among the LLT-like solvers predictions and BEM predictions is also observed for the $E_{21}$ and $E_{22}$ transfer functions, as the wing aspect ratio decreases (see Figure 9 and Figure 12).

3.3. Reduced-Order Modelling

As stated in Section 2.3, from the rational matrix approximation of the aerodynamic transfer function, it is possible to derive a reduced-order model of the aerodynamic operator to be suitably coupled with the structural dynamics one in order to define the aeroelastic operator.

As an example of this outcome of the application of the proposed frequency-domain LLT-like solvers, the rational matrix approximation is applied to the matrix related to the NMV shape functions, obtained by the U-T algorithm, for the wing with $l=5\mathrm{m}$.

The comparisons between the transfer functions directly given by the LLT-like solver, and those analytically provided by the rational matrix expression of Section 2.3 with the introduction of four poles, are shown in Figure 13, Figure 14 and Figure 15 in terms of the magnitude and phase.

These results demonstrate the excellent accuracy of the approximation procedure, and hence of the reduced-order unsteady aerodynamic model that can correspondingly be derived for aeroelastic applications.

4. Conclusions

A novel frequency-domain LLT-like approach for wing unsteady aerodynamics has been presented. It includes two different versions of the Kutta–Joukowski theorem for the definition of the relationship between the bound circulation and circulatory lift and two different approaches for the evaluation of the sectional loads (namely the Theodorsen theory and the Glauert theory). The magnitude and phase of the aerodynamic transfer functions provided by the proposed solvers have been calculated for two sets of shape functions describing the elastic displacements and three wings with different aspect ratios. The results from the LLT-like tools have been compared to those from a three-dimensional BEM solver for potential flows, extensively validated in the past. The predictions provided by the U-T and S-T solvers (based on the sectional Theodorsen theory), for which only trailed vortices affect the normalwash distribution along the wingspan, are very close and in very good agreement with the three-dimensional BEM solution for a slender wing. Their accuracy decreases as the aspect ratio of the wing decreases, and significant discrepancies with respect to the BEM simulation appear for a very small aspect ratio. The solvers based on the Glauert theory and wake normalwash given by both the shed and trailed vortices (namely the U-G and S-G solvers) provide similar, but different predictions for any wing aspect ratio, particularly in terms of the transfer function magnitude. Their accuracy is satisfactory for high values of the aspect ratio and, akin to the U-T and S-T tools, decreases significantly for a very low aspect ratio. For these solvers, the presence of the shed vorticity in the evaluation of the wake effects on the normalwash produces evident differences between the results based on the steady and the unsteady versions of the Kutta–Joukowski theorem. In addition, it is interesting to observe that the presence of the shed vortices to be considered for the evaluation of the wake-induced velocity makes the U-G and S-G solvers computationally more expensive than the U-T and S-T ones. Nonetheless, all the LLT-like computational tools examined are significantly faster than the higher-fidelity BEM solver applied for comparison, and hence they are particularly suitable for the aerodynamic matrix frequency sampling aimed at aeroelastic investigations for preliminary design and aeroservoelastic purposes. Finally, for the slender wing solved by the U-T approach, the aerodynamic matrix has been successfully approximated through a rational matrix form, from which the aerodynamic reduced-order model to be effectively coupled with the structural dynamics one can be directly derived.