26 Time-Dependent Flows


In fluid dynamics and aerodynamics, a time-dependent or unsteady flow refers to the conditions where the properties (such as velocity, pressure, temperature, and density) at a particular point in space will change with time. This behavior contrasts with steady flow, where the fluid properties at a given point remain constant over time. Unsteady flows occur, and indeed are expected, in many real-world situations, such as when a fluid flow is subjected to changes in operating conditions, the initiation of transient events, or when time-varying external events produce non-steady forces on flows or flow systems.

In aerospace engineering, unsteady flow effects occur frequently. For example, they are found and must be considered to be able to predict the fluid behavior in the fields of aeroelasticity and flutter, rotating machinery, including jet engines and helicopter rotors, mixing processes, hydraulic and pneumatic control systems, and combustion processes in air-breathing and rocket engines. Analyzing and predicting unsteady flow can be more challenging than steady flow because it requires careful consideration of all relevant time-dependent factors. Closed-form solutions to most unsteady flow problems are generally only realizable in cases with specific assumptions. Therefore, experimental techniques and computational fluid dynamics (CFD) are often used to study unsteady flows in various engineering applications.

Learning Objectives

  • Understand how to differentiate between steady and unsteady flows and recognize situations in which unsteady flow properties are relevant to problem-solving.
  • Know about reduced frequency and Strouhal number and why they are essential for characterizing unsteady flows and the degree of unsteadiness.
  • Be able to apply the conservation principles of fluid dynamics to solving some exemplar unsteady flow problems.

Classification of Time-Dependent Flows

Steady flows are flows with properties that do not change as a function of time. A steady-state flow refers to the condition where the macroscopic flow properties, such as the velocity and pressure at a point, do not change with respect to time, as shown in the figure below on the left. Mathematically, for steady flows, then

(1)   \begin{equation*} \frac{\partial P}{\partial t} \equiv 0 \end{equation*}

where P is any property such as pressure, temperature, velocity, density, etc. However, in a time-dependent flow, also known as an unsteady flow, the flow properties at a point will change in time, as shown in the figure below on the right. In this case, all of the unsteady terms in any equations used to describe the flow must be retained and solved, i.e.,

(2)   \begin{equation*} \frac{\partial P}{\partial t} \ne 0 \end{equation*}


Examples of the laminar flow at a point in space show the difference between a steady and an unsteady flow.

Unsteady flow phenomena are encountered in many engineering applications. Examples include the flows in turbomachinery and other internal combustion engines, helicopter aerodynamics, wind turbines, and numerous problems in aeroacoustics where creating time-varying aerodynamic loads produces sound (noise). A turbulent flow is an unsteady flow, by definition. However, a turbulent flow can still be statistically steady. This definition means that the average flow velocity and other quantities are constant with respect to time, and all the statistically varying properties, such as the fluctuating velocity component, are constant in time.

The figure below shows the difference between a statistically steady turbulent flow and a statistically unsteady flow. Turbulence enhances the mixing of fluids, promoting the transport of momentum, heat, and other properties within the fluid. The irregular motion of fluid particles makes it difficult to predict the behavior of the flow over time precisely; this randomness, or mathematically stochastic behavior, is a fundamental characteristic of turbulence. However, statistically a flow property P can be decomposed into a mean or average part, \overline{P}, and a mean fluctuating part, P', i.e.,

(3)   \begin{equation*} P = \overline{P} + P' \end{equation*}

This latter process is called a Reynolds decomposition, which is a helpful concept in modeling turbulent flows that will be explained later in this chapter.


Examples of the turbulent flow at a point showing a statistically steady flow versus a statistically unsteady one.

One reason it is helpful to distinguish between steady and unsteady flows is that the former is often more tractable to understand and predict. Eliminating time from the solution of the equations that govern fluid flow problems usually results in a significant simplification of the governing equations as well as the mathematical and/or numerical techniques needed to solve these equations. Characterizing steady, quasi-steady, or unsteady flows is often done using parameters such as reduced frequency or Strouhal number.

Reduced Frequency

One nondimensional parameter used to help categorize whether a flow is steady or unsteady is the reduced frequency, given the symbol k. The reduced frequency is defined, in general, as

(4)   \begin{equation*} k = \frac{\omega \, L}{V} \end{equation*}

where \omega is a characteristic physical frequency of the unsteady flow (in units of radians per second), L is a characteristic length scale, and V is a reference flow velocity. Base units should always be used to obtain the correct numerical value for the reduced frequency. For a wing or airfoil, such as one oscillating in angle of attack or an oscillatory vertical (heaving) motion as shown in the figure below, the reduced frequency is often defined in terms of its semi-chord, i.e., L = c/2, and the free-stream velocity, V_{\infty}, i.e., in this case, it is defined as

(5)   \begin{equation*} k = \frac{\omega \, c}{2 \, V_{\infty}} \end{equation*}

The aerodynamics of oscillating airfoil sections and wings must be addressed in the field of aeroelasticity and flutter.


Unsteady airfoil motion produces different aerodynamic characteristics than a steady flow; these issues are encountered in aeroelasticity.

For k = 0, the flow is steady, and all the usual aerodynamic results apply regarding the relationships between the aerodynamic quantities and the angle of attack. For 0 \le k \le 0.05, the flow can be considered quasi-steady; that is, unsteady effects are generally minor, and for some problems, they may be neglected completely. However, such reduced frequency bounds may be somewhat subjective.

Flows with characteristic reduced frequencies of 0.05 and above are usually considered unsteady, so all unsteady terms must be retained in the governing flow equations. Such problems are more challenging to determine for two reasons. First, the local flow properties also depend on the previous time, i.e., what has happened in the time history of the lift and other aerodynamic forces. Second, the effects of compressibility may be significant even if the free-stream flow velocity is low. This means that the time scales of the unsteady motion become of the order of propagating pressure disturbances in the flow, which occur at the speed of sound.

