3.1 Background and main concept

Brown & Daniels (Reference Brown and Daniels1975) were the first to extend the steady triple-deck theory to the case of a high-frequency ( $\unicode[STIX]{x1D714}$ ), small-amplitude oscillatory pitching flat plate. Unlike the steady case, there is a Stokes layer near the wall that is of order $\sqrt{\unicode[STIX]{x1D708}/\unicode[STIX]{x1D714}}$ where the viscous term is balanced by the time-derivative term in the equations. Brown & Daniels assumed that the Stokes layer and the lower deck have the same thickness, which results in a reduced frequency $k=O(R^{1/4})$ , where $k=\unicode[STIX]{x1D714}b/U$ . This range is too large for engineering applications: $k\simeq 5{-}15$ for $R=10^{4}{-}10^{6}$ . Then, the matching between the adverse pressure gradient due to oscillation and the triple-deck favourable pressure gradient results in a pitching amplitude of the order of $O(R^{-9/16})$ , which is also impractically small for engineering applications: $\simeq 0.02^{\circ }{-}0.32^{\circ }$ for $R=10^{4}{-}10^{6}$ . Indeed, their work is for very-high-frequency, very-small-amplitude oscillations. Note that, over this range, in contrast to the engineering-relevant problem considered in this work, the triple-deck problem becomes mathematically interesting as the time-derivative term appears in the lower-deck equations.

Brown & Cheng (Reference Brown and Cheng1981) extended the work of Brown & Daniels (Reference Brown and Daniels1975) to a more practical range of parameters $0<k\ll R^{1/4}$ . They provided a solution for the case of a flat plate pitching about its mid-chord at $k=1/2$ . In this section, we extend their work to an arbitrarily deforming thin airfoil, arbitrary $k$ in the range $0<k\ll R^{1/4}$ , and correct for a few minor mistakes in their derivation. More importantly, we use the developed theory, within a describing function formulation (Krylov & Bogoliubov Reference Krylov and Bogoliubov1943) assuming weakly nonlinear dynamics, to provide a viscous extension of the classical Theodorsen’s lift frequency response, which was not provided by Brown & Cheng (Reference Brown and Cheng1981).

A key element in the theoretical development below is the vanishing of the time-derivative term from the triple-deck equations over the range $0<k\ll R^{1/4}$ . That is, the lower-deck equations are quite similar to those of the steady case at a non-zero $\unicode[STIX]{x1D6FC}_{s}$ studied by Brown & Stewartson (Reference Brown and Stewartson1970) with a proper definition for the equivalent steady angle of attack. However, we emphasize that this approach is not a quasi-steady solution; although the time-derivative term does not show up in the lower-deck equations, the correspondence with the steady equations implies an equivalent angle of attack that is dependent on the oscillation frequency, as shown below. Therefore, the lower-deck system is dynamical (i.e. possesses a non-trivial frequency response). In fact, even with no time-derivative term in the lower-deck equations, it is not obvious how the steady results of Brown & Stewartson (Reference Brown and Stewartson1970) can be readily applied because the upstream flow is unsteady with Stokes layer in the perturbed Blasius layer. Brown & Daniels (Reference Brown and Daniels1975) encountered a similar problem: how to match the solution of the perturbed Blasius boundary layer (with its inner Stokes layer) with the main deck of the triple-deck structure? They resolved this issue by introducing a transition region, whose length is $O(k^{-1})$ , between the perturbed Blasius boundary layer and the triple deck, called the fore deck. It has similar structure to that upstream of the triple deck: outer potential flow, main boundary layer and an inner Stokes layer, as shown in figure 5. Therefore, in the low-frequency problem where the triple-deck equations are void of the time-derivative term, to match the unsteady flow upstream of the triple deck with the ‘quasi-steady’ flow in the triple deck, Brown & Cheng (Reference Brown and Cheng1981) inserted a second fore deck between the first fore deck and the triple deck, as shown in figure 5. As such, the numerical results of Chow & Melnik (Reference Chow and Melnik1976) to the steady problem of Brown & Stewartson (Reference Brown and Stewartson1970) could be readily used with an equivalent angle of attack. Since the equivalent steady angle of attack $\unicode[STIX]{x1D6FC}_{s}$ is proportional to $A_{\unicode[STIX]{x1D6FC}}k^{2}$ (where $A_{\unicode[STIX]{x1D6FC}}$ is the amplitude of oscillation), and the steady triple-deck formulation of Brown & Stewartson (Reference Brown and Stewartson1970) is valid for $\unicode[STIX]{x1D6FC}_{s}=O(R^{-1/16})$ , the current unsteady formulation is valid for $A_{\unicode[STIX]{x1D6FC}}k^{2}=O(R^{-1/16})$ , which is quite relevant to engineering applications (e.g. $R\simeq 10^{6}$ , $A_{\unicode[STIX]{x1D6FC}}\simeq 5^{\circ }$ and $k<5$ ).

Figure 5. Low-frequency triple-deck structure and flow regimes.

3.2 Theoretical development

3.2.1 Potential-flow set-up

Figure 6. A flexible/deformable thin airfoil defined by the time-varying camber function $y_{c}(x,t)$ .

Consider an arbitrarily deforming thin airfoil (i.e. of time-varying camber) in the presence of a uniform stream $U$ , as shown in figure 6. In classical thin airfoil theory (e.g. Robinson & Laurmann Reference Robinson and Laurmann1956; Schlichting & Truckenbrodt Reference Schlichting and Truckenbrodt1979; Bisplinghoff, Ashley & Halfman Reference Bisplinghoff, Ashley and Halfman1996), it is typical to assume the following series solution for the pressure distribution over the upper surface, which automatically satisfies the Kutta condition (zero loading at the trailing edge):

(3.1) $$\begin{eqnarray}P(\unicode[STIX]{x1D703},t)-P_{\infty }=\unicode[STIX]{x1D70C}\left[\frac{1}{2}a_{0}(t)\tan \frac{\unicode[STIX]{x1D703}}{2}+\mathop{\sum }_{n=1}^{\infty }a_{n}(t)\sin n\unicode[STIX]{x1D703}\right],\end{eqnarray}$$

where $\unicode[STIX]{x1D703}$ is related to $x$ via $x=b\cos \unicode[STIX]{x1D703}$ and $a_{0}$ represents the leading-edge singularity. The pressure on the lower side is given by the negative of (3.1). Moreover, if the plate’s normal velocity $v_{p}$ is written as

(3.2) $$\begin{eqnarray}v_{p}(\unicode[STIX]{x1D703},t)=\frac{1}{2}b_{0}(t)+\mathop{\sum }_{n=1}^{\infty }b_{n}(t)\cos n\unicode[STIX]{x1D703},\end{eqnarray}$$

then the no-penetration boundary condition will provide a means to determine all the coefficients $a_{n}$ (except $a_{0}$ ) in terms of the plate motion kinematics ( $b_{n}$ ) as (Robinson & Laurmann Reference Robinson and Laurmann1956, p. 491)

(3.3) $$\begin{eqnarray}a_{n}(t)=\frac{b}{2n}{\dot{b}}_{n-1}(t)+Ub_{n}(t)-\frac{b}{2n}{\dot{b}}_{n+1}(t),\quad \forall ~n\geqslant 1.\end{eqnarray}$$

The determination of $a_{0}$ is more involved in the sense that it requires solving an integral equation, which cannot be solved analytically for arbitrarily time-varying wing motion. It has been solved for some common inputs; e.g. step change in the angle of attack resulting in Wagner’s response (Wagner Reference Wagner1925), simple harmonic motion resulting in Theodorsen’s frequency response (Theodorsen Reference Theodorsen1935) and a sharp-edged gust (Küssner Reference Küssner1929). Since the focus of this work is to provide a viscous extension of Theodorsen’s frequency response, consider the simple harmonic motion

(3.4a,b ) $$\begin{eqnarray}v_{p}(\unicode[STIX]{x1D703},t)=V_{p}(\unicode[STIX]{x1D703})\text{e}^{\text{i}\unicode[STIX]{x1D714}t},\quad V_{p}(\unicode[STIX]{x1D703})=\frac{1}{2}B_{0}+\mathop{\sum }_{n=1}^{\infty }B_{n}\cos n\unicode[STIX]{x1D703},\end{eqnarray}$$

where the spatially varying amplitude $V_{p}(\unicode[STIX]{x1D703})$ may be complex and $\unicode[STIX]{x1D714}$ is the oscillation frequency. Then, $a_{0}$ is written as (Robinson & Laurmann Reference Robinson and Laurmann1956, p. 496)

(3.5) $$\begin{eqnarray}a_{0}(t)=U[(B_{0}+B_{1})C(k)-B_{1}]\text{e}^{\text{i}\unicode[STIX]{x1D714}t}.\end{eqnarray}$$

(Note that the presentation of Robinson & Laurmann has been adapted to a more common and modern notation.) In equation (3.5), $C(k)$ is Theodorsen’s frequency response function, which depends on the reduced frequency $k=\unicode[STIX]{x1D714}b/U$ via

