2.1. Oscillation of a thin aerofoil in uniform flow

The aerodynamic forces generated by the harmonic motion of an aerofoil in uniform flow were analysed by Theodorsen (Reference Theodorsen1935). Theodorsen developed an exact solution using the potential flow and the Kutta condition. The unsteady lift ( $L$ ) can be expressed as (Bisplinghoff, Ashley & Halfman Reference Bisplinghoff, Ashley and Halfman2013)

(2.1) \begin{equation} L = \pi \rho {b}^2 \left [\ddot {h} + U_{\infty } \dot {\alpha _g} - {b} a \ddot \alpha _g\right ] + 2 \pi \rho U_{\infty } b C(k)\left [\dot h + U_{\infty } \alpha _g + b\left (\frac {1}{2} - a\right ) \dot \alpha _g\right ], \end{equation}

where $\rho$ is the fluid density, $b$ is the half-chord length, $h$ is the vertical displacement, $U_{\infty }$ is the free stream velocity, $\alpha _g$ is the geometric angle of attack, $k$ is the reduced frequency based on semi-chord length ( $k = \varOmega c / 2 U_{\infty }$ ), $a$ is the distance between the rotation axis and the mid-chord, and $C(k)$ is Theodorsen’s function related to the Hankel functions,

(2.2) \begin{equation} C(k) = \frac {H_1^{(2)}(k)}{H_1^{(2)}(k)+iH_0^{(2)}(k)}, \end{equation}

where $H_n^{(2)}$ is a combination of Bessel function of first and second kind, $H_n^{(2)}=J_n - iY_n$ . This study focuses on the second term of (2.1), which corresponds to the circulatory component. Due to the absence of an explicit surface in ALM simulations, we do not expect significant effects from the the first term (apparent mass term).

The lift coefficient generated by a quasi-steady effect is related to the unsteady lift in (2.1) (considering only the circulatory term) by

(2.3) \begin{equation} C_l(\dot h, \alpha _g, \dot \alpha _g) = C(k) {C_l}_{\textit{QS}}(\dot h, \alpha _g, \dot \alpha _g), \end{equation}

where ${C_l}_{\textit{QS}}$ is the quasi-steady lift simplification, which is calculated as

(2.4) \begin{equation} {C_l}_{\textit{QS}} = 2 \pi \left [\dot h^* + \alpha _g + \left (\frac {1}{2} - a \right ) \dot \alpha _g^* \right ]\!, \end{equation}

where $\dot h^* = \dot h / U_{\infty }$ and $\dot \alpha _g^*$ = $\dot \alpha _g$ $b/U_{\infty }$ are the non-dimensional plunge and pitch rate, respectively. The relation in (2.3) indicates that the ratio of the unsteady lift and the quasi-steady lift is governed by Theodorsen’s function $C(k)$ , which reflects the wake’s influence. However, there is an additional simplification, which neglects the pitch rate in (2.4),

(2.5) \begin{equation} {C_l}_{\textit{SS}} = 2 \pi \Big [\dot h^* + \alpha _g^*\Big ], \end{equation}

commonly used in ALM simulations considering the steady-state conditions, although the neglect of the pitch rate in (2.5) results in a difference in the lift calculation. Later, we will focus on the effects of applying these two simplifications to unsteady aerodynamics.

2.2. ALM unsteady effects

2.2.1. Vorticity created by an unsteady lift force convoluted with a three-dimensional Gaussian function

The proposed theoretical model is based on the vorticity generated by actuator lines, extending the work of Forsythe et al. (Reference Forsythe, Lynch, Polsky and Spalart2015) and Martínez-Tossas & Meneveau (Reference Martínez-Tossas and Meneveau2019) to consider unsteady effects. For completeness, this section reproduces the derivation initially presented by Kleine (Reference Kleine2022). The inviscid vorticity equation, considering unsteady effects, is

(2.6) \begin{equation} \frac {\partial \boldsymbol{\omega }}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\omega } = \boldsymbol{\omega } \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u} + \boldsymbol{\nabla }\times \boldsymbol{f} . \end{equation}

