2.1 Inviscid methods and the discrete vortex method

We first consider briefly some methods of modelling vortex shedding in flows in which viscous effects can be neglected in a first analysis.

The first of these models of shed vorticity is limited to bodies which induce only small perturbations in a uniform incident flow. For a streamlined body such as a plate or aerofoil at small angle-of-attack (AoA) the wake perturbations including any roll-up are small enough to allow representation by a plane sheet. The vorticity shed from the trailing edges of the body is then assumed to be convected by the free stream alone, all other effects being negligible. This thin aerofoil theory devised for steady flows by Munk [Reference Munk12] was extended to unsteady flow to analyse flutter [Reference Theodorsen13], impulsively started motion [14], effects of periodic convected gusts [Reference Sears15] and impulsive gusts [Reference Kussner16] on aerofoils. Useful analytic results can be obtained for many applications. For larger perturbations this representation may be extended by means of lifting surface panel methods to unsteady flows around more complicated 3D bodies and wakes [Reference Katz and Plotkin6] still widely used in industry for aeroelastic calculations.

The second group of methods was developed for flows with shed vortex sheets exhibiting strong roll-up. The rolled-up trailing vortices shed from lifting wings were represented by concentrated, streamwise line vortices close to the centroids of the rolled-up vortex sheets and aligned with the free stream. For high angles of attack and stronger roll-up closer to the wing as in the case of vortices separating from the leading edges of a slender delta wing this representation was improved by including the local flow around the wing in the vortex alignment but still gave significant errors [Reference Legendre17]. The analogous 2D unsteady flow model for the shedding of a rolling-up vortex from an edge, such as for starting flow around the trailing edge of an aerofoil, was similarly inaccurate [Reference Rott18].

Further improvement of the concentrated line vortex model of the leading edge vortices (LEVs) of a slender wing at high AoA [Reference Edwards19] and similarly for unsteady 2D vortex shedding and roll-up (for shock diffraction around a sharp edge [Reference Rott18]) recognised that the part of the vortex sheet joining the growing spiral vortex to its separation edge needed to be included in the model. The effect of this part was represented by a plane ‘feeder’ sheet (FS in Fig. 1) contributing to the condition of zero overall force on the shedding vortex system. With this the equation for ${\boldsymbol{r}_{\textbf{0}}}(\boldsymbol{x}) = \left( {{y_0}(x),{z_0}(x)} \right)$ the vector position in the cross-flow $\left( {y,z} \right)$ plane of the LEVs above the slender wing becomes:

(7) \begin{align}\frac{{d{{\boldsymbol{r}}_{\textbf{0}}}}}{{dx}} = \frac{{{\boldsymbol{V}}\!\left( {{\boldsymbol{r}}_0} \right)}}{{{U_\infty }}} - \frac{{{{\boldsymbol{r}}_0} - {{\boldsymbol{r}}_{\bf{e}}}}}{{{{\rm{\Gamma }}_0}}}\frac{{d{{\rm{\Gamma }}_0}}}{{dx}}\end{align}

where ${{\boldsymbol{r}}_{{\boldsymbol{e}}}}$ is the coordinate of the edge. The second term on the RHS represents the effect of the feeder sheet on the force balance. ${\boldsymbol{V}}\!\left( {{r}_0} \right) = \left( {V\!\left( {{\boldsymbol{r}_\textbf{{0}}}} \right),W\!\left( {{\boldsymbol{r}_\textbf{{0}}}} \right)} \right)$ is the velocity in the cross-flow plane induced at the vortex position and ${U_\infty }$ the free stream velocity in the $x$ direction. This line vortex model then gives the positions of the vortex centres and the additional lift due to the LEVs more accurately [Reference Brown and Michael20].

Figure 1. Line vortex + feeder sheet (FS) model of delta wing LEVs (from Ref. [Reference Brown and Michael20]).

The corresponding improved equation for the unsteady 2D separating vortex [Reference Rott18] is:

(8) \begin{align}\frac{{d{\boldsymbol{r}_\textbf{{0}}}}}{{dt}} = {\boldsymbol{V}}\!\left( {{\boldsymbol{r}_\textbf{{0}}}} \right) - \frac{{{\boldsymbol{r}_\textbf{{0}}} - {\boldsymbol{r}{_\boldsymbol{e}}}}}{{{\Gamma _0}}}\frac{{d{{\rm{\Gamma }}_0}}}{{dt}}\end{align}

or:

(9) \begin{align}\frac{{d\!\left[ {{{\rm{\Gamma }}_0}\!\left( {{\boldsymbol{r}_\textbf{{0}}} - {\boldsymbol{r}{_\boldsymbol{e}}}} \right)} \right]}}{{dt}} = {{\rm{\Gamma }}_0}{\boldsymbol{V}}\!\left( {{\boldsymbol{r}_\textbf{{0}}}} \right)\end{align}

showing that the addition of the second term on the RHS of Equation (8) equates the vortex induced force to the rate of change of the vortex moment formed by a counter-rotating vortex pair formed by the vortex and its virtual image in the body.

More detailed representation of the rolling-up sheet requires, at least, a multi-element discretisation. The Lagrangian method of solution of the vorticity transport equation:

(10) \begin{align}\frac{{\partial \boldsymbol\omega }}{{\partial t}} + \boldsymbol{U}.\nabla \boldsymbol\omega = \boldsymbol\omega .\nabla \boldsymbol{U} + \nu {\nabla ^2}\boldsymbol\omega \end{align}

tracks elements of the vorticity field $\boldsymbol\omega $ simulating the LHS convection of Equation (10), while preserving the vortex structure. The first term on the RHS of Equation (10) $\boldsymbol\omega .\nabla \boldsymbol{U}$ is zero in 2D flow and hence omitted.

The velocity field, $U$ , is computed simultaneously from the vorticity field as the solution of the Poisson equation:

(11) \begin{align}{\nabla ^2}\boldsymbol{U} = - \nabla \times \boldsymbol\omega .\end{align}

Equations (10) and (11) are sufficient with boundary conditions to specify the flow kinematically. The pressure and body forces may be determined retrospectively when required. Many different discretisations of the circulation density in a vortex sheet are possible to increase accuracy but the simplest, known as the discrete vortex method (DVM) or point vortex method, is a lumped vortex discretisation of the sheet, each short element being represented by a concentrated point vortex at the centre of the element,

The DVM was first developed [Reference Rosenhead21] as an inviscid 2-D model of rolling-up of vortex sheets in an isolated wake. It was then adapted further [Reference Clements and Maull22] to study 2D shedding of a vortex street wake from a rectangular-base body.

With externally specified separation points the inviscid DVM was found to predict the Strouhal frequency and spacing of the von Karman vortex street quite well. Alternating rolled-up vortices arose naturally in the wake due to inherent numerical instability of the wake.

The DVM is a time-stepping method, new vortices being shed into the flow at each step. The solution of Equation (10) at each time step with RHS zero for inviscid 2D flow, requires the solution of Equation (11) for the velocity field from all vortices in the flow field, most easily calculated using the 2D Biot-Savart integral solution:

(12) \begin{align}{\boldsymbol{U}_m} = \mathop \sum \limits_j \frac{{{{\rm{\Gamma }}_j}}}{{2\pi }}\frac{{\boldsymbol{k} \times \left( {{\boldsymbol{r}_m} - {\boldsymbol{r}_j}} \right)}}{{|\left( {{\boldsymbol{r}_m} - {\boldsymbol{r}_j}} \right){|^2}}}\end{align}

excluding the local singularity at $j = m$ , where $\boldsymbol{k}$ is the unit vector $\left( {0,0,1} \right)$ . Thus the operation count of the method is $O\!\left( {N_v^2} \right)$ per time step, where ${N_v}$ is the number of vortices shed which increases in proportion to the number of time-steps. While the method has attractions of simplicity, apparently reducing diffusion, long flow evolution times, with an operation count $O\!\left( {{N^3}} \right)$ for $N$ time-steps, have only become possible with the advent of large computers. Faster methods of solution of the velocity from Equation (11) have been developed to obviate this problem (eg multipole [Reference Greengard and Rohklin23], vortex-in-cell [Reference Baker24, Reference Stansby and Dixon25] and other mesh methods) which have lead to the possibility of much longer and better resolved flow computations.

Point vortex discretisation of a separating vortex sheet in inviscid flow, including the concentrated vortex models, fails locally at the separation edge since without continuous sheet representation there the only possible flows are either a continuous flow around the edge, with a singularity, or a full stagnation separation point, both incorrect. More accurate discretisation methods have been developed to overcome this using continuous elements and integration along the sheet, which are more appropriate [Reference Giesing26, Reference Pullin27].

Numerical instability is inherent increasing with resolution, see Refs [Reference Moore28, Reference Van der Vooren29], but often controlled by a desingularising addition in the denominator of the influence kernel in Equation (12), see e.g. Ref. [Reference Krasny30]. It increases rapidly in the close-packed inner turns of a vortex spiral but can be averted by amalgamating the inner turns progressively into a single central core vortex. An example of this for a rolling up spiral vortex sheet shedding from the sharp edge of a plate in a starting flow, is shown with core amalgamation (red curve) and without (black curve) in Fig. 2. The two curves agree well until the third arm of the spiral where the one without amalgamation becomes unstable and breaks up. Similar strong roll-up occurs at all free ends of vortex sheets unless the vorticity gradient is weak.

Figure 2. DVM computation of a rolling up vortex sheet shed from the edge of a flat plate Reference Graham[31].

2.2 Modelling of viscous flows by the DVM

Modelling of flows with appreciable viscous effects is more usually (and virtually always for 3-dimensional flows) carried out using the Eulerian mesh methods of computational fluid dynamics (CFD). Numerical simulations using Lagrangian discrete vortex methods for the convection, have been carried out for 2D unsteady flows where there is some advantage that only two flow variables need be solved for, compared with four for CFD. With the point vortices desingularised by a finite support viscous diffusion has been modelled by the Random Walk Method [Reference Chorin32], see Ref. [Reference Roberts33] regarding the high number density required. With sufficient vortices to ensure overlapping cores the vorticity exchange method [Reference Leonard34] has proved accurate for laminar flows, hybrid methods in which the Lagrangian vortices are projected onto a Eulerian mesh to carry out the diffusion have also been developed. A result using this type of hybrid method for flow around a circular cylinder at a Reynolds number of 200 is shown in Fig. 3 The hybrid method can be extended with subgrid modelling of turbulence to simulate vortical flows at higher Reynolds numbers [Reference Willden and Graham35].

Figure 3. DVM simulation of the wake of a circular cylinder $Re = 200,{S_t} = 0.196,{C_D} = 1.22,{G_A}rkel{l_Z}$ .

2.3 Scaling of 2D vortex shedding from a sharp edge

The following scaling analysis was first derived for a local edge flow [Reference Graham36] and later specifically for an aerofoil in starting flow [Reference Graham37]. The edge region of a 2D body such as an aerofoil, formed by two locally plane sections of its surface meeting with included angle $\delta $ , is subject to a starting flow at time $t = 0$ , as shown in Fig. 4. The local flow near the edge can be divided into two velocity components: ${U_0}$ symmetric with respect to the bisector of the edge angle and ${V_0}$ normal to it. Assuming effectively inviscid flow and $\delta \lt \pi $ the ${V_0}$ component generates an infinite velocity at the edge which causes separation and shedding of a vortex sheet. There are two cases to consider: (1) ${U_0}(t) = {\hat U_0}{t^\beta }$ and ${V_0}(t) = {\hat V_0}{t^\beta }$ , $\beta \geqslant 0$ , representing starting motion of a body in an initially stationary fluid, and (2) ${U_0}(t) = {\hat U_0}$ and ${V_0}(t) = {\hat V_0}{t^\beta }$ , $\beta \geqslant 0$ representing the body starting pitching or transverse motion in a pre-existing uniform stream.

Figure 4. Diagram of vortex ${\rm{\Gamma }}$ shed from a sharp edge of internal angle $\delta $ . Imposed velocity field ${U_0}$ symmetric, ${V_0}$ antisymmetric, ${q_e}\!\left( r \right)$ surface velocity at distance $r$ from edge.

While the vortex sheet which is shed as a result of flow separation at the edge is still close to it the flow there has no length scale given explicitly by the body. The length-scale ${L_e}$ must be evaluated from the product of a local velocity-scale and the time $t$ . The local velocity close to the edge, ${q_e}\!\left( r \right) = {Q_e}{\boldsymbol{r}^{\left( {1/\kappa } \right) - 1}}$ , where ${Q_e}$ is a constant which does depend implicitly on the overall body scale, is singular at the edge since $\kappa = 2 - \delta /\pi $ . Therefore, it must be evaluated at a non-zero distance $r \gt 0$ from the edge to define a velocity-scale and the evaluation distance must itself scale on ${L_e}$ . For simplicity let $r = {L_e}$ , but it could be any multiple.

Therefore ${L_e} = {Q_e}L_e^{\left( {1/\kappa } \right) - 1}t$ and hence ${L_e} = {({Q_e}t)^{\kappa /\left( {2\kappa - 1} \right)}}$ .

The parameter ${Q_e}$ which defines the edge flow locally, is proportional to the imposed velocity ${V_0}$ . It can be calculated in practice from values of ${q_e}\!\left( r \right)$ at points on the body surface close to the edge, averaging the product ${q_e}\!\left( r \right){\boldsymbol{r}^{\left( {1 - 1/\kappa } \right)}}$ .

The development of the shedding vortex sheet can only scale on this local length scale. ${L_e}$ , grows with time as:

(13) \begin{align}{L_e} = {({Q_e}t)^{\kappa /\left( {2\kappa - 1} \right)}}\end{align}

and the sheet therefore rolls up initially into a growing, self-similar spiral.

Let the total circulation shed since the start depend on time as:

(14) \begin{align}{\rm{\Gamma }}(t) = {{\rm{\Gamma }}_0}{t^\mu }\end{align}

The local velocity induced by the vortex and its image decays as ${\boldsymbol{r}^{ - 2}}$ but the potential on which the force depends as an integral decays in 2D flow as ${\boldsymbol{r}^{ - 1}}$ and is therefore not local but must depend on a global body length-scale parameter such as a mean diameter. However, this parameter can be derived from ${Q_e}$ which, although evaluated locally, contains the required body-scale information.

