3.1. Induced velocity due to normal forcing

The normal induced velocity along the wake centreline $(x_*, y_*=0)^\top$ is obtained by solving the normal component of the Biot–Savart integral given in (3.1) for a generic normal forcing $C_y(t_*)$ . First, we focus here on the contribution from the second term of the unsteady vorticity solution $\omega ^{\textit{us}}_*$ (the time integral in (2.5)). By switching the order of temporal and spatial integration, we may obtain

(3.2) \begin{align} v^{\textit{us}}_*(x_*,0,t_*)=\int \limits _0^{t_*}\!\! \frac {1}{2\pi }\int \limits _{-\infty }^{+\infty } \int \limits _{-\infty }^{+\infty } \!\! \frac {(x_*-x_*^\prime ) (x_*^\prime +s-t_*)}{(x_*-x_*^\prime )^2+(-y_*^\prime )^2} \frac {C_y(s)}{\pi \epsilon _*^4} \text{e}^{-((x_*^\prime +s-t_*)^2+y_*^{\prime 2})/\epsilon _*^2} \: \mathrm{d}x_*^\prime \mathrm{d} y_*^\prime \mathrm{d} s, \end{align}

which is the normal perturbation velocity induced by the vorticity shed as a result of the unsteady forcing in the time interval $t_*\in [0,t_*]$ . Integration in the two spatial directions leads to

(3.3) \begin{align} v^{\textit{us}}_*(x_*,0,t_*)= \int \limits _0^{t_*} \frac {-C_y(s)}{2\pi } \bigg [ \frac {\text{e}^{-(x_*+s-t_*)^2/\epsilon _*^2}}{\epsilon ^2} + \frac {\text{e}^{-(x_*+s-t_*)^2/\epsilon _*^2} -1}{2(x_*+s-t_*)^2} \bigg ] \: \mathrm{d} s. \end{align}

We note that the integrand in the limit $(x_*+s-t_*)\rightarrow 0$ tends to $-C_y(s)/(4\pi \epsilon _*^2)$ , i.e. it is not singular. The derived solution shows the identical structure as the formula obtained by Martínez-Tossas et al. (Reference Martínez-Tossas, Sakievich, Churchfield and Meneveau2024) for their generalised steady three-dimensional filtered lifting line theory. This equivalence stems from the fact that the two problems are similar upon applying a rotation of $90^\circ$ to their coordinate system and replacing the spatial convolution integral in their solution with a spatio-temporal convolution as in (3.3). So far, only the contribution from the second term of the vorticity solution $\omega ^{\textit{us}}_*$ has been considered. Thus, the induced velocity solution in (3.3) has to be complemented with a term $v^{\textit{IC}}_*(x_*,0,t_*)$ , accounting for the velocity induced by the initial vorticity field in case $\omega _*^{\textit{IC}} \neq 0$ in (2.5). If we assume that the initial condition is given by the steady-state solution of the flow over a Gaussian body force, one can take the initial condition derived by Martínez-Tossas et al. (Reference Martínez-Tossas, Churchfield and Meneveau2017), i.e. the velocity field induced by a Lamb–Oseen vortex with core size $\epsilon _*$ advected downstream by the background flow:

(3.4) \begin{align} v^{\textit{IC}}_*(x_*-t_*,0) = \frac {-C_y(t_*=0)}{4\pi } \frac {1 - \text{e}^{-((x_*-t_*)/\epsilon _*)^2}}{x_*-t_*}. \end{align}

The complete solution for the normal induced velocity along the wake centreline follows by combining (3.3) and (3.4) as

(3.5) \begin{align} v_*(x_*,0,t_*)=v^{\textit{IC}}_*(x_*-t_*,0)+v^{\textit{us}}_*(x_*,0,t_*). \end{align}

While for numerical implementation purposes the structure of the solution in (3.3) is beneficial since it does not require the calculation of the temporal gradient of $C_y(t_*)$ , it is insightful to rewrite the time integral using integration by parts, i.e.

(3.6) \begin{align} v^{\textit{us}}_*(x_*,0,t_*) & = -\bigg [ \frac {C_y(s)}{4\pi } \frac {1 - \text{e}^{-((x_*+s-t_*)/\epsilon _*)^2}}{x_*+s-t_*} \bigg ]_{s=0}^{s=t_*} \nonumber \\ & \quad + \int \limits _0^{t_*} \frac {1}{4\pi } \frac {\mathrm{d} C_y(s)}{\mathrm{d}s} \frac {1 - \text{e}^{-((x_*+s-t_*)/\epsilon _*)^2}}{x_*+s-t_*} \: \mathrm{d} s, \end{align}

where we define the indicial response function as

(3.7) \begin{align} \varphi (x_*+s-t_*) = \frac {1 - \text{e}^{-((x_*+s-t_*)/\epsilon _*)^2}}{x_*+s-t_*}. \end{align}

In the steady case ( $\mathrm{d}C_y/\mathrm{d}s=0$ ), the time integral vanishes, and the term in the first brackets evaluated for the lower bound $s=0$ cancels the initial condition in (3.4). Thus, the only remaining term is the first term evaluated at the upper integration bound, which represents the induced velocity of the bound vortex, e.g. one recovers the solution of the steady problem.

When combining the contributions of the vorticity initial condition (3.4) and the vorticity due to unsteady forcing (3.6) to the normal induced velocity at the actuator point as $v_*(0,0,t_*)=v^{\textit{IC}}_*(0-t_*,0)+v^{\textit{us}}_*(0,0,t_*)$ , we obtain the solution

(3.8) \begin{align} v_*(0,0,t_*) = - \int _0^{t_*} \frac {1}{4\pi } \frac {\mathrm{d}C_y(s)}{\mathrm{ds}} \varphi (t_*-s) \: \mathrm{d}s, \end{align}

where we exploit the fact that $\varphi (x_*+s-t_*)$ for $x_*=0$ is an anti-symmetric function and $\varphi (0)=0$ . This solution is a Duhamel’s integral with the indicial response function $\varphi$ (Leishman Reference Leishman2002b ). The induced velocity at the actuator point at time $t_*$ due to the unsteady Gaussian forcing is then obtained as the superposition of all indicial responses to the forcing in the time interval $s\in [0,t_*]$ . Wagner (Reference Wagner1925) derived an indicial response function for the response of a thin aerofoil to a step change of the angle of attack, which has been shown by Garrick (Reference Garrick1938) to be directly related to Theodorsen’s function $C(k)$ in the frequency domain (Theodorsen Reference Theodorsen1935). In our case, the indicial response function for an unsteady Gaussian forcing is parameterised with the Gaussian kernel width $\epsilon _*$ , and from (3.7) it can be seen that the indicial response function is given by the shape representative for the induced velocity of a Lamb–Oseen vortex with core size $\epsilon _*$ located at the streamwise location $x_*=t_*-s$ .

3.2. Induced velocity due to streamwise forcing

The streamwise induced velocity along the wake centreline $(x_*, y_*=0)^\top$ is obtained by solving the streamwise component of the Biot–Savart integral given in (3.1) for a generic streamwise forcing $C_x(t)$ . We first focus again on the contribution from the second term of the unsteady vorticity solution $\omega ^{\textit{us}}_*$ (the time integral in (2.5)). By switching the order of temporal and spatial integration, it follows that

(3.9) \begin{align} & u^{\textit{us}}_*(x_*,0,t_*)= \nonumber \\ & \quad -\int \limits _0^{t_*} \frac {1}{2\pi }\int \limits _{-\infty }^{+\infty } \int \limits _{-\infty }^{+\infty } \frac {(-y_*^\prime )}{(x_*-x_*^\prime )^2+(-y_*^\prime )^2} \frac {-y^\prime _* C_x(s)}{\pi \epsilon _*^4} \text{e}^{-((x_*^\prime +s-t_*)^2+y_*^{\prime 2})/\epsilon _*^2} \mathrm{d}x_*^\prime \mathrm{d} y_*^\prime \mathrm{d} s \end{align}

and the two integrals in $x_*^\prime$ and $y_*^\prime$ can be solved analytically to obtain

(3.10) \begin{align} u^{\textit{us}}_*(x_*,0,t_*)= \int \limits _0^{t_*} \frac {C_x(s)}{4\pi } \bigg [ \frac {\text{e}^{-(x_*+s-t_*)^2/\epsilon _*^2}-1}{(x_*+s-t_*)^2} \bigg ] \: \mathrm{d} s. \end{align}

It should be noted that in the limit case of a vanishing spatio-temporal shift $x_*+s -t_*\rightarrow 0$ , the integrand tends to $-C_x(s)/(4\pi \epsilon _*^2)$ .

Furthermore, given the assumption of constant streamwise forcing $C_x \neq C_x(t_*)$ , one can recover the induced velocity solution of the steady two-dimensional problem derived by Martínez-Tossas et al. (Reference Martínez-Tossas, Churchfield and Meneveau2017). We derive here the solution valid for any $x_*$ , but restricted to $y_*=0$ . To this end, (3.10) is integrated for general $t_*$ but constant $C_x$ , which yields

