NOMENCLATURE

MAV

Micro-Air-Vehicle

LEV

Leading Edge Vortex

UVLM

Unsteady Vortex Lattice Method

CFD

Computational Fluid Dynamics

LE,TE

Leading edge, trailing edge

LEV

Leading edge vortex

SS4

Model skill score

QS

Quasi-Steady

Symbol

$\mathcal{O}(\cdot)$

Order of magnitude

$\mathcal{E}$ , $\partial\mathcal{E}$ , $\mathcal{C}$ , $\partial\mathcal{C}$

Set of points within and at the boundary of an elliptical and circular cross-sectional shaped body, respectively

${\boldsymbol{\omega}}$ , $\omega$

Angular velocity and its magnitude, respectively

$\hat{\textbf{\textit{n}}}_1, \hat{\textbf{\textit{n}}}_2, \hat{\textbf{\textit{n}}}_3$

Orthonormal set of unit vectors defining an inertial reference frame

$\hat{\boldsymbol{e}}_1, \hat{\boldsymbol{e}}_2, \hat{\boldsymbol{e}}_3$

Orthonormal set of unit vectors defining a body-fixed reference frame

O

Pivot point of the wing section

t

Time

${\textbf{\textit{r}}}$

Position vector

${\textbf{\textit{v}}}({\textbf{\textit{r}}}_p,t)$

Velocity of an arbitrary point P within the body

${\textbf{\textit{u}}}({\textbf{\textit{r}}},t)$

Fluid velocity at location ${\textbf{\textit{r}}}$

${\textbf{\textit{r}}}^{P/O}$ , ${\textbf{\textit{r}}}^{C/O}$

Relative position of points P and C, respectively, with respect to point O

${\textbf{\textit{r}}}^{O/C}$ , $z_O$

Relative position of point O with respect to point C, and its complex counterpart, respectively

${\textbf{\textit{v}}}_O$ , $V_O$

Velocity of point O and its magnitude, respectively

${\textbf{\textit{v}}}_C$ , $\mathscr{V}_C$

Velocity of point C and its complex counterpart, respectively

${\textbf{\textit{u}}}_\infty$ , $U_\infty$

Velocity of the free flow and its magnitude, respectively

$\varphi$ , $\vartheta$

Orientation of the ${\textbf{\textit{v}}}_C$ and ${\textbf{\textit{u}}}_\infty$ with respect to $\hat{\textbf{\textit{b}}}_1$ , respectively

$\hat{\textbf{\textit{s}}}_\varphi$ , $\hat{\textbf{\textit{s}}}_\vartheta$

Unit vectors defining the direction of ${\textbf{\textit{v}}}_C$ and ${\textbf{\textit{u}}}_\infty$ , respectively

$\Phi$

Velocity Potential

$\hat{\textbf{\textit{n}}}({\textbf{\textit{r}}})$

Unit vector normal to the surface of the body $\mathcal{E}$ at location ${\textbf{\textit{r}}}$

$p({\textbf{\textit{r}}})$ , $p_\infty$

Pressure at location ${\textbf{\textit{r}}}$ and at infinity, respectively

$\rho$

Fluid density

$g(\cdot)$

Conformal map

$\lambda$

Conformal map parameter

j

Imaginary unit

$\zeta, z$

Complex coordinates in the circle and physical plane, respectively, defined in the body-fixed frame

$\mathscr{F}$

Complex potential

$\Gamma$ , $\Gamma_j$

Point vortex intensity

$\mu_j$ , $\beta_j$

Doublet intensity and orientation, respectively

R

Circle radius

$\mathcal{F}$

Complex force

${\textbf{\textit{F}}}$ , ${\textbf{\textit{F}}}_{\text{2D}}$ , ${\textbf{\textit{F}}}_{\text{3D}}$

Vector force per se, two-dimensional, and three-dimensional, respectively

${\textbf{\textit{F}}}_{\text{LE}}$ , ${\textbf{\textit{F}}}_{\text{LEV}}$

Leading edge force, Leading edge vortex force

$\varrho$

Leading edge vortex force parameter

$\gamma_\mathcal{E}$ , $\gamma_\mathcal{W}$

Vortex sheet intensity distribution

$A_{ij}$

Influence coefficient of vortex j over control point i

$C_{d\_{qs}}$ , $C_{l\_{qs}}$

Quasi-steady translation drag and lift coefficient

$\alpha$

Flat plate angle-of-attack

$t_vV/c$

Nondimensional time for computation of quasi-steady translation loads

$C_d$ , $C_l$

Normalised drag and lift coefficient