Two equations define the flow development. The first is the flow separation condition. In any single or multi-point vortex representation of the shedding it is taken as the Kutta (full stagnation point) condition, since the alternative is a velocity singularity at the edge. Using it does introduce some small errors [Reference Graham38]. Higher fidelity vortex shedding models based on sheet-like elements give better accuracy with a flow stagnation point only on the separated side [Reference Pullin27]. All lead to the same inviscid scaling.

The second equation is for the convection of the vorticity in the physical plane. There are two possibilities initially depending on which term in the equation dominates. If ${V_0}$ dominates, as often the case, the circulation ${\rm{\Gamma }}$ of the shed vorticity (a single vortex at ${\boldsymbol{r}_\textbf{{0}}} = \left( {{x_0},{y_0}} \right)$ is shown) increases as ${t^\mu }$ where:

(15) \begin{align}\mu = \frac{{2\kappa \beta + 1}}{{2\kappa - 1}},\left[ {39} \right]\end{align}

For impulsively started flow $\beta = 0$ and $\mu = \frac{1}{{2\kappa - 1}}$ , for constant acceleration $\beta = 1$ and $\mu = \frac{{2\kappa + 1}}{{2\kappa - 1}}$ .

The second possibility occurs when the symmetric component ${U_0}$ dominates, for example when a body immersed in a continuous flow symmetrically off the edge is subject to a ${V_0}$ component around the edge which starts to accelerate gradually from zero. Then

(16) \begin{align}\mu = \frac{1}{{2\kappa - 2}} + \beta, \left[ {39} \right]\end{align}

For constant acceleration of ${V_0}$ $\left( {\beta = 1} \right)$ , $\mu = \frac{{2\kappa - 1}}{{2\kappa - 2}}$ . So for an aerofoil commencing heave motion with constant acceleration in a mean free-stream $\mu \sim 3/2$ since $\kappa \sim 2$ .

The single point vortex model applied to a starting flow about a sharp edge has a small-time analytic solution. The result is shown for a Karman-Trefftz aerofoil in impulsively started flow compared with longer time DVM simulation and the Wagner linearised theory [14] in Figs 5, 6 taken from [Reference Graham37]. For small times the shed concentrated vortex and the DVM spiral centre leave the trailing edge in a direction normal to the edge angle bisector before their paths bend towards the free-stream direction. More details of the scaling of the flow are given in Graham [Reference Graham37]. For impulsively started aerofoil-free-stream flows these isolated edge results are of more theoretical than practical use since the starting vortex is rapidly carried away from the locality of the trailing edge, but have greater validity for oscillatory and slender body flows.

Figure 5. DVM computation of the TEV formed above the trailing edge of an aerofoil at two times after an impulsive start [Reference Lamb40].

Figure 6. DVM computation of the lift on a Karman-Trefftz aerofoil after an impulsive start [Reference Lamb40].

The scaling can be applied equally to leading edge vortex (LEV) as to trailing edge vortex (TEV) shedding from an aerofoil. Viscous effects are more important in the LEV case since the shed vortex may remain or move downstream close to the aerofoil surface. For effectively sharp leading-edge separation the analysis is as for the TEV with the sign of ${U_\infty }$ changed to towards the edge rather than away from it. For a flat plate or thin aerofoil in a starting flow at a moderate to large AoA, both ${U_\infty }$ and ${V_\infty }$ have the same dependence on $t$ . The shedding from both edges depends only on the ${V_\infty }$ component of free stream normal to the plate and is independent of ${U_\infty }$ . Therefore reversing the sign of ${U_\infty }$ in the local edge coordinate system has no effect on the initial shedding path of the vortex core, the roll-up is the same sign, but the geometry is reversed. The LEV core initially leaves in a direction normal to the plate at the leading edge with scaling as for the TEV. For longer times the path of the vortex moves from the normal direction to a downstream direction now over the plate’s surface, but its image effect is therefore stronger causing a stronger opposing back-flow which slows its downstream convection. This LEV flow has been studied extensively with PIV and computer simulation [Reference Pitt Ford and Babinsky39]. For more gradually starting heave or pitch motion the concentrated vortex model predicts a rearwards core vortex path leaving the surface at the edge tangentially with ${y_0}\sim x_0^{3/2}$ . Leading edge separations due to rapid changes in effective AoA are a major component of dynamic stall, see Section 3.3.2.

2.4 Time-dependent slender wing flow

In steady flow, the LEVs shed by a plane delta wing, as shown in Fig. 1, retain conical similarity. The locations ${y_0}(x),{z_0}(x)$ of the vortex cores in the cross-flow plane $\left( x \right)$ predicted by the concentrated model move gradually inwards from the edge ${y_e}(x)$ following the path ${y_e} - {y_0}\sim z_0^{2/3}$ as AoA is increased [Reference Brown and Michael20].

For time-dependent flow such as following an impulsive start, sudden heave or roll motion or oscillatory motion the time dependence combines with rearward propagation of the shedding vorticity feeding the LEV. The streamwise derivative in Equation (7) is replaced by the convective derivative $\frac{\partial }{{\partial t}} + {U_\infty }\frac{\partial }{{\partial x}}$ in the corresponding vortex force balance equation:

(17) \begin{align}\left( {\frac{\partial }{{\partial t}} + {U_\infty }\frac{\partial }{{\partial x}}} \right)\left[ {{{\rm{\Gamma }}_0}\left( {{\boldsymbol{r}_\textbf{{0}}} - {\boldsymbol{r}{_\boldsymbol{e}}}} \right)} \right] = {{\rm{\Gamma }}_0}{U_0}\!\left( {{\boldsymbol{r}_\textbf{{0}}}} \right)\end{align}

Similar vortex shedding occurs from the bilges and keels of long ship hulls and is important in determining vortex damping of time-dependent roll, and other sea-keeping motions. For small amplitudes around a slender body of constant section the vortex shedding can be considered to be as from an isolated slender edge in a free stream ${U_\infty }$ with a cross-flow ${V_0}$ , due to AoA or lee-way, starting impulsively. Equation (17) is hyperbolic and solved along characteristics $t + x/{U_\infty } = $ constant. Assuming the vortex shedding to develop from the upstream end of the edge, eg the start of a long slender rectangular planform wing or of a bilge-keel without other disturbances, the solution shows the effect of the upstream origin of shedding spreading rearwards. Considering an impulsively started flow ${U_\infty }$ and defining a convection length, $X = {U_\infty }T$ , which the influence of the upstream end has just reached at a specific time $T$ since the start, upstream of $X$ the growth of the vortex depends on distance $x({\lt} X)$ from the upstream origin. In this region its normal distance from the plate, ${z_0}(x)$ follows the envelope $z/s = 0.25{(x/s)^{2/3}}{({V_0}/{U_\infty })^{2/3}}$ . Downstream of $X$ vortex growth is dependent on time $t$ since the start only and is independent of lengthwise distance $x({\gt} X)$ . Here ${z_0}$ grows as ${t^{2/3}}$ until $t = x/{U_\infty }$ , Fig. 7.

Figure 7. Vertical distance of the growing vortex above a streamwise edge of a slender, rectangular wing at incidence in an impulsively started flow. The horizontal lines are at equal steps in time.

If the cross-section geometry of the shedding body changes with streamwise distance $x$ , streamwise growth for the same conditions occurs over the whole range of $x$ but at different rates on either side of the convection length $X$ . For a slender, plane, delta wing, semi-span $s(x) = \varepsilon x$ , at AoA $ = \alpha $ , in the upstream region $x \lt X = {U_\infty }T$ a growing length of steady conical flow develops from the apex, Fig. 8 in agreement to first order in $\alpha /\varepsilon $ with the results given by Ref. [Reference Brown and Michael20]. In the downstream region $x \gt X$ , ${z_0}\sim {({U_\infty }t)^{2/3}}{s^{1/3}}$ , which while independent of the apex still increases here with lengthwise distance $x$ since $\frac{{ds}}{{dx}} \gt 0$ .

Figure 8. Vertical distance of the growing vortex above the leading edge of a slender delta wing at incidence in an impulsively started flow. The transverse lines are each for constant time separated by equal time-steps.

Therefore for gust or operational frequencies quasi-steady aerodynamics as for $x \lt {U_\infty }T$ are adequate because aircraft speeds are high and the convective distances associated with the time-scales, with the exception of those for flutter, are much larger than the wing chord. In contrast for sea-keeping and manoeuvring motions of ships the axial velocities are lower and the bilge lengths longer so that time dependency effects are important.

3.1 The impulse method

The conventional method of calculating forces in CFD by integrating the pressure and the shear stress on the surface of a body can be difficult to apply because of the need to fully resolve the boundary layers on the body where the velocity and vorticity gradients are very large. In many unsteady flow applications, direct measurement of pressures or forces is also difficult because of structural resonances, time lags and body geometry problems. For unsteady flows the concept of the fluid dynamic impulse applied to the entire vorticity field, called the Impulse Method, is a useful alternative for the calculation of aerodynamic or hydrodynamic forces. It has gained the interest of the CFD community, especially in the context of vorticity-based algorithms where pressure is not explicitly evaluated. In experimental work, the impulse method is useful for evaluation of forces when combined with techniques such as Particle Image Velocimetry (PIV) which provides time-sequenced velocity fields but not pressure data.

3.1.1 The original method

The impulse method applies Newton’s third law directly to the flow field through the momentum equation in integral form without involving the Navier-Stokes equations explicitly. In classical steady flow aerodynamics, Prandtl’s steady vortex-force theory, the Kutta-Joukowski theorem in 2D flow, shows the advantage of linking force directly to shed vorticity. More generally the impulse method [Reference Lamb40] applies to unsteady flow, the force on a body being expressed as the time derivative of an impulse over the infinite flow domain, which is always unsteady. It is valid for inviscid and viscous flow at any Reynolds number and can be evaluated as a linear superposition of the impulses of all vortical structures present. In Lamb’s original work, the flow fields of point vortices were superposed onto a background potential flow. The fluid impulse on the body was shown to be generated by the growth and motion of counter-rotating vortex pairs of equal and opposite strength.

Lamb’s impulse method has been developed into an aerodynamic lift theory for an aerofoil in unsteady motion [Reference Von Karman and Sears41]. The unsteady part of the force was divided into two components: added-mass force caused by the acceleration of the fluid in response to the aerofoil’s acceleration and circulatory force which could be evaluated from the vortices shed into the flow field.

Fluctuating lift and drag formulae have been derived [Reference Phillips42] for a 2D body by integration of the acceleration field of the wake over the wake domain. More general formulae for vorticity-induced aerodynamic forces and moments on multiple bodies in unsteady, viscous flows were proposed by Wu [Reference Wu9]. These, related the forces on one or more immersed solid bodies of arbitrary shape and motion in 2- or 3D flows to the time variation of vorticity-moment integrals, by applying the momentum theorem to the infinite control volume ${R_\infty }$ , including both the fluid region and the solid body regions, where the induced flow appears as that of a time-dependent vortex dipole in Equation (5) [Reference Wu9]. Because of this, the impulse method does not separate the individual body forces in a multi-body group. This work was one of the first to present a rigorous analytical model for calculating aerodynamic forces and moments in viscous flows.

(18) \begin{align}\boldsymbol{F} = - \frac{1}{{\mathbb{D} - 1}}\rho \frac{{d\boldsymbol\alpha }}{{dt}} + \rho \mathop \sum \limits_{i = 1}^N \frac{d}{{dt}}\mathop \int \nolimits_{{R_{iB}}} {\boldsymbol{U}_{iB}}dR,\end{align}

where $\mathbb{D} = 2$ in 2D flows, $\mathbb{D} = 3$ in 3D flows, and $\boldsymbol\alpha $ the first moment of the vorticity field defined in Equation (6). In 2D flows $\boldsymbol\alpha = \left( {{\alpha _x},{\alpha _y}} \right)$ , where

(19) \begin{align}{\alpha _x} = \mathop \int \nolimits_{{R_\infty }} y\omega dR,{\rm{\;\;\;\;}}{\alpha _y} = - \mathop \int \nolimits_{{R_\infty }} x\omega dR.\end{align}

The first term in Equation (18) is proportional to the time rate of change of the first vorticity moment in the total region occupied by the fluids and solid bodies. The second term in Equation (18) is the inertial force due to the acceleration of the mass of fluid around the bodies. A comprehensive and rigorous derivation of this equation was also given [Reference Noca43]. Equation (18) is the general form of the unsteady force acting on a body assuming incompressible flow, a no-penetration boundary condition at the fluid-body surface, and ${R_\infty }$ enclosing all the bodies and shed vorticity.

The moment is obtained similarly and has the same form for 2- and 3D flows [Reference Wu9],

(20) \begin{align}\boldsymbol{l} \times \boldsymbol{F} = \frac{1}{2}\rho \frac{{d\boldsymbol\beta }}{{dt}} + \rho \mathop \sum \limits_{i = 1}^N \frac{d}{{dt}}\mathop \int \nolimits_{{R_{iB}}} \boldsymbol{r} \times {\boldsymbol{U}_{iB}}dR,\end{align}

where $\boldsymbol\beta $ is defined as the second moment of the vorticity field [Reference Wu9],

(21) \begin{align}\beta = \mathop \int \nolimits_{{R_\infty }} {\boldsymbol{r}^2}\boldsymbol\omega dR.\end{align}

The impulse method has been applied in many unsteady flow cases [Reference Lin and Rockwell44–Reference Devoria and Ringuette46], which show load estimations agreeing well with direct measurements. It also permits considerably simplified formulae for forces induced by vortices while they remain local to a body shedding edge.

3.1.2 Impulse theory for a finite region

The integral for the aerodynamic force as in Equation (18) and similarly for the aerodynamic moment as in Equation (20) require description of the vorticity over the whole flow field, when the far field should be irrelevant. Noca [Reference Noca, Shield and Jeon47] presented an accurate expression for the evaluation of instantaneous forces on a body in an incompressible cross-flow from knowledge of the velocity and vorticity field in a finite, arbitrarily chosen domain ${R_m}$ surrounding the body and closed by a surface ${S_m}$ over which residual stress is integrated.

(22) \begin{align}\boldsymbol{F} = - \frac{d}{{dt}}\mathop \int \nolimits_{{R_m}} \rho \boldsymbol{U}dR + \mathop \oint \nolimits_{{S_m}} \boldsymbol{n} \cdot {\boldsymbol{\Sigma }}dS,\end{align}

