3.1. Vortex lift response model

The Anderson and DiPerna VRM will serve as the kernel function in the new lift spectrum model that is convolved with the vorticity wavenumber spectrum of the turbulent flow. In the Anderson and DiPerna VRM, a vortex with circulation, $\varGamma$ , drifts along a path, $Z_{s} (h,t)=X_{s} (h,t)+iY_{s} (h,t)$ , in the two-dimensional complex plane that follows a selected streamline with a vertical standoff position, $h$ , relative to the stagnation streamline well upstream of the aerofoil, as pictured in figure 7(a). The foil surface, $B$ , and the flow field are mapped into a domain, z, where the foil surface lies on a circle of radius, b, to efficiently achieve the impenetrable boundary condition by invoking the Milne-Thomson circle theorem (Milne-Thomsan Reference Milne-Thomsan1940), wherein an image vortex with circulation, $\varGamma _{I}=-\varGamma$ , is placed inside the foil at the image location, $z_{I} (t)=b^2/z_s^* (t)$ . To map any desired foil shape to the circle of radius, $b$ , as shown in figure 7(b), Anderson and DiPerna proposed the following generalised Joukowsky mapping function, according to

(3.1) \begin{equation} Z=z-s+\frac {(b \cos\beta -s_x)^2}{z-s}+\sum _{n=-1} ^{N}q_n^-(z-s)^{-n}, \end{equation}

where $s=s_x+is_y$ is the origin of the circle in the z-plane, $\beta =\arcsin {s_y/b}$ and $q_n^-$ are coefficients needed to match the shape of any desired foil shape. The real and imaginary components of $s$ provide thickness and camber to the foil, respectively, and the method to determine the coefficients, $q_n^-$ , is presented in the appendix of Anderson & DiPerna (Reference Anderson and DiPerna2022). If the coefficients, $q_n^-$ , are all set to zero, then (3.1) reduces to the traditional Joukowsky mapping function that is capable of producing three foil shapes: the teardrop foil with cusped trailing edge, the perfectly thin foil and the elliptically shaped foil. The last kinematic attribute of the model is the introduction of the foil wake vortex sheet, with circulation per unit length, $\gamma _w (Z_w (t))$ , and its image pairs, $\gamma _{wI} (Z_{wI} (t))$ , that allow for enforcement of the unsteady Kutta condition.

Figure 7.

Geometry of the vortex passing an aerofoil, the shed vortex sheet and the Blasius force equation contour of integration in (a) the physical $Z$ -plane and in (b) the cylinder $z$ -plane.

Through the conformal mapping approach, Anderson and DiPerna supplied expressions for the vortex kinematics to the Blasius unsteady force equation for the force per unit span, $\boldsymbol{F} _{\varGamma } (h,t)=d_{\varGamma } (h,t)+il_{\varGamma } (h,t)$ , to arrive at the following expression in the two-dimensional complex z-plane:

(3.2) \begin{align} \boldsymbol{F}_{\varGamma } (h,t)&=\left ( \frac {i\rho }{2}\int _{\mathcal{C}_z} \left ( \frac {\overline {W}^2(z)}{{\rm d}Z/{\rm d}z}\right ) {\rm d}z\right )^*+\left ( i\rho \int _{\mathcal{C}_z} \left ( W_v^{\prime }(z,t)+W_w^{\prime }(z,t)\right ) \left ( \frac {\overline {W}(z)}{{\rm d}Z/{\rm d}z}\right ) {\rm d}z\right )^* \nonumber\\[4pt] &\quad - i\rho \frac {\partial }{\partial t}\int _{\mathcal{C}_z}Re [ \tilde {F}(z,t)] \frac {{\rm d}Z}{{\rm d}z}{\rm d}z. \end{align}

The complex velocity potential is denoted as, $\tilde {F}(z,t)$ , and the complex velocity was decomposed according to

(3.3) \begin{equation} W(z,t)=\overline {W}(z)+W_v^{\prime }(z,t)+W_w^{\prime }(z,t), \end{equation}

where $\overline {W}(z)$ is the time-mean complex velocity, $W_v^{\prime }(z,t)$ is the unsteady complex velocity due to the passing vortex and $W_w^{\prime }(z,t)$ is the unsteady complex velocity induced by the vortices shed into the foil’s wake. Only the first-order velocity perturbation terms were retained in the $W^2 (z,t)$ term in producing equation (3.2), so that this expression could be used to superimpose the forces caused by a field of turbulent vortices. For turbulent flow applications, where the fluctuating-root-mean-square velocities are generally less than 10 % of the mean flow velocity, this first-order approximation of the unsteady force expression is appropriate.

Through extensive application of the residue theorem, Anderson and DiPerna solved the contour integral in (3.2), arriving at the following expression for the vortex-induced force acting on the aerofoil:

(3.4) \begin{align} \boldsymbol{F}_{\varGamma } (h,t)&=\overline {\boldsymbol{F}}_{\varGamma }+\boldsymbol{F}_{\!\varGamma ,qs}^{\prime}(h,t)+\boldsymbol{F}_{\!\varGamma ,mw}^{\prime }(h,t)+\boldsymbol{F}_{\varGamma ,mv}^{\prime }(h,t) \nonumber\\[3pt] & =-\,i\frac {\rho U\left (\varGamma _{m}+\varGamma \right )e^{i\alpha }}{\left (1+q_{-1}^{-}\right )^{*}}+i\rho \varGamma \frac {\partial Z_s}{\partial t} \nonumber\\ &\quad + i\rho \left (1+q_{-1}^{-}\right )\int _{-\infty }^{t}\dot {\gamma }_w\left (b,\tau \right )\left [1+\frac {b^2}{x_w^{2}} \right ]\frac {\partial x_w}{\partial t}{\rm d}\tau \nonumber\\[5pt] &\quad +i\rho \varGamma \left [\frac {b^2 \left (1+q_{-1}^{-}\right )}{\left (z_s^2\left (t\right )\right )^*}\frac {\partial z_s^*\left (t\right )}{\partial t}+\left [ \frac {\left (b-s\right )^2}{\left (z_s\left (t\right )-s\right )^2}\right . \right . \nonumber \\[5pt] &\quad +\left . \left . \sum _{n=1} ^{N}\frac {-n q_n^{-}\left (-1\right )^n}{\left (z_s\left (t\right )-z\right )^{\left (n+1\right )}}\right ]\frac {\partial z_s\left (t\right )}{\partial t}\right ]\!. \end{align}