Strouhal Number

The Strouhal number, given the symbol St, is another dimensionless number used in fluid dynamics to characterize the behavior of oscillating or vibrating objects within a fluid flow. It is beneficial for analyzing phenomena such as vortex shedding and sound generation. The Strouhal number is defined as the ratio of the oscillation frequency to the product of the object’s characteristic length (or diameter) of the object and the reference or free-stream velocity of the fluid, e.g., V_{\infty}. The equation is given by

(6)   \begin{equation*} St = \frac{f \, L}{V_{\infty}} \end{equation*}

Notice that the physical frequency {\scriptstyle{f}} (in units of per second or Hz) is used to determine the Strouhal number rather than the circular frequency in the reduced frequency. The Strouhal number is named after the Czech physicist and engineer Vincenc Strouhal, who significantly contributed to the study of oscillating flows in the late 19th and early 20th centuries.

The Strouhal number is often used to analyze various fluid flow phenomena, such as the shedding of vortices behind a cylinder or another type of bluff (non-streamlined) body. It helps characterize the vortex shedding frequency relative to the fluid flow velocity and the object’s size. For different types of flows and body shapes, there are often characteristic Strouhal number ranges corresponding to specific flow regimes and behaviors, i.e., the behavior locks into a particular shedding frequency.


The frequency of vortex shedding from a bluff body, such as a circular cylinder, is often characterized by the Strouhal number.

Check Your Understanding #1 – Calculating the value of the Strouhal number

A circular cylinder with a diameter d of 0.1 m is immersed in a fluid flow with a velocity V of 1.0 m/s. The shedding frequency {\scriptstyle{f}} of vortex shedding behind the cylinder is measured to be 2 Hz. What is the Strouhal number in this case?

Show solution/hide solution.

The Strouhal number is defined by

    \[ St = \frac{f \, L}{V} \]

Using the numerical values gives

    \[ St = \frac{ 2.0 \times 0.1}{1.0} = 0.2 \]

This value of 0.2 is often associated with the onset of vortex shedding behind a cylindrical object in a fluid flow. For example, suppose a particular combination of flow speeds and length scales results in a Strouhal number of 0.2. In that case, vortex shedding can be expected, and vortex-induced vibrations on a structure are possible. However, it is essential to note that specific Strouhal number values may vary for different flow conditions and other body shapes.

Time-Dependent Fluid Flows

Time-dependent flows are often more challenging to deal with in fluid dynamics, but the same principles of conservation of mass, momentum, and energy will apply. Studying exemplars is an excellent way to learn how to solve unsteady-flow problems. In this regard, several classic exemplars in fluid dynamics can be used to lay down the solution principles when considering the additional dimension of time.

Liquid Draining from a Tank

Consider the problem of a tank of circular cross-section, A, that lets a liquid leave through a drain valve at the bottom, as shown below. The discharge velocity, {V_d}, of the flow from the drain varies with the height, h, of the fluid level above the drain according to the relationship V_d = \sqrt{ 2 g h}. The fluid is a liquid, so \varrho = constant. Notice that the height, h, will decrease with time as the liquid leaves the tank.

Illustration of a time-dependent flow problem where a fluid continuously discharges from a drain.

The general form of the continuity equation is

(7)   \begin{equation*} \frac{\partial}{\partial t}\oiiint \varrho \, d{\cal{V}} + \oiint_S \varrho \, \vec{V} \bigcdot d\vec{S} = 0 \end{equation*}

In this case, \varrho = constant so that it can be reduced to

(8)   \begin{equation*} -\frac{d{\cal{V}}}{dt} +  a \, V_d  = 0 \end{equation*}

where a is the area of discharge of the liquid from the drain. The volume of liquid \cal{V} in the tank for any height h is

(9)   \begin{equation*} {\cal{V}} = A h \end{equation*}

where A is the cross-sectional area of the tank.

It is given that

(10)   \begin{equation*} V_d = \sqrt{ 2 g h} \end{equation*}

and this flow rate depends on the instantaneous height of the liquid in the tank, h. This latter result is called Torricelli’s theorem, and it applies to a liquid that flows out of an orifice from the effects of a hydrostatic pressure head. Therefore, in a short time interval dt, then the decrease in the volume \cal{V} of liquid in the tank will be

(11)   \begin{equation*} \frac{d \cal{V}}{dt} = - A \frac{dh}{dt} \end{equation*}

where the minus sign denotes that h decreases as the tank empties. Notice that the units in the preceding equation are in terms of a volume flow rate.

By continuity considerations (conservation of mass), this flow rate must be equal to the flow rate of liquid coming out of the drain, i.e., the problem is time-dependent, so

(12)   \begin{equation*} \frac{d \cal{V}}{dt} = Q \end{equation*}

where Q is the flow rate. Therefore, using the conservation of mass gives

(13)   \begin{equation*} - A \frac{dh}{dt} = a V_d = a \sqrt{ 2 g h} \end{equation*}

Rearranging this latter equation gives

(14)   \begin{equation*} \frac{dh}{dt} = - \frac{a}{A} \sqrt{ 2 g h} \end{equation*}

which is an ordinary differential equation relating the height of the level of the liquid to time.

Separating the variables and integrating gives

(15)   \begin{equation*} t = \int_{0}^{t} dt = - \frac{A}{a} \int_{h_0}^0 \frac{dh}{\sqrt{ 2 g h}} = - \frac{A}{a \sqrt{ 2 g}} \int_{h_0}^0 \frac{dh}{\sqrt{h}} = \frac{A}{a \sqrt{ 2 g}} \int_0^{h_0} \frac{dh}{\sqrt{h}} \end{equation*}

where the limits of integration are: At time t = 0, then h = h_0, and when h = 0, the tank is considered empty. Performing the integration gives

(16)   \begin{equation*} t = \frac{2A}{a \sqrt{ 2 g}} \sqrt{h_0} \end{equation*}

