Viscosity Notes 3

Viscosity Notes 3#

These are some supporting notes for the document Effective Viscosity.

Saint-Venant Equations#

Here we consider the following 1D version of the Navier-Stokes equations called the Saint-Venant equations:

\[\begin{gather*} \dot{n} + (nu)' = 0, \qquad D_{t}u = \dot{u} + uu_{,x} = -\frac{\Bigl(V(x) + \mathcal{E}'(n)\Bigr)_{,x}}{m} + \frac{\Bigl(\nu(n)u_{,x}\Bigr)_{,x}}{n}. \end{gather*}\]

A Physical Analogy: Shallow Water

These equations follow from the flow of an incompressible fluid of mass density \(\rho\) like water flowing along a channel of width \(w\) whose bottom height \(b(x)\) varies, but not enough to excite appreciable vertical dynamics. If the fluid height is \(\eta(x)\), then the internal energy density is given by the gravitational potential energy (here \(a_g\) is the acceleration due to gravity)

\[\begin{gather*} \int_{b(x)}^{b(x) + \eta(x)} a_g \rho w z\;\d{z} = a_g w\rho \frac{(b+\eta)^2 - b^2}{2} = a_g w\rho \frac{\eta^2}{2} + a_g w b \eta \end{gather*}\]

This suggests identifying the following 1D number density and coupling constant to match the GPE.

\[\begin{gather*} n_{1}(x) \equiv \frac{\rho w \eta}{m}, \qquad g \equiv \frac{m^2 a_g}{w \rho^2}, \qquad V(x) \equiv \frac{m a_g b(x)}{\rho}. \end{gather*}\]

To help with intuition, we will sometimes plot our flow \(n(x, t)\) and \(u(x, t)\) in terms of this channel flow:

\[\begin{gather*} \eta(x) \equiv \frac{m}{\rho w}n(x), \qquad b(x) \equiv \frac{\rho}{ma_g} V(x) = \frac{m}{\rho w} \frac{V(x)}{g} \end{gather*}\]

For the purposes of plotting, we will choose the width of the channel \(w\) and fluid density \(\rho\) such that \(m/\rho w = 1\).

Comparison with Wikipedia

To compare with the conservative form of the shallow water equations, note that, in 2D, the following are equivalent, after application of the continuity equation:

\[\begin{gather*} \rho \eta D_t u = \pdiff{(\rho \eta u)}{t} + \pdiff{\rho\eta u^2}{x} + \pdiff{\rho\eta u v}{y}. \end{gather*}\]

Their form does not include the potential \(V(x)\) (they consider a flat bed). Explicitly

\[\begin{gather*} \dot{n} = - \sum_{i}\nabla_{i}(nv_{i}), \qquad D_t v_i = \dot{v}_{i} + \sum_{j} v_j \nabla_{i}v_i,\\ \pdiff{n v_i}{t} + \sum_{j}\nabla_j(nv_iv_j) = \underbrace{n \dot{v}_i + n\sum_{j}v_j\nabla_j(v_i)}_{nD_tv_i} + \underbrace{\dot{n}v_i + v_i\sum_{j}\nabla_j(nv_j)}_{0}. \end{gather*}\]

Numerical Implementation#

We start with a simple numerical implementation of these equations on a periodic lattice. This should work well when everything is smooth, but will likely fail when the density drops to zero which will produce a kink. Another potential failure mode appears if the local flow velocity exceeds the speed of sound. Finite viscosity allows for mild supersonic flow

Show code cell content

Hide code cell content

from scipy.integrate import solve_ivp
from gpe.contexts import FPS

%pylab inline
from scipy.integrate import solve_ivp
from gpe.contexts import FPS

Lx = 10.0
Nx = 128
dx = Lx/Nx
x = np.arange(Nx) * dx - Lx/2
k = 2*np.pi * np.fft.fftfreq(Nx, dx)

# Fourier-space representation of derivatives
D = {1: None, 2: None}
D[1] = 1j*k
assert Nx % 2 == 0
D[1][Nx//2] = 0  # Kill highest unmatched frequency
D[2] = -k**2

def diff(f, d=1):
    """Compute the derivative of f."""
    global D
    return np.fft.ifft(D[d] * np.fft.fft(f)).real

def pack(n, u):
    return np.concatenate([n, u])
    
def unpack(y):
    return np.reshape(y, (2, Nx))

def get_V_m(t, x, d=0):
    """Return the dth derivative of V/m."""
    global V0_m, sigma, dV_m_dx
    V_m = V0_m * np.exp(-(x/sigma)**2/2)
    if d == 0:
        dV = dV_m_dx*x
        return V_m + dV - np.min(dV)  # Constant to make V_min = 0
    elif d == 1:
        return - x * V_m / sigma**2 + dV_m_dx

def get_E_m(n, d=0):
    """Return the dth derivative of E/m where E is the energy density."""
    global g_m
    if d == 0:
        return g_m * n**2 / 2 
    elif d == 1:
        return g_m * n
    elif d == 2:
        return g_m

def get_nu_n(n, d=0):
    """Return the dth derivative of the vicosity diff(nu(n), n, d)/n, divided by n."""
    if d == 0:
        return nu # * n / n
    elif d == 1:
        return nu / n
    else:
        return 0
        
def compute_dy_dt(t, y):
    global get_V_m, get_E_m, get_nu_n
    n, u = unpack(y)
    j = u*n    
    dn_dt = -diff(j)
    n_x = diff(n)
    u_x = diff(u)
    u_xx = diff(u, d=2)
    du_dt = -u*u_x
    du_dt -= get_V_m(t, x, d=1) + get_E_m(n, d=2)*n_x
    du_dt += get_nu_n(n, d=1)*n_x*u_x + get_nu_n(n) * u_xx / n
    
    # Correction needed to keep finite j = j_avg
    du = j_avg/n - u
    du_dt += j_corr_factor * du
    return pack(dn_dt, du_dt)
    
# Quick test of derivatives
k = 2*np.pi / Lx
s, c = np.sin(k*x), np.cos(k*x)
f = np.exp(s)
df = k*c*f
ddf = k**2*f*(c**2 - s)
assert np.allclose(diff(f), df)
assert np.allclose(diff(f, 2), ddf)

%pylab is deprecated, use %matplotlib inline and import the required libraries.
Populating the interactive namespace from numpy and matplotlib

In this example, we add a gradient to our potential, noting that only \(V'(x)\) appears in our equations of motion. Thus, although the potential is not periodic, the force is. This constant is needed to ensure that a finite flow remains.

Viscosity Notes 3

Contents

Viscosity Notes 3#

Saint-Venant Equations#

Numerical Implementation#