kh.c

Lock exchange (Kelvin–Helmoltz shear instability)
- References
- Code

Lock exchange (Kelvin–Helmoltz shear instability)

A rectangular horizontal box is filled with a fluid with a non-uniform density. The left half of the domain is “warmer” and thus lighter than the right half. The fluid is viscous and the density field is diffusive with a Péclet number of unity (i.e. the kinematic viscosity of the fluid is identical to the diffusivity of the density).

This is a similar test case to that proposed by Berntsen et al. 2006 to illustrate the importance of non-hydrostatic effects. The main difference is the symmetric slip boundary conditions on the top and bottom boundaries.

Animation of the density field

The results can be directly compared with the solution obtained using the (completely different) centered Navier–Stokes solver. The agreement is excellent as illustrated below.

set term svg enhanced size 1000,170 font ",10"
set size ratio -1
unset key
unset ytics
set xrange [-4:4]
plot 'log' u 1:2 w l lc rgbcolor "black"

Ten equally-spaced contours of density at t= 21. (script)

plot '../kh-ns/log' u 1:2 w l lc rgbcolor "black"

Results with the Navier–Stokes solver. (script)

The hydrostatic solution is different.

plot '../kh-hydro/log' u 1:2 w l lc rgbcolor "black"

Hydrostatic solution. (script)

References

[berntsen2006]

Jarle Berntsen, Jiuxing Xing, and Guttorm Alendal. Assessment of non-hydrostatic ocean models using laboratory scale problems. Continental Shelf Research, 26(12-13):1433–1447, 2006. [ DOI ]

Code

We use an x-z grid and the non-hydrostatic multilayer solver.

#include "grid/multigrid1D.h"
#include "layered/hydro.h"
#if HYDRO
# include "layered/implicit.h"
#else
# include "layered/nh.h"
#endif

The relative density difference is 10^{-2}.

#define drho(T) (- 0.01*(T))
#include "layered/dr.h"

#include "layered/remap.h"
#include "layered/perfs.h"

The domain is 8 \times 1.

int main()
{
  L0 = 8.;
  X0 = -L0/2.;
  G = 9.81;
  nl = 64;
  N = nl*L0;
  nu = 2e-4;

  theta_H = 0.51;
  cell_lim = mono_limit;
  DT = 0.03;

We use a slip boundary condition at the bottom, rather than the default no-slip (the top boundary is slip by default).

  const vector slip[] = {HUGE};
  lambda_b = slip;
  
  run();

We generate movie and graphics at the end of the run.

  system ("for f in plot-*.png; do convert $f ppm:- && rm -f $f; done | "
	  "ppm2mp4 T.mp4");
#if HYDRO
  system ("gnuplot -e 'set table' kh-hydro.plot | sed '/^#.*/d' > log");
#else
  system ("gnuplot -e 'set table' kh.plot | sed '/^#.*/d' > log");
#endif
}

The initial density profile is a smooth hyperbolic tangent function.

event init (i = 0)
{
  foreach() {
    double z = zb[];
    foreach_layer() {
      h[] = 1./nl;
      z += h[]/2.;
      T[] = - 0.5*tanh((x + 0.1*cos(pi*z/2.))/0.04);
      z += h[]/2.;
    }
  }
}

Horizontal viscosity and diffusion of density (with the same viscosity/diffusion coefficient).

event viscous_term (i++)
{
#if NH
  horizontal_diffusion ({u,w,T}, nu, dt);
#else
  horizontal_diffusion ({u,T}, nu, dt);
#endif // NH
  
  foreach()
    vertical_diffusion (point, h, T, dt, nu, 0, 0, HUGE);
}

Density field at t= 21.

event density (t = 21)
{
  foreach() {
    double z = zb[];
    foreach_layer() {
      z += h[]/2.;
      printf ("%g %g %g\n", x, z, T[]);
      z += h[]/2.;
    }
    printf ("\n");
  }
  fflush (stdout);
}

We generate an animation of the density field.

event gnuplot (t += 1; t <= 30)
{
  static FILE * fp = popen ("gnuplot 2> /dev/null", "w");
  if (i == 0) {
    system ("rm -f plot*.png");
    fprintf (fp,
	     "set term png font ',9' size 1000,250\n"
	     "set pm3d map\n" // interpolate 10,4\n"
	     "set size ratio -1\n"
	     "unset key\n"
	     "unset colorbox\n"
	     "unset ytics\n"
	     "# jet colormap\n"
	     "set palette defined ( 0 0 0 0.5647, 0.125 0 0.05882 1, 0.25"
	     " 0 0.5647 1, 0.375 0.05882 1 0.9333, 0.5 0.5647 1 0.4392, "
	     "0.625 1 0.9333 0, 0.75 1 0.4392 0, 0.875 0.9333 0 0,"
	     " 1 0.498 0 0 )\n"
	     "set cbrange [-0.5:0.5]\n");
  }
  fprintf (fp,
	   "set output 'plot-%02g.png'\n"
	   "sp [%g:%g][0:1.001]'-' u 1:2:3\n",
	   t, X0, X0 + L0);
  foreach (serial) {
    double z = zb[];
    fprintf (fp, "%g %g %g\n", x, z, T[]);
    foreach_layer() {
      z += h[];
      fprintf (fp, "%g %g %g\n", x, z, T[]);
    }
    fprintf (fp, "\n");
  }
  fprintf (fp, "e\n\n");
  //  fprintf (fp, "pause 0.05\n");
  fflush (fp);
}