In terms of the diameter of the tank, D, and the outlet diameter, d, of the drain, then

(17)   \begin{equation*} \frac{A}{a} = \frac{D^2}{d^2} \end{equation*}

Therefore, the time required to empty the tank will be

(18)   \begin{equation*} t = \frac{2D^2}{d^2} \sqrt{\frac{h_0}{2 g}} \end{equation*}

where h_0 is the initial height of the liquid in the tank when the drain is first opened.

Filling an Air Tank

Consider now a rigid tank of volume {\cal{V}} with air pumped into it at a constant mass flow rate, as shown in the figure below. In this problem, the effect of compressibility must be considered. Still, it can be assumed that the process is slow enough that it is isothermal, and the heat generated by compressing the air is dissipated through radiation by the tank.

Filling an air tank is an unsteady flow process that also involves the compressibility of the air.

This is also an unsteady flow problem because a mass of air is pumped into a fixed volume, and for mass conservation, the air density must increase in time. The initial density and pressure are \varrho_0 and p_0, respectively. The general form of the continuity equation is

(19)   \begin{equation*} \frac{\partial}{\partial t}\oiiint \varrho \ d {\cal{V}} + \oiint_S \varrho \, \vec{V} \bigcdot d\vec{S} = 0 \end{equation*}

so, in this case

(20)   \begin{equation*} \oiint_S \varrho \, \vec{V} \bigcdot d\vec{S} = -\overbigdot{m}_{\rm in} \end{equation*}

and so

(21)   \begin{equation*} \frac{\partial}{\partial t}\oiiint \varrho \, d{\cal{V}} = {\cal{V}}\frac{d \varrho}{dt} \end{equation*}

Assume a uniform mixing of air so the density inside the volume is uniform. Therefore,

(22)   \begin{equation*} {\cal{V}} \frac{d \varrho}{dt} = \overbigdot{m}_{\rm in} \end{equation*}

Integrating using the separation of variables gives

(23)   \begin{equation*} \int_{\varrho_0}^{\varrho} d \varrho = \frac{ \overbigdot{m}_{\rm in}}{{\cal{V}}} \int_{t_0}^{t} dt \end{equation*}

so that

(24)   \begin{equation*} \varrho = \varrho_0 + \frac{ \overbigdot{m}_{\rm in}}{{\cal{V}}} \left( t - t_0 \right) \end{equation*}

which shows that the density of the air increases linearly with time.

The corresponding pressure can be obtained from the equation of state, i.e., p = \varrho \, R  \, T. If the process is assumed to be isothermal, then T = T_0, and so

(25)   \begin{equation*} p = \varrho_0 \,  R \,  T_0 + \frac{ R \,  T_0 \, \overbigdot{m}_{\rm in}}{{\cal{V}}} \left( t - t_0 \right) = p_0 + \frac{R \,  T_0 \, \overbigdot{m}_{\rm in}}{{\cal{V}}} \left( t - t_0 \right) \end{equation*}

Hydraulic Shock

Hydraulic shock, also known as “water hammer” or “hydraulic hammer,” is a time-dependent flow phenomenon that occurs in fluid flow systems, such as pipelines, when there is a sudden change in the flow velocity. This abrupt change can result in strong pressure waves propagating through the fluid, causing transitory surges in pressure inside the pipeline, as shown in the schematic below. Water hammer is a common occurrence in fluid flow systems where valves and regulators are rapidly opened or closed; the rapid, time-dependent changes in pressure can generate loud banging or knocking sounds, often audible throughout the piping system.

Hydraulic shock or “water hammer” involves the generation of large pressure pulses in a flowing fluid after the rapid closure of a valve.

The fundamental cause of water hammer is the fluid’s inertia and the need for energy conservation when a sudden stoppage of the flow occurs. When a fluid is in motion and experiences a sudden change in velocity, its kinetic energy must be either absorbed or dissipated. If a valve is closed suddenly, the fluid comes to an abrupt stop, all kinetic energy is converted into pressure (or the equivalent in terms of pressure), and the flow stagnates in the pipe.

This process creates a pressure wave that travels back upstream through the pipe; this pressure wave travels through the fluid with speed, a, known as the wave speed. When the pressure wave reaches the reservoir, the higher pressure head starts the flow downstream again, with continuing oscillations until the energy of the pressure wave is dissipated through frictional effects.

The resulting pressure pulses can be several times higher than the normal operating pressure. The change in pressure from the water hammer depends on how quickly the value is closed, the type of fluid, the length of the pipe, its elastic properties, and how the pipe is mounted, e.g., supported at its ends and/or clamped along its length. Prolonged or severe water hammer can lead to structural damage to the pipes, their fittings, and other system components. The phenomenon can cause fatigue, leaks, and pipe bursts and accelerate the wear and tear on valves and pumps, increasing maintenance requirements. Preventing or mitigating water hammer effects involves engineering solutions and using various devices, such as surge tanks, air chambers, and water hammer arrestors, which are designed to absorb the excess pressure and prevent it from causing damage to the system.


Unsteady effects, called water hammer, are produced by closing a valve in a pipeline.

Consider the analysis of this problem. A liquid stored in a tank flows steadily through a pipe of length L, as shown in the figure below. At the time t = 0, the valve at the downstream end is quickly closed, producing the classic pressure pulse of a water hammer. Water hammer transients are typically axisymmetric because the axial mass, momentum, and energy changes are significantly more significant than their radial counterparts.

Nikolay Joukowsky laid the foundation of the water hammer theory[1] where the pressure amplitude, \Delta p, in the liquid in the pipe is related to the change in flow velocity, \Delta V, using

(26)   \begin{equation*} \Delta p = \pm \, \varrho \, a \, \Delta V \end{equation*}

