Hydrodynamic Stability: 2D Stability Analysis in Primitive Variables

Welcome back. I hope everyone is having a wonderful first week after spring break and is doing well going into the final stretch of the semester. In today’s blog post, I want to discuss how we could perform linear stability analysis using primitive velocity and pressure variables instead of intricately-derived streamfunction or vorticity quantities. One downside of primitive variables is that they can lead to larger matrices, which take longer to decompose into eigenvalues and eigenvectors. However, the primitive variable formulation can be beneficial because it includes only first and second-order derivatives, while the streamfunction and more complicated formulations can have higher-order derivatives that can introduce noise to the results. Therefore, we will provide an overview of how to formulate a 2D spatial stability problem in primitive variables. In this case, we will fix a real temporal frequency 𝜔 and compute a complex streamwise wavenumber 𝛼 = 𝛼_r + i𝛼_i, where -𝛼_i is the spatial growth rate of a temporally harmonic perturbation. 

If our base flow is a 2D parallel shear flow U(y), then the nondimensionalized, linearized Navier-Stokes equations for a 2D perturbation [u, v, p](x,y,t) are:

∂_tu + U ∂_xu + v dU/dy = -∂_xp + 1/Re (∂_x² + ∂_y²)u

∂_tv + U ∂_xv = -∂_yp + 1/Re (∂_x² + ∂_y²)v

0 = ∂_xu + ∂_yv. 

Now suppose that [u, v, p](x,y,t) = [uˆ, vˆ, pˆ](y) exp[i(𝛼x – 𝜔)] for real 𝜔 and complex 𝛼. Then, the linearized equations become:

-i𝜔uˆ = -i𝛼Uuˆ- i𝛼pˆ – dU/dy vˆ+ 1/Re (∂_y² – 𝛼²)uˆ

-i𝜔vˆ = -i𝛼Uvˆ- ∂_ypˆ+ 1/Re (∂_y² – 𝛼²)vˆ

0 = i𝛼uˆ + ∂_yvˆ. 

The idea with primitive variable analysis is that we compute our eigenvalues and eigenmodes directly from these equations. If we choose our state variable to be X = [uˆ, vˆ, pˆ]^T, then the three equations above can be cast into block matrix form, yielding the polynomial eigenvalue problem:

0 = (A₀ + 𝛼A₁ + 𝛼²A₂) X, 

where 

A₀ = [i𝜔I + ∂_y²/Re,   -dU/dy,  0 ;

          0,   i𝜔I + ∂_y²/Re,   -∂_y;

          0,   ∂_y,   0]

A₁ = [-iU,   0,   iI;

          0,   -iU,   0;

          iI,   0,   0]

A₂ = [-I/Re,   0,   0;

          0,   -I/Re,   0;

          0,   0,   0].

If U(y) is a boundary layer profile, then the perturbation velocity [u, v] is expected to vanish at the wall and in the far field. Therefore, the boundary conditions are 

uˆ = vˆ  = 0 at  y = 0 and as y → ∞. 

(Interestingly, we don’t need to impose boundary conditions for pˆ.) For boundary layers and channel flow, accurate results can be obtained if the derivative operator ∂_y is discretized using Chebyshev collocation methods. 

Next week, we will explore a new topic in hydrodynamic stability. Until then, please take care.

---