$A_f$ , $\Theta$

Amplitude of the translational and rotational motion, respectively

b, $\mathcal{S}$

Wing span and planform area of a three-dimensional planar wing

$r_{11}$

Radius of gyration abut the axis along the $\hat{\textbf{\textit{e}}}_1$ direction

2.0 MATERIALS, METHODS AND THEORY

Insects and small birds fly at low-to-moderate chord based Reynolds numbers: $\mathcal{O}(10^1)$ for small insects such as Chalcid wasps, $\mathcal{O}(10^2)$ for fruit flies, $\mathcal{O}(10^3)$ for bees and dragonflies, and $\mathcal{O}(10^4)$ for moths and hummingbirds (Table 1.1, p. 15 of ref. (Reference Altenbuchner and Hubbard30)). Experimental evidence^{(Reference Dickinson, Lehmann and Sane7,Reference Wang, Birch and Dickinson24)} has shown that the aerodynamic forces generated by flapping wings are dominated by pressure and that the viscous effects play a secondary role, even for Reynolds numbers as low as 75. Hence, the observed aerodynamic forces are primarily normal to the solid surfaces. These observations support the hypothesis that the main features of the flow and forces generated by a flapping wing can be described to certain extent through inviscid fluid models.

The system considered is a two-dimensional wing section undergoing flapping-like motions (translation and rotation) within still or moving air. In general, the wing section has an elliptical shape, which can degenerate to flat plate. Only the flat plate section is studied in this article. The resulting airflow is assumed to be incompressible and irrotational except at the interface (vortex sheet) between the solid surfaces and the air where infinitesimally thin boundary layers are present.

Figure 1.

Two-dimensional elliptic body $\mathcal{E}$ undergoing prescribed translation and rotation. (a) Conformal map of a circle into an ellipse. (b) ULVM discretization.

The motion of a body $\mathcal{E}$ in the x-y plane is described by an angular velocity $\boldsymbol{\omega}=\omega\hat{\textbf{\textit{n}}}_3$ about a pivot point O and a translation velocity ${{\textbf{\textit{v}}}_O=V_O \hat{\textbf{\textit{s}}}_\varphi}$ of point O, with $\hat{\textbf{\textit{s}}}_\varphi$ a unit vector with direction given by the angle $\varphi$ (Fig. 1(a)). The fluid (air) may have a uniform velocity ${\textbf{\textit{u}}}_\infty =U_\infty \hat{\textbf{\textit{s}}}_\vartheta$ , with $\hat{\textbf{\textit{s}}}_\vartheta$ a unit vector with direction given by the angle $\vartheta$ . The velocity of an arbitrary point $P \in\mathcal{E}$ of the body is given by:

(1) \begin{align}{\textbf{\textit{v}}}({\textbf{\textit{r}}}_P,t)= {\textbf{\textit{v}}}_O+ \boldsymbol{\omega} \times{\textbf{\textit{r}}}^{P/O},\end{align}

where ${\textbf{\textit{r}}}^{P/O}$ is the relative position of point P with respect to O. Of particular interest is the velocity of the centroid of the section C, which according to (1) is given by:

(2) \begin{align}{\textbf{\textit{v}}}_C(t):={\textbf{\textit{v}}}({\textbf{\textit{r}}}_C,t)= {\textbf{\textit{v}}}_O+ \boldsymbol{\omega} \times{\textbf{\textit{r}}}^{C/O}={\textbf{\textit{v}}}_O- \boldsymbol{\omega} \times{\textbf{\textit{r}}}^{O/C}.\end{align}

Assuming that the flow is irrotational outside the vortex sheets, the continuity equation for irrotational and incompressible flows becomes the Laplace Equation (3) for the velocity potential $\Phi$ :

(3) \begin{align}\nabla^2 \Phi =0.\end{align}

The flow field ${\textbf{\textit{u}}}({\textbf{\textit{r}}},t)=\nabla \Phi$ must satisfy the no-penetration condition on the surface of the body $\partial \mathcal{E}$ :

(4) \begin{align}\left({\textbf{\textit{u}}}({\textbf{\textit{r}}},t) - {\textbf{\textit{v}}}({\textbf{\textit{r}}},t)\right)\cdot\hat{\textbf{\textit{n}}}({\textbf{\textit{r}}}) = 0&& \forall\ t\in[0,\infty);\ {\textbf{\textit{r}}}\in\partial \mathcal{E}.\end{align}