where ${\boldsymbol{\Sigma }}$ is the stress tensor. If the domain is kept relatively small Equation (22) becomes

(23) \begin{align}\boldsymbol{F} = - \frac{1}{{\mathbb{D} - 1}}\rho \frac{{{d_{{R_m}}}}}{{dt}} - \frac{1}{{\mathbb{D} - 1}}\mathop \sum \limits_{j = 1}^N \frac{d}{{dt}}\mathop \oint \nolimits_{{S_j}} \boldsymbol{r} \times \left( {\boldsymbol{n} \times \boldsymbol{U}} \right)dS + \mathop \oint \nolimits_{{S_m}} \boldsymbol{n} \cdot {\boldsymbol{\Theta }}dS.\end{align}

Here ${\boldsymbol\alpha _{{R_m}}}$ is the first moment of the vorticity field in the control volume ${R_m}$ (including the fluid region and the solid bodies) bounded by the control surface ${S_m}$ .

(24) \begin{align}{\boldsymbol\alpha _{{R_m}}} = \mathop \int \nolimits_{{R_m}} \boldsymbol{r} \times \boldsymbol\omega dR\end{align}

The detailed derivation is given in Noca, where ${\boldsymbol{\Theta }}$ is a tensor related to the velocity, vorticity, and viscous stress [Reference Noca43].

The term $\mathop \oint \nolimits_{{S_m}} \boldsymbol{n} \cdot {\rm{\Theta }}dS$ vanishes when the surface ${S_m}$ is taken to infinity and the system is started from rest. For a non-rotating cylinder, the last term in Equation (23) becomes:

(25) \begin{align}\mathop \oint \nolimits_{{S_m}} \boldsymbol{n} \cdot {\boldsymbol{\Theta }}dS = \mathop \oint \nolimits_{{S_m}} \left[ {\frac{1}{2}{U^2}\boldsymbol{n} - \left( {\boldsymbol{n} \cdot \boldsymbol{U}} \right)\boldsymbol{U}} \right]dS = \mathop \int \nolimits_{{R_m}} \boldsymbol{U} \times \boldsymbol\omega dR.\end{align}

Noca et al. [Reference Noca48] applied the above equation for the unsteady force to a transverse oscillating circular cylinder placed in a free stream and compared the results with those obtained by the ‘momentum equation’ and the ‘flux equation’ choosing different arbitrary and finite control volumes enclosing the body. The ‘flux equation’ for the unsteady force $\boldsymbol{F}$ contains only the control-surface integral of the velocity and its derivatives, which is useful as in most experiments the flow adjacent to a body surface in the boundary layer is highly under-resolved.

(26) \begin{align}\boldsymbol{F} = \mathop \oint \nolimits_{{S_m}} \hat{\boldsymbol{n}} \cdot {\boldsymbol\gamma _{flux}}dS - \mathop \sum \limits_{j = 1}^N \mathop \oint \nolimits_{{S_j}} \hat{\boldsymbol{n}} \cdot \left( {\boldsymbol{U} - {\boldsymbol{U}_{Sj}}} \right)\boldsymbol{U}dS - \frac{d}{{dt}}\mathop \sum \limits_{j = 1}^N \mathop \oint \nolimits_{{S_j}} \hat{\boldsymbol{n}} \cdot \left( {\boldsymbol{U}\boldsymbol{r}} \right)dS,\end{align}

where

(27) \begin{align}{\boldsymbol\gamma _{flux}} & = \frac{1}{2}{U^2}\boldsymbol{I} - \boldsymbol{U}\boldsymbol{U} - \frac{1}{{\mathbb{D} - 1}}\boldsymbol{U}\!\left( {\boldsymbol{r} \times \boldsymbol\omega } \right) + \frac{1}{{\mathbb{D} - 1}}\boldsymbol\omega \!\left( {\boldsymbol{r} \times \boldsymbol{U}} \right)\nonumber\\[4pt]& - \frac{1}{{\mathbb{D} - 1}}\left[ {\left( {\boldsymbol{r} \cdot \frac{{\partial \boldsymbol{U}}}{{\partial t}}} \right)\boldsymbol{I} - \boldsymbol{r}\frac{{\partial \boldsymbol{U}}}{{\partial t}} + \left( {\mathbb{D} - 1} \right)\frac{{\partial \boldsymbol{U}}}{{\partial t}}\boldsymbol{r}} \right]\\[4pt]& + \frac{1}{{\mathbb{D} - 1}}\left[ {\boldsymbol{r} \cdot \left( {\nabla \cdot \boldsymbol{T}} \right)\boldsymbol{I} - \boldsymbol{r}\!\left( {\nabla \cdot \boldsymbol{T}} \right)} \right] + \boldsymbol{T},\nonumber\end{align}

and ${\boldsymbol{U}_{Sj}}$ is the ${j^{{\rm{th}}}}$ -body wall velocity. Their work shows the ability of these formulae to predict unsteady force using truncated domains for relatively high Reynolds numbers. The intrinsic three-dimensionality of a nominally 2D flow is not considered. For low Reynolds number flows, in which the force coefficients are small, a change of domain size or relocation of the coordinates doesn’t provide completely self-consistent results though the comparison is still satisfactory. This may be caused by the asymmetric end conditions and oblique vortex shedding or the ’moment arm’ $\boldsymbol{r}$ [Reference Noca48].

The unsteady loads exerted on bodies moving in a fluid strongly depend on the local kinematics and vortex flow structures. To calculate the force using a finite domain enclosing the body, Wu [Reference Wu, Ma and Zhou49], proposed the ‘derivative-moment transformations’ (DMT), a formulation for force and moment, which still includes the pressure term. Wu et al. [Reference Wu, Lu and Zhuang50] used three variants of the DMT-based equations to study the dynamic forces on circular cylinders due to local flow structures and their evolution in unsteady viscous flow. A unified explanation of the influence on the body forces of the evolution of each time-varying flow structure was given. The rich forms and wide range of applications of the DMT-based force and moment expression make it popular.

Although not as popular as Wu’s force formula [Reference Wu9], Noca’s impulse equation finds application in many fields, including force prediction from a PIV field [Reference Siala and Liburdy51, Reference Siala, Fard and Liburdy52], discussed later. Li and Lu [Reference Li and Lu53] discovered when studying flapping locomotion using the DMT formula, that force and power are influenced strongly by the vortical structures near the body.

3.1.3 Force decomposition in impulse theory

It is better, if possible, to derive a procedure to obtain forces which only requires flow data within a finite domain instead of an infinitely large domain. However even so obtaining the complete description of the velocity/vorticity data throughout the domain can be difficult in some practical cases. For example, in experimental measurements such as PIV, the crucial area within the boundary-layer often cannot be adequately resolved at high Reynolds number. If the diffusion of the vorticity generated at the body surface cannot be captured by the experimental measurements, the missing vorticity information leads to a deficit in force prediction. Graham et al. [Reference Graham, Pitt Ford and Babinsky54] showed that the vortex sheet required to fulfil the non-slip boundary condition represents the incomplete measurement of vorticity on the body surface. The vortical flow field is decomposed into a vortex sheet associated with the body motion, and the vorticity field measured from the experiment as shown in Fig. 9. This surface vortex sheet arises naturally in inviscid flow, replacing the boundary-layer vorticity as a result of the slip boundary condition. It appears similar to the bound vortex sheet representations in some numerical vortex panel methods for inviscid flows. The impulse defined in Equation (6) can then be evaluated from the contribution from the measured vorticity in the fluid ${{\rm{\Gamma }}_k}$ and the contribution from the surface vortex sheet $\gamma $ :

(28) \begin{align}{\alpha _{{R_m},x}} = \alpha _{{R_m},x}^{kv} + \alpha _{{R_m},x}^{vs},{\rm{\;\;\;\;}}{\alpha _{{R_m},y}} = \alpha _{{R_m},y}^{kv} + \alpha _{{R_m},y}^{vs},\end{align}

where

(29) \begin{align}\alpha _{{R_m},x}^{kv} = \mathop \sum \limits_k {y_k}{{\rm{\Gamma }}_k},{\rm{\;\;\;\;}}\alpha _{{R_m},y}^{kv} = - \mathop \sum \limits_k {x_k}{{\rm{\Gamma }}_k},\end{align}

and

(30) \begin{align}\alpha _{{R_m},x}^{vs} = \oint y\gamma dl,{\rm{\;\;\;\;}}\alpha _{{R_m},y}^{vs} = - \oint x\gamma dl.\end{align}

Figure 9. Schematic representation of the decomposition of the vortical flow field into a vortex sheet associated with the body motion and the vorticity field measured from the experiment [Reference Wang56].

The closed-loop integrals are around the body surface.

The vortex sheet impulse can be further broken down into a circulatory part induced by the shed vorticity in the field and an added mass part due to the body movement

(31) \begin{align}\alpha _{{R_m},x}^{vs} = \alpha _{{R_m},x}^c + \alpha _{{R_m},x}^{am},{\rm{\;\;\;\;}}\alpha _{{R_m},y}^{vs} = \alpha _{{R_m},y}^c + \alpha _{{R_m},y}^{am}.\end{align}

The details on how to calculate these terms can be found in Ref. [Reference Graham, Pitt Ford and Babinsky54]. Implementing this method for a 2D flat-plate flow and translating wing experiments, they found the vortex-sheet contribution to be significant and the correction needed to be added whenever the data was incomplete for the flow around the body surface.

3.2 Impulse integral: Vortex-induced forces due to vortices close to a body edge

3.2.1 Force calculation and aerofoil starting flows

The force induced on a body in 2D flow by a vortex ${\rm{\Gamma }}(t)$ at $\left( {{x_0}(t),{y_0}(t)} \right)$ shed from and still close to an edge of internal angle $\delta $ , Fig. 4, is shown in the Appendix A to be:

(32) \begin{align}\boldsymbol{F} = - \rho \frac{d}{{dt}}\left\{ {{\rm{\Gamma }}\boldsymbol{k} \times \left( {{\boldsymbol\zeta _0} - {\boldsymbol\zeta _{I0}}} \right)} \right\}\end{align}

where ${\boldsymbol\zeta _0}$ is the coordinate of the transformed vortex position under the transformation $\boldsymbol{z}\!\left( { = \left( {r,\theta } \right)} \right) \to \boldsymbol\zeta $ , of the body to a circle of radius $R$ , and ${\boldsymbol\zeta _{I0}}$ is its image point in the circle. Evaluated close to the edge, $\boldsymbol{z} = 0$ , $\boldsymbol\zeta = 0$ , the transformation can be written as $\boldsymbol\zeta = {(Cr)^{1/\kappa }}{R^{1 - 1/\kappa }}[{\rm{cos}}( {\theta /\kappa }),{\rm{sin}} ( {\theta /\kappa }) ]$ , where C is a constant dependant on the body geometry. $\theta $ is measured from the external bisector of the edge angle and $\kappa = 2 - \delta /\pi $ . If the vortex is close to the edge ${\boldsymbol\zeta _0} - {\boldsymbol\zeta _{I0}} = [2\!\left( {C{\boldsymbol{r}_\textbf{{0}}}{)^{1/\kappa }}{R^{1 - 1/\kappa }}{\rm{cos}}( {\theta /\kappa }),0} \right]$ . The above result and the scaling analysis were given in complex notation in Ref. [Reference Graham36].

Therefore the force amplitude on a body due to an array of shed vortices, ${{\rm{\Gamma }}_j}$ , is:

(33) \begin{align}F = \rho \frac{{\kappa {Q_e}}}{{{V_\infty }}}\frac{d}{{dt}}\sum \left[ {2{{\rm{\Gamma }}_j}r_j^{1/\kappa }{\rm{cos}}\!\left( {{\theta _j}/\kappa } \right)} \right]\end{align}

with ${Q_e}$ , the edge velocity singularity parameter defined in Section 2.3 providing the influence of the overall body length-scale. The force $F$ acts in the direction normal to the bisector of the edge angle.

If a rolling-up shedding vortex can be represented by a single concentrated vortex ${{\rm{\Gamma }}_0}$ at ${\underline z _0}$ then the Kutta condition for separation at the edge allows further simplification of Equation (33) to:

(34) \begin{align}F = \pi \rho {\kappa ^2}\frac{d}{{dt}}\frac{{Q_e^2r_0^{2/\kappa }}}{{{V_\infty }}}\end{align}

The advantage of using these two formulations for the vortex force is that they can be evaluated solely from quantities which can be calculated or measured locally near the edge without reference to the rest of the body shape nor the implicit transformation of its perimeter to a circle. The body-scale factor required for the force is contained in the parameter $\frac{{{Q_e}}}{{{V_\infty }}}$ which is obtainable from local conditions at the edge. Applying Equation (A6) for the vortex shed from the trailing edge of a Karman-Trefftz aerofoil in impulsively started flow shows that this formula predicts the results given in Ref. [Reference Graham37], Equation (14).

Force calculations involving contributions from vortices further away which cannot be considered local to the edge can be computed using the time derivative of the vortex dipole formula given in the appendix, Equation (A10).

For an impulsively started flow $\beta = 0$ past an aerofoil for which $\kappa \sim 2$ :

(35) \begin{align}F\sim {({U_\infty }t/c)^{ - 1/3}}\end{align}

where c is the chord of the aerofoil. After the initial singularity $\sim {t^{ - 1/3}}$ in the vortex lift due to the shedding of the starting vortex, the lift decreases rapidly, and a DVM simulation shows that it asymptotes towards the curve predicted by the linearised Wagner theory [14], Fig. 6 [Reference Graham37]. The singularity in the vortex induced lift should be distinguished from the singularity in the potential flow added mass of the aerofoil due to infinite acceleration of the flow at the start.

If however the velocity around the aerofoil trailing edge causing the separation grows more gradually, ${V_0},\sim {t^\mu }$ with $\mu \gt 1/2$ , in the presence of an on-going free-stream the shed vortex convects downstream immediately, closer to the stream direction, with a weaker roll-up. There is no singularity in the vortex lift which remains closer to linear theory [Reference Graham37].

3.3 Strong vortex development due to rapid change of geometry

