Contents

Theodorsen Theory

Theodorsen Theory

“…why linear systems are so important. The answer is simple: because we can solve them!”

Richard Feynman

Theodore Theodorsen was first to publish a complete and detailed solution to the problem of aerodynamic loads on an oscillating lifting surface. Chapter 5 in [BAH96] is a great resource for understanding the theoretical foundations of the Theodorsen model (even better than the original paper by Dr Theodorsen that describes the model [The40]). That said, the complete derivation is mathematically involved and requires meticulously solving a lot of complex-valued integrals. Here we will start with the final results proposed by Theodorsen and work our way to backwards with the motivation to understand the expressions involved intuitively and take assummptions involved at face value.

The Theodorsen theory, or alternately the Theodorsen model, has been developed in the frequency domain i.e. the variable of interest is the frequency of input. The input in this case refers to the frequency at which the airfoil is moving. The Theodorsen theory only accounts for the pitching and heaving motions of the airfoil. The lead-lag motion of the airfoil, or alternately the variation in the oncoming flow velocity is not accounted for.

Aside: Linear time-invariant systems (LTI)

The behavior of a dynamical system ccan be mathematically represented using the input forcing \(x(t)\), the output response \(y(t)\) and the operator \(h\) that operates on the input.

If we just focus on the cases where the operator $h$ is given by only a combination of the following -
In this case, the system is expected to have the following characteristics-

The discussion above and all the corresponding figures have been obtained from

[source]

.

Aside: Frequency Response

In the analysis of dynamical systems there has been sufficient progress so much so that there are multiple domains of modelling strategies that are available for one to choose from - frequency response, state-space, time-domain and impulse response. All these strategies apply to linear systems only and are equivalent in that any one form can be obtained from any one of the other forms. Dynamical systems are often studied in the frequency domain because the linearity assumption means that solving for the system response becomes especially straightforward once the transfer function \(G(s)\) of the system is known. \(G(s)\) relates the input forcing of the system at a given frequency to the output response of the system at the same frequency but now at a different amplitude as well as a phase shift.

Bode plot [source]

\[s=i\omega\]
\[\frac{B}{A}=\left | G(s) \right | \textrm{ dB}\]
\[\phi=\text{arg}(G(s)) \textrm{ deg}\]

Mathematical formulation

Airfoil+trailing-edge flap model in Theodorsen model [source]

The schematic model of the airfoil with trailing-edge used in the Theodorsen model is shown above. The airfoil is placed such that the origin coincides with the mid-chord of the airfoil. The total chord length is \(2b\) where \(b\) is referred to as the semi-chord. This is done because quite often it is the semi-chord that keeps popping up in expressions rather than the total chord itself. The symbol \(c\) is reserved for the distance of the flap hinge from the mid-chord. The oncoming flow velocity is assumed to be constant within Theodorsen model and only steady-state aerodynamic pressure distribution over the airfoil due harmonic motion can be solved - no transients allowed. The centers of gravity of the airfoil and the trailing-edge flap are relevant here because Theodorsen went on to use the results of the model to analyse flutter.

Final equations

The total lift acting on the complete airfoil (i.e. together with flap) is

\[\begin{equation*} \begin{split} L= \underbrace{\rho b^2\left(U \pi \dot{\alpha}+\pi \ddot{h}-\pi b a \ddot{\alpha}-U T_4 \dot{\beta}-T_1 b \ddot{\beta}\right)}_\text{non-circulatory} + \underbrace{2 \pi \rho U b C(k)\left\{U \alpha+\dot{h}+ b\left(\frac{1}{2}-a\right) \dot{\alpha}+\frac{1}{\pi} T_{10} U \beta\right. \left.+b \frac{1}{2 \pi} T_{11} \dot{\beta}\right\}}_\text{circulatory} \end{split} \end{equation*}\]

Moment acting about the airfoil pitching location \(x=ba\) is

\[\begin{equation*} \begin{split} M = & \underbrace{-\rho b^2\left[\pi\left(\frac{1}{2}-a\right) U b \dot{\alpha}+\pi b^2\left(\frac{1}{8}+a^2\right) \ddot{\alpha} +\left(T_4+T_{10}\right) U^2 \beta +\left(T_1-T_8-(c-a) T_4+\frac{1}{2} T_{11}\right) U b \dot{\beta} -\left(T_7+(c-a) T_1\right) b^2 \ddot{\beta}-a \pi b \ddot{h}\right]}_\text{non-circulatory} \\ & \underbrace{+2 \rho U b^2 \pi\left(a+\frac{1}{2}\right) C(k)\left[U \alpha+\dot{h}+b\left(\frac{1}{2}-a\right) \dot{\alpha}\right. \left.+\frac{1}{\pi} T_{10} U \beta+b \frac{1}{2} T_{11} \dot{\beta}\right]}_\text{circulatory} \end{split} \end{equation*}\]

