NOMENCLATURE

${A_i}$

area of a wing stirp

$b$

wingspan

${c_i}$

local wing chord

${\rm{C}}(u)$

parameterised form of an aerofoil curve

${C_{{d_i}}}$

drag coefficient of two-dimensional aerofoil located at wing section $i$

${C_D}$

wing drag coefficient

${C_{{D_0}}}$

wing profile (form) drag

${C_{{l_i}}}$

lift coefficient of two-dimensional aerofoil located at wing section $i$

${C_L}$

wing lift coefficient

$CP$

pressure coefficient

$CP_i^{visc}$

pressure coefficient for wing strip $i$ as determined from viscous analysis

$\Delta CP$

pressure coefficient difference

${\textbf{c}}{{\textbf{s}}_i}$

wing strip unit chordwise vector

${\textbf{d}}{{\textbf{l}}_i}$

geometric segment coinciding with a bound vortex filament

${\textbf{d}}{{\textbf{F}}_i}$

aerodynamic force acting on a surface panel or bound vortex segment

${\textbf{d}}{{\textbf{M}}_i}$

two-dimensional pitching moment of the aerofoil located at wing section $i$

${\textbf{F}}$

total aerodynamic force

$\textbf{\textit{G}}\!\left({\textbf{\textit{w}},\boldsymbol\alpha } \right)$

equality and inequality constraints

$J\!\left( {w,\alpha } \right)$

objective functional

$k$

number of NURBS control points

${\textbf{M}}$

total aerodynamic moment

${\textbf{n}_i}$

unit normal vector

${\textbf{n}}{{\textbf{s}}_i}$

wing strip unit normal vector

$N$

number of horseshoe or ring vortices on a lifting surface

${N_{i,n}}$

NURBS basis functions

${{\textbf{P}}_i}$

NURBS control point coordinates

${Q_\infty }$

freestream dynamic pressure

${{\textbf{r}}_i}$

position vector

${\textbf{R}}$

residual vector

$\textbf{\textit{R}}\left( {\textbf{\textit{w}},\boldsymbol\alpha } \right)$

nonlinear system of equations

$Re$

Reynolds number

$S$

wing area

${t_i}$

NURBS parameterisation knots

${\rm{\Delta }}t$

time step

$u$

aerofoil curve parameter

${{\textbf{v}}_{ij}}$

velocity induced by vortex $j$ at control point $i$

${{\textbf{v}}_{rel}}_i$

local relative velocity of a point on the wing

${{\textbf{V}}_\infty }$

freestream velocity

${{\textbf{V}}_i}$

local velocity vector at control point $i$

${\textbf{V}}_i^T$

surface transpiration velocity

${\textbf{VS}}_i^T$

average transpiration velocity for wing strip $i$

${w_i}$

weights associated with NURBS control points

$\textbf{\textit{w}}$

vector of system-dependent variables

$\textbf{\textit{X}}$

coordinates of wake point

$\Delta{y_i}$

width of a wing strip

${\alpha _i}$

effective angle-of-attack

$\boldsymbol\alpha $

vector of system design parameters

${{\boldsymbol\gamma }}$

geometric vector coinciding with one vortex ring side

${\Gamma _j}$

strength of vortex $j$

$\Delta{\Gamma _j}$

correction to strength of vortex $j$

$\varepsilon $

error

$\rho $

air density

${\boldsymbol\psi _{\textbf{1}}}$,${\boldsymbol\psi _{\textbf{2}}}$

Lagrange multipliers

${\rm{\Delta \tau }}$

fictitious time step

${{\boldsymbol\Omega }}$

angular velocity of wing-fixed reference system

2.0 NON-LINEAR LIFTING LINE METHOD

The LLT has seen extensive usage for the analysis and design of straight lifting surfaces with moderate to high aspect ratio, with applications in analysis domains including low-speed aircraft wings, boat sails, propellers and wind turbine blades. Owing to its relatively good accuracy in the range of linear aerodynamic behaviour and its minimal computational costs, various authors have developed alternative formulations to the classical LLT, thus increasing its range of applicability and/or the accuracy of the predicted results. All these formulations, however, retain the fundamental idea of a concentrated distribution of vorticity bound to the lifting surface quarter-chord line. Applications of modified, steady-state and unsteady LLT have included lifting surfaces with arbitrary camber, sweep and dihedral angle^{(Reference Spall, Phillips and Pincock2)}, complex multi-element wings in take-off and landing configurations^{(Reference Gallay and Laurendeau3,Reference Gallay and Laurendeau4)} , analysis of flapping bird wings in forward flight^{(Reference Phlips, East and Pratt5)} and the design and optimisation of wind turbine blades^{(Reference Cline and Crawford6–Reference Fluck8)}.

Various non-linear LLT models have been proposed in literature. In ref. [Reference Chreim, Esteves, Pimenta, Assi, Dantas and Kogishi44], a model was developed to design and analyse the performance of ship propellers. The helical wake created by a constant-pitch propeller was modelled by splitting the typical horseshoe vortices into segments whose position is updated as the wake evolves. The Pistolesi boundary condition was enforced to guarantee the flow tangency condition at a control point whose coordinates were changed according to local sectional aerofoil data read from lookup tables. The model developed in ref. [Reference Gallay and Laurendeau3] used the effective angle-of-attack updated strategy of van Dam and a lifting-line model based on the works of Weissinger, but the non-linear viscous correction was linearised and embedded into a single coupled set of linear equations, thus solving simultaneously for both the circulation and the effective angle-of-attack correction. The results showed good predictions for the stall and post-stall flow regime but were restricted to wings with an elliptical planform. A morphing wing design was introduced in ref. [Reference Gamble, Pankonien and Inman46] with the aim of increased stall recovery, its performance being numerically predicted with a non-liner LLT model. It used two-dimensional Reynolds-Averaged Navier–Stokes (RANS) results to predict the viscous aerodynamic characteristics of the aerofoil sections and a linearisation of Prandtl’s original integrodifferential equation into which the RANS result for the lift coefficient was directly introduced. The method used an artificial viscosity approach to stabilise the iterative solution. Good agreement was obtained with low-speed wind-tunnel results, while the morphing technique proved effective up to maximum lift but not in the post-stall regime. Other non-linear lifting line models can be found in refs. [Reference Owens52–Reference Vernengo, Bonfiglio and Brizzolara54]. These are based on modifying Weissinger’s original method by including viscous aerofoil results^{(Reference Owens52)}, determining a circulation strength correction based on lookup tables of viscous aerofoil results and iterating until the difference between inviscid and viscous circulation values becomes negligible^{(Reference van Garrel53)}, or by updating the coordinates of the control point according to the viscous lift curve slope followed by an application of a Pistolesi-type boundary condition to build a linear system of equations for the circulation values^{(Reference Vernengo, Bonfiglio and Brizzolara54)}.

A non-linear LLT model was developed and presented in refs. [Reference Sugar Gabor, Koreanschi and Botez9,Reference Sugar Gabor, Koreanschi and Botez10]. The model uses the fully three-dimensional vortex lifting law as initially proposed in ref. [Reference Phillips and Snyder1] but also reformulates the equations to allow coupling with two-dimensional experimental aerofoil results provided via a pre-built database file.

