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:

1. (i) a collection of $N_p$ panels with linear strength located on $\mathcal {S}_b$ to model the infinitely thin boundary layer;

2. (ii) a unique panel in the very near wake $\mathcal {S}_s$ with length $b_s$ and uniform strength $\gamma _s$;

3. (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)\begin{equation} \boldsymbol{K}^*(\boldsymbol{x}) =- \frac{\boldsymbol{x}}{2{\rm \pi} (|\boldsymbol{x}|^2 + \delta^2)}, \end{equation}
   where $\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

(2.2)\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))

(2.3)\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:

1. (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);

2. (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}

3. (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)\begin{equation} \gamma_s = \gamma_1 \cos{(\theta_+)} + \gamma_{N_p+1} \cos{(\theta_{TE} - \theta_+)}, \end{equation}
   where $\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:

(2.7) \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

(2.8)\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:

(2.9)\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:

(2.10)\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:

(2.11)\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}

(2.12)\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}

(2.13)\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)$.

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

(4.1)\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:

(4.2)\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

(4.3)\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$).