Hydrodynamic Stability: Three-Dimensional Linear Stability and Input-Output Analysis of Parallel Shear Flows

Introduction:

Welcome back. Today, I am excited to return to our series on hydrodynamic stability. In previous blogs, we considered the stability of two-dimensional parallel shear flows. Now, we will derive the governing equations for stability analysis of three-dimensional parallel shear flows. Our mathematical work here will give us the tools to consider oblique modes and three-dimensional instabilities that are characteristic of channel flows. We will also have the opportunity to derive the Orr-Sommerfeld equation in physical space, rather than Fourier space, to get more practice manipulating the Navier-Stokes equations using vector calculus operations. This Orr-Sommerfeld equation will be coupled with the Squire equation for the wall-normal vorticity, which will allow us to study three-dimensional effects. The resulting system of the Orr-Sommerfeld and Squire equations will form a dynamical system for the wall-normal velocity and vorticity, which will allow us to compute instability modes of the unforced flow, and also map input disturbances to output amplifications in velocity [1]. The resulting flow structures identified from this analysis can be beneficial in designing control strategies that encourage or suppress turbulence in channel flows. 

Orr-Sommerfeld Equation in Physical Coordinates:

Consider a base flow of the form U(y) = U(y) e_x, where y is the wall-normal coordinate and e_x is the unit vector in the streamwise direction. Any base flow of this form will be called a parallel shear flow, because the only velocity component is in the streamwise direction, and varies in the cross-stream direction. Now, consider infinitesimal disturbances of the form u(x, y, z, t) = [u e_x + v e_y + w e_z](x, y, z, t), where x is the streamwise coordinate and z is the spanwise coordinate. Assume that all quantities are non-dimensionalized. If we assume very small disturbances, then the non-dimensionalized Navier-Stokes equations can be linearized around the base flow to yield the following: 

∂_tu + U ∂_xu + U’ v = -∂_x p + 1/Re ∆u + F_x       (x-momentum)

∂_tv + U ∂_xv = -∂_y p + 1/Re ∆v + F_y       (y-momentum)

∂_tw + U ∂_xw = -∂_z p + 1/Re ∆w + F_z       (z-momentum)

∂_xu + ∂_yv + ∂_zw = 0,        (continuity)

where U’ = dU/dy, Re is the Reynolds number, ∆ = ∂_xx² + ∂_yy² + ∂_zz² is the Laplacian, and F_x, F_y, and F_z are input forcing terms in the x, y, and z directions, respectively. To reduce our problem from four variables (u, v, w, p) to three variables (u, v, w), we can eliminate pressure from the above equations. We do this by taking the divergence of the (vector) momentum equation, which results in the following equation for the Laplacian of pressure:

-∆p + ∂_xF_x + ∂_yF_y + ∂_zF_z = ∂_x [U ∂_xu + U’ v] + ∂_y [U ∂_xv] + ∂_z [U ∂_xw] =

= 2U’ ∂_xv + U ∂_x [∂_xu + ∂_yv + ∂_zw].

Recognizing that ∂_xu + ∂_yv + ∂_zw = 0 yields the pressure Poisson equation: 

∆p = -2U’ ∂_xv + ∂_xF_x + ∂_yF_y + ∂_zF_z.

For computational efficiency, it is useful to cast the three equations for u, v, and w into just two equations for wall-normal velocity v and wall-normal vorticity 𝜔_y. To do this, let’s first derive the Orr-Sommerfeld equation for v by taking the Laplacian of the y-momentum equation:

∂_t ∆v + ∆ [U ∂_xv] = -∂_y ∆p + 1/Re ∆²v + ∆F_y.

The term ∆ [U ∂_xv] can be simplified as follows: 

∆ (U ∂_xv) = U ∂_x( ∂_xx² + ∂_zz²)v + ∂_yy² [U ∂_xv] = 

 = U ∂_x( ∂_xx² + ∂_zz²)v + ∂_y [U’ ∂_xv + U ∂_xy² v] = 

 = U’’ ∂_xv + 2 U’ ∂_xy² v + U ∂_x ∆v.

Therefore, the Laplacian of the y-momentum equation becomes:

∂_t ∆v + U’’ ∂_xv + 2 U’ ∂_xy² v + U ∂_x ∆v = 

-∂_y [-2U’ ∂_xv + ∂_xF_x + ∂_yF_y + ∂_zF_z] + 1/Re ∆²v + ∆F_y = 