Here $\hat{\textbf{\textit{n}}}({\textbf{\textit{r}}})$ is the unit vector normal to the surface of the body at point ${\textbf{\textit{r}}}\in\partial\mathcal{E}$ . Additionally, it is required that the perturbations introduced by the body decay to zero far away from it. This regularity condition at infinity is satisfied by the fundamental solutions (vortex, source/sink, dipole) of the Laplace Equation (3). For incompressible, irrotational flows, the unsteady Bernoulli’s equation is employed to obtain the pressure distribution over the body surface once the velocity potential $\Phi$ is known. This equation is given by

(5) \begin{align}p({\textbf{\textit{r}}},t)-p_\infty=&\frac{\rho}{2}\left({\textbf{\textit{u}}}_\infty\cdot{\textbf{\textit{u}}}_\infty-{\textbf{\textit{u}}}({\textbf{\textit{r}}},t)\cdot{\textbf{\textit{u}}}({\textbf{\textit{r}}},t)\right)-\frac{\partial}{\partial t}\left(\Phi({\textbf{\textit{r}}},t)-\Phi_\infty\right)\ \ \text{for}\ {\textbf{\textit{r}}}\in \partial\mathcal{E}.\end{align}

2.1 Analytical potential flow formulation

The circulatory flow around a moving body $\mathcal{E}$ is obtained through conformal mapping of the equivalent flow around a moving circle $\mathcal{C}$ (Fig. 1(a)). A Joukouwsky transformation ${g:\mathbb{C}\ni \zeta\mapsto z\in \mathbb{C}}$ ( ${z=x+jy}$ , ${\zeta=\xi+j\eta}$ ) is employed to transform the fluid domain outside a circle $\mathcal{C}$ of radius R (defined in the $\zeta$ -complex plane) into the fluid domain outside an elliptical body $\mathcal{E}$ (defined in the z-complex plane). This conformal map is given by:

(6) \begin{align}g(\zeta):=a\left(\zeta+\frac{\lambda^2}{\zeta}\right).\end{align}

The coordinates x-y are referred to a reference frame $\left\{\hat{\textbf{\textit{e}}}_1,\hat{\textbf{\textit{e}}}_2\right\}$ which is attached to the centroid C of the body $\mathcal{E}$ . Similarly, the coordinate system $\zeta \hbox{-} \eta$ is referred to the frame $\left\{\hat{\textbf{\textit{e}}}_1,\hat{\textbf{\textit{e}}}_2\right\}$ and its origin is located at the centroid C of the circle $\mathcal{C}$ .

It can be shown that the complex potential $\mathscr{F}$ of the equivalent flow around the circle $\mathcal{C}$ is given by:

(7) \begin{align}\mathscr{F}(\zeta,t){}={}& U_\infty(t) g(\zeta) e^{-j\vartheta(t)}-j\frac{\Gamma(t)}{2\pi} \ln\left(\frac{\zeta}{R}\right)-\frac{\mu_1(t)}{2\pi}\frac{e^{\,j\beta_1(t)}}{\zeta}-\frac{\mu_2(t)}{2\pi}\frac{e^{\,j\beta_2(t)}}{\zeta}\nonumber\\ &+\frac{1}{\zeta}a\left(\lambda^2\bar{\mathscr{V}}_C(t)-R^2\mathscr{V}_C(t)\right)-j\omega(t) \frac{\lambda^2a^2R^2}{\zeta^2},\end{align}

where the point vortex circulation $\Gamma$ is taken positive in the counter-clockwise sense, the overline indicates complex conjugation, $\mu_1e^{\,j\beta_1}$ and $\mu_2e^{\,j\beta_2}$ are two dipoles defined in the $\zeta$ -plane and centered at C, $\mathscr{V}_C\in\mathbb{C}$ is the complex counterpart of ${\textbf{\textit{v}}}_C$ (Equation (2)) such that ${\textbf{\textit{v}}}_C=\mathscr{V}_C\hat{\textbf{\textit{e}}}_1$ Footnote ¹. If $\mathscr{V}_O=V_Oe^{\,j\varphi}\in\mathbb{C}$ and $z_O\in \mathbb{C}$ are the complex counter parts of ${\textbf{\textit{v}}}_O$ ( $=\mathscr{V}_O\hat{\textbf{\textit{e}}}_1$ ) and ${\textbf{\textit{r}}}^{O/C}$ ( $=z_O\hat{\textbf{\textit{e}}}_1$ ), respectively, then it can be shown that

(8) \begin{align}{\textbf{\textit{v}}}_C=\mathscr{V}_C\hat{\textbf{\textit{e}}}_1=\left(\mathscr{V}_O-j\omega z_O\right)\hat{\textbf{\textit{e}}}_1.\end{align}