The final equations are a result of the airfoil kinematics and the complex valued function \(C(k)\), commonly referred to as the Theodorsen function in literature. The former invloves reactionary loads experienced by the moving airfoil from pushing the fluid around and the latter manifests the time history of the oscillation, thereby incorporating hysteresis (effects of the history of motion) in the aerodynamic loads experienced by an airfoil. While the Theordorsen model involves an airfoil with trailing-edge flap, for simplicity we would look at the derivation below only for the case of an airfoil, i.e. no \(\beta\) terms. Note that the \(T_n\) terms are all attached to the \(\beta\) terms so we need not be focussed on what those are but for the curious minded these are different integrals that are laid down in Ref [The40].

Relevant building blocks

Understanding the Theodorsen theory requires a refresher of a number of theoretical aerodynamics concepts.

  • Potential flow theory

  • Kutta condition

  • No-penetration boundary condition

  • Conservation of circulation

  • Conformal mapping

    • Joukowski transformation

  • Unsteady Bernoulli’s equation

All these concepts have already been discussed in sufficient detail so that we can start using them now for the practical problem of an oscillating airfoil.

Non-circulatory lift

So the problem boils down to solving the Laplacian for the flow around a circular cross-section and then using the conformal mapping to obtain the flow around the thin airfoil profile. Thereafter, once the potential function in the z-plane is obtained only the corresponding potential function value at the airfoil surface is of interest since it, and its derivates, are directly related to the aerodynamic pressure forces using the unsteady Bernoulli’s equation. The pressure forces can be integrated over the entire airfoil to obtain lift force and moment about the quarter chord.

Potential flow theory and conformal mapping

The Theodorsen’s model begins with the conformal representation of a flat plate (airfoil) as a circle. The idea is to place sources/sinks around the circle so that a streamline is obtained that imitates flow around a circle (2D cylinder) placed in a fluid. Thereafter, using conformal mapping techniques, the fluid velocity distribution around the airfoil can be obtained.

Theodorsen begins by placing sources in the upper half of the cylinder and sinks of equal strength in the lower half placed symmetrically about X-axis as shown below. This maps onto the \(z\)-plane at a single location (\(x\),0). Given the nature of the Joukowski transformation, and without going into further details, the source/sink pair doesn’t simply cancel the effect of each other.

Conformal representation of sources/sinks in the $\xi$-plane and the $z$-plane.

Recall that the complex function \(F(\xi)\) can be written down for the case of a source/sink. In case the source/sink is not exactly as the origin, an offset can be applied as

\[F(\xi) = \frac{\frac{b}{2}q(\varphi)d\varphi}{2\pi}\text{ln} (\xi-\frac{b}{2}e^{i\varphi}) - \frac{\frac{b}{2}q(\varphi)d\varphi}{2\pi}\text{ln} (\xi-\frac{b}{2}e^{-i\varphi})\]

The corresponding velocity in the complex plane at some location \(re^{i\theta}\) off the cylinder is

\[\begin{equation*} \begin{split} \bar{W}=\frac{d F}{d \xi}&=\frac{b q(\varphi) d \varphi}{4 \pi}\left[\frac{1}{\xi-\frac{b}{2} e^{i\varphi}}-\frac{1}{\xi-\frac{b}{2} e^{-i\varphi}}\right]\\ &=\frac{ b^2 q(\varphi) d \varphi}{4 \pi} \frac{2 i \sin \varphi}{r^2 e^{2 i \theta}-\frac{b}{2} r e^{i \theta}(2 \cos \varphi)+\frac{b^2}{4}} \end{split} \end{equation*}\]
\[\bar{W}=V_{\text{real}}+iV_{\text{imag}} = (dV_{r}-idV_{\theta})e^{-i\theta}\]

where \(\xi=re^{i\theta}\), \(dV_{r}\) and \(dV_{\theta}\) denote the contribution to \(V_r\) and \(V_\theta\) at location \(\xi\) due to the source/sink pair on the circle at \(\frac{b}{2}e^{i\varphi}\) and \(\frac{b}{2}e^{-i\varphi}\). An integration over all such source/sink pairs over the entire cylinder will give the net radial and tangential flow velocities at any point in the domain. In particular, the tangential flow velocity at the cylinder surface i.e. at \(r=\frac{b}{2}\) simplifies to

\[V_\theta(\frac{b}{2},\theta) = -\frac{1}{2\pi }\int_{0}^{\pi }\frac{q(\varphi) \sin\varphi d\varphi}{\cos \theta -\cos \varphi}\]

Similarly, the expression for \(V_r(\frac{b}{2},\theta)\) can be obtained but some additional treatment is required so only the final expression will be stated here. Those interested in the details can check out Ref [BAH96] (Section 5-6).

\[V_r(\frac{b}{2},\theta) = \frac{q(\theta)}{2}\]

So, based on som unknown distribution of sources and sinks \(q(\theta)\) on the cylinder surface, velocity distribution in the entire domain can be obtained including the surface of the cylinder itself. Much like the case of a doublet in a flow, we can obtain \(q(\theta)\) by satisfying the condition that normal flow at the boundary of the cylinder is zero.