where c is the wave propagation speed (assumed to be constant) and \varrho is the density of the liquid; this equation is basically what is referred to today as piston theory. The negative sign in Eq. 26 describes a water hammer wave moving downstream, while the positive sign describes the wave moving upstream. Because pressure head is often used in the field of hydraulics, then Eq. 26 can also be written as

(27)   \begin{equation*} \Delta h = \frac{\Delta p}{\varrho \,  g} = \pm \left( \frac{ a \, \Delta V}{g} \right) \end{equation*}

There are three cases of interest when it comes to predicting water hammer effects, all of which involve the wave propagation speed:

  1. Gradual closure of the valve.
  2. Sudden closure of the valve.
  3. Closure of the valve allowing for the elasticity of the pipe, which will affect the wave propagation speed.

The transit time, also known as reflection time, taken for the pressure wave to propagate from the valve to the tank and then back to the valve will be

(28)   \begin{equation*} t_t = \frac{2 L}{a} \end{equation*}

If the time it takes to close the valve is t_c, then if t_c > 2L/a is longer than the reflection time, the valve closure is considered gradual. If t_c < 2L/a is shorter than time, the reflection of the valve closure is referred to as being sudden.

If t_c > 2L/a, which is a gradual closure of the valve, then the increase in pressure from the water hammer is given by

(29)   \begin{equation*} \Delta p = \frac{\varrho \, L \, V_0}{t_c} \end{equation*}

V_0 is the average flow velocity inside the pipe before the valve closes. The equivalent pressure head of the water hammer is

(30)   \begin{equation*} \Delta h = \frac{\Delta p}{\varrho \,  g} = \frac{\varrho \, L \, V_0}{\varrho \, g \, t_c} = \frac{L \, V_0}{g t} \end{equation*}

If t_c < 2L/c, which is a sudden closure of the valve, then the increase in pressure from the water hammer is given by

(31)   \begin{equation*} \Delta p = \varrho \, a \, V_0 \end{equation*}

and the equivalent pressure head is

(32)   \begin{equation*} \Delta h = \frac{a  \, V_0}{g} \end{equation*}

Notice that in both cases, the propagation speed of the pressure wave, a, is needed. If the pipe is rigid, then

(33)   \begin{equation*} a = \sqrt{ \frac{E_b}{\varrho} } \end{equation*}

where E_b is called the bulk modulus of the liquid, for which values for various liquids are available. This latter equation is often referred to as the Newton-Laplace equation.

If the pipe is elastic, for which most will be to a lesser or greater degree, then

(34)   \begin{equation*} a = \sqrt{ \frac{E_c}{\varrho} } \end{equation*}

where E_c is called the effective bulk modulus of the liquid in the pipe. The value of E_c can be obtained using

(35)   \begin{equation*} \frac{1}{E_c} = \frac{1}{E_b} + \frac{D \,  K}{E_p \, e} \end{equation*}

where E_p is the modulus of elasticity of the pipe, D is the pipe diameter, and e is the wall thickness of the pipe. The value of K depends on exactly how the pipe is mounted and anchored; usually K \approx 1.0. The modulus of elasticity of the pipe is a material property for which values are also widely available.

Finally, consider the valve’s sudden closure when the pipe’s elasticity is considered. In this case

(36)   \begin{equation*} \Delta p = \varrho \,  a \, V_0 = \varrho \sqrt{ \frac{E_c}{\varrho} } \, V_0 = V_0 \sqrt{ \frac{\varrho}{ \displaystyle{ \left( \frac{1}{E_b} + \frac{D \, K}{E_p \, e} \right) }} } \end{equation*}

Check Your Understanding #2 – Determining hydraulic shock pressure in a pipeline

A steel pipe that is 1,000 ft long has a 15-inch diameter and a 1.5-inch wall thickness. The pipe carries water from a reservoir, and a valve is located downstream. When the valve is fully open, the flow rate through the pipe is 23.3 ft{^3}/s. If the valve is closed completely in 0.15 seconds, determine the water hammer pressure in the pipe. Assume E_b = 3.0 \times 10^5 lb/in^2, and E_p = 2.8 \times 10^7 lb/in^2, and K = 1.0. The density of water is 1.94 slugs/ft{^{3}}.

Show solution/hide solution.

The relevant equation, in this case, for the water hammer pressure in an elastic pipe is

    \[ \Delta p = \varrho \,  a \, V_0 = \varrho \sqrt{ \frac{E_c}{\varrho} } V_0 \]

where E_c is the effective bulk modulus of the water in the pipe. The average velocity through the pipe, V_0, is related to flow rate, Q, using

    \[ V_0 = \frac{Q}{A} = \frac{4 Q}{\pi \,  D^2} = \frac{4 \times 23.3}{\pi (15.0/12.0)^2}  = 18.99~\mbox{ft/s} \]

The effective bulk modulus is given by

    \[ \frac{1}{E_c} = \frac{1}{E_b} + \frac{D \,  K}{E_p \, e} = \frac{1}{3.0 \times 10^5} + \frac{(15.0/12.0) \times 1.0}{2.8 \times 10^7 \times (2.0/12.)} = 3.75 \times 10^{-6} \]

so that the value of E_c is 2.67 \times 10^5 lb/in{^2}. Therefore, the wave propagation speed is

    \[ a = \sqrt{ \frac{E_c}{\varrho} } = \sqrt{ \frac{144 \times 2.67 \times 10^5}{1.94} } = 4,451.8~\mbox{ft/s} \]

and the wave reflection time is

    \[ t_t = \frac{2 L}{a} = \frac{2 \times 1,000}{4,451.8} = 0.45~\mbox{s} \]

And so with a valve closure time of 0.15 s, then t_c < t_c, this is a sudden closure. Therefore, the water hammer change in pressure will be

    \[ \Delta p = \varrho \,  V_0 \,  a = 1.94 \times 18.99 \times 4,451.8 = 1.64 \times 10^5~\mbox{lb/ft$^2$} = 1,138.94~\mbox{lb/in$^2$} \]