The last two terms in Equation (7), related to the body kinematics, arise from the boundary function introduced by Milne-Thompson^{(Reference Milne-Thompson31)} (Section 9.63). It can be shown that the no-penetration Equation (4) on the circle yields:

(9) \begin{align}\mu_1e^{\,j\beta_1}&= -2\pi aR^2 U_\infty e^{\,j\vartheta}\nonumber\\\mu_2e^{\,j\beta_2}&=2\pi a\lambda^2 U_\infty e^{-j\vartheta}.\end{align}

2.1.1 Fluid force

The resultant (complex) force $\mathcal{F}$ is obtained through the unsteady version of Blasius’ formula^{(Reference Milne-Thompson31)}. Following Section 6.41 of ref. (Reference Milne-Thompson31), the complex conjugate of the force is given by

(10) \begin{align}\overline{\mathcal{F}}&=-j\rho \oint_{\partial\mathcal{E}} P(z)d\bar{z}=-j\rho \oint_{\partial\mathcal{C}} {\mathcal{P}}(\zeta)\frac{d\bar g}{d\bar{\zeta}}d\bar{\zeta},\end{align}

where $\mathcal{P}(\zeta)=P(g(\zeta))$ is the pressure distribution on the surface of the circle $\partial\mathcal{C}$ . Upon substituting Equations (7), (9) and (5) (with vector operations replaced by the equivalent complex number products) into (10) and carrying out the integration, the following result is obtained:

(11) \begin{align}\overline{\mathcal{F}}{}={}&\rho j \left(\Gamma-2\pi a^2\omega \frac{\lambda^4}{R^2}-\omega \mathcal{A}_\mathcal{E}\right)\left( U_\infty e^{-j\vartheta}-\bar{\mathscr{V}}_C \right)+2\pi j \omega \rho a^2\lambda^2\left(U_\infty e^{\,j\vartheta}-\mathscr{V}_C\right)\nonumber\\&-2\pi \rho a^2\lambda^2 \frac{d}{dt}\left(U_\infty e^{\,j\vartheta}-\mathscr{V}_C\right)+j\rho \bar{g}(R)\frac{d\Gamma}{dt}\nonumber\\&+\rho \pi a^2\left(\frac{\lambda^4}{R^2}+R^2\right)\frac{d}{dt}\left(U_\infty e^{-j\vartheta}-\bar{\mathscr{V}}_C\right).\end{align}

where $\mathcal{A}_\mathcal{E}$ is the cross-section of the elliptical body. This expression of the force is valid for any elliptical two-dimensional body $\mathcal{E}$ with arbitrary circulation $\Gamma$ .

2.1.2 Flat plate with smooth flow off of the trailing edge

Choosing $\lambda=R$ in (6) degenerates the elliptical body $\mathcal{E}$ into a flat plate ( $\mathcal{A}_\mathcal{E}=0$ ) with chord ${c=4aR}$ . It can be shown that the quasi-steady Kutta condition requiring zero fluid velocity relative to the solid at trailing edge point is given by:

(12) \begin{equation}\Gamma_K =\pi c \left(V_O \sin \varphi- U_{\infty} \sin\vartheta\right) +c^2 \pi \omega\left(\frac{1}{4}- \frac{x_O}{c}\right).\end{equation}

Here, $V_O$ is the magnitude of ${\textbf{\textit{v}}}_O$ and $\Re\left[z_O\right]=x_O\in\mathbb{R}$ is the position of the pivot point with respect to the centroid C. It is noted that $\Gamma_K$ is in general a function of time. However, the quasi-steady Kutta condition also implies that the rate of generation of circulation from the trailing edge is zero $d\Gamma_K/dt=0$ (Equation (27)). This contradiction is eliminated if a wake model is introduced, as it is done in the UVLM approach. Following Minotti^{(Reference Minotti17)}, $d\Gamma/dt=0$ is adopted in the analytical models for all the cases analyzed in this article.

From Equation (11) the resultant (vector) fluid force on the flat plate is given by:

(13) \begin{align}{\textbf{\textit{F}}}{}={}&\rho \Gamma\left( U_\infty \sin\vartheta-V_O \sin\varphi\right)\hat{\textbf{\textit{t}}}-\rho \Gamma\left( U_\infty \cos\vartheta-V_O \cos\varphi\right)\hat{\textbf{\textit{n}}}\nonumber\\&-\pi \rho \omega \frac{c^2}{4}\left(U_\infty \sin\vartheta-V_O \sin \varphi\right)\hat{\textbf{\textit{t}}}+\rho \Gamma \omega x_O \hat{\textbf{\textit{t}}} - \frac{1}{4}\pi \rho \omega c^2x_O\hat{\textbf{\textit{t}}}\nonumber\\&-\rho \frac{c}{2}\frac{d\Gamma}{dt}\hat{\textbf{\textit{n}}}+\pi \rho \frac{c^2}{4}x_O\frac{d\omega}{dt}\hat{\textbf{\textit{n}}}+\pi \rho \frac{c^2}{4} \frac{d}{dt}\left(U_\infty \sin \vartheta-V_O \sin\varphi\right)\hat{\textbf{\textit{n}}},\end{align}