(3.11) \begin{align} u^{\textit{us}}_*(x_*,0,t_*) &= \frac {C_x}{4\pi } \bigg [ \frac {1-\text{e}^{-(x_*+s-t_*)^2/\epsilon _*^2}}{(x_*+s-t_*)} -\frac {\sqrt {\pi }}{\epsilon _*} \textrm {erf}{\left (\frac {(x_*+s-t_*)}{\epsilon _*}\right )} \bigg ]_{s=0}^{s=t_*} . \end{align}

The steady-state solution along the curve $(x_*,y_*=0)^\top$ is obtained in the limit $t_*\rightarrow \infty$ using the rule of l’Hôpital:

(3.12) \begin{align} u^{st}_*(x_*,0) = \frac {C_x}{4\pi } \bigg [ \frac {1-\text{e}^{-x_*^2/\epsilon _*^2}}{x_*} -\frac {\sqrt {\pi }}{\epsilon _*} \textrm {erf}{\left (\frac {x_*}{\epsilon _*}\right )} - \frac {\sqrt {\pi }}{\epsilon _*} \bigg ]. \end{align}

In the case of $x_*\rightarrow \infty$ , the first term vanishes and the error function tends to one that leads to a limit value of $u^{st}_*(y_*=0)=-C_x/(2\sqrt {\pi }\epsilon _*)$ far downstream, which is exactly the result derived by Martínez-Tossas et al. (Reference Martínez-Tossas, Churchfield and Meneveau2017). Furthermore, at the actuator point the solution is exactly half of the limit value far downstream, i.e. $u^{st}_*(y_*=0)=-C_x/(4\sqrt {\pi }\epsilon _*)$ . Similar to the normal forcing case, the solution in (3.12) can also be used as an initial condition in order to account for a non-zero initial vorticity field $\omega ^{\textit{IC}}_*(x_*-t_*,y_*)$ due to streamwise forcing, i.e.

(3.13) \begin{align} u^{\textit{IC}}_*(x_*-t_*,0) = \frac {C_x(t_*=0)}{4\pi } \bigg [ \frac {1-\text{e}^{-(x_*-t_*)^2/\epsilon _*^2}}{x_*-t_*} -\frac {\sqrt {\pi }}{\epsilon _*} \textrm {erf}{\left (\frac {x_*-t_*}{\epsilon _*}\right )} - \frac {\sqrt {\pi }}{\epsilon _*} \bigg ]. \end{align}

The complete solution for the normal induced velocity along the wake centreline then follows by combining (3.10) and (3.13) as

(3.14) \begin{align} u_*(x_*,0,t_*)=u^{\textit{IC}}_*(x_*-t_*,0)+u^{\textit{us}}_*(x_*,0,t_*). \end{align}

4.1. Time domain solution: the root-finding problem

We now seek to obtain the time history of both the forcing and the velocity at the actuator point as a function of time using (3.5) and (3.14). We write the time history of the velocity at the actuator point with a single equation in terms of the flow angle, as done by Ning (Reference Ning2014) and Martínez-Tossas et al. (Reference Martínez-Tossas, Allaerts, Branlard and Churchfield2024). The tangent of the flow angle at the actuator point is defined by

(4.1) \begin{align} \tan (\phi (t_*))=\frac {v_*(0,0,t_*)}{1+u_*(0,0,t_*)}, \end{align}

where $u_*$ and $v_*$ can both be written in terms of $\phi$ ; thus, this is an implicit equation where the time history, $\phi (t_*)$ , is the only unknown. Rearranging (4.1) yields

(4.2) \begin{align} R(\phi )= v_*(0,0,t_*) \cos (\phi ) - (1+u_*(0,0,t_*)) \sin (\phi )=0. \end{align}

Given a guess for $\phi (t_*)$ as input, the following algorithm provides the steps to compute (4.2).

1. (i) Compute the angle of attack using the known pitch input, $\beta (t_*)$ :

   (4.3) \begin{align} \alpha (t_*)=\phi (t_*)+\beta (t_*). \end{align}

2. (ii) This effective angle of attack, which incorporates the effect of the shed vorticity, allows for the evaluation of the lift, $C_L(\alpha (t_*))$ , and drag, $C_D(\alpha (t_*))$ , coefficients from tabulated aerofoil data, which can be subsequently converted into the corresponding streamwise and normal forcing coefficients using the projection based on the flow angle given in (2.4):

   (4.4) \begin{align} C_x(t_*) & = -C_L(t_*) \sin (\phi (t_*)) + C_D(t_*) \cos (\phi (t_*)) , \\[-12pt] \nonumber \end{align}
   (4.5) \begin{align} C_y(t_*) & = C_L(t_*) \cos (\phi (t_*)) + C_D(t_*) \sin (\phi (t_*)). \\[12pt] \nonumber \end{align}

3. (iii) The integral equations for the induced velocities are solved via numerical integration:

   (4.6) \begin{align} u_*(0,0,t_*) & = u^{\textit{IC}}_*(0-t_*,0) + \int \limits _0^{t_*} \frac {C_x(s)}{4\pi } \bigg [ \frac {\text{e}^{-(s-t_*)^2/\epsilon _*^2}-1}{(s-t_*)^2} \bigg ] \: \mathrm{d} s, \\[-12pt] \nonumber \end{align}
   (4.7) \begin{align} v_*(0,0,t_*) & = v^{\textit{IC}}_*(0-t_*,0) + \int \limits _0^{t_*} -\frac {C_y(s)}{2\pi } \bigg [ \frac {\text{e}^{-(s-t_*)^2/\epsilon _*^2}}{\epsilon _*^2} + \frac {\text{e}^{-(s-t_*)^2/\epsilon _*^2}-1}{2(s-t_*)^2} \bigg ] \: \mathrm{d} s. \\[6pt] \nonumber \end{align}

4. (iv) Finally, the residual is computed:

   (4.8) \begin{align} R(\phi )= v_*(0,0,t_*) \cos (\phi ) - (1+u_*(0,0,t_*)) \sin (\phi ). \end{align}

A multidimensional root-finding algorithm then iteratively evaluates steps (i)–(iv) to compute the $\phi (t_*)$ that solves $R(\phi (t_*))=0$ . The resulting $\phi (t_*)$ is the full solution to the problem that provides the time history of forcing at the actuator point. The algorithm above is implemented based on the python scipy.optimize.root() function (Virtanen et al. Reference Virtanen2020) using the derivative-free spectral algorithm for nonlinear equations (DF-SANE) by La Cruz et al. (Reference La Cruz, Martínez and Raydan2006).

4.2. Frequency domain solution: the closed-loop transfer function

The time domain solution from § 4.1 can be compared with results from ALM-LES. However, it is insightful to also formulate the feedback problem in the frequency domain since it allows for a concise analysis of the Gaussian kernel width’s influence on the phase and magnitude modulation of the unsteady loading through a closed-loop transfer function. This connection is analogous to the one that can be obtained between Theodorsen’s frequency domain solution for unsteady aerofoil motion (Theodorsen Reference Theodorsen1935) and Wagner’s time domain solution (Wagner Reference Wagner1925) as shown by Garrick (Reference Garrick1938).

In order to derive the transfer function from the quasi-steady to the unsteady lift, we explicitly need to linearise the projection of the lift and drag coefficients onto the streamwise and normal directions (2.4), which is not necessary for the time-domain solution. When doing so, one obtains $C_y \approx C_L + C_D \phi$ , and since $C_D$ is at least an order of magnitude smaller than $C_L$ in the operating regime of interest for this study, it is consistent to neglect the drag because $C_D \phi$ is a higher-order term. The linearised lift coefficient is then described by its value at the linearisation point, $C_L(\beta ^0)$ , and the corresponding lift slope $\mathrm{d}C_L/\mathrm{d}\alpha (\beta ^0)$ . For a given sinusoidal pitch input signal $\beta (t_*)=\beta ^0+\Delta \beta \sin (2kt_*)$ with amplitude $\Delta \beta$ , the quasi-steady and unsteady lift coefficients are then, after the decay of any transient, given by

(4.9) \begin{align} C_L^{qs}(t_*;k) &= C_L(\beta ^0) + \Delta \beta \sin (2kt_*) \left .\frac {\mathrm{d}C_L}{\mathrm{d}\alpha }\right \rvert _{\beta ^0}, \\[-12pt] \nonumber \end{align}

(4.10) \begin{align} C_L^{\textit{us}}(t_*;k) &= C_L(\beta ^0) + \Delta \beta |G(k)| \sin (2kt_*+\angle G(k)) \left .\frac {\mathrm{d}C_L}{\mathrm{d}\alpha }\right \rvert _{\beta ^0}, \\[6pt] \nonumber \end{align}