Why is the speed of sound so much higher in a liquid than in air?

The speed of sound in a medium, whether it’s a gas, liquid, or solid, depends on the properties of that medium. Sound is faster in liquids and solids than in gases because the molecules in liquids and solids are much closer together, allowing for more rapid transmission of mechanical waves. The speed of sound in a medium is determined by its density and elastic properties, specifically its bulk modulus. The bulk modulus measures a substance’s resistance to compression or expansion when subjected to pressure. Liquids and solids generally have higher bulk moduli than gases, making them less compressible. In a liquid, the molecules are closer together than in a gas, and they can still transmit mechanical waves. This closeness of molecules and the resistance to compression contribute to a higher sound speed in liquids than in gases.

 Producing Jet Thrust

A monopropellant thruster is a very basic rocket, often used for satellite propulsion and control, with a certain mass of propellant being stored in a pressurized storage tank. The propellant is released over time by opening the valve to flow over a catalyst bed, causing thermal decomposition and the liberation of energy from the propellant, as shown in the figure below. The flow then expands through a nozzle to produce an exit velocity and, hence, a thrust from the time rate of change of momentum of the expanding gases. The idea is that the thrust boosts the satellite to change its orbital velocity and altitude.

This is a time-varying flow problem because a propellant mass is continuously discharged from a tank, so the mass of the propellant inside the control volume is decreasing. The most general form of the continuity equation is

(37)   \begin{equation*} \frac{\partial}{\partial t}\oiiint \varrho \, d {\cal{V}} + \oiint_S \varrho \, \vec{V} \bigcdot d\vec{S} = 0 \end{equation*}

and in this case, it can be reduced to

(38)   \begin{equation*} - \frac{d M_p}{d t} + \overbigdot{m}_p = 0 \end{equation*}

where M_p is the mass of the propellant and \overbigdot{m}_p is the propellant flow rate. Notice that the minus sign occurs because the total propellant mass is decreasing.

Let the initial mass of the spacecraft be M_0. Conservation of mass gives the current mass of the spacecraft, M, at time t later as the propellant is discharged as

(39)   \begin{equation*} M  = M_0 - \left( \frac{dM_p}{dt} \right)  t = M_0 -\overbigdot{m}_p \, t \end{equation*}

where \overbigdot{m}_p = constant. The thrust force from the propulsion system is given by Newton’s second law (force equals time rate of change of momentum), i.e.,

(40)   \begin{equation*} T = \overbigdot{m}_p V_e \end{equation*}

and so the acceleration of the spacecraft, which is not constant because of its continuously decreasing mass, is given by

(41)   \begin{equation*} a = \frac{T}{M} = \frac{\overbigdot{m}_p \, V_e}{M} = \frac{dV}{dt} \end{equation*}


(42)   \begin{equation*} \frac{dV}{dt} = \frac{\overbigdot{m}_p \, V_e}{M_0 -\overbigdot{m}_p \, t} \end{equation*}

Integrating the equation gives

(43)   \begin{equation*} \Delta V = V_e \ln \left( \frac{M_0}{M_0 -\overbigdot{m}_p \, t} \right) \end{equation*}

which has a special name called the rocket equation.

Check Your Understanding #3 – Discharge of a monopropellant thruster

An orbiting satellite is propelled using a monopropellant thruster. The satellite has an initial mass of 5,000 kg, including its propellant mass. To give a slight corrective boost to its orbit, it opens the control valve and ejects the propellant at a constant rate of 0.1 kg/s at an effective exit velocity of 2,000 m/s. Assume that the satellite operates in a vacuum and no gravity forces are acting on it. Determine the thrust, the change in velocity of the satellite, and the new mass after 90 seconds from the start of the propulsive burn.

Show solution/hide solution.

The thrust produced can be determined from the time rate of change of momentum of the exhaust gases, i.e.,

    \[ T = \overbigdot{m}_p \, V_e = 0.1 \times 2,000 = 200~\mbox{N} \]

After 90 seconds of thrusting, i.e., t = 90, then using the rocket equation gives

    \[ \Delta V = 2,000 \ln \left( \frac{5,000}{5,000 - 0.1 \times 90.0 } \right) = 3.6~\mbox{m/s} \]

The final mass after the burn will be

    \[ M  = M_0 -\overbigdot{m}_p \, t = 5,000 - 0.1 \times 90 = 4,991~\mbox{kg} \]

Primer on Unsteady Airfoil Behavior

Unsteady airfoil theory deals with the aerodynamic behavior of airfoils and wings when subjected to time-varying boundary conditions, such as oscillations in angle of attack, various vertical and horizontal gusts, or other transient phenomena. Understanding unsteady aerodynamics is crucial for applications involving aircraft maneuvering, rotorcraft aerodynamics, and wind turbine blades.

In an unsteady flow, the wake behind the airfoil is time-varying in structure with the shedding of vorticity and circulation, thereby influencing the aerodynamic forces experienced by the airfoil. These effects change the magnitude of the lift and the phase of the lift response relative to the time-varying boundary conditions. At high angles of attack, the flow may separate from the airfoil surface, causing a rapid loss of lift and increased drag, a phenomenon called dynamic stall.

There are many ways of calculating these unsteady effects. One approach that is very instructive for those first learning the subject of unsteady airfoil behavior is the classic unsteady aerodynamics theories, including Theordorsen’s theory, which is discussed next. Unsteady panel methods discretize the airfoil surface into panels on the surface and extend into the wake. CFD methods, which are the most computationally intensive, solve the Navier-Stokes equations for unsteady flow, which can capture detailed flow phenomena such as vortex shedding and dynamic stall.

Quasi-Steady Aerodynamics

In quasi-steady flows, the behavior can be evaluated by applying the principles of steady flow under the instantaneous boundary conditions of flow tangency to the airfoil surface.  For example, for a thin airfoil, the changes in the angle of attack of a harmonically oscillating airfoil can be expressed as