We assume an undisturbed velocity in the $z$ -direction, $\boldsymbol{U}=(0,0,U_{\infty })$ . If only the lift force, $\boldsymbol{f}=(0,f_y,0)$ , along a line located at $(y,z)=(0,0)$ is considered, the linearisation of the vorticity equation considering small disturbances to $\boldsymbol{U}$ gives rise to (Kleine Reference Kleine2022)

(2.7) \begin{equation} \frac {\partial \boldsymbol{\omega }}{\partial t} + U \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\omega } = \boldsymbol{\nabla }\times \boldsymbol{f} , \end{equation}

which can be written as

(2.8) \begin{align} \frac {\partial \omega _x}{\partial t} + U_{\infty } \frac {\partial \omega _x}{\partial z} &= -\frac {\partial f_y}{\partial z} , \end{align}

(2.9) \begin{align} \frac {\partial \omega _z}{\partial t} + U_{\infty } \frac {\partial \omega _z}{\partial z} &= \frac {\partial f_y}{\partial x} , \end{align}

with the boundary condition $\boldsymbol{\omega } \rightarrow 0$ as $z \rightarrow -\infty$ . Assuming a periodic body force with the format $f_y(x,y,z,t)=\widehat {f_y}(x,y,z) e^{i \varOmega t}$ , the solutions of the first-order linear differential equations are

(2.10) \begin{align} \omega _x & = \frac {\exp {\left (i \varOmega \left (t - \frac {z}{U_{\infty }} \right )\right )}}{U_{\infty }} \int _{-\infty }^{z} -\frac {\partial \widehat {f_y}}{\partial z} \exp {\left (\frac {i \varOmega z}{U_{\infty }}\right )}\, {\rm d}z \nonumber\\ & = - \frac {\widehat {f_y}}{U_{\infty }} \exp {(i \varOmega t}) + \frac {\exp {\left (i \varOmega \left (t - \frac {z}{U_{\infty }} \right )\right )}}{U_{\infty }} \int _{-\infty }^{z} \frac {i \varOmega }{U_{\infty }} \widehat {f_y} \exp {\left (\frac {i \varOmega z}{U_{\infty }}\right )} \,{\rm d}z , \end{align}

(2.11) \begin{align} \omega _z &= \frac {\exp {\left (i \varOmega \left (t - \frac {z}{U_{\infty }} \right )\right )}}{U_{\infty }} \int _{-\infty }^{z} \frac {\partial \widehat {f_y}}{\partial x} \exp {\left (\frac {i \varOmega z}{U_{\infty }}\right )} \,{\rm d}z . \end{align}

We consider an oscillatory force convoluted with the Gaussian kernel such as

(2.12) \begin{equation} \widehat {f_y}(x,y,z)= -\frac {1}{\rho }(\widehat {F_l}\eta )(x) \, \eta (y) \, \eta (z) , \end{equation}

where

(2.13) \begin{equation} \eta (x) := \frac {1}{\pi ^{1/2}\varepsilon } \exp {\left ({-}\frac {x^{2}}{\varepsilon ^{2}}\right )} = \frac {1}{\varepsilon } \eta \left (\frac {x}{\varepsilon } \right ) \quad \mathrm{and}\quad \eta (x^*) := \frac {1}{\pi ^{1/2}} \exp {\big ({-}x^{*2}\big )}, \end{equation}

and $F_l(x,t)=\widehat {F_l}(x) e^{i \varOmega t}$ is the 2-D lift force. Here, $\varepsilon$ denotes the smearing parameter, which in numerical simulations should be properly resolved by the chosen grid, and superscript $*$ indicates non-dimensional quantities. This force is written as $F_l = \rho U_{\infty } \varGamma$ , which can be interpreted as a linearised form of the Kutta–Joukowski theorem, in the absence of apparent mass effects.

The solutions of (2.10) and (2.11) become

