Welcome back. In this post, we will discuss the celebrated Orr-Sommerfeld equation, which is an eigenvalue problem governing the 2D stability of viscous parallel shear flows. Although this equation might initially seem narrow in scope since it only accounts for two spatial dimensions, it can actually be applied to a substantial range of 3D flows. This works because of Squire’s theorem, which states that 3D parallel shear flows first become unstable from 2D disturbances, which start to grow at lower Reynolds numbers than 3D disturbances. This powerful result has allowed the Orr-Sommerfeld equation to be applied to boundary layers on 3D airplane wings to reveal how they might become turbulent. The equation has also been applied to channel and pipe flows, and has typically found them to become unstable at lower Reynolds numbers than what simpler, inviscid equations have predicted. The pipe flow investigations have been particularly insightful because inviscid equations suggest that they are linearly stable at all Reynolds numbers, which cannot be true in light of Osborne Reynolds’ experiments of 1883. Since the Orr-Sommerfeld equation accounts for viscosity, it has shown pipe flows to become unstable at a finite Reynolds number. This is significant because it suggests that viscosity can *destabilize* a fluid, an unintuitive result that goes against the traditional wisdom that viscosity dampens disturbances. For the remainder of this post, we will use the general form of the linearized Navier-Stokes equations to derive the Orr-Sommerfeld equation for 2D parallel flows. We will then cast the Orr-Sommerfeld equation in terms of a non-normal differential operator acting on a perturbation to a stream function. The consequences of non-normality will be explored in more detail next week.

In a previous blog post, we linearized the Navier-Stokes equations around a base flow with velocity **U** and pressure P to determine how infinitesimal velocity perturbations, **u**’, and pressure perturbations, p’, might grow or decay. The resulting disturbance equations were

∂**u**/∂t** **+ **U** • ∇**u** + **u** • ∇**U** – 1/Re ∇^{2}**u = **– ∇p,

∇ • **u** = 0.

In 2D parallel shear flows, the base flow velocity can be expressed as the vector **U** = U(y) **e**_{x}, where **e**_{x} is the unit vector in the x-direction. The velocity perturbation can also be expressed as the 2D vector **u**’ = u(x,y,t) **e**_{x} + v(x,y,t) **e**_{y}. With these simplifications, the disturbance equations can be expressed as follows:

*x-momentum:* ∂u/∂t + U ∂u/∂x + v dU/dy – 1/Re ∇^{2}u = -∂p/∂x,

*y-momentum:* ∂v/∂t + U ∂v/∂x – 1/Re ∇^{2}v = -∂p/∂y,

*continuity:* ∂u/∂x + ∂v/∂y = 0.

Since 2D parallel shear flows are homogeneous in the streamwise x-direction and in time, we can express u, v, and p in terms of wave functions of x and t. These wave functions take the form of the complex exponential exp[ik(x-ct)], where k∊ℝis a streamwise *wavenumber* and c = c_{r} + ic_{i} ∊ℂis a complex wave speed:

[u, v, p](x,y,t) = [û, v̂, p̂](y) exp[ik(x-ct)].

The amplitudes û, v̂, and p̂ are unknown functions of y, and must be determined using the linearized disturbance equations. If the assumed form of u, v, and p are substituted into these equations, the resulting expressions are:

*x-momentum:* ik (U-c) û + v̂ U’ – 1/Re (D^{2}-k^{2}) û = -ik p̂,

*y-momentum:* ik (U-c) v̂ – 1/Re (D^{2}-k^{2}) v̂ = -dp̂/dy,

*continuity:* ik û + dv̂/dy = 0.

In these equations dU/dy has been substituted with U’ and the operator D = d/dy represents differentiation with respect to y. These relations above are three coupled equations in three unknowns, [û, v̂, p̂], which are still not easy to deal with. To create a simpler expression with just one known, we can introduce a stream function, 𝜓. This stream function acts as a Hamiltonian for the flow, and the velocity components are expressed in terms of partial derivatives of 𝜓:

u = ∂𝜓/∂y,

v = -∂𝜓/∂x.

If u and v are expressed using a stream function, the continuity constraint is satisfied automatically and the remaining analysis can focus on the x and y-momentum equations. These momentum equations can be further simplified by first assuming that

𝜓(x,y,t) = 𝜑(y) exp[ik(x-ct)].

This decomposition of the stream function implies that û = d𝜑/dy = D𝜑 and v̂ = -ik𝜑. Therefore, the following relations hold:

*x-momentum:* ik (U-c) D𝜑 – ik U’𝜑 – D/Re (D^{2}-k^{2}) 𝜑 = -ik p̂,

*y-momentum:* k^{2} (U-c) 𝜑 + ik/Re (D^{2}-k^{2}) = -Dp̂.

If the x-momentum equation is multiplied on the left by D and the y-momentum equation is multiplied on the right by ik, then the right-hand-sides of both equations become equal to -ik Dp̂. This ultimately means that

ik (U-c) D^{2}𝜑 + ik U’ D𝜑 – ik U’ D𝜑 – ik U’’𝜑 – D^{2}/Re (D^{2}-k^{2}) 𝜑 =

ik k^{2} (U-c) 𝜑 – k^{2}/Re (D^{2}-k^{2}).

With a little more rearrangement, we obtain the Orr-Sommerfeld equation:

[(U-c)(D^{2}-k^{2}) – U’’] 𝜑 = -i/(k Re) (D^{2}-k^{2})^{2} 𝜑.

This is a 4th order, variable coefficient, linear differential equation that determines a stream function amplitude 𝜑(y) and complex wave speed c. If c_{i} > 0, then the stream function blows up. This means that u, v, and p also blow up, and the perturbation determined by 𝜑(y) is asymptotically unstable. If c_{i} = 0, then the flow is neutrally stable and perturbations neither grow nor decay asymptotically. In the case that c_{i} < 0, the flow is asymptotically stable, and the perturbation determined by 𝜑(y) decays as t → ∞.

It turns out that there are infinitely many amplitudes 𝜑(y) and complex wave speeds c that satisfy the Orr-Sommerfeld equation. Each of these pairs (𝜑, c) is an eigenfunction-eigenvalue pair of the Orr-Sommerfeld equation. The eigenfunctions 𝜑 determine characteristic spatial structures in the flow, the real part of c determines the temporal frequencies at which these spatial structures oscillate, and the imaginary part of c determines the rate at which these spatial structures grow or decay in time. These eigenfunctions and eigenvalues can be computed numerically by turning the linear differential equation into a matrix equation. This can be done by approximating the infinite-dimensional operator D = d/dy with a finite-dimensional differentiation matrix, and then using that finite-dimensional matrix to approximate the operators

A = i/(k Re) (D^{2}-k^{2})^{2} + U (D^{2}-k^{2}) – U’’,

B = D^{2}-k^{2}.

If these operators are discretized, then the Orr-Sommerfeld equation can be solved as the generalized eigenvalue problem,

A𝜑 = cB𝜑.

Since B is invertible, this generalized eigenvalue problem can be turned into the traditional eigenvalue problem,

L𝜑 = c𝜑,

where L = B^{-1}A is the Orr-Sommerfeld operator. The eigenfunctions of L belong to an inner product space equipped with an energy norm. Within this inner product space, it can be shown that L does not commute with its Hermitian adjoint, L^{ad}. This means that L is non-normal, and hence has non-orthogonal eigenfunctions. These non-orthogonal eigenfunctions can combine together to produce surprising results, such as significant transient energy growth, even if all of the eigenvalues are stable. In the next post, we will explore this transient energy growth and its implications on the transition to turbulence. Until then, please take care.

## Leave a Reply