1. Introduction
Flying animals display impressive abilities in controlling the flow surrounding them. By flapping and clapping their wings, insects and small birds demonstrate high lift performance and can undertake manoeuvres whose mechanics are not fully understood yet. Such agility is mainly attributed to a leading edge vortex attached to the suction side of the flyer's wings (see Ellington Reference Ellington1984; Dickinson & Gotz Reference Dickinson and Gotz1993) even though other transient phenomena are thought to play a role in their manoeuvrability (see Dickinson, Lehmann & Sane Reference Dickinson, Lehmann and Sane1999). Despite the vast amount of experimental and computational studies on the effect of wake vorticity on the aerodynamic loads of wings and airfoils in unsteady motion (e.g. see Birch & Dickinson Reference Birch and Dickinson2001; Wang, Birch & Dickinson Reference Wang, Birch and Dickinson2004; Taira et al. Reference Taira, Dickson, Colonius, Dickinson and Rowley2007), there is still a need for inexpensive models able to produce results in real time.
In that context, vortex methods have been successfully applied to the prediction of the flow separation from an airfoil at high Reynolds number. Their very formulation indeed makes their use as reduced order models quite seamless for such flows, enabling substantial computational savings with respect to a classical Navier–Stokes solver. This is especially due to the underlying flow discretization of the method: a set of discrete compact vortex elements, whose evolution is governed by the Biot–Savart law (see Giesing Reference Giesing1968; Katz Reference Katz1981; Pullin & Wang Reference Pullin and Wang2004; Li & Wu Reference Li and Wu2015; Xia & Mohseni Reference Xia and Mohseni2017).
As time progresses though, the flow is represented by a growing number of vortex elements, which requires the resolution of increasingly larger 
$n$-body problems. This growth in dimensionality makes the investigation of long-term flow behaviour impractical.
To alleviate the increase in complexity, Tchieu & Leonard (Reference Tchieu and Leonard2011) and Wang & Eldredge (Reference Wang and Eldredge2013) proposed to avoid the production of a vortex element at each time step by a process called impulse matching, i.e. feeding vorticity into already existing point vortices, keeping the most coherent vortical structures only. Because the natural instability of shear layers favours the formation of larger scale structures, such models fail at representing flows at long horizon times with accuracy. Another technique relies on merging a pair of neighbouring vortices into one another. In the work of Xia & Mohseni (Reference Xia and Mohseni2013), the merging operation is carried out under the condition that the Taylor expansions of both near- and far-field are not altered by the amalgamation. This double condition translates into a modification of both position and velocity of the point vortex resulting from the merging operation. More recently, Darakananda & Eldredge (Reference Darakananda and Eldredge2019) corrected the velocity of the resulting vortex to exactly cancel the spurious force induced by the local change of circulation around the vortex pair. They qualified their model as hybrid since the proportion of the vortex sheet to lump into the roll-up point vortex depends on a parameter left at the discretion of the user. Even though this model seems more attractive because it does not rely on Taylor expansions of the velocity field, a spurious moment on the flat plate appears. With this framework, spurious force and moment cannot be cancelled simultaneously.
All of the aforementioned efforts have led to lower-order models that are applicable to the case of a flat plate. To the best of the authors’ knowledge, the sole effort towards the handling of a generic airfoil geometry is the derivation by Eldredge (Reference Eldredge2019) of a complex analysis expression for vorticity amalgamation (see (7.45) of aforementioned reference). The present work aims at filling this gap through a work on the real-valued expressions of impulse matching and from there, the construction of a numerical scheme capable of handling generic airfoil geometries with a finite angle trailing edge. For the sake of straightforwardness, we only treat the case of a single shedding source, i.e. at the trailing edge; this allows us to base our shedding scheme on the unsteady Kutta condition derived by Xia & Mohseni (Reference Xia and Mohseni2017). Finally, as the target applications lie in biolocomotion and unsteady aerodynamics, we derive control-volume based expressions for the forces that are more numerically stable than formulations involving the budget of fluid impulse.
The remainder of the paper is organized as follows. In § 2, we present the methodology used for the unsteady panel method. Section 3 addresses the transfer of circulation between the wake vortex sheet and the roll-up vortex generalized for flows around airfoils with a wedge-shaped trailing edge. Section 4 shows the validation of our framework through comparisons with computational and experimental results: cases of impulsively started and heaving and pitching airfoils are studied. We summarize and discuss the results in § 5.
2. Unsteady panel method
We develop the mathematical tools necessary to build the flow model for wake-shedding airfoils. The present model is closely related to the framework proposed by Xia & Mohseni (Reference Xia and Mohseni2017).
2.1. Model formulation
Let us consider a region of inviscid fluid 
$\mathcal {V}_f$ internally bounded by a body surface 
$\mathcal {S}_b$. The body is moving in this volume with a velocity 
$\boldsymbol {u}_b$ and thereby generates vorticity: an infinitely thin boundary layer forms at the body surface and a wake is shed from its trailing edge of finite angle 
$\theta _{TE}$. Vorticity is compact and captured by a set of three types of vortex elements. They are:
(i) a collection of
$N_p$ panels with linear strength located on $\mathcal {S}_b$ to model the infinitely thin boundary layer;(ii) a unique panel in the very near wake
$\mathcal {S}_s$ with length $b_s$ and uniform strength $\gamma _s$;(iii) a series of
$N_v$ point vortices in the wake shed by the airfoil. These point vortices are regularized with the low-order algebraic kernel (see Rosenhead Reference Rosenhead1930)
(2.1)where\begin{equation} \boldsymbol{K}^*(\boldsymbol{x}) =- \frac{\boldsymbol{x}}{2{\rm \pi} (|\boldsymbol{x}|^2 + \delta^2)}, \end{equation}$\delta$ is the blob radius. Figure 1 summarizes this flow discretization where the uniform shed vorticity is unknown and so are the $N_p +1$ bound vortex strengths.
Figure 1. Discretization of an airfoil in an inviscid fluid. The trailing edge and shed panels are cropped by an infinitesimal 
$\varepsilon$ to avoid the singularity that would result from the discontinuity of the vortex sheet strength.
This vorticity field enables the computation of the velocity everywhere in 
$\mathcal {V}_f \,\backslash \,\mathcal {S}_b$ thanks to the Biot–Savart integral which becomes
\begin{align} \boldsymbol{v}(\boldsymbol{x}) &= \int_{\mathcal{S}_b+\mathcal{S}_s} \boldsymbol{K}(\boldsymbol{x} - \boldsymbol{y}) \times \boldsymbol{\gamma}(\boldsymbol{y}) \,\text{d} s(\boldsymbol{y}) + \sum_{k=1}^{N_v}\boldsymbol{K}^*(\boldsymbol{x} - \boldsymbol{x}_k) \times \boldsymbol{\varGamma}_k\nonumber\\ &\quad + \int_{\mathcal{S}_b} [\boldsymbol{K}(\boldsymbol{x}-\boldsymbol{y}) \times (\boldsymbol{n}(\boldsymbol{y})\times \boldsymbol{u}_b(\boldsymbol{y}) - \boldsymbol{K}(\boldsymbol{x}-\boldsymbol{y})(\boldsymbol{n}(\boldsymbol{y})\boldsymbol{\cdot} \boldsymbol{u}_b(\boldsymbol{y}) ) ] \,\text{d} s(\boldsymbol{y}), \end{align}
where 
$\boldsymbol {x}_k$ and 
$\boldsymbol {\varGamma }_k$ are the position and circulation of point vortex 
$k$, respectively; 
$\boldsymbol {K}$ is the singular planar velocity kernel and 
$\boldsymbol {n}$ is the normal pointing towards 
$\mathcal {V}_f$. The limit of this equation as we approach 
$\mathcal {S}_b$ is (see e.g. Eldredge (Reference Eldredge2019))
\begin{align} &- \frac{1}{2} \boldsymbol{n}(\boldsymbol{x})\times\boldsymbol{\gamma}(\boldsymbol{x}) - {\int\hskip -1,05em -\,}_{\mathcal{S}_b} \boldsymbol{K}(\boldsymbol{x}-\boldsymbol{y}) \times \boldsymbol{\gamma}(\boldsymbol{y})\,\text{d} s(\boldsymbol{y}) - \int_{\mathcal{S}_s} \boldsymbol{K}(\boldsymbol{x}-\boldsymbol{y}) \times \boldsymbol{\gamma}_s\,\text{d} s(\boldsymbol{y})\nonumber\\ &\quad = \sum_{k=1}^{N_v} \boldsymbol{K}^*(\boldsymbol{x} - \boldsymbol{x}_k) \times \boldsymbol{\varGamma}_k - \frac{1}{2}\boldsymbol{u}_b(\boldsymbol{x})\nonumber\\ &\qquad+ {\int\hskip -1,05em -\,}_{\mathcal{S}_b} [\boldsymbol{K}(\boldsymbol{x}-\boldsymbol{y}) \times (\boldsymbol{n}(\boldsymbol{y})\times \boldsymbol{u}_b(\boldsymbol{y})) - \boldsymbol{K}(\boldsymbol{x}-\boldsymbol{y})(\boldsymbol{n}(\boldsymbol{y})\boldsymbol{\cdot} \boldsymbol{u}_b(\boldsymbol{y}) ) ] \,\text{d} s(\boldsymbol{y}), \end{align}
where the symbol 
${\int\hskip -1,05em -\,}$ denotes the Cauchy principal value. In the particular case of uniform flows, the Cauchy integral of the right-hand side of (2.3) simplifies to 
$- \frac {1}{2} \boldsymbol {u}_b(\boldsymbol {x})$.
We then can impose at the same time:
(i) the no-flow-through condition at the centre (i.e. collocation point) of the
$N_p$ bound panels. This condition corresponds to the normal component of (2.3);(ii) Kelvin's theorem for conservation of circulation
(2.4)\begin{equation} - \int_{\mathcal{S}_b} \boldsymbol{\gamma} \,\text{d} s = \sum_{k=1}^{N_v} \boldsymbol{\varGamma}_k + \int_{\mathcal{S}_b} \boldsymbol{n}\times\boldsymbol{u}_b \,\text{d} s; \end{equation}(iii) the unsteady Kutta condition (see Xia & Mohseni Reference Xia and Mohseni2017) linking the shedding angle of the shed panel and its strength,
(2.5)where\begin{equation} \gamma_s = \gamma_1 \cos{(\theta_+)} + \gamma_{N_p+1} \cos{(\theta_{TE} - \theta_+)}, \end{equation}$\theta _+$ is the angle between the upper trailing edge panel and the shedding direction given by
(2.6)\begin{equation} \theta_+ = \arccos{\left(\frac{u_{+}^2 + u_3^2 - u_-^2}{2u_{+}u_3}\right)}, \end{equation}$u_+$ and $u_-$ being the fluid velocities in the reference frame of the body above and below the trailing edge, respectively, and $u_3 = -\sqrt {u_{+}^2 + u_-^2 + 2u_{+}u_-\cos {(\theta _{TE}})}$. It is worth noting that vorticity is always shed in the sector defined by the trailing edge panels, as was originally observed by Poling & Telionis (Reference Poling and Telionis1986). In the cases where $u_+=0$ or $u_-=0$, the shedding angle matches the angle of the lower or upper trailing edge side, respectively. In the case of a steady flow, the shedding direction bisects the trailing edge angle since no vorticity is shed in this particular case.
These three conditions lead to a linear system of size 
$N_p+2$ with as many unknowns and can be represented by block matrices:
\begin{equation}
\renewcommand{\arraystretch}{1.5} \left[ \begin{array}{@{}c|c@{}}
\boldsymbol{\mathsf{A}} & \boldsymbol{a}
\\ \hline \boldsymbol{k}^T & 2b_s \\ \hline \boldsymbol{\kappa}^T & -1
\end{array}\right] \left[\begin{array}{@{}c@{}}
\boldsymbol{\gamma} \\ \hline \gamma_s \end{array}\right]
=\left[\begin{array}{@{}c@{}} \boldsymbol{\chi} \\ \hline -
\varGamma \\ \hline 0 \end{array}\right].
\end{equation} The first line enforces the no-flow-through condition with 
$A_{ij}$ the influence of vortex strength 
$j$ and 
$a_i$ the influence of the shed panel at panel 
$i$. On the right-hand side, 
$\chi _i$ holds the normal component of the right-hand side of (2.3) computed at the collocation point of panel 
$i$. The second and third lines impose the Kelvin theorem and Kutta condition, respectively. Additionally, 
$\varGamma$ is the third component of the right-hand side of (2.4). The solution of (2.7) provides the distribution of circulation on the profile 
$\boldsymbol {\gamma }$ as well as the strength 
$\gamma _s$ of the shed panel. Once the system is solved, we have all the elements at our disposal to compute forces on the body, update its location and convect the wake point vortices.
At the end of a time step, the shed panel is transformed into a point-vortex at its centre and thus has zero length. In this specific case, the last column of the left-hand side matrix is filled with zeros except for its last element. In other words, the strength 
$\gamma _s$ has no contribution over the bound vortex sheet whatsoever and the Kutta condition is trivially satisfied by the existing set of point vortices. Hence, the system is reduced to
\begin{equation} \renewcommand{\arraystretch}{1.5} \underbrace{\left[ \begin{array}{@{}c@{}} \boldsymbol{\mathsf{A}} \\ \hline \boldsymbol{k}^T \end{array}\right]}_{\tilde{\boldsymbol{\mathsf{A}}}} [\,\,\boldsymbol{\gamma}\,] =\left[\begin{array}{@{}c@{}} \boldsymbol{\chi} \\ \hline - \varGamma \end{array}\right]. \end{equation} The present approach differs from the model of Xia & Mohseni (Reference Xia and Mohseni2017) in two ways. First, the quantity 
$\boldsymbol {\gamma }$ in our model balances vorticity associated with non-uniform body motion, which explicitly appears as a surface integral in the right-hand side of (2.3). In spite of the equivalence of both approaches in the end result, this difference results in a dissimilar formulation for the total impulse in the flow and also yields an alternative flow inside of the body volume. The nature of this internal flow is however irrelevant to the simulation purposes. The second minor difference is the use of a fourth-order Runge–Kutta numerical scheme for the time marching procedure in lieu of a forward Euler method.
2.2. Aerodynamic force and torque
Following a straightforward approach, the force applied on the airfoil can be retrieved from the rate of change of total linear impulse in the fluid. One can readily derive a vorticity-based expression for this quantity:
\begin{equation} \boldsymbol{F} =- \frac{\text{d}}{\text{d} t}\left(\int_{\mathcal{S}_b} \boldsymbol{x}\times (\boldsymbol{n}\times\boldsymbol{u})\,\text{d} s + \sum_{k=1}^{N_v} \boldsymbol{x}_k\times\boldsymbol{\varGamma}_k\right), \end{equation}
where 
$\boldsymbol {u} = \boldsymbol {\gamma }\times \boldsymbol {n} + \boldsymbol {u}_b$ consequentially to the no-flow-through condition. The torque can be obtained in a similar way using the rate of change of angular impulse in the inertial frame of reference:
\begin{equation} \boldsymbol{T} =- \frac{1}{2} \frac{\text{d}}{\text{d} t}\left(\int_{\mathcal{S}_b}\boldsymbol{x}\times [\boldsymbol{x} \times (\boldsymbol{n}\times\boldsymbol{u})] \,\text{d} s+\sum_{k=1}^{N_v} \boldsymbol{x}_k\times(\boldsymbol{x}_k\times\boldsymbol{\varGamma}_k) \right) . \end{equation}However, using the budget of total linear (respectively angular) impulse in the flow not only leads to the computation of the force (respectively moment) applied on all the immersed bodies, thereby precluding estimations for individual bodies, but it also leads to numerical complexity and difficulties. Indeed, this formulation involves the numerical integration of moments of vorticity over its entire support. It is readily seen that this will be increasingly expensive as the flow develops and will also make the calculation sensitive to round-off errors.
We address both those shortcomings through a control volume approach involving velocity and vorticity fields only, as proposed by Noca, Shiels & Jeon (Reference Noca, Shiels and Jeon1999). We further extend this formulation to retrieve the hydrodynamic moment on the body as well. For each body, the control volume encloses the infinitely thin attached boundary layer and undergoes a well-defined flux of vorticity with strength 
$\boldsymbol {\gamma }_s$ across its surface at position 
$\boldsymbol {x}_s$ just beyond the trailing edge of the body:
\begin{align} \boldsymbol{F} &=- \frac{\text{d}}{\text{d} t}\int_{\mathcal{S}_b}\boldsymbol{x}\times(\boldsymbol{n}\times\boldsymbol{u})\,\text{d} s \end{align}
\begin{align} &\quad +\int_{\mathcal{S}_b} \left(\frac{1}{2} u^2\boldsymbol{n}- \boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{u}\boldsymbol{u}_b\right)\text{d} s - \boldsymbol{n}_s\boldsymbol{\cdot}(\boldsymbol{u}_s-\boldsymbol{u}_{bs})(\boldsymbol{x}_s\times \boldsymbol{\gamma}_s), \end{align}
\begin{align} \boldsymbol{T} &=- \frac{1}{2} \frac{\text{d}}{\text{d} t}\int_{\mathcal{S}_b}\boldsymbol{x}\times [ \boldsymbol{x} \times (\boldsymbol{n}\times\boldsymbol{u})] \,\text{d} s\nonumber\\ &\quad+ \int_{\mathcal{S}_b} \boldsymbol{x}\times\left(\frac{1}{2} u^2\boldsymbol{n}- \boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{u}\boldsymbol{u}_b\right)\text{d} s - \frac{1}{2}\boldsymbol{n}_s \boldsymbol{\cdot}(\boldsymbol{u}_s-\boldsymbol{u}_{bs})[\boldsymbol{x}_s \times(\boldsymbol{x}_s\times \boldsymbol{\gamma}_s)], \end{align}
where 
$\boldsymbol {n}_s=\boldsymbol {n}(\boldsymbol {x}_s)$, 
$\boldsymbol {u}_s=\boldsymbol {u}(\boldsymbol {x}_s)$ and 
$\boldsymbol {u}_{bs}=\boldsymbol {u}_b(\boldsymbol {x}_s)$.
3. Vortex lumping for general foil geometries
In previous studies, Wang & Eldredge (Reference Wang and Eldredge2013) and Darakananda & Eldredge (Reference Darakananda and Eldredge2019) have developed low-order models for the roll-up of the vortex sheets shed by a moving flat plate. The second of these studies has shown that transferring vorticity from source point 
$\boldsymbol {x}_s$ to target point 
$\boldsymbol {x}_t$ induces a spurious force due to the local violation of Kelvin's conservation of circulation around both points. They proposed to cancel this parasitic force by correcting the velocity of the target particle:
\begin{equation} \Delta\dot{\boldsymbol{x}}_t = \frac{\dot{\varGamma}_t}{\varGamma_t}\left(\frac{\text{d}\boldsymbol{p}(\boldsymbol{x}_t)}{\text{d}\boldsymbol{x}_t}\right)^{-1} [\boldsymbol{p}(\boldsymbol{x}_s) - \boldsymbol{p}(\boldsymbol{x}_t)], \end{equation}
where 
$\boldsymbol {p}(\boldsymbol {\xi })$ is introduced as the impulse due to a point vortex of unit circulation located at 
$\boldsymbol {\xi }$:
\begin{equation} \boldsymbol{p}(\boldsymbol{\xi}) = \boldsymbol{\xi} \times {\boldsymbol{e}}_3 + \int_{\mathcal{S}_b} \boldsymbol{x}(s) \times \hat{\boldsymbol{\gamma}}(s) \,\text{d} s, \end{equation}
$\hat {\boldsymbol {\gamma }}(s)$ being the so-called unit image vortex sheet caused by this point vortex on 
$\mathcal {S}_b$.
The main drawback of this strategy is the impossibility to cancel both the spurious force and moment caused by the displacement of vorticity. Hence, cancelling the spurious force will inevitably result in a spurious moment given by
\begin{equation} \Delta \boldsymbol{T} =-\frac{\rho\dot{\varGamma}_t}{2}\Bigg[ \frac{\text{d} \boldsymbol{m}(\boldsymbol{x}_t)}{\text{d} \boldsymbol{x}_t}\left(\frac{\text{d} \boldsymbol{p}(\boldsymbol{x}_t)}{\text{d} \boldsymbol{x}_t}\right)^{-1} [\boldsymbol{p}(\boldsymbol{x}_s) - \boldsymbol{p}(\boldsymbol{x}_t)] + [\boldsymbol{m}(\boldsymbol{x}_t) - \boldsymbol{m}(\boldsymbol{x}_s)]\Bigg], \end{equation}
where the unit angular impulse due to a vortex located at 
$\boldsymbol {\xi }$ is
\begin{equation} \boldsymbol{m}(\boldsymbol{\xi}) = \boldsymbol{\xi} \times (\boldsymbol{\xi} \times {\boldsymbol{e}}_3) + \int_{\mathcal{S}_b} \boldsymbol{x}(s) \times (\boldsymbol{x}(s) \times \hat{\boldsymbol{\gamma}}(s) ) \,\text{d} s. \end{equation}Darakananda & Eldredge (Reference Darakananda and Eldredge2019) applied the thin airfoil theory to derive the analytical velocity correction in the case of the flat plate. However, in the panel-based model, we propose that such an analytical derivation is impossible because the Jacobian in (3.1) must be obtained numerically. In two dimensions, the Jacobian of the unit impulse is the following:
\begin{equation} \frac{{\rm d}\boldsymbol{p}(\boldsymbol{\xi})}{{\rm d}\boldsymbol{\xi}} = \begin{bmatrix} 0 & 1\\ -1 & 0 \end{bmatrix}+\int_{\mathcal{S}_b} \begin{bmatrix} y(s)\\ -x(s) \end{bmatrix} \nabla_{\boldsymbol{\xi}}\hat{\boldsymbol{\gamma}}(s)\,\text{d} s. \end{equation}
The only missing piece in this expression is the gradient of the bound sheet circulation 
$\hat {\boldsymbol {\gamma }}(s)$ on a contour 
$\mathcal {S}_b$ of general shape. It represents the sensitivity of the bound vortex sheet to a change in position of a free point vortex of unit circulation.
Because the lumping operation is performed at the end of a time step (i.e. when the shed panel is replaced by a point vortex), we can take the gradient of (2.8) to obtain 
$\nabla _{\boldsymbol {\xi }}\hat {\boldsymbol {\gamma }}$. Since the total circulation is conserved throughout the merging of vortices and because matrix 
$\tilde {\boldsymbol{\mathsf{A}}}$ does not depend on 
$\boldsymbol {\xi }$ but solely on the relative positions of the bound panels, the gradient on the body surface reduces to
\begin{equation} \renewcommand{\arraystretch}{1.5} [\nabla_{\boldsymbol{\xi}}\hat{\boldsymbol{\gamma}}] = \tilde{\boldsymbol{\mathsf{A}}}^{-1} \Bigg[\begin{array}{@{}c@{}} \nabla_{\boldsymbol{\xi}}\hat{\boldsymbol{\chi}} \\ \hline 0 \end{array}\Bigg], \end{equation}
where 
$\tilde {\boldsymbol{\mathsf{A}}}^{-1}$ has already been computed in the no-flow-through condition determining the bound vortex strength and where the induced normal velocity gradient 
$\nabla _{\boldsymbol {\xi }}\hat {\boldsymbol {\chi }}$ has to be obtained at the collocation point of each panel. For panel 
$j$, we have
\begin{equation} \nabla_{\boldsymbol{\xi}}\hat{\chi}_{j} = \nabla_{\boldsymbol{\xi}}(\boldsymbol{n}_j \boldsymbol{\cdot} (\boldsymbol{K}^*(\boldsymbol{x}_j - \boldsymbol{\xi}) \times {\boldsymbol{e}}_3 )), \end{equation}
which is obtained analytically in the frame of reference of panel 
$j$.
First of all, the lumping operation is only allowed between vortex elements of the same sign. Furthermore, even though the merging of point vortices itself is carried out in a way to cancel spurious forces, the main structure of the flow is altered in the vicinity of the transfer of vorticity. This change in flow topology results in a different solution for the motion of vortex elements computed over a time step. This in turn causes both the global impulse and the force on the airfoil to drift away from the values which would result from an unperturbed vorticity field. To control this discrepancy and in line with the work of Darakananda & Eldredge (Reference Darakananda and Eldredge2019), we use a parameter 
$B_F$ that serves as a threshold on the force discrepancy observed one time step after the lumping operation has occurred. This discrepancy is readily evaluated by integrating both the original vortex sheet and the lumping-affected vortex sheet over a time step. The choice in 
$B_F$ balances the representational accuracy of the flow with the computational efficiency of the method. Indeed, low values of 
$B_F$ will result in more resolved vortex sheets while large 
$B_F$ values are associated with coarser wake structures. One can readily identify a lower bound for 
$N_v$: the number of lumped vortices can go as low as the number of sign changes in the shed vorticity plus an initial vortex. Once this bound is reached at a particular 
$B_F$ level, further increasing the threshold no longer has an effect on the structure of the wake. It follows that the behaviour of 
$N_v$ with respect to 
$B_F$ will be flow specific.
This model embeds two other parameters: the minimal vortex sheet length, described by the minimum number of shed point vortices 
$N_{min}$ before which no lumping can occur whatsoever, and the minimum time interval 
$T_{min}$, within which no new vortex can be shed from the tip of the vortex sheet. Parameter 
$N_{min}$ allows to keep unaltered the coherent structure of the velocity close to the airfoil while 
$T_{min}$ aims at avoiding instabilities in the vicinity of the vortex sheet tip.
4. Experiments and validation
The model presented above is now validated on: (i) an impulsively started NACA0012 and (ii) a heaving and pitching NACA0013. The results obtained are then compared with the simulations of Xia & Mohseni (Reference Xia and Mohseni2017) and the experimental results from Izraelevitz & Triantafyllou (Reference Izraelevitz and Triantafyllou2014). In all simulations, the airfoil of chord 
$c$, towed at horizontal velocity 
$U_\infty$, is uniformly discretized with 200 panels of linear vortex strength and a time step of 
$10^{-2} \times c/U_\infty$ is used. The blob radius 
$\delta$ is set to 
$10^{-2} c$.
4.1. Impulsively started NACA0012
We first compare the influence of the parameter 
$B_F$ of the lumping model on the aerodynamic loads of the impulsively started NACA0012 at a 
$10^\circ$ angle of attack (AoA). Snapshots of simulations at two and five chord lengths of travel are presented in figure 2. The first case in which 
$B_F = 0$ corresponds to the standard unsteady panel method, later referred to as the baseline model. Increasing this parameter yields a dramatic reduction in the number of point vortices in the wake. At 
$B_F=10^{-3}$, the wake is only captured by means of 
$L_{min}+3$ point vortices. Increasing this parameter beyond 
$B_F=10^{-2}$ results in the formation of a unique point vortex corresponding to the starting vortex of the airfoil.
Figure 2. Comparison of the vorticity distribution predicted for an impulsively started NACA0012 airfoil at 
$10^\circ$ after two and five chord lengths of travel between the proposed model with different values of the error threshold 
$B_F$. The other parameters are fixed to 
$T_{min} = L_{min} = 25$.
Figure 3 shows the aerodynamic coefficients up to ten chord lengths of travel computed with these different values of 
$B_F$. The aerodynamic coefficients are defined as 
$C_D = -2F_{\boldsymbol {e}_1}(\rho U^2_\infty c)^{-1}$, 
$C_L = 2F_{\boldsymbol {e}_2}(\rho U^2_\infty c)^{-1}$ and 
$C_{M,1/4} = 2T_{\boldsymbol {e}_3}(\rho U^2_\infty c^2)^{-1}$. Because the lumping operation affects the angular momentum budget in the computational domain, 
$\boldsymbol {F}$ and 
$\boldsymbol {T}$ are computed using (2.11) and (2.13) following the control volume approach. When parameter 
$B_F$ is equal to or higher than 
$10^{-2}$, the error in drag coefficient is smaller than 
$10\,\%$. The model predicts a lift coefficient that is higher by less than 
$2\,\%$ of its final value while the moment coefficient is not accurately predicted at the early stage of the simulation.
Figure 3. Impulsively started NACA0012 with 
$\alpha =10^\circ$. (a) Drag, (b) lift and (c) moment coefficients obtained with the baseline model (blue, —), and the lumping model with 
$T_{min} = L_{min} = 25$, 
$B_F = 10^{-3}$ (green, -
$\cdot$-
$\cdot$) and 
$B_F = 10^{-2}$ (orange, - - -).
The small loss of accuracy in lift and drag coefficients comes with a significant speedup in computational efficiency. Since the time complexity per time step scales quadratically with the number of vortex particles, limiting their number to 
$L_{min}+1$ (instead of a linear increase up to 1000 at the end of the baseline simulation) keeps the cost constant throughout the whole simulation. The total speedup in this case is close to 3.5.
Finally, results of computations at various AoAs with 
$B_F = 10^{-2}$ are compared with the simulations of Xia & Mohseni (Reference Xia and Mohseni2017) in figure 4. For any configuration, the lift coefficient agrees well with their model with a unique point vortex shed at the trailing edge. As expected, the shedding angle converges towards the bisector of the trailing edge. However, this occurs at a slower rate than what Xia & Mohseni (Reference Xia and Mohseni2017) observed.
Figure 4. Impulsively started NACA0012 airfoil at 
$\alpha =2^\circ$ (green), 
$\alpha =6^\circ$ (orange) and 
$\alpha =10^\circ$ (blue). Comparison of (a) the lift coefficient and (b) the shedding angle between the lumping model with 
$T_{min} = L_{min} = 25$, 
$B_F=10^{-2}$ (—) and computational results of Xia & Mohseni (Reference Xia and Mohseni2017) (- - -) along with steady values obtained with a potential flow solver (
$\cdot \cdot \cdot \cdot \cdot \,\cdot$). Shedding angle 
$\theta _s$ is defined with respect to the bisector of the trailing edge; 
$\theta _s \in [-\theta _{TE}/2, \theta _{TE}/2]$.
4.2. Heaving and pitching NACA0013
In this section, we consider a NACA0013 airfoil evolving with constant horizontal velocity 
$U_\infty$. Its trajectory is characterized by a sinusoidal heaving motion 
$y(t)$ perpendicular to 
$U_\infty$ given by
\begin{equation} y(t) = hc\cos(\omega t), \end{equation}
and an angular motion about its quarter chord 
$\theta (t)$. The pitch motion is such that the relative AoA is the sinusoid:
\begin{equation} \alpha(t) = \theta(t) - \arctan\left(\frac{\dot{y}(t)}{U_\infty}\right) = \alpha_{max} \sin(\omega t). \end{equation}Finally, the dimensionless number linking the frequency and amplitude of the motion is the Strouhal number defined as
\begin{equation} \textit{St}= \frac{\omega hc}{{\rm \pi} U_\infty}. \end{equation}
The motion can now be fully described with the three fixed parameters: 
$h = 1$, 
$\alpha _{\text {max}} = 25^\circ$ and 
$\textit {St} = 0.3$. In this situation, Izraelevitz & Triantafyllou (Reference Izraelevitz and Triantafyllou2014) have experimentally shown the absence of leading edge separation at Reynolds number of 11 000. Therefore, this choice of motion parameters suits the proposed model.
First, we examine the effect of applying the lumping model on vorticity distribution in the wake of the airfoil. Figure 5 shows snapshots of simulations where the error parameter 
$B_F$ is set to either 
$0$, 
$10^{-3}$ or 
$10^{-1}$, captured after one and a half motion periods. Past the highest tolerance level, the wake can be fully reduced to a pair of vortices of equal and opposite strength for each motion period. Their arrangement fits the characteristic structure of a thrust wake (see Oskouei & Kanso Reference Oskouei and Kanso2013). As a point of reference, the original model would generate 
$666$ vortices per oscillation period with the same set of simulation parameters (and 27 vortices per period with 
$B_F=10^{-3}$). This allows a maximal reduction in the number of vortices by a factor larger than 300 compared to the baseline model, at the cost of slightly affecting the aerodynamic coefficients. This is shown in figure 6 where the lightest model is compared with the baseline, experimental results of Izraelevitz & Triantafyllou (Reference Izraelevitz and Triantafyllou2014) and computations of Xia & Mohseni (Reference Xia and Mohseni2017). The thrust coefficient 
$C_T$ is defined here as the opposite of the drag coefficient 
$C_D$. The error due to the most severe lumping procedure does not exceed 
$10\,\%$ of the amplitude of the curves obtained with the baseline. Hence, independently of the choice of the force error threshold 
$B_F$ and given the specific motion under examination, our lumping-enabled model is able to capture accurately the dynamics of the airfoil and produce results complying with the literature (see figure 5).
Figure 5. Heaving and pitching NACA0013 airfoil at 
$\textit {St}=0.3$, 
$h=c$, 
$\alpha _{max}=25^\circ$, 
$T_{min} = L_{min} = 25$ and (a) 
$B_F=0$, (b) 
$B_F=10^{-3}$ and (c) 
$B_F=10^{-1}$. Comparison of vorticity field.
Figure 6. Heaving and pitching NACA0013 airfoil with 
$\textit {St}=0.3$, 
$h=c$ and 
$\alpha _{max}=25^\circ$. Comparison of (a) thrust, (b) lift and (c) moment coefficients between the model with 
$T_{min} = L_{min} = 25$, 
$B_F = 10^{-1}$ (blue, —), the baseline model (orange, - - -), computational results of Xia & Mohseni (Reference Xia and Mohseni2017) (
$\cdot \cdot \cdot \cdot \cdot \,\cdot$) and the experimental results of Izraelevitz & Triantafyllou (Reference Izraelevitz and Triantafyllou2014) (
$\triangle$).
5. Conclusions
In this paper, we bring two contributions to the modelling of unsteady airfoil flows. As a first contribution, we have developed and validated the generalization of the hybrid model of Darakananda & Eldredge (Reference Darakananda and Eldredge2019) for modelling flows separating from the trailing edge of general airfoils in unsteady motions. The model enables a user-controlled reduction of computational complexity by lumping circulation from one vortex element to another without introducing any instantaneous force imbalance on the airfoil. This transfer of circulation could be applied to richer models embedding multiple bodies or/and LEV separation. To achieve the model extension, we have determined the net effect of the shape of the airfoil on the impulse of a point vortex. This contribution is computed through a recycling of the inverse of the matrix used in the standard unsteady panel method to determine the vortex strength distribution on the profile. This recycling allows a fast computation of the velocity correction the target vortex has to undergo to avoid spurious forces on the body. The remaining computational effort lies in calculation of the force error caused by the transfer of circulation between vortices over a time step.
In addition, we proposed the application and extension of a control volume approach relying on velocity and its derivatives only (see Noca et al. Reference Noca, Shiels and Jeon1999) to compute aerodynamic force and moment undergone by a body with an explicit treatment of the singularity in the wall-vorticity flux in an inviscid flow. To the best of our knowledge, this is the first time this method is involved in the computation of the aerodynamic moment. This compact formulation also opens the field to low-order simulations of multiple immersed bodies and avoids the cost and sensitivity associated with largely populated far-wake vortex sheets.
We have validated this model on two benchmarks: an impulsively started NACA0012 and a heaving and pitching NACA0013 airfoil. In both cases, we showed the effect of lumping vortices on the vorticity distribution in the wake of the airfoils. The main structure of the wake was conserved even though far fewer vortex elements were required to describe it. In the limit of acceptable accuracy, one single point vortex per roll-up core of the wake vortex sheet is necessary. Such an important dimensionality reduction accompanies a significant speedup allowing for real-time simulations. In addition, aerodynamic forces corroborate results obtained with the unsteady panel method of Xia & Mohseni (Reference Xia and Mohseni2017) as well as experimental results of Izraelevitz & Triantafyllou (Reference Izraelevitz and Triantafyllou2014). The lumping operation has shown to only degrade the accuracy in the moment coefficient at early stages of a simulation if a very large transfer of circulation occurs between the body and the starting vortex close to the body surface. This significant degradation was not observed in the case of the heaving and pitching airfoil because the airfoil is initially parallel to its towing motion and the generation of circulation in the flow is much more progressive than in the case of the impulsively started airfoil at 
$10^{\circ }$ of AoA. Therefore, in cases of unsteady motion of an airfoil shedding vorticity continuously from its trailing edge and without leading edge separation, our model is fully capable of capturing essential dynamics of the flow and predict the aerodynamic loads on the immersed body while offering the user the freedom to adjust the balance between accuracy and computational time.
Acknowledgements
Thomas Gillis and Denis-Gabriel Caprace deserve special thanks. Their guidance was remarkable and they contributed greatly, not only to D.D.'s mathematical development but also to building a motivating working environment.
Funding
This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement no. 725627).
Declaration of interests
The authors report no conflict of interest.
Open access