(3.6) $$\begin{eqnarray}C(k)=\frac{\text{H}_{1}^{(2)}(k)}{\text{H}_{1}^{(2)}(k)+\text{i}\text{H}_{0}^{(2)}(k)},\end{eqnarray}$$

where $\text{H}_{n}^{(m)}$ is the Hankel function of the $m$ th kind of order $n$ . Finally, the potential-flow coefficients of lift and pitching moment (positive pitching up) at the mid-chord point are written as

(3.7a,b ) $$\begin{eqnarray}C_{L_{P}}=-\frac{\unicode[STIX]{x03C0}}{U^{2}}(a_{0}+a_{1})\quad \text{and}\quad C_{M_{0P}}=\frac{\unicode[STIX]{x03C0}}{4U^{2}}(a_{2}-a_{0}).\end{eqnarray}$$

3.2.2 Viscous correction

Following the approach of Brown & Stewartson (Reference Brown and Stewartson1970) in the steady problem (described above), we relax the Kutta condition in the unsteady inviscid pressure distribution (3.1) by introducing a correction $\unicode[STIX]{x1D6E4}_{v}$ to the Kutta circulation. This additional circulation will naturally introduce a singularity at the trailing edge. Similar to the steady case, there is no means within potential flow to determine the strength of such a singularity (additional circulation); the essence behind the Kutta condition is to remove such a singularity (condition of least singularity (Crighton Reference Crighton1985)). This additional circulation will have two effects on the unsteady inviscid pressure distribution (3.1): (i) a steady-like effect with two singularities at the leading and trailing edges, similar to the $B_{s}$ term in (2.2); and (ii) an unsteady effect from the interaction with the wake. The latter has a singularity only at the leading edge. As such, the modified pressure distribution can be written as

(3.8) $$\begin{eqnarray}P(\unicode[STIX]{x1D703},t)-P_{\infty }=\unicode[STIX]{x1D70C}\left[\frac{1}{2}a_{0}(t)\tan \frac{\unicode[STIX]{x1D703}}{2}+\mathop{\sum }_{n=1}^{\infty }a_{n}(t)\sin n\unicode[STIX]{x1D703}+\frac{1}{2}B_{v}(t)\left(\cot \frac{\unicode[STIX]{x1D703}}{2}+a_{0_{v}}(t)\tan \frac{\unicode[STIX]{x1D703}}{2}\right)\right],\end{eqnarray}$$

where the correction $B_{v}$ is related to the additional circulation as $B_{v}=U\unicode[STIX]{x1D6E4}_{v}/(2\unicode[STIX]{x03C0}b)$ and $a_{0_{v}}$ is the total leading-edge singularity effect from the two contributions of $\unicode[STIX]{x1D6E4}_{v}$ , mentioned above. This term has a non-trivial dynamics (there is a non-trivial transfer function from $\unicode[STIX]{x1D6E4}_{v}$ to $a_{0_{v}}$ ). It can be determined from potential-flow considerations: it is the $a_{0}$ term in the unsteady inviscid pressure distribution (3.1) over the plate due to a bound circulation $\unicode[STIX]{x1D6E4}_{v}$ , ignoring the quasi-steady contribution (i.e. the wake effects only). Therefore, similar to the general $a_{0}$ term, it cannot be determined analytically for arbitrary kinematics; there is an analytical expression in the special case of harmonic motion ( $a_{0_{v}}=2C(k)-1$ (Brown & Cheng Reference Brown and Cheng1981)).

To determine the viscous correction $B_{v}$ , without resorting to the Kutta condition ( $B_{v}=0$ ), we use the unsteady triple-deck theory. Approaching the trailing edge ( $\unicode[STIX]{x1D703}\rightarrow 0$ or $\tilde{x}\rightarrow 1$ ), the inviscid pressure (with the $B_{v}$ term) is written as

(3.9) $$\begin{eqnarray}P(\tilde{x}\rightarrow 1;t)-P_{\infty }=\unicode[STIX]{x1D70C}\left[\left(\frac{1}{2}a_{0}(t)+2\mathop{\sum }_{n=1}^{\infty }na_{n}(t)+\frac{1}{2}B_{v}(t)a_{0_{v}}(t)\right)\sqrt{\frac{1-\tilde{x}}{2}}+\frac{B_{v}(t)/2}{\sqrt{{\displaystyle \frac{1-\tilde{x}}{2}}}}\right].\end{eqnarray}$$

(Equation (3.9) reduces to (2.2) in the work of Brown & Stewartson (Reference Brown and Stewartson1970) for the steady case ( $k=0$ ).) In this case, $a_{n}=0$ for all $n\geqslant 1$ , $a_{0}=-2U^{2}\unicode[STIX]{x1D6FC}_{s}$ and $a_{0_{v}}=0$ , which yields (2.3) in this work. It also reduces to (2.7) in the work of Brown & Daniels (Reference Brown and Daniels1975) for their case of a harmonically pitching flat plate about its mid-chord point at very high frequency ( $k\gg 1$ ). In this case, $a_{0}$ and $a_{1}$ (proportional to $\dot{\unicode[STIX]{x1D6FC}}$ ) are neglected with respect to $a_{2}=-b^{2}\ddot{\unicode[STIX]{x1D6FC}}/4$ , and $a_{n}=0$ for all $n>2$ .

The unsteady inviscid pressure (3.9) has the same form as the steady one (2.3) with

(3.10) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FC}_{s}(t)\equiv \frac{1}{U^{2}}\left|\frac{1}{2}a_{0}(t)+2\mathop{\sum }_{n=1}^{\infty }na_{n}(t)+\frac{1}{2}B_{v}(t)a_{0_{v}}(t)\right|, & \displaystyle\end{eqnarray}$$

(3.11) $$\begin{eqnarray}\displaystyle & \displaystyle B_{v}(t)\equiv B_{s}=-2\unicode[STIX]{x1D716}^{3}\unicode[STIX]{x1D706}^{-5/4}\left(\frac{1}{2}a_{0}(t)+2\mathop{\sum }_{n=1}^{\infty }na_{n}(t)+\frac{1}{2}B_{v}(t)a_{0_{v}}(t)\right)B_{e}(\unicode[STIX]{x1D6FC}_{e}). & \displaystyle\end{eqnarray}$$

This comparison, along with the fact that the time-derivative term does not enter the triple-deck equations, points to the possibility of directly using the steady solution by Chow & Melnik (Reference Chow and Melnik1976) of the inner-deck equations for the unsteady case with the equivalence shown above, valid in the range $0<k\ll O(R^{1/4})$ . In the above equivalence, if the term $(1/2)a_{0}+2\sum _{n=1}^{\infty }na_{n}+(1/2)B_{v}a_{0_{v}}$ is negative, then the top of the oscillating thin airfoil will correspond to the top of the steady plate; and if it is positive, then the top of the oscillating thin airfoil should correspond to the bottom of the steady plate. In either case, the equivalent steady angle of attack $\unicode[STIX]{x1D6FC}_{s}$ would be positive. In fact, this correspondence has led to the following interesting behaviour: while there is always a significant lift decrease at the trailing-edge stall angle in the steady case as shown in figure 4, there can be either an increase or a decrease in the unsteady lift when $\unicode[STIX]{x1D6FC}_{s}$ reaches the trailing-edge stall value, as shown by Brown & Cheng (Reference Brown and Cheng1981). The above equivalence was mistakenly performed in Brown & Cheng (Reference Brown and Cheng1981) – see equations (2.2), (2.9) and (2.12) in their work.

The application procedure will be as follows. The airfoil kinematics will be used to determine the instantaneous values of $a_{n}(t)$ via the no-penetration boundary condition (3.3) and (3.5). These coefficients will define the equivalent steady angle of attack $\unicode[STIX]{x1D6FC}_{s}$ according to (3.10), which defines $\unicode[STIX]{x1D6FC}_{e}$ according to (2.4), resulting in $B_{e}$ via the numerical solution of Chow & Melnik (Reference Chow and Melnik1976). Finally, $B_{v}$ will be determined from $B_{e}$ and the $a_{n}$ values according to (3.11), which represents the viscous correction to the Kutta condition and consequently to the lift and moment as

(3.12a,b ) $$\begin{eqnarray}C_{L}=-\frac{\unicode[STIX]{x03C0}}{U^{2}}[a_{0}+a_{1}+B_{v}(1+a_{0_{v}})]\quad \text{and}\quad C_{M_{0}}=\frac{\unicode[STIX]{x03C0}}{4U^{2}}[a_{2}-a_{0}+B_{v}(1-a_{0_{v}})].\end{eqnarray}$$

This procedure is iterative at each time step because the input ( $\unicode[STIX]{x1D6FC}_{s}$ ) depends on the output ( $B_{v}$ ): to determine $\unicode[STIX]{x1D6FC}_{s}$ , one needs $B_{v}$ , which would not be determined until $\unicode[STIX]{x1D6FC}_{s}$ is known. However, our computational results show that the $B_{v}$ contribution to $\unicode[STIX]{x1D6FC}_{s}$ is quite negligible. As such, it is fair to consider the following equivalence instead of (3.10) and (3.11)