For more rapid changes in effective AoA the strong roll up previously described occurs and induces an initial increase in lift force which can be well above the linearised prediction. Such a flow occurs when an upper wing-surface spoiler is rapidly deployed during forward flight. This is intended to dump lift rapidly and add drag if the aircraft encounters a sudden gust or needs to make a rapid manoeuvre. However it has been shown e.g. [Reference Bearman, Graham and Kalkanis7] that in rapid deployment the growth of the vortex at the tip of the spoiler generates a strong ’adverse’ effect. This initially increases the lift on the wing significantly before it falls with oscillations over a time for the aircraft to travel $\sim 4c$ (wing chords) to a new value below the previous mean, appropriate to the final angle of deployment and the intended goal, Fig. 10.

Figure 10. Lift Coefficient time history due to rapid deployment to ${90^{\rm{o}}}$ of a $0.1c$ spoiler, at $0.7c$ for ${12^o}$ AoA [Reference Bearman, Graham and Kalkanis7].

3.3.1 Flapping flight

The high lift-to-drag ratio which a large bird generates in forward flight by fixing its wing at a small angle-of-attack, within the streamlined flow regime, has led to conventional fixed-wing aerodynamics. However, in many biological flight regimes which cover a wide Reynolds number range ( $Re = {10^0}$ to ${10^6}$ ), conventional wing aerodynamics are insufficient and vortex-dominated unsteady separated flows are essential for the excellent flight stability and maneuverability which many flying creatures achieve [Reference Liu, Wang and Liu55]. In particular insects [Reference Wang56], birds [Reference Wu57], and bats [Reference Muijres, Johansson, Barfield, Wolf, Spedding and Hedenström58] all leverage the large-scale coherent vortices generated through the flapping wing mechanism [Reference Chen and Lentink59] to generate the lift and thrust which enables the larger birds to take off from the ground or perch and to generate lift which enables smaller birds and insects to hover or to slow forward-flight and execute soft landings. The impulsive wing opening motion and resulting vortex formation generating additional lift, seen in this ’clap and fling’ wing motion [Reference Weis-Fogh60] have some similarities with the rapidly deployed spoiler, considered above. Edwards and Cheng [Reference Edwards and Cheng61] have studied the effect of the inviscid formation of the vortex at each wing-tip of a pair of wings in a rapidly opening ‘V’ configuration without forward speed. These vortices were shown to significantly increase the lift in this phase of wing motion. More detailed studies of insect flight, e.g. the Hawk Moth [Reference Ellington62], have suggested a more complex flow may be involved. A LEV running along the wing develops due to rapid forward rotation of the wing about the ‘shoulder’ during typical ‘figure-of-eight’ wing motion. This type of vortical flow may be similar to the Himmelskamp delayed stall effect observed on wide-chord inboard sections of rotor blades [Reference Himmelskamp63].

These unsteady separated flows and vortex dynamics are not yet fully understood but play an important role in the high lift and thrust generation in flapping flight. The present understanding of their flow physics has developed to a large extent from examination of CFD simulations and experimental results, quasi-steady models, and evaluation of relevant vortex force models [Reference Liu, Wang and Liu55, Reference Eldredge and Jones64]. The vortex force map (VFM) method, which will be introduced in Section 5, should be very suitable to contribute to this.

Vortex formation and evolution have had not only a great impact on lift and net thrust generation but because of the variation in these effects for different wing types (Fig. 11) have influenced the wings which have evolved. In the low $Re$ regime the vortex dynamics involved in insect flight is better understood [Reference Ellington, van den Berg, Willmott and Thomas68] while at somewhat higher $Re$ the vortex dynamics for bird and bat flight is less well understood due to the complex morphology and dynamic variation of their wings [Reference Muijres, Johansson, Barfield, Wolf, Spedding and Hedenström58].

Figure 11. (a) Vortex structures in the near-field of hawk moth hovering flight (up) and in the revolving bumble bee wing at the onset of and long after the LEV breakdown. Panel adapted from Liu et al. [Reference Liu, Wang and Liu55]. (b) PIV flow fields of mosquitoes forward flight during the end of pronation, late downstroke, end of supination and late upstroke. Panel adapted from Bomphrey et al. [Reference Bomphrey, Nakata, Phillips and Walker65]. (c) Side view of the vortex wake structure during the transition from hovering to forward acceleration and during forward cruising of a hummingbird hawkmoth (Macroglossum stellatarum). Panel adapted from Yao and Yeo [Reference Yao and Yeo66], (d) Vortical structures around a bat wing identified by the ${\lambda _2}$ -criterion and coloured by the normalised spanwise vorticity. Panel adapted from Wang et al. [Reference Wang, Hang, He and Liu67]. LEV, leading-edge vortex; TEV, trailing-edge vortex; TV, tip vortex; RV, root vortex; DVR, downstroke vortex ring; UVR, upstroke vortex ring.

This high-lift LEV aerodynamic mechanism is observed in flapping flight over a wide range of $Re$ for insects [Reference Ellington, van den Berg, Willmott and Thomas68], birds [Reference Videler, Stamhuis and Povel69, Reference Warrick, Tobalske and Powers70], and bats [Reference Muijres, Johansson, Barfield, Wolf, Spedding and Hedenström58], as shown in Fig. 11. Many insects move their wings with a translation motion that combines with forward flight, flapping up- and down-strokes and rotation, executing figure of eight wing motions. In the translational phase of the wing, a leading-edge vortex (LEV) is developed during the second half of the down-stroke. The LEV grows larger as the flight speed increases, covering the whole wing chord. Dynamic stall (or possibly delayed stall) mechanisms explain the formation and stabilisation of the LEV before it stalls. The LEV formation can also be attributed to the wing’s rotational motion at the end of each stroke, at which point an LEV forms during the flip-over of the wing in preparation for the next half-stroke, and may stay attached to the wing during the downstroke. Alternate shedding of the LEV and TEV at large angles of attack of the wing motion is common in a quasi-2D flow regime. The vortex shedding frequency is proportional to the wing frequency in hovering flight. The characteristics of the vortices vary with the $Re$ , the wing reduced frequency, Strouhal number of the shedding, and the wing kinematics (including wing deformation). In the 3-dimensional flow regime, spanwise flow along the wing stabilises the LEV, which spirals out towards the wingtip in a similar manner to the conical leading-edge vortex on a delta wing, see also delayed stall [Reference Himmelskamp63]. For high aspect-ratio wings a Kelvin-Helmholtz instability occurs on the LEV feeding sheet at relatively high Reynolds numbers, causing the formation of a secondary LEV. Additionally, vortex breakdown (bursting) may occur at 60–70% of the wing span due to a fall in axial velocity in the core. However, both of these phenomena have very little impact on the aerodynamic forces, which are generated [Reference Liu, Wang and Liu55].

During the wing rotation, vortices are also found to augment force production in mainly three ways: the additional circulation produced through the Framer effect which is dependent on the timing and duration of rotation, the clap and fling of paired flapping wings [Reference Weis-Fogh60, Reference Wakeling and Ellington71], and wing-wake interactions including wake capture mechanisms during stroke reversal [Reference Birch and Dickinson72]. This has been suggested as the cause of force production by TEVs in high-frequency high-aspect-ratio mosquito flight [Reference Bomphrey, Nakata, Phillips and Walker65].

In the flapping and gliding flight of larger birds, complex near-body and wake-vortex structures, such as vortex rings are generated for the lift and thrust/drag production. The slotted wing-tip feathers in both gliding and flapping flight generate multiple-core vortices, and vortex spreading, which are responsible for the induced drag reduction [Reference Wu57]. The bird’s V-formation flight also benefits from the energy in the wake vortices. Bats utilise vortices in another way. The tube-like TVs trail in the path of the wing-tip and shed continously throughout the wing stroke, forming incomplete ring-like vortices that travel downstream. A pair of root vortices sometimes can be observed [Reference Muijres, Johansson, Barfield, Wolf, Spedding and Hedenström58] (Fig. 11). The dynamically changing wingspan significantly intensifies the LEV thus augmenting time-averaged lift [Reference Wang, Hang, He and Liu67].

3.3.2 Dynamic stall

Dynamic stall occurs when the effective AoA of an aerofoil section increases rapidly to beyond the mean flow stall angle. Depending on the Reynolds number, wing section shape and rapidity of change of angle, separation can start some way back from the leading edge, even near the trailing edge, then move quickly forward to the leading edge where a strong vortex separates and rolls-up. This vortex induces additional lift increasing it to well above the steady flow value for the current AoA. This is followed by further increase in lift and the development of a strong nose-down pitching moment as the vortex stops growing and is convected over the surface of the aerofoil. When the vortex passes the trailing edge the lift collapses, Fig. 12, and the pitching moment recovers [Reference Leishman and Beddoes73]. The time scales involved are rather similar to those during rapid spoiler deployment but with variations due to the more complex movements of the separation. Because of the important applications for helicopter and wind turbine blades, the phenomenon has been widely investigated, e.g. [Reference McCroskey and Fisher74] and empirical methods of prediction developed [Reference Leishman and Beddoes73, Reference Petot75, Reference Gangwani76]. Flow separation is assumed to start when the normal force coefficient reaches a prescribed value depending on the aerofoil section, followed by LEV growth and downstream convection. The time history of the normal force on the aerofoil is obtained from a combination of linear unsteady lift theory [14] and a separated flow model, often based on Kirchhoff’s formulae, described by first order differential equations with time lags and empirical coefficients. The whole subject is covered extensively in [Reference Leishman77].

Figure 12. Sketch of an aerofoil section lift coefficient ${C_L}$ through a cycle of AoA $\alpha $ for a typical dynamic stall. (a) Flow separates and separating vortex sheet starts to roll up. (b) Vortex convects over upper surface of aerofoil. (c) Vortex passes over trailing edge. (d) Flow is separated over the whole upper surface. (e) Upper surface flow reattaches.

Stages in the process are similar to those of sudden separation from a sharp leading-edge plate discussed earlier. As the LEV grows, back-flow induced on it by its virtual image in the aerofoil approximately balances main flow convection so that the vortex remains close to the leading edge, the lift continues to increase and the centre of the vortex as it grows starts to move very slowly downstream. Eventually back-flow weakens the shedding until LEV growth ceases. At this point streamwise convection dominates, the vortex moves downstream past the trailing edge into the wake and the lift collapses.

3.3.3 3D flows: slender wings in impulsive flow

The formulation for the vortex induced force at an isolated edge can be applied to low aspect ratio wings, AoA $ = \alpha $ small but sufficient to cause leading edge separation, with the LEV distance above the leading edge at any x-section small compared with the wing span at that section. Applying the isolated edge scaling to a cross-section at $x$ of the delta wing, with semi-span $s = \varepsilon x$ where $\varepsilon $ is the semi-apex angle, replacing time $t$ by the equivalent time $x/{U_\infty }$ the concentrated vortex model, Equation (A6) gives:

(36) \begin{align}{C_L} = {2^{2/3}}\pi {\varepsilon ^2}{\left(\frac{\alpha }{\varepsilon }\right)^{5/3}}\end{align}

where ${C_L} = \frac{L}{{\frac{1}{2}\rho U_\infty ^2s}}$ is the lift coefficient for the lift $L$ on a half-wing section.

This is the same to first order in $\frac{\alpha }{\varepsilon }$ as that derived by Brown and Michael [Reference Brown and Michael20] using the concentrated vortex model for the full wing with both leading-edge vortices accounted for in the cross-flow. If the slender wing is subject to a sudden impulsive motion such as heave, the vortex induced forces can be evaluated similarly. LEVs on slender wings are much more stable than those on high aspect ratio wings because of the stabilising effect of the strong axial flow in the core. If the AoA of a slender wing is increased gradually the axial velocity in the vortex core weakens, then reverses locally and vortex ‘bursting’ occurs [Reference Wedermeyer78], initially in the trailing regions of the LEVs downstream of the wing. As the AoA increases further bursting moves forward until it reaches the apex when the wing is then fully stalled. If the rate of change of AoA is faster the development of the vortex burst lags the AoA considerably and exhibits considerable hysteresis if caused by oscillatory AoA [Reference Atta and Rockwell79].

4.1 2D oscillatory flow with zero mean

Oscillatory 2D flows around a body with one or more edges shed strongly rolled up vortex structures. Shedding is typically of one pair of opposite sign vortices propagating outwards from each shedding edge per flow cycle as shown in Fig. 13. The oscillatory flow amplitude ratio parameter is the Keulegan-Carpenter number ${K_c} = \frac{{{V_0}T}}{D}$ , where ${V_0}$ is the velocity amplitude, $T$ is the cycle time and $D$ is a length-scale of the body, usually its diameter. ${K_c}$ is proportional to the ratio of the amplitude of the spatial oscillations of the fluid (or the body) to the diameter of the body. After the oscillatory free stream velocity (the ‘driving flow’) has passed through its maximum during a flow cycle and started decreasing the vortex currently being shed from a given edge continues to grow in strength until the amplitude of the opposing back-flow which it induces at the edge reaches and then exceeds the velocity induced by the free stream. At this point the vorticity being shed from the edge into the rolling-up sheet changes sign. Shortly after this the sheet breaks and opposite signed vorticity commences a new roll-up in the opposite sense to the previous vortex. When the second vortex has reached approximately its maximum (and the first of a new pair has started to form) the two contra-rotating vortices of the first pair convect away together from the edge under their mutually induced velocities in a direction at an acute angle to the edge bisector [Reference Graham36, Reference Bearman81, Reference Obasaju, Bearman and Graham82]. A heuristic argument is suggested to explain why the force induced on the body in small amplitude oscillatory flow which cannot itself convect vortices far from the shedding edge, is dominated by the development of the vortices close to the edge and relatively unaffected by those further away. Provided the strength of one is approximately equal and opposite to the strength of the other in a pair their self-induced path is nearly straight and due to mutually induced convection they rapidly leave the vicinity of the edge. The growth of the vortex moments of each of a nearly equal and opposite strength counter-rotating pair with their effective images in the body, tend to nearly cancel leaving the combined moment approximately constant and hence making little contribution to the force on the body. If on the other hand the two vortices are far from being of equal but opposite strength they follow a more highly curved path and propagation away from the body is restricted. If close to the body surface but far from the edge vortices are convected along the surface by the induced velocity of their images in the surface again preserving the constancy of moment. Therefore in nearly all cases of vortices far from the edge vortex moment does not change appreciably and there is little contribution to the force. For relatively low Keulegan Carpenter numbers, ${K_c} \le 10$ , the vortex force scaling given by the isolated edge theory can be applied.