2U’’ ∂_xv + 2U’ ∂_xy² v – ∂_xy² F_x – ∂_yy² F_y – ∂_yz² F_z + 1/Re ∆²v + ∆F_y.

With a little more simplification, we arrive at the Orr-Sommerfeld equation in physical space: 

∂_t ∆v = -U ∂_x ∆v + U’’ ∂_xv + 1/Re ∆²v – ∂_xy² F_x + (∂_xx² + ∂_zz²) F_y – ∂_yz² F_z.        (Orr-Sommerfeld)

Squire Equation for Wall-Normal Vorticity: 

The next step in formulating our stability problem is to derive an equation for the wall-normal vorticity, defined as: 

𝜔_y = ∂_zu – ∂_xw.

To form our equation for 𝜔_y, we first take the z-derivative of the x-momentum equation:

∂_t ∂_z u + U ∂_x ∂_z u + U’ ∂_zv = -∂_x ∂_z p + 1/Re ∆ ∂_zu + ∂_zF_x.

Now, take the x-derivative of the z-momentum equation:

∂_t ∂_x w + U ∂_x ∂_x w = -∂_z ∂_x p + 1/Re ∆ ∂_x w + ∂_x F_z.

Subtracting the x-derivative of the z-momentum equation from the z-derivative of the x-momentum equation yields:

∂_t (∂_z u – ∂_x w) + U ∂_x (∂_z u – ∂_x w) + U’ ∂_z v = 1/Re ∆ (∂_z u – ∂_x w) + ∂_z F_x – ∂_x F_z ,

which simplifies to the Squire equation: 

∂_t 𝜔_y = – U’ ∂_z v – U ∂_x 𝜔_y + 1/Re ∆𝜔_y + ∂_z F_x – ∂_x F_z ,        (Squire)

Output Equations:

Using the Orr-Sommerfeld and Squire equations, we can compute instability modes and input-output relationships in terms of v and 𝜔_y. However, it is also useful to relate these quantities back to the streamwise and spanwise velocity components u and w. To form expressions for u and w in terms of v and 𝜔_y, we can compute the x and z-derivatives of 𝜔_y and apply the continuity equation as follows: 

∂_x𝜔_y =  ∂_x (∂_zu – ∂_xw) = ∂_z (- ∂_yv – ∂_zw) – ∂_xx²w = -∂_zy²v – (∂_xx² + ∂_zz²)w,

∂_z𝜔_y = ∂_z (∂_zu – ∂_xw) = ∂_zz²u – ∂_x (-∂_xu – ∂_yv) = ∂_xy²v + (∂_xx² + ∂_zz²)u.

Therefore, output equations for u and w are as follows: 

u = – (∂_xx² + ∂_zz²)^-1 (∂_xy²v – ∂_z𝜔_y),        (u-output)

w = – (∂_xx² + ∂_zz²)^-1 (∂_yz²v + ∂_x𝜔_y),        (w-output)

Linear Control System Formulation:

The Orr-Sommerfeld, Squire, and output equations can then be formally cast into the following state-space form, where F = [F_x   F_y   F_z]^T is the input vector, 𝜓 = [v   𝜔_y]^T is the state vector, and 𝜉 = [u   v   w]^T is the output vector: 

∂_t 𝜓(x, y, z, t) = [А 𝜓(t)](x, y, z) + [B F(t)](x, y, z)

𝜉(x, y, z, t) = [C 𝜓(t)](x, y, z).

The operator A is defined as: 

А = [L_OS,   0;

         L_C,   L_SQ],

where the Orr-Sommerfeld, coupling, and Squire operators are respectively:

L_OS = ∆^-1 (-U ∂_x ∆ + U’’ ∂_x + 1/Re ∆²)

L_C = -U’ ∂_z

L_SQ = -U ∂_x + 1/Re ∆.

The input matrix is:

B = [- ∆^-1 ∂_xy²,  ∆^-1 (∂_xx² + ∂_zz²),    – ∆^-1 ∂_yz²;

        ∂_z ,   0,   – ∂_x].

Lastly, the output matrix is: 

C = [- (∂_xx² + ∂_zz²)^-1 ∂_xy²,   (∂_xx² + ∂_zz²)^-1 ∂_z;

        1,   0;

        – (∂_xx² + ∂_zz²)^-1 ∂_zy², – (∂_xx² + ∂_zz²)^-1 ∂_x].

