CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutSign UpSign In

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place. Commercial Alternative to JupyterHub.

| Download

Bálint Kaszás and George Haller, 'Capturing the Edge of Chaos as a Spectral Submanifold in Pipe Flows ', Journal of Fluid Mechanics, submitted, Figure 5 key words: spectral submanifold, reduced-order model, edge of chaos, pipe flow

Views: 4
Visibility: Unlisted (only visible to those who know the link)
Image: ubuntu2204
Tweet
Share
Share
Kernel: Python 3 (system-wide)

logo

Spectral analysis of the laminar flow and the lower branch traveling wave

In this notebook we present calculations related to the spectrum of the two fixed point solutions, the lower branch traveling wave and the laminar flow. The eigenvalues are obtained using the Newton-Krylov algorithm built into Openpipeflow (developed by Ashley Willis, https://openpipeflow.org/).

  • We show the leading eigenvalues of the lower branch traveling wave and check that the nonresonance conditions required by results of Haller et al. 2023 hold.

  • We show the calculations necessary to explicitly compute the analytic (primary) slow SSM of the laminar flow

Check for resonances in the spectrum of the lower branch traveling wave

import numpy as np

eig2LB = np.loadtxt('data/arnoldi.dat', skiprows = 3)

The existence results for the mixed-mode SSM apply if Sternberg's linearization theorem applies. Let us denote the eigenvalues as λi\lambda_i and ordering them according to descending real part. We then require that there is no resonance among eigenvalues, i.e

λj≠∑imiλi∑imi≥2\lambda_j \neq \sum_{i} m_i \lambda_i \quad \sum_i m_i \geq 2

We select the leading 1515 eigenvalues, but this corresponds to 2424 eigenvalues in total, due to the complex-conjugate pairs. 

ms = np.arange(0, 5, 1)

eigenvalues = eig2LB[np.abs(eig2LB[:,1])>0.001, 1] + 1j*eig2LB[np.abs(eig2LB[:,1])>0.001,2]

eigenvalues = np.unique(np.concatenate((eigenvalues, np.conjugate(eigenvalues))))

sort_ind = np.argsort(-np.real(eigenvalues), axis=0)

eigenvalues = np.take_along_axis(eigenvalues, sort_ind, axis = 0)

from sklearn.preprocessing import PolynomialFeatures

polyfeat = PolynomialFeatures(degree = 6)


polyfeat.fit(np.ones(eigenvalues.shape).reshape(1,-1)) # we use PolynomialFeatures to get all linear combinations

We create an array that enumerates all possible integer linear combinations of the eigenvalues up to order 5. 

poww = polyfeat.powers_

poww.shape
(593775, 24)
lincombs = []

for i in range(poww.shape[0]):
    if np.sum(poww[i,:]>1):
        summ = 0.
        for j in range(poww.shape[1]):
            summ = summ + poww[i,j]*eigenvalues[j]
        lincombs.append(summ)
    else:
        lincombs.append(1e10)

From the linear combinations ∑imiλi\sum_{i} m_i \lambda_i we select the one that is closest to every eigenvalue, that is we take

D(λj)=minm∣λj−∑imiλi∣D(\lambda_j) = \text{min}_m |\lambda_j - \sum_{i} m_i \lambda_i |
mindist = []
minlincomb = []
for i in range(len(eigenvalues)):
    distances = []
    for j in range(len(lincombs)):
        distances.append(np.abs(eigenvalues[i] - lincombs[j]))
    amax = np.argmin(distances)
    mindist.append(distances[amax])
    minlincomb.append(poww[amax,:])

We then compute a relative metric of resonance-nearness, D(λj)∣λj∣ \frac{D(\lambda_j)}{|\lambda_j|} 

relative_dist = mindist/np.abs(eigenvalues)

The eigenvalue closest to resonance is

closest_index = np.argmin(relative_dist)

closest_index
20
relative_dist[closest_index]
0.0024162315256366787

The closest near-resonance observed is at 

from sympy import latex
from IPython.display import display_latex



def disp(idx, symObj):
    eqn = '$' + idx + ' ' + latex(symObj) + '$'
    display_latex(eqn,raw=True)
    return

disp('\lambda_{20}=', eigenvalues[closest_index])

[\lambda_{20}= \mathtt{\text{(-0.065289236+0.0082860247j)}}]

and the closest linear combination is

disp('\lambda_{1} + \lambda_{4} + 3 \lambda_{8}=', eigenvalues[0] + eigenvalues[3] + 3*eigenvalues[7])

[\lambda_{1} + \lambda_{4} + 3 \lambda_{8}= \mathtt{\text{(-0.065141413+0.0082274116j)}}]

eigenvalues[(minlincomb[closest_index])>0]
array([ 0.01976398-0.j , -0.01511716+0.00822741j, -0.02326274-0.j ])
f = plt.figure()
import matplotlib
ratio = 8/5
ax = f.add_subplot(122)
ax2 = f.add_subplot(121)

f.set_size_inches(3.38*2*1.7, (3.38*2)/ratio, forward=True)
ax.loglog(np.abs(eigenvalues), mindist/np.abs(eigenvalues), '.', c='black')
ax.text(.96*np.abs(eigenvalues)[20], 1.05*(mindist/np.abs(eigenvalues))[20], '$\\lambda_{20}$', fontsize = 17)
ax.loglog(np.abs(eigenvalues)[20], (mindist/np.abs(eigenvalues))[20], 'o', c='crimson')

ax.set_xlabel(r'$|\lambda|$', fontsize = 25)
ax.set_ylabel(r'$\frac{D(\lambda)}{|\lambda|}$', rotation=0, labelpad = 20, fontsize = 20)
ax2.plot([-0.075, 0.025], [0,0], '--', c='black', alpha = 0.5, linewidth =0.8)
ax2.plot([0, 0.0], [-0.035,0.035], '--', c='black', alpha = 0.5, linewidth =0.8)

ax2.plot(np.real(eigenvalues), np.imag(eigenvalues), 'x', c='black')
ax2.plot(np.real(eigenvalues)[20], np.imag(eigenvalues)[20], 'o', c='crimson')
ax2.text(1.05*np.real(eigenvalues)[20], 1.18*np.imag(eigenvalues)[20], '$\\lambda_{20}$', fontsize = 17)

ax2.plot(np.real(eigenvalues)[0], np.imag(eigenvalues)[0], 'o', c='orange')
ax2.text(.95*np.real(eigenvalues)[0], 0.0015 + np.imag(eigenvalues)[0], '$\\lambda_{1}$', fontsize = 17)

ax2.plot(np.real(eigenvalues)[3], np.imag(eigenvalues)[3], 'o', c='orange')
ax2.text(1.05*np.real(eigenvalues)[3], 1.18*np.imag(eigenvalues)[3], '$\\lambda_{4}$', fontsize = 17)


ax2.plot(np.real(eigenvalues)[7], np.imag(eigenvalues)[7], 'o', c='orange')
ax2.text(1.15*np.real(eigenvalues)[7], 0.0019+np.imag(eigenvalues)[7], '$3\\lambda_{8}$', fontsize = 17)


ax2.set_xlim(-0.075, 0.025)
ax2.set_ylim(-0.035, 0.035)

ax2.set_xlabel(r'real$(\lambda)$', fontsize = 20)
ax2.set_ylabel(r'imag$(\lambda)$', fontsize = 20)

from matplotlib import ticker

ax.set_xlim(0.006, 0.1)
#ax.set_ylim(0.005, 0.06)
f.tight_layout()
Image in a Jupyter notebook

Calculating the spectrum of the laminar flow

import sympy as sy

from sympy.vector import CoordSys3D, Del
from sympy.vector import divergence, gradient, curl, matrix_to_vector, dot, cross


r, phi, z, t = sy.symbols('r \\varphi z t')
ur = sy.Function('u_r')(r, phi, z, t)
uphi = sy.Function('u_\\varphi')(r, phi, z, t)
uz = sy.Function('u_z')(r, phi, z, t)
p = sy.Function('p')(r, phi, z, t)
nu = sy.symbols('\\nu')

R = CoordSys3D('R', transformation='cylindrical', variable_names=("r", "phi", "z"))
to_cyl_sub = {r:R.r, phi:R.phi, z:R.z}
from_cyl_sub = {R.r:r, R.phi:phi, R.z:z}

u = ur.subs(to_cyl_sub)*R.i + uphi.subs(to_cyl_sub)*R.j + uz.subs(to_cyl_sub)*R.k

umatrix = u.to_matrix(R).subs(from_cyl_sub)

laplacian = sy.simplify(curl(curl(u)) - gradient(divergence(u))).to_matrix(R).subs(from_cyl_sub)

ugradu = sy.simplify((gradient(dot(u,u)) - 2* cross(u, curl(u)))/2).to_matrix(R)

gradp = gradient(p.subs(to_cyl_sub)).to_matrix(R).subs(from_cyl_sub)

nonlinearity = ugradu.subs(from_cyl_sub)

navierstokes = umatrix.diff(t) + ugradu.subs(from_cyl_sub) + gradp - nu * laplacian

from sympy import latex
from IPython.display import display_latex



def disp(idx, symObj):
    eqn = '\\[' + idx + ' ' + latex(symObj) + '\\]'
    display_latex(eqn,raw=True)
    return

disp('\\mathcal{N}(\mathbf{u})=',navierstokes)

N(u)=[−ν(−r2∂2∂r2ur⁡(r,φ,z,t)−r2∂2∂z2ur⁡(r,φ,z,t)−r∂∂rur⁡(r,φ,z,t)+ur⁡(r,φ,z,t)+2∂∂φuφ⁡(r,φ,z,t)−∂2∂φ2ur⁡(r,φ,z,t))r2+ur⁡(r,φ,z,t)∂∂rur⁡(r,φ,z,t)+uz⁡(r,φ,z,t)∂∂zur⁡(r,φ,z,t)+∂∂rp(r,φ,z,t)+∂∂tur⁡(r,φ,z,t)−uφ⁡2(r,φ,z,t)r+uφ⁡(r,φ,z,t)∂∂φur⁡(r,φ,z,t)r−ν(−r2∂2∂r2uφ⁡(r,φ,z,t)−r2∂2∂z2uφ⁡(r,φ,z,t)−r∂∂ruφ⁡(r,φ,z,t)+uφ⁡(r,φ,z,t)−∂2∂φ2uφ⁡(r,φ,z,t)−2∂∂φur⁡(r,φ,z,t))r2+ur⁡(r,φ,z,t)∂∂ruφ⁡(r,φ,z,t)+uz⁡(r,φ,z,t)∂∂zuφ⁡(r,φ,z,t)+∂∂tuφ⁡(r,φ,z,t)+uφ⁡(r,φ,z,t)ur⁡(r,φ,z,t)r+uφ⁡(r,φ,z,t)∂∂φuφ⁡(r,φ,z,t)r+∂∂φp(r,φ,z,t)r−ν(−∂2∂r2uz⁡(r,φ,z,t)−∂2∂z2uz⁡(r,φ,z,t)−∂∂ruz⁡(r,φ,z,t)r−∂2∂φ2uz⁡(r,φ,z,t)r2)+ur⁡(r,φ,z,t)∂∂ruz⁡(r,φ,z,t)+uz⁡(r,φ,z,t)∂∂zuz⁡(r,φ,z,t)+∂∂zp(r,φ,z,t)+∂∂tuz⁡(r,φ,z,t)+uφ⁡(r,φ,z,t)∂∂φuz⁡(r,φ,z,t)r]\mathcal{N}(\mathbf{u})= \left[\begin{matrix}- \frac{\nu \left(- r^{2} \frac{\partial^{2}}{\partial r^{2}} \operatorname{u_{r}}{\left(r,\varphi,z,t \right)} - r^{2} \frac{\partial^{2}}{\partial z^{2}} \operatorname{u_{r}}{\left(r,\varphi,z,t \right)} - r \frac{\partial}{\partial r} \operatorname{u_{r}}{\left(r,\varphi,z,t \right)} + \operatorname{u_{r}}{\left(r,\varphi,z,t \right)} + 2 \frac{\partial}{\partial \varphi} \operatorname{u_{\varphi}}{\left(r,\varphi,z,t \right)} - \frac{\partial^{2}}{\partial \varphi^{2}} \operatorname{u_{r}}{\left(r,\varphi,z,t \right)}\right)}{r^{2}} + \operatorname{u_{r}}{\left(r,\varphi,z,t \right)} \frac{\partial}{\partial r} \operatorname{u_{r}}{\left(r,\varphi,z,t \right)} + \operatorname{u_{z}}{\left(r,\varphi,z,t \right)} \frac{\partial}{\partial z} \operatorname{u_{r}}{\left(r,\varphi,z,t \right)} + \frac{\partial}{\partial r} p{\left(r,\varphi,z,t \right)} + \frac{\partial}{\partial t} \operatorname{u_{r}}{\left(r,\varphi,z,t \right)} - \frac{\operatorname{u_{\varphi}}^{2}{\left(r,\varphi,z,t \right)}}{r} + \frac{\operatorname{u_{\varphi}}{\left(r,\varphi,z,t \right)} \frac{\partial}{\partial \varphi} \operatorname{u_{r}}{\left(r,\varphi,z,t \right)}}{r}\\- \frac{\nu \left(- r^{2} \frac{\partial^{2}}{\partial r^{2}} \operatorname{u_{\varphi}}{\left(r,\varphi,z,t \right)} - r^{2} \frac{\partial^{2}}{\partial z^{2}} \operatorname{u_{\varphi}}{\left(r,\varphi,z,t \right)} - r \frac{\partial}{\partial r} \operatorname{u_{\varphi}}{\left(r,\varphi,z,t \right)} + \operatorname{u_{\varphi}}{\left(r,\varphi,z,t \right)} - \frac{\partial^{2}}{\partial \varphi^{2}} \operatorname{u_{\varphi}}{\left(r,\varphi,z,t \right)} - 2 \frac{\partial}{\partial \varphi} \operatorname{u_{r}}{\left(r,\varphi,z,t \right)}\right)}{r^{2}} + \operatorname{u_{r}}{\left(r,\varphi,z,t \right)} \frac{\partial}{\partial r} \operatorname{u_{\varphi}}{\left(r,\varphi,z,t \right)} + \operatorname{u_{z}}{\left(r,\varphi,z,t \right)} \frac{\partial}{\partial z} \operatorname{u_{\varphi}}{\left(r,\varphi,z,t \right)} + \frac{\partial}{\partial t} \operatorname{u_{\varphi}}{\left(r,\varphi,z,t \right)} + \frac{\operatorname{u_{\varphi}}{\left(r,\varphi,z,t \right)} \operatorname{u_{r}}{\left(r,\varphi,z,t \right)}}{r} + \frac{\operatorname{u_{\varphi}}{\left(r,\varphi,z,t \right)} \frac{\partial}{\partial \varphi} \operatorname{u_{\varphi}}{\left(r,\varphi,z,t \right)}}{r} + \frac{\frac{\partial}{\partial \varphi} p{\left(r,\varphi,z,t \right)}}{r}\\- \nu \left(- \frac{\partial^{2}}{\partial r^{2}} \operatorname{u_{z}}{\left(r,\varphi,z,t \right)} - \frac{\partial^{2}}{\partial z^{2}} \operatorname{u_{z}}{\left(r,\varphi,z,t \right)} - \frac{\frac{\partial}{\partial r} \operatorname{u_{z}}{\left(r,\varphi,z,t \right)}}{r} - \frac{\frac{\partial^{2}}{\partial \varphi^{2}} \operatorname{u_{z}}{\left(r,\varphi,z,t \right)}}{r^{2}}\right) + \operatorname{u_{r}}{\left(r,\varphi,z,t \right)} \frac{\partial}{\partial r} \operatorname{u_{z}}{\left(r,\varphi,z,t \right)} + \operatorname{u_{z}}{\left(r,\varphi,z,t \right)} \frac{\partial}{\partial z} \operatorname{u_{z}}{\left(r,\varphi,z,t \right)} + \frac{\partial}{\partial z} p{\left(r,\varphi,z,t \right)} + \frac{\partial}{\partial t} \operatorname{u_{z}}{\left(r,\varphi,z,t \right)} + \frac{\operatorname{u_{\varphi}}{\left(r,\varphi,z,t \right)} \frac{\partial}{\partial \varphi} \operatorname{u_{z}}{\left(r,\varphi,z,t \right)}}{r}\end{matrix}\right]

uhp = sy.Matrix([0,0,1-r**2])

We now linearize around the parabolic profile

uHP=(001−r2)\mathbf{u}_{HP} = \begin{pmatrix}0 \\ 0 \\ 1-r^2 \end{pmatrix}
epsilon = sy.symbols('\\varepsilon')


Linearized_Ns = navierstokes.subs([(ur, epsilon*ur),
                                   (uphi, epsilon*uphi),
                                   (uz, 1-r**2 + epsilon * uz),
                                   (p, 4*nu*z + epsilon*p)])


Linearized_Ns = Linearized_Ns.diff(epsilon).subs(epsilon, 0).doit()

disp('D\mathcal{N}(\mathbf{u}_{HP})=',Linearized_Ns)

DN(uHP)=[−ν(−r2∂2∂r2ur⁡(r,φ,z,t)−r2∂2∂z2ur⁡(r,φ,z,t)−r∂∂rur⁡(r,φ,z,t)+ur⁡(r,φ,z,t)+2∂∂φuφ⁡(r,φ,z,t)−∂2∂φ2ur⁡(r,φ,z,t))r2+(1−r2)∂∂zur⁡(r,φ,z,t)+∂∂rp(r,φ,z,t)+∂∂tur⁡(r,φ,z,t)−ν(−r2∂2∂r2uφ⁡(r,φ,z,t)−r2∂2∂z2uφ⁡(r,φ,z,t)−r∂∂ruφ⁡(r,φ,z,t)+uφ⁡(r,φ,z,t)−∂2∂φ2uφ⁡(r,φ,z,t)−2∂∂φur⁡(r,φ,z,t))r2+(1−r2)∂∂zuφ⁡(r,φ,z,t)+∂∂tuφ⁡(r,φ,z,t)+∂∂φp(r,φ,z,t)r−ν(−∂2∂r2uz⁡(r,φ,z,t)−∂2∂z2uz⁡(r,φ,z,t)−∂∂ruz⁡(r,φ,z,t)r−∂2∂φ2uz⁡(r,φ,z,t)r2)−2rur⁡(r,φ,z,t)+(1−r2)∂∂zuz⁡(r,φ,z,t)+∂∂zp(r,φ,z,t)+∂∂tuz⁡(r,φ,z,t)]D\mathcal{N}(\mathbf{u}_{HP})= \left[\begin{matrix}- \frac{\nu \left(- r^{2} \frac{\partial^{2}}{\partial r^{2}} \operatorname{u_{r}}{\left(r,\varphi,z,t \right)} - r^{2} \frac{\partial^{2}}{\partial z^{2}} \operatorname{u_{r}}{\left(r,\varphi,z,t \right)} - r \frac{\partial}{\partial r} \operatorname{u_{r}}{\left(r,\varphi,z,t \right)} + \operatorname{u_{r}}{\left(r,\varphi,z,t \right)} + 2 \frac{\partial}{\partial \varphi} \operatorname{u_{\varphi}}{\left(r,\varphi,z,t \right)} - \frac{\partial^{2}}{\partial \varphi^{2}} \operatorname{u_{r}}{\left(r,\varphi,z,t \right)}\right)}{r^{2}} + \left(1 - r^{2}\right) \frac{\partial}{\partial z} \operatorname{u_{r}}{\left(r,\varphi,z,t \right)} + \frac{\partial}{\partial r} p{\left(r,\varphi,z,t \right)} + \frac{\partial}{\partial t} \operatorname{u_{r}}{\left(r,\varphi,z,t \right)}\\- \frac{\nu \left(- r^{2} \frac{\partial^{2}}{\partial r^{2}} \operatorname{u_{\varphi}}{\left(r,\varphi,z,t \right)} - r^{2} \frac{\partial^{2}}{\partial z^{2}} \operatorname{u_{\varphi}}{\left(r,\varphi,z,t \right)} - r \frac{\partial}{\partial r} \operatorname{u_{\varphi}}{\left(r,\varphi,z,t \right)} + \operatorname{u_{\varphi}}{\left(r,\varphi,z,t \right)} - \frac{\partial^{2}}{\partial \varphi^{2}} \operatorname{u_{\varphi}}{\left(r,\varphi,z,t \right)} - 2 \frac{\partial}{\partial \varphi} \operatorname{u_{r}}{\left(r,\varphi,z,t \right)}\right)}{r^{2}} + \left(1 - r^{2}\right) \frac{\partial}{\partial z} \operatorname{u_{\varphi}}{\left(r,\varphi,z,t \right)} + \frac{\partial}{\partial t} \operatorname{u_{\varphi}}{\left(r,\varphi,z,t \right)} + \frac{\frac{\partial}{\partial \varphi} p{\left(r,\varphi,z,t \right)}}{r}\\- \nu \left(- \frac{\partial^{2}}{\partial r^{2}} \operatorname{u_{z}}{\left(r,\varphi,z,t \right)} - \frac{\partial^{2}}{\partial z^{2}} \operatorname{u_{z}}{\left(r,\varphi,z,t \right)} - \frac{\frac{\partial}{\partial r} \operatorname{u_{z}}{\left(r,\varphi,z,t \right)}}{r} - \frac{\frac{\partial^{2}}{\partial \varphi^{2}} \operatorname{u_{z}}{\left(r,\varphi,z,t \right)}}{r^{2}}\right) - 2 r \operatorname{u_{r}}{\left(r,\varphi,z,t \right)} + \left(1 - r^{2}\right) \frac{\partial}{\partial z} \operatorname{u_{z}}{\left(r,\varphi,z,t \right)} + \frac{\partial}{\partial z} p{\left(r,\varphi,z,t \right)} + \frac{\partial}{\partial t} \operatorname{u_{z}}{\left(r,\varphi,z,t \right)}\end{matrix}\right]

k, m= sy.symbols('k m', integers=True)

alpha = sy.symbols('\\alpha')
lamb = sy.symbols('\\lambda')

ii = sy.sqrt(-1)

ar = sy.Function('a_r')(r)
aphi = sy.Function('a_\\varphi')(r)
az = sy.Function('a_z')(r)
ap = sy.Function('a_p')(r)

This operator can be decomposed to normal modes of the form

u=aeikαzeikmφeλt\mathbf{u} = ae^{ik\alpha z} e^{ikm\varphi} e^{\lambda t}

We expect that the least stable mode is obtained with m=0m=0, in which case the equation simplifies to the eigenvalue equation

Normalmode_ns = Linearized_Ns.subs([(ur, ar*sy.exp(ii*m*phi + ii*k*alpha*z)*sy.exp(lamb*t)),
                   (uphi,aphi*sy.exp(ii*m*phi + ii*k*alpha*z)*sy.exp(lamb*t)),
                   (uz,az*sy.exp(ii*m*phi + ii*k*alpha*z)*sy.exp(lamb*t)),
                   (p, ap*sy.exp(ii*m*phi + ii*k*alpha*z)*sy.exp(lamb*t))]).doit()

Normalmode_ns = sy.simplify(Normalmode_ns/(sy.exp(ii*m*phi + ii*k*alpha*z)*sy.exp(lamb*t)))

equation = sy.Matrix([eq for eq in Normalmode_ns])

simplified_eq = equation.subs(k, 0).subs(ar,0).subs(aphi, 0).subs(ap,0).doit()

disp('0=',simplified_eq[2])

0=λaz⁡(r)−νm2az⁡(r)r2+νd2dr2az⁡(r)+νddraz⁡(r)r0= \lambda \operatorname{a_{z}}{\left(r \right)} - \frac{\nu m^{2} \operatorname{a_{z}}{\left(r \right)}}{r^{2}} + \nu \frac{d^{2}}{d r^{2}} \operatorname{a_{z}}{\left(r \right)} + \frac{\nu \frac{d}{d r} \operatorname{a_{z}}{\left(r \right)}}{r}