This force expression is composed of four terms: the mean force, $\overline {\boldsymbol{F}}_{\varGamma }$ , the quasi-steady force due to the vortex, $\boldsymbol{F}_{\!\varGamma ,qs}^{\prime }(h,t,)$ , which is sometimes referred to as the impulsive force, the inertia force due to the shed vortex sheet, $\boldsymbol{F}_{\!\varGamma ,mw}^{\prime }(h,t)$ , and the force due to the inertia of the passing vortex, $\boldsymbol{F}_{\!\varGamma ,mv}^{\prime }(h,t)$ . Referring to (3.4), the mean flow approaches the foil with relative speed, $U$ , at an angle of attack, $\alpha$ . In the event that the foil is either cambered, or at a non-zero angle of attack, the mean force is related to the foil’s mean bound circulation, $\varGamma _m$ . For the quasi-steady force term, $\boldsymbol{F}_{\!\varGamma ,qs}^{\prime }(h,t,)$ , ${\partial Z_s}/{\partial t}$ is the velocity of the source vortex in the Z-plane. Lastly, in the shed wake inertia term, $\boldsymbol{F}_{\!\varGamma ,mw}^{\prime }(h,t)$ , $\dot {\gamma }_w$ is the time rate of change of vorticity shed into the cylinder wake due to imposition of the Kutta condition, and $x_W$ are the positions of the shed vorticities drifting in the cylinder wake.

The unsteady lift component of the force, $\boldsymbol{F}_{\varGamma } (h,t)$ , will be the focus of this study because it is considerably stronger than the unsteady drag component for aerofoils immersed in high-Reynolds-number turbulent flows. Based on (3.4), the unsteady coefficient of lift, $C_{l,\varGamma }^{\prime } (h,t)$ , defined by

(3.5) \begin{equation} C_{l,\varGamma }^{\prime } (h,t) = \frac {\mathrm{Im} \left [\boldsymbol{F}_{\!\varGamma }^{\prime } (h,t) \right ]}{\rho U \varGamma } , \end{equation}

is expressed as

(3.6) \begin{equation} \begin{aligned} C_{l,\varGamma }^{\prime } (h,t) &= \textit{Re} \left [ \frac {1}{U} \frac {\partial Z_s}{\partial t} \right ] + \frac {\left (1+q_{-1}^{-}\right )}{U \varGamma }\int _{-\infty }^{t}\dot {\gamma }_w\left (b,\tau \right )\left [1+\frac {b^2}{x_w^{2}} \right ]\frac {\partial x_w}{\partial t}{\rm d}\tau \\ &\quad + \frac {1}{U} \textit{Re} \left [\frac {b^2 \left (1+q_{-1}^{-}\right )}{\left (z_s^2\left (t\right )\right )^*}\frac {\partial z_s^*\left (t\right )}{\partial t}+\left [ \frac {\left (b-s\right )^2}{\left (z_s\left (t\right )-z\right )^2}\right . \right . \nonumber \\ &\quad +\left . \left .\sum _{n=1} ^{N}\frac {-n q_n^{-}\left (-1\right )^n}{\left (z_s\left (t\right )-z\right )^{\left (n+1\right )}}\right ]\frac {\partial z_s\left (t\right )}{\partial t} \right ]\! . \end{aligned} \end{equation}

The expression in (3.6) is the form of the VRM that will be convolved with the voticity wavenumber spectrum to compute unsteady lift due to turbulent flow, which will be described in the next section.

3.2. Lift spectrum due to turbulence

Toward deriving an expression for the lift on a foil immersed in a turbulent flow field, which utilises the VRM from (3.6), we begin by imagining that the turbulence vorticity in the cross-stream direction, $\omega _3$ , is a continuum distribution of the point vortices, $\varGamma$ , that were considered for the VRM. The turbulence is assumed here to drift along the streamlines in a frozen manner, according to Taylor’s hypothesis (Taylor Reference Taylor1938), in the moving reference frame, $X_m^{\prime }-h$ , where the $X_m^{\prime }$ -coordinate moves along the streamlines well upstream of the foil at the free-stream velocity, $U$ . The stationary reference frame, $X^{\prime }-h$ , as shown in figure 8, is centred on the stagnation streamline at a selected position well upstream of the foil. With the upstream, flow-aligned reference frames being related by, $X^{\prime }=X_m^{\prime }+Ut$ , the vorticity is related to the point vortices of the VRM according to

(3.7) \begin{align} \varGamma \left(t+\frac {X_m^{\prime }}{U},h\right) = \omega _3\left(t+\frac {X_m^{\prime }}{U},h\right){\rm d}X_m^{\prime }{\rm d}h . \end{align}

Figure 8.

Placement of point vortices, $ \varGamma (h,t)$ , to compute the VRM lift responses, $C_{l,\varGamma }^{\prime } (h,t)$ , needed for computing the foil lift response to an incident turbulent flow field in the physical $Z$ -plane.

Thus, to compute the lift from the entire field of turbulent vorticity, one can superimpose the VRM results for a field of point vortices, $\varGamma (t+({X_m^{\prime }}/{U}),h)$ , that have been placed on each streamline, defined by $h$ , with varying circulation strengths as a function of time delay, $ {X_m^{\prime }}/{U}$ . Invoking this superposition of the VRM from (3.6) with the field of vorticity from (3.7), gives the following relationship for the unsteady lift per unit span:

(3.8) \begin{equation} l^{\prime }(t) = \rho U \int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\omega _3\left(X_m^{\prime },h\right)C_{l,\varGamma }^{\prime } \left(h,t+\frac {X_m^{\prime }}{U}\right){\rm d}X_m^{\prime }{\rm d}h . \end{equation}

When computing the lift coefficients, $C_{l,\varGamma }^{\prime } (h,t+({X_m^{\prime }}/{U}))$ , for each streamline, the quadrature for $h$ must be chosen to avoid the stagnation streamline, where $h=0$ , as it results in a force singularity at the stagnation point on the foil. Using this two-dimensional lift response function in a strip theory manner, the unsteady lift acting on a wing of span, $S$ , can then be estimated for a three-dimensional cross-stream vorticity field, $\omega _3 (X_m^{\prime },h,s)$ , according to

(3.9) \begin{equation} L^{\prime }(t) = \int _{\frac {-S}{2}}^{\frac {S}{2}}l^{\prime }(s,t){\rm d}s=\rho U\int _{\frac {-S}{2}}^{\frac {S}{2}}\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\omega _3(X_m^{\prime },h,s)C_{l,\varGamma }^{\prime } \left(h,t+\frac {X_m^{\prime }}{U}\right){\rm d}X_m^{\prime }{\rm d}h {\rm d}s , \end{equation}

and in the temporal frequency domain according to

(3.10) \begin{equation} \tilde {L}^{\prime }(\omega ) = 2\pi \rho U\int _{\frac {-S}{2}}^{\frac {S}{2}}\int _{-\infty }^{\infty }\tilde {\omega }_3\left(\frac {\omega }{U},h,s\right)\tilde {C}_{l,\varGamma }^{\prime } (h,\omega ){\rm d}X_m^{\prime }{\rm d}h {\rm d}s . \end{equation}

The vorticity field, $\omega _3 (x,h,s)$ , is not generally known for a turbulent flow field because space–time-accurate measurements or numerical simulations over large three-dimensional domains for high-Reynolds-number flows are costly in terms of time and computational/measurement capabilities. For these ubiquitous situations, a stochastic modelling approach can adequately suffice. For the stochastic approach, the turbulent vorticity field, assumed spatially homogeneous in this case, is described in a reference frame, $X_m^{\prime }-h$ , which moves along with the flow in the $X^{\prime }$ -direction, in terms of the wavenumber cross-stream vorticity spectrum, $\varOmega _{33}(\boldsymbol{k})$ , where $\boldsymbol{k}$ is the wavenumber vector. Characterising the turbulent flow in terms of the vorticity wavenumber spectrum necessitates computation of the foil lift spectrum

(3.11) \begin{align} \varPhi _L(\omega ) &= \left ( 2\pi \rho U \right )^2 \int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\int _{-\infty }^{\infty } \varOmega _{33} \left (\frac {\omega }{U},k_2,k_3 \right ) \tilde {C}_{l,\varGamma }^{\prime }(h,\omega ) \tilde {C}_{l,\varGamma }^{\prime *}(h^{\prime },\omega ) \notag\\ &\quad \times \frac {\sin^2\left (k_3 \dfrac {S}{2}\right )}{\left(\dfrac {k_3 }{2}\right)^2} e^{i k_2 \left (h^{\prime }-h \right )} {\rm d}k_2 {\rm d}k_3 {\rm d}h {\rm d}h^{\prime } . \end{align}

If desired, the integrals over $h$ and $h^{\prime }$ can be evaluated to yield, $\mathcal{C}_{l,\varGamma }^{\prime }(k_2,\omega )$ , which simplifies (3.11) to

(3.12) \begin{equation} \varPhi _L(\omega ) = \left (2\pi \right )^4 \left (\rho U\right )^2 \int _{-\infty }^{\infty } \int _{-\infty }^{\infty } \varOmega _{33} \left (\frac {\omega }{U},k_2,k_3 \right ) \mathcal{C}_{l,\varGamma }^{\prime }(k_2,\omega ) \frac {\sin^2\left (k_3 \dfrac {S}{2}\right )}{\left( \dfrac {k_3}{2}\right)^2} {\rm d}k_2 {\rm d}k_3 . \end{equation}

At this point in the analysis, one is left to supply the model of the cross-stream vorticity wavenumber spectrum, and to implement the VRM to solve for the lift responses, $C_{l,\varGamma }^{\prime } (h,t)$ , to a series of point vortices set adrift on selected streamlines defined by $h$ , as shown in figure 8.

3.3. Lift spectrum for turbulence modelled using the von Kármán energy spectrum

A spatially homogenous turbulent flow has been assumed for the model presented in this paper, which only truly exists for flows in an infinite domain, but is a practical approximation for many flows. The cross-stream vorticity wavenumber transform for a flow over an infinite domain can be written in terms of the velocities according to

(3.13) \begin{equation} \tilde {\omega }_3(\boldsymbol{k}) = ik_2\tilde {u}_1(\boldsymbol{k}) - ik_1\tilde {u}_2(\boldsymbol{k}) , \end{equation}

which leads to the following relationship between the cross-stream vorticity spectrum, $\varOmega _{33} (\boldsymbol{k})$ , and the velocity wavenumber spectrum tensor, $\varPhi _{ij} (\boldsymbol{k})$ :

(3.14) \begin{equation} \varOmega _{33} (\boldsymbol{k}) = k_2^2\varPhi _{11}(\boldsymbol{k}) + k_1^2\varPhi _{22}(\boldsymbol{k}) - k_2k_1\left (\varPhi _{12}(\boldsymbol{k}) + \varPhi _{21}(\boldsymbol{k}) \right ). \end{equation}

Assuming isotropic turbulence, the velocity tensor spectrum is described by

(3.15) \begin{equation} \varPhi _{ij} (\boldsymbol{k}) = \frac {E(k)}{4\pi k^4}\left (k^2\delta _{ij} - k_i k_j \right ) , \end{equation}