No-penetration boundary condition

The velocity field at the cylinder surface is a function of the source/sink distribution strength \(q(\theta)\). This is an unknown quantity but knowledge about the motion of the airfoil can help arrive at a unique distribution using the no-penetration boundary condition.

Representative airfoil motion involving pitching and plunging together with forward motion.

In the airfoil domain, the velocity normal to the airfoil surface due to its motion is given by

\[v_{\text {motion }}(x, y=0)= v_{\text {airfoil }}=-U_\infty \sin\alpha-\dot{h} \cos \alpha-\dot{\alpha}(x-a b) \]

Theodorsen theory is concerned with only small motions of the airfoil such that the small angle approximation can be applied.

\[\Rightarrow v_{\text {motion }}(x, y=0)= v_{\text {airfoil }}\approx-U_\infty \alpha-\dot{h} -\dot{\alpha}(x-a b) \]

This is related to the radial velocity in the \(\xi\) domain on the surface of the cylinder by the relation-

\[V_r(\frac{b}{2},\theta) = 2\sin \theta v_{motion}\]

Putting it all together

\[V_r(\frac{b}{2},\theta) = \frac{q(\theta)}{2}\]
\[V_r(\frac{b}{2},\theta) = 2\sin \theta v_{motion}\]
\[\Rightarrow q(\theta) = 4\sin\theta v_{motion}\]

Once the distribution of source/sinks is known around the cylinder, based on the motion of the airfoil, then the tangential flow velocity at the cylinder surface \(V_\theta(\frac{b}{2},\theta)\) can be evaluated

\[\begin{equation*} \begin{split} V_\theta(\frac{b}{2},\theta) &= -\frac{1}{2\pi }\int_{\theta}^{\pi }\frac{q(\varphi) \sin\varphi d\varphi}{\cos \theta -\cos \varphi}\\ &=-\frac{1}{2\pi }\int_{\theta}^{\pi }\frac{4\sin\varphi v_{motion} \sin\varphi d\varphi}{\cos \theta -\cos \varphi} \end{split} \end{equation*}\]

Substituting for \(v_motion\) brings the integral in the form of the so called Galuert’s integral \(I_n\), commonly used in thin airfoil theory results, provided in [Gla36](Pg 92-93).

\[I_n=\int_0^\pi \frac{\cos n \theta}{\cos \theta-\cos \phi} d \theta=\pi \frac{\sin n \phi}{\sin \phi} .\]

Again skipping the intermediate steps, the final result obtained as a consequence of carrying out the above integration is as follows-

\[V_\theta(\frac{b}{2},\theta)=2\left[(U \alpha+\dot{h}-\dot{\alpha} a b) \cos \theta+\frac{\dot{\alpha}b}{2} \cos 2 \theta\right]\]

Recall from the last chapter that,

\[\phi_{cyl}(\theta) =\frac{-b}{2}\int_{\theta}^{\pi}V_\theta d\varphi \]
\[\begin{equation*} \begin{split} \Rightarrow \phi_{cyl}(\theta) &=\frac{-b}{2}\int_{\theta}^{\pi}2\left[(U \alpha+\dot{h}-\dot{\alpha} a b) \cos \varphi+\frac{\dot{\alpha}b}{2} \cos 2 \varphi\right] d\varphi \\ &=b\left[(U \alpha+\dot{h}-\dot{\alpha} a b) \sin \theta+\frac{\dot{\alpha}b}{4} \sin 2 \theta\right] \end{split} \end{equation*}\]

We now have all the elements for obtaining lift generated by the airfoil undergoing \(V_{motion}\). But before solving for the expression of total lift we need \(\phi\) at the trailing-edge and \(\frac{\partial \phi}{\partial t}\)

\[\phi_{TE} = \phi_{cyl}(\theta =0) = 0\]
\[\frac{\partial \phi(\theta)}{\partial t} = b\left[(U \dot\alpha+\ddot{h}-\ddot{\alpha} a b) \sin \theta+\frac{\ddot{\alpha}b}{2} \sin 2 \theta\right]\]

Note that the Theodorsen theory assumes a constant oncoming flow velocity \(U\), however if the oncoming flow velocity is also time-varying then the above expression for \(\frac{\partial \phi}{\partial t}\) changes accordingly. Recall that a general airfoil section on rotor blade undergoes pitching, plunging and lead-lag motion so a general unsteady model that incorporates time-varying oncoming flow velocity is more relevant for rotocraft applications. Additionally, during forward flight motion there is a large time-varying component of the oncoming flow velocity too. Such a model exists and was proposed by Greenberg [Gre47] a few years after the Theodorsen’s model. However we continue to focus on the Theodorsen model for its relative simplicity.

From the last chapter, recall that

\[p-p_\infty = -\rho\left [ \frac{\partial \phi}{\partial t}+U_\infty \frac{\partial \phi}{\partial x} \right ] \]