(2.14) \begin{align} \omega _x &= ( \widehat {\varGamma } * \eta )(x) \, \eta (y) \, \eta (z) \exp {(i \varOmega t})\nonumber \\ &\quad - \frac {i \varOmega }{U_{\infty }} (\widehat {\varGamma } * \eta )(x) \, \eta (y) H_{\varepsilon }\left (\frac {z}{\varepsilon }-i \frac {\varOmega \varepsilon }{2 U_{\infty }}\right ) \exp {\left (i \varOmega \left (t - \frac {z}{U_{\infty }} \right )\right )} \exp {\left ({-}\frac {\varOmega ^{2} \varepsilon ^{2}}{4 U_{\infty }^{2}} \right )} , \end{align}

(2.15) \begin{align} \omega _z &= - \left (\frac {\partial \widehat {\varGamma }}{\partial x} * \eta \right )(x) \, \eta (y) H_{\varepsilon }\left (\frac {z}{\varepsilon }-i \frac {\varOmega \varepsilon }{2 U_{\infty }}\right ) \exp {\left (i \varOmega \left (t - \frac {z}{U_{\infty }} \right )\right )} \exp {\left ({-}\frac {\varOmega ^{2} \varepsilon ^{2}}{4 U_{\infty }^{2}} \right )} , \end{align}

where

(2.16) \begin{equation} H_{\varepsilon }(z)=\frac {\textrm{erf}(z)+1}{2} . \end{equation}

From these equations, it is possible to note that a relevant parameter of the problem is the reduced frequency, $k_{\varepsilon }$ , defined by $\varepsilon$ :

(2.17) \begin{equation} k_{\varepsilon } = \frac {\varOmega \varepsilon }{2 U_{\infty }} , \end{equation}

where factor $2$ is used not only because of its explicit occurrence in the formulae, but also due to the convention in the definition of reduced frequency (Bisplinghoff et al. Reference Bisplinghoff, Ashley and Halfman2013). We can define $\omega _x(x,y,z,t)=\widehat {\omega _x}(x,y,z) e^{i \varOmega t}$ and $\omega _z(x,y,z,t)=\widehat {\omega _z}(x,y,z) e^{i \varOmega t}$ . Non-dimensionalising by $\varepsilon$ and $U_{\infty }$ , the non-dimensional variables are defined as $x^*=x/\varepsilon$ , $t^*=tU_{\infty }/\varepsilon$ , $\omega ^*=\omega \varepsilon /U_{\infty }$ , $\varGamma ^*=\varGamma /(U_{\infty }\varepsilon )$ ,

(2.18) \begin{align} \widehat {\omega _x^*} &= ( \widehat {\varGamma ^*} * \eta )(x^*) \, \eta (y^*) \, \eta (z^*) \nonumber\\&\quad - i 2 k_{\varepsilon } (\widehat {\varGamma ^*} * \eta )(x^*) \, \eta (y^*) H_{\varepsilon ^*}\big (z^*-i k_{\varepsilon } \big ) \exp {\big ({-} i 2 k_{\varepsilon } z^* \big )} \exp {\big ({-} k_{\varepsilon }^2 \big )}, \end{align}

(2.19) \begin{align} \widehat {\omega _z^*} &= - \left (\frac {\partial \widehat {\varGamma ^*}}{\partial x^*} * \eta \right )(x^*) \, \eta (y^*) H_{\varepsilon }\big (z^*-i k_{\varepsilon } \big ) \exp {\big ({-} i 2 k_{\varepsilon } z^*\big )} \exp {\big ({-}k_{\varepsilon }^2 \big )} . \end{align}

2.2.2. Induced velocity in 2-D flow

Now, we restrict ourselves to the case of a 2-D flow. From Helmholtz formula (Helmholtz Reference Helmholtz1858), the velocity induced by vorticity can be written as