where $|G|$ and $\angle G$ are the magnitude and phase of the corresponding closed-loop transfer function $G$ mapping from the quasi-steady to the unsteady lift. This transfer function can be derived from the system’s block diagram shown in figure 2, where we invoke the small-angle approximation for the flow angle $\phi \approx v_*$ . Calculating the Laplace transform of the individual blocks, the transfer function follows from block diagram algebra. It should be noted that the Laplace transform has the convenient property that the convolution of two functions in the time domain reduces simply to a multiplication in the Laplace domain. Denoting the Laplace transform of the indicial response function $\varphi$ for the normal induced velocity defined in (3.7) as $\varPhi (ik)$ , the closed-loop transfer function for the feedback problem then reads as

(4.11) \begin{align} G(k) = \frac {\Delta C^{\textit{us}}_L(k)}{\Delta C^{qs}_L(k)} = \frac {1}{1-2ki \left .\frac {\mathrm{d}C_L}{\mathrm{d}\alpha }\right \rvert _{\beta ^0} \varPhi (ik)}, \end{align}

with $i$ denoting the imaginary unit (not to be confused with the streamwise unit vector $\boldsymbol{i}$ ). The full expression for the Laplace transform of the indicial response function $\varphi$ is given by

(4.12) \begin{align} \varPhi (ik) =& \frac {1}{16\pi } \bigg [ 2 \mathbb{\gamma } - 2\pi \,\textrm {erfi}(ik \epsilon _*) - 2\log (1/\epsilon _*^2) \nonumber \\ &+4\log (2ki) + (2\epsilon _*ki)^2 \,_2F_2((1,1);(3/2,2);(ik\epsilon _*)^2) \bigg ], \\[6pt] \nonumber \end{align}

where $\mathbb{\gamma }$ is Euler’s constant and $_pF_q(a;b;z)$ is the generalised hypergeometric function. Unsteady aerodynamics problems are commonly classified using the reduced frequency $k=\pi f c/U_\infty$ , which is based on the semi-chord length (Leishman Reference Leishman2002b ). Recalling that the non-dimensional variables in this work are based on the chord length and the free-stream velocity, it follows that the angular frequency is given by $\sigma _*=2k$ , which explains the added factor of 2 when writing the frequency in terms of the reduced frequency $k$ .

Figure 2. Schematic block diagram of the feedback problem for the unsteady lift experienced by a pitching aerofoil with $C_D=0$ at a linearised operating point $\beta ^0$ . The convolution operator is denoted by $\ast$ and $\varphi$ is the indicial response function for the normal induced velocity (defined in (3.7)).

Figure 3. The closed-loop transfer function, $G(k)$ , for a pitching flat plate, ( $\mathrm{d}C_L/\mathrm{d}\alpha =2\pi$ ), in the complex plane. The function is shown for fifteen different Gaussian kernel widths. The Theodorsen function, $T^{\beta }_C(k)=C(k)$ , the circulatory component of Theodorsen’s transfer function, $T_C(k)$ , and Theodorsen’s complete transfer function, $T(k)$ , are shown for reference, where the aerofoil pitches around the quarter-chord point ( $a=-1/2$ ).

4.3. Relation between $G(k)$ and Theodorsen theory

Theodorsen (Reference Theodorsen1935) derived a transfer function $T(k)$ from quasi-steady to unsteady lift for a flat plate aerofoil modelled as a line distribution of point vortices. For an aerofoil in pure pitching motion, it consists of two non-circulatory (apparent mass) $T_{\textit{NC}}(k)$ and two circulatory $T_C(k)$ contributions given by

(4.13) \begin{align} T(k) = \frac {\Delta C^{\textit{us}}_L(k)}{\Delta C^{qs}_L(k)} = T_{\textit{NC}}(k) + T_{C}(k) &= \left [T^{\dot {\beta }}_{\textit{NC}}(k) + T^{\ddot {\beta }}_{\textit{NC}}(k) \right ] + \left [T^{\beta }_C(k) + T^{\dot {\beta }}_C(k) \right ]\nonumber \\ &= \frac {\left [ \pi i k + \pi a k^2 \right ] + 2\pi C(k) \left [ 1 + ik (\frac {1}{2}-a) \right ]}{2\pi }, \\[6pt] \nonumber \end{align}

where $a\in [-1,1]$ (from leading to trailing edge) is the location of the rotation axis of the pitching motion, and the complex Theodorsen function $C(k)$ is given by

(4.14) \begin{align} C(k) = \frac {K_1(ik)}{K_0(ik) + K_1(ik)}. \end{align}

Here, $K_0$ and $K_1$ denote modified Bessel functions of the second kind. The superscripts $(\boldsymbol{\cdot })^{\beta }$ , $(\boldsymbol{\cdot })^{\dot {\beta }}$ and $(\boldsymbol{\cdot })^{\ddot {\beta }}$ denote terms depending on the pitch angle, pitch rate and pitch acceleration, respectively. In the frequency domain, these dependencies translate to terms scaling proportionally to $1$ , $k$ and $k^2$ , respectively. The classical ALM, as it was proposed by Sørensen & Shen (Reference Sørensen and Shen2002) and as it is also used in this work, only models the $T^{\beta }_C(k)$ term. The remaining three terms, which are proportional to either $k$ or $k^2$ , are not captured by the ALM. It follows that the derived transfer function $G(k)$ is the equivalent of the Theodorsen function $C(k)$ consistent with the mathematical formulation of the ALM.

The ALM (for $\epsilon _* \in [0.125,16.0]$ ) and Theodorsen frequency domain solutions are compared in the complex plane in figure 3 and as a Bode plot in figure 4 up to a reduced frequency of $k=0.6$ . The magnitude of the transfer function is always $|G|\leqslant 1$ , whereas the phase can be both positive or negative depending on the considered kernel width and the reduced frequency. In general, it can be seen that smaller kernel widths lead to larger damping. The phase of $G(k)$ (see figure 4 b) first decreases with increasing $k$ before a characteristic phase inversion point is approached where the phase switches from negative to positive. This phase inversion point occurs for larger $k$ as the kernel width decreases. The magnitude $|C|(k)$ closely follows $|G|(k;\epsilon _*=0.375)$ for $k\lt 0.4$ . Adding the second circulatory term $T^{\dot {\beta }}_C(k)$ , which is proportional to the pitch rate, shows that the circulatory Theodorsen solution $|T_C|(k)$ departs from $|G|(k;\epsilon _*=0.375)$ for $k\gt 0.2$ . Finally, also adding the two non-circulatory terms leads to an additional deviation between the ALM and the complete Theodorsen solution $|T|(k)$ for $k\gt 0.4$ . The phase $\angle C(k)$ closely follows $\angle G(k;\epsilon _*=0.375)$ for $k\lt 0.4$ . In contrast to the magnitude, adding the second circulatory and eventually the two non-circulatory terms both leads to immediate deviations in phase compared with $\angle G(k;\epsilon _*=0.375)$ for any $k$ rather than the two distinct points of departure of $k=0.2$ and $k=0.4$ observed for the magnitude.

Figure 4. Magnitude (a) and phase (b) of $G(k)$ for a pitching flat plate, ( $\mathrm{d}C_L/\mathrm{d}\alpha =2\pi$ ), for fifteen different Gaussian kernel widths. The Theodorsen function, $T^{\beta }_C(k)=C(k)$ , the circulatory component of Theodorsen’s transfer function, $T_C(k)$ , and Theodorsen’s complete transfer function, $T(k)$ , are shown for reference, where the aerofoil pitches around the quarter-chord point ( $a=-1/2$ ).

These results suggest that when taking Theodorsen theory as a reference, the classical ALM can capture the unsteady damping of the lift accurately up to a reduced frequency of $k=0.2$ , whereas errors in phase are immediately encountered as soon as $k\gt 0$ . However, these phase errors are relatively small and at $k=0.2$ the difference in phase predicted by the ALM and the complete Theodorsen theory is $T(k=0.2)-G(k=0.2)\approx 19^\circ$ . The model validation shown later in § 7 will further support these identified bounds of the ability of the ALM to capture unsteady aerodynamic effects. For the purpose of further elucidating in § 7 the impact of the ALM not capturing the terms $T^{\dot {\beta }}_C$ , $T^{\dot {\beta }}_{\textit{NC}}$ and $T^{\ddot {\beta }}_{\textit{NC}}$ , we also define the extended closed-loop transfer function $G_{\textit{Ext}}(k)$ for the ALM. The extended ALM transfer function incorporates the three missing terms compared with Theodorsen theory in an ad-hoc approach and then is defined accordingly as

(4.15) \begin{align} G_{\textit{Ext}}(k) = \frac {\Delta C^{\textit{us}}_L(k)}{\Delta C^{qs}_L(k)} = \frac {\left [ \pi i k + \pi a k^2 \right ] + \dfrac {\mathrm{d}C_L}{\mathrm{d}\alpha } G(k) \left [ 1 + ik \left(\dfrac {1}{2}-a \right) \right ]}{\dfrac {\mathrm{d}C_L}{\mathrm{d}\alpha }}, \end{align}