From the expression for \(\phi(\theta)\), we can see that for points on the surface of the cylinder that form a mirror image about the \(X\) axis, the potential functions at the two points are related by

\[\phi(\theta) = -\phi(-\theta)\]

Therefore the total lift generated over the airfoil surface spanning from \(-b\) to \(+b\) is given by-

\[\begin{equation*} \begin{split} L &= \int_{-b}^{b}p_l-p_u dx \\ &= \int_{-b}^{b}2\rho\left [ \frac{\partial \phi}{\partial t}+U_\infty \frac{\partial \phi}{\partial x} \right ] dx \end{split} \end{equation*}\]

Again, assuming \(\phi_{LE}=0\) the lift expression simplifies to

\[L=2\rho\int_{-b}^{b} \frac{\partial \phi}{\partial t} dx +U_\infty \phi_{TE}\]

what we’ve ended up with is the non-circulatory component of lift i.e. aerodynamic force perpendicular to the oncoming flow that is not due to circulation. Here we start differentiating this from the total lift \(L\) by denoting this using \(L_{NC}\)

\[L_{NC} = \pi \rho b^2 \underbrace{\left[(U \dot{\alpha}+\ddot{h}-\ddot{\alpha} a b) \right]}_\textrm{acceleration at mid-chord}\]

Salient features of \(L_{NC}\)

  • Circulation not needed for lift when unsteadiness involved

    • only sources/sinks distribution lead to \(L_{NC}\)

    • \(L_{NC}=0\) if \(\frac{\partial \phi}{\partial t}=0\)

  • Also refered to as lift due to added mass
    \(L_{NC} = \underbrace{\pi \rho b^2}_\text{added mass} \underbrace{\left[(U \dot{\alpha}+\ddot{h}-\ddot{\alpha} a b) \right]}_\text{acceleration at mid-chord}\)

    • due to acceleration of surrounding fluid as the airfoil moves

    • the amount of fluid influenced depends on size of airfoil

  • \(L_{NC}\) is the reactionary inertial force acted on the airfoil by the fluid

    • airfoil lift acts opposite to the direction of airfoil acceleration

  • Net \(L_{NC}\) over one actuation cycle is 0

There are some additional results worth noting based on the derivation so far. The velocity in the cylinder domain at the point corresponding to the airfoil trailing-edge is given by

\[\begin{equation*} \begin{split} V_\theta(\frac{b}{2},0)&=2\left[(U \alpha+\dot{h}-\dot{\alpha} a b) +\frac{\dot{\alpha}b}{2} \right] \\ &=2\left[-v_{c/2}+\frac{\dot{\alpha}b}{2}\right] \\ &=-2v_{3c/4} \end{split} \end{equation*}\]

Alternately,

\[V_\theta(\frac{b}{2},0)_{NC}=-2v_{3c/4}\]
\[ \left. u(x,y) \right|_{TE} = u(b,0)=u_{TE} = -\frac{V_\theta(\frac{b}{2},0)}{2\sin 0} = \frac{v_{3c/4}}{0} \rightarrow \infty\]

This is a problem because we’ve ended up with an unphysical result at the trailing-edge. But this is no cause for concern yet because the lift generated from the circulatory component of flow has not been accounted for yet and this ‘issue’ can perhaps be fixed afterall.

Circulatory lift

Next we quantify the lift generated as a consequence of circulation around the airfoil.

Conformal representation of vortices in the $\xi$-plane and the $z$-plane.

Conservation of circulation

That the dominant mechanism of lift generation by an airfoil is by introducing circulation has already been discussed. Any change in this circulation proceeds with shedding of an equivalent vortex of opposite strength in the wake that then moves away from the airfoil with the freestream. So clearly one needs to solve the flowfield for the situation when there is a vortex in the wake. We do this in the cylinder domain but now care needs to be taken because the addition of vortices that do not lie at the center of the cylinder as a streamline already established via placement of sources/sinks.

No-penetration boundary condition

The no-penetration boundary condition was already satisfied when solving for the non-circulatory flow. In the Theodorsen’s model this boundary condition is continued to be satisfied by judicious placement of \(-\Gamma\) such that the surface of the cylinder continues to be a streamline i.e. there is no flow perpendicular to it. The distance from the center of cylinder where the vortex needs to be placed is

\[d = \frac{b^2}{4X}\]

Again, the detailed derivation is skipped in the interest of conciseness.

Kutta condition

The vortex strength \(\Gamma\) is unknown but in order for the Kutta condition to be satisfied at the trailing edge, i.e. \(u_{TE}\) is finite, the only possibility is that \(V_\theta(\frac{b}{2},0) = 0\) i.e.

\[V_\theta(\frac{b}{2},0)_{NC} + V_\theta(\frac{b}{2},0)_{C} =0\]

Putting it all together

\[F=\frac{\Gamma}{2 \pi i}[\ln (\xi-X)-\ln (\xi-d)]\]
\[\bar{W}=\frac{d F}{d \xi}=\frac{\Gamma}{2 \pi i}\left[\frac{1}{\xi-X}-\frac{1}{\xi-d}\right]=\frac{1}{2 \pi} \frac{X-d}{\xi^2-\xi(X+d)+X d} = V_r -i V_\theta\]