(2.20) \begin{equation} \boldsymbol{u}(\boldsymbol{x})= \frac {1}{4 \pi } \int \frac {\boldsymbol{\omega }(\boldsymbol{x'}) \times (\boldsymbol{x}-\boldsymbol{x'})}{|\boldsymbol{x}-\boldsymbol{x'}|^3}\, {\rm d}\boldsymbol{x'} = \frac {1}{4 \pi } \int \frac {\widehat {\boldsymbol{\omega }}(\boldsymbol{x'}) \times (\boldsymbol{x}-\boldsymbol{x'})}{|\boldsymbol{x}-\boldsymbol{x'}|^3}\, {\rm d}\boldsymbol{x'} e^{i \varOmega _f t} , \end{equation}

hence,

(2.21) \begin{equation} \widehat {\boldsymbol{u}}(\boldsymbol{x}) =\frac {1}{4 \pi } \int \frac {\widehat {\boldsymbol{\omega }}(\boldsymbol{x'}) \times (\boldsymbol{x}-\boldsymbol{x'})}{|\boldsymbol{x}-\boldsymbol{x'}|^3} \,{\rm d}\boldsymbol{x'} . \end{equation}

For a 2-D flow, restricted to the $y{-}z$ -plane, the vorticity is

(2.22) \begin{align} \widehat {\omega _x^*} = \widehat {\varGamma ^*} \, \eta (y^*) \, \eta (z^*) - i 2 k_{\varepsilon } \widehat {\varGamma ^*} \, \eta (y^*) H_{\varepsilon }\big (z^*-i k_{\varepsilon } \big ) \exp {\big ({-} i 2 k_{\varepsilon } z^* \big )} \exp {\big ({-} k_{\varepsilon }^2 \big )}. \end{align}

The first term of (2.22) is the bound vortex, which does not induce velocity onto a line located at $(y,z)=(0,0)$ , due to symmetry. The induced velocity at the actuator line is then

(2.23) \begin{align} \widehat {u_y}(0,0,0) & = \frac {1}{4 \pi } \int \limits _{-\infty }^{+\infty } \int \limits _{-\infty }^{+\infty } \int \limits _{-\infty }^{+\infty } \frac {\widehat {\omega _x^*} {z^*}'}{\big ({x^*}'^2 + {y^*}'^2 + {z^*}'^2 \big )^{3/2}} \,{\rm d}{x^*}' \,{\rm d}{y^*}'\, {\rm d}{z^*}' \nonumber\\& = - \frac {i 2 k_{\varepsilon } \widehat {{\varGamma ^*}} \exp {\big ({-} k_{\varepsilon }^2 \big )}}{4 \pi } \int \limits _{-\infty }^{+\infty } \int \limits _{-\infty }^{+\infty } \int \limits _{-\infty }^{+\infty } \eta ( {y^*}') H_{\varepsilon }\big ( {z^*}'-i k_{\varepsilon } \big )\nonumber\\&\quad \times \exp {\big ({-} i 2 k_{\varepsilon } {z^*}' \big )} \frac { {z^*}'}{\big ( {x^*}'^2 + {y^*}'^2 + {z^*}'^2 \big )^{3/2}}\, {\rm d} {x^*}' \,{\rm d} {y^*}'\, {\rm d} {z^*}' \nonumber\\& = \kappa (k_{\varepsilon }) \widehat { {\varGamma ^*}} , \end{align}

where

(2.24) \begin{equation} \kappa (k_{\varepsilon }) = - \frac {i 2k_{\varepsilon } \exp ({-}k_{\varepsilon }^2) }{2 \sqrt {\pi }} \int \limits _{-\infty }^{+\infty } H_{\varepsilon }\big ( {z^*}-i k_{\varepsilon } \big ) \exp ({-}i 2 k_{\varepsilon } {z^*}) \textrm {sgn}({z^*}) \textrm {erfcx}(| {z^*}|) {\rm d} {z^*} , \end{equation}

where $\textrm{erfcx}(z) = \exp (z^2) ( 1 - \textrm {erf}(z) )$ is the scaled complementary error function (Johnson Reference Johnson2017). The complex function $\kappa$ , which is the ratio between the non-dimensional induced velocity and the non-dimensional circulation, is only a function of the reduced frequency $k_{\varepsilon }$ . The values of $\kappa (k_{\varepsilon })$ , shown in figure 1, were calculated by integrating (2.24) in Matlab, using the Faddeeva package (Johnson Reference Johnson2017) for the error functions.