(3.13a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}_{s}(t)\equiv \frac{1}{U^{2}}\left|\frac{1}{2}a_{0}(t)+2\mathop{\sum }_{n=1}^{\infty }na_{n}(t)\right|,\quad B_{v}(t)\equiv -2\unicode[STIX]{x1D716}^{3}\unicode[STIX]{x1D706}^{-5/4}\left(\frac{1}{2}a_{0}(t)+2\mathop{\sum }_{n=1}^{\infty }na_{n}(t)\right)B_{e}(\unicode[STIX]{x1D6FC}_{e}),\end{eqnarray}$$

which eliminates the need for iteration at each time step.

It should be noted that this procedure admits arbitrary time variation of the airfoil camber (not necessarily harmonic); only $a_{0}$ should be modified accordingly instead of using (3.5). Nevertheless, because there might not be exact closed-form expressions for $a_{0}(t)$ due to other kinematics (e.g. step input), we recommend using (3.5) to construct a viscous frequency response (describing function (Krylov & Bogoliubov Reference Krylov and Bogoliubov1943)), assuming a weakly nonlinear system, and then using the Fourier transform to obtain the viscous lift force and pitching moment due to an arbitrarily time-varying camber (as shown in Garrick (Reference Garrick1938) and Bisplinghoff et al. (Reference Bisplinghoff, Ashley and Halfman1996, pp. 282–283)). This procedure is demonstrated in more detail in the next section.

4.1 Set-up of the frequency response (describing function)

The above approach can be used to construct a viscous extension of Theodorsen’s function at a given Reynolds number. For practical use, we opt to show such an extension for a pitching–plunging flat plate. In this case, the normal velocity of the plate (assuming small disturbances ${\dot{h}}$ and $\unicode[STIX]{x1D6FC}$ ) is written as

(4.1) $$\begin{eqnarray}v_{p}(x,t)={\dot{h}}(t)-\dot{\unicode[STIX]{x1D6FC}}(t)(x-ab)-U\unicode[STIX]{x1D6FC},\quad -b\leqslant x\leqslant b,\end{eqnarray}$$

where $h$ is the plunging displacement (positive upwards), $\unicode[STIX]{x1D6FC}$ is the pitching angle (angle of attack, positive clockwise) and $ab$ represents the chordwise distance from the mid-chord point to the hinge point, as shown in figure 6. This kinematics results in

(4.2a-c ) $$\begin{eqnarray}b_{0}(t)=2[{\dot{h}}(t)+ab\dot{\unicode[STIX]{x1D6FC}}(t)-U\unicode[STIX]{x1D6FC}]=2v_{1/2}(t),\quad b_{1}(t)=-b\dot{\unicode[STIX]{x1D6FC}}(t)\quad \text{and}\quad b_{n}=0\quad \forall ~n>1,\end{eqnarray}$$

where $v_{1/2}$ is the normal velocity of the mid-chord point. As such, for the harmonic motion

(4.3a,b ) $$\begin{eqnarray}h(t)=Hb\text{e}^{\text{i}\unicode[STIX]{x1D714}t}\quad \text{and}\quad \unicode[STIX]{x1D6FC}(t)=A_{\unicode[STIX]{x1D6FC}}\text{e}^{\text{i}\unicode[STIX]{x1D714}t},\end{eqnarray}$$

equations (3.3) and (3.5) result in the following coefficients:

(4.4a-c ) $$\begin{eqnarray}a_{0}(t)=U[2V_{3/4}C(k)\text{e}^{\text{i}\unicode[STIX]{x1D714}t}+b\dot{\unicode[STIX]{x1D6FC}}(t)],\quad a_{1}(t)=b\dot{v}_{1/2}-bU\dot{\unicode[STIX]{x1D6FC}}(t),\quad a_{2}(t)=-\frac{b^{2}\ddot{\unicode[STIX]{x1D6FC}}(t)}{4},\end{eqnarray}$$

where $v_{3/4}(t)=V_{3/4}\text{e}^{\text{i}\unicode[STIX]{x1D714}t}$ is the normal velocity at the three-quarter-chord point. Also, the coefficient $a_{0_{v}}$ is given by (Brown & Cheng Reference Brown and Cheng1981)

(4.5) $$\begin{eqnarray}a_{0_{v}}=2C(k)-1.\end{eqnarray}$$

Note that in this harmonic formulation, equation (4.4) may yield complex values for the coefficients $a_{n}$ ; the actual coefficients, to be used in the series (3.8), are determined by taking the real parts of those in (4.4). In the common classification proposed by Theodorsen (Reference Theodorsen1935), the coefficient $a_{0}$ represents the circulatory contribution while the other two coefficients ( $a_{1},a_{2}$ ) represent the non-circulatory contribution. The lift coefficient is then written as

(4.6) $$\begin{eqnarray}C_{L}(t)=\underbrace{\underbrace{-\unicode[STIX]{x03C0}\frac{b\dot{v}_{1/2}(t)}{U^{2}}}_{non\text{-}circulatory}+\underbrace{2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FC}_{3/4}(t)C(k)}_{circulatory}}_{potential\;flow\;solution}-\underbrace{2\unicode[STIX]{x03C0}\tilde{B}_{v}(t)C(k)}_{viscous\;correction},\end{eqnarray}$$

where $\tilde{B}_{v}=B_{v}/U^{2}$ , $\unicode[STIX]{x1D6FC}_{3/4}$ is the local angle of attack at the three-quarter-chord point (as recommended by the Pistolesi theorem (Schlichting & Truckenbrodt Reference Schlichting and Truckenbrodt1979, p. 80)), and the multiplication $\unicode[STIX]{x1D6FC}_{3/4}(t)C(k)$ is interpreted after writing $\unicode[STIX]{x1D6FC}_{3/4}(t)=\overline{\unicode[STIX]{x1D6FC}_{3/4}}\text{e}^{\text{i}\unicode[STIX]{x1D714}t}$ as

(4.7) $$\begin{eqnarray}(\unicode[STIX]{x1D6FC}_{3/4}C(k))(t)=\text{Re}(\overline{\unicode[STIX]{x1D6FC}_{3/4}}C(k)\text{e}^{\text{i}\unicode[STIX]{x1D714}t}),\end{eqnarray}$$

where $\overline{\unicode[STIX]{x1D6FC}_{3/4}}$ may be a complex number and $\text{Re}(\cdot )$ denotes the real part of its complex argument.

Recall that if $u(t)=A\text{e}^{\text{i}\unicode[STIX]{x1D714}t}$ is the input to a linear dynamical system whose frequency response is $G(\text{i}\unicode[STIX]{x1D714})$ , then the steady-state output is simply written as $y(t)=A|G(\text{i}\unicode[STIX]{x1D714})|\text{e}^{\text{i}\unicode[STIX]{x1D714}t+\text{arg}\,G(\text{i}\unicode[STIX]{x1D714})}$ (Ogata & Yang Reference Ogata and Yang1970). The describing function technique represents an extension of the frequency response concept for weakly nonlinear systems (Krylov & Bogoliubov Reference Krylov and Bogoliubov1943). In this technique, only the response at the fundamental frequency is considered and the higher harmonics are neglected. As such, the response of a weakly nonlinear system to the input $u(t)=A\text{e}^{\text{i}\unicode[STIX]{x1D714}t}$ is approximated as $y(t)=Y(A,\unicode[STIX]{x1D714})\text{e}^{\text{i}\unicode[STIX]{x1D714}t+\unicode[STIX]{x1D719}(A,\unicode[STIX]{x1D714})}$ . That is, unlike linear systems, the magnitude and phase of the transfer function depend on the input amplitude. Using such a technique, we provide below a viscous extension of Theodorsen’s frequency response, i.e. the frequency response between the quasi-steady lift (input) and the viscous circulatory lift (output).

The system possesses two nonlinearities as shown in figure 8: a multiplicative nonlinearity and the triple-deck viscous nonlinearity. The former is due to interactions between the airfoil motion (represented by the $a_{n}$ ) and the trailing-edge singular behaviour (represented by $B_{s}$ ) while the latter is due to the steady triple-deck nonlinear characteristics, determined numerically by Chow & Melnik (Reference Chow and Melnik1976), and shown here in figure 3(b). The effect of these nonlinearities is minimal with respect to the main linear contribution at small angles of attack as evident from the power spectra (fast Fourier transform, FFT) of the total circulatory lift coefficient $C_{L_{C}}$ shown in figure 7(a) for the case of pitching around the mid-chord point with $A_{\unicode[STIX]{x1D6FC}}=1^{\circ }$ at $k=0.8$ and $R=10^{5}$ . The FFT of $C_{L_{C}}$ has a single distinct peak at the operating frequency ( $k=0.8$ ). In fact, even the viscous contribution (the $B_{v}$ term) is mostly linear despite the existence of a weak cubic nonlinearity as evident from its FFT shown in figure 7(b). This weakly nonlinear behaviour of the system justifies the use of the describing function approach (akin to linearization) to construct a frequency response. It is interesting to note that the triple-deck theory confirms the common expectation that the most significant term in the power series expansion of lift in terms of the angle of attack after the linear term is the cubic one (Librescu, Chiocchia & Marzocca Reference Librescu, Chiocchia and Marzocca2003; Ding & Wang Reference Ding and Wang2006). It is also interesting to point out that even at this small amplitude ( $A_{\unicode[STIX]{x1D6FC}}=1^{\circ }$ ) and relatively large Reynolds number ( $R=10^{5}$ ), the viscous contribution is approximately 19 %, as shown in figure 7(a).