From the complex velocity conjugate result above, the required \(d\) that maintains the cylinder as a streamline is obtained using \(V_r=0\). However, we’ve already stated this result to be \(d = \frac{b^2}{4X}\). This gives the tangential flow velocity due to the pair of vortices

\[V_\theta(\frac{b}{2},\theta) = \frac{-\Gamma}{b \pi} \frac{X^2-b^2}{X^2+\frac{b^2}{4}-b X \cos \theta}\]
\[\Rightarrow V_{\theta,C}(\frac{b}{2},0) = \frac{-\Gamma}{b \pi} \frac{X^2-b^2}{X^2+\frac{b^2}{4}-b X } = \frac{-\Gamma}{\pi b} \frac{X+\frac{b}{2}}{X-\frac{b}{2}}\]

Using the Joukowski transformation expression the corresponding expression for tangential flow velocity using the coordinates of the vortex in the airfoil domain becomes

\[ V_{\theta,C}(\frac{b}{2},0) = \frac{-\Gamma}{\pi b}\sqrt{ \frac{z+b}{z-b}}\]

further considering the fact that in case of unsteady motion of the airfoil instead of a discrete starting vortex-like wake there would exist a distribution of vorticity \(\gamma(z,t)\),

\[ V_{\theta,C}(\frac{b}{2},0) = -\frac{1}{\pi b}\int_{b}^{\infty} \gamma_w(z,t)\sqrt{ \frac{z+b}{z-b}}dz\]

Conformal representation of a flat plate (airfoil) as a circle within the Theodorsen model.

Much like the steady aerodynamics case, the Kutta condition in the unsteady aerodynamics case needs to be satisfied which takes the form

\[-2v_{3c/4} -\frac{1}{\pi b}\int_{b}^{\infty} \gamma_w(z,t)\sqrt{ \frac{z+b}{z-b}}dz = 0\]

This can be alternately written as

\[ \frac{1}{ U b}\int_{b}^{\infty} \gamma_w(z,t)\sqrt{ \frac{z+b}{z-b}}dz = 2\pi\alpha_{3c/4}(t)\]

This is the integral circulation equation. The quantity on the RHS is known (since the motion of the airfoil is known). Solving for \(\gamma_w\), however, is not straightforward since it is an integral equation. Note also that the RHS corresponds to quasi-steady lift expression i.e. \(L_{QS}(t)=2\pi\alpha_{3c/4}(t)\)

Alternately, we can carry out the same exercise as in the case of the non-circulatory part and evaluate \(\phi_{cyl,C}(\theta)\), i.e. the circulatory component of \(\phi\) at the cylinder surface. Thereafter, once \(\phi_{cyl,C}\) at the trailing-edge and \(\frac{\partial \phi_{cyl,C}(\theta)}{\partial t}\) are known, the lift force due to the presence of circulation can be obtained. Here we again go back to considering only a single vortex \(\Gamma\) in the wake. Note that \(\Gamma\) is the strength of the vortex (which is constant since there is no viscous dissipation in potential flow theory) in the wake which moves with the flow, so it is the position of the vortex \(z\) which is a function of time \(\rightarrow dz/dt = U\).

\[\begin{equation*} \begin{split} \phi_{cyl,C}(\theta) &=\frac{-b}{2}\int_{\theta}^{\pi}V_{\theta,C} d\varphi \\ &=\frac{-b}{2}\int_{\theta}^{\pi}\frac{-\Gamma}{b \pi} \frac{X^2-b^2}{X^2+\frac{b^2}{4}-b X \cos \varphi} d\varphi \\ & =\frac{\Gamma}{\pi } \tan ^{-1}\left(\sqrt{\frac{z+b}{z-b}} \tan \frac{\theta}{2}\right)_\theta^\pi \\ & =\frac{\Gamma}{\pi }\left[\frac{\pi}{2}-\tan ^{-1} \sqrt{\frac{z+b}{z-b}}\sqrt{\frac{1-\cos \theta}{1+\cos \theta}}\right]. \end{split} \end{equation*}\]
\[\phi_{TE,C} = \phi_{cyl,C}(0) = \frac{\Gamma}{2}\]
\[\begin{equation*} \begin{split} \frac{\partial \phi_{cyl,C}(\theta)}{\partial t} &=-\frac{\Gamma}{\pi} \frac{1}{1+\frac{z+b}{z-b} \frac{1-\cos\theta}{1+\cos\theta}} \sqrt{\frac{1-\cos\theta}{1+\cos \theta}} \frac{((z-b)U-(z+b)U)}{2 \sqrt{\frac{z+b}{z-b}}}(z-b)^2 \\ & =\frac{\Gamma U b}{\pi} \frac{\sqrt{(1+\cos\theta)(1-\cos\theta)} }{(z-b)(1+\cos\theta)+(z+b)(1-\cos\theta)} \frac{1}{\sqrt{z^2-b^2}} \\ &=\frac{\Gamma U b}{2 \pi} \frac{\sin \theta}{(z-b \cos \theta) \sqrt{z^2-b^2}} \\ & \end{split} \end{equation*}\]
\[\begin{equation*} \begin{split} \Rightarrow L_C &=2\rho\left[\int_{-b}^{b} \frac{\partial \phi}{\partial t} dx +U_\infty \phi_{TE}\right] \\ &=2\rho\left[\int_{\pi}^{0} \frac{\Gamma U b}{2 \pi} \frac{\sin \varphi}{(z-b \cos \varphi) \sqrt{z^2-b^2}} (-2b\sin \varphi d\varphi) +U_\infty \frac{\Gamma}{2}\right] \\ &=\rho U \Gamma \left[1+\frac{b^2}{2 \pi \sqrt{z^2-b^2}} \int_0^{\pi} \frac{\sin^2 \varphi}{z-b\cos \varphi} d \varphi \right] \\ &=\rho U \Gamma \frac{z}{ \sqrt{z^2-b^2}} \\ \end{split} \end{equation*}\]