(44)   \begin{equation*} \alpha (t) = \alpha_0 + \bar{\alpha} \sin \omega t \end{equation*}

where \alpha_0 is the mean angle of attack and \bar{\alpha} is the angle of attack amplitude of the oscillation. If the flow were quasi-steady, in the sense that flow adjustments were to take place instantaneously, the lift coefficient would be given by

(45)   \begin{equation*} C_l (t) = C_{l_{\alpha}} \left( \alpha + \frac{\overbigdot{\alpha} \, c}{2 \, V_{\infty}} \right) \end{equation*}

recognizing that the pitch rate, i.e., \overbigdot{\alpha}, must affect the lift coefficient because it changes the angle of attack. This equation then leads to

(46)   \begin{equation*} C_l (t) = C_{l_{\alpha}} \bigg( \alpha_0 + \bar{\alpha} \sin \omega t + k \,  \bar{\alpha} \cos \omega t \bigg) = C_{l_{\alpha}} \sqrt{1 + k^2} \bigg( \alpha_0 + \bar{\alpha} \sin (\omega t + \phi ) \bigg) \end{equation*}

where \phi is a frequency-dependent phase angle given by

(47)   \begin{equation*} \phi = \tan^{-1} \big( k \big) \end{equation*}

Therefore, even based on quasi-steady arguments, the lift response will no longer be in phase with the angle of attack (it will lead by the angle \phi in this case) nor have the same amplitude (it is larger than C_{l_{\alpha}} by a factor \sqrt{1 + k^2}). The upshot is that the aerodynamic forces will have different values depending on whether the angle of attack is increasing or decreasing with respect to time.

Indeed, measurements suggest that the unsteady lift on an airfoil or wing is attenuated (diminished) by unsteady effects not amplified as the quasi-steady theory would suggest, and there is a lag (not a lead) in the aerodynamic response with respect to the angle of attack. These effects are formally embodied within the theory of unsteady aerodynamic behavior, for which extensive literature exists. One of these theories is called Theordorsen’s theory.

Theodorsen’s Theory

Theodorsen’s theory[2] provides a linearized, small-disturbance solution for oscillating airfoils in an incompressible flow. The theory is based on the assumption of small oscillation amplitudes. It gives an exact analytic theory (within the assumptions and limitations) to describe the unsteady air loads on a thin airfoil. Theodorsen’s problem was obtaining the solution for the loading (bound circulation) \gamma_b on the airfoil surface under harmonic forcing conditions, the flow model being shown in the figure below.

Flow model used for Theordorsen’s theory. The objective is to find the unsteady loading on a thin airfoil accounting for the effects of circulation shed into the downstream wake.

The governing integral equation is

(48)   \begin{equation*} w(x,t) = \frac{1}{2\pi} \int_0^c \frac{\gamma_b(x,t)}{(x-x_0)} dx + \frac{1}{2\pi}\int_c^{\infty} \frac{\gamma_w(x,t)}{(x-x_0)} dx \end{equation*}

where w is the downwash on the airfoil surface. This equation must be solved subject to invoking the Kutta condition at the trailing edge, i.e., \gamma_b (c,t)=0. There is also a connection to be drawn between the change in circulation about the airfoil and the circulation shed into the wake. Conservation of circulation requires that

(49)   \begin{equation*} \gamma_w(c,t) dx = -\frac{d\Gamma(t)}{dt} dt. \end{equation*}

Assuming that the vorticity in the shed wake is convected at the free stream velocity, V_{\infty}, this gives dx = V_{\infty} dt, and so

(50)   \begin{equation*} V_{\infty} \gamma_w(c,t) = -\frac{d\Gamma(t)}{dt} , \end{equation*}

where the total circulation about the airfoil, \Gamma(t) is given by

(51)   \begin{equation*} \Gamma(t) = \int_0^c \gamma_b(x,t) dx. \end{equation*}

This wake vorticity changes the downwash velocity over the airfoil, and the loads on the airfoil are also affected. Suppose the net circulation about the airfoil is changing with respect to time. In that case, circulation will be continuously shed into the wake and will continuously affect the aerodynamic loads on the airfoil. In the limit, as the changing boundary conditions become constant, the shed wake vorticity cast off the trailing edge of the airfoil becomes zero, and the remaining circulation in the wake convects downstream to infinity.

The unsteady problem posed above is certainly not trivial to solve, but for simple harmonic motion, the solution is given by Theodorsen in a form that represents a transfer function between the forcing (angle of attack) and the aerodynamic response (i.e., the chordwise loading distribution, lift, and pitching moment). The Theodorsen function, denoted by C(k), is a complex-valued function that captures the instantaneous effect of the circulation shed into the wake on the unsteady aerodynamic response and is given by

(52)   \begin{equation*} C(k) = F(k) + i G(k) = \frac{H_1^{(2)}(k)}{H_1^{(2)}(k) + i H_0^{(2)}(k)} \end{equation*}

where H_0^{(2)} and H_1^{(2)} are Hankel functions of the second kind, and k is the reduced frequency. The Hankel function is defined as H_{\rm v}^{(2)} = J_{\rm v} - i Y_{\rm v}, with J_{\rm v} and Y_{\rm v} being Bessel functions of the first and second kind, respectively. The C(k) function is an aerodynamic transfer function that relates the input (forcing) to the response.

Recognizing that each Bessel function has an argument k, then the real or in-phase (\Re) part and imaginary or out-of-phase part (\Im) of the Theodorsen function can be written as

(53)   \begin{equation*} \Re C(k) = F = \frac{J_1 ( J_1 + Y_0 ) + Y_1 (Y_1 - J_0)}{(J_1 + Y_0)^2 + (J_0 - Y_1)^2} \end{equation*}


(54)   \begin{equation*} \Im C(k) = G = -\frac{Y_1 Y_0 + J_1 J_0}{(J_1 + Y_0)^2 + (J_0 - Y_1)^2} \end{equation*}