where $\hat{\textbf{\textit{t}}}$ is the tangent unit vector to the plate directed from C towards the trailing edge and $\hat{\textbf{\textit{n}}}$ the normal unit vector such that $\hat{\textbf{\textit{t}}}\times\hat{\textbf{\textit{n}}}=\hat{\textbf{\textit{n}}}_1\times\hat{\textbf{\textit{n}}}_2=\hat{\textbf{\textit{n}}}_3$ (Fig. 1(b)).

Analysing the flow around the leading and trailing edges (LE and TE, respectively), the following forces develop to allow the flow to turn around the sharp edges

(14) \begin{equation}\hspace*{-2pt}{\textbf{\textit{F}}}_{\text{LE}}=-\frac{\rho \pi c}{4}\left(\frac{\Gamma}{\pi c}-c\omega\left(\frac{1}{4}+\frac{x_O}{c}\right)+V_O\sin\varphi -U_{\infty}\sin\vartheta \right)^2\hat{\textbf{\textit{t}}}\end{equation}

(15) \begin{align}{\textbf{\textit{F}}}_{\text{TE}}=\frac{\rho \pi c}{4}\left(-\frac{\Gamma}{\pi c}+c\omega\left(\frac{1}{4}-\frac{x_O}{c}\right)+V_O\sin\varphi-U_{\infty}\sin\vartheta \right)^2\hat{\textbf{\textit{t}}}.\end{align}

It is observed that ${\textbf{\textit{F}}}_{\text{TE}}\equiv{\textbf{0}}$ if $\Gamma \equiv \Gamma_K $ ; that is to say, the suction force on the trailing edge vanishes when the flow leaves the edge smoothly and the Kutta condition is satisfied. In that situation, however, the flow still turns around the leading edge and the (suction) force there becomes

(16) \begin{equation}\left.{\textbf{\textit{F}}}_{\text{LE}}\right|_{\Gamma=\Gamma_K}=-\rho\pi c\left(V_O\sin\varphi -U_{\infty}\sin\vartheta-\omega x_O\right)^2\hat{\textbf{\textit{t}}}.\end{equation}

2.1.3 Modified force model for leading edge vortex effect

A modification to the aerodynamic force predicted by the potential flow model (13) on a flat plate is proposed. By analogy to the situation at the trailing edge, it is argued that if the flow were to separate from the leading edge, the suction force there (Equation (16)) vanishes. With this mind, a two-step modification of the theoretical aerodynamic force (Equation (13)) is proposed to account for an “attached” leading edge vortex. The first step is to suppress the leading edge force (16). The second step is to add a force ${\textbf{\textit{F}}}_{\text{LEV}}$ normal to the leeward side of the section with magnitude proportional to $F_{\text{LE}}$ . This normal force is thought to prevent the leading edge vortex from detaching. This reasoning is consistent with the the well-known leading edge suction analogy proposed by Polhamus^{(Reference Polhamus32)} to compute the forces on delta wings.

The proposed resultant modified fluid force on the flat plate is given by:

(17) \begin{equation}{\textbf{\textit{F}}}={\textbf{\textit{F}}}_{\text{Pot}} -{\textbf{\textit{F}}}_{\text{LE}} +{\textbf{\textit{F}}}_{\text{LEV}},\end{equation}

where ${\textbf{\textit{F}}}_{\text{Pot}}$ is given by Equation (13), ${\textbf{\textit{F}}}_{\text{LE}}$ by Equation (16) and

(18) \begin{align}\equationnumber{\normalfont(18)}{\textbf{\textit{F}}}_{\text{LEV}} &= \kappa F_{\text{LE}}\hat{\textbf{\textit{n}}}\end{align}

(19) \begin{align}\equationnumber{\normalfont(19)}\kappa &= \varrho \text{sign} (\left({\textbf{\textit{V}}}_\infty-{\textbf{\textit{v}}}_{O}\right)\cdot\hat{\textbf{\textit{n}}}).\\[6pt]\nonumber\end{align}