Figure 1.

Function $\kappa$ . Ratio between induced velocity and circulation.

2.2.3. Lift, circulation and angle of attack

We consider a linear relation between the angle of attack and the aerofoil lift coefficient:

(2.25) \begin{equation} C_l = a_0 (\alpha - \alpha _0) , \end{equation}

where $\alpha _0$ is the angle for zero lift and $a_0$ is the value of ${\rm d} C_l/{\rm d}\alpha$ of the aerofoil. In a linear harmonic analysis, $\alpha _0$ only contributes to the steady component of lift. Therefore, the harmonic component can be separated from the steady component:

(2.26) \begin{equation} \widehat {C_l} = a_0 \widehat {\alpha } . \end{equation}

The effective local angle of attack $\alpha$ can be calculated from

(2.27) \begin{equation} \alpha = \alpha _g + \arctan {\left ( \frac {u_{\textit{ry}}}{u_{rz}}\right )}, \end{equation}

or including the pitch rate term (see § 2.3)

(2.28) \begin{equation} \alpha = \alpha _g + \arctan {\left ( \frac {u_{\textit{ry}}}{u_{rz}}\right )} + b\left (\frac {1}{2}-a \right )\dot \alpha _g, \end{equation}

where $u_{\textit{ry}}$ and $u_{rz}$ are the relative velocities. Linearising the second term,

(2.29) \begin{equation} \arctan {\left ( \frac {u_{\textit{ry}}}{u_{rz}} \right )} = \arctan {\left ( \frac {u_y + \dot h}{u_{rz}}\right )} \approx \frac {u_y}{U_{\infty }}+\frac {\dot h}{U_{\infty }}, \end{equation}

and linearising $\alpha$ ,

(2.30) \begin{equation} \alpha = \alpha _{\textit{QS}} + u_y^*. \end{equation}

Assuming small induced velocities,

(2.31) \begin{equation} \widehat {\alpha } = \widehat {\alpha }_{\textit{QS}} + \widehat { {u_y*}} = \widehat {\alpha }_{\textit{QS}} + \kappa (k_{\varepsilon }) \widehat { {\varGamma ^*}} . \end{equation}

Therefore, the lift coefficient can be rewritten as

(2.32) \begin{equation} \widehat {C_l} = a_0 \widehat {\alpha }_{\textit{QS}} + a_0\kappa (k_{\varepsilon }) \widehat { {\varGamma ^*}} = \widehat {C_l}_{\textit{QS}} + a_0\kappa (k_{\varepsilon }) \widehat { {\varGamma ^*}} , \end{equation}

where $\widehat {C_l}_{\textit{QS}} = a_0 \widehat {\alpha }_{\textit{QS}}$ is the quasi-steady component of the lift coefficient. The lift force is calculated from

(2.33) \begin{equation} {F_l} = \frac {1}{2} {\rho } {u_{r}}^{2} {c} C_l , \end{equation}

where the relative velocity $ {u_{r}}=\sqrt { {u_y}^{2}+ {U_{\infty }}^{2}}$ can be linearised to $ {U_{\infty }}$ . After linearisation,

(2.34) \begin{equation} {F_l} = {\rho } {U_{\infty } \varGamma } \approx \frac {1}{2} {\rho } {U_{\infty }}^{2} {c} C_l, \end{equation}

the relation between circulation and lift coefficient in non-dimensional form becomes

(2.35) \begin{equation} {\varGamma ^*} = \frac { {c^*} C_l}{2} , \end{equation}

where $ {c^*}= {c}/ {\varepsilon }$ . Combining (2.32) and (2.35),

(2.36) \begin{equation} \widehat {C_l} = \frac {1}{1-\frac {a_0 {c^*}}{2}\kappa }\widehat {C_l}_{\textit{QS}} , \end{equation}