Figure 7. Power spectra (FFT) of (a) the total circulatory lift coefficient $C_{L_{C}}$ and (b) the viscous correction $2\unicode[STIX]{x03C0}\tilde{B}_{v}C(k)$ , both normalized by $2\unicode[STIX]{x03C0}A_{\unicode[STIX]{x1D6FC}}$ , for the case of a flat plate pitching around the mid-chord point with $A_{\unicode[STIX]{x1D6FC}}=1^{\circ }$ at $k=0.8$ and $R=10^{5}$ . The behaviour is almost linear with a single distinct peak at this small amplitude.

Figure 8. A block diagram showing the different components constituting the dynamics of the viscous circulatory lift. The airfoil motion dictates the angle of attack $\unicode[STIX]{x1D6FC}_{3/4}$ at the three-quarter-chord point, which is the main input to Theodorsen’s inviscid linear dynamics, resulting in the potential-flow circulatory lift. The upper branch represents the viscous correction developed in this work. The correction term $\tilde{B}_{v}$ represents a singularity in the inviscid pressure distribution at the trailing edge. It should be set to zero according to the Kutta condition. Rather, it is obtained here from the triple-deck boundary layer theory. The airfoil motion goes into some linear dynamics (that includes the Theodorsen function) to obtain an equivalent steady angle of attack, which will be used in the nonlinear triple-deck theory to obtain the viscous correction to the circulatory lift.

4.2 Computation procedure

Let $k$ and $R$ be given. Then, the quasi-steady lift coefficient (input to the sought lift dynamical system) is written as

(4.8) $$\begin{eqnarray}C_{L_{QS}}(t)=2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FC}_{3/4}(t).\end{eqnarray}$$

Also, the coefficients $a_{0}$ , $a_{1}$ and $a_{2}$ are given from (4.4). Thus, $\unicode[STIX]{x1D6FC}_{s}$ can be obtained accordingly from (3.13); however, care should be taken when applying (3.13). It should be applied instantaneously as

(4.9) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}_{s}(t)=\frac{1}{U^{2}}\left|\frac{1}{2}a_{0}(t)+2a_{1}(t)+4a_{2}(t)\right|.\end{eqnarray}$$

The scaled angle of attack $\unicode[STIX]{x1D6FC}_{e}(t)$ for the numerical solution of Chow & Melnik (Reference Chow and Melnik1976) is obtained from (2.4) with $\unicode[STIX]{x1D716}=R^{-1/8}$ . Using figure 3(b), one can obtain $B_{e}(t)$ , which in turn is substituted in (3.13) to determine the viscous correction $B_{v}(t)$ . The unsteady viscous circulatory lift coefficient is determined from (4.6) by excluding the first term, i.e.

(4.10) $$\begin{eqnarray}C_{L_{C}}(t)=2\unicode[STIX]{x03C0}\,\text{Re}[(\unicode[STIX]{x1D6FC}_{3/4}(t)-\tilde{B}_{v}(t))C(k)].\end{eqnarray}$$

Finally, a spectral analysis (e.g. FFT) is applied to $C_{L_{C}}(t)$ to extract its relative amplitude and phase shift with respect to $C_{L_{QS}}(t)$ . That is, the circulatory lift viscous transfer function $C_{v}$ is defined as

(4.11) $$\begin{eqnarray}C_{v}(k;R)\triangleq \frac{C_{L_{C}}(k;R)}{C_{L_{QS}}(k)}.\end{eqnarray}$$

Figure 8 shows a block diagram for the dynamics of the unsteady viscous circulatory lift.

Note that if $\unicode[STIX]{x1D6FC}_{e}(t)$ exceeds 0.47, then the simulation should be terminated because such a value implies trailing-edge stall beyond which the current analysis is not valid. This limitation defines the region of applicability of the developed model. For example, the model can handle an oscillation about the mid-chord point with $3^{\circ }$ amplitude at $0.4$ reduced frequency and $10\,000$ Reynolds number. However, it cannot handle the same situation when the amplitude is increased to $4^{\circ }$ as $\unicode[STIX]{x1D6FC}_{e}(t)$ would exceed 0.47 during the course of the simulation.

Following the above procedure, we construct frequency responses of the unsteady, viscous, circulatory lift coefficient $C_{L_{C}}$ at different Reynolds numbers, which are shown in figure 9 in comparison to Theodorsen’s. Intuitively, as $R$ increases, the viscous response approaches the inviscid Theodorsen’s response and vice versa. The viscous lift transfer function has a smaller magnitude than the inviscid Theodorsen’s solution as $R$ decreases. Moreover, it is found that viscosity induces a significant phase lag beyond Theodorsen’s; the larger the oscillation frequency and the smaller the Reynolds number, the larger is the discrepancy between Theodorsen’s phase lag and the viscous results. Interestingly, this finding supports the conclusions of some of the earlier experimental efforts (Chu & Abramson Reference Chu and Abramson1959; Bass et al. Reference Bass, Johnson and Unruh1982): Chu & Abramson (Reference Chu and Abramson1959) suggested adding a phase lag of $-10^{\circ }$ to the Theodorsen function for a better estimate of the unsteady lift and flutter boundary when $k\simeq 0.5$ . Bass et al. (Reference Bass, Johnson and Unruh1982) conducted a water tunnel experiment for a NACA 16-012 undergoing pitching oscillations around its quarter-chord point in the range of $0.5<k<10$ and $R=6500{-}26\,500$ . They compared their force measurements to Theodorsen’s potential-flow frequency response. They found bad agreement in the range $0.5<k<2$ where the most pronounced boundary layer activity is observed and the flow near the trailing edge being separated and alternating around the trailing edge. They concluded that adding a phase lag of $-30^{\circ }$ to the Theodorsen’s $C(k)$ would make the predicted lift from the classical theory of unsteady aerodynamics match their experimental measurements over this range, which supports the current results shown in figure 9.

Figure 9. Comparison between the frequency responses of the unsteady, viscous, circulatory lift coefficient $C_{L_{C}}$ at different Reynolds numbers and that of the potential-flow circulatory lift coefficient (i.e. Theodorsen’s): (a) magnitude variation; (b) phase variation. The larger the frequency and the lower the Reynolds number, the larger is the discrepancy in phase between Theodorsen’s function and the viscous frequency response.

This finding is particularly important for the determination of the flutter boundary (e.g. Alben Reference Alben2008; Mandre & Mahadevan Reference Mandre and Mahadevan2010; Zakaria, Al-Haik & Hajj Reference Zakaria, Al-Haik and Hajj2015; Hussein & Canfield Reference Hussein and Canfield2017). Note that the flutter instability, similar to any typical Hopf bifurcation, is mainly dictated by when energy is added/subtracted during the cycle. That is, the phase difference between the applied loads (aerodynamic loads) and the system motion (e.g. angle of attack) plays a crucial role in dictating the stability boundary (Bisplinghoff et al. Reference Bisplinghoff, Ashley and Halfman1996, p. 280). Therefore, if the Theodorsen function does not capture such a phase difference correctly, it will typically lead to an erroneous flutter stability boundary. Hence, if the current model better captures the phase lag, it may enhance our flutter predictability, if occurring at high reduced frequencies. Based on this discussion, we suggest using the obtained viscous frequency responses in place of Theodorsen’s for a more accurate, yet efficient, estimate of the flutter boundary. This point will be discussed further below from an added-mass point of view in the last section.

5.1 Computational set-up

The O-type far field located $25c$ away from the solid body has been implemented for grid generation around the standard NACA 0012 airfoil with sharp trailing edge. In return of closing the blunt trailing edge of the original NACA 0012, the thickness of the airfoil altered to $11.9\,\%$ . To construct the dynamic mesh due to the airfoil motion, the computational domain is divided into three rings as shown in figure 10. The inner ring (red), which encloses the airfoil, has the radius of $6c$ . In this region, a hybrid mesh is used such that a boundary layer structure dense mesh near the airfoil guarantees $y+<1$ , in conjunction with an unstructured tri-mesh attached to it. The distance of the first layer of the mesh was set to be $10^{-5}c$ with $1.1$ growth factor, which guarantees that the triple-deck structure, $O(R^{-3/4})$ , is well resolved at the highest Reynolds number ( $10^{6}$ ) in this simulation. A total of $300$ mesh points were used on each side of the airfoil. A size function has been used to ensure that the unstructured mesh in the inner ring is dense enough to capture the shed vortices if needed. The whole inner ring including the airfoil undergos a rigid-body pitching motion. No dynamic mesh is used in this region to ensure that the grids near the airfoil maintain their fine configuration and quality as they were before the motion.

Figure 10. O-type mesh around the airfoil: the outer ring is fixed, the inner ring moves rigidly with the airfoil, and the intermediate ring represents the deforming dynamic mesh.