Aside

The above result suggests that the influence of the starting vortex, on the lift generated by an airfoil in case of a steady angle of attack, reduces as the vortex moves downstream with the flow. In fact as \(t\rightarrow \infty\) so does \(z\rightarrow \infty\), making the total lift generated by the airfoil reduce to the familiar Kutta-Joukowski expression. However, when the starting vortex is close to the airfoil trailing-edge the above result suggests that the amount of lift generated can get very large and even \(L_C \rightarrow \infty\) when \(z\rightarrow b\). One would presume that the situation modelled by the above expression relates closely to the step increase of angle of attack where a single discrete starting vortex might be formed in the flow corresponding to an equivalent circulation established around the airfoil. This is something that is of course not what is observed in reality as has been mentioned earlier.

The disconnect comes from the fact that \(\left. \Gamma_{bound} \right |_{t=0+} \neq \left. \Gamma_{bound} \right |_{t\rightarrow \infty}\). In fact \(\Gamma_{bound}\) develops progressively around the airfoil as time progresses and this is sometimes referred to as the ‘Wagner effect’ [San03] after Herbert Wagner who was the first to propose analytical expressions for the development of this unsteady lift on an impulsively pitched airfoil [Wag25] as shown above.

Airfoil lift developement after step-input [source]

The preceding expressions were for the case of an isolated vortex of strength \(\Gamma\) in the wake of the airfoil. If there exists a continuous distribution of vorticity in the wake \(\gamma_w\) then the circulatory lift is obtained by accounting for the cumulative effect of this vorticity distribution.

\[L_C =\rho U \int_{b}^{\infty}\frac{\gamma_w z}{ \sqrt{z^2-b^2}} dz \]

Now we are at an impasse because the solution of \(\gamma_w\) using the integral equation from earlier is not possible so the net circulatory lift cannot be obtained. At this point in the derivation, Theodorsen introduced the following-

\[L_C =L_{QS} \frac{\int_{b}^{\infty}\frac{\gamma_w z}{ \sqrt{z^2-b^2}} dz}{\int_{b}^{\infty} \gamma_w\sqrt{ \frac{z+b}{z-b}}dz} \]

where,

\[L_{QS}= \rho U^2 b 2\pi\alpha_{3c/4}(t)=\rho U\int_{b}^{\infty} \gamma_w(z,t)\sqrt{ \frac{z+b}{z-b}}dz \]

So, the circulatory component of lift (i.e. lift due to circulation around the airfoil) is equal to the quasi-steady lift times a factor that is a function of the wake. Since \(\gamma_w\) is unknown, the integrals cannot be solved for general airfoil motion. Therefore, while the non-circulatory component of lift acting on an airfoil can be obtained for any random motion, the circulatory component can only be obtained for specific types of motion for which the integrals involving \(\gamma_w\) CAN be solved. Theodorsen solved them assuming simple harmonic motion of the airfoil i.e. \(\alpha_{3c/4}=\alpha_0 e^{i\omega t}\). The crucial bit of realisation necessary for a solution comes from the equation for \(L_{QS}\) (or the integral circulation equation mentioned earlier) which is linear, a not-so-easily solvable integral equation but linear nonetheless, and the property of linear systems that respond harmonically to harmonic excitation with certain amplitude and a phase shift. With that in mind, the ‘response’ of our system \(\gamma_w\) would be of the form

\[\gamma_w(z,t) = \gamma_w(z)e^{i\omega t +\phi} = \overline{\gamma}_w(z)e^{i\omega t}\]

Not just that, because any shed vortices in the wake move with the flow speed \(U\), the following relation also exists

\[\gamma_w(z,t) = \gamma_w\left ( b,t-\frac{z-b}{U} \right )= \overline{\gamma}_w(b)e^{i\omega \left (t-\frac{z-b}{U}\right )} = \overline{\gamma}_w(b)e^{i(\omega t-k(z-b))} = \overline{\overline{\gamma}}_w(b)e^{i(\omega t-kz)}\]