Figure 13. Vortex shedding from a ${30^{\rm{o}}}$ edge in oscillatory flow. anticlockwise vortices (blue), clockwise vortices (red) [Reference Graham80].

The process of shedding counter-rotating pairs, as shown for a ${30^{\rm{o}}}$ edge in Fig. 13, induces an irregular oscillatory force ${F_v}(t)$ on the shedding body with the same averaged period as the oscillatory free stream but with varying phase difference as in Fig. 15 due mainly to the strength of induced back-flow.

Morison’s equation [Reference Morison, O’Brien, Johnson and Schaaf83],

(37) \begin{align}F(t) = \rho A\frac{{dU}}{{dt}}{C_m} + \frac{1}{2}\rho DU\!\left| U \right|{C_D}\end{align}

is the standard equation used to express the flow direction force $F(t)$ , in this case sectional, per unit transverse length of a body with diameter $D$ and cross-sectional area A in oscillatory flow:

(38) \begin{align}U(t) = {U_0}{\rm{sin}}\!\left( {2\pi t/T} \right)\end{align}

in terms of two cycle averaged coefficients: ${C_m}$ , the inertia coefficient and ${C_D}$ , the drag coefficient.

If $\psi $ is the phase lead of the fundamental of ${F_v}(t)$ with respect to the velocity $U(t)$ then

(39) \begin{align}{C_m} = {C_{m0}} + {C_{mv}}\end{align}

with

(40) \begin{align}{C_{mv}} \propto \frac{{{K_c}}}{{2\pi }}{\rm{sin}}\!\left( \psi \right)\end{align}

and

(41) \begin{align}{C_D} \propto \frac{{3\pi }}{2}{\rm{cos}}\!\left( \psi \right)\end{align}

where the unseparated flow inertia coefficient ${C_{m0}}$ associated with the potential flow induced by the free-stream around the body is the dominant contribution to ${C_m}$ . ${C_{mv}}$ due to the vortex shedding is usually much smaller and less important except when it effects resonant frequencies in cases where the body/fluid density ratio is small, rarely in air but quite often in water [Reference Khalak and Williamson84]. The drag, although often the smaller part of the overall force, can be very important since it determines the vortex damping of the oscillatory motion.

Morison’s equation [Reference Morison, O’Brien, Johnson and Schaaf83] was originally derived to describe the wave-direction force induced on long structural elements of circular cross-section exposed to wave action. Offshore oil and gas jacket structures have many components of this type, subject to large wave forces. Floating and highly compliant structures and ship hulls can perform large oscillatory motions. Baffles and other plates may be used to increase damping (e.g. heave plates on spar buoys, bilge keels on ship hulls). Many fundamental studies of sectional forces induced on bluff bodies transverse to planar oscillatory flow have been carried out to establish values of the force coefficients, particularly for circular, square and flat plate sections [Reference Sarpkaya85, Reference Singh86] rounded square pontoon sections [Reference Bearman, Graham, Obasaju and Drossopoulos87] and ship hull cross-sections [Reference Downie, Bearman and Graham88, Reference Cozens89].

For the common symmetric sections with typically just two separation points, one either side, the convection paths of shed vortex pairs may be symmetric or antisymmetric with respect to the flow direction and often, apparently randomly, flip between the different states which may also exist simultaneously at different transverse sections of a long body [Reference Obasaju, Bearman and Graham82]. Spanwise (lengthwise) mode jumping is a frequent occurrence in more complex oscillatory cylinder flows. Since the vortices which are shed as counter-rotating pairs in oscillatory flow contribute little to the dipole force moment provided the value of ${K_c}$ is not too large the isolated edge scaling analysis of Section 2.3 should apply throughout the cycles of flow. Therefore changing the characteristic time to the cycle time $T$ , gives the following dependence of the Morison force coefficients on ${K_c}$ [Reference Graham36]:

(42) \begin{align}{C_D} = AK_c^{\frac{{3 - 2\kappa }}{{2\kappa - 1}}}\end{align}

(43) \begin{align}{C_m} = {C_{m0}} + BK_c^{\frac{2}{{2\kappa - 1}}}\end{align}

where $\kappa = 2 - \delta /\pi $ , is the included angle of the shedding edge and A and B are dimensionless coefficients which depend on the body section. Viscous effects are more usually expressed for these flows by the Stokes parameter $\beta = \frac{{{D^2}}}{{\nu T}} = \frac{{Re}}{{{K_c}}}$ than the Reynolds number $Re$ . For the three most common edge angles: $\delta \simeq 0$ (plates and bilge keels),

(44) \begin{align}{C_D} = AK_c^{ - 1/3},{C_m} = {C_{m0}} + BK_c^{2/3},\end{align}

$\delta = \pi /2$ (box sections, ship and pontoon hulls)

(45) \begin{align}{C_D} = A,{C_m} = {C_{m0}} + B{K_c},\end{align}

and $\delta = \pi $ referring to body sections without sharp or high curvature edges such as circular cylinder. Although for separation from continuous surfaces with mobility of the separation point the analysis is less justifiable, nevertheless the scaling appears to be reflected in the measured data [Reference Graham36]:

(46) \begin{align}{C_D} = A{K_c},{C_m} = {C_{m0}} + BK_c^2\end{align}

These predictions for ${C_D}$ as powers of ${K_c}$ are compared with experimental data in Fig. 14 showing good agreement. Normalising the force cycle by the above powers of Keulegan-Carpenter number for the three sections, flat plate, diagonally aligned square and circular section, similarly gives a very good collapse of the data for force cycles over the low ${K_c}$ range as shown in Fig. 15. They also predict well the vortex shedding generated by small amplitude roll motion from the bilges of ship hulls with and without bilge keels [Reference Downie, Bearman and Graham88, Reference Cozens89]. Rectangular hulls are predicted by this scaling to have a constant Morison drag coefficient, independent of flow amplitude (i.e. ${K_c}$ ), which is the usual assumption (damping force proportional to roll velocity squared) made by the industry, see e.g. [Reference Mathisen and Price90]. But it should be noted that the result for thin plates predicts an increasingly large drag coefficient as the amplitude of motion decreases which supports the effectiveness of bilge keels and baffle plates as damping devices. In the usual ’mixed’ cases where right-angle bilges are either rounded at the corner or have a bilge-keel added there (or both), the behaviour of the drag coefficient follows the expected behaviour for the local geometry, rounding or plate, at small amplitude of oscillation but reverts to the constant ${C_D}$ as the amplitude becomes large [Reference Xi, Milne, Draper and Chen91, Reference Xi and Milne92].

Figure 14. Drag coefficient at low ${K_c}$ for $ + $ Normal Flat plate, $\diamondsuit $ diamond square, $o$ circular cylinder, [Reference Sarpkaya85, Reference Xi, Milne, Draper and Chen91].

Figure 15. Normalised streamwise force versus time over one cycle [Reference Pitt Ford and Babinsky39].

The circular section is the one which has been most widely studied because it is the section used in so many structures. For very small amplitudes of oscillatory motion, ${K_c} \lt 1$ , the flow is unseparated but boundary-layer displacement and skin friction generate a small additional streamwise force. These boundarylayers also have a spanwise instability which develops into alternating Honji vortices running round the cylinder [Reference Honji93]. For ${K_c} \gt \sim 3$ separation occurs and vortices shed in pairs antisymmetrically and obliquely to the flow direction, Fig. 16. At medium and large values of the Keulegan Carpenter number the flows exhibit spanwise correlated shedding regimes similar to those for sharp sections with two fixed shedding edges [Reference Bearman81]. The process exhibits periods of departure from these states, single vortices circling around the body, flows showing abrupt jumps between alternative states in the spanwise direction [Reference Obasaju, Bearman and Graham82]. Averaged over many cycles the shedding process for a symmetric section doesn’t induce any net mean force but the irregularity of the strong cyclic in-line force can produce non-zero averages over a few flow cycles. Because the shedding does not form a conventional streamwise wake there is no strong periodic transverse force generated in this ${K_c}$ range.

Figure 16. Vortex shedding from a circular cylinder at ${K_c} = 4.0$ . DVM simulation with passive particle visualisation [Reference Lin, Bearman and Graham94].

At high Keulegan Carpenter numbers, $\left( { \geqslant 10} \right)$ , the effect of the free-stream $U(t) = {U_0}{\rm{sin}}\frac{{2\pi t}}{T}$ on the flow tends towards quasi-steady behaviour [Reference Sarpkaya95]. The in-line force in Morison’s equation is dominated by drag as the ratio of $D\frac{{dV}}{{dt}}$ to ${V^2}$ decreases and the drag coefficient tends towards its value in steady flow. A periodic transverse force develops together with the formation of a finite length ‘von Karman type’ vortex wake containing an antisymmetric array of vortices downstream of the body over each half-cycle of the free-stream. Remarkably, even when these wakes contain only very few vortices the principle of a constant Strouhal number and hence constant streamwise vortex spacing, defines the shedding process and wake, e.g. [Reference Verley96]. A frequency modulated transverse force with constant (lift) coefficient ${C_L}$ and Strouhal number $S$ follows closely the equation for the time dependent phase $\phi (t)$ :

(47) \begin{align}\phi = \frac{S}{D}\int U(t)dt\end{align}

where $D$ is the cylinder diameter and $S$ is its Strouhal number. Therefore

(48) \begin{align}L(t) = \frac{1}{2}U_0^2D{\rm{si}}{{\rm{n}}^2}\frac{{2\pi t}}{T}{C_L}cos\!\left( {\frac{{{K_c}S}}{{2\pi }}{\rm{cos}}\frac{{2\pi t}}{T} + {\phi _0}} \right)\end{align}

where ${\phi _0}$ is the starting phase of the shedding in each half-cycle. Figure 17 shows the excellent agreement between measured data for a circular cylinder in oscillatory flow with ${K_c} = 43.4$ and the quasi-steady Equation (48) with a standard value of Strouhal number, $S = 0.20$ for steady incident flow at similar Reynolds number. In this case there are only about five vortices altogether in the wake. Surprisingly the new wake forming in each half-cycle appears to accommodate the previous returning wake without disruption. Results for Keulegan Carpenter numbers down to ${K_c} \le 20$ continue to show a strong measure of agreement with this theory [Reference Bearman, Graham and Obasaju97].

Figure 17. Transverse force coefficient for a circular cylinder over a half-cycle of oscillatory incident flow: ${K_c} = 43.4$ , $S = 0.2$ , ${\hat C_L} = 1.33$ , $\psi = {351.0^{\rm{O}}}$ , x experimental data, - quasi-steady theory Equation (48), from Ref. [Reference Bearman, Graham and Obasaju97].

4.2 Oscillatory flows with non-zero mean

4.2.2 3-Dimensional flows: vortex induced forces on slender bodies in oscillatory cross-flow with forward motion

In the case of roll motion of slender wing aircraft and the calculation of roll damping the parameter ${U_\infty }T/c$ is small and application of quasi-steady analysis is reasonable. However, in the case of ship hulls, this assumption is less certain. Ikeda et al. [Reference Ikeda, Himeno and Tanaka100] assessed the effect of forward motion on eddy making damping on a rolling hull by experimental subtraction of other contributions to the damping, and concluded that it led to a decrease in the damping coefficient. On the other hand Al Hukail [Reference Al-Hukail, Bearman, Downie, Graham and Zhao101] carried out computations of vortex shedding from a rolling slender rectangular block section hull. Some vortex pairing for each cycle of roll was observed in the computer simulations. An example is shown in Fig. 18. It was found in these computations that the forward motion led to a significant increase in the roll damping.

Figure 18. Vortices shed from a rolling ship hull of rectangular cross-section in forward motion, from Ref. [Reference Ramesh, Gopalarathnam, Granlund, Ol and Edwards102].

5.1 Element-based methods

The general expressions for the force and moment on a rigid body within a general framework and a unified perspective, which involve the velocity field instead of the pressure at the body surface was first proposed by Quartapelle and Napolitano [Reference Quartapelle and Napolitano105]. Chang [Reference Chang106] expressed the force acting on a body as the scalar product of a potential flow and the Lamb vector. Howe [Reference Howe107] expresses the force and moment exerted on a rigid body in unsteady motion in a uniform incompressible flow in terms of the velocity and vorticity in the flow and the body velocity. These methods together with the vortex force/moment map method are intrinsically the same.

In the element-based theory, the Lamb vector depends on the choice of frame of reference, thus the projection of Lamb’s vector, i.e. the elementary force contribution, is not Galilean invariant [Reference Gao and Wu108]. A Galilean invariant representation of the force is important, especially in those arbitrary multiple-body dynamic problems where no universal frame of reference exists. To address this, a modified force-element theory was proposed by Gao et al. [Reference Gao, Zou, Shi and Wu109] based on the weighted integral of the second invariant of the velocity gradient tensor $Q = \frac{1}{{2\rho }}{\nabla ^2}p$ . This modified force-element method satisfies Galilean invariance. It has been applied in analysing the 3D force of low Reynolds number flows past a plunging aerofoil [Reference Gao, Cantwell, Son and Sherwin110], where a curve of the accumulated volume fraction with respect to the downstream position was drawn to quantify the contribution of local structures on the force. More applications of the modified force-element theory include the jellyfish flow [Reference Kang, Gao, Han, Cui and Lu111], where a negative $Q$ is also shown to be the predominant propulsion mechanism and a self-propelled flexible plate in a uniform shear flow where the background vorticity is non zero [Reference Luo, Gao and Lu112]. It has been further explored by Menon and Mittal [Reference Menon and Mittal113], and the derivation of the equations for the moment has been given. The element-based methods have also been explained as the virtual power principle in fluid dynamics by Yu [Reference Yu114], where the integral term in the equation needs to satisfy the Galilean invariance requirement.

5.2 VFM for 2D inviscid flow