In general, the parameter $\varrho$ can assume any value. It is observed that $\varrho=1$ results in Polhamus^{(Reference Polhamus32)} leading edge suction analogy model. The case ${\varrho=0}$ , simply the suppression of the leading edge suction force, represents a situation where leading edge flow separation occurs, but the generated vortex does not remain “attached” to the body. The case of $\varrho<0$ can be thought to represent a situation when the leading edge vortex is on the windward side of the plate, as it may occur after a rapid flapping stroke reversal. Since intermediate values of $\varrho$ are allowed, the proposed modification is considered a generalisation of Polhamus’ leading edge suction analogy.

After substituting Equations (13), (16) and (18), and $\Gamma=\Gamma_K$ into Equation (17), the modified resultant force is obtained:

(20) \begin{align}{\textbf{\textit{F}}}{}={}&\rho\pi c \left( U_{\infty} \sin\vartheta-V_O \sin \varphi\right)\left( U_\infty \cos\vartheta-V_O \cos\varphi\right)\hat{\textbf{\textit{n}}}\nonumber\\&-\rho c^2 \pi \omega\left(\frac{1}{4}- \frac{x_O}{c}\right)\left( U_\infty \cos\vartheta-V_O \cos\varphi\right)\hat{\textbf{\textit{n}}}-\rho \frac{c}{2}\frac{d\Gamma}{dt}\hat{\textbf{\textit{n}}}\nonumber\\&+\pi \rho \frac{c^2}{4}x_O\frac{d\omega}{dt}\hat{\textbf{\textit{n}}}+\pi \rho \frac{c^2}{4} \frac{d}{dt}\left(U_\infty \sin \vartheta-V_O \sin\varphi\right)\hat{\textbf{\textit{n}}}\nonumber\\&+\kappa\rho\pi c\left(V_O\sin\varphi -U_{\infty}\sin\vartheta-\omega x_O\right)^2\hat{\textbf{\textit{n}}}.\end{align}

2.2 Unsteady Vortex Lattice Method (UVLM) formulation

The Unsteady Vortex Lattice Method (UVLM) provides a way for describing the wake motions. Here a succinct description of the method is presented; for more details of the formulation, the reader is referred to Katz and Plotkin^{(Reference Katz and Plotkin33)}.

Two types of vortex sheets are recognised: bound vortex sheets modeling the boundary layers attached to the solid surfaces, and free vortex sheets representing the wakes. The bound vortex sheets move with the bodies; as a result, a finite pressure jump develops across them. Conversely, the free vortex sheets grow and deform freely to attain force free positions. The two types of vortex sheets merge at the points of the bodies where flow separation occurs.

In the UVLM formulation, the continuous bound vortex sheets are discretised with collections of $N_\mathcal{E}$ vortices of circulation $\Gamma_i$ (Fig. 1(b)). As a result, the no-penetration condition (4) can only be satisfied at a finite number $M\le N_\mathcal{E}$ of control (collocation) points on the solid surface. By virtue of this discretisation, the contribution to the normal component of the fluid velocity at the location ${\textbf{\textit{r}}}_i\in\partial \mathcal{E}$ of the control point i due to the bound vortex sheet of intensity $\gamma_\mathcal{E}$ can be approximated as follows:

(21) \begin{align}{\textbf{\textit{u}}}_{\gamma_\mathcal{E}}({\textbf{\textit{r}}}_i)\cdot\hat{\textbf{\textit{n}}}({\textbf{\textit{r}}}_i)&=\int_{\partial \mathcal{E}}\frac{\gamma_\mathcal{E}(s)\hat{\textbf{\textit{k}}}\times\left({\textbf{\textit{r}}}_i -{\textbf{\textit{r}}}(s)\right)}{2\pi\Vert{\textbf{\textit{r}}}_i - {\textbf{\textit{r}}}(s)\Vert^2}ds\cdot\hat{\textbf{\textit{n}}}({\textbf{\textit{r}}}_i)\nonumber\\& \approx \sum_{j=1}^{N_\mathcal{E}}\Gamma_j\frac{\hat{\textbf{\textit{k}}}\times\left({\textbf{\textit{r}}}_i-{\textbf{\textit{r}}}_j \right)\cdot\hat{\textbf{\textit{n}}}({\textbf{\textit{r}}}_i)}{2\pi\Vert{\textbf{\textit{r}}}_i - {\textbf{\textit{r}}}_j\Vert^2} \nonumber\\&\approx\sum_{i=1}^{N_\mathcal{E}}\Gamma_j A_{ij},\\[3pt] \forall i&=1,2,...,M.\nonumber\end{align}