where, \(k= \frac{\omega b}{U}\) is called the reduced frequency and is an important parameter that decides the degree of unsteadiness in the flow; for e.g. if the airfoil is oscillating such that \(k<0.03\) then the unsteady effects are so small that they can be ignored and quasi-steady results suffice. Of course, the corresponding ‘unsteadiness’ is actually coming from \(\omega\) or harmonic osciallation of \(\alpha_{3c/4}\) but \(k\) tells us that the speed with which the vortices in the wake will convect also relevant. The faster the vortices convect the lower is the unsteadiness which is physically intuitive in the sense the \(\omega\) might be large but if the vortices shed quickly move away from the airfoil would resemble the steady state scenario more closely. The circulatory lift expression leads to the following simplification

\[\begin{equation*} \begin{split} L_C &=L_{QS} \frac{\int_b^{\infty} \frac{z}{\sqrt{z^2-b^2}} \gamma_w(z, t) d z}{\int_b^{\infty} \sqrt{\frac{z+b}{z-b}} \gamma_w(z, t) d z} \\ & =L_{QS} \frac{\bar{\gamma}_w e^{i \omega t} b \int_1^{\infty} \frac{z^*}{\sqrt{z^{* 2}-1}} e^{-i k z^*} d z^*}{\bar{\gamma}_w e^{i \omega t} b \int_1^{\infty} \sqrt{\frac{z^*+1}{z^*-1}} e^{-i k z^*} d z^*} \quad \textrm{ where } z^*=\frac{z}{b} \\ & =L_{QS} \frac{\int_1^{\infty} \frac{z^*}{\sqrt{z^*-1}} e^{-i k z^*} d z^*}{\int_1^{\infty} \sqrt{\frac{z^*+1}{z^*-1}} e^{-i k z^*} d z^*} \\ &=L_{QS} C(k) \end{split} \end{equation*}\]

where, the assumption of harmonic motion reduces the integral into a relatively intimidating-looking form referred to as the Theodorsen function \(C(k)\). However, this integral is actually solvable! Without going into details, \(C(k)\) can be written as a function of Hankel functions of the second kind \(H_n^{(2)}\), which in turn can be written using Bessel functions of the first and second kind \(J_n\) and \(Y_n\), respectively. These are analytical functions that often occur in engineering problems, so much so that most standard programming languages have one line commands implementing them.

\[\begin{equation*} \begin{split} C(k) &=\frac{\int_1^{\infty} \frac{z^*}{\sqrt{z^*-1}} e^{-i k z^*} d z^*}{\int_1^{\infty} \sqrt{\frac{z^*+1}{z^*-1}} e^{-i k z^*} d z^*} \\ &= F(k)+i G(k) \\ &=\frac{H_1^{(2)}(k)}{H_1^{(2)}(k)+i H_0^{(2)}(k)} \\ &=\frac{J_1(k) -i Y_1(k)}{J_1(k) -i Y_1(k) + i(J_0(k) -i Y_0(k))} \\ \end{split} \end{equation*}\]

Which represents as incredible achievement because \(C(k)\) is a complex function of the reduced frequency \(k\). If the frequency of harmonic pitching of the airfoil, say, and the constant oncoming flow velocity \(U\) is known using the half-chord value \(b\), \(C(\frac{\omega b}{U})\) can be evaluated in a straightforward manner. Once the quasi-steady three quarter-chord angle of attack of the airfoil is known then \(L_{QS}\) can be evaluated and thereafter the circulatory lift. Remember, all of this can only be carried out only for harmonic motion of the airfoil!

Final Results

Finally, the complete expression for unsteady lift on an airfoil can be written as a sum of the circulatory and the non-circulatory expressions derived above

\[\begin{equation*} \begin{split} L &= L_{NC}+L_C \\ &= L_{NC}+L_{QS}C(k) \\ &= \pi \rho b^2 \left[(U \dot{\alpha}+\ddot{h}-\ddot{\alpha} a b) \right] + \rho U^2 b 2 \pi \alpha_{3c/4} C(k)\\ &= \pi \rho b^2 \left[(U \dot{\alpha}+\ddot{h}-\ddot{\alpha} a b) \right] + 2 \rho U^2 b \pi C(k) \left [\alpha +\frac{\dot{h}}{V}+b \left ( \frac{1}{2}-a \right ) \frac{\dot{\alpha}}{V} \right] \end{split} \end{equation*}\]

Looking at the behavior of \(C(k)\) as a function of \(k\) a clear trend emerges - unsteady circulatory lift is ALWAYS less than its quasi-steady counterpart! Note that the total lift can be higher or lower based on the contribution from the inertial forces (non-circulatory term).

Theodorsen function $C(k)$ as a function of $k$ [source]

Aside