The outer ring located at $25c$ away from the airfoil is stationary as if no motion is taking place inside the domain. This fixed mesh near the outer boundaries certifies that the far-field boundary conditions are applied correctly. The intermediate ring plays the main role of the dynamic mesh. The inner radius of this ring is $6c$ and the outer radius is $18c$ ; it occupies a large region inside the domain. Both remeshing and deforming techniques are utilized to damp the deformations in the region caused by the motion of the inner ring. A user-defined function is attached to the solver to impose an arbitrary motion to the airfoil and prevent high skewness in the dynamic mesh zone. The advantage of this method may not be sensible when deflections are small, yet it demonstrates its ability in damping mesh deformations at large motion amplitude. The large size of the intermediate region gives enough room to handle excessive deflections. It has to be pointed out that the position of each ring has been chosen based on numerous simulations. The total number of grids is roughly $2\times 10^{5}$ . For more information about the computational set-up, the reader is referred to our earlier effort (Rezaei & Taha Reference Rezaei and Taha2017). A mesh independence study has been performed by running another case in which the grids were twice as dense, and no alteration in the results has been observed. Thus, the aforementioned mesh configuration was utilized in all of the forthcoming simulations.

Since the flow is assumed to be incompressible, the velocity inlet and pressure outlet, corresponding to the left semicircle and right semicircle, respectively, were set as the far-field boundary conditions. The chord length of the airfoil is $18~\text{cm}$ and the magnitude of the velocity at the inlet boundary is $10~\text{m}~\text{s}^{-1}$ and $1~\text{m}~\text{s}^{-1}$ , corresponding to $R=10^{5}$ and $R=10^{4}$ , respectively. The turbulent intensity of the flow was set to $0.1\,\%$ and the gauge pressure at the outlet boundary condition was set to zero. Note that the free-stream turbulence intensity is not expected to significantly affect the results because the $k{-}\unicode[STIX]{x1D714}$ SST turbulence model utilizes the $k{-}\unicode[STIX]{x1D716}$ model in the free-stream region, which is robust to changes in the inlet conditions. Our simulations with changing the inlet turbulence intensity 10-fold did not show an appreciable change in the lift dynamics. The no-slip boundary condition is imposed on the airfoil which is undergoing the harmonic pitching motion $\unicode[STIX]{x1D6FC}=A_{\unicode[STIX]{x1D6FC}}\sin \unicode[STIX]{x1D714}t$ with a pitching amplitude of $A_{\unicode[STIX]{x1D6FC}}=3^{\circ }$ to ensure that the airfoil is in the pre-stall regime (Schlichting & Truckenbrodt Reference Schlichting and Truckenbrodt1979).

Figure 11. Comparison of the pressure distribution over the flat plate from the inviscid theory, the current viscous theory and computational simulations in the case of a harmonically pitching airfoil about its quarter-chord with $3^{\circ }$ amplitude at $k=1$ and (a)  $R=10^{5}$ and (b)  $R=10^{4}$ . In the former case, a URANS solver is used whereas a laminar solver is used in the latter case. The panels show the pressure distributions at the instant of zero pitching angle $\unicode[STIX]{x1D6FC}(t)=0$ and maximum upward pitching velocity $\dot{\unicode[STIX]{x1D6FC}}$ .

Figure 11 shows the distribution of the pressure difference (i.e. lift distribution) over the flat plate for the case of pitching about the quarter-chord with $3^{\circ }$ amplitude at $k=1$ for two Reynolds numbers, $R=10^{5}$ and $R=10^{4}$ . Figure 11 shows results from the inviscid theory (i.e. twice the result of (3.1)), which is insensitive with respect to $R$ ; the developed viscous theory (i.e. twice the result of (3.8)); and the computational simulations described above. URANS is used in the case of the relatively high $R=10^{5}$ and the laminar solver is used at the lower Reynolds number $R=10^{4}$ . Very good matching between the computational simulations and the developed theory is found. It should be noted that the selected instant ( $\unicode[STIX]{x1D6FC}(t)=0$ and $\dot{\unicode[STIX]{x1D6FC}}$ maximum) is a critical instant during the cycle where the flat plate is on the verge of trailing-edge stall, as shown in figure 12: when $\unicode[STIX]{x1D6FC}(t)=0$ and $\dot{\unicode[STIX]{x1D6FC}}$ is maximum, at this relatively high $k$ , the equivalent $\unicode[STIX]{x1D6FC}_{e}$ becomes very close to the trailing-edge stall value 0.47. Both results from computational simulations and the current boundary layer theory indicate that, as $R$ decreases, the pressure distribution at this critical moment deviates more from the inviscid one; the pressure distribution decreases and shifts to the left (i.e. the pressure attains its maximum earlier on the airfoil). The developed theory captures this behaviour by adding a singularity at the trailing edge, which bends the pressure distribution curve downwards and to the left. However, it is precisely this trailing-edge singularity that causes the discrepancy between the predicted singular pressure and the actual non-singular pressure at the trailing edge. Nevertheless, the strength of this singularity decreases as $R$ increases and the discrepancy becomes more confined to the immediate vicinity of the trailing edge (i.e. agreement with the computational simulations over a wider range of the airfoil). We emphasize that the developed theory should be valid only for high Reynolds numbers. In other words, the results of the developed theory should be interpreted only as providing a first-order correction to the inviscid results for finite Reynolds numbers; from this point of view, it is quite satisfactory. One more point that is noteworthy, though it does not seem to be significant here, is that the developed theory is for an infinitely thin airfoil whereas the computational simulations are performed for a finite-thickness airfoil (NACA 0012).

Figure 12. Variation of the equivalent angle of attack $\unicode[STIX]{x1D6FC}_{e}$ over the cycle along with the actual angle of attack $\unicode[STIX]{x1D6FC}$ for the case of a harmonically pitching flat plate about its quarter-chord with $3^{\circ }$ amplitude at $k=1$ and $R=10^{4}$ . At this relatively high $k$ , the instant of maximum $\dot{\unicode[STIX]{x1D6FC}}$ renders the airfoil on the verge of trailing-edge stall in comparison to the instant of the maximum $\unicode[STIX]{x1D6FC}$ .

5.2 Computation of the viscous lift frequency response function

In this study, we show how the viscous lift frequency response (describing function) is determined from computational simulations at a given Reynolds number. Similar to Theodorsen (Reference Theodorsen1935), this transfer function is defined as the ratio between the circulatory lift coefficient and the quasi-steady lift coefficient. As such, given a combination of $k$ and $R$ , our computational set-up is simulated to result in a time history of the total lift coefficient $C_{L_{tot}}$ . According to the describing function approach (Krylov & Bogoliubov Reference Krylov and Bogoliubov1943), the Fourier transform ${\hat{C}}_{L_{tot}}$ of the total lift coefficient at the fundamental frequency $k$ is considered; it is a complex number, as shown in figure 13. To extract the circulatory contribution ${\hat{C}}_{L_{C}}$ from ${\hat{C}}_{L_{tot}}$ , the non-circulatory contribution must be subtracted. Adopting Theodorsen’s estimate for the non-circulatory loads in the case of a pitching airfoil around the quarter-chord point, we obtain

(5.1) $$\begin{eqnarray}C_{L_{NC}}=\unicode[STIX]{x03C0}\frac{b}{U}\left(\dot{\unicode[STIX]{x1D6FC}}+\frac{b\ddot{\unicode[STIX]{x1D6FC}}}{2U}\right)\quad \Longrightarrow \quad {\hat{C}}_{L_{NC}}(k)=\unicode[STIX]{x03C0}A_{\unicode[STIX]{x1D6FC}}(\text{i}k-k^{2}/2).\end{eqnarray}$$

(Theodorsen’s estimate for the non-circulatory loads may not be accurate, as will be discussed below.) As such, the circulatory lift frequency response is obtained as ${\hat{C}}_{L_{C}}={\hat{C}}_{L_{tot}}-{\hat{C}}_{L_{NC}}$ , as shown graphically in figure 13. Then, the complex number ${\hat{C}}_{L_{C}}$ is divided by the quasi-steady lift

(5.2) $$\begin{eqnarray}C_{L_{QS}}=2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FC}_{3/4}=2\unicode[STIX]{x03C0}\left(\unicode[STIX]{x1D6FC}+\frac{\dot{\unicode[STIX]{x1D6FC}}b}{U}\right)\quad \Longrightarrow \quad {\hat{C}}_{L_{QS}}=2\unicode[STIX]{x03C0}A_{\unicode[STIX]{x1D6FC}}\left(1+\text{i}k\right)\end{eqnarray}$$

to obtain the viscous lift transfer function $C_{v}(k;R)={\hat{C}}_{L_{C}}(k;R)/{\hat{C}}_{L_{QS}}(k)$ .

Figure 13. Complex plane showing the different lift components.