The concept of the 2D method was first proposed by Li and Wu [Reference Li and Wu115] in studying the flow around a flat plate at a moderate AoA ( $ \lt {20^{\rm{o}}}$ ). It amends the inviscid Wagner lift model for starting flow to account for the influence of separated leading edge and trailing edge vortices. The separated leading-edge vortices are always presented as discrete vortices released sequentially from the leading edge. Assuming the trailing edge vortices as either a planar vortex sheet or a series of discrete vortices released sequentially gives two sets of lift formulae. Then the lift models are coupled with a discrete vortex simulation with the results validated against numerical computation by CFD. The vortex system of the starting flow around a flat plate at moderate AoA with an arbitrary number of vortices shed from the leading and trailing edges, plus a quasi-planar vortex sheet, (VS, which may also be represented by discrete vortices) is shown in Fig. 19(a). The vortex system of the starting flow around a flat plate at arbitrarily high AoA is shown in Fig. 19(b). An explicit lift expression in terms of modified Wagner functions and the properties of the separated vortices has been derived.

Figure 19. Schematic representation of the Vortex flow around a flat plate [Reference Li and Wu115, Reference Li and Wu116] at (a) moderate and (b) high AoA. LEV, leading edge vortex; TEV, trailing edge vortex; VS, vortex sheet; RVS, residual vortex sheet.

The non-dimensional time is defined by $\tau = {V_\infty }t/c$ (equivalent to the number of chord lengths traveled by the aerofoil at constant speed in a ground-based frame of reference). In the vortex structures shown in Fig. 19, the vortex sheets are presented by $C({s,t})$ and have the strength of $k({s,t})$ , where $s$ is the location and $t$ is the physical time, the separated discrete vortices with strength ${\rm{\Gamma }}_v^{(i)}$ located at $\left( {x_v^{(i)},y_v^{(i)}} \right)$ move with the local speed $\boldsymbol{U}_v^{(i)} = \left( {dx_v^{(i)}/dt,dy_v^{(i)}/dt} \right)$ . In these cases, the body speed once started is constant, thus the added mass force is zero, and the impulse Equation (18) takes the form of the Phillips equation [Reference Phillips42]. Here, the vorticity consists of three parts, the bound vorticity ${{\rm{\Gamma }}_b}$ , the point vortices ${\rm{\Gamma }}_v^{(i)}$ , and the vortex sheet $k({s,t})$ , and they satisfy the Kelvin theorem, Equation (49)

(49) \begin{align}{{\rm{\Gamma }}_b} + \mathop \sum \limits_{i = 1}^I {\rm{\Gamma }}_v^{(i)} + \mathop \int \nolimits_C k({s,t})ds = 0\end{align}

Accordingly, the total, i.e., normal force formula (there is no leading-edge suction) can be divided into three parts and using Equation (49) we have

