Classical Hydrodynamics

Hide code cell source

 
import mmf_setup; mmf_setup.nbinit(quiet=True)
import os
from pathlib import Path
FIG_DIR = Path(mmf_setup.ROOT) / '../Docs/_build/figures/'
os.makedirs(FIG_DIR, exist_ok=True)
import logging; logging.getLogger("matplotlib").setLevel(logging.CRITICAL)
import numpy as np, matplotlib.pyplot as plt
try: from myst_nb import glue
except: glue = None

Classical Hydrodynamics#

This notebook describes the feature of the code to solve simple classical hydrodynamics problems. We use these to match to coarse-grained superfluid hydrodynamics simulations, for example, to extract the effective viscosity.

Example#

As a sample problem, we consider here a harmonically trapped gas oscillating back and forth in a harmonic trap with a small barrier potential inserted. The idea is that the barrier will produce some drag, slowing the oscillations as energy is dissipated. For comparison, a superfluid will not lose energy, but will generate phonons (wave turbulence) resulting in a slow decay of the center of mass oscillations.

from gpe.Examples.hydro_1d import Experiment, u

expt = Experiment()
state = expt.get_state()
state.plot()
state.evolve(steps=200, t__final=10*u.ms);

[I 03:37:43 numexpr.utils] NumExpr defaulting to 2 threads.
[I 03:37:44 root] Patching zope.interface.document.asReStructuredText to format code

/home/docs/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/hydro_1d.py:122: UserWarning: Coupling constant g should have the correct dimensions.
  warnings.warn("Coupling constant g should have the correct dimensions.")

Linearization of Equation#

Linearization of the above equation is little difficult in the current form. So let’s do some redefinitions: $\( h=-\frac{V}{g} \qquad \eta = n-h \qquad g' = \frac{g}{m} \)\( Then the shallow water equation in 1D takes the form \)\( \dot{\eta} = -\frac{\partial (\eta+h)u}{\partial x} \\ \dot{u} = -\frac{\partial (\frac{u^2}{2}+g'\eta)}{\partial x} + \nu \frac{\partial^2 u}{\partial x^2} \)$

The jacobian matrix for transforming fluxes for \(X = (n, nu)\)is given by:

(16)#\[\begin{align} \mathbf{J} = \begin{bmatrix} u & h\\ g' & u\\ \end{bmatrix} \end{align}\]

The eigenvalues are \(u \pm \sqrt{g'h}\)

Now, we can linearize this equation to get $\( \dot{\eta} = -\frac{\partial hu}{\partial x} \\ \dot{u} = -\frac{\partial (g\eta)}{\partial x} + \nu \frac{\partial^2 u}{\partial x^2} \)\( I am little unclear on whether to include the double derivative in linear version. For simpler version of stability analysis, we drop the double derivative \)\( \dot{\eta} = -\frac{\partial hu}{\partial x} \\ \dot{u} = -\frac{\partial (g\eta)}{\partial x} \)\( i.e. \)\( \dot{Y} = -\mathbf{C} \frac{\partial Y}{\partial x} \)$ where

(17)#\[\begin{align} \mathbf{C} = \begin{bmatrix} 0 & h\\ g & 0\\ \end{bmatrix} \end{align}\]
(18)#\[\begin{align} \mathbf{Y} = \begin{bmatrix} \eta\\ u\\ \end{bmatrix} \end{align}\]

Von-Neumann Stability Analysis: Convection Equation#

Consider the Finite Difference discretization $\( \frac{Y_{i}^{r+1} - Y_{i}^{r}}{\delta t} = - C \frac{Y_{i+1}^{r} - Y_{i}^{r}}{\delta x} \)\( Putting, \)Y_{i}^{r} = A^r (1, 1) e^{ikx}\( \)\( \frac{A-1}{\delta t} = -C \frac{e^{ik\delta x} - 1}{\delta x} \\ A = 1 - \frac{\delta t}{\delta x}C\left(e^{ik\delta x} - 1\right) \)\( \)\( |A|^2 = (1-\frac{\delta t}{\delta x}C (\cos(k \delta x)-1))^2 + (\frac{\delta t}{\delta x}C \sin(k \delta x))^2 \)$

\[ |A|^2 = (1 + \frac{|C| \delta t}{\delta x})^2 - 2 \frac{\delta t}{\delta x}|C| \cos(k \delta x) (1 +\frac{\delta t}{\delta x}|C|) + (\frac{\delta t}{\delta x}C)^2 \]

Now, \(|A|^2 \leq 1\), for stability $\( (1+\frac{\delta t}{\delta x}|C|)(1-\cos(k \delta x)) \leq 0 \)$ Hence it is unconditionally unstable

The Formulation (Finite Volume Method):#

(19)#\[\begin{gather} \dot{\rho} = -\frac{\partial F}{\partial x} + \mathcal{S} \\ \frac{d \rho}{dt} = - \frac{1}{\bigtriangleup x} (F_{xR} - F_{xL}) + \mathcal{S}_{x} \\ F_{xR} = \frac{1}{2} \left((F_{x+dx} + F_{x}) - |\alpha_{R}|(\rho_{x+dx} - \rho_{x}) \right) \\ F_{xL} = \frac{1}{2} ((F_{x} + F_{x-dx}) - |\alpha_{L}|(\rho_{x} - \rho_{x-dx})) \\ \alpha_I = \frac{\bigtriangleup F}{\bigtriangleup \rho} \quad \bigtriangleup \rho > 1e-6 \\ = \frac{\partial F}{\partial \rho} \quad otherwise \end{gather}\]
(20)#\[\begin{gather} \frac{\partial j}{\partial t} = -\frac{1}{2*\bigtriangleup x} \left(F_{x+dx} - F_{x-dx} - (\alpha_R j_{x+dx} - (\alpha_R + \alpha_L)j_{x} + \alpha_L j_{x-dx}) \right) + \mathcal{S_x} \end{gather}\]