where $E(k)$ is the energy spectrum. The cross-stream vorticity spectrum for isotropic turbulence can then be expressed in terms of the turbulence energy spectrum as

(3.16) \begin{equation} \varOmega _{33}(\boldsymbol{k}) = E(k) \frac {k_1^2 + k_2^2}{4\pi k^2} . \end{equation}

For the example cases to be considered later in this paper, we use the von Kármán energy spectrum model

(3.17) \begin{equation} E(k) = \frac {55 \varGamma (5/6)}{9\sqrt {\pi }\varGamma (1/3)}u_{\textit{rms}}^2 L_e \frac {(kL_e)^4}{\left [ 1 +(kL_e)^2 \right ]^{17/6}} , \end{equation}

where the length scale of the energy-containing flow structures, $L_e$ , is related to the integral length scale, $\varLambda$ , according to

(3.18) \begin{equation} L_e = \frac {\varGamma (1/3)}{\sqrt {\pi }\varGamma (5/6)}\varLambda . \end{equation}

Using the von Kármán energy spectrum model, the unsteady lift spectrum from (3.12) becomes

(3.19) \begin{align} \varPhi _L(\omega ) &= \dfrac {55\pi ^{1/2} \varGamma (5/6) \left ( \rho U u_{\textit{rms}}\right )^2 L_e}{9 \varGamma (1/3)} \int _{-\infty }^{\infty } \int _{-\infty }^{\infty } \dfrac {\left (\left (\dfrac {\omega }{u}\right )^2 + k_2^2\right )\left (\left (\dfrac {\omega }{u}\right )^2 + k_2^2 + k_3^2\right ) L_e^4 }{\left (1+ \left (\left (\dfrac {\omega }{u}\right )^2 + k_2^2 + k_3^2\right )L_e^2\right )^{17/6}} \notag\\ &\quad \times \mathcal{C}_{l,\varGamma }^{\prime }(k_2,\omega ) \dfrac {\sin^2\left (k_3 \dfrac {S}{2}\right )}{\left(k_3 \dfrac {S}{2}\right)^2} {\rm d}k_2 {\rm d}k_3 . \end{align}

To utilise this VRM, (3.19), to compute unsteady lift on a foil buffeted by a homogeneous turbulent flow, one must first compute the unsteady lift coefficients, $C_{l,\varGamma }^{\prime } (h,t)$ , according to (3.6) for a set of vortex standoff distances, $h$ . It should be noted that the choices of spatial resolution, $\Delta h$ , for the vortex standoff distances and the spatial extent over which to place the vortices when calculating the lift coefficients, $C_{l,\varGamma }^{\prime } (h,t)$ , will determine the frequency range over which the unsteady lift calculation will be accurate, and this will be discussed further when presenting example calculations in the next section. Then, a wavenumber transform of these lift coefficients will be performed along the $h$ -dimension to obtain the lift vortex lift response coefficient wavenumber–frequency spectrum, $\mathcal{C}_{l\varGamma } (k_2,\omega )$ . Then, to carry out the calculation of (3.19), one only needs knowledge of the foil span, the mean flow speed and an estimate of the turbulence fluctuating root-mean-square velocity and the length scale of the energy-containing structures, which is related to the integral length scale according to (3.18).

3.4. Verification of vortex response model for a thin foil

To verify the VRM for turbulence-induced unsteady foil lift, (3.19), we begin by benchmarking this model for the case of a thin foil immersed in turbulent flow, a case for which analytical models exist. The VRM unsteady lift spectrum due to turbulence for a thin foil will utilise the thin-foil VRM lift coefficient, which was shown by Anderson and DiPerna to be

(3.20) \begin{equation} {C}_{l,\varGamma ,\textit{thin}}^{\prime } (h,t) = -b \textit{Re}\left [\frac {1}{z_s+b}\right ] + \frac {1}{2\varGamma }\int _{-\infty }^{t}\dot {\gamma }_w\left (b,h,\tau \right )\left [\frac {x_W(h,t-\tau )-b}{x_W(h,t-\tau )-b} \right ] {\rm d}\tau . \end{equation}

This lift coefficient was computed over a selected range of $h$ -values, and then the temporal and spatial Fourier transforms were computed to arrive at $\mathcal{C}_{l,\varGamma ,\textit{thin}} (k_2,\omega )$ .

Here, the VRM is benchmarked against the analytical modelling approach described by Liepmann (Reference Liepmann1955), which used an upwash velocity lift response function convolved with a model of the statistically expected upwash velocity wavenumber–frequency spectrum, $\varPhi _{v_g} (k_3,\omega )$ , where $k_3$ is the foil spanwise wavenumber, which is modelled here using the von Kármán isotropic turbulence energy spectrum model. The Liepmann lift spectrum is given by

(3.21) \begin{align} \varPhi _{L,Thin,nD}(\omega ) &= \left (\dfrac {55\varGamma (7/3)\pi ^{1/2}}{18\varGamma (17/6)} \right )\left ( \rho ^2C^2u_{\textit{rms}}^2U\varLambda L_e \right ) \notag\\ &\quad\times \int _{-\infty }^{\infty }\dfrac {L_e^2\left (\left (\dfrac {\omega }{u}\right )^2+k_3^2 \right ) G_{nD}\left (\dfrac {\omega C}{2U},k_3\right )}{\left [1 + L_e^2\left (\left (\dfrac {\omega }{u}\right )^2+k_3^2 \right ) \right ]^{14/6}} \dfrac {\sin^2\left (k_3 \dfrac {S}{2}\right )}{\left(\dfrac {k_3}{2}\right)^2} {\rm d}k_3 , \end{align}

where $G_{nD} (({\omega C}/{2U}),k_3 )$ is the upwash velocity lift response function. Sears (Reference Sears1941) provided the most commonly used two-dimensional upwash velocity response function, which is