which implies analogous relations for the lift and circulation. We observe that the unsteady lift in the present formulation is a function of the reduced frequency $k_\varepsilon$ and is thus a function of the smearing parameter $\varepsilon$ of the Gaussian kernel; different choices of kernel thus modify the unsteady lift. The absolute value of $\widehat {C_l}/\widehat {C_l}_{\textit{QS}}$ and its phase difference are shown in figure 3, which employs the reduced frequency based on the semi-chord, $ {k}$ , for direct comparison with Theodorsen’s function.

2.3. Numerical methods

To implement the actuator line method, we employ the Dedalus code (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020), an open-source framework that enables solving differential equations using spectral methods. We use ALM to compute the lift generated by the harmonic transverse oscillation of a thin aerofoil, which may be compared with the present theoretical results and also with the reference theory by Theodorsen. The ALM approach is solved using direct numerical simulation (DNS) of the governing incompressible Navier–Stokes equations, resulting in the following continuity equation:

(2.37) \begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0, \end{equation}

and the momentum equation in primitive variables (velocity-pressure),

(2.38) \begin{equation} \frac {\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u} = - \frac {1}{\rho } \boldsymbol{\nabla }\boldsymbol{p} + \nu {\nabla} ^2 \boldsymbol{u} + \boldsymbol{f}, \end{equation}

where $\rho$ is the fluid density and $\nu$ is the kinematic viscosity.

In this study, we calculate the body forces $\boldsymbol{f}$ using the thin aerofoil theory, which employs the lift coefficient ( $c_l = 2 \pi \alpha$ ), the free stream velocity, the fluid density and the chord length. In the ALM, the angle of attack ( $\alpha$ ) is usually calculated from the local velocity and the geometric angle of attack, $\alpha _g$ , according to the steady-state approximation,

(2.39) \begin{equation} \alpha = \alpha _g + \arctan \left (\frac {u_{\textit{ry}}}{ {u_{rz}}}\right )\!, \end{equation}

where $u_{\textit{ry}}$ and $u_{rz}$ represent the vertical and horizontal components of the relative velocity. In our simulations, we sample the velocities at the centre of the Gaussian distribution. However, to maintain consistency with Theodorsen’s theory, a modified angle of attack

(2.40) \begin{equation} \alpha _m = \alpha _g + \arctan \left (\frac {u_{\textit{ry}}}{ {u_{rz}}}\right ) + b\left (\frac {1}{2} - a\right ) \dot \alpha _g, \end{equation}

with $a = -1/2$ , is employed. The differences between cases with and without the pitch rate term ( $\dot \alpha_g$ ) will be discussed in § 3, where it is demonstrated that this term is necessary to achieve results consistent with (2.3). It should be noted that the theory developed in § 2.2 produces the same ratio of $\widehat {C_l}$ (2.36) for either (2.39) or (2.40), as long as $\widehat {\alpha }_{\textit{QS}}$ is defined consistently. The body forces are obtained by convolving the two-dimensional force vector ( $\boldsymbol{F}_{\text{2-D}}$ ) term (Kleine et al. Reference Kleine, Hanifi and Henningson2023a )

(2.41) \begin{equation} \boldsymbol{f}^i = \frac {1}{\rho } \boldsymbol{F}_{\text{2-D}} \delta (y) \delta ( {z}) , \end{equation}

with a two-dimensional Gaussian kernel,

(2.42) \begin{align}&\qquad\qquad \boldsymbol{f} = \boldsymbol{f}^i * \eta _2 , \end{align}

(2.43) \begin{align}& \eta _{2} ( {y,z}) := \frac {1}{\varepsilon ^2 \pi } \exp \left ({-}\frac { {y^2+z^2}}{\varepsilon ^2} \right )\!, \end{align}

where $y$ and $z$ denote position, $\eta _{2}$ is the Gaussian kernel, and $\varepsilon$ is the smearing parameter, defining the width parameter of Gaussian function. The value of mesh discretisation $\Delta y = \Delta z= 4$ is selected according to Shives & Crawford (Reference Shives and Crawford2013) (see also Forsting & Troldborg Reference Forsting and Troldborg2020 and Kleine et al. Reference Kleine, Hanifi and Henningson2023a ), who concluded that an $\varepsilon /\Delta y = \varepsilon /\Delta z$ ratio of approximately 4 minimises significant errors.