In this model, vorticity is distributed in a finite number $N$ of horseshoe vortices, each having its own strength ${\Gamma _i}$. The three-dimensional vortex lifting law^{(Reference Saffman11)} is applied to express the inviscid force ${\bf{d}}{{\bf{F}}_i}$ acting on the bound segment ${\bf{d}}{{\bf{l}}_i}$ of each horseshoe vortex:

(1) \begin{equation}{\bf{d}}{{\bf{F}}_i} = \rho {\Gamma _i}\left( {{{\bf{V}}_\infty } + \mathop \sum \nolimits_{j=1}^N {\Gamma _j}{{\bf{v}}_{ij}}} \right) \times {\bf{d}}{{\bf{l}}_i}\end{equation}

In Equation (1), $\rho $ is the fluid density, ${{\bf{V}}_\infty }$ is the freestream velocity, and ${{\bf{v}}_{ij}}$ is the velocity induced by horseshoe vortex $j$ at a suitably chosen control point associated with horseshoe vortex $i$. Various control point locations can be considered, including on the bound vortex itself or the three-quarter-chord point of the wing’s local chord.

The magnitude of the force acting on a wing strip of area ${A_i}$ and having a local airfoil lift coefficient ${C_{{l_i}}}$ is given by

(2) \begin{equation}\|{{\bf{F}}_i}\| = \frac{1}{2}\rho V_\infty ^2{A_i}{C_{{l_i}}}\end{equation}

Pre-built databases of experimental aerofoil results (or, equivalently, numerical results provided by either RANS-based CFD or panel-boundary layer codes such as XFOIL) are used to provide the $C_{l{i}}$ values. These are pre-calculated at a given number of flow conditions (Reynolds number, angle-of-attack), while for any other flow condition not in the database, a simple linear interpolation is used between the two closest data points available. If the wing strips are taken such that each bound horseshoe vortex segment corresponds to one strip, then the modulus of the force given by Equation (1) can be set equal to that given by Equation (2) and the following non-linear equation is obtained for each of the $N$ horseshoe vortices:

(3) \begin{equation}\left\|\rho {\Gamma _i}\left( {{{\bf{V}}_\infty } + \mathop \sum \nolimits_{j=1}^N {\Gamma _j}{{\bf{v}}_{ij}}} \right) \times {\bf{d}}{{\bf{l}}_i}\right\| - \frac{1}{2}\rho V_\infty ^2{A_i}{C_{{l_i}}}=0,\quad i=1,2, \ldots ,N\end{equation}

The model can be solved for the strengths ${\Gamma _i}$, after which the total aerodynamic force and moment are determined as follows:

(4) \begin{equation}{\bf{F}} = \rho \mathop \sum \nolimits_{i=1}^N \left[ {\left( {{{\bf{V}}_\infty } + \mathop \sum \nolimits_{j=1}^N {\Gamma _j}{{\bf{v}}_{ij}}} \right){\Gamma _i} \times {\bf{d}}{{\bf{l}}_i}} \right]\end{equation}

(5) \begin{equation}{\bf{M}} = \rho \mathop \sum \nolimits_{i=1}^N {{\bf{r}}_i} \times \left[ {\left( {{{\bf{V}}_\infty } + \mathop \sum \nolimits_{j=1}^N {\Gamma _j}{{\bf{v}}_{ij}}} \right){\Gamma _i} \times {\bf{d}}{{\bf{l}}_i}} \right] + {\bf{d}}{{\bf{M}}_i}\end{equation}

where ${\bf{d}}{{\bf{M}}_i}$ is the local, two-dimensional pitching moment of the aerofoil section (again obtained from the pre-built database) and ${{\bf{r}}_i}$ is a vector from the position of the chosen moment reference point to the control point.

The non-linear system of Equations (3) is solved using Newton’s method and requires no under-relaxation to achieve convergence. Typical non-linear LLT models found in literature (see for example refs. [Reference Chreim, Esteves, Pimenta, Assi, Dantas and Kogishi44,Reference Gamble, Pankonien and Inman46,Reference van Garrel53]) require significant under-relaxation if the variable being updated from iteration to iteration is the circulation or the local lift coefficient value. The much better convergence properties of the current model are attributed to the strong coupling between the circulation and the term representing the viscous correction. Other non-linear LLT models first solve for linear circulation values and then iterate through a non-linear correction step until a certain convergence criterion (typically linked to the effective angle-of-attack or the difference between inviscid and viscous sectional lift coefficients) is satisfied. This type of approach is considered to be loosely coupled, and significant under-relaxation is required to guarantee convergence. The model proposed here integrates the viscous data into the non-linear equations and solves directly for the corrected (final) circulation values. It has been shown by other authors^{(Reference Gallay and Laurendeau3)} that such a strong coupling resulted in an increased computation effort per iteration but achieved significantly faster convergence with minimal under-relaxation for a non-linear LLT utilising an effective angle-of-attack correction approach.

It must be noted that the method does not use Trefftz plane analysis for the induced drag but calculates the overall aerodynamic force vector using on-body velocities only. The overall accuracy in estimating induced drag with the near-field approach remains good. The profile drag can also be approximated by

(6) \begin{equation}{C_{{D_0}}} = \frac{1}{S}\mathop \smallint \nolimits_{ - b/2}^{b/2} {C_d}\left( y \right)c\left( y \right)dy \cong \frac{1}{S}\mathop \sum \nolimits_{i=1}^N {C_{{d_i}}}{c_i}\Delta{y_i}\end{equation}

where $S$ is the wing area, ${C_{{d_i}}}$ is the two-dimensional aerofoil drag coefficient, ${c_i}$ is the local chord and $\Delta{y_i}$ is the width of the wing strip (which coincides with the distance between the two trailing segments of a horseshoe vortex).

A first verification and validation test is done using geometrical and experimental data from the NACA 1270 Technical Note^{(Reference Neely, Bollech, Westrick and Graham12)}. The wing is a straight wing having an aspect ratio of 12, a taper ratio of 0.285 and a twist of $3^\circ$. The aerofoil section progressively changes from a NACA 4422 at the root section to a NACA 4412 at the wing tip. The experimental results were obtained for an airspeed of 65m/s and a Reynolds number equal to $4 \times {10^6}$, as calculated with the mean aerodynamic chord value. The numerical results are obtained with 35 strips per semi-span, while the sectional aerofoil database is generated using the two-dimensional XFOIL solver^{(Reference Drela13)}. It is important to point out that, while the strongly coupled integral boundary layer method implemented in XFOIL allows for limited regions of flow separation, its applicability to high-${C_L}$ conditions and the determination of ${C_{{L_{max}}}}$ are not always dependable, leading to a loss of accuracy. However, for the current set of calculations, it was deemed sufficient and chosen over RANS-based CFD calculations due to the considerable time savings obtained while generating the database.

Figure 1 shows a comparison between the numerical and experimental results for the lift, drag and pitching moment coefficients. The estimation of the lift coefficient for angles of attack up to $12.5{^\circ}$ is very accurate, as expected for a straight, high aspect ratio wing. The stalling angle and ${C_{{L_{max}}}}$ are overestimated, although this is due, at least in part, to the overestimation of the same parameters by XFOIL for the NACA 44-series aerofoils. The drag polar prediction is good overall, with underestimations at both very low and very high ${C_L}$ values. Since no data are provided in ref. [Reference Neely, Bollech, Westrick and Graham12] on the turbulence intensity levels, the drag underestimation at low lift might be due to an extended laminar boundary layer in the XFOIL analysis compared with the actual wind-tunnel conditions.