Figure 14 shows the computed viscous lift frequency response function at two different Reynolds numbers, $R=10^{5}$ and $R=10^{4}$ ; the former is obtained using URANS and the latter is obtained using the laminar solver. Figure 14 also shows the inviscid Theodorsen’s lift frequency response function and the developed viscous extension for comparison. For the case of $R=10^{4}$ , the convergence properties of the laminar solver at lower $k$ values were not satisfactory and are therefore omitted. It is found that the magnitude of the transfer function decreases, as $R$ decreases and $k$ increases, more than what is predicted by the triple-deck theory (conforming with the results of Zakaria, Taha & Hajj (Reference Zakaria, Taha and Hajj2017)). However, the computational phase results corroborate the theoretical findings discussed above. That is, at lower Reynolds numbers and higher frequencies, there is a significant deviation from Theodorsen’s phase prediction. In fact, there is a satisfactory quantitative agreement between the theoretical phase lag predictions and computational results. This additional phase lag may significantly affect the prediction of an instability boundary (e.g. flutter) as discussed above. Note that Bisplinghoff et al. (Reference Bisplinghoff, Ashley and Halfman1996, p. 280) emphasized the importance of unsteady phase lag in dictating the flutter boundary even at low reduced frequencies (e.g. $k=0.1$ ) where the phase lag is already very small. Therefore, since the developed theory provides a better estimate of the unsteady phase lag than the Theodorsen function, particularly at large $k$ and low $R$ , it is expected to enhance our retarding capability in predicting flutter, which is discussed further below from an added-mass point of view. It should be noted that most of the earlier experimental efforts that reported failure in predicting flutter (or unsteady loads) lie in the high-frequency range: $k\simeq 0.5$ (Chu & Abramson Reference Chu and Abramson1959), $k\simeq 0.6{-}1.4$ (Henry Reference Henry1961), $k\simeq 0.7$ (Abramson & Ransleben Reference Abramson and Ransleben1965) and $k\simeq 0.5{-}10.0$ (Bass et al. Reference Bass, Johnson and Unruh1982). Therefore, it is expected that the developed theory may help to reconcile the concerns raised in these efforts; a quantitative assessment of the effect of the predicted additional phase lag on the flutter boundary will be addressed in future work. Having said that, one should emphasize that the flutter frequency of conventional airplane wings is usually of the order of $k=0.1$ (Bisplinghoff et al. Reference Bisplinghoff, Ashley and Halfman1996, p. 280), for which the current theory results in a phase close to Theodorsen’s. However, it is expected that the developed theory will be useful for the flutter prediction of the next-generation unconventional designs with highly flexible wings (typically with higher flutter frequencies).

Figure 14. Computational results of the frequency responses of the unsteady, viscous, circulatory lift coefficient $C_{L_{C}}$ at different Reynolds numbers. The computational results support the theoretical finding that the larger the frequency and the lower the Reynolds number, the larger is the discrepancy in phase between Theodorsen function and the viscous frequency response.

It is noteworthy to mention that a better matching between the theoretical phase lag predictions and computational results is obtained when using the eddy viscosity in the developed boundary layer theory, as shown in figure 15 for the case of $R=10^{6}$ . In this case, a Reynolds number based on the average value of the eddy viscosity (obtained from the URANS simulations over a domain close enough and surrounding the airfoil) is used in the developed theoretical model. That is, even when operating at a high Reynolds number where the deviations between Theodorsen results and the developed theory and computational simulations are minimal, the effective Reynolds number is actually less, implying that the phase predictions of the Theodorsen function are quite off at high frequencies. Note that, as $R$ decreases, the turbulent viscosity ratio decreases and its effect on the developed theory may be neglected.

Figure 15. Phase of the lift frequency response at $R=10^{6}$ when using molecular and eddy viscosity. Using eddy viscosity enhances the matching between the theoretical phase lag predictions and computational simulations. A turbulent viscosity ratio of 10 is used based on URANS simulations.

6.1 Viscosity-induced lag

Figure 16. Stokes second problem: flow above an oscillating infinite plate.

The fact that viscosity induces phase lag in the flow response is well known from classical fluid problems. For example, the laminar viscous flow in a pipe due to an oscillatory pressure gradient shows a phase lag between the input pressure gradient and the flow response (e.g. velocity distribution, wall shear or vorticity (Langlois & Deville Reference Langlois and Deville2014, pp. 113–116)). Also, recall the Stokes second problem: the flow above an oscillating infinite plate, shown in figure 16. This problem is one of the few simple problems where an analytical solution of the Navier–Stokes equations is available (Lamb Reference Lamb1932, pp. 619–623; Batchelor Reference Batchelor2000, pp. 191–193; Langlois & Deville Reference Langlois and Deville2014, pp. 109–111), which results in the following velocity distribution:

(6.1) $$\begin{eqnarray}u(y,t)=U\text{e}^{-y/\unicode[STIX]{x1D6FF}}\cos (\unicode[STIX]{x1D714}t-y/\unicode[STIX]{x1D6FF}),\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}=\sqrt{2\unicode[STIX]{x1D708}/\unicode[STIX]{x1D714}}$ is the thickness of the boundary layer (Stokes layer) and $\unicode[STIX]{x1D708}$ is the fluid’s kinematic viscosity. Equation (6.1) clearly shows the phase lag between the input (plate motion) to the flow dynamics and the flow response and that this phase lag increases with the fluid viscosity. Moreover, the vorticity in the boundary layer experiences even more lag in development with respect to the plate motion than the fluid velocity as shown in the vorticity response:

(6.2) $$\begin{eqnarray}\unicode[STIX]{x1D701}(y,t)=-\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}y}=\frac{\sqrt{2}U}{\unicode[STIX]{x1D6FF}}\text{e}^{-y/\unicode[STIX]{x1D6FF}}\cos (\unicode[STIX]{x1D714}t-y/\unicode[STIX]{x1D6FF}-\unicode[STIX]{x03C0}/4).\end{eqnarray}$$

The generation of vorticity in the boundary layer is particularly important for the explanation of the observed lag in the lift frequency response of an oscillating airfoil, for the lift evolution is intimately related to vorticity generation and circulation development.

6.2 Viscous damping and lag in circulation development

To show that the observed phase lag in the lift frequency response is due to lag in the development of the bound circulation, we computed the latter by performing a line integral of the tangential velocity along a closed contour around the airfoil. Then, we followed a similar procedure to the one presented in the last section to construct a frequency response (describing function) between the quasi-steady circulation $\unicode[STIX]{x1D6E4}_{QS}$ as an input and the viscous unsteady bound circulation $\unicode[STIX]{x1D6E4}$ as an output. The former is determined as

(6.3) $$\begin{eqnarray}\unicode[STIX]{x1D6E4}_{QS}=UbC_{L_{QS}}=2\unicode[STIX]{x03C0}bU\left(\unicode[STIX]{x1D6FC}+\frac{\dot{\unicode[STIX]{x1D6FC}}b}{U}\right)\quad \Longrightarrow \quad \hat{\unicode[STIX]{x1D6E4}}_{QS}=2\unicode[STIX]{x03C0}UbA_{\unicode[STIX]{x1D6FC}}(1+\text{i}k).\end{eqnarray}$$

Figure 17 shows a comparison between the viscous transfer function $\hat{\unicode[STIX]{x1D6E4}}/\hat{\unicode[STIX]{x1D6E4}}_{QS}$ of the circulation response from computational simulations and the corresponding potential-flow one, which is different from the Theodorsen function $C(k)$ . Rather, it is given as (Bisplinghoff et al. Reference Bisplinghoff, Ashley and Halfman1996, pp. 275–276)

(6.4) $$\begin{eqnarray}\frac{\hat{\unicode[STIX]{x1D6E4}}_{P}}{\hat{\unicode[STIX]{x1D6E4}}_{QS}}(k)=\frac{-2\text{e}^{-\text{i}k}}{\text{i}k\unicode[STIX]{x03C0}(\text{H}_{1}^{(2)}(k)+\text{i}\text{H}_{0}^{(2)}(k))},\end{eqnarray}$$

where $\unicode[STIX]{x1D6E4}_{P}$ denotes the potential-flow unsteady bound circulation. Similar to the lift transfer function $C_{v}(k;R)={\hat{C}}_{L_{C}}(k;R)/{\hat{C}}_{L_{QS}}(k)$ , the circulation transfer function $\hat{\unicode[STIX]{x1D6E4}}/\hat{\unicode[STIX]{x1D6E4}}_{QS}$ experiences more phase lag due to viscosity at high frequencies than its potential-flow counterpart.

Figure 17. Computational results of the frequency responses of the unsteady, viscous, bound circulation at different Reynolds numbers: (a) magnitude variation and (b) phase variation. A behaviour similar to the circulatory lift frequency response is observed: the larger the frequency, the larger is the discrepancy in phase between inviscid and viscous responses.

