Stability of Poiseuille flow¶
Orr-Sommerfeld equation¶
The Orr-Sommerfeld (OS) equation is a fourth order differential equation for the modal small perturbation $v'$ over the incompressible "parallel" flow $U(y)$. The perturbations including $v'$ is supposed to be decomposed in several modes $v_i'(x,y,t)=\hat{v_i}(y) \exp(j(\omega_i t - \alpha_i x))$ where $\alpha_i$ are the streamwise wave number (real), and $\omega_i$ are the (complex) associated pulsations. This leads to the following eigenvalue problem.
\begin{equation*} \Bigg[ \frac{1}{i Re} (D^2 - \alpha^2)^2 - (\alpha U - \omega)(D^2 - \alpha^2) + \alpha U'' \Bigg] \hat{v} = 0 \end{equation*} where $D = \frac{d}{dy}$
The OS problem is discretized along y (x in the python class) and solved.
Channel incompressible flow¶
The OS problem is derived for the specific case of a laminar incompressible flow in a channel, named Poiseuille flow. The base state is $U(y)=U_{\max{}}(1-y^2)$. The problem is solved for a given pair of parameters $(\alpha, R_e)$ close to the critical point. Eigenvalues and eigenmodes are then selected and sorted by imaginary value (from most unstable to most stable). Darker blue circles show pair of eigenvalues.
import aerokit.stability.OrrSommerfeld as OS
import numpy as np
model = OS.Poiseuille(n=101, alpha=1.0, Reynolds=6000.0)
model.solve_eig()
vals, _, _ = model.select_and_sort(realmin=0., imagmin=-1.) # default is real order
import matplotlib.pyplot as plt
plt.plot(vals.real, vals.imag, "ob", markersize=5, alpha=0.3)
plt.xlabel("$\omega_r$")
plt.ylabel("$\omega_i$")
plt.grid()
48/101 remaining eigenvalues
Grid convergence and eigenmodes¶
Two models are run with respective grids 81 and 161. Spectra are superimposed. The coarsest grid provides a numerically sensitive branch of modes for most stables eigenvalues (as blue and red modes are not superimposed).
model = OS.Poiseuille(n=81, alpha=1.0, Reynolds=6000.0)
model.solve_eig()
vals0, _, _ = model.select_and_sort(realmin=0., imagmin=-1., sort="imag") # default is real order
model = OS.Poiseuille(n=161, alpha=1.0, Reynolds=6000.0)
model.solve_eig()
vals, vects, order = model.select_and_sort(realmin=0., imagmin=-1., sort="imag") # default is real order
51/81 remaining eigenvalues 48/161 remaining eigenvalues
plt.title("Spectrum")
plt.plot(vals0.real, vals0.imag, "ob", markersize=5, alpha=0.3)
plt.plot(vals.real, vals.imag, "or", markersize=3, alpha=0.9)
for i, v in enumerate(vals[order[:20]]):
plt.text(v.real, v.imag, ' '+str(i), fontsize='small')
plt.xlabel("$\omega_r$")
plt.ylabel("$\omega_i$")
plt.grid()
nv = 4
fig, ax = plt.subplots(nv, 2)
ax[0,0].set_title("u")
ax[0,1].set_title("v")
for i in range(nv):
ax[i,0].plot(model.x, model._diffop.matder(1) @ vects[:, order[i]].real)
ax[i,1].plot(model.x, vects[:, order[i]].real)