Figure 1. Comparison of lift, drag and pitching moment coefficient values obtained with the non-linear LLT method with experimental data published for the NACA TN1270 high-aspect-ratio straight wing geometry.

Accurate prediction of pitching moment coefficient values is known to be generally very difficult. The model captures the variation trend across the lift coefficient range, but the prediction is offset by a relatively constant value.

A second verification and validation test is done using geometrical and experimental data from the NACA L50F16 Research Memorandum^{(Reference Cahill and Gottlieb14)}. The wing has a moderate sweep angle of $30^\circ $, an aspect ratio of 4 and a taper ratio of 0.6. The wing has a constant NACA 65A006 aerofoil section from root to tip. The experimental results were obtained for an airspeed of 80m/s and a Reynolds number equal to $3 \times {10^6}$, as calculated with the mean aerodynamic chord value. The numerical results are obtained with 50 strips per semi-span, while the sectional aerofoil database is constructed using the experimental results provided in ref. [Reference Loftin15].

A comparison between the numerical and experimental results for the lift, drag and pitching moment coefficients is presented in Fig. 2. An excellent agreement exists for both the lift and drag coefficient results for the ${C_L}$ range below 0.60. This shows that the non-linear LLT model can be used for predicting the aerodynamic behaviour of moderately swept wings. Pitching moment values are accurately predicted only for the low ${C_L}$ range. The lift coefficient plateau occurring around ${C_{{L_{max}}}}$ and the high-lift pitching moment behaviour observed in the experimental results could result from boundary-layer separation in the region close to the tip of the swept wing, with the possible formation of localised quasi-steady leading-edge vortices. This would delay the local loss of lift, accompanied by a significant variation in drag and pitching moment. This highly nonlinear phenomenon cannot be captured by potential flow models such as the lifting line, but the prediction quality for low to moderate angle-of-attack values is very good.

Figure 2. Comparison of lift, drag and pitching moment coefficient values obtained with the non-linear LLT method with experimental data published for the NACA L50F16 moderate-aspect-ratio moderate-sweep wing geometry.

The lack of true 3D interactions is one of the main drawbacks of lifting-line methods. This can become a major source of error, especially for wings having some degree of sweep while close to the maximum lift conditions (as observed in both Figs 1 and 2).

While it is difficult to address this drawback in the simple mathematical framework of the lifting-line method, a promising alternative was proposed recently in refs. [Reference Goitia and Llamas47,Reference Parenteau, Sermeus and Laurendeau55]. In this so-called 2.5D approach, the usual 2D viscous lift curves are corrected either by analytically applying sweep theory to the results (as in ref. [Reference Goitia and Llamas47]) or by utilising 3D RANS results for a swept wing with infinite span (as in ref. [Reference Parenteau, Sermeus and Laurendeau55]). While leaving the mathematics of the lifting-line model unchanged, the approach has been reported to significantly improve the 3D prediction accuracy for high angles of attack and for the wing tip region of swept wings.

3.0 NON-LINEAR VORTEX LATTICE METHOD

The VLM^{(Reference Hedman16,Reference Katz and Plotkin17)} represents a powerful tool for preliminary wing design and optimisation. It has been used in a very wide range of applications, from multi-objective optimisation studies for existing commercial aircraft^{(Reference Leifsson, Slawomir and Bekasiewicz19,Reference Mariens, Elham and van Tooren25)} , the development of morphing wings^{(Reference Smith, Ajaj, Isikveren and Friswell20)}, Unmanned Aerial Vehicles (UAV) aerodynamic performance optimisations^{(Reference Tianyuan and Xiongqing21)}, the design of non-conventional blended wing–body aircraft geometries^{(Reference Peifeng, Zhang, Yingchun, Changsheng and Yu22)}, up to unsteady variants of the method (UVLM) used to calculate aerodynamic loads for aeroelasticity and flight dynamics simulations^{(Reference Murua, Palacios and Graham18)}.

Non-linear extensions to the VLM have been proposed in literature. It must be noted that the focus of this paper (and thus of the papers reviewed here) is on a VLM model based on a non-linear correction of the circulation so as to improve the prediction accuracy and not non-linearities introduced due to wake relaxation models (such as in ref. [Reference Lee and Park50]). In ref. [Reference Dos Santos and Marques45], a model was developed to predict the stall and post-stall characteristics of aircraft wings. It used two-dimensional RANS simulations to predict the boundary-layer separation point and the Kirchhoff flow approach to model the non-linear lift variation of the RANS solution in the VLM. The model obtained improved results compared with other approaches such as iterative de-cambering. The stall behaviour of a horizontal tail was investigated in ref. [Reference Goitia and Llamas47] using a non-linear VLM. The approach used van Dam’s loosely coupled algorithm based on the effective angle-of-attack correction but incorporated 2.5D aerofoil characteristics based on the XFOIL solver and wing sweep theory.

Recently, it was shown in ref. [Reference Jo, Jardin, Gojon, Jacob and Moschetta48] than nonlinear VLM results can be successfully used to predict the base flow field with sufficient accuracy for aero-acoustic noise calculations for micro air vehicle rotors. The method used an iterative correction of the sectional circulation values based on a lookup table of aerofoil viscous XFOIL or CFD results, the wake being modelled as vortex particles rather than rings to reduce the computational effort. An unsteady non-linear VLM was developed and applied to the analysis of wind turbine rotors in ref. [Reference Kim, Lee and Lee49]. At each time step, the linear lift force values of each wing strip were first determined, then the corrections to the circulation values were obtained by solving the non-linear system of equations resulting from the difference between the linear lift forces and the non-linear lift forces extracted from lookup table results. This method can be considered as more coupled compared with other approaches, as the viscous corrections are determined by solving an implicit system of equations. Another model was proposed recently in ref. [Reference Lee and Lee51]. The model also utilises a loosely coupled approach in which an initially determined set of linear circulation values for the vortex rings is iteratively corrected based on lookup table lift coefficient values, from which a lookup table circulation was determined and compared with the linear circulation to determine the correction magnitude. The importance of carefully defining the control point location is highlighted, and a method for updating the location is proposed based on a strip-wise averaging of the circulation followed by an enforcement of the flow tangency condition in a manner similar to the classic thin-aerofoil theory. Results showed very good agreement with experimental data for the MEXICO wind turbine.

Other non-linear VLM approaches were developed in refs. [Reference Parenteau, Sermeus and Laurendeau55,Reference Paul and Gopalarathnam56] based on the effective angle-of-attack updated strategy of van Dam and utilising 2.5D RANS sectional data to increase the accuracy for swept wings (as in ref. [Reference Parenteau, Sermeus and Laurendeau55]) or utilising an iterative de-cambering approach implemented in a non-linear Newton–Raphson iteration scheme so as to drive the potential flow results towards viscous data saved in lookup tables^{(Reference Paul and Gopalarathnam56)}.