(3.22) \begin{equation} G_{2D}\left (\frac {\omega C}{2U},k_3=0\right )=S_e\left (\dfrac {\omega C}{2U}\right ) = \dfrac {2}{\pi \dfrac {\omega C}{2U}\left [ H_0^{(2)}\left (\dfrac {\omega C}{2U}\right ) - iH_1^{(2)}\left (\dfrac {\omega C}{2U}\right )\right ]} , \end{equation}

where $H_0^{(2)}$ and $H_1^{(2)}$ are Hankel functions. To assess the effects of using a two-dimensional lift response function, we will compare these results with those computed with the three-dimensional upwash lift response model by Filotas (Reference Filotas1969)

(3.23) \begin{equation} G_{3D}\left (\dfrac {\omega C}{2U},k_3\right ) = \dfrac {e^{-i(2k/C)}\left [\sin\beta - \dfrac {\pi \beta (1+1/2 \cos\beta )}{1+\pi (2k/C)(1+1/2 \cos\beta )} \right ]}{\sqrt {1+\pi (2k/C)(1+\sin^2\beta +\pi (2k/C)\cos\beta )}} , \end{equation}

where $\beta = a\tan( ({\omega C/2U})/{k_3})$ .

The goal for this model verification exercise was to achieve agreement between the VRM and the Sears model for frequencies, $1\lt \omega C/2U\lt 50$ . To help ensure achieving the low-frequency goal for the VRM approach, a criteria was set to include all streamlines from where a vortex contributes to the foil upwash velocity to within 30 dB of the vortex with infinitesimal offset distance from the foil. Since the contribution to the upwash from a vortex drifting by with standoff distance, $h$ , is proportional to $\exp { [-\omega h/U ]}$ , the vertical standoff distance where this value is −30 dB at a reduced frequency, $\omega C/2U = 1$ , occurs at $h = 1.72 C$ . The maximum standoff distance of $h = 2C$ was chosen to meet this criteria. To confirm if this criteria was actually met, we examine the vortex lift response coefficient plotted versus vortex offset distance, $h/C$ , and reduced frequency, $\omega C/2U$ , in figure 9. By examining figure 9 at $\omega C/2U = 1$ , it can be seen that $\tilde {C}_{l,\varGamma } (h = 0,\pi fC/U = 1)\approx$ −17 dB and that $\tilde {C}_{l,\varGamma } (h = 1.5C,\pi fC/U = 1)\approx -50 \,\rm dB$ , thus confirming that the goal of including all vortices that will contribute to the lift at $\pi fC/U =1$ to within 30 dB was achieved using a vertical offset range of $\pm 2C$ . To achieve accuracy in estimating the lift spectrum at the upper frequency limit, $\omega C/2U = 50$ , it was determined that the vortex lift response function wavenumber–frequency spectrum, $\mathcal{C}_{l,\varGamma }(k_2,\omega )$ , should cover a wavenumber span that captures amplitudes within 20 decibels of the maximum value at a given frequency. From the plot of $\mathcal{C}_{l,\varGamma }(k_2,\omega )$ for the thin foil, shown in figure 10, it can be seen that this criteria is satisfied at all frequencies, up to $\omega C/2U = 50$ , at a wavenumber, $k_2 C = 250$ radians. To adequately resolve wavenumbers, $k_2 C \leqslant 250$ radians, it was decided to sample at a wavenumber, $k_{2,sample} C=1250$ , which equated to solving the vortex lift response function at a streamline spacing of $\Delta h/C = 0.005$ for $-2\lt \Delta h/C \lt 2$ .

Figure 9.

Vortex lift response function spectrum, $\tilde {C}_{l,\varGamma } (h,\omega )$ , for a thin foil and turbulent flow parameters listed in table 2.

Figure 10.

Vortex lift response wavenumber–frequency spectrum, $\mathcal{C}_{l,\varGamma } (k_2,\omega )$ , for a thin foil and turbulent flow parameters listed in table 2.

The lift spectra were computed using the VRM approach, by evaluating equation (3.19) at the streamline spacing and range identified to achieve accuracy for frequencies, $1\lt \omega C/2U \lt 50$ , and using the Sears model approach, (3.21), for a thin foil in a spatially homogenous turbulent flow described by the parameters listed in table 2. The computed lift spectra using the two approaches agree almost identically for the targeted frequency range, $1\lt \omega C/2U\lt 50$ , as can be seen in figure 11, thus confirming the accuracy of the VRM approach for computing the turbulence-induced lift spectrum on a foil. As an additional step, the error of using a two-dimensional lift response model for this case was studied by computing the lift spectrum with the Filotas three-dimensional gust response model, (3.23), combined with (3.21). As seen in figure 11, the Filotas three-dimensional model lift spectrum nearly identically matches the two-dimensional Sears spectrum in this case for $\omega C/2U\gt 2$ , and the difference is less than 0.5 dB for $\omega C/2U\lt 2$ , which shows that two-dimensional lift response models can be quite accurate for cases when $S\gt \varLambda$ .

Figure 11.

Comparison of VRM and Sears model for lift spectrum on a thin foil due to turbulence for thin foil and turbulent flow parameters listed in table 2.

5.1. Numerical study of varied leading-edge thicknesses

We utilise the VRM to analyse three foil shapes with 8 % thickness-to-chord ratios and the same foil chord, span and turbulence parameters as the Lysak experiment. The foil shapes to be studied here, the NACA $65_{1A}$ -0008, the NACA 0008-65 and the NACA 0008-61 aerofoils, all previously shown in figure 2, have progressively thicker, more rounded leading edges while maintaining the same maximum thickness. The mapping coefficients used in the mapping function, (3.1), for these three foil shapes are listed in table 4. The computed lift spectra for all three aforementioned foil shapes are plotted in figure 18, along with the result for a thin foil as a reference. As seen in figure 18, the NACA 0008-61 foil produces greatly attenuated high-frequency lift compared with the other foil shapes, particularly the NACA $65_{1A}$ -0008 with its sharp leading edge. For example, the lift spectral density level for the NACA 0008-61 foil at $\omega C/2U=30$ is approximately 8 dB lower than the NACA $65_{1A}$ -0008 foil, and the difference is nearly 15 dB at $\omega C/2U=50$ . Relative to the NACA 0008-65 foil, the NACA 0008-61 produced lift levels that were 5 and 9 dB lower at frequencies, $\omega C/2U$ of 30 and 40, respectively. This large attenuation for merely moving the maximum thickness location closer to the leading edge motivated a closer look at the physics responsible for this effect, described in § 5.2, and the conduct of a wind tunnel experiment, described in § 2.3, to verify this impressive effect of thickening the leading edge.

