# Pipe flow (disk EVP)

## Script

"""
Dedalus script solving the linear stability eigenvalue problem for pipe flow.
This script demonstrates solving an eigenvalue problem in the periodic cylinder
using the disk basis and a parametrized axial wavenumber. It should take just
a few seconds to run (serial only).

The radius of the pipe is R = 1, and the problem is non-dimensionalized using
the radius and laminar velocity, such that the background flow is w0 = 1 - r**2.

No-slip boundary conditions are implemented on the velocity perturbations.
For incompressible hydro with one boundary, we need one tau term each for the
scalar axial velocity and vector horizontal (in-disk) velocity. Here we choose
to left the tau terms to the original (k=0) basis.

The eigenvalues are compared to the results of Vasil et al. (2016) [1] in Table 3.

To run, print, and plot the slowest decaying mode:
$python3 pipe_flow.py References: [1]: G. M. Vasil, K. J. Burns, D. Lecoanet, S. Olver, B. P. Brown, J. S. Oishi, "Tensor calculus in polar coordinates using Jacobi polynomials," Journal of Computational Physics (2016). """ import numpy as np import dedalus.public as d3 import matplotlib.pyplot as plt import logging logger = logging.getLogger(__name__) # Parameters Re = 1e4 kz = 1 m = 5 Nphi = 2 * m + 2 Nr = 64 dtype = np.complex128 # Bases coords = d3.PolarCoordinates('phi', 'r') dist = d3.Distributor(coords, dtype=dtype) disk = d3.DiskBasis(coords, shape=(Nphi, Nr), radius=1, dtype=dtype) phi, r = dist.local_grids(disk) # Fields s = dist.Field(name='s') u = dist.VectorField(coords, name='u', bases=disk) w = dist.Field(name='w', bases=disk) p = dist.Field(name='p', bases=disk) tau_u = dist.VectorField(coords, name='tau_u', bases=disk.edge) tau_w = dist.Field(name='tau_w', bases=disk.edge) tau_p = dist.Field(name='tau_p') # Substitutions dt = lambda A: s*A dz = lambda A: 1j*kz*A lift_basis = disk.derivative_basis(2) lift = lambda A: d3.Lift(A, lift_basis, -1) # Background w0 = dist.Field(name='w0', bases=disk.radial_basis) w0['g'] = 1 - r**2 # Problem problem = d3.EVP([u, w, p, tau_u, tau_w, tau_p], eigenvalue=s, namespace=locals()) problem.add_equation("div(u) + dz(w) = 0") problem.add_equation("dt(u) + w0*dz(u) + grad(p) - (1/Re)*(lap(u)+dz(dz(u))) + lift(tau_u) = 0") problem.add_equation("dt(w) + w0*dz(w) + u@grad(w0) + dz(p) - (1/Re)*(lap(w)+dz(dz(w))) + lift(tau_w) = 0") problem.add_equation("u(r=1) = 0") problem.add_equation("w(r=1) = 0") # Solver solver = problem.build_solver() sp = solver.subproblems_by_group[(m, None)] solver.solve_dense(sp) evals = solver.eigenvalues[np.isfinite(solver.eigenvalues)] evals = evals[np.argsort(-evals.real)] print(f"Slowest decaying mode: λ = {evals[0]}") solver.set_state(np.argmin(np.abs(solver.eigenvalues - evals[0])), sp.subsystems[0]) # Plot eigenfunction scales = (32, 4) ω = d3.div(d3.skew(u)).evaluate() ω.change_scales(scales) u.change_scales(scales) w.change_scales(scales) p.change_scales(scales) phi, r = dist.local_grids(disk, scales=scales) x, y = coords.cartesian(phi, r) cmap = 'RdBu_r' fig, ax = plt.subplots(2, 2, figsize=(6, 6)) ax[0,0].pcolormesh(x, y, u['g'][0].real, cmap=cmap) ax[0,0].set_title(r"$u_\phi$") ax[0,1].pcolormesh(x, y, u['g'][1].real, cmap=cmap) ax[0,1].set_title(r"$u_r$") ax[1,0].pcolormesh(x, y, w['g'].real, cmap=cmap) ax[1,0].set_title(r"$w$") ax[1,1].pcolormesh(x, y, p['g'].real, cmap=cmap) ax[1,1].set_title(r"$p\$")
for axi in ax.flatten():
axi.set_aspect('equal')
axi.set_axis_off()
fig.tight_layout()
fig.savefig("pipe_eigenfunctions.png", dpi=200)