Here, ${\Gamma_j=\Delta s_j\gamma_\mathcal{E}(s_j)}$ is the intensity of vortex j and $A_{ij}$ its influence coefficient over the control point i. In a similar fashion, the free vortex sheets are discretised with collections of $N_\mathcal{W}$ vortices. Therefore, the contribution to the normal component of the fluid velocity at the location ${\textbf{\textit{r}}}_i\in\partial \mathcal{E}$ of the control point i due to the free vortex sheet of intensity $\gamma_\mathcal{W}$ can be approximated as follows:

(22) \begin{align}{\textbf{\textit{u}}}_{\gamma_\mathcal{W}}({\textbf{\textit{r}}}_i)\cdot\hat{\textbf{\textit{n}}}({\textbf{\textit{r}}}_i)&=\int_{\partial \mathcal{W}}\frac{\gamma_\mathcal{W}(s)\hat{\textbf{\textit{k}}}\times\left({\textbf{\textit{r}}}_i-{\textbf{\textit{r}}}(s) \right)}{2\pi\Vert{\textbf{\textit{r}}}_i - {\textbf{\textit{r}}}(s)\Vert^2}ds\cdot\hat{\textbf{\textit{n}}}({\textbf{\textit{r}}}_i)\nonumber\\& \approx \sum_{j=1}^{N_\mathcal{W}}\Gamma_{\mathcal{W}j}\frac{\hat{\textbf{\textit{k}}}\times\left({\textbf{\textit{r}}}_i-{\textbf{\textit{r}}}_j \right)\cdot\hat{\textbf{\textit{n}}}({\textbf{\textit{r}}}_i)}{2\pi\Vert{\textbf{\textit{r}}}_i - {\textbf{\textit{r}}}_j\Vert^2}.\\[3pt] \forall i&=1,2,...,M.\nonumber\end{align}

Once the vortices in the wakes are created, they move freely in the fluid domain with constant circulation. Their motions are determined by the local fluid velocities at the location of the free vortices. Therefore, at a given instant of time, the contribution given by Equation (22) is known. Considering the attached and free vortex sheets contributions (Equations (21) and (22), respectively), and the contribution of the free stream, the no-penetration condition (4) at each control point i yields:

(23) \begin{align}\begin{split}\sum_{j=1}^{N_\mathcal{E}}\Gamma_j A_{ij}=& -\sum_{j=1}^{N_\mathcal{W}}\Gamma_{\mathcal{W}j}\frac{\hat{\textbf{\textit{k}}}\times\left({\textbf{\textit{r}}}_i-{\textbf{\textit{r}}}_j \right)\cdot\hat{\textbf{\textit{n}}}({\textbf{\textit{r}}}_i)}{2\pi\Vert{\textbf{\textit{r}}}_i - {\textbf{\textit{r}}}_j\Vert^2}-\left({\textbf{\textit{u}}}_\infty-{\textbf{\textit{v}}}({\textbf{\textit{r}}}_i)\right)\cdot \hat{\textbf{\textit{n}}}({\textbf{\textit{r}}}_i),\end{split}\end{align}

\begin{align}&\forall i=1,2,...,M.\nonumber\end{align}

Additionally, the Kelvin Condition of conservation of circulation requires that

(24) \begin{align}\sum_{i=1}^{N_\mathcal{E}}\Gamma_j (t)+\sum_{i=1}^{N_\mathcal{W}}\Gamma_{\mathcal{W}j}(t) =\sum_{i=1}^{N_\mathcal{E}}\Gamma_j (t=0).\end{align}

These $M+1$ equations (Equations (23) and (24)) can be cast into a linear algebraic problem as follows:

(25) \begin{align}\left[\text{A}\right]\left\{\Gamma(t)\right\} = \left\{\text{RHS}(t)\right\}\!.\end{align}

Here, $\left[\text{A}\right]\in\mathbb{R}^{\left(M+ 1\right)\times N_\mathcal{E}}$ , $\left\{\Gamma(t)\right\}\in\mathbb{R}^{N_\mathcal{E}}$ , $\left\{\text{RHS}(t)\right\}\in\mathbb{R}^{M+ 1}$ . Solving this system at each time step yields the intensities of the bound vortices $\Gamma_i$ . Once the circulation of the bound vortices is known, the velocity field becomes known and the pressure distribution on the surface of the body can be computed by means of Bernoulli’s Equation (5).

2.2.1 Formation of vortex sheets at the sharp edges

It can be shown through Equation (5) that the pressure difference at a separation point S is given by: 4