Figure 18.

Comparison of lift spectra computed using the VRM for different foils, all with 8 % thickness-to-chord ratio and chord, span and turbulence parameters the same as the Lysak experiment, but each having different section shapes.

5.2. Explanation for unsteady lift attenuation due to thick leading edges

Toward developing a physics-based understanding of how thick leading edges cause attenuation of unsteady lift response to turbulent flow, we begin by revisiting the vortex lift response coefficient model in (3.6), where the unsteady lift is ascribed to three primary mechanisms: the inertia of the passing vortex, the inertia of the shed wake vorticity and the quasi-steady term. The three contributions to the unsteady lift are decomposed in figure 19 for the case when a unit strength vortex passes by a NACA 0007-65 aerofoil at an undistorted offset distance of $h/C$ = 0.0025. The total lift, shown as the black curve in figure 19, reaches peak levels just before the vortex passes the leading edge at $Ut/C\approx$ −1/2, and during this time it is almost entirely attributed to the vortex inertia lift, making this term the focus of our quest for a physical explanation for how the leading-edge shape affects the leading-edge pulse in lift, which will control the high-frequency levels.

Figure 19.

Mechanisms comprising the lift on a NACA 0007-65 aerofoil due to a vortex passing at $h/C$ = 0.0025.

If the baseline Joukowsky mapping coefficient, $s$ , is selected to closely approximate the actual foil thickness, then the mapping perturbation coefficients, $q_{n}^-$ , become small relative to the circle radius, $b$ , which allows the vortex inertia lift coefficient to be approximated by

(5.1) \begin{equation} \begin{aligned} C_{l,\varGamma ,mv}^{\prime } (h,t) \approx & \frac {1}{U} Re \left [ \frac {b^2 }{\left (z_s^2\left (t\right )\right )^*}\frac {\partial z_s^*\left (t\right )}{\partial t}+\frac {\left (b-s\right )^2}{\left (z_s\left (t\right )-s\right )^2}\frac {\partial z_s\left (t\right )}{\partial t} \right ] . \end{aligned} \end{equation}

To attenuate the peak lift, based on (5.1), the foil must be shaped in a way that either reduces the vortex’s drift velocity, ${\partial z_s (t )}/{\partial t}$ , or increases the position of the vortex from the origin, as it drifts near the foil leading edge. For the scenario of a vortex at an offset distance of $ {h}/{C}=0.0025$ passing by a thin foil compared with a NACA 0007-65 foil, the vortex drift velocity in the circle plane is significantly attenuated for the thick foil as the vortex approaches the leading edge at time, $ {Ut}/{C}\approx -0.5$ , as shown in figure 20(a), while the vortex distance from the origin in the circle plane is essentially the same for the thin foil and the thick foil, as shown in figure 20(b). In the physical plane of the foil, the Z-plane, the foils with a thicker leading edge reduce the real component of the vortex drift velocity as the vortex approaches the leading edge at times just before ${Ut}/{C}\approx -0.5$ , as shown in figure 21 for the case of the NACA 0007-61 compared with the NACA 0007-65. Since lift is proportional to $\rho \varGamma Re [{\partial Z_s}/{\partial t} ]$ at the moments when the vortex is closely approaching the foil leading edge, the ability of a foil with a thick leading edge to reduce the real component of the vortex drift velocity as it approaches the foil leading edge has a profound effect on reducing the peak lift and, consequently, the high-frequency portion of the lift spectrum. This observation that the leading-edge bluntness serves to reduce high-frequency unsteady lift because it slows down the convection speed of the vorticity approaching the foil nose is the same mechanism cited by Ayton & Chaitanya (Reference Ayton and Chaitanya2017) that serves to also attenuate leading-edge sound due to turbulence–foil interactions. Also, the foil with a thicker leading edge will quickly accelerate the vortex just after it passes the leading edge, as seen in figure 21, but this does not affect the total lift because the vortex inertia is no longer the single dominant contribution to the lift once the vortex has drifted approximately 5 %–10 % of the chord downstream from the leading edge, as was shown in figure 19.

Figure 20.

(a) Normalised vortex drift velocity magnitude in the circle plane and (b) the source vortex inverse distance from the origin in the circle plane for a vortex passing by a thin foil and a NACA 0007-65 foil at a vortex offset distance, $ {h}/{C}=0.0025$ .

Figure 21.

Normalised real component of the vortex drift velocity in the foil plane for a vortex passing by a thin foil, a NACA 0007-65, and a NACA 0007-61 foil at a vortex offset distance, $h/C$ = 0.0025.

5.3. Simplified high-frequency lift model

As was described in § 5.2, the amplitude of the peak lift pulse, which translates to high-frequency lift, is controlled by the vortex inertia and, more specifically, the speed of the vortices as they approach the foil leading edge. This indicates that only the vortex inertia in the lift vortex response expression, (3.6), needs to be considered when estimating the high-frequency portion of the lift spectrum, which will expedite the computation time of the many vortex responses over the range of distances, $h$ , that are needed to arrive at the wavenumber–frequency lift response function, $\mathcal{C}_{l,\varGamma }(k_2,\omega )$ .

