# Viscous friction between layers

Boundary conditions on the top and bottom layers need to be added to close the system for the viscous stresses. We chose to impose a Neumann condition on the free-surface i.e. \displaystyle \partial_z u |_t = \dot{u}_t and a Navier slip condition on the bottom i.e. \displaystyle u|_b = u_b + \lambda_b \partial_z u|_b By default the viscosity is zero and we impose free-slip on the free-surface and no-slip on the bottom boundary i.e. \dot{u}_t = 0, \lambda_b = 0, u_b = 0.

double nu = 0.;
(const) scalar lambda_b = zeroc, dut = zeroc, u_b = zeroc;

For stability, we discretise the viscous friction term implicitly as \displaystyle \frac{(hu_l)^{n + 1} - (hu_l)^{\star}}{\Delta t} = \nu \left( \frac{u_{l + 1} - u_l}{h_{l + 1 / 2}} - \frac{u_l - u_{l - 1}}{h_{l - 1 / 2}} \right)^{n + 1} which can be expressed as the linear system \displaystyle \mathbf{Mu}^{n + 1} = \mathrm{rhs} where \mathbf{M} is a tridiagonal matrix. The lower, principal and upper diagonals are a, b and c respectively.

void vertical_viscosity (Point point, scalar h, scalar s, double dt)
{
double a[nl], b[nl], c[nl], rhs[nl];

The rhs of the tridiagonal system is h_lu_l.

  foreach_layer()
rhs[_layer] = s[]*h[];

The lower, principal and upper diagonals a, b and c are given by \displaystyle a_{l > 0} = - \left( \frac{\nu \Delta t}{h_{l - 1 / 2}} \right)^{n + 1} \displaystyle c_{l < \mathrm{nl} - 1} = - \left( \frac{\nu \Delta t}{h_{l + 1 / 2}} \right)^{n + 1} \displaystyle b_{0 < l < \mathrm{nl} - 1} = h_l^{n + 1} - a_l - c_l

  for (int l = 1; l < nl - 1; l++) {
a[l] = - 2.*nu*dt/(h[0,0,l-1] + h[0,0,l]);
c[l] = - 2.*nu*dt/(h[0,0,l] + h[0,0,l+1]);
b[l] = h[0,0,l] - a[l] - c[l];
}

For the top layer the boundary conditions give the (ghost) boundary value \displaystyle u_{\mathrm{nl}} = u_{\mathrm{nl} - 1} + \dot{u}_t h_{\mathrm{nl} - 1}, which gives the diagonal coefficient and right-hand-side \displaystyle b_{\mathrm{nl} - 1} = h_{\mathrm{nl} - 1}^{n + 1} - a_{\mathrm{nl} - 1} \displaystyle \mathrm{rhs}_{\mathrm{nl} - 1} = (hu)_{\mathrm{nl} - 1}^{\star} + \nu \Delta t \dot{u}_t

  a[nl-1] = - 2.*nu*dt/(h[0,0,nl-2] + h[0,0,nl-1]);
b[nl-1] = h[0,0,nl-1] - a[nl-1];
rhs[nl-1] += nu*dt*dut[];

For the bottom layer a third-order discretisation of the Navier slip condition gives \displaystyle \begin{aligned} b_0 & = h_0 + 2 \Delta t \nu \left( \frac{1}{h_0 + h_1} + \frac{h^2_1 + 3 h_0 h_1 + 3 h^2_0}{\det} \right),\\ c_0 & = - 2 \Delta t \nu \left( \frac{1}{h_0 + h_1} + \frac{h^2_0}{\det} \right),\\ \text{rhs}_0 & = (hu_0)^{\star} + 2 \Delta t \nu u_b \frac{h^2_1 + 3 h_0 h_1 + 2 h^2_0}{\det},\\ \det & = h_0 (h_0 + h_1)^2 + 2\lambda (3\,h_0 h_1 + 2\,h_0^2 + h_1^2), \end{aligned}

  double den = h[]*sq(h[] + h[0,0,1])
+ 2.*lambda_b[]*(3.*h[]*h[0,0,1] + 2.*sq(h[]) + sq(h[0,0,1]));
b[0] = h[] + 2.*dt*nu*(1./(h[] + h[0,0,1]) +
(sq(h[0,0,1]) + 3.*h[]*h[0,0,1] + 3.*sq(h[]))/den);
c[0] = - 2.*dt*nu*(1./(h[] + h[0,0,1]) + sq(h[])/den);
rhs[0] += 2.*dt*nu*u_b[]*(sq(h[0,0,1]) + 3.*h[]*h[0,0,1] + 2.*sq(h[0]))/den;

if (nl == 1) {
b[0] += c[0];
rhs[0] += (- c[0]*h[] - nu*dt) * dut[];
}

We can now solve the tridiagonal system using the Thomas algorithm.

  for (int l = 1; l < nl; l++) {
b[l] -= a[l]*c[l-1]/b[l-1];
rhs[l] -= a[l]*rhs[l-1]/b[l-1];
}
a[nl-1] = rhs[nl-1]/b[nl-1];
s[0,0,nl-1] = a[nl-1];
for (int l = nl - 2; l >= 0; l--)
s[0,0,l] = a[l] = (rhs[l] - c[l]*a[l+1])/b[l];
}

In the layered solver, vertical viscosity is applied to the velocity field just after advection, but before the pressure gradient/acceleration term is applied. To take the pressure gradient into account, we first apply the acceleration of the previous timestep, apply vertical viscosity and then substract the previous acceleration.

event viscous_term (i++,last)
{
if (nu > 0.) {
foreach() {
foreach_layer()
foreach_dimension()
u.x[] += dt*(ha.x[] + ha.x[1])/(hf.x[] + hf.x[1] + dry);
foreach_dimension()
vertical_viscosity (point, h, u.x, dt);
foreach_layer()
foreach_dimension()
u.x[] -= dt*(ha.x[] + ha.x[1])/(hf.x[] + hf.x[1] + dry);
}
boundary ((scalar *) {u});
}
}

## Horizontal diffusion

This approximates \displaystyle h \partial_t s = D \nabla \cdot (h \nabla s) with D the diffusion coefficient.

void horizontal_diffusion (scalar * list, double D, double dt)
{
if (D > 0.) {
scalar * d2sl = list_clone (list);
foreach_layer() {
foreach() {
scalar s, d2s;
for (s,d2s in list,d2sl) {
double a = 0.;
foreach_dimension()
a += (hf.x[]*fm.x[]/(cm[-1] + cm[])*(s[-1] - s[]) +
hf.x[1]*fm.x[1]/(cm[1] + cm[])*(s[1] - s[]));
d2s[] = 2.*a/(cm[]*sq(Delta));
}
}
foreach()
if (h[] > dry) {
scalar s, d2s;
for (s,d2s in list,d2sl)
s[] += dt*D*d2s[]/h[];
}
}
boundary (list);
delete (d2sl);
free (d2sl);
}
}

## References

 [popinet2020] Stéphane Popinet. A vertically-Lagrangian, non-hydrostatic, multilayer model for multiscale free-surface flows. Journal of Computational Physics, 418:109609, May 2020. [ DOI | http | .pdf ] [devita2019] Francesco De Vita, Pierre-Yves Lagrée, Sergio Chibbaro, and Stéphane Popinet. Beyond Shallow Water: appraisal of a numerical approach to hydraulic jumps based upon the Boundary Layer Theory. European Journal of Mechanics - B/Fluids, 79:233–246, 2019. [ DOI | http | .pdf ]