Figure 18 shows the vorticity contours during the cycle of a harmonically pitching NACA 0012 airfoil about its quarter-chord with $3^{\circ }$ amplitude at $k=1$ and two different Reynolds numbers, $R=10^{5}$ and $R=10^{4}$ . First, figure 18(a) shows that for this relatively high $k$ , the wake is more deformed at the instants of zero pitching angle ( $\unicode[STIX]{x1D6FC}=0$ ) and maximum angular velocity $\dot{\unicode[STIX]{x1D6FC}}$ in comparison to the instants of maximum $\unicode[STIX]{x1D6FC}$ with $\dot{\unicode[STIX]{x1D6FC}}=0$ . This fact, similar to the theoretical findings above, implies that the airfoil would be more prone to trailing-edge stall at the instants of maximum $\dot{\unicode[STIX]{x1D6FC}}$ when oscillating at high frequencies. Second, the comparison between the two sets of figures at the two values of $R$ indicates a much smaller wake activity (deformation) for lower $R$ , which is intuitively expected due to viscous damping. From a dynamical system perspective, this damping of wake vorticity will be typically associated with lag in its development which, via conservation of circulation, points to a lag in the bound circulation development. Therefore, it may explain the larger phase lag of the lift frequency response found at lower Reynolds numbers and higher frequencies.

Figure 18. Vorticity contours $(1~\text{s}^{-1})$ during the cycle of a harmonically pitching NACA 0012 airfoil about its quarter-chord with $3^{\circ }$ amplitude at $k=1$ and (a) $R=10^{5}$ and (b)  $R=10^{4}$ . As $R$ decreases, viscosity damps the deformation of wake vorticity.

6.3 Lag in circulation development and the Kutta condition

One may be able to relate the observed lag in the circulatory lift frequency response (due to lag in the circulation development) to the Kutta condition. Note that the trailing-edge singularity term in the pressure distribution (3.8) is the main modification introduced to the inviscid pressure distribution (3.1). Therefore, the observed additional lag in the lift response may be related to the pressure near the trailing edge, which motivates the following analysis.

Figure 19. A zoom at the trailing edge and its boundary layer. The blue lines represent the edge of the boundary layers and the red dots (points 1 and 2) represent the edge of the boundary layers at the trailing-edge $x$ station. The potential-flow theory, ignoring the boundary layers, assumes that points 1 and 2 and the trailing edge all lie on top of each other. Hence, the Kutta condition (continuous pressure at the trailing edge) would dictate $P_{1}=P_{2}$ neglecting the pressure rise across the boundary layer and its effect on the circulation development over the airfoil.

The Kutta condition at the sharp trailing edge can be stated in several ways, such as (i) finite velocity, (ii) zero loading or (iii) continuous pressure, among others (see Sears Reference Sears1956). In fact, some of these representations are, indeed, exact. For example, clearly, the pressure must be continuous at the trailing edge (TE). That is,

(6.5) $$\begin{eqnarray}\lim _{y\rightarrow 0^{+}}P(TE,y)=\lim _{y\rightarrow 0^{-}}P(TE,y).\end{eqnarray}$$

However, the inviscid pressure distribution over the plate represents the distribution at the edge of the boundary layer. That is, applying the condition (6.5) within the framework of potential-flow results in $P_{1}=P_{2}$ , where the points 1 and 2 lie on the edge of the boundary layer at the trailing-edge station as shown in figure 19. While Prandtl’s boundary layer assumption (pressure is constant along the boundary layer thickness) is valid over the majority of the airfoil length, it is not necessarily valid in the singular trailing-edge region. As such, if $\unicode[STIX]{x0394}P$ is the pressure rise across the boundary layer, then the condition (6.5) results in

(6.6) $$\begin{eqnarray}P_{1}-\unicode[STIX]{x0394}P_{1}=P_{2}-\unicode[STIX]{x0394}P_{2},\end{eqnarray}$$

which is also suggested by Preston (Reference Preston1943) and Spence (Reference Spence1954) as a modification of the classical Kutta condition ( $P_{1}=P_{2}$ ). Since the points 1 and 2 lie on the edge of the boundary layer, one can use the unsteady Bernoulli’s equation to relate $P_{1}$ and $P_{2}$ as (Bisplinghoff et al. Reference Bisplinghoff, Ashley and Halfman1996; Katz & Plotkin Reference Katz and Plotkin2001)

(6.7) $$\begin{eqnarray}\frac{P_{1}}{\unicode[STIX]{x1D70C}}+\frac{1}{2}V_{1}^{2}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}t}=\frac{P_{2}}{\unicode[STIX]{x1D70C}}+\frac{1}{2}V_{2}^{2}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x2202}t},\end{eqnarray}$$

where $V$ is the potential-flow velocity at the edge of the boundary layer and $\unicode[STIX]{x1D719}$ is the corresponding velocity potential. Combining (6.6) and (6.7) and realizing that $\unicode[STIX]{x1D719}_{1}-\unicode[STIX]{x1D719}_{2}=\unicode[STIX]{x1D6E4}$ , one obtains

(6.8) $$\begin{eqnarray}\dot{\unicode[STIX]{x1D6E4}}=\frac{1}{2}(V_{2}^{2}-V_{1}^{2})+\frac{\unicode[STIX]{x0394}P_{2}-\unicode[STIX]{x0394}P_{1}}{\unicode[STIX]{x1D70C}}.\end{eqnarray}$$

Equation (6.8) represents an exact (derived) version of the hypothesized Kutta condition. In particular, it governs the evolution of the bound circulation over the airfoil, i.e. it provides the dynamics of the bound circulation. Interestingly, it can be derived from a completely different point of view than the continuity argument (6.5) which we opt to show below. The underpinning concept is that the circulation is instantaneously conserved. That is, the instantaneous rate of change of circulation $\dot{\unicode[STIX]{x1D6E4}}$ is related to the total vorticity flux at separation (Sears Reference Sears1976). As such, assuming that separation occurs at the trailing edge (complying with the triple-deck theory used in this paper), we write

(6.9) $$\begin{eqnarray}\dot{\unicode[STIX]{x1D6E4}}=-\left[\int _{0}^{\unicode[STIX]{x1D6FF}_{1}}\unicode[STIX]{x1D701}(y)u(y)\,\text{d}y+\int _{-\unicode[STIX]{x1D6FF}_{2}}^{0}\unicode[STIX]{x1D701}(y)u(y)\,\text{d}y\right],\end{eqnarray}$$

where $\unicode[STIX]{x1D701}$ is the clockwise vorticity, $u$ is the velocity parallel to the wall inside the boundary layer and $\unicode[STIX]{x1D6FF}$ is the boundary layer thickness. Also, $\unicode[STIX]{x1D6E4}$ is assumed clockwise positive in the entire paper. Then, one can use the boundary layer theory along a curved surface (Goldstein Reference Goldstein1938, pp. 119–120; Sears Reference Sears1956) to write

(6.10) $$\begin{eqnarray}\int _{0}^{\unicode[STIX]{x1D6FF}_{1}}\unicode[STIX]{x1D701}u\,\text{d}y=\int _{0}^{\unicode[STIX]{x1D6FF}_{1}}\left(\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}y}+\unicode[STIX]{x1D705}u\right)u\,\text{d}y=\frac{V_{1}^{2}}{2}+\frac{\unicode[STIX]{x0394}P_{1}}{\unicode[STIX]{x1D70C}},\end{eqnarray}$$

where $\unicode[STIX]{x1D705}$ is the curvature of the wall. Writing an expression for the vorticity flux out of the boundary layer on the lower surface similar to (6.10) and substituting both in (6.9), one immediately arrives at the condition (6.8).

Setting $\unicode[STIX]{x0394}P_{1}=\unicode[STIX]{x0394}P_{2}=0$ along with $V_{1,2}=U\pm (1/2)\unicode[STIX]{x1D6FE}_{TE}$ as typically done in the classical theory of unsteady aerodynamics, the condition (6.8) results in

(6.11) $$\begin{eqnarray}\dot{\unicode[STIX]{x1D6E4}}_{Kutta}(t)=-U\unicode[STIX]{x1D6FE}_{TE}(t),\end{eqnarray}$$

where $\unicode[STIX]{x1D6FE}_{TE}$ is the circulation distribution at the trailing edge (instantaneous strength of the shed vortex sheet per unit length at the shedding time). Equation (6.11) is equivalent to the classical Kutta condition ( $P_{1}=P_{2}$ ) and it is ubiquitously used in the classical theory of unsteady aerodynamics (Wagner Reference Wagner1925; Loewy Reference Loewy1957; Bisplinghoff et al. Reference Bisplinghoff, Ashley and Halfman1996, equations (5)–(318); Peters Reference Peters2008, equation (11c); among others). Note that the main difference between the exact condition (6.8) and the classical Kutta condition (6.11) is two assumptions: (i) linearization ( $V_{1,2}=U\pm (1/2)\unicode[STIX]{x1D6FE}_{TE}$ ); and (ii) neglect of the viscous terms $\unicode[STIX]{x0394}P_{1}$ and $\unicode[STIX]{x0394}P_{2}$ . The first assumption may be valid for small disturbance (small angle of attack). Using higher-fidelity computational simulations, we assess the validity of the second assumption. Figure 20 shows a comparison between the Kutta’s rate of circulation development $\dot{\unicode[STIX]{x1D6E4}}_{Kutta}$ and the viscous contribution proportional to $\unicode[STIX]{x0394}P_{TE}=P_{2}-P_{1}$ , at different frequencies and Reynolds numbers. It is found that the viscous contribution $\dot{\unicode[STIX]{x1D6E4}}_{\unicode[STIX]{x0394}P}$ relative to the inviscid $\dot{\unicode[STIX]{x1D6E4}}_{Kutta}$ is not sensitive to frequency; it depends mainly on Reynolds number. This ratio is approximately 18 %, irrespective of the frequency $k$ , at the lower Reynolds number $R=10^{4}$ versus 6–7 % at $R=10^{5}$ . Moreover, a higher frequency, though it does not significantly affect the magnitude, causes a significant phase shift for the viscous contribution, which will in turn affect the phase of the total lift force at high frequencies.