In the classical VLM approach, the unknown intensities of all the vortex rings distributed over the wing surface are determined by requiring that the flow tangency condition be satisfied for all collocation points, leading to the following linear system of equations:

(7) \begin{equation}\mathop \sum \nolimits_{j=1}^N {\textbf{\textit{v}}_{ij}}\cdot{\textbf{\textit{n}}_i}{\Gamma _j} = - {\textbf{\textit{V}}_\infty }\cdot{\textbf{\textit{n}}_i}\quad i=1,2, \ldots ,N\end{equation}

In Equation (7), ${{\bf{V}}_\infty }$ is the freestream velocity, $N$ is the total number of vortex rings over the wing surface, ${{\bf{v}}_{ij}}$ is the velocity induced by the unit strength vortex ring $j$ at the ${i^{{\rm{th}}}}$ panel collocation point and ${{\bf{n}}_i}$ is the surface normal vector calculated at the ${i^{{\rm{th}}}}$ panel collocation point.

Like the non-linear LLT method presented in the previous section, the predictive capabilities of the VLM can be enhanced by introducing, to a limited extent, viscous effects^{(Reference Sugar Gabor, Koreanschi and Botez23,Reference Sugar Gabor, Koreanschi, Botez, Mamou and Mebarki24)} . For each vortex ring, a correction $\Delta\Gamma $ is defined, so that the final values of the vortex intensities become

(8) \begin{equation}{\Gamma _j} \to {\Gamma _j} + \Delta{\Gamma _j}\quad j=1,2, \ldots ,N\end{equation}

The approach developed is conceptually similar to that presented in ref. [Reference Kim, Lee and Lee49]. The corrections are calculated by solving a non-linear system of equations constructed based on the available 2D viscous results. Unlike the non-linear LLT, the inviscid–viscous coupling is looser, but the results showed good convergence without the need for significant under-relaxation. This is attributed to the fact that the correction values are determined in a highly coupled approach, the specific $\Delta{\Gamma _j}$ value for each ring being determined as a nonlinear function of all other corrections and of the 2D viscous results for all the wing strips. This approach alleviates the need for significant under-relaxation (such as in ref. [Reference Goitia and Llamas47] or ref. [Reference Lee and Lee51]), where the vortex strength corrections are determined independently for each wing strip.

The flow tangency boundary condition of Equation (7) is modified through the introduction of a surface transpiration velocity ${\bf{V}}_i^T$:

(9) \begin{equation}\left( {{{\bf{V}}_\infty } + \mathop \sum \nolimits_{j=1}^N \left( {{\Gamma _j} + \Delta{\Gamma _j}} \right){{\bf{v}}_{ij}} + {\bf{V}}_i^T} \right)\cdot{{\bf{n}}_i}=0\end{equation}

The effective angle-of-attack seen by each wing strip is determined as

(10) \begin{equation}{\alpha _i} = {\tan ^{ - 1}}\left( {\frac{{{{\bf{V}}_i}\cdot{\bf{n}}{{\bf{s}}_i}}}{{{{\bf{V}}_i}\cdot{\bf{c}}{{\bf{s}}_i}}}} \right) = {\tan ^{ - 1}}\left( {\frac{{\left( {{{\bf{V}}_\infty } + \mathop \sum \nolimits_{j=1}^N \left( {{\Gamma _j} + \Delta{\Gamma _j}} \right){{\bf{v}}_{ij}} + {\bf{VS}}_i^T} \right)\cdot{\bf{n}}{{\bf{s}}_i}}}{{\left( {{{\bf{V}}_\infty } + \mathop \sum \nolimits_{j=1}^N \left( {{\Gamma _j} + \Delta{\Gamma _j}} \right){{\bf{v}}_{ij}} + {\bf{VS}}_i^T} \right)\cdot{\bf{c}}{{\bf{s}}_i}}}} \right)\end{equation}

where ${\bf{n}}{{\bf{s}}_i}$ and ${\bf{c}}{{\bf{s}}_i}$ are the wing strip unit normal and unit chord vectors, ${{\bf{v}}_{ij}}$ is the velocity induced by the unit strength vortex ring $j$ at the ${i^{{\rm{th}}}}$ strip collocation point (located at the strip three-quarter chord) and ${\bf{VS}}_i^T$ is the average transpiration velocity for the strip, calculated as a simple arithmetic average of the transpiration velocities of all panels of that strip, ${\bf{VS}}_i^T = \frac{1}{{{N_x}}}\mathop \sum \nolimits_{j=1}^{{N_x}} {\bf{V}}_j^T$. It is interesting to observe in Equation (10) that the proposed approach amounts to two corrections to the effective angle-of-attack: one due to the viscous circulation corrections, and the second due to the transpiration velocities which ensure that the flow tangency condition is always enforced at the panel collocation points.

Pre-built databases of experimental or numerical aerofoil results are again used, but this time containing pressure coefficient distributions for several Reynolds number values ${\rm{R}}{{\rm{e}}_i}$ and angles of attack ${\alpha _i}$:

(11) \begin{equation}CP_i^{visc} = f\!\left( {air\ foi{l_i},{\rm{R}}{{\rm{e}}_i},{\alpha _i}} \right)\end{equation}

These pressure coefficient distributions are pre-calculated at a given number of flow conditions, while for any other flow condition not in the database, a simple linear interpolation is used between the two closest data points available.

The equations needed to calculate the vortex rings’ intensity corrections are constructed based on the assumption that, for all $N$ panels on the wing surface, the pressure coefficient variation ($\Delta CP = C{P_{upper - surface}} - C{P_{lower - surface}})$ obtained from the vortex rings’ intensities is equal to the nonlinear viscous pressure coefficient variation $\Delta C{P^{visc}}$ obtained from the database. For all panels, the following equality is written:

(12) \begin{equation} - {{\bf{F}}_i}\cdot{{\bf{n}}_i} + {A_i}{Q_\infty }\Delta CP_i^{visc}=0\quad i=1,2, \ldots ,N\end{equation}

where ${{\bf{F}}_i}$ is the aerodynamic force generated by all the vortex lines placed on the panel, ${{\bf{n}}_i}$ is the surface normal vector calculated at the panel collocation point, ${A_i}$ is the panel area and ${Q_\infty }$ is the freestream dynamic pressure.

Coupling Equations (9) and (12) leads to a nonlinear system of $2N$ equations for the vortex strength corrections $\Delta{\Gamma _j}$ and the surface transpiration velocities ${\bf{V}}_i^T$, shown in Equation (13) below, which can be solved using Newton’s method (since the Jacobian matrix can be analytically determined without much difficulty):

(13) \begin{equation}{\bf{R}} =\left\{ {\begin{array}{c}{\begin{array}{{c}} \vdots \\{ - {{\bf{F}}_i}\cdot{{\bf{n}}_i} + {A_i}{Q_\infty }\Delta CP_i^{visc}}\\ \vdots \end{array}}\\{ - - - - - - - - - - }\\{\begin{array}{{c}} \vdots \\{\mathop \sum \nolimits_{j=1}^N {{\bf{v}}_{ij}}\cdot{{\bf{n}}_i}\Delta{\Gamma _j} + V_i^T}\\ \vdots \end{array}}\end{array}} \right\} ={\bf 0}\end{equation}