Despite its apparent mathematical complexity, the Theodorsen function is readily calculated numerically using MATLAB and is plotted in the figure below.

Theodorsen’s function plotted in terms of its real and imaginary parts.

Another way of interpreting the Theodorsen function, which is perhaps more tangible in the first instance, is in terms of the amplitude, \| C(k) \|, and phase angle, \phi, of the response, as also shown in the figure above, which are given by

(55)   \begin{equation*} \| C(k) \| = \sqrt{ F^2 + G^2} \quad \mbox{and} \quad \phi = \tan^{-1} \left( \frac{G}{F} \right) \end{equation*}

MATLAB code to calculate the Theodorsen function

Show the code/hide the code

function theodorsen_function_plot_and_save
% Define the range of reduced frequencies
k = linspace(0.01, 10, 500); % Example range from 0.01 to 10 with 500 points

% Preallocate arrays for real and imaginary parts
Ck_real = zeros(size(k));
Ck_imag = zeros(size(k));

% Calculate the Theodorsen function for each k
for i = 1:length(k)
[Ck_real(i), Ck_imag(i)] = theodorsen_function(k(i));

% Plot the real part
subplot(2, 1, 1);
plot(k, Ck_real, ‘b-‘, ‘LineWidth’, 1.5);
xlabel(‘Reduced Frequency k’);
ylabel(‘Real Part of C(k)’);
title(‘Real Part of Theodorsen Function’);
grid on;

% Plot the imaginary part
subplot(2, 1, 2);
plot(k, Ck_imag, ‘r-‘, ‘LineWidth’, 1.5);
xlabel(‘Reduced Frequency k’);
ylabel(‘Imaginary Part of C(k)’);
title(‘Imaginary Part of Theodorsen Function’);
grid on;

% Save the data to a file
data = [k.’, Ck_real.’, Ck_imag.’];
save(‘theodorsen_function_data.txt’, ‘data’, ‘-ascii’, ‘-double’);

function [Ck_real, Ck_imag] = theodorsen_function(k)
% This function calculates the Theodorsen function C(k)
% and returns its real and imaginary parts.
% k is the reduced frequency

% Bessel functions of the first kind
J0 = besselj(0, k);
J1 = besselj(1, k);

% Bessel functions of the second kind
Y0 = bessely(0, k);
Y1 = bessely(1, k);

% Hankel functions of the second kind
H1 = J1 – 1i * Y1;
H0 = J0 – 1i * Y0;

% Theodorsen function C(k)
Ck = H1 ./ (H1 + 1i * H0);

% Separate into real and imaginary parts
Ck_real = real(Ck);
Ck_imag = imag(Ck);

Theodoresen’s theory is a frequency domain theory; even though the results for the aerodynamic response may be expressed as a function of time, the unsteady forcing must be oscillatory, which is expressed in terms of the reduced frequency. Recall that the reduced frequency is a dimensionless parameter representing the unsteadiness of the flow, defined by

(56)   \begin{equation*} k = \frac{\omega c}{2\, V_\infty} \end{equation*}

Physically, the reduced frequency can be interpreted as the ratio of the oscillation frequency to the convective frequency of the flow. Airfoils respond differently to different values of k. Low k values imply a quasi-steady behavior, while high k values indicate highly unsteady conditions. Quasi-steady means that the flow responds almost immediately to any changes in the angle of attack; therefore, the aerodynamic response will be in phase with the airfoil motion.

The effect of the Theodorsen function, which emulates a physical wake lag response, can be appreciated if a pure oscillatory variation in the angle of attack is considered, as before, i.e., \alpha (t) = \alpha_0 + \bar{\alpha} \sin \omega t,[3] where \alpha_0 is the mean angle of attack and \bar{\alpha} is the angle of attack amplitude of the oscillation.  If the flow were quasi-steady, in the sense that flow adjustments were to take place instantaneously, the lift coefficient would be given by

(57)   \begin{equation*} C_l (t) = C_{l_{\alpha}} \alpha = C_{l_{\alpha}} \left( \alpha_0 + \bar{\alpha} \sin \omega t \right) \end{equation*}

where C_{l_{\alpha}}  = 2 \pi for a thin-airfoil in incompressible flow. In the unsteady case, this circulatory part of the airfoil lift coefficient becomes

(58)   \begin{equation*} C_l(t) = C_{l_{\alpha}} \| C(k) \|  \bigg( \alpha_0 + \bar{\alpha} \sin ( \omega t + \phi ) \bigg) \end{equation*}

or in terms of the real and imaginary parts, then

(59)   \begin{equation*} \Re C_l(t) = C_{l_{\alpha}} \left( \alpha_0 + \bar{\alpha} \sin \omega t \right) \, F(k)\quad \mbox{and} \quad \Im C_l(t) = C_{l_{\alpha}} \left( \alpha_0 + \bar{\alpha} \sin \omega t \right) \, G(k) \end{equation*}

The amplitude and phase lag of the lift response is then

(60)   \begin{equation*} \| C_l(t) \| \sqrt{ \Re C_l^2 + \Im C_l^2 }\quad \mbox{and} \quad \phi = \tan^{-1} \left( \dfrac{\Im C_l}{\Re C_l} \right) \end{equation*}

Representative results from the Theodorsen theory are shown in the figure below in terms of lift coefficient C_l versus angle of attack, \alpha. Notice that the Theodorsen function introduces an amplitude reduction and phase lag effect on the circulatory part of the lift response compared to the result obtained under quasi-steady conditions.

Representative unsteady lift coefficient versus angle of attack on an airfoil when oscillating at different reduced frequencies.