In summary, the lower the Reynolds number, the larger is the viscous contribution to the bound circulation development relative to the inviscid one; and the higher the frequency, the larger is the phase shift of this viscous contribution. That is, at higher frequencies and lower Reynolds numbers, the viscous contribution to the bound circulation rate of development is of relatively larger magnitude and phase shift. This point, in addition to the wake viscous damping discussed above, may help to explain the physical reasons behind the additional phase lag in the lift response at high $k$ and low $R$ that could not be captured by the inviscid theory even at very small amplitudes. Since this viscous contribution is essentially neglected in the classical potential-flow framework by virtue of the classical Kutta condition (6.11), it is inferred that the Kutta condition is one of the reasons behind the inaccurate phase prediction of Theodorsen’s lift frequency response function.

Figure 20. Inviscid and viscous contributions ( $\dot{\unicode[STIX]{x1D6E4}}_{Kutta}$ , $\dot{\unicode[STIX]{x1D6E4}}_{\unicode[STIX]{x0394}P}$ ) to the rate of bound circulation development over NACA 0012 undergoing a pitching oscillation about the quarter-chord point at different reduced frequencies and Reynolds numbers: (a)  $R=10^{5}$ and $k=0.1$ ; (b)  $R=10^{5}$ and $k=1$ ; (c)  $R=10^{4}$ and $k=0.3$ ; (d)  $R=10^{4}$ and $k=1$ .

6.4 Viscous reduction in virtual mass

A closer look at the viscous correction $B_{v}$ in (3.13) and its contribution to the circulatory lift in (4.10) implies that the viscous lift contribution is proportional to (in phase with) the ‘effective’ angle of attack

(6.12) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}_{eff}=\frac{1}{U^{2}}\left[\frac{1}{2}a_{0}(t)+2\mathop{\sum }_{n=1}^{\infty }na_{n}(t)\right].\end{eqnarray}$$

Note that this $\unicode[STIX]{x1D6FC}_{eff}$ is different from the common notion of the effective angle of attack in potential flow. The former is a term special to the developed theory while the latter is simply given by the angle of attack $\unicode[STIX]{x1D6FC}_{3/4}$ at the three-quarter-chord point (Schlichting & Truckenbrodt Reference Schlichting and Truckenbrodt1979, p. 80). Equations (3.13) and (4.10) imply that the viscous contribution $\unicode[STIX]{x0394}_{v}C_{L}$ to the lift coefficient is simply proportional to $\unicode[STIX]{x1D6FC}_{eff}$ :

(6.13) $$\begin{eqnarray}\unicode[STIX]{x0394}_{v}C_{L}=-2\unicode[STIX]{x03C0}\tilde{B}_{v}(t)C(k)=4\unicode[STIX]{x03C0}\unicode[STIX]{x1D716}^{3}\unicode[STIX]{x1D706}^{-5/4}B_{e}(\unicode[STIX]{x1D716}^{-1/2}\unicode[STIX]{x1D706}^{-9/8}|\unicode[STIX]{x1D6FC}_{eff}(t)|)\unicode[STIX]{x1D6FC}_{eff}(t)C(k),\end{eqnarray}$$

where $B_{e}(\cdot )$ is a nonlinear function coming from the numerical solution of Chow & Melnik (Reference Chow and Melnik1976) to the triple-deck problem, specifically from figure 3(b). This nonlinear function mainly affects the magnitude of $\tilde{B}_{v}$ with a weak effect on its phase, as shown in figures 21 and 22: the viscous contribution is in phase with  $\unicode[STIX]{x1D6FC}_{eff}$ .

Figure 21. A comparison between the time history of the inviscid circulatory lift $\unicode[STIX]{x1D6FC}_{3/4}C(k)$ , the effective angle of attack $\unicode[STIX]{x1D6FC}_{eff}$ for the developed viscous theory, the weighted effective angle of attack or the viscous contribution $-\tilde{B}_{v}C(k)$ , and the inviscid added-mass lift. All lift coefficients are represented as effective angles of attack (i.e. normalized by $2\unicode[STIX]{x03C0}$ ). Simulation of the developed viscous model is performed for a flat plate pitching about its quarter-chord point with amplitude $1^{\circ }$ at $k=1$ and $R=10^{4}$ . The viscous contribution $-\tilde{B}_{v}C(k)$ is opposite to the inviscid non-circulatory lift, decreasing the added-mass effect by 70 %.

Considering the studied case of a pitching flat plate about its quarter-chord point, the above definition of the effective angle of attack can be manipulated to write its Fourier transform as

(6.14) $$\begin{eqnarray}\hat{\unicode[STIX]{x1D6FC}}_{eff}=\hat{\unicode[STIX]{x1D6FC}}_{3/4}C(k)-A_{\unicode[STIX]{x1D6FC}}(3\text{i}k-2k^{2}).\end{eqnarray}$$

Since the viscous contribution is proportional to $\unicode[STIX]{x1D6FC}_{eff}$ , it is fair to write

(6.15) $$\begin{eqnarray}\unicode[STIX]{x0394}_{v}C_{L}=A[\hat{\unicode[STIX]{x1D6FC}}_{3/4}C(k)-A_{\unicode[STIX]{x1D6FC}}(3\text{i}k-2k^{2})]C(k),\end{eqnarray}$$

where $A=4\unicode[STIX]{x03C0}\unicode[STIX]{x1D716}^{3}\unicode[STIX]{x1D706}^{-5/4}B_{e}(\unicode[STIX]{x1D716}^{-1/2}\unicode[STIX]{x1D706}^{-9/8}|\unicode[STIX]{x1D6FC}_{eff}(t)|)$ . Written this way, equation (6.15) shows how the viscous contribution induces more lag to the inviscid circulatory lift $\hat{\unicode[STIX]{x1D6FC}}_{3/4}C(k)$ : it provides a similar contribution ( $\hat{\unicode[STIX]{x1D6FC}}_{3/4}C(k)$ ) that is multiplied by $C(k)$ even further (i.e. possesses double the lag of Theodorsen’s function). Also, $A$ provides more lag because it depends algebraically on $\unicode[STIX]{x1D6FC}_{eff}$ , which includes $C(k)$ in it (i.e. possesses lag).

Recalling the inviscid non-circulatory lift in this case,

(6.16) $$\begin{eqnarray}{\hat{C}}_{L_{NC}}=\unicode[STIX]{x03C0}A_{\unicode[STIX]{x1D6FC}}(\text{i}k-k^{2}/2),\end{eqnarray}$$

equation (6.15) implies that the second component $-A_{\unicode[STIX]{x1D6FC}}(3\text{i}k-2k^{2})$ of the viscous contribution is opposite to the non-circulatory lift (added mass). This fact is clearly seen in figures 21 and 22. Therefore, adding $\unicode[STIX]{x0394}_{v}C_{L}$ to the inviscid lift coefficient would decrease the added mass and cause a phase lag for its contribution.

Figure 22. Argand diagram showing different components of lift for a pitching flat plate about its quarter-chord point at $k=1$ and two different Reynolds numbers, (a)  $R=10^{5}$ and (b)  $R=10^{4}$ . All lift coefficients are represented as effective angles of attack (i.e. normalized by $2\unicode[STIX]{x03C0}$ ). The term $\unicode[STIX]{x1D6FC}_{eff}$ is scaled down to one-quarter size to enhance visualization. The viscous contribution $\tilde{B}_{v}$ increases as $R$ decreases, resulting in a larger phase difference between the inviscid circulatory contribution $\unicode[STIX]{x1D6FC}_{3/4}C(k)$ and the total (viscous) circulatory component $C_{L_{c}}$ , or a larger decrease in the added-mass effect.

It is noteworthy to comment on the results shown in figure 21 for a pitching flat plate about its quarter-chord point with amplitude $1^{\circ }$ at $k=1$ and $R=10^{4}$ . In addition to the points addressed above ( $-\tilde{B}_{v}C(k)$ is in phase with $\unicode[STIX]{x1D6FC}_{eff}$ and out of phase with $C_{L_{NC}}$ ), it is interesting to see the apparent nonlinear response of the viscous contribution even at this very small amplitude. Moreover, the phase lag between the inviscid circulatory lift and viscous contribution is very clear also in the Argand diagrams in figure 22. More importantly, the viscous contribution is as strong as the inviscid one at this low $R$ and high $k$ : its magnitude is 60 % of the inviscid circulatory lift or 82 % of the non-circulatory lift, which would have a significant effect on flutter (Bisplinghoff et al. Reference Bisplinghoff, Ashley and Halfman1996), if it happens at this low $R$ and high  $k$ .