(50) \begin{align}\begin{array}{l}F = {F_b} + {F_{PV}} + {F_{VS}}\\[3pt]\left\{ {\begin{array}{*{20}{l}}{{F_b} = - \rho {V_\infty }{{\rm{\Gamma }}_b} + \rho {{\dot{\Theta}}_b}}\\[3pt]{{F_{PV}} = \rho \mathop \sum \limits_{i = 1}^I \left( { - {V_\infty }{\rm{\Gamma }}_v^{(i)} + \dfrac{{d{\rm{\Gamma }}_v^{(i)}{x_i}}}{{dt}}} \right)}\\[9pt]{{F_{VS}} = - \rho {V_\infty }\mathop \int \nolimits_{\rm{C}} k({s,t})ds + \rho \dfrac{d}{{dt}}\mathop \int \nolimits_C k({s,t})C({s,t})ds}\end{array}} \right.\\[3pt]Or\\[3pt]\left\{ {\begin{array}{*{20}{l}}{{F_b} = \rho {{\dot{\Theta}}_b}}\\[5pt]{{F_{PV}} = \rho \mathop \sum \limits_{i = 1}^I \dfrac{{d{\rm{\Gamma }}_v^{(i)}{x_i}}}{{dt}}}\\[9pt]{{F_{VS}} = \rho \dfrac{d}{{dt}}\mathop \int \nolimits_C k({s,t})C({s,t})ds}\end{array}} \right.\end{array}\end{align}

Here ${{\rm{\Gamma }}_b} = \mathop \int \nolimits_0^c {\gamma _b}\!\left( {x,t} \right)dx$ is the bound vorticity and ${{\rm{\Theta }}_b} = \mathop \int \nolimits_0^c {\gamma _b}\!\left( {x,t} \right)xdx$ is the first moment of the bound vorticity. For a flat plate or a thin aerofoil, the bound vorticity ${\gamma _b}\!\left( {x,t} \right)$ can be obtained by thin aerofoil theory with the added separated vortices [Reference Bisplinghoff, Ashley and Halfman1].

(51) \begin{align}\frac{1}{{2\pi }}\mathop \int \nolimits_0^c \frac{{{\gamma _b}\!\left( \xi \right)}}{{x - \xi }}d\xi = {v^{\left( g \right)}} + {v^{\left( {PV} \right)}} + {v^{\left( {VS} \right)}}\end{align}

Here ${v^{\left( g \right)}} = - {V_\infty }sin\alpha $ , ${v^{\left( {PV} \right)}}$ and ${v^{\left( {VS} \right)}}$ are given by the Biot–Savart relation [Reference Li and Wu115]. Accordingly, the bound vorticity can be decomposed into three parts ${\gamma _b}\left( {x,t} \right) = \gamma _b^{\left( g \right)}\left( {x,t} \right) + \gamma _b^{\left( {PV} \right)}\left( {x,t} \right) + \gamma _b^{\left( {VS} \right)}\left( {x,t} \right)$ , ${{\rm{\Gamma }}_b} = {\rm{\Gamma }}_b^{\left( g \right)} + {\rm{\Gamma }}_b^{\left( {PV} \right)} + {\rm{\Gamma }}_b^{\left( {VS} \right)}$ , and ${{\rm{\Theta }}_b} = {\rm{\Theta }}_b^{\left( g \right)} + {\rm{\Theta }}_b^{\left( {PV} \right)} + {\rm{\Theta }}_b^{\left( {VS} \right)}$ . Where ${\rm{\Gamma }}_b^{\left( g \right)} = - \pi c{V_\infty }sin\alpha $ . For the flow around a moderate AoA flat plate in Fig. 19(a), by solving Equation (51) under assumption of planar trailing edge vortex sheets and $sin\alpha \approx \alpha $ , the lift formula can be obtained by projection of the normal force (50) in the lift direction [Reference Li and Wu115]

(52) \begin{align}L = {\phi _A}{L_\infty } + \rho \mathop \sum \limits_{i = 1}^I \left( {{O_A}{V_\infty } + {{\boldsymbol{\Lambda }}_A} \cdot \boldsymbol{U}_v^{(i)}} \right){\rm{\Gamma }}_v^{(i)}\end{align}

where ${\phi _A}$ is the modified Wagner function, ${O_A}$ and ${{\boldsymbol{\Lambda }}_A}$ are problem-independent factors [Reference Li and Wu115]. ${L_\infty } = \rho \pi cV_\infty ^2{\rm{sin}}\alpha $ is the steady asymptotic lift which can be obtained from the Kutta–Joukowski theorem. Equation (52) decomposes the lift as a modified Wagner lift with separated vortex contributions. It can be proven that giving the strength of a self-contained stationary vortex above a flat plate, it is precisely the same as Saffman’s lift formula [Reference Saffman117].

For the flow around a high AoA flat plate in Fig. 19(a), by solving Equation (51) under the assumption of vanishing RVS and substituting the solution into Equation (52), the force formula for a flat plate can be reorganised as a sum of the vortex movement ( $vm$ ) term and the circulation production $\left( {cp} \right)$ term [Reference Li and Wu116]

(53) \begin{align}\begin{array}{*{20}{l}}{F = {F^{\left( {vm} \right)}} + {F^{\left( {cp} \right)}}}\\[4pt]{\left\{ {\begin{array}{*{20}{l}}{{F^{\left( {vm} \right)}} = \rho \mathop \sum \limits_{i = 1}^I \left( {{{\boldsymbol{\Lambda }}_i} \cdot \boldsymbol{U}_v^{(i)}} \right){\rm{\Gamma }}_v^{(i)}}\\[10pt]{{F^{\left( {cp} \right)}} = \dfrac{1}{2}\rho c\!\left( {{R_{LE}}{{\dot{\Gamma}}_{LE}} + {R_{TE}}{{\dot{\Gamma}}_{TE}}} \right)}\end{array}} \right.}\end{array}\end{align}

where ${{\boldsymbol{\Lambda }}_i}$ , ${R_{LE}}$ , and ${R_{TE}}$ are problem-independent factors [Reference Li and Wu116]. It is proved in Li and Wu [Reference Li and Wu116] that, ${F^{\left( {cp} \right)}}$ which is due to the newly shed element of vortex sheet, is small relative to ${F^{\left( {vm} \right)}}$ , tending to zero as the time step tends to zero.

In Equations (52) and (53), the flow field ( $\boldsymbol{U}_v^{(i)}$ and ${\rm{\Gamma }}_v^{(i)}$ ) can be obtained either from experiment (e.g., PIV) or from numerical methods (e.g. discrete vortex simulation or CFD). The factors ${O_A}$ , ${{\boldsymbol{\Lambda }}_A}$ , ${{\boldsymbol{\Lambda }}_i}$ , ${R_{LE}}$ , and ${R_{TE}}$ can be pre-obtained independently once the geometry of the body is given and used to draw the vortex force (line) maps. These explicitly indicate the lift-enhancing and -reducing directions, and the effect of the sign and position of vortices on the force.

The vortex force line map for the vector ${{\boldsymbol{\Lambda }}_i}$ in Equation (53) is shown in Fig. 20. The explanation and application of this can be found in Ref. [Reference Li and Wu116].

Figure 20. (a) Vortex force line map, panels (b) and (c) are the amplification of zones near the leading and trailing edges [Reference Li and Wu116].

5.3 Vortex force map for 2D viscous flows

5.3.1 Single-body flow vortex force map

Li and Wu [Reference Li and Wu118] adopted Howe’s integral force formulae [Reference Howe107] to extend the above inviscid vortex force line map approach for the normal force on a flat plate to viscous flow around ^arbitrary 2D bodies. The normal and axial force components for 2D flow are expressed as

(54) \begin{align}\left\{ {\begin{array}{*{20}{l}}{{F_N} = \rho \mathop \int\!\!\!\int \nolimits_{\rm{\Omega }} {\boldsymbol{\Lambda }}_N^S \cdot \boldsymbol{U}{\omega _z}d{\rm{\Omega }}}\\[7pt]{{F_A} = \rho \mathop \int\!\!\!\int \nolimits_{\rm{\Omega }} {\boldsymbol{\Lambda }}_A^S \cdot \boldsymbol{U}{\omega _z}d{\rm{\Omega }}}\end{array}} \right.\end{align}

Here ${\boldsymbol{\Lambda }}_{\boldsymbol{N}}^{\boldsymbol{S}} = \left( {\partial {\phi _{NS}}/\partial y, - \partial \phi _N^S/\partial x} \right)$ and ${\boldsymbol{\Lambda }}_{\boldsymbol{A}}^{\boldsymbol{S}} = \left( {\partial \phi _A^S/\partial y, - \partial \phi _A^S/\partial x} \right)$ are standard VFM factors for normal and axial forces respectively. Note that the integral over the whole fluid region ${\rm{\Omega }}$ should include all shedding vortices in viscous flow, and it replaces the summation of the potential vortices in Equations (52) and (53). $\phi _N^S$ and $\phi _A^S$ are the standard hypothetical potentials for a fixed body with unit inflow from infinity in the normal and axial directions respectively, as defined in Li and Wu [Reference Li and Wu118].

The lift and drag can be obtained by a combination of the normal and axial forces in Equation (54). Note that, when $({x,y}) \to \infty $ , ${\boldsymbol{\Lambda }}_{\boldsymbol{N}}^{\boldsymbol{S}}$ and ${\boldsymbol{\Lambda }}_{\boldsymbol{A}}^{\boldsymbol{S}}$ are non-zero, which means that vortices far away from the body have a significant impact on the force. This is physically unlikely and inconvenient practically as correction terms should be introduced when applying to truncated domains [Reference Li and Wu118]. Thus, the force formula where only near-body vortices have important effects through coordinate transformations and Howe’s formulae [Reference Howe107] was derived [Reference Li, Wang, Graham and Zhao119] as

(55) \begin{align}\left\{ {\begin{array}{*{20}{l}}{{F_k} = F_k^{\left( {dd} \right)} + F_k^{\left( {vp} \right)} + F_k^{\left( {vis - p} \right)} + F_k^{\left( {vis - f} \right)}}\\[4pt]{F_k^{\left( {add} \right)} = \mathop \oint \nolimits_{{l_B}} {\phi _k}\frac{{d{\boldsymbol{U}_{\boldsymbol{B}}}}}{{dt}} \cdot {\boldsymbol{n}}dl}\\[4pt]{F_k^{\left( {vp} \right)} = \rho \mathop \int\!\!\!\int \nolimits_{\rm{\Omega }} {{\boldsymbol{\Lambda }}_k} \cdot {\boldsymbol{U}}{\omega _z}d{\rm{\Omega }}}\\[4pt]{\left. {\begin{array}{*{20}{l}}{F_k^{\left( {vis - p} \right)} = \mu \mathop \oint \nolimits_{{l_B}} {\omega _z}d{\phi _k}}\\[4pt]{F_k^{\left( {vis - f} \right)} = \mu \mathop \oint \nolimits_{{l_B}} {\omega _z}{\boldsymbol{k}} \cdot {\boldsymbol{dl}}}\end{array}} \right\}F_k^{\left( {correction} \right)}}\end{array}} \right.\end{align}

where ${\boldsymbol{U}_{\boldsymbol{B}}} = {\boldsymbol{U}_{\boldsymbol{B}}}\!\left( {{x_B},{y_B},t} \right)$ is the velocity of the body surface ( ${l_B}$ ),

(56) \begin{align}{{\boldsymbol{\Lambda }}_{\boldsymbol{k}}} = \left( {\partial {\phi _k}/\partial y, - \partial {\phi _k}/\partial x} \right)\end{align}

is the vortex pressure force vector, and ${\phi _k}$ is the hypothetical potential defined as the velocity potential for hypothetical fluid motion induced by the translational motion of the body surface at unit speed in the negative $k$ th direction. According to the definition of ${\phi _k}$ , we have

(57) \begin{align}\begin{array}{*{20}{l}}{\left\{ {\begin{array}{*{20}{l}}{{\nabla ^2}{\phi _k} = 0}\\[4pt]{ - \nabla {\phi _k} \cdot {\boldsymbol{n}_{\boldsymbol{B}}} = \boldsymbol{k} \cdot {\boldsymbol{n}_{\boldsymbol{B}}} = {n_{B,k}}}\\[4pt]{\nabla {\phi _k} = 0{\rm{\;\;\;\;\;\;\;\;}}({x,y}) \to \infty }\end{array}\!,} \right.}\end{array}\end{align}

where ${\boldsymbol{n}_B}$ is the normal vector pointing inward from the body surface. Note that these hypothetical potentials are different from the standard vortex force factors $\phi _N^S$ and $\phi _A^S$ . ${\phi _k}$ ensures that only near-body vortices have an impact on the vortex-induced body force. Equation (55) includes the added mass force $F_k^{\left( {add} \right)}$ caused by body motion, the viscous pressure force $F_k^{\left( {vis - p} \right)}$ caused by vorticity diffusion on the body surface, and the skin friction force $F_k^{\left( {vis - f} \right)}$ caused by the shear stress. The unsteady effects are included in the force formula implicitlythrough time-dependent $\omega $ and $\boldsymbol{U}$ . Since time $t$ does not occur explicitly in the force formula, temporal resolution of the flowfield data does not affect the accuracy of the method and snapshots of the field can be used. However, it is essential to precisely charact erise the body acceleration for an accurate prediction of the added mass force.

The Laplace equations and the corresponding boundary conditions (57) with $k$ in the lift and drag directions are solved to give the potential ${\phi _L}$ and ${\phi _L}D$ , thus ${{\rm{\Lambda }}_L}$ and ${{\rm{\Lambda }}_D}$ [Reference Li, Wang, Graham and Zhao119]. Analytical solutions of (57) exist for simple bodies, such as a circular cylinder (Fig. 21(a) and (b)). For a general shape of body, such as a NACA0012 aerofoil at $\alpha = {45^o}$ , (57) can be solved numerically, and the results are plotted in Fig. 21(c) and (d). Figure 21(a)–(d) are called the vortex pressure force maps, where each line with arrows defines the direction for a counterclockwise vortex (having positive strength) to give a positive lift/drag or the direction for a clockwise vortex (having negative strength) to give a negative lift/drag. These VFMs can be employed to evaluate the force exerted by a specific vortex and to extract forces from the information of the flow field enclosing the body [Reference Li, Wang, Graham and Zhao119]. E.g., for a circular cylinder, the symmetrical vortex pressure lift map about the $y = 0$ axis explains that a pair of counter-rotating vortices moving symmetrically downstream of a cylinder contribute zero lift. The vortex force approach has been validated against CFD results of NACA0012 at several angles of attack and Reynolds numbers by Li et al. [Reference Li, Wang, Graham and Zhao119]. The dominant force is the vortex pressure force. The viscous forces are negligible when the Reynolds number is large enough. Moreover, the computation of vortex forces was demonstrated using coarsely sampled data on a small, truncated domain, and it was found that the method produces acceptable accuracy, which facilitates its application to computing body forces from PIV data.

Figure 21. Vortex pressure force maps for the (a) lift of a cylinder, (b) drag of a cylinder, (c) lift of NACA0012 aerofoil at ${\rm{\alpha }} = {45^o}$ , and (d) drag of NACA0012 aerofoil at $\alpha = {45^{\rm{o}}}$ [Reference Wang, Zhao, Graham and Li120].

5.3.2 Multi-body flow vortex force map

We consider a number $M$ of solid bodies, relatively stationary to each other, in the viscous flows of constant density $\rho $ and viscosity $\mu $ in a control volume ${\rm{\Omega }}$ bounded by ${S_\infty }$ at infinity. Each body has a volume ${{\rm{\Omega }}_{mB}}$ $\left( {m = 1,2, \ldots, M} \right)$ , bounded by closed surfaces ${S_{mB}}$ $\left( {m = 1,2, \ldots, M} \right)$ (degrade to ${l_{mB}}$ in a 2D case). The free-stream velocity is ${V_\infty }$ incident at an angle $\alpha $ to the $i$ th-body axis. The force ${\boldsymbol{F}_{\boldsymbol{i}}}$ acting on the $ith$ -body has been derived rigorously in Wang et al. [Reference Wang, Zhao, Graham and Li120] in the body-fixed frame $({x,y,z})$ of the $i{\rm{th}}$ -body. The schematic of the flow and the force components are shown in Fig. 22(a). Similar to a single body, the $kth$ component of the force (with a unit vector $\boldsymbol{k}$ ) of the $i{\rm{th}}$ -body is decomposed into the added mass force $F_{ik}^{\left( {Add} \right)}$ , the vortex-pressure force $F_{ik}^{\left( {vor - P} \right)}$ , the viscous-pressure force $F_{ik}^{\left( {vis - P} \right)}$ and the skin-friction force $F_{ik}^{\left( {friction} \right)}$ , and the first three make up the pressure force $F_{ik}^{\left( {pressure} \right)}$ :

(58a) \begin{align}{F_{{{\rm{\;}}^{ik}}}} = \underbrace {F_{ik}^{\left( {Add} \right)} + F_{ik}^{\left( {vor - P} \right)} + F_{ik}^{\left( {vis - P} \right)}}_{F_{ik}^{\left( {pressure} \right)}} + F_{ik}^{\left( {friction} \right)}\end{align}

(58b) \begin{align}F_{ik}^{\left( {Add} \right)} = - \rho \mathop \sum \limits_{i = 1,2, \ldots, M} \mathop \int\!\!\!\int \nolimits_{{l_{iB}}} {\phi _{ik}}\frac{{d{\boldsymbol{U}}}}{{dt}} \cdot {\boldsymbol{n}}dl\end{align}

(58c) \begin{align}F_{ik}^{\left( {vor - P} \right)} = \rho \mathop \int\!\!\!\int \nolimits_{\rm{\Omega }} {{\boldsymbol{\Lambda }}_{ik}} \cdot {\boldsymbol{U}}{\omega _z}d{\rm{\Omega }}\end{align}

(58d) \begin{align}F_{ik}^{\left( {vis - P} \right)} = \mu \mathop \sum \limits_{i = 1,2, \ldots, M} \mathop \oint \nolimits_{{l_{iB}}} {\omega _z}d{\phi _{ik}}\end{align}

(58e) \begin{align}F_i^{\left( {friction} \right)} = \mu \mathop \oint \nolimits_{{l_{iB}}} {\omega _z}{\boldsymbol{k}} \cdot {\boldsymbol{dl}}\end{align}

where the vortex-pressure force factor is expressed as

(59) \begin{align}{{\boldsymbol{\Lambda }}_{ik}} = \left( {\frac{{\partial {\phi _{ik}}}}{{\partial y}}, - \frac{{\partial {\phi _{ik}}}}{{\partial x}}} \right).\end{align}

where ${\phi _{ik}}$ is defined as

(60) \begin{align}\left\{ {\begin{array}{*{20}{l}}{{\nabla ^2}{\phi _{ik}} = 0}\\[4pt]{ - \nabla {\phi _{ik}} \cdot {{\boldsymbol{n}}_{{\boldsymbol{iB}}}} = {\boldsymbol{k}} \cdot {{\boldsymbol{n}}_{{\boldsymbol{iB}}}},{\rm{\;\;\;\;}}({x,y,z}) \to {l_{iB}}}\\[4pt]{ - \nabla {\phi _{ik}} \cdot {{\boldsymbol{n}}_{{\boldsymbol{mB}}}} = 0,{\rm{\;\;\;\;}}({x,y,z}) \to {l_{mB}} \cup {l_\infty }\!\left( {m \ne i} \right)}\end{array}} \right.\end{align}

Figure 22. The 3D method for the multi-body configuration. (a) The schematic of the flow around multi-body configuration. (b) The 3D for the main aerofoil and for the flap in a wing-flap configuration. (c) The lift coefficient on the main aerofoil and on the flap in the wing-flap configuration at $\alpha = {45^{\rm{o}}}$ . (d) Contours of vortex lift distribution displaying lift coefficients acting on the main aerofoil (left column) and on the flap (right column) contributed by local vortices, for a wing-flap configuration starting flow at $\alpha = {45^{\rm{o}}}$ and at instants $\tau $ = 7 and 8, with streamlines [Reference Li, Wang, Graham and Zhao121].

Here ${\boldsymbol{n}_{\boldsymbol{iB}}}$ and ${\boldsymbol{n}_{\boldsymbol{mB}}}$ are the normal vectors pointing inwards from each body surface. ${l_{iB}}$ and ${l_{mB}}$ denote 2D body surfaces. Taking lift as an example, solving the Laplace Equation (60) in the direction of $\boldsymbol{k} = {\boldsymbol{k}_{\boldsymbol{L}}} = \left({-}{\rm{sin}}\alpha, {\rm{cos}}\alpha \right)$ for $i = 1$ and $i = 2$ in the wing-flap configuration, one can obtain the vortex force vectors for the main aerofoil ( ${{\boldsymbol{\Lambda }}_{1L}}$ ) and for the flap ( ${{\boldsymbol{\Lambda }}_{2L}}$ ) respectively. The 3D for the main aerofoil and for the flap are shown in Fig. 22(b). The concept of vortex pressure maps still exists for individual bodies when the flow field contains multiple bodies and should not be understood as a superposition of the individual contributions. The maps separate the body forces and provide links to the flow features (velocity and vorticity).

Vortex lift evolution for the main aerofoil and for the flap in an impulsively started wing-flap configuration with zero flap deflection angle and for $\alpha = {45^{\rm{o}}}$ at Re = 1,000 are shown in Fig. 22(c). The total force and other force components are also shown. Figure 22(d) is the spatial distribution of local vortex pressure lift due to the local vorticity at typical instants. The pressure force contributed by vorticity in the flow field was found to be the dominant force and its variation can be explained by the motion of the vortices relative to each body.

Figure 23. The vortex moment map for the NACA0012 aerofoil with a pitching axis (a) LE, (b) $c/2$ , (c) $c/4$ , (d) TE. (e) The comparison between the pitching moment obtained by the VMM method and CFD with a pitching axis $c/4$ at $\alpha = {60^o}$ and $Re = 1,000$ . (f) Contours of the vortex moment coefficient per unit area (left) and the vorticity (right) at a typical instant $\tau = 2.0$ with streamlines [Reference Li, Wang, Graham and Zhao121].

5.3.3 Single-body flow vortex moment map

A vortex moment map (VMM) method has also been derived to predict the pitching moment on a single body [Reference Wang, Zhao, Graham and Li120]:

(61) \begin{align}{M_p} = \rho \zeta \mathop \int\!\!\!\int \nolimits_{\rm{\Omega }} {\textsf{F}_{\boldsymbol{p}}} \cdot \boldsymbol{U}{\omega _z}d{\rm{\Omega }}.\end{align}

where the VMM vector

(62) \begin{align}{\textsf{F}_{\boldsymbol{p}}} = \frac{1}{\zeta }\left( {\frac{{\partial {\chi _p}}}{{\partial y}}, - \frac{{\partial {\chi _p}}}{{\partial x}}} \right)\end{align}

is normalised by a unit angular velocity multiplied by the characteristic length. The flow-independent hypothetical potential ${\chi _p}$ is defined as the velocity potential for hypothetical fluid motion induced by the rotation of the body ${{\rm{\Omega }}_B}$ at unit angular velocity about an axis that passes through the reference point ${{\boldsymbol{x}}_{\boldsymbol{p}}}$ and perpendicular to the coordinate plane, satisfying the following Laplace equation and boundary conditions

(63) \begin{align}\left\{ {\begin{array}{*{20}{l}}{{\nabla ^2}{\chi _p} = 0}\\[9pt]{\dfrac{{\partial {\chi _p}}}{{\partial {n_B}}} = \left( {\boldsymbol{x} - {\boldsymbol{x}_{\boldsymbol{p}}}} \right) \times \boldsymbol{p} \cdot {\boldsymbol{n}_{\boldsymbol{B}}}\;\;{\rm{\;\;\;\;}}({x,y}) \to {l_B}}\\[12pt]{\nabla {\chi _p} = 0{\rm{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;}}({x,y}) \to \infty }\end{array}} \right.\!,\end{align}

where $\boldsymbol{p}$ is the unit vector along the moment axis. The hypothetical potential ${\chi _p}$ vanishes at infinity and is made unique by requiring no circulation about any irreducible path. The flow-independent VMMs, shown in Fig. 23 for a NACA0012 aerofoil with different pitching axis ${\boldsymbol{x}_{\boldsymbol{p}}}$ , were constructed to identify the moment contribution of each given vortex in the flow field [Reference Li, Wang, Graham and Zhao121]. It is shown in the maps that an LEV moving away from the reference point or a TEV moving towards the reference point contributes to a nose-down pitching moment, and vice versa. Moreover, for any reference points located on the aerofoil except for the LE and TE, an LEV moving away from the edges or a TEV moving towards the edges contributes to a nose-up pitching moment, and vice versa. An impulsively started NACA0012 aerofoil has been used as an example. Good comparisons have been found between the moment given by CFD itself and by the precomputed VMM, together with the vortices obtained by CFD. As expected, the contribution from viscous forces on the pitching moment is negligible at the Reynolds numbers considered ( $50 \le Re \le 1{e^6}$ ), and the corresponding VMM can accurately reflect the total force contribution of the vortices in the viscous flow field. Similar to the VFM, the VMM method is insensitive to vortices far away from the body, which reflects the fact that pressure loads on the aerofoil are mainly due to near-body vortices in accordance with the Biot-Savart law. It is found that, due to the rolling up of LEVs and TEVs, the unsteady nose-down moment about the quarter chord is higher than the steady-state value. Both LEVs and TEVs can generate both positive and negative moments depending on their location. As an LEV grows and moves away from the body, its net contribution of the moment changes from positive to negative, while a TEV always contributes a net positive moment. The time-variation of the total moment is the overall effect of both LEVs and TEVs.

5.3.4 Multi-body flow vortex moment map

In a similar way the above VMM can be extended to the multi-body assembly shown in Fig. 22(a). The moment ${M_{ip}}$ acting on the $ith$ -body is, without completely rigorous derivation:

(64) \begin{align}{M_{ip}} = \rho {\zeta _i}\mathop \int\!\!\!\int \nolimits_{\rm{\Omega }} {\textsf{F}_{\boldsymbol{ip}}} \cdot \boldsymbol{U}{\omega _z}d{\rm{\Omega }}.\end{align}

where the VMM vector for the $i{\rm{th}}$ -body

(65) \begin{align}{\textsf{F}_{\boldsymbol{ip}}} = \frac{1}{{{\zeta _i}}}\left( {\frac{{\partial {\chi _{ip}}}}{{\partial y}}, - \frac{{\partial {\chi _{ip}}}}{{\partial x}}} \right)\end{align}

is normalised by a unit angular velocity multiplied by the characteristic length of the $ith$ -body. The flow-independent hypothetical potential ${\chi _{ip}}$ is defined as the velocity potential for hypothetical fluid motion induced by the rotation of the body ${{\rm{\Omega }}_{iB}}$ at unit angular velocity about an axis that passes through the reference point ${\boldsymbol{x}_{\boldsymbol{p}}}$ and perpendicular to the coordinate plane, which satisfies the following Laplace equation and boundary conditions

(66) \begin{align}\left\{ {\begin{array}{*{20}{l}}{{\nabla ^2}{\chi _{ip}} = 0}\\[12pt]{\dfrac{{\partial {\chi _{ip}}}}{{\partial {n_{iB}}}} = \left( {\boldsymbol{x} - {\boldsymbol{x}_{\boldsymbol{p}}}} \right) \times \boldsymbol{p} \cdot {\boldsymbol{n}_{\boldsymbol{iB}}}\;\;\;{\rm{\;\;\;}}({x,y}) \to {l_B}}\\[12pt]{\dfrac{{\partial {\chi _{ip}}}}{{\partial {n_{mB}}}} = 0{\rm{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;}}({x,y}) \to {l_{mB}} \cup {l_\infty }\left( {m \ne i} \right)}\end{array}} \right..\end{align}

A more rigorous derivation of the moment formula for multi-body assembly is still required. However, this method can separate the forces for individual bodies although high resolution is needed in narrow gaps.

5.4 Vortex force map method for 3D viscous flows

The 3D VFM method [Reference Li, Zhao and Graham122] described in this section has been derived and used to obtain body forces from the existing velocity/vorticity field near the body. The starting flow of a delta wing at high AoAs was used as an example to demonstrate the force method.

In the 3D body-fixed frame where the velocity $\boldsymbol{U} = \left( {u,v,w} \right)$ , and the vorticity $\boldsymbol\omega = \left( {{\omega _x},{\omega _y},{\omega _z}} \right) = \left( {\frac{{\partial w}}{{\partial y}} - \frac{{\partial v}}{{\partial z}},\frac{{\partial u}}{{\partial z}} - \frac{{\partial w}}{{\partial x}},\frac{{\partial v}}{{\partial x}} - \frac{{\partial u}}{{\partial y}}} \right)$ , the vortex force on the body can be written as

(67) \begin{align}{F_k} & = \iiint_{{\Omega }}\delta {F_k}d{\rm{\Omega }}\nonumber\\[4pt]\delta {F_k} & = \rho \!\left[ { - {\boldsymbol{\Lambda }}_{\boldsymbol{k},\textbf{3}}^{\boldsymbol{S}} \cdot \left( {u,v,0} \right){\omega _z} - {\boldsymbol{\Lambda }}_{\boldsymbol{k},\textbf{1}}^{\boldsymbol{S}} \cdot \left( {0,v,w} \right){\omega _x} - {\boldsymbol{\Lambda }}_{\boldsymbol{k},\textbf{2}}^{\boldsymbol{S}} \cdot \left( {u,0,w} \right){\omega _y}} \right] - \rho \boldsymbol{k} \cdot \left( {\boldsymbol\omega \times {\boldsymbol{V}_\infty }} \right),\end{align}

where the standard 3D VFM vectors are

(68) \begin{align}{\boldsymbol{\Lambda }}_{\boldsymbol{F},\textbf{1}}^{\boldsymbol{S}} = \left( {0, - \frac{{\partial \phi _k^S}}{{\partial z}},\frac{{\partial \phi _k^S}}{{\partial y}}} \right),{\boldsymbol{\Lambda }}_{\boldsymbol{k},\textbf{2}}^{\boldsymbol{S}} = \left( {\frac{{\partial \phi _k^S}}{{\partial z}},0, - \frac{{\partial \phi _k^S}}{{\partial x}}} \right),{\boldsymbol{\Lambda }}_{\boldsymbol{k},\textbf{3}}^{\boldsymbol{S}} = \left( { - \frac{{\partial \phi _k^S}}{{\partial y}},\frac{{\partial \phi _x^S}}{{\partial x}},0} \right).\end{align}

Here the standard hypothetical potential $\phi _k^S$ is the potential of an ideal flow with a unit incident velocity in the kth-direction about the body and governed by a 3D Laplace equation

(69) \begin{align}\begin{array}{*{20}{l}}{\left\{ {\begin{array}{*{20}{l}}{{\nabla ^2}\phi _k^S = 0}\\[4pt]{ - \nabla \phi _k^S \cdot {\boldsymbol{n}_{\boldsymbol{B}}} = 0{\rm{\;\;\;\;\;\;\;\;}}({x,y,z}) \to {S_B}}\\[4pt]{\nabla \phi _k^S = \boldsymbol{k} = \left( {{k_1},{k_2},{k_3}} \right){\rm{\;\;\;\;\;\;\;\;}}({x,y,z}) \to \infty }\end{array}} \right.}\end{array}.\end{align}

Note that at infinity, substituting $\nabla \phi _k^S{|_\infty } = \left( {{k_1},{k_2},{k_3}} \right)$ and $\left( {{k_1},{k_2},{k_3}} \right){|_\infty } = {{\boldsymbol{V}}_\infty }$ into Equation (67) we have $\delta {F_k}{|_\infty } = 0$ . This is compatible with the physical understanding that only near-body vortices induce pressure load on the body and affect the lift and drag. The vorticity at infinity appears to give a zero force contribution.

Li et al. [Reference Li, Zhao and Graham122] have discussed the case of a slender delta wing with ${\boldsymbol{V}_\infty } = \left( {{V_\infty }{\rm{cos}}\alpha, 0,{V_\infty }{\rm{sin}}\alpha } \right)$ (Fig. 24(a)). The force formula (67) together with (68) and (69) have been used to study the components of the force coefficient contributed by streamwise, spanwise, vertical vorticity ( ${\omega _x}$ , ${\omega _y}$ , ${\omega _z}$ ), and the correction term $ - \rho {\iiint _{\rm{\Omega }}}\boldsymbol{k} \cdot \left( {\boldsymbol\omega \times {\boldsymbol{V}_\infty }} \right)d{\rm{\Omega }}$ . The force from the LEV, mainly ${\omega _x}$ , is found to be dominant.

Figure 24. (a) A schematic diagram of vortex flow and force components for a delta wing at an incidence of $\alpha $ . (b) The analogy between the normal force of a slender delta wing and that of a 2D flat plate with steadily growing width and the analysis of the force contribution from a given vortex according to Wu’s [Reference Wu9] force formula. (c) Contours of vortex lift coefficient distribution (left column) and contours of streamwise vorticity (right column) at the mid-cross-section $x = 0$ of a delta wing at typical instants $\tau = 1$ (top line), $2.5$ (middle line), and $11.5$ (bottom line) after it impulsively started from rest at $\alpha = {60^o}$ and $Re = 5 \times {10^4}$ [Reference Gehlert and Babinsky123].

The VFM for the lift contribution of streamwise vorticity has been compared with that for ${F_N}cos\alpha $ on a 2D flat plate normal (at ${90^{\rm{o}}}$ incidence) to the flow to demonstrate the analogy between the normal force on a slender delta wing and that of a flat plate with a steadily growing span (see Fig. 24(b)). It is found that the analogy exists in the cross-flow sections close to the apex of the delta wing, even for relatively large AOAs, while the similarity doesn’t extend to the rear part, e.g., close to the TE. The precomputed VFMs, together with the vorticity field obtained by CFD shown in Fig. 24(c) (right column) can be used to calculate the force (here lift) contribution of the vorticity field and the results are shown in Fig. 24(c) (left column). Part of the LEV contributes a large positive lift and part of it contributes a small negative lift, leading to a positive total lift. The result is consistent with the conclusion from Wu’s theory that a pair of streamwise vortices contribute to positive lift when separating laterally and negative lift when moving towards each other (Fig. 24(b)). Moreover, the force contribution is mainly due to the conical vortex sheet (with a lateral flow speed of ${V_\infty }{\rm{cos}}\alpha {\rm{tan}}\beta $ ) rather than the central core. A quantitative understanding of the influence of vortical structures in four different evolution regimes on the body force has also been given in Ref. [Reference Li, Zhao and Graham122].

Note that Equation (67) can be rewritten in an alternative form

(70) \begin{align}{F_k} = \rho {\iiint _{\rm{\Omega }}}\left[ { - {{\boldsymbol{\Lambda }}_{\boldsymbol{k},\textbf{3}}} \cdot \left( {u,v,0} \right){\omega _z} - {{\boldsymbol{\Lambda }}_{\boldsymbol{k},\textbf{1}}} \cdot \left( {0,v,w} \right){\omega _x} - {{\boldsymbol{\Lambda }}_{\boldsymbol{k},\textbf{2}}} \cdot \left( {u,0,w} \right){\omega _y}} \right]d{\rm{\Omega }},\end{align}

with VFM vectors ${{\boldsymbol{\Lambda }}_{\boldsymbol{k},\textbf{1}}} = \left( {0, - \frac{{\partial {\phi _k}}}{{\partial z}},\frac{{\partial {\phi _k}}}{{\partial y}}} \right),{{\boldsymbol{\Lambda }}_{\boldsymbol{k},\textbf{2}}} = \left( {\frac{{\partial {\phi _k}}}{{\partial z}},0, - \frac{{\partial {\phi _k}}}{{\partial x}}} \right),{{\boldsymbol{\Lambda }}_{\boldsymbol{k},\textbf{3}}} = \left( { - \frac{{\partial {\phi _k}}}{{\partial y}},\frac{{\partial {\phi _x}}}{{\partial x}},0} \right)$ and hypothetical potential satisfying

(71) \begin{align}\begin{array}{*{20}{l}}{\left\{ {\begin{array}{*{20}{l}}{{\nabla ^2}{\phi _k} = 0}\\[4pt]{ - \nabla {\phi _k} \cdot {\boldsymbol{n}_{\boldsymbol{B}}} = \boldsymbol{k} \cdot {\boldsymbol{n}_{\boldsymbol{B}}}\;{\rm{\;\;\;\;\;\;\;}}({x,y,z}) \to {S_B}}\\[4pt]{\nabla {\phi _k} = 0{\rm{\;\;\;\;\;\;\;\;}}({x,y,z}) \to \infty }\end{array}} \right.}\end{array}\end{align}

Whether it is better to use Equation (67) or Equation (70) depends on the situation.

5.5 Application of VFM to practical force evaluation from PIV data

The following is an outline of how VFM may be used with PIV of the flow-field around a body to provide values of force on the body.

The typical approach to obtain fluid dynamic loads from the integration of surface pressures and skin friction in CFD has its difficulties in unsteady PIV data since resolving the entire boundary layer to an adequate resolution near the body surface is impractical. Additionally, though accurate experimental measurements on fluid dynamic loads in separated flows are possible with the development of high-precision experimental techniques, they are not easy to execute for low Reynolds number flow where the loads tend to be small and are subject to significant measurement errors and can be significantly contaminated by the resonance effect. Thus, the above-mentioned vortex methods can be very useful for extracting fluid dynamic forces in a non-intrusive way from the experimental measurements, such as PIV, of the flow fields.

The impulse force formulae have been extensively explored in extracting force from the unsteady PIV field. For instance, Noca et al. [Reference Noca, Shield and Jeon47] tested the expression for evaluating the instantaneous forces on a body in an incompressible cross-flow from the PIV measurements of the velocity and vorticity field in a finite and arbitrarily chosen region enclosing the circular cylinder. The results showed that the impulse force formula was accurate for fully resolved computational results, and captured the trend for the under-resolved experimental data. Later, Graham et al. [Reference Graham, Pitt Ford and Babinsky54] got good agreement between impulse formulae and force-balance measurements on a 2D translating wing flow by adding the influence of the missing vorticity via a vortex sheet, whose strength was obtained by the no-slip boundary condition. Following the same procedure by implicitly accounting for the incomplete PIV data close to the cylinder, applying the rotational and translational boundary condition of the cylinder kinematics to the cylinder surface, Gehlert and Babinsky [Reference Gehlert and Babinsky123, Reference Gehlert and Babinsky124] recovered the lift and drag from the PIV data on accelerating and rotating cylinder flows and validated it with force balance measurements. Ōtomo et al. [Reference Otomo, Henne, Mulleners, Ramesh and Viola125] successfully estimated unsteady lift due to a pair of LEV and TEV from PIV data.

The use of impulse methods in PIV data analysis requires a time sequence of data and a filtering process to align well with force balance measurements. These methods are unable to extract instantaneous forces from a single PIV data snapshot. Conversely, element-based techniques, such as auxiliary potential-based methods and the VFM method, are often underutilised. These methods hold significant potential for analysing unsteady flows without explicitly incorporating time derivatives. Li et al. [Reference Li, Wang, Graham and Zhao119] have facilitated the application of VFM method to computing body forces from PIV data by demonstrating a satisfactory level of force prediction accuracy using the VFM method on CFD data, even when relying on coarsely sampled data from a small, truncated domain. Building on this, Ōtomo et al. [Reference Otomo126] have further investigated the application of the VFM method to PIV data. This is the first application of the VFM method to PIV velocity data, where forces are extracted from snapshot PIV velocity and vorticity fields. A remarkable agreement between the direct force measurements and the VFM method has been observed for three different kinematics for surging flat plates and pitching NACA0018 aerofoils at Reynolds numbers of $O\!\left( {{{10}^4}} \right)$ . It has been found that the VFM method is highly robust to noise compared to the impulse method.