The expression for the normal projection of the aerodynamic force appearing in Equations (12) and (13) is determined using the three-dimensional vortex lifting law^{(Reference Saffman11)} and follows the notations introduced in Fig. 3.

(14) \begin{align}&{{\mathbf F}}_i\cdot {{\mathbf n}}_i=\rho\! \left[\!{{\mathbf n}}_i\times \!\left(\!{{\mathbf V}}_{\infty }{+}\sum^N_{j=1}{\left({\Gamma }_j+{\triangle \Gamma }_j\right)\!{{\mathbf v}}_{ij}}\!\right)\!\right]\cdot [\left({\Gamma }_i-{\Gamma }_U\right){\boldsymbol{ \gamma }}_{12}+\left({\Gamma }_i-{\Gamma }_R\right)\!{\boldsymbol{\gamma }}_{23}+\left({\Gamma }_i-{\Gamma }_L\right)\!{\boldsymbol{\gamma }}_{61}\nonumber\\[2pt] &\quad +\left({\Gamma }_U-{\Gamma }_{UR}\right){\boldsymbol{\gamma }}_{34}+\left({\Gamma }_U-{\Gamma }_{UL}\right){\boldsymbol{\gamma }}_{56}+\left({\triangle \Gamma }_i-{\triangle \Gamma }_U\right){\boldsymbol{\gamma }}_{12}+\left({\triangle \Gamma }_i-{\triangle \Gamma }_R\right){\boldsymbol{\gamma }}_{23}\nonumber\\[3pt] &\quad +\left({\triangle \Gamma }_i-{\triangle \Gamma }_L\right){\boldsymbol{\gamma }}_{61}+\left(\triangle {\Gamma }_U-{\triangle \Gamma }_{UR}\right){\boldsymbol{\gamma }}_{34}+\left(\triangle {\Gamma }_U-{\triangle \Gamma }_{UL}\right){\boldsymbol{\gamma }}_{56}]\ \ \ \ \ i=1,2,\dots ,N\end{align}

where $\Gamma $ is the strength of a vortex ring as obtained in an initial, purely inviscid solution from solving Equation (7) and ${\bf{\gamma }}$ is the supporting geometric segment of any given vortex line.

Figure 3. Illustration of the vortex segments on a typical panel.

Once the nonlinear system (13) is solved and the corrected vortex ring strengths are determined from Equation (8), the aerodynamic force and moment acting on each vortex segment on the wing can be determined using equations very similar to Equations (4) and (5), and the profile drag coefficient can be determined using Equation (6).

A verification and validation test is done using geometrical and experimental data from the NACA Technical Note 1208^{(Reference Schneider26)}. The wing has a sweep angle of $45^\circ $, an aspect ratio of 8 and a taper ratio of 0.4. The experimental results were obtained for an airspeed of 65m/s and a Reynolds number equal to $4 \times {10^6}$, as calculated with the mean aerodynamic chord value. The numerical results are obtained with a mesh of 18 chordwise panels and 35 spanwise panels per wing semi-span, while the sectional aerofoil database is generated by XFOIL. For reference purposes, linear VLM results obtained with the widely used XFLR5 code are also included in the comparison.

In Fig. 4, the results for the wing lift coefficient and quarter chord pitching moment coefficient are compared with the experimental data. The nonlinear VLM predicts the lift curve characteristics very well, with a slight overestimation in predicted ${C_L}$ values for angles of attack higher than $10^\circ $. A very good agreement exists for (1.01 for the experiment, versus 1.04 in the numerical results), but there is an underestimation of the stall angle by about $1.5^\circ $. The linear variation of the pitching moment coefficient is very well captured, but there are some differences for the nonlinear higher lift conditions, where the swept-back wing experiences an early tip stall phenomenon which is difficult to capture. Figure 5 shows a comparison of the span-wise wing loading for an angle of attack of $4.7^\circ $, the agreement being good.

Figure 4. Comparison of lift and pitching moment coefficient values obtained with the non-linear VLM method and the XFLR5 VLM code with experimental data published for the NACA TN1208 high-aspect-ratio high-sweep wing geometry.

Figure 5. Comparison of span-wise loading obtained with the non-linear VLM method with experimental data published for the NACA TN1208 high-aspect-ratio high-sweep wing geometry.

4.0 CASE STUDY: AERODYNAMIC SHAPE OPTIMISATION OF A GENERIC UAV WING

The aerodynamic shape optimisation of a generic UAV is performed using the nonlinear VLM as solver. The wing is chosen to be representative for a tactical UAV with an all-up mass under 20kg, having an aspect ratio of 8, a span of just over 4 m, a taper ratio of 0.5 and a NACA 2412 aerofoil section. The optimisation is done using the Artificial Bee Colony (ABC) algorithm, as described in ref. [Reference Sugar Gabor, Simon, Koreanschi and Botez27], with the objective of improving the wing lift-to-drag ratio $L/D$ over a range of fixed angle-of-attack values, at an airspeed of 50m/s and with a Reynolds number of $2.13 \times {10^6}$, as calculated with the mean aerodynamic chord of the original wing geometry. The optimisation procedure is focused on aerodynamic performance, so no structural or weight aspects are considered.

Two optimisation cases are considered. The first case uses the wingspan, taper ratio, and quarter-chord sweep angle as the optimisation variables. The wingspan is constrained between 3.5m and 5m, the quarter-chord sweep angle between $0^\circ $ and $30^\circ $ and the taper ratio between 0.3 and 1. The second case changes the aerofoil shape in addition to the wing planform variables indicated for the first case. To achieve this, a Non-Uniform Rational B-Splines (NURBS)^{(Reference Piegl and Tiller28)} parameterisation of the NACA 2412 is introduced:

(15) \begin{align}&\hspace*{-7pt}{\bf{C}}\!\left(u\right) = \mathop \sum \nolimits_{i=1}^k \frac{{{N_{i,n}}{w_i}}}{{\mathop \sum \nolimits_{j=1}^k {N_{j,n}}{w_j}}}{{\bf{P}}_i}\nonumber\\[3pt] &{N_{i,1}} = \left\{ {\begin{array}{*{20}{l}}{1,\quad if{t_i} \le u \le {t_{i+1}}}\\{0,\quad otherwise}\end{array}} \right.\nonumber\\[3pt] &{N_{i,n}} = \frac{{u - {t_i}}}{{{t_{i + n}} - {t_i}}}{N_{i,n - 1}} + \frac{{{t_{i + n+1}} - u}}{{{t_{i + n+1}} - {t_i}}}{N_{i+1,n - 1}}\end{align}

In Equation (15), ${\bf{C}}\left( u \right)$ is the parametrised form of the aerofoil curve, $u$ is the curve parameter, ranging from 0 (the start of the curve) to 1 (the end of the curve), $k$ is the number of control points, $n$ is the order of the curve, ${w_i}$ are the weights associated with the control points, ${t_i}$ are the knots, ${N_{i,n}}$ are the basis functions and ${{\bf{P}}_i} = \left[ {{x_i},{y_i},{z_i}} \right]$ are the control points.

The initial coordinates of the control points ${{\bf{P}}_i}$ are determined through a non-linear regression in which the NURBS aerofoil shape generated with Equation (14) is iteratively driven towards the actual NACA 2412 aerofoil shape. In the numerical optimisation, the change of the aerofoil shape is achieved by changing the coordinates of the NURBS control points, which act as optimisation variables. Only the upper surface of the aerofoil between the leading edge and $0.5c$ can change, and the control points are constrained to move only on the y-axis and by a maximum of $0.05c$. To capture the effects of the aerofoil shape change on the wing performance, pre-built databases of 2D results are not used in this scenario. Rather, the XFOIL solver is called while the optimisation procedure runs during each iteration of the non-linear VLM system solution and for each span-wise wing strip.

Figure 6 presents a comparison between the original and redesigned wing shapes, as well as between the original and optimised (morphed) airfoil, while Fig. 7 presents a comparison between the lift curve and drag polar for the original wing, the planform-only optimised wing and the planform and aerofoil optimised wing.

Figure 6. Comparison between the original and the redesigned wing and airfoil shapes.

It must be noted that the planform shape of the optimised wing in the two cases was approximately the same, regardless of whether the aerofoil optimisation variables were included in the process or not. The results hold no surprises, as the case study is meant to be more of a test of the sensitivity of the nonlinear VLM with respect to aerofoil shape changes rather than a rigorous optimisation scenario. As expected, the optimal values of the planform variables provide a higher aspect ratio, lower sweep, lower taper wing, which achieves a higher $L/D$ than the original design at any ${C_L}$ value. This result is visible in both the increased lift curve slope and the reduced drag at any fixed lift coefficient value. The second optimisation case, in which the parameterised aerofoil shape is added to the optimisation variables, provides more significant results. As can be seen from Fig. 7, the nonlinear VLM is sensitive enough to capture the effects of small changes in the aerofoil thickness on the aerodynamic performance of the wing. This represents a net advantage over the classic VLM, which considers the wing to be a zero-thickness surface. The small decrease in aerofoil thickness between the original and optimised (morphed) aerofoil shapes over $0.3c$ amounts to a further reduction in drag over the entire range of angles of attack considered. This further reduction is only due to changes in the wing profile drag, as the planform shape obtained is essentially unchanged between the two optimisation cases.

Figure 7. Lift curve and drag polar for the original wing, the redesigned wing using only wing planform parameters and the redesigned wing including aerofoil shape optimisation.

The wing profile drag is calculated using Equation (6), which means that drag coefficient reductions could potentially be obtained only due to the local drag coefficient reductions for the aerofoil sections, as calculated with the XFOIL solver. In order to show that the circulation values calculated by the VLM are sensitive to the changes in the aerofoil shape, a plot of the local lift coefficient along the wingspan at $\alpha =2^\circ $ is depicted in Fig. 8 for both the planform-only optimised wing and the planform and aerofoil optimised wing. The local lift is determined for each spanwise wing strip by summing the contributions of all the vortex rings distributed chordwise along the strip. Differences can be observed in the calculated local lift coefficient, indicating that the circulation strength of the vortex rings does change in response to aerofoil shape changes.

Figure 8. Local lift coefficient calculated using the corrected vortex ring strengths for an angle-of-attack of 2°.

5.0 CASE STUDY: INTEGRATION INTO A SIMULTANEOUS ANALYSIS AND DESIGN METHODOLOGY

Using global optimisation procedures such as the ABC algorithm or Genetic Algorithm (GA) can provide the global optimum point in the design space, however doing so requires a very large number of objective function evaluations and thus runs of the aerodynamic solver. Additionally, any increase in the number of design variables leads to further, significant increases in the number of evaluations, to a point where even using relatively simple methods such as the nonlinear LLT or VLM becomes much too time-consuming for conceptual design purposes.

Gradient-based optimisation algorithms can arrive at the global optimum point in a deterministic way and with a much smaller number of evaluations, provided the objective function is differentiable and a suitable method for evaluating the gradient is found. The adjoint equation approach has become particularly popular in the field of aerodynamic shape optimisation because it allows the evaluation of the objective function gradient with a computational cost roughly equal to one additional aerodynamic solver run, regardless of the dimension of the design space. This approach has been extensively used over recent years, for a wide variety of problems across various industries. These include helicopter rotors^{(Reference Fabiano, Mishra, Mavriplis and Mani29)}, automobiles^{(Reference Blacha, Gregersen, Islam and Bensler30)}, wide-body transport aircraft^{(Reference Kenway and Martins31)}, supersonic aircraft configurations^{(Reference Reuther, Alonso, Rimlinger and Jameson32)}, tidal turbines^{(Reference Funke, Farrell and Piggott33)} and ship hulls^{(Reference Ragab34)}.

The nonlinear LLT and VLM methods presented above are particularly well suited for a gradient-based adjoint equation implementation since all the required matrices can be analytically calculated and efficiently implemented. Details of a methodology for further reducing the computational cost of a discrete adjoint method, with the construction of an adjoint-based simultaneous analysis and design (SAND) technique, are presented in detail in ref. [Reference Sugar Gabor35].

To minimise the objective functional $J(\textbf{\textit{w,}}\ \boldsymbol\alpha),J:\mathbb{R}^{N}\times\mathbb{R}^{D}\rightarrow\mathbb{R}$ subject to the equality and inequality constraints given by $\textbf{\textit{G}}(\textbf{\textit{W}},\boldsymbol{\alpha})\leq {\bf 0}, \textbf{\textit{G}}:\mathbb{R}^{N}\times\mathbb{R}^{D}\rightarrow\mathbb{R}^{K}$ and to satisfy the system of equations $\textbf{\textit{R}}(\textbf{\textit{W}},\boldsymbol{\alpha})= {\bf 0}, \textbf{\textit{R}}:\mathbb{R}^{N}\times\mathbb{R}^{D}\rightarrow\mathbb{R}^{N}$ that models the aerodynamic problem, a Lagrange functional can be defined as^{(Reference Cacuci36)}

(16) \begin{equation}L\left( {\textbf{\textit{w}},{\boldsymbol\psi _{\textbf{1}}},\boldsymbol\alpha ,{\boldsymbol\psi _{\textbf{2}}}} \right) = J\left( {\textbf{\textit{w}},\boldsymbol\alpha } \right) + {\left( {{\boldsymbol\psi _{\textbf{1}}},\textbf{\textit{R}}\left( {\textbf{\textit{w}},\boldsymbol\alpha } \right)} \right)_N} + {\left( {{\boldsymbol\psi _{\textbf{1}}},\textbf{\textit{G}}\left( {\textbf{\textit{w}},\boldsymbol\alpha } \right)} \right)_K}\end{equation}

where $\textbf{\textit{w}}\in\mathbb{R}^{N}$ are the system-dependent variables and $\boldsymbol\alpha\in\mathbb{R}^{D}$ are the system design parameters, whose values must be given to determine a solution for the system of equations, while $\boldsymbol\psi_{\textbf{1}}\in\mathbb{R}^{N}$ and $\boldsymbol\psi_{\textbf{2}}\in\mathbb{R}^{K}$ are two sets of Lagrange multipliers which will also play the role of the adjoint variables, and ${\left( {^\circ ,^\circ } \right)_N}$ signifies an appropriately defined inner product in the $\mathbb{R}^{N}$ space.

As shown in ref. [Reference Sugar Gabor35], introducing the adjoint ${\textbf{\textit{R}}^ + }$ for the model equations and ${\textbf{\textit{G}}^ + }$ for the set of equality and inequality constraints leads to the following set of SAND equations:

(17) \begin{equation}\frac{{\partial J}}{{\partial \textbf{\textit{w}}}} + {\frac{{\partial \textbf{\textit{R}}}}{{\partial \textbf{\textit{w}}}}^ + }{\boldsymbol\psi _{\textbf{1}}} + {\frac{{\partial \textbf{\textit{G}}}}{{\partial \textbf{\textit{w}}}}^ + }{\boldsymbol\psi _{\textbf{2}}}={\bf 0}\end{equation}

(18) \begin{equation}\textbf{\textit{R}}=\textbf{0}\end{equation}

(19) \begin{equation}\frac{{\partial J}}{{\partial \boldsymbol\alpha }} + {\frac{{\partial \textbf{\textit{R}}}}{{\partial \boldsymbol\alpha }}^ + }{\boldsymbol\psi _{\textbf{1}}} + {\frac{{\partial \textbf{\textit{G}}}}{{\partial \boldsymbol\alpha }}^ + }{\boldsymbol\psi _{\textbf{2}}}=\textbf{0}\end{equation}

(20) \begin{equation}{\left( {\textbf{\textit{G}},\boldsymbol\delta {\boldsymbol\psi _{\textbf{2}}}} \right)_K}=0\end{equation}

The solution of the above system, ${\left( {\textbf{\textit{w}}},{\boldsymbol\psi _{\textbf{1}}},\boldsymbol\alpha ,{\boldsymbol\psi _{\textbf{2}}} \right)^T}$, includes the constrained optimal values of the design parameters as required to minimise the objective functional $J$, the solution of the system of equations $\textbf{\textit{R}}$ modelling the aerodynamic problem as obtained with the optimal design parameters, and the values for the two vectors of adjoint variables. The Jacobian matrix of the system of Equations (17–20) can become ill conditioned or singular at points, thus requiring a trust-region method^{(Reference Byrd, Gilbert and Nocedal37)} for its solution.

The nonlinear LLT has been integrated into the SAND approach and applied to the problem of optimising the geometry of a winglet fitted to a generic UAV wing geometry whose details are presented in Table 1. Morphing wing tip devices^{(Reference Gomes and Suleman38,Reference Smith, Ajaj, Isikveren and Friswell39)} have attracted much research as a promising solution that could efficiently increase a wing’s lift-to-drag ratio. The optimisation variable is the winglet toe angle, and the objective function is improving the wing lift-to-drag ratio ${C_L}/{C_D}$ for three angle-of-attack values, namely –3, 0$^\circ $ and 3$^\circ $, at an airspeed of 10m/s. The toe angle is constrained between a lower limit of ${\tau _{min}} = - 10^\circ $ and an upper limit of ${\tau _{max}}=10^\circ $. The wing is modelled using 60 horseshoe vortices, clustered towards the wing tips, while 20 horseshoe vortices are used for each winglet, clustered towards both the winglet tip and the junction with the wing tip. The variation of the two-dimensional lift coefficient is limited to the linear region, with the lift curve slope and the zero-lift angle-of-attack for both NACA 4412 and NACA 4409 aerofoils being estimated based on experimental data. The drag computations are limited to the induced drag component only, since 2D profile drag data are not considered.

Table 1 Wing and winglet geometry parameters

Table 2 presents a comparison between the lift-to-drag ratio for the original design and the optimised design, together with the toe angle values that the morphing winglets take at each different angle-of-attack to achieve the indicated performance increase. A successful increase in ${C_L}/{C_D}$ has been achieved for all three angles of attack. The trust-region algorithm required 7, 12 and 10 iterations to achieve convergence in the three conditions considered. In terms of computational effort, only a very small number of evaluations is required with the SAND approach, demonstrating the high efficiency sought in the conceptual design stage, where hundreds of different configurations might be tested.

Table 2 Performance improvements obtained using morphing winglet

6.0 EXTENSIONS TO UNSTEADY FLOWS

The nonlinear lifting line model summarised earlier and presented in detail in refs. [Reference Sugar Gabor, Koreanschi and Botez9,Reference Sugar Gabor, Koreanschi and Botez10] was extended to unsteady flows in ref. [Reference Sugar Gabor40]. In the context of unsteady flows, the continuous distribution of bound vorticity over the lifting surface and of trailing vorticity in the wake are approximated using a finite number $N$ of ring vortices bound to the geometry, and at each time step, a new row of $N$ vortex rings are shed into the wake, as illustrated in Fig. 9.

Figure 9. Sketch of the unsteady trailing vortex system.

To determine the force acting on a bound vortex segment, the following unsteady form of the vector Kutta–Joukowski theorem is introduced (a full derivation can be found in ref. [Reference Sugar Gabor41]):

(21) \begin{equation}{\bf{dF}} = \rho \Gamma \left( {{\bf{V}} \times {\bf{dl}}} \right) + \rho c\frac{{\partial \Gamma }}{{\partial t}}{\bf{n}} + \rho c\Gamma \frac{{d{\bf{n}}}}{{dt}}\end{equation}

Equation (21) is written for the quarter-chord vortex segment of all vortex rings placed over the lifting surface. The magnitude of the aerodynamic force acting on a span-wise strip was already presented in Equation (3). Following the same idea as for the steady-state case (obtaining the 2D aerofoil aerodynamic characteristics using either experimental data or numerical solvers), the modulus of the force in Equation (21) can be set equal to the modulus of Equation (2), and after some algebraic manipulation, the following equation is obtained:

(22) \begin{align}&\rho {\Gamma _i}\sqrt {{{\left[ {\left( {{{\bf{V}}_i} \times {\bf{d}}{{\bf{l}}_i}} \right)\cdot{{\bf{n}}_i}} \right]}^2} + {{\left[ {\left( {{{\bf{V}}_i} \times {\bf{d}}{{\bf{l}}_i}} \right)\cdot{{\bf{c}}_i}} \right]}^2}} + \rho {c_i}{\left( {\frac{{\partial \Gamma }}{{\partial t}}} \right)_i} + \rho {c_i}{\Gamma _i}\sqrt {{{\left[ {\frac{{d{{\bf{n}}_i}}}{{dt}}\cdot{{\bf{n}}_i}} \right]}^2} + {{\left[ {\frac{{d{{\bf{n}}_i}}}{{dt}}\cdot{{\bf{c}}_i}} \right]}^2}}\nonumber\\ &\quad = \frac{1}{2}\rho \left[ {{{\left( {{{\bf{V}}_i}\cdot{{\bf{n}}_i}} \right)}^2} + {{\left( {{{\bf{V}}_i}\cdot{{\bf{c}}_i}} \right)}^2}} \right]d{A_i}{C_l}_i\quad i=1, \ldots ,N\end{align}

The local airspeed ${{\bf{V}}_i}$ includes contributions from the local kinematic velocity, the velocities induced by the vortex segments on the wing surface and by the vortex rings shed in the wake and is given by

(23) \begin{equation}{{\bf{V}}_i} = - \left( {{{\bf{V}}_\infty } + {{\bf{v}}_{rel}}_{\textbf{\textit{i}}} + {\boldsymbol{\Omega }} \times {{\bf{r}}_i}} \right) + \mathop \sum \nolimits_{j=1}^N \Gamma _j^n{{\bf{w}}_{ij}} + \mathop \sum \nolimits_{k=2}^M \mathop \sum \nolimits_{j=1}^N \Gamma _j^{n - k+1}{{\bf{w}}_{ikj}}\end{equation}

where ${{\bf{v}}_{rel}}_{\textbf{\textit{i}}}$ is the local relative velocity of a point on the wing with respect to the wing-fixed reference system (non-zero in the case of flapping wings, for example), ${\boldsymbol{\Omega }}$ is the angular velocity of the wing-fixed reference system with respect to the ground-fixed reference system, ${{\bf{r}}_i}$ is the position vector of a point on the wing and $M$ is the number of rows of vortex rings shed in the wake.

By inserting Equation (23) into Equation (22) and estimating the time derivative using a second-order backwards difference, the following time-dependent nonlinear system of equations is determined:

(24) \begin{equation}{R_i}\left( {{{\rm{\Gamma }}^n}} \right) = \left( {{E_i}\left( {{\Gamma ^{\textbf{n}}}} \right) + \frac{{3{G_i}}}{{2{\rm{\Delta }}t}}} \right)\Gamma _i^n - {F_i}\left( {{\Gamma ^n}} \right) - \frac{{4{c_i}}}{{2{\rm{\Delta }}t}}\Gamma _i^{n - 1} + \frac{{{c_i}}}{{2{\rm{\Delta }}t}}\Gamma _i^{n - 2}=0,i=1, \ldots ,N\end{equation}

The various terms appearing in the above equation are lengthy and will not be presented here but can be found in full in ref. [Reference Sugar Gabor40].

The unsteady calculation procedure begins with updating the wing geometry at the new timestep (if required, in the case of flapping wings, for example). Next, the nonlinear system in Equation (24) is iteratively solved, assuming that the wake shape is frozen in time. Once the new vortex strength values ${{\rm{\Gamma }}^n}$ are converged to a desired precision, the final step requires the shape of the wake to be updated accordingly. This requires an iterative relaxation of the wake surface to ensure a physically representative force-free wake surface, as the current position ${{\bf{X}}^n}$ of each wake point depends on the current position of all other points via the wake-induced velocities ${{\bf{W}}_{kj}}$. To handle the inherent non-linearity of the wake relaxation process, the following scheme is proposed:

(25) \begin{align}{{\bf{X}}^{\textbf{\textit{t}}}} &= {{\bf{X}}^{{\textbf{\textit{t}}} - {\bf 1}}}\nonumber\\[3pt] {{\bf{X}}^{{\textbf{\textit{t}}}+{\bf 1}}} &= {{\bf{X}}^{\textbf{\textit{t}}}} + \left[ {\frac{{{{\bf{X}}^{\textbf{\textit{t}}}} - {{\bf{X}}^{{\textbf{\textit{n}}} - {\bf 1}}}}}{{\Delta {\textbf{\textit{t}}}}} - \left( {\mathop \sum \nolimits_{j=1}^N \Gamma _j^n{{\bf{W}}_j}\left( {{{\bf{X}}^{\textbf{\textit{t}}}}} \right) + \mathop \sum \nolimits_{k=2}^M \mathop \sum \nolimits_{j=1}^N \Gamma _j^{n - k+1}{{\bf{W}}_{kj}}\left( {{{\bf{X}}^{\textbf{\textit{t}}}}} \right)} \right)} \right]{\rm{\Delta \tau }}\nonumber\\[3pt] &\qquad\qquad\qquad\qquad\qquad\qquad\quad\qquad\qquad if\left\|{{\bf{X}}^{{{\textbf{\textit{t}}}+{\bf 1}}}} - {{\bf{X}}^{\bf t}}\right\| \le \varepsilon ,{{\bf{X}}^{\bf n}} = {{\bf{X}}^{{\textbf{\textit{t}}}+{\bf 1}}}\end{align}

Here, ${\rm{\Delta \tau }}$ represents a fictitious time step and $\varepsilon $ is a desired convergence criterion. The scheme in Equation (25) is inspired by the dual time-stepping approach^{(Reference Jameson57)} used in many unsteady CFD solvers. The explicit marching in the fictitious time is set up to guarantee an implicit non-linear approximation at the current physical time. A small number of only three iterations in the fictitious time was found to be sufficient to ensure reasonable convergence of the wake relaxation process to $\varepsilon/ $$b = {10^{ - 2}}$, where $b$ is the wingspan.

To test the capabilities of the method, a rectangular wing undergoing harmonic flapping is analysed. A comparison is made with the Unsteady Vortex Lattice Method (UVLM), which has been repeatedly proven to provide unsteady lift and thrust predictions with relatively high accuracy and at low computational cost for flapping flight^{(Reference Fritz and Long43)}. The geometry chosen has an aspect ratio of 8 and is generated using a highly cambered aerofoil from the NACA 83 series^{(Reference Fritz and Long43)}. UVLM results are available for two reduced frequencies, ${k_w}=0.08$ and ${k_w}=1$, with ${k_w} = \left( {4l{\beta _0}n} \right)/{V_\infty }$, where ${\beta _0}$ is the maximum amplitude of the flapping motion, $l$ is the wing half-span and $n$ is the flapping frequency. The lower-frequency case is representative of a bird the size of a pigeon, while the high-frequency case is representative of insect flight. Research has shown that various flapping regimes each have their associated modelling challenges^{(Reference Ho, Nassef, Pornsinsirirak, Tai and Ho42)}. The numerical simulation of insect flapping is particularly demanding, due to the very complex flow behaviour and the development of highly nonlinear lift and thrust generation mechanisms such as the clap and fling mechanism, rotational lift, wake capture (especially low advance ratio, hovering flight), laminar boundary-layer separation, unsteady leading-edge vortex formation and strong aeroelastic coupling.

It is thus very important to note that the focal point of the comparison is not to show accurate physical modelling but to demonstrate the ability of the unsteady lifting line model to predict the same aerodynamic behaviour as the vortex lattice in a field where it has been only very rarely used. Figures 10 and 11 present the variation of the steady and unsteady lift components during the flapping motion as calculated by the unsteady lifting line and by the UVLM. The low-frequency case is dominated by the steady ${C_L}$ component, while in the high-frequency case, both components are significant, and out of phase. These results agree with the observation that unsteady flapping effects contribute to lift generation only if ${k_w} \ge 0.66$^{(Reference Ho, Nassef, Pornsinsirirak, Tai and Ho42)}. Figures 12 and 13 indicate how the wake development differs qualitatively between the two cases.

Figure 10. Comparison of steady and unsteady lift contributions for the flapping-wing case having a reduced frequency of 0.08.

Figure 11. Comparison of steady and unsteady lift contributions for the flapping-wing case having a reduced frequency of 1.00.

Figure 12. Wake development for the flapping-wing case having a reduced frequency of 0.08.

Figure 13. Wake development for the flapping-wing case having a reduced frequency of 1.00.