(26) \begin{align}\frac{p({\textbf{\textit{r}}}_S^-,t)-p({\textbf{\textit{r}}}_S^+,t)}{\rho}{}={}&\frac{1}{2}\left(u_{\text{slip}}^+{}^2-u_{\text{slip}}^-{}^2\right)-\left.\frac{d\Gamma}{dt}\right|_{{\textbf{\textit{r}}}_S},\end{align}

where $u_{\text{slip}}$ is the slip speed on the surface, $+$ and $-$ indicates the merging surfaces “above” and “below” the separation point, respectively, and $d\Gamma/{dt}$ is the rate of generation of circulation at the separation point. From this equation, the Kutta Conditon of zero pressure jump on the trailing edge implies that the rate of generation of vorticity off of a sharp edge S is given by

(27) \begin{align}\frac{d\Gamma_\mathcal{W}({\textbf{\textit{r}}}_S)}{dt}=\frac{1}{2}\left(u_{\text{slip}}^+{}^2-u_{\text{slip}}^-{}^2\right).\end{align}

This equation was also used by Nitsche and Krasny^{(Reference Nitsche and Krasny34)} and by Jones^{(Reference Jones35)}, among others. In the numerical scheme, an equivalent point vortex is generated at time t at the separation point whose intensity is given by

(28) \begin{align}\Gamma_S^{t}&=\int_{t-\Delta t}^{t}\left.\frac{d\Gamma}{dt}\right|_{{\textbf{\textit{r}}}_S}dt\approx \left.\frac{d\Gamma}{dt}\right|_{{\textbf{\textit{r}}}_S}\Delta t+\mathcal{O}\left(\Delta t^2\right),\end{align}

and its position at a later time by

(29) \begin{align}{\textbf{\textit{r}}}_{\Gamma_S^t} &={\textbf{\textit{r}}}_S(t)+\int_t^{t+\Delta t}{\textbf{\textit{u}}}_S dt\nonumber\\&\approx{\textbf{\textit{r}}}_S(t)+{\textbf{\textit{u}}}_S(t) \Delta t+\frac{1}{2}\frac{d{\textbf{\textit{u}}}_S(t)}{dt}\Delta t^2+\mathcal{O}\left(\Delta t^3\right).\end{align}

2.2.2 Fluid force

The fluid force acting on the body is determined by integrating the pressure distribution (Equation (5)) over the surface. This integral is carried out numerically, and for a flat plate yields a force ${\textbf{\textit{F}}}_p$ that is always normal.

To reconcile the theory with the numerical method, the leading edge suction force has to be computed and added to the pressure-based force. The added contribution to the leading edge suction force due to the bound circulation is obtained from Equation (14) by replacing $\Gamma(t)$ by $\sum \Gamma_j(t)$ . It can be shown^{(Reference Katz and Plotkin33)}, Section 6.9 that the contribution due to the free vorticity ${\textbf{\textit{F}}}_{\text{LE}\mathcal{W}}$ can be computed as follows

(30) \begin{equation}{\textbf{\textit{F}}}_{\text{LE}\mathcal{W}}=\rho\sum_{i=1}^{N_\mathcal{E}}\Gamma_i\left(\sum_{j=1}^{N_\mathcal{W}}\frac{\Gamma_{\mathcal{W}j}}{2\pi}\frac{\left({\textbf{\textit{r}}}_i-{\textbf{\textit{r}}}_j\right)\times\hat{{\textbf{\textit{n}}}_3}}{\Vert{\textbf{\textit{r}}}_i - {\textbf{\textit{r}}}_j\Vert^2}\cdot\hat{\textbf{\textit{n}}}({\textbf{\textit{r}}}_i)\right)\hat{\textbf{\textit{t}}},\end{equation}

where the term in parenthesis is the total downwash (normal component of velocity) due to the wake vortices at the location of the bound vortices i. The total force on the flat wing section is computed as follows:

(31) \begin{equation}{\textbf{\textit{F}}}={\textbf{\textit{F}}}_p+{\textbf{\textit{F}}}_{\text{LEV}},\end{equation}

where ${\textbf{\textit{F}}}_{\text{LEV}}$ is given by:

(32) \begin{equation}{\textbf{\textit{F}}}_{\text{LEV}}=\kappa \Vert{\textbf{\textit{F}}}_{\text{LE}}+{\textbf{\textit{F}}}_{\text{LE}\mathcal{W}}\Vert\hat{\textbf{\textit{n}}},\end{equation}

with $\kappa$ defined in Equation (19).

This hybrid UVLM-theoretical approach to compute the fluid force is a novel contribution of this article.