With a fixed value of $\varepsilon = 1$ , the chord length of the thin aerofoil varies based on $\varepsilon /c$ for each case. The Reynolds number, based on $\varepsilon$ for all simulations, is $Re = \varepsilon \rho _{\infty } U_{\infty }/{\mu } = 10^5$ , where $\rho _{\infty } =1$ is the fluid density, $U_{\infty } = 1$ is the free stream velocity and $\mu$ is the dynamic viscosity. A vertical displacement and pitching around a quarter of the chord ( $a = -1/2$ ) are imposed on the thin aerofoil, which oscillates according to the following equations:

(2.44) \begin{equation} h = h_0 \sin {\varOmega t} \quad \mathrm{and}\quad \alpha _g = \alpha _0 \sin {\varOmega t}, \end{equation}

where $h_0 = 0.01 \varepsilon$ and $\alpha _0 = 1^{\circ }$ represent a small amplitude to enable comparison with the linear theory of § 2.2. Simulations were conducted over a range of frequencies and the results are presented in terms of reduced frequency based on the semi-chord ( $ {k}$ ).

Figure 2.

Illustration of the computational mesh resolution for specific parameters: $\varepsilon$ = 1.0 (smearing length) and $c=2.5$ (chord length), leading to a ratio of $\varepsilon /c=0.4$ . (a) Localised view of the mesh, emphasising the constant discretisation and the spatial influence of the Gaussian kernel (represented by the colour scale) in the vicinity of the force application point. (b) Overview of the full computational domain, demonstrating the global mesh setup.

All parameters are non-dimensionalised by the width parameter of the Gaussian function $\varepsilon$ and the free stream velocity $U_{\infty }$ . The domain dimensions, $L_y$ and $L_{ {z}}$ , are chosen to facilitate the calculation of induced velocity for distant wakes at low frequencies. Given the low-frequency values used for thin aerofoil oscillation, sufficient space is required for detached vortices to develop. The aerofoil is located at the origin of the planar domain ( $h,0$ ) for plunging motion and for pitching motion is located at ( $0,0$ ). For the periodic dimensions, a Fourier basis is used on the z-axis and a sine/cosine basis on the y-axis, which implicitly impose free-slip boundary conditions at the top and bottom domain boundaries. The distance between the inlet ( $ {z}_{in}=-102.4$ ) and the outlet ( $z_{out}=409.6$ ) is $L_{ {z}} = 512$ , and the vertical boundary distance is $L_y = 64$ , equidistant from the z-axis. The mesh resolution has a constant discretisation $\Delta z = \Delta y = \varepsilon /4 = 0.25$ . Figure 2(a) illustrates this resolution within a limited domain, facilitating visualisation of the Gaussian function positioned at (0,0) for pitching simulations. Figure 2(b) displays the mesh resolution across the entire domain.

A fringe zone, designed to absorb the wake, is added around $3/4$ of $L_{ {z}}$ with a width of $L_{ {z}}/16$ . At the outlet, the sponge geometry is carefully defined to prevent wake reflections, enforcing the $y$ -component of the velocity to zero and setting the $z$ -component to the initial free stream velocity. We add the fringe as a forcing term $\sigma (\boldsymbol{u -u_0})$ to the Navier–Stoke equations solver, where $\boldsymbol{u_0}$ denotes the free stream velocity and $\sigma$ represents the damping coefficient, which is smoothly varied along z from zero to its maximal value with a tanh function. The accuracy of the fringe is verified by observing that negligible vorticity appears in the left boundary of the domain. For the direct numerical simulation, no subgrid was employed. However, to mitigate numerical instability arising from large wavenumbers produced by the highest frequencies in the simulations, a filter provided by the Dedalus framework was used. This tool acts as a spatial filter, setting to zero all Fourier coefficients of the velocity solutions associated with wavenumbers larger than a chosen threshold. This effectively introduces artificial dissipation for the energy generated by these larger wavenumbers.