To simply focus on the vortex inertia term to estimate the high-frequency lift is not adequate, because this term also introduces a large pulse as the vortex passes the foil trailing edge, as shown in figure 19, which is nullified by the wake inertia. Therefore, the aim here is to mathematically remove the trailing-edge pulse from the vortex inertia term. We begin by approximating the vortex velocity, also referred to as the vortex drift velocity, according to

(5.2) \begin{equation} \begin{aligned} \frac {\partial z_s(t)}{\partial t} \approx & \frac {\partial z_s}{\partial Z_s}\left |\frac {\partial z_s}{\partial Z_s}\right |^2, \end{aligned} \end{equation}

where Routh’s correction (Clements Reference Clements1973) was ignored because it generally has a small effect. This allows the unsteady lift due to the inertia of the vortex to be approximated as

(5.3) \begin{equation} \begin{aligned} C_{l,\varGamma ,mv}^{\prime } (h,t) \approx & Re \left [ \frac {b^2 }{\left (z_s^2\right )^*}\left (\frac {z_s^2-b^2}{z_s^2}\right )+\frac {\left (b-s\right )^2}{\left (z_s\left (t\right )-s\right )^2}\left (\frac {z_s^2-b^2}{z_s^2}\right )^* \right ] \left |\frac {\partial z_s}{\partial Z_s}\right |^2. \end{aligned} \end{equation}

The expression in (5.3) is a product of two terms, the first consisting of $Re[b^2/(z_s^* )^2 (((z_s)^2-b^2)/(z_s)^2 ) +(b-s)^2/(s-z_s )^2 (((z_s)^2-b^2)/(z_s)^2 )^* ]$ and the second, $|(\partial z_s)/(\partial Z_s )|^2$ . The first term is related to the mean drift path, $z_s$ , in the circle plane, and the second, being a derivative of the mapping function, quantifies the distortion of the coordinates when mapped from the circle plane to the physical plane. This distortion is greatest at locations of smallest radius of curvature of the foil shape, typically the leading and trailing edges, which is confirmed by the large pulses in the curve for $|(\partial z_s)/(\partial Z_s )|^2$ , shown in figure 22, resulting from a vortex at $h/C=0.0025$ drifting by a NACA 0007-65 foil. The $z_s$ drift term typically exhibits a more relatively muted lift response near the leading and trailing edges, as shown in figure 22. Therefore, the mapping derivative term is the likely indicator of foil leading-edge shape effects on the peak lift amplitude as vortex passes the leading edge.

Figure 22.

Decomposition of vortex inertia force on the NACA 0007-65 aerofoil due to a vortex passing at $h/C=0.0025$ .

To focus on the vortex inertia leading-edge force pulse, and to remove the mathematical terms associated with the trailing-edge pulse, the mapping function derivative is recast in terms of a partial fraction expansion, and only the term associated with the pole, $p_{\textit{LE}}$ , nearest the leading edge is retained to provide the following approximation:

(5.4) \begin{equation} \begin{aligned} \left (\frac {\partial z_s}{\partial Z_s}\right )_{\textit{LE}} \approx & \frac {r_{\textit{LE}}}{z_s(t)-p_{\textit{LE}}}+1. \end{aligned} \end{equation}

When using the general mapping function, (3.1), the polynomials in the mapping function derivative make it exceedingly difficult to analytically determine the poles and residues of the partial fraction expansion, so a computational tool will generally be required; however, the partial fraction for the thin foil is readily determined to be

(5.5) \begin{equation} \begin{aligned} \left (\frac {\partial z_s}{\partial Z_s}\right )_{{\textit{LE}},\textit{thin}} \approx & \frac {\frac {-b}{2}}{z_s(t)+b}+1, \end{aligned} \end{equation}

where $b$ is the radius of the circle representing the foil surface in the z-plane. Equation (5.5) shows that the leading-edge pole for a thin foil exists at the foil leading edge, $z=-b$ . For foils with non-zero thickness, the leading-edge pole will move away from the leading edge along the real z-axis, when the foil is uncambered, and the residue will become a slightly smaller negative value. As an example, the leading-edge pole and residue for a traditional Joukowsky foil shape are, $p_{{\textit{LE}},Joukowsky}=2s-b$ and $r_{{\textit{LE}},Joukowsky}=(s-b)/2$ , respectively, which shows that foil thickness, accounted for by the real component of $s$ , moves the pole at a rate twice that of the residue. Therefore, increasing foil thickness will have the mathematical effect of moving the leading-edge pole increasingly further away from the leading edge, thereby increasing the minimum distance between the vortex path, $z_s (t)$ , and the pole, $p_{\textit{LE}}$ , thus decreasing the leading-edge force pulse. The locations of the mapping derivative poles for the NACA 0007-65 and 0007-61 foil shapes, plotted in figure 23(a), show migration of the poles, with the exception of the trailing-edge pole, toward the origin with increased foil leading-edge thickness. The relative distances between the leading-edge poles and the closest position of approach (CPA) for the passing vortex with offset distance, $h/C=0.0025$ , plotted in figure 23(b), show that the added leading-edge foil thickness for the 0007-61 foil shape increases the distance, $\left |z_{s,\textit{CPA}}-p_{\textit{LE}}\right |$ , relative to the 0007-65 foil.

Figure 23.

(a) Mapping function derivative poles for the NACA 0007-65 and 0007-61 foil shapes, and (b) a close-up view of the leading-edge poles.

The mapping derivative poles and residues are independent of the passing vortex paths, $z_s (h,t)$ , so $r_{\textit{LE}}$ and $p_{\textit{LE}}$ need only be computed once for a given foil shape, yielding a comparatively quick calculation relative to the full expression in (3.4). It is proposed here that the high-frequency, leading-edge vortex lift response function, $\tilde {C}_{l,\varGamma ,LE}^{\prime } (h,f)$ , can be approximated by relating it to the thin-foil high-frequency, leading-edge votex lift response function, $\tilde {C}_{l,\varGamma ,LE,\textit{thin}}^{\prime } (h,f)$ , according to