A common set of boundary conditions for these equations, arising from the no-slip and continuity equations at the walls at y = ± 1 are: 

v = ∂_yv = 𝜔_y = 0 at y = ± 1 for all x and z.

If the flow is homogeneous in the x and z directions, then periodic boundary conditions can be used in x and z. This state-space model can serve several purposes. Since A is in block-diagonal form, the eigenvalues of A are determined by the eigenvalues of L_OS and L_SQ. Therefore, the three-dimensional instability modes of channel flows can be separated between modes arising from the Orr-Sommerfeld equation and modes arising from the Squire equation. Furthermore, the impact of disturbances on the velocity field can be assessed using the transfer function, 

H(𝜔) = C (i𝜔 – A)^-1 B.

A singular value decomposition of the transfer function reveals that the most amplified disturbances come in the form of streamwise vortices and streaks, rather than the unstable eigenmodes. This is a reflection of the non-normality of the operator A, which arises both from the non-normality of L_OS and L_SQ, and from the off-diagonal coupling term L_C.

Following [2], the transfer function H(𝜔) can be split into six different components, which map each input (F_x, F_y, F_z) to each velocity output (u,v,w):

H(𝜔) = [H_ux(𝜔),   H_uy(𝜔),   H_uz(𝜔);

              H_vx(𝜔),   H_vy(𝜔),   H_vz(𝜔);

              H_wx(𝜔),   H_wy(𝜔),   H_wz(𝜔)].

Analysis of each of these components can provide significant insight into the precise directions of disturbances that amplify each of the velocity components, and therefore provide a more nuanced understanding than just analyzing H(𝜔) alone. 

Linear Control System in Fourier Space:

Finally, we mention that the Orr-Sommerfeld and Squire equations derived above can be Fourier transformed with respect to the homogeneous coordinates x and z. This yields a dynamical system in terms of the streamwise wavenumber k_x and spanwise wavenumber k_z. In this case, the input is the Fourier-transformed forcing vector F(k_x, k_z) = [F_x   F_y   F_z]^T(k_x, k_z), the state is the Fourier-transformed wall-normal velocity and vorticity 𝜓(k_x, k_z) = [v   𝜔_y]^T(k_x, k_z), and the output is the Fourier-transformed velocity field 𝜉(k_x, k_z) = [u   v   w]^T(k_x, y_z) is the output vector. The matrix equation is again: 

∂_t 𝜓(k_x, y, k_z, t) = [А(k_x, k_z) 𝜓(k_x, k_z, t)](y) + [B(k_x, k_z) F(k_x, k_z, t)](y)

𝜉(k_x, y, k_z, t) = [C(k_x, k_z) 𝜓(k_x, k_z, t)](y).

However, the operator A is defined as: 

А(k_x, k_z) = [L_OS(k_x, k_z),   0;

         L_C(k_x, k_z),   L_SQ(k_x, k_z)],

where the Orr-Sommerfeld, coupling, and Squire operators are respectively:

L_OS(k_x, k_z) = ∆^-1 (-U ik_x ∆ + U’’ ik_x + 1/Re ∆²)

L_C(k_x, k_z) = -U’ ik_z

L_SQ(k_x, k_z) = -U ik_x + 1/Re ∆,

where the Laplacian is now ∆ = ∂_yy² – k_x² – k_z².

The input matrix is:

B(k_x, k_z) = [- ∆^-1 ik_x∂_y,  – ∆^-1 (k_x² + k_z²),    – ∆^-1 ik_z∂_y;

        ik_z ,   0,   – ik_x].

Lastly, the output matrix is: 

C(k_x, k_z) = [(k_x² + k_z²)^-1 ik_x∂_y,   – (k_x² + k_z²)^-1 ik_z;

        1,   0;

        – (k_x² + k_z²)^-1 ik_z∂_y, (k_x² + k_z²)^-1 ik_x].

In future blog posts, we will explore the discoveries made from these equations in more detail. Until then, please take care. 

References:

[1] M. R. Jovanovic and B. Bamieh. The spatio-temporal impulse response of the linearized Navier-Stokes equations. In Proceedings of the 2001 American Control Conference, Arlington, VA, pages 1948-1953, 2001. 
[2] M. R. Jovanovic and B. Bamieh. Componentwise energy amplification in channel flows. J. Fluid Mech., 534:145-183, July 2005.

---