where $G(k)$ effectively replaces the Theodorsen function $C(k)$ . It should be stressed that the current comparison between Theodorsen theory and the ALM is based on the classical ALM relying on point velocity sampling. An integral velocity sampling approach as introduced by Churchfield et al. (Reference Churchfield, Schreck, Martínez-Tossas, Meneveau and Spalart2017) could influence the transfer function and should be explored in future work.

5.1. Investigated cases

The range of investigated cases is chosen such that a range of operating points, types of actuation, reduced frequencies and the interplay between streamwise and normal forcing are explored. All cases are conducted using the NACA64-A17 aerofoil, which is, for example, also used for the outer part of the blade of the NREL 5 MW reference turbine (Jonkman et al. Reference Jonkman, Butterfield, Musial and Scott2009). The aerofoil’s tabulated lift and drag coefficients are shown in figure 5. Each investigated case is analysed for the set of Gaussian kernel widths $\epsilon _* \in \{0.25, 0.5, 1.0, 2.0, 4.0\}$ . The lower limit is chosen here on the order of the steady-state optimal kernel width (Martínez-Tossas et al. Reference Martínez-Tossas, Churchfield and Meneveau2017), while the upper limit corresponds to coarse grid ALM-LES simulations conducted for wind energy purposes. An overview of all cases and their acronyms is provided in table 1. In the following, the different cases are explained in more detail.

Table 1. Overview of the investigated cases. The initial condition for all cases is $\omega ^{\textit{IC}}=0$ , corresponding to a situation without initial forcing. The step magnitude is defined as the difference between the pitch angle after the step (given by the operating point $\beta ^0$ ) and the angle of attack, which results in approximately zero lift force for the cambered aerofoil $\Delta \alpha = \beta ^0-\alpha ^{C_L=0}$ ( $\alpha ^{C_L=0}\approx -4^\circ$ ). The periodic component of the actuation is of the form $\beta (t_*)=\Delta \beta \sin (2kt_*)$ . The step components contain a continuous frequency spectrum where the magnitude of the Fourier transform of the step scales as $|\hat {\beta }| \propto 1/k$ .

Figure 5. Tabulated lift, $C_L(\alpha )$ , and drag, $C_D(\alpha )$ , coefficients of the NACA64-A17 aerofoil.

In order to characterise the obtained theoretical solutions and their scaling with the kernel width, we begin by considering two cases with a step in pitch angle. Both cases start from an initial condition of zero vorticity, $\omega ^{\textit{IC}}=0$ , i.e. zero forcing. This choice ensures that there is no mismatch between the initial conditions used for the model and LES. The simulation start thus corresponds to a sudden step actuation to reach the operating pitch angle $\beta ^0$ for $t_*\geqslant 0$ . The considered operating points are $\beta ^0=0^\circ$ (A0-Cxy-S4) and $\beta ^0=8^\circ$ (A8-Cxy-S12). The same cases are also repeated without applying any normal forcing (A0-Cx-S4 and A8-Cx-S12). It should be noted that the frequency content of a step function is proportional to the inverse of the frequency ( $|\hat {\beta }|\propto 1/k$ ) and, thus, tests the model simultaneously across the entire frequency range. The step case A14-Cxy-S18 is only conducted for $\epsilon _*=0.25$ to obtain the largest forcing magnitudes. This case is used to perform a grid and domain size convergence study for the LES set-up.

The second set of cases considers the same initial step but additionally superposes continuous periodic pitch actuation with an amplitude of $\Delta \beta =3^\circ$ and reduced frequencies of $k=0.1$ , $k=0.2$ , $k=0.3$ and $k=0.6$ (indicated by the endings -k01, -k02, -k03 and -k06 of the case acronyms). The periodic pitch actuation is of the form $\beta (t_*)=\Delta \beta \sin (2kt_*)$ . It is applied at the operating point $\beta ^0=0^\circ$ (A0-Cxy-P3-k…). This second set of cases is well suited to test the predictions of the derived closed-loop transfer function for an isolated frequency component after the initial transient decayed.

5.2. Employed numerical code and simulation set-up