(5.6) \begin{equation} \begin{aligned} \tilde {C}_{l,\varGamma ,LE}^{\prime } (h,f) \approx & \tilde {C}_{l,\varGamma ,\textit{LE},\textit{thin}}^{\prime } (h,f)\frac {\left |\dfrac {\partial z_s}{\partial Z_s}(h,f)\right |^2_{\textit{LE}}}{\left |\dfrac {\partial z_s}{\partial Z_s}(h,f)\right |^2_{{\textit{LE}},\textit{thin}}}, \end{aligned} \end{equation}

where $\tilde {C}_{l,\varGamma ,LE,\textit{thin}}^{\prime } (h,f)$ can be computed by performing either the Fourier transform of (3.20), or by using the following expression in terms of the Sears function:

(5.7) \begin{equation} \begin{aligned} \tilde {C}_{l,\varGamma ,LE,\textit{thin}}^{\prime } (h,f) = & i\pi \frac {Ce^{-2\pi fh/U}}{2U}S_e\left (2\pi fC/2U\right ). \end{aligned} \end{equation}

The accuracy of the high-frequency leading-edge lift model, (5.6), will be demonstrated in the next section, where it is utilised to predict the unsteady lift measurements performed in the wind tunnel experiment described in § 2.3.

5.4. Comparison of model predictions and experimental measurements

To physically confirm the model-based predictions from § 5.1 that showed significant high-frequency lift attenuation for foils with thickened leading edges, a wind tunnel test was performed, as described in § 2.3, wherein the unsteady lift was measured for a NACA 0007-65 and a NACA 0007-61m aerofoil immersed in a homogeneous turbulent flow characterised by the parameters listed in table 1. The lift was measured using the sum of four force gauges, labelled as FG1 to FG4 in figure 5. Unavoidable structural vibrations prevented direct measurement of the turbulence-induced aerodynamic lift force, so a novel speed-scaling technique, described in Appendix A, was employed to extract the aerodynamic lift force from the measured force gauge signals.

The measured unsteady lift spectrum for the NACA 0007-65 aerofoil due to the turbulent inflow, inferred from the speed-scaling analysis of the force gauge measurement data, is compared with the VRM and high-frequency, leading-edge force model predictions in figure 24. The VRM shows excellent agreement with the measured lift spectrum throughout the entire frequency range, and the high-frequency model becomes very accurate for frequencies, $\omega C/2U\gt 10$ . The NACA 0007-65 aerofoil shape is thus confirmed to provide only modest attenuation, approximately 3–6 dB, relative to an infinitesimally thin foil at frequencies, $\omega C/2U\gt 10$ .

Figure 24.

Comparison of full VRM and high-frequency, leading-edge model predictions with measured lift spectra for the NACA 0007-65 aerofoil, with a calculated thin-foil lift spectrum for reference.

Similar excellent agreement between model-based predictions and lift measurements was also achieved for the NACA 0007-61m, as shown in figure 25. Both models accurately predicted the increased attenuation levels for high frequencies due to the increased leading-edge thickness of the NACA 0007-61m, which was approximately 15 dB lower than an infinitesimally thin foil at the frequency, $\omega C/2U = 30$ . It is also worth noting that the high-frequency, leading-edge lift model maintained excellent accuracy down to a rather low frequency of $\omega C/2U = 5$ . These experimental results confirmed the theoretical prediction that increasing the foil leading-edge thickness can significantly reduce high-frequency lift for aerofoils immersed in turbulent flow.

Figure 25.

Comparison of full VRM and high-frequency, leading-edge model predictions with measured lift spectra for the NACA 0007-61m aerofoil, with a calculated thin-foil lift spectrum for reference.

5.5. Discussion of model uncertainty

While this new model is shown to accurately match two measurement data sets, it is wise to understand where the accuracy of this model might be challenged. Some error may be introduced at lower frequencies because the model utilises a two-dimensional vorticity lift response function. As was shown by Manar & Anderson (Reference Manar and Anderson2022), this will only cause small errors of less than 2 dB at low frequencies, $\omega C/2U\lt 1$ , and minimal errors of approximately 1 dB in the mid-frequency range, $1\lt \omega C/2U\lt 10$ . Manar and Anderson also showed that these two-dimensional lift response effects are further minimised when the turbulence integral length scale is much smaller than the foil span. The lift spectrum calculated for the flow case considered in the wind tunnel experiment using the two- and three-dimensional upwash lift response functions, plotted in figure 11, showed close agreement over most of the frequency range, with the largest difference of 0.5 dB found only at the lowest frequencies, $\omega C/2U\lt 2$ .

The second model limitation arises from use of Taylor’s frozen turbulence assumption, which could introduce error in modelling the turbulence velocity coherence for large streamwise separations at high frequencies. Tong (Reference Tong1996) showed that invoking Taylor’s hypothesis most prominently affects the coherence of the streamwise velocity component for streamwise separation distances at high frequencies. Their results showed that the coherence is underestimated by approximately 0.36, 1 and 2 dB at frequencies, $\omega \eta _1/U$ , of 1, 2 and 3, respectively, where $\eta _1$ is the streamwise separation distance. As was seen in § 5.2, the peak lift in the response to a vortex passing near the leading edge is impulsive, occurring over the short time duration when the approaching vortex moves an approximate streamwise distance, $\eta _{1c} \approx C/20$ , and this impulsive portion of the lift response controls the lift spectrum levels at high frequencies. Thus, Taylor’s hypothesis needs to induce minimal error for streamwise velocity coherence over this critical streamwise distance, $\eta _{1c}$ . Applying this critical streamwise distance to the streamwise velocity coherence error levels observed by Tong, gives possible errors in the computed lift spectra of 0.36, 1 and 2 dB at frequencies, $\omega C/2U$ , of 10, 20 and 30, respectively. Therefore, Taylor’s frozen turbulence assumption could introduce a possible 2 dB error at the upper end of the frequency range studied in this work.