# Breaking Stokes wave

We solve the two-phase Navier–Stokes equations using a momentum-conserving transport of each phase. Gravity is taken into account using the “reduced gravity approach”.

#include "navier-stokes/centered.h"
#include "two-phase.h"
#include "navier-stokes/conserving.h"
#include "reduced.h"

The primary parameters are the wave steepness ak and the Reynolds number.

double ak = 0.35;
double RE = 40000.;
int LEVEL = 9;

The error on the components of the velocity field used for adaptive refinement.

double uemax = 0.005;

The density and viscosity ratios are those of air and water.

#define RATIO (1./850.)
#define MURATIO (17.4e-6/8.9e-4)

The wave number, fluid depth and acceleration of gravity are set to these values. T0 is the wave period.

#define k_  (2.*pi)
#define h_   0.5
#define g_   1.
#define T0  (k_/sqrt(g_*k_))

The program takes optional arguments which are the level of refinement, steepness and Reynolds numbers.

int main (int argc, char * argv[])
{
if (argc > 1)
LEVEL = atoi (argv);
if (argc > 2)
ak = atof(argv);

The domain is a cubic box centered on the origin and of length L0=1, periodic in the x-direction.

  origin (-L0/2, -L0/2);
periodic (right);

Here we set the densities and viscosities corresponding to the parameters above.

  rho1 = 1.;
rho2 = RATIO;
mu1 = 1.0/RE; //using wavelength as length scale
mu2 = 1.0/RE*MURATIO;
G.y = -g_;

When we use adaptive refinement, we start with a coarse mesh which will be refined as required when initialising the wave.

#if TREE
N = 32;
#else
N = 1 << LEVEL;
#endif
DT = 1e-2;
run();
}

## Initial conditions

We either restart (if a “restart” file exists), or initialise the wave using the third-order Stokes wave solution.

#include "stokes.h"

event init (i = 0)
{
if (!restore ("restart")) {

We need to make sure that fields are properly initialised before refinement below, otherwise a -catch exception will be triggered when debugging.

    event ("properties");

do {
fraction (f, wave(x,y));
foreach()
foreach_dimension()
u.x[] = u_x(x,y) * f[];
boundary ((scalar *){u});
}

On trees, we repeat this initialisation until mesh adaptation does not refine the mesh anymore.

#if TREE
(double[]){0.01,uemax,uemax,uemax}, LEVEL, 5).nf);
#else
while (0);
#endif
}
}

event profiles (t += T0/4.; t <= 2.5*T0) {
output_facets (f, stderr);
fprintf (stderr, "\n");
}

event logfile (i++)
{
double ke = 0., gpe = 0.;
foreach (reduction(+:ke) reduction(+:gpe)) {
double norm2 = 0.;
foreach_dimension()
norm2 += sq(u.x[]);
ke += norm2*f[]*dv();
gpe += y*f[]*dv();
}
printf ("%g %g %g\n", t/(k_/sqrt(g_*k_)), rho1*ke/2., rho1*g_*gpe + 0.125);
}

#if TREE
#endif