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.

In [1]:

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()

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).

In [2]:

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

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

In [3]:

plt.title("Spectrum")
plt.plot(vals0.real, vals0.imag, "ob", markersize=5, alpha=0.3, label="$n=81$")
plt.plot(vals.real, vals.imag, "or", markersize=3, alpha=0.9, label="$n=161$")
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() ; plt.legend()

plt.title("Spectrum") plt.plot(vals0.real, vals0.imag, "ob", markersize=5, alpha=0.3, label="$n=81$") plt.plot(vals.real, vals.imag, "or", markersize=3, alpha=0.9, label="$n=161$") 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() ; plt.legend()

Out[3]:

<matplotlib.legend.Legend at 0x75105cf5dc60>

In [4]:

nv = 4
fig, ax = plt.subplots(nv, 2, figsize=(8, 2*nv), sharex=True)
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)

nv = 4 fig, ax = plt.subplots(nv, 2, figsize=(8, 2*nv), sharex=True) 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)

Maximum amplification¶

In [5]:

def max_mode(alpha, Rey, n=81):
    model = OS.Poiseuille(n=n, alpha=alpha, Reynolds=Rey)
    model.solve_eig()
    vals0, _, _ = model.select_and_sort(realmin=0., imagmin=-1., verbose=False)
    return np.max(np.imag(vals0))

Alpha = np.linspace(0.9, 1.2, 10)
Rey = np.linspace(5000, 7000, 3)
for i, Rey in enumerate(Rey):
    growthrate = np.zeros(Alpha.shape)
    for j, alpha in enumerate(Alpha):
        growthrate[j] = max_mode(alpha, Rey)
    plt.plot(Alpha, growthrate, label="Re = "+str(Rey))
plt.xlabel("$\\alpha$")
plt.ylabel("$\\omega_i$")
plt.grid(); plt.legend()

def max_mode(alpha, Rey, n=81): model = OS.Poiseuille(n=n, alpha=alpha, Reynolds=Rey) model.solve_eig() vals0, _, _ = model.select_and_sort(realmin=0., imagmin=-1., verbose=False) return np.max(np.imag(vals0)) Alpha = np.linspace(0.9, 1.2, 10) Rey = np.linspace(5000, 7000, 3) for i, Rey in enumerate(Rey): growthrate = np.zeros(Alpha.shape) for j, alpha in enumerate(Alpha): growthrate[j] = max_mode(alpha, Rey) plt.plot(Alpha, growthrate, label="Re = "+str(Rey)) plt.xlabel("$\\alpha$") plt.ylabel("$\\omega_i$") plt.grid(); plt.legend()

Out[5]:

<matplotlib.legend.Legend at 0x75105c86e950>

In [6]:

import scipy.optimize as opt

def max_growthrate(Rey, n=81):
    # Maximize f over alpha for given Re
    res = opt.minimize_scalar(lambda alpha: -max_mode(alpha, Rey, n), bounds=(.9, 1.2), method='bounded')
    return -res.fun, res.x  # return max value

Alpha = np.linspace(0.9, 1.2, 20)
Rey = np.linspace(5000, 7000, 10)
max_growth = np.zeros((len(Rey), 2))
grid_alpha, grid_Rey = np.meshgrid(Alpha, Rey)
alphamax = np.zeros_like(Rey)
for i, Re in enumerate(Rey):
    alphamax[i] = max_growthrate(Re, n=41)[1]
    for j, alpha in enumerate(Alpha):
        growthrate[i, j] = max_mode(alpha, Re, n=41)
plt.contourf(grid_alpha, grid_Rey, growthrate, levels=41, cmap='OrRd')
plt.contour(grid_alpha, grid_Rey, growthrate, levels=[0.], colors='k', linewidths=2)
plt.plot(alphamax, Rey, 'b--', label="max growth rate")
plt.xlabel("$\\alpha$")
plt.ylabel("Re")
plt.legend()

import scipy.optimize as opt def max_growthrate(Rey, n=81): # Maximize f over alpha for given Re res = opt.minimize_scalar(lambda alpha: -max_mode(alpha, Rey, n), bounds=(.9, 1.2), method='bounded') return -res.fun, res.x # return max value Alpha = np.linspace(0.9, 1.2, 20) Rey = np.linspace(5000, 7000, 10) max_growth = np.zeros((len(Rey), 2)) grid_alpha, grid_Rey = np.meshgrid(Alpha, Rey) alphamax = np.zeros_like(Rey) for i, Re in enumerate(Rey): alphamax[i] = max_growthrate(Re, n=41)[1] for j, alpha in enumerate(Alpha): growthrate[i, j] = max_mode(alpha, Re, n=41) plt.contourf(grid_alpha, grid_Rey, growthrate, levels=41, cmap='OrRd') plt.contour(grid_alpha, grid_Rey, growthrate, levels=[0.], colors='k', linewidths=2) plt.plot(alphamax, Rey, 'b--', label="max growth rate") plt.xlabel("$\\alpha$") plt.ylabel("Re") plt.legend()

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[6], line 16
     14     alphamax[i] = max_growthrate(Re, n=41)[1]
     15     for j, alpha in enumerate(Alpha):
---> 16         growthrate[i, j] = max_mode(alpha, Re, n=41)
     17 plt.contourf(grid_alpha, grid_Rey, growthrate, levels=41, cmap='OrRd')
     18 plt.contour(grid_alpha, grid_Rey, growthrate, levels=[0.], colors='k', linewidths=2)

IndexError: too many indices for array: array is 1-dimensional, but 2 were indexed

Critical Reynolds number¶

In [7]:

# Root-finding for g(Re) = 1
Re_crit, info = opt.newton(lambda Rey: max_growthrate(Rey)[0], 6000, rtol=1e-6, full_output=True)
alpha_crit, niter = max_growthrate(Re_crit, n=41)[1], info.iterations  # return max value and number of iterations
#
print(f"Critical Re = {Re_crit:.1f} with alpha = {alpha_crit:.3f} within {niter} iterations")

# Root-finding for g(Re) = 1 Re_crit, info = opt.newton(lambda Rey: max_growthrate(Rey)[0], 6000, rtol=1e-6, full_output=True) alpha_crit, niter = max_growthrate(Re_crit, n=41)[1], info.iterations # return max value and number of iterations # print(f"Critical Re = {Re_crit:.1f} with alpha = {alpha_crit:.3f} within {niter} iterations")

Critical Re = 6000.6 with alpha = 1.016 within 3 iterations