For k = 0, the steady-state lift behavior is obtained, i.e., C_l is linearly proportional to \alpha. However, as k is increased, the lift plots develop into hysteresis loops, and these loops rotate such that the amplitude of the lift response (slope of the lift curve) decreases with increasing values of reduced frequency. Notice also that these loops are circumvented in a counterclockwise direction such that the lift is lower than the steady value when \alpha is increasing with time and higher than the steady value when \alpha is decreasing with time, i.e., there is a phase lag accounted for by \phi. Notice that for infinite reduced frequency, the circulatory part of the lift amplitude is half that at k=0, and there is no phase lag angle.

These concepts are readily extended to more complicated boundary conditions, including the effects of pitch rate and combinations of frequencies. The results are obtained by linear superposition.

Apparent Mass Terms

It should also be appreciated that Theodorsen’s theory applies only to the “circulatory” component of the lift. Additional non-circulatory forces (and moments) are produced by inertial effects or “apparent mass” in an unsteady flow. This must be included to round out the unsteady aerodynamic theory if it is to be used in any predictive sense.

The apparent mass terms account for the pressure forces required to accelerate the fluid near the airfoil. For example, for a thin airfoil of chord c=2b moving normal to its surface at velocity w(t), the noncirculatory fluid force, F^{\rm nc}, acting on the surface is

(61)   \begin{equation*}F^{\rm nc} = -M_a \frac{dw}{dt}. \end{equation*}

The term M_a is known as the apparent mass and, in this case, is given by

(62)   \begin{equation*} M_a = \pi \varrho b^2 = \pi \varrho \frac{c^2}{4} \end{equation*}

Therefore, the non-circulatory lift for a planning motion where the airfoil moves normal to its surface with velocity w = \overbigdot{h} is

(63)   \begin{equation*} L_{\rm nc} = \pi \varrho \dfrac{c^2}{4} \ddot{h} \end{equation*}

or in terms of lift coefficient, then

(64)   \begin{equation*} C_{L_{\rm nc}} = \dfrac{\pi \varrho \dfrac{c^2}{4}}{\dfrac{1}{2} \varrho V_{\infty}^2 c} \ddot{h} = \pi b \dfrac{\ddot{h} }{V_{\infty}^2} \end{equation*}

Unsteady (Dynamic) Stall

Under unsteady conditions, airfoils and wings behave differently than in steady flow. In particular, their stall characteristics are very different, and the behavior is called a dynamic stall. This phenomenon has received much research interest because of the unsteady aerodynamic effects produced, which include the shedding of a leading-edge vortex, as shown in the simulation below.

A CFD animation of a dynamic stall on an oscillating airfoil. (With permission of Dr. George Barakos.)

The vortex shedding increases lift and drag, giving powerful nose-down (negative) pitching moments because of the aft-moving center of pressure induced by the shedding vortex, as shown in the figure below, based on measurements and flow visualization. Dynamic stall is known for exhibiting higher levels of maximum lift, drag, and pitching moments and hysteresis effects that can lead to flutter phenomena. Dynamic stall is particularly relevant in designing and analyzing helicopter rotors, wind turbine blades, and specific aircraft configurations.

Dynamic stall characteristics are very different from what happens under steady conditions.

Engineers and researchers study the phenomenon of dynamic stall to understand the effects better and develop strategies to mitigate any adverse impact on the performance of aerospace systems. In particular, helicopter rotor blades often experience dynamic stall during higher-speed forward flight or during maneuvers. Indeed, the onset of dynamic stall limits the helicopter’s operational flight envelope. Wind turbine blades may encounter dynamic stall during sudden changes in wind conditions, leading to variations in power output and high blade loads that may produce concerning vibration levels. Aircraft wings can also undergo dynamic stall during aggressive maneuvers, which is a behavior that must be considered when establishing the maneuvering flight envelope for military combat aircraft.

Summary & Closure

Unsteady flows occur whenever fluid properties fluctuate with respect to time. Understanding unsteady flows is crucial in various engineering and scientific disciplines, as they tend to be more prevalent in real-world problems. The practical implications of unsteady flows extend to advanced aerodynamic design and propulsion systems, necessitating advanced mathematical models and computational methods. Being able to categorize an unsteady flow is the first step in its analysis, which can often be quantified using a reduced frequency parameter. This chapter has shown that learning how to solve unsteady flow problems by following some classic exemplars helps expose the solution principles when considering the dimension of time.

5-Question Self-Assessment Quickquiz

For Further Thought or Discussion

  • Can you provide real-world situations where unsteady flow is crucial in fluid systems or engineering applications?
  • Why is reduced frequency an essential parameter for analyzing fluid flow oscillatory or vibratory motion?
  • In what scenarios might the Strouhal number be a relevant parameter to consider?
  • How does turbulence contribute to the complexity of fluid dynamics, and why is it often associated with unsteadiness?
  • Discuss the practical implications of unsteady flows in engineering design. How might engineers account for unsteady conditions in the design of systems like pipelines, aircraft, or water turbines?
  • How might unsteady flows impact the efficiency and performance of propulsion systems, such as those in aircraft or marine vehicles?


Other Useful Online Resources

To learn more about unsteady flows, take a look at some of these online resources:

  • A good video explaining the differences between steady and unsteady flows.
  • An explanation of unsteady flows, including the effects of thermodynamics.
  • An explanation of the conservation of mass in unsteady flow.
  • A video presentation on the application of the continuity equation to unsteady flows.

  1. Joukowsky, N.,  “Über den hydraulischen Stoss in Wasserleitungsröhren.” (“On the hydraulic hammer in water supply pipes.”) Mémoires de l'Académie Impériale des Sciences de St.-Pétersbourg, Series 8, 9(5), 1-71, 1900 (in German).
  2. Theodorsen, T., "General Theory for Aerodynamic Instability and the Mechanism of Flutter," NACA TR 496, 1935.
  3. Notice that this case does not represent an oscillating airfoil because the pitch-rate terms are not included.


Icon for the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License

Introduction to Aerospace Flight Vehicles Copyright © 2022 – 2024 by J. Gordon Leishman is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, except where otherwise noted.