An alternate form of looking at \(C(k)\) is that it is a transfer function of the linear fluid system, since fluid has been modelled using potential flow theory which is based on the linear Laplace’s equation, such that when the input to the system is the motion the final output is a response at the same frequency but with a phase delay and change in amplitude.

Theodorsen function as a response transfer function

For an exemplary case of an airfoil of chord \(c\) pitching as \(\alpha = \alpha_0 \sin (\omega t)\) in an oncoming flow with freestream velocity of \(U_\infty\), the corresponding unsteady lift then is given by

\[L =\pi \rho b^2 \alpha_0 \left[(U\omega \cos (\omega t)+\omega^2 \cos (\omega t) a b) \right] + 2 \rho U^2 b \pi \left | C(k) \right | \left [\sin (\omega t+\phi) + b \left ( \frac{1}{2}-a \right ) \frac{\omega \cos (\omega t+\phi)}{V} \right]\]

where \(\phi\) is the phase of the complex \(C(k)\) at \(k=\omega b/U_\infty\).

Summary and Final Thoughts

In summary, unsteadiness comes from accounting for the vorticity in the wake and when wake vorticity (as in the case of the starting vortex) is not acocunted for the flow is considered to be steady. In addition, a steady airfoil generates lift which is purely due to the circulation of the flow but an airfoil undergoing unsteady motion will generate an added lift that is due to the displacement of the surrounding fluid.

The Theodorsen model is a solution of the unsteady potential flow so the problem of fluid flow represents a linear problem. Potential flow theory is based on the inviscid assumption so clearly the predicted drag force (i.e. force along the direction of flow) is zero leading to the d’Alembert’s paradox. Of course this is not much of a paradox now given we modelled the fluid as an ideal fluid is some sense by ignoring viscosity altogether. Therefore, the Theodorsen model does not make any prediction about the nature of drag.

From the discussion till now, it is clear that there is a significant difference between the wake of an airfoil when it maintains a constant angle of attack versus when the airfoil is pitching. The difference arises due to the shed vorticity in the airfoil wake. In case of steady airfoil, the starting vortex convects far away but in case of the pitching airfoil vorticity is constantly being shed due to the changing bound circulation around the airfoil. These shed vortices in the vicinity of the airfoil have the potential to influence the aerodynamics of the airfoil and need to be accounted for. The criterion used to decide the threshold when unsteady aerodynamic effects might be significant enough, that they cannot be ignored, is called reduced frequency.

Finally, in general, in order to evaluate the aerodynamic load (steady or unsteady) acting on rotating wings, the set of non-linear PDE called NSE need to be solved. This can be a daunting excercise using CFD methods that is both computationally intensive, i.e. requiring significant computing infrastructure, as well as time consuming, i.e. it would be a few days before you would see the final result! The lifting line theory presents a great simplification in that what happens to the flow in the rotor wake can be modelled separate from the effect of flow on the blade itself. The effect of the rotor wake is modelled as an inflow that can be obtained via a hierarchy of inflow models of varying levels of fidelity but we discussed the simplest ones in this course - the momentum theory and the blade element momentum theory. The blade element momentum theory requires prior knowledge of the airfoil characterisicts (i.e. \(C_l\), \(C_d\) and \(C_m\)) but presents a vast improvement over the simple momentum theory in that the effect of using different airfoils making up the blade sections can be accounted for. These airfoil characteristics are based on steady state airfoil aerodynamics. This solution strategy has proven to provide acceptable results in hover and very low speed forward flight but not for high speed forward flight cases.

When unsteady motion of the airfoil is involved, like in forward flight conditions, even an airfoil has a wake whose effect on the airfoil section itself needs to be accounted for in order to evaluate lift, drag and moment on the airfoil accurately. So if we consider the general case of a rotor where the blades are undergoing unsteady motion then one component of the rotor wake vorticity is due to the finiteness of the wings and a secondary component now exists due to the unsteady motion of the blade sections that make up the rotor. Instead of solving for the influence of this combined complex wake on the blade sections, i.e. evaluating the induced inflow, we use the same strategy as used in the case of steady rotor aerodynamics. We can do this be evaluating the effect of the unsteady airfoil wake on the airfoil loads, this is where the Theodorsen theory comes in handy when harmonic motion is involved. Once the effect of the airfoil wake has been accounted for, the effect of the rest of the wake (due to finite wing effects) can be treated using blade element momentum theory. This is the standard process used for modelling rotor aerodynamics using ‘non-CFD’ methods. There exist a number of unsteady airfoil aerodynamics analysis models (particularly in the field of helicopter aerodynamics) one can choose from - Peters model, indicial response methods, state-space modelling etc. - for quantifying the airfoil wake effects on blade section loads and once they are known the rotor inflow can be obtained using BEMT. Here too a number of models exist in addition to the annular momentum theory - for e.g. dynamic inflow, free-wake, lifting surface, vortex particles etc.- but, irrespective of the specific models used, the overarching strategy towards a quick solution follows the standard strategy convexed through this course.

previous

Unsteady Aerodynamics

next

References