Numerical nonlinear reference results for the assessment of the developed model are obtained by means of LES using the open-source code AMR-Wind (available at https://github.com/Exawind/amr-wind). The AMR-Wind code solves the three-dimensional incompressible Navier–Stokes equations on block-structured Cartesian grids and is specifically designed for wind energy applications. For details on the numerical schemes employed in AMR-Wind, the reader is referred to Almgren et al. (Reference Almgren, Bell, Colella, Howell and Welcome1998), Sharma et al. (Reference Sharma2024) and Kuhn et al. (Reference Kuhn2025). For the purpose of verifying the developed model, which is based on the Euler equations, the governing equations are solved in the limit of infinite Reynolds number by setting the molecular viscosity to zero. Throughout the simulations, the inflow and the wake flow behind the actuator points remain laminar, removing the need for explicit turbulence modelling.

This work utilises AMR-Wind’s implementation of the ALM for straight translating wings. The AMR-Wind code does not allow for strictly two-dimensional simulations; rather, we simulate the aerofoil in a three-dimensional domain as an infinite wing, i.e. periodic boundary conditions are imposed in the spanwise direction. The lower and upper boundary condition in the normal direction is a slip wall ( $y$ -direction). In the streamwise direction ( $x$ direction), a constant inflow velocity of $\boldsymbol{u}=U_{\infty }\boldsymbol{i}$ is imposed at the upstream boundary (which in non-dimensional units corresponds to $\boldsymbol{u}_*=\boldsymbol{i}$ ), and a pressure outflow boundary condition is applied at the downstream boundary. The wing spans two actuator points located on the spanwise boundaries of the domain at which the actuator forces are applied. The forces are regularised with a two-dimensional isotropic Gaussian kernel of width $\epsilon _*$ . The velocity used for the actuator force calculation is the known inflow velocity (see Appendix C).

The convergence of LES results with domain size and grid resolution is studied in Appendix D for the A14-Cxy-S18 case using $\epsilon _*=0.25$ . Based on this convergence study, all simulations are conducted with a semi-domain height $L_{y_*}=256$ and an upstream fetch of $L^{\textit{upstr}}_{x_*}=256$ between the inflow boundary and the actuator point. The streamwise distance from the actuator point to the outlet is $L^{\textit{do}w\textit{nstr}}_{x_*}=T_*+L_{y_*}$ , where $T_*$ is the maximum simulated time. The grid resolution on the coarsest grid level is constant for all simulations; however, depending on the studied kernel width, four ( $\epsilon _*=4$ ) to eight ( $\epsilon _*=0.25$ ) refinement levels are added such that on the finest grid around the actuator point the resolution is always $\epsilon _*/\Delta x_*=8$ . The grid resolution between two adjacent grid levels varies with a constant factor of 2. The size of the finest grid level is dependent on the kernel width. Its extent in the $x$ and $y$ directions is $-4\epsilon _*\lt x_* \lt 16+4\epsilon _*$ and $-4\epsilon _*\lt y_* \lt 4\epsilon _*$ , respectively. Letting the index $m$ denote the grid level (the finest being $m=0$ ), the extent of all other refinement levels ( $m\gt 0$ ) in $x$ and $y$ is given by $-4\epsilon _* 2^m\lt x_* \lt T_*+4\epsilon _* 2^m$ and $-4\epsilon _* 2^m\lt y_* \lt 4\epsilon _* 2^m$ . The time-step size $\Delta t_{*,{\textit{LES}}}$ is chosen such that the Courant–Friedrichs–Levy (CFL) number defined on the finest grid level is $CFL\lt 0.5$ at all times. Step cases are simulated up to $T_*=128$ , whereas the step cases followed by periodic actuation are simulated up to $T_*=256$ . These long simulation times provide enough time for the actuator loading to reach the new steady state (step cases) or periodic limit cycle (step plus periodic cases) after the decay of the initial transient. The simulation time-dependent streamwise domain length ensures that no vorticity is truncated at the outlet. The simulation time for the periodic cases is chosen twice as long in order to collect a large number of actuation periods, even for the smallest frequency $k=0.1$ . The LES set-up is visualised in figure 6.

Figure 6. The LES domain and grid refinement structure for the Gaussian kernel width $\epsilon _*=0.25$ with the resolution on the finest grid level being $\epsilon _*/\Delta x=8$ . Note that the domain here is configured for a simulation time of $T_*=64$ in order to keep the aspect ratio limited for visual clarity. The actual simulations are performed for $T_*=128$ (steps) and $T_*=256$ (periodic). The red dot marks the actuator point. The orange dot marks the approximate downstream location of the start-up vortex at time $t_*=T_*$ .

5.3. Response to pitch step actuation

We start by probing the developed model for two pitch step actuation cases within the linear regime of the lift curve. Their respective operating points are $\beta ^0=0^\circ$ (A0-Cxy-S4) and $\beta ^0=8^\circ$ (A8-Cxy-S12). The initial condition is chosen as undisturbed free-stream flow, i.e. no forcing is active. Consequently, the two operating points correspond to effective angle of attack steps of $\Delta \alpha =4^\circ$ and $\Delta \alpha =12^\circ$ upon simulation start (the lift coefficient of the cambered aerofoil is approximately zero at $\alpha \approx -4^\circ$ ).

The vorticity distribution around the aerofoil and in the wake as predicted by the model (see Appendix B) with $\epsilon _*=0.25$ is shown in figure 7 for the larger pitch step case A8-Cxy-S12. Around the aerofoil, a bound vortex forms in response to the non-zero loading after $t_*\gt 0$ . The bound vortex is accompanied by the formation of a start-up vortex, which ensures conservation of the overall vorticity. Since the model assumes the advection speed in the wake to match the free-stream velocity, the start-up vortex is located at $x_*=12$ in the wake at the shown time instance $t_*=12$ . Furthermore, it is seen that the wake centreline $(x_*, y_*=0)^\top$ is the symmetry axis for both the vorticity due to streamwise and the vorticity due to normal forcing (figure 7 a,b). The vorticity due to streamwise forcing exhibits an anti-symmetric distribution around this symmetry axis, whereas the normal forcing contribution is distributed symmetrically around the wake centreline. These distributions stem from the symmetry properties of the forcing terms in the linearised vorticity transport equation(2.3).

Figure 7. Model solution for the vorticity field of the A8-Cxy-S12 case with the kernel width $\epsilon _*=0.25$ at time instance $t_*=12$ . The vorticity created by streamwise $C_x$ and normal $C_y$ forcing is shown in (a) and (b), respectively. The total resultant vorticity field is shown in (c). The aerofoil modelled by the Gaussian body force is shown in white with a black outline. The aerofoil’s quarter-chord point is located at the actuator point $(x^{\textit{Act}}_*, y^{\textit{Act}}_*)^\top =(0,0)^\top$ .

The formation of the start-up vortex and its subsequent advection downstream causes a dynamic angle of attack response. A comparison of the model and the LES solution for this response is shown in figure 8 for the small pitch step case A0-Cxy-S4. The start-up vortex induces velocity at the actuator point, which in turn alters the angle of attack. As it is advected further downstream, its influence on the angle of attack fades, which leads the angle of attack to converge towards the new steady-state operating point $\alpha =\beta ^0=0^\circ$ (figure 8 a). The initial modulation of the angle of attack depends on the kernel width since it sets the local strength of the applied forcing and, thus, also the strength of the start-up vortex. Furthermore, the core size of the start-up vortex is also set by the kernel width. The core size determines the time period between $t_*=0$ and the time instance when the angle of attack reaches its minimum. The minimum angle of attack is attained faster with a smaller kernel width and, thus, core size.

The model error compared with the LES reference is presented in figure 8(b). The very initial relative angle of attack error is as large as $3.4\,\%$ and then subsequently drops to the order of tenths of a percent. The relative error starts approaching zero after the actuator point no longer lies within the vortex core region of the start-up vortex. The relative error is larger for smaller kernel widths as long as there is overlap between the vortex core and the actuator point. This is the case because the velocity induced by the bound and start-up vortex both reach their absolute maximum values in proximity to their respective core boundary, where the maximum values increase with decreasing kernel width. Thus, after the time instance $t_* \gt |x^{\textit{Act}}_*-t_*| \approx |0-\epsilon _*|=\epsilon _*$ the impact of the bound vortex on the trajectory of the start-up vortex as well as the impact of the start-up vortex on the angle of attack at the actuator point start to decrease. Hence, the effect of modelling errors also decreases as soon as $t_*\approx \epsilon _*$ . The larger pitch step case A8-Cxy-S12 shown in figure 9 follows the previous observations made apart from the relative error evolution for $\epsilon _*=0.25$ . In this case, the error also persists after the vortex core no longer overlaps with the actuator point and only fades at $t_*\approx 10$ . The larger pitch step combined with the smallest kernel width makes for a stronger bound vortex, increasingly challenging the model, which is based on a linearisation around the vorticity-free background flow. The detailed impact of the error due to the linearisation will be explored in the following § 6.

Figure 8. Comparison of the model solution and LES reference for the initial time evolution of the angle of attack (a) and the corresponding relative error (b) for the pitch step actuation case A0-Cxy-S4. Reference data points are only shown in steps of $4\Delta t_{*,{\textit{LES}}}/\epsilon _*$ .

Figure 9. Comparison of the model solution and LES reference for the initial time evolution of the angle of attack (a) and the corresponding relative error (b) for the pitch step actuation case A8-Cxy-S12. Reference data points are only shown in steps of $4\Delta t_{*,{\textit{LES}}}/\epsilon _*$ .

Figure 10. Comparison of the model solution and LES reference for the initial time evolution of the streamwise (a) and normal (c) induced velocity and the streamwise (b) and normal (d) force coefficients for the pitch step actuation case A0-Cxy-S4. Reference data points are only shown in steps of $4\Delta t_{*,{\textit{LES}}}/\epsilon _*$ .

Figure 11. Comparison of the model solution and LES reference for the initial time evolution of the streamwise (a) and normal (c) induced velocity and the streamwise (b) and normal (d) force coefficients for the pitch step actuation case A8-Cxy-S12. Reference data points are only shown in steps of $4\Delta t_{*,{\textit{LES}}}/\epsilon _*$ .

We explore the two pitch step cases further by turning our attention towards the loading, which the model ultimately should predict. The model is able to accurately reproduce the time evolution of the streamwise, $C_x$ , and normal, $C_y$ , force coefficients for both the small (figure 10 b,d) and large (figure 11 b,d) pitch step cases. In particular, the transient drop in normal forcing and peak in streamwise forcing caused by the start-up vortex are captured. It is also worthwhile to study the induced velocities at the actuator point, which induce the change in flow angle and, thus, angle of attack (panels aand c of figures 10 and 11). The model matches the LES prediction for the normal induced velocity, $v_*$ , for both pitch step magnitudes, but discrepancies are apparent for the streamwise induced velocity, $u_*$ , which become more pronounced for smaller kernel widths and the larger pitch step (compare figures 10 a and 11 a). This discrepancy also matches the previous observation that the angle of attack error of the A8-Cxy-S12 case for $\epsilon _*=0.25$ persists for a longer time. Since we deal with small perturbation velocities, the flow angle is small and, thus, approximately given by $\phi \approx v_*$ . Thus, the error in streamwise induced velocity does not affect the quality of the streamwise and normal force coefficients as long as the small perturbation assumption is justified. In fact, $u_*$ is a measure of the system’s nonlinearity and is relevant for the flow angle calculation given larger induced velocities. However, larger induced velocities start to challenge the linearisation of the vorticity transport equation for the model derivation. Thus, the mismatch for $u_*$ and its increase for stronger forcing is consistent. In § 6 we show that the observed streamwise velocity residual is explained by the vorticity residual obtained from the difference of the model (linear) and the LES reference (nonlinear) solution. However, the two pitch step cases show that the developed model can predict the loading on the aerofoil in the linear regime of the lift curve despite the mismatch in $u_*$ . The cases also illustrate the importance of the kernel width. The largest studied kernel width of $\epsilon _*=4.0$ barely modifies the normal force coefficient, while the kernel width $\epsilon _*=0.25$ , which is of the order of the steady-state optimal value, results in a transient drop of more than $50\,\%$ (see figure 11 d).

Motivated by the mismatch of the streamwise induced velocity, we also conduct the pitch step cases with an absent normal force $C_y=0$ (the cases labelled A0-Cx-S4 and A8-Cx-S12). The resulting evolution of $u_*$ and the corresponding absolute model error are shown in figure 12 for the large pitch step case. In the absence of normal forcing, the previously observed mismatch between the model and the LES reference vanishes, and the model error after the initial time steps is 2 orders of magnitude smaller than $u_*$ . This highlights the impact of the assumptions made when linearising the vorticity transport equation for the model derivation. Neglecting the nonlinear transport terms, which include the induced velocity, causes the model to miss the normal displacement of the start-up vortex as it is advected in the induced normal velocity field of the bound vortex, which is maximal at the wake centreline $(x_*,y_*=0)$ . This lack of vertical displacement is also apparent in the previously shown model solution for the vorticity field (figure 7). Switching off the normal forcing, $C_y$ , removes its anti-symmetric contribution to the induced velocity field. Thus, as long as the streamwise induced velocity remains small, the model can now correctly predict $u_*$ too. It can also be noted that the induced streamwise velocity after the decay of the transient tends towards the steady-state limit given by $u_*=-C_x/(4\sqrt {\pi }\epsilon _*)$ , which was also derived by Martínez-Tossas et al. (Reference Martínez-Tossas, Churchfield and Meneveau2017).

Figure 12. Comparison of the model solution and LES reference for the time evolution of the streamwise induced velocity (a) and the corresponding absolute error (b) for the pitch step actuation case A8-Cx-S12. Reference data points are only shown in steps of $64\Delta t_{*,{\textit{LES}}}/\epsilon _*$ .

5.4. Response to periodic pitch actuation

After characterising the developed model with pitch steps, we now move on towards signals that are more present in reality, i.e. smooth signals and, in particular, periodic signals with sinusoidal shape. The solution to the latter can also be used to construct the solution for more complex forcing signals via their representation in terms of a Fourier series. In the following, we consider the four cases A0-Cxy-P3-k01, A0-Cxy-P3-k02, A0-Cxy-P3-k03 and A0-Cxy-P3-k06, which initially possess the same pitch step as the A0-Cxy-S4 case but then proceed with a continuous sinusoidal pitch actuation of amplitude $\Delta \beta =3^\circ$ . The four different cases vary in terms of the actuation frequency; namely, the four reduced frequencies $k=0.1$ , $k=0.2$ , $k=0.3$ and $k=0.6$ are considered.

Representative angle of attack evolutions obtained for these four reduced frequencies are shown for the case of $k=0.3$ in figure 13. The angle of attack now varies periodically, where its initial evolution shows a transient due to the influence of the start-up vortex, which initially is still present. After the transient, the angle of attack attains a limit cycle. For reference, figure 13 shows horizontal lines indicating the angle of attack amplitude predicted by the model closed-loop transfer function (given by $|G(k)|\Delta \beta$ ). The frequency domain solution in terms of the transfer function is able to predict the limit cycle amplitude modulation since the peaks of the angle of attack time series converge towards the horizontal lines. Therefore, the effect of increasingly smaller kernel widths is a successive reduction of the angle of attack amplitude compared with the quasi-steady reference ( $\Delta \beta =3^\circ$ ) and a successive increase in phase lag. For the given aerofoil and reduced frequency of $k=0.3$ , the reduction in effective angle of attack amplitude is as large as $35\,\%$ for the smallest studied kernel width $\epsilon _*=0.25$ , which is of the order of the steady-state optimal kernel width. The evolution of the induced velocities and the force coefficients corresponding to the previous angle of attack evolution of the A0-Cxy-P3-k03 case are shown in figure 14. Similar to the step cases, the model cannot predict the correct streamwise induced velocity, $u_*$ , in the presence of combined streamwise and normal forcing, where the error increases with decreasing kernel width.

Figure 13. Comparison of the model solution and LES reference for the initial time evolution of the angle of attack for the periodic pitch actuation case A0-Cxy-P3-k03. Reference data points are only shown in steps of $4\Delta t_{*,{\textit{LES}}}/\epsilon _*$ . The horizontal coloured dashed lines indicate the unsteady angle of attack amplitudes $|G(k)|\Delta \beta$ predicted by the model closed-loop transfer function $G(k)$ .

Figure 14. Comparison of the model solution and LES reference for the initial time evolution of the streamwise (a) and normal (c) induced velocity and the streamwise (b) and normal (d) force coefficients for the periodic pitch actuation case A0-Cxy-P3-k03. Reference data points are only shown in steps of $4\Delta t_{*,{\textit{LES}}}/\epsilon _*$ .

After the decay of the initial transient, it is convenient to represent the force coefficients in terms of hysteresis plots (figure 15). The grey reference lines in this plot show the quasi-steady force coefficients, i.e. the forcing as it would be obtained when the instantaneous angle of attack matches the instantaneous pitch angle. However, the velocity induced by the shed vorticity in the wake causes a non-zero flow angle, and it is $\alpha \neq \beta$ , which causes the unsteady force coefficients to differ from the quasi-steady ones both in terms of magnitude and phase. These differences between unsteady and quasi-steady force coefficients depend on the kernel width. For the shown case of $k=0.3$ in figure 15(b), it can be seen that smaller kernel widths cause stronger hysteresis, e.g. the minor axis of the ellipse for $C_y$ is larger (phase modulation). Furthermore, the major axis of the ellipse becomes increasingly tilted towards the horizontal axis, i.e. the force amplitude is damped compared with the quasi-steady case (amplitude modulation). For the streamwise forcing, $C_x$ , smaller kernel widths also lead to more pronounced hysteresis as observed for $C_y$ (figure 15 a). However, the trend for the amplitude modulation is different compared with $C_y$ in the sense that the unsteady solution for $C_x$ , compared with the quasi-steady reference, shows an increase in amplitude rather than a decrease. This amplitude increase is larger for smaller kernel widths because for these cases, the normal induced velocity reaches larger values (figure 14 c), and thus, the same follows for the flow angle $\phi$ . Since the streamwise force coefficient in the linear regime of the lift curve is dominated by the $C_L$ contribution, it attains its maximum when the product $-C_L(\alpha ) \sin (\phi )$ becomes maximal, which in addition to large $C_L$ also requires flow angles close to their maximal negative value since $C_L$ and $v_*$ are approximately half a period out of phase (compare panels c and d of figure 14).

Figure 15. Comparison of the model solution and LES reference for the hysteresis of the streamwise (a) and normal (b) force coefficient (periodic pitch actuation case A0-Cxy-P3-k03). The quasi-steady references $C_x(\beta )$ and $C_y(\beta )$ are shown in grey.

Given the dominance of the normal forcing in the linear regime of the lift curve for which the present model is developed, the observations regarding the phase and amplitude modulations made from the time series and hysteresis curves can be concisely summarised in a Bode plot. Figure 16(a,b) shows the Bode magnitude and Bode phase plot of the model closed-loop transfer function $G(k)$ , which was derived assuming that $C_D=0$ . Here $G(k)$ is computed using the linearised value of the lift slope of the NACA64-A17 aerofoil at the operating point $\beta ^0=0^\circ$ , matching the operating point of the periodic actuation cases. The model solution is presented for the set of kernel widths $\epsilon _* \in \{0.25, 0.5, 1.0, 2.0, 4.0\}$ in the range $k\in [0,0.75]$ . In addition, the phase and amplitude modulation are extracted from the last period of the LES reference data for the cases A0-Cxy-P3-k01 ( $k=0.1$ ), A0-Cxy-P3-k02 ( $k=0.2$ ), A0-Cxy-P3-k03 ( $k=0.3$ ) and A0-Cxy-P3-k06 ( $k=0.6$ ), which are shown as dots in the plot. The transfer function solution matches the LES reference for all four reduced frequencies.

Figure 16. The magnitude (a) and phase (b) of the model closed-loop transfer function $G(k)$ for the NACA64-A17 aerofoil at the linearisation point $\beta ^0=0^\circ$ . Dots indicate the corresponding amplitude and phase modulations extracted from the last period of the LES reference cases A0-Cxy-P3-k01, A0-Cxy-P3-k02, A0-Cxy-P3-k03 and Cxy-P3-k06.

For reduced frequencies $k\leqslant 0.3$ , magnitude and phase trends are monotone for varying kernel widths at a given constant frequency. However, the magnitude and phase variation with frequency for a constant kernel width are not monotone. The magnitude reaches a minimum, after which it again approaches $|G|=1$ for further increasing $k$ . The convergence to $|G|=1$ happens for smaller $k$ , the larger the kernel width, and for the shown relevant frequency range, it is only visible for $\epsilon _*=4.0$ . In this limit, the unsteady lift magnitude again approaches the quasi-steady solution. Also, the phase tends towards a common limit of $\angle G=0^\circ$ for all kernel widths as $k\rightarrow \infty$ , although this convergence happens for $k\gt 0.75$ except for the largest studied kernel width. Consequently, for a chosen kernel width representing the aerofoil, the maximum reduction of the quasi-steady lift amplitude $\text{min}(|G(k; \epsilon _*)|)$ and maximum phase lag $\text{min}(\angle G(k; \epsilon _*))$ occur for an intermediate reduced frequency $0\lt k\lt \infty$ . Furthermore, the maximum phase lag occurs for smaller $k$ than does the maximum amplitude reduction. Within the range of relevant $k$ , smaller kernel widths lead to larger damping. Comparing $\epsilon _*=0.25$ , which is of the order of the steady-state optimal kernel width, with $\epsilon _*=4$ , which corresponds to a value encountered for coarse grid ALM-LES for wind energy purposes, one would miss about $\Delta (|G|)\approx 0.3$ of damping and $\Delta (\angle G) \approx 23^\circ$ of phase lag at a reduced frequency of $k=0.3$ . This difference is the unsteady equivalent of the steady-state lift error observed for suboptimal large kernel widths, which motivated the development of smearing corrections (Martínez-Tossas & Meneveau Reference Martínez-Tossas and Meneveau2019; Meyer Forsting et al. Reference Meyer Forsting, Pirrung and Ramos-García2019a ; Dağ & Sørensen Reference Dağ and Sørensen2020; Kleine et al. Reference Kleine, Hanifi and Henningson2023a ).

Figure 17. Model solution for the total resultant vorticity field of the A0-Cxy-P3-k03 case at time instance $t_*=32$ . Panels (a–e) show the solution for the Gaussian kernel widths $\epsilon \in \{4.0, 2.0, 1.0, 0.5, 0.25\}$ . The aerofoil modelled by the Gaussian body force is shown in white with a black outline. The aerofoil’s quarter-chord point is located at the actuator point $(x^{\textit{Act}}_*, y^{\textit{Act}}_*)^\top =(0,0)^\top$ .

The source of the unsteady error due to large kernel widths can also be made apparent by comparing the topology of the shed vorticity in the wake for varying $\epsilon _*$ (figure 17). The regularisation with the Gaussian kernel conserves both the integrated force as well as the integrated vorticity. That is assuming the magnitude of the point force would not be affected by the induced velocities, which in reality it is. But even neglecting this effect, the conservation on an integral level can still be achieved by locally very different vorticity distributions. The vorticity for small kernel widths is more localised, and thus, integral conservation directly implies locally higher vorticity magnitudes in proximity of the aerofoil and the wake centreline. In turn, large kernel widths spread the vorticity across a larger area, reducing its local magnitude. Since the induced velocity scales with the inverse of the distance outside of the vortex core, this directly implies smaller induced velocities at the actuator point. Finally, the phase error is also apparent from the vorticity distributions. To this end, one can track the vertical band of zero vorticity, which is located between the bound vortex (blue) and the neighbouring positive vorticity patch (red) directly downstream of the aerofoil. Its location is closest to the aerofoil for $\epsilon _*=0.25$ and then shifts continuously farther downstream for larger kernel widths, i.e. the larger the kernel width, the further ahead the phase.

Summarising the model verification presented in this section, we conclude that the developed model matches the ALM-LES predictions for the unsteady loading on the aerofoil within the considered operating regime, which is the linear part of the lift curve. The step and periodic actuation cases demonstrated the frequency-dependent impact of the Gaussian kernel width on the unsteady loading, both in terms of lift reduction and added phase lag.

7.1. Validation with inviscid (Euler) CFD

We first compare the model with numerical reference data from Motta, Guardone & Quaranta (Reference Motta, Guardone and Quaranta2015). They employ a finite-volume Euler solver for moving, overset, multi-block grids to compute the loads on symmetric aerofoils in response to sinusoidal pitching motion around the quarter-chord point ( $a=-0.5$ ). The pitching motion is characterised by an operating point of $\beta ^0=0^\circ$ , a pitch amplitude of $\Delta \beta =1^\circ$ and considers the range of reduced frequencies $k\in [0.01, 0.75]$ . The flow is approximately incompressible with a Mach number of $Ma=0.117$ . Four different NACA aerofoils (0004, 0012, 0018, 0024) were simulated to study the impact of aerofoil thickness on the unsteady loads and the results were compared with Theodorsen theory.

We compare in figure 22 the ALM transfer function, $G(k)$ , for the unsteady lift with the results from Motta et al. (Reference Motta, Guardone and Quaranta2015) obtained for their thinnest studied aerofoil (NACA 0004). In addition, the Theodorsen transfer function, $T(k)$ , is shown. It should be noted that Motta et al. (Reference Motta, Guardone and Quaranta2015) chose to present their results in terms of the unsteady lift coefficient normalised with the pitch amplitude ( $\Delta C^{\textit{us}}_L(k)/\Delta \beta$ ) and, thus, the transfer function $G(k)$ and $T(k)$ need to be multiplied with the aerofoil lift slope ( $\mathrm{d}C_L/\mathrm{d}\alpha =2\pi$ for the flat plate) in order to obtain the appropriate scaling. Furthermore, it should be noted that the transfer function $G(k)$ for the ALM is also computed based on the flat plate lift slope to be consistent with Theodorsen theory. The reference data for the NACA 0004 aerofoil closely match the magnitude and phase of $T(k)$ . Because of the good match between the thin-aerofoil CFD results and $T(k)$ , the conclusions reached in § 4.3 regarding the limits of the ALM’s applicability to model unsteady aerodynamics apply to the validation with the CFD data. A Gaussian kernel width of $\epsilon _*\approx 0.4$ closely matches the magnitude of the unsteady lift up to $k\lt 0.2$ , whereas deviations in phase occur as soon as $k\gt 0$ .

Figure 22. Magnitude (a) and phase (b) of the unsteady lift coefficient normalised with the actuation pitch amplitude $\Delta \beta$ as function of reduced frequency $k$ . The predictions of the model frequency domain solution assuming the lift slope of a flat plate are compared with the unsteady lift predicted by Theodorsen theory and inviscid CFD results from Motta et al. (Reference Motta, Guardone and Quaranta2015) for a NACA 0004 aerofoil.

Figure 23. Magnitude (a) and phase (b) of the unsteady lift coefficient normalised with the actuation pitch amplitude $\Delta \beta$ as function of reduced frequency $k$ . The predictions of the extended model frequency domain solution assuming the lift slope of a flat plate are compared with the unsteady lift predicted by Theodorsen theory and inviscid CFD results from Motta et al. (Reference Motta, Guardone and Quaranta2015) for a NACA 0004 aerofoil.

In § 4.3 we also showed that the three remaining terms modelled by Theodorsen theory, but neglected by the ALM, can be incorporated into an extended ALM transfer function $G_{\textit{Ext}}(k)$ (4.15). In figure 23, the previous comparison with the reference data from Motta et al. (Reference Motta, Guardone and Quaranta2015) for the NACA 0004 aerofoil is repeated for $G_{\textit{Ext}}(k)$ . It can be seen that $G_{\textit{Ext}}(k;\epsilon _*=0.4)$ closely reproduces the prediction from Theodorsen theory and, thus, also leads to an accurate magnitude and phase prediction compared with the CFD reference. While the objective of the present study is to reproduce the behaviour of the classical ALM with a simple model, this result highlights avenues for future work to incorporate the non-circulatory terms, $T^{\dot {\beta }}_{\textit{NC}}$ and $T^{\ddot {\beta }}_{\textit{NC}}$ , and the pitch rate dependent circulatory term, $T^{\dot {\beta }}_C$ , into the ALM formulation in order to increase the bounds of the ALM’s applicability to unsteady aerodynamics beyond $k=0.2$ .

The presented comparison with inviscid CFD and Theodorsen theory leads to recommendations regarding the optimal kernel width in order to accurately capture unsteady aerodynamic effects with the classical ALM based on point velocity sampling. We find that the reference data are closely matched up to $k\lt 0.2$ for a kernel width in the range $\epsilon _*\in [0.3,0.5]$ . Based on the validation of the extended ALM transfer function $G_{\textit{Ext}}(k)$ , a kernel width of $\epsilon _*\approx 0.4$ provides a good match with Theodorsen theory and CFD throughout the frequency range $k\in [0,0.75]$ . These identified values for the optimal unsteady kernel width are slightly higher than the range of the optimal kernel width found for steady conditions $\epsilon _{*,Opt}\in [1/8,1/4]$ (Shives & Crawford Reference Shives and Crawford2013; Martínez-Tossas et al. Reference Martínez-Tossas, Churchfield and Meneveau2017). However, we note that when using a Gaussian kernel with different widths in the chord and thickness directions, Martínez-Tossas et al. (Reference Martínez-Tossas, Churchfield and Meneveau2017) found that the optimal width along the chord direction is approximately $\epsilon _* \approx 0.3$ for thin aerofoils. Since the chord is the principal direction in the unsteady problem, it is reasonable for the optimal width to scale accordingly. Given that practical ALM-LES for wind turbine/farm modelling typically operate with kernel widths that are an order of magnitude larger, the difference between the identified optimal kernel width for steady and unsteady aerodynamics is rather small. Nevertheless, future work should explore if this difference persists when considering an advanced ALM formulation based on integral velocity sampling as proposed by Churchfield et al. (Reference Churchfield, Schreck, Martínez-Tossas, Meneveau and Spalart2017).

Figure 24. Magnitude (a) and phase (b) of the unsteady lift coefficient normalised with the actuation pitch amplitude $\Delta \beta$ as function of reduced frequency $k$ . The predictions of the extended model frequency domain solution are shown for $\epsilon _*=0.4$ . The employed lift slope varies with the aerofoil thickness according to $\mathrm{d}C_L/\mathrm{d}\alpha =2\pi (1+0.77 d/c)$ . Results are compared with inviscid CFD results from Motta et al. (Reference Motta, Guardone and Quaranta2015) for the NACA 0004, NACA 0012, NACA 0018 and NACA 0024 aerofoils.

Lastly, Motta et al. (Reference Motta, Guardone and Quaranta2015) also studied the unsteady loads on NACA 0012, 0018 and 0024 aerofoils and showed that the increasing non-dimensional aerofoil thickness of up to $d/c=24\,\%$ challenges the flat plate assumption made by Theodorsen theory. Figure 24 shows that this equally applies to the ALM. Figure 24(a) shows that the steady lift magnitude, $|\Delta C^{\textit{us}}_L|(k=0)$ , continuously increases with the aerofoil thickness. In fact the term $\Delta C^{\textit{us}}_L(k=0)/\Delta \beta$ is the lift slope of the aerofoil, which can be approximated with the relation $\mathrm{d}C_L/\mathrm{d}\alpha =2\pi (1+0.77 d/c)$ (Motta et al. Reference Motta, Guardone and Quaranta2015). The only way the aerofoil thickness enters the ALM is via the lift(-slope) input. Thus, we use the relation between thickness and lift slope given above to determine varying lift slope inputs for $G_{\textit{Ext}}(k)$ in order to model the impact of increasing aerofoil thickness. Comparing the unsteady lift magnitude in figure 24(a) shows that the $G_{\textit{Ext}}(k)$ now more closely matches the CFD results in proximity of $k=0$ than when using the flat plate lift slope to ensure consistency with Theodorsen theory (compare with figure 23). However, these deviations are small compared with the deviations between the CFD and ALM results for larger aerofoil thicknesses of . The increased damping of the lift magnitude predicted by CFD is not captured by the ALM, in particular, as soon as $k\gt 0.05$ . Similar conclusions are found for the phase in figure 24(b). The CFD results predict a shift of the phase inversion point to larger $k$ as the thickness increases. This trend is only qualitatively captured by the ALM, but both the maximum phase lag and the reduced frequency corresponding to the phase inversion point are quantitatively underestimated.

7.2. Validation with experiment at $Re=32\,000$

Lastly, we validate the ALM predictions with the experimental data set from Ōtomo et al. (Reference tomo, Henne, Mulleners, Ramesh and Viola2020). They studied a NACA 0018 aerofoil with periodic smoothed triangular pitching kinematics at $Re=32\,000$ . Hence, this validation data set also probes the model assumption of inviscid flow, and furthermore, the model is tested for actuation signals that are more complex than mono-frequency sinusoids. The triangular pitching kinematics are characterised by an operating point of $\beta ^0=0^\circ$ , pitch amplitudes of $\Delta \beta \in \{4^\circ , 8^\circ \}$ and periods corresponding to the reduced frequencies $k\in \{0.22, 0.44, 0.66\}$ . The location of the pitching axis is the quarter-chord point ( $a=-0.5$ ). Ōtomo et al. (Reference tomo, Henne, Mulleners, Ramesh and Viola2020) compared their experimental results also to predictions from Theodorsen theory by considering the superposition of the individual Theodorsen solutions corresponding to the twenty leading Fourier harmonics.

The temporal evolution of the unsteady lift coefficient in response to the smoothed triangular pitching kinematics is shown for one period $T_{p*}$ in figure 25 and the corresponding maxima of the lift coefficient are summarised in figure 26. The experimental data and the Theodorsen solution are compared with the time-domain solution of the developed model. The model solution is computed for the kernel width $\epsilon _*=0.4$ , which was determined as the recommended kernel width for unsteady conditions in the previous validation section. We start by focusing on the case $k=0.22$ / $\Delta \beta = 4^\circ$ (figure 25 a), which can be seen as the limit case for the ALM’s applicability to unsteady aerodynamics given the preceding comparison with Theodorsen theory and the CFD results from Motta et al. (Reference Motta, Guardone and Quaranta2015) in § 7.1. The Theodorsen solution is also shown split into the non-circulatory and circulatory contributions to visualise their relative importance. The non-circulatory term proportional to the pitch acceleration is zero within the linear parts of the pitch actuation signal due to the constant pitch rate and, thus, only contributes close to the extreme values of the pitch. It can be seen that the ALM predicts a slightly smaller unsteady lift magnitude than Theodorsen theory, which can be attributed to the observation that the non-circulatory terms neglected by the ALM reduce the damping of the unsteady lift (also recall figure 4). Due to the inclusion of the non-circulatory terms, the Theodorsen theory predicts a positive phase shift with respect to the pitch actuation signal with the fundamental frequency of $k=0.22$ , whereas the ALM still predicts a small phase lag. The experimental data shows the overall largest damping and phase lag. Motta et al. (Reference Motta, Guardone and Quaranta2015) also studied the effect of including viscosity in their CFD study such that $Re=10^6$ for a NACA 0018 aerofoil. They observed a reduced maximum lift coefficient for the viscous results compared with the inviscid reference. Furthermore, experiments conducted by Mackowski & Williamson (Reference Mackowski and Williamson2015) for a NACA 0012 aerofoil subject to sinusoidal pitching at $Re=16\,600$ also showed stronger damping, i.e. smaller peak lift coefficient values, compared withTheodorsen theory for the pitch amplitudes $\Delta \beta \in \{2,4,8\}^\circ$ in the range $k\in [0,1]$ . The increased deviation between the inviscid ALM and Theodorsen solutions compared with the experimental data from Ōtomo et al. (Reference tomo, Henne, Mulleners, Ramesh and Viola2020) thus is in agreement with the observations made in the two aforementioned studies. It should be further pointed out that one should not conclude from the comparison with the experiment that the ALM prediction is more accurate than Theodorsen theory. The closer match between the ALM and the experiment is rather likely due to fortunate cancellation of the modelling errors due to neglecting viscosity, but in addition also the three terms $T^{\dot {\beta }}_{\textit{NC}}$ , $T^{\ddot {\beta }}_{\textit{NC}}$ and $T^{\dot {\beta }}_C$ that in turn are captured by Theodorsen theory.

Doubling the pitch amplitude (the case $k=0.22$ / $\Delta \beta = 8^\circ$ shown in figure 25 c) leads to doubled lift magnitudes predicted by ALM and Theodorsen theory, which is consistent with their linear nature. In contrast, the increase in lift magnitude seen in the experiment is only $60\,\%$ . Furthermore, phase predictions of the linear models are unaltered, whereas the experiment shows a gain in phase such that the unsteady lift is almost in phase with the pitch actuation. Further increasing the fundamental frequency of the pitch actuation to $k=0.44$ increases the importance of the non-circulatory terms, where the unsteady lift signals from Theodorsen theory and experiment now are ahead of the pitch actuation (figure 25 b,d).

Figure 25. One period of the unsteady lift coefficient in response to smoothed triangular pitching kinematics with $k=0.22$ (a,c) and $k=0.44$ (b,d) with actuation amplitudes of $\Delta \beta =4^\circ$ and $\Delta \beta =8^\circ$ . The predictions of the model time domain solution are shown for $\epsilon _*=0.4$ using the lift coefficient of a flat plate. For reference, the experimental results from Ōtomo et al. (Reference tomo, Henne, Mulleners, Ramesh and Viola2020) are shown (the shaded grey area denotes the 95 % confidence interval). In addition, the predictions of Theodorsen theory are shown, where the non-circulatory and circulatory contributions are also shown separately. The thick semi-transparent line indicates the time series of the pitch actuation signal. Dots mark the maxima and minima of the signals.

As a summary of the observations made above, the maximum lift coefficients for the ALM, Theodorsen theory and the experiment are shown as a function of the reduced frequency $k$ in figure 26. It can be concluded that for the case of $k=0.22$ / $\Delta \beta = 4^\circ$ both the ALM and Theodorsen theory still provide reasonable lift predictions despite the low Reynolds number $Re=32\,000$ of the experiment. In this light, the observed errors can be seen as an upper bound, given that for wind energy purposes, where the ALM finds its most use, chord-based Reynolds numbers are two orders of magnitude larger. It should be again stressed that the apparent better match between the ALM and experiment compared with Theodorsen theory is expected to stem from the lack of modelling the terms $T^{\dot {\beta }}_{\textit{NC}}$ , $T^{\ddot {\beta }}_{\textit{NC}}$ and $T^{\dot {\beta }}_C$ . Thus, the experimental validation affirms the conclusions drawn from the ALM-Theodorsen comparison (§ 4.3) and from the validation with inviscid CFD results (§ 7.1), namely that the ALM can capture unsteady aerodynamic effects up to $k\approx 0.2$ given that the Gaussian kernel width is close to the determined optimum $\epsilon _* \approx 0.4$ ; the aerofoil thickness is limited; and errors in phase of approximately twenty degrees can be accepted.

Figure 26. Maximum unsteady lift coefficient resulting from the smoothed triangular pitching kinematics with amplitude $\Delta \beta =4^\circ$ and $\Delta \beta =8^\circ$ . The predictions of the model time domain solution are shown for $\epsilon _*=0.4$ using the lift coefficient of a flat plate. For reference, the experimental results of Ōtomo et al. (Reference tomo, Henne, Mulleners, Ramesh and Viola2020) are shown (error bars denote the 95 % confidence interval). In addition, the predictions of Theodorsen theory are shown.