sandbox/M1EMN/BASIC/disperse.c
Wave equation with dispersion due to capillarity
Explicit resolution of linearized Saint Venant with dispersion \displaystyle \frac{\partial h}{\partial t}+ \frac{\partial hu}{\partial x}=0 \displaystyle \frac{\partial hu}{\partial t}= - \frac{\partial h}{\partial x} + \sigma \frac{\partial^3 h}{\partial x^3} Note that there are no nonlinear terms \frac{\partial hu^2}{\partial x} in the total derivative d_t hu of momentum for the moment (hu is the flux, h is the height of water, \sigma the surface tension parameter).
//#include "grid/cartesian1D.h"
#include "grid/multigrid1D.h"
#include "run.h"
scalar h[];
face vector hu[];
double dt, sigma;
FILE * out;
Outflow boundary conditions.
hu.n[left] = neumann(0);
h[left] = neumann(0);
hu.n[right] = neumann(0);
h[right] = neumann(0);
Parameters, siez of domain, origin, number of points, time step
int main() {
L0 = 50.;
X0 = -5;
N = 512;
DT = 5e-5;
Loop on surface tension \sigma. Results are written in a file indexed by surface tension.
for (sigma = 0.; sigma <= 0.1; sigma += 0.1) {
char name[80];
sprintf (name, "out-%g", sigma);
out = fopen (name, "w");
run();
}
}
The initial field is a linear wave moving to the right.
event init (t = 0) {
#define EPS 0.01
foreach()
h[] = 1. + EPS*exp(-x*x);
foreach_face()
hu.x[] = EPS*exp(-x*x);
boundary ({h,hu});
}
Output fields for plots.
event printdata (t += .1; t <= 25) {
foreach()
fprintf (out, "%g %g %g %g\n", x, h[], hu.x[], t);
fprintf (out, "\n");
}
Integration.
event integration (i++) {
double dt = DT;
scalar u[];
foreach_face() {
u[] = 2.*hu.x[]/(h[] + h[-1,0]);
double un = fabs(u[]);
if (un > 0. && CFL*Delta/un < dt)
dt = CFL*Delta/un;
}
dt = dtnext ( dt);
Compute the curvature K = \frac{\partial^2 h}{\partial x^2}
Curvature contribution to pressure and hydro \displaystyle \frac{\partial hu}{\partial t}= - \frac{\partial h}{\partial x} + \sigma \frac{\partial K}{\partial x}
foreach_face()
hu.x[] += dt*(- (h[] - h[-1,0]) + sigma*(K[] - K[-1,0]))/Delta;
Variation of h \displaystyle \frac{\partial h}{\partial t}+ \frac{\partial hu}{\partial x}=0
scalar dh = K; // we just reuse K as dh
foreach()
dh[] = ((u[] < 0. ? h[] : h[-1,0])*u[] -
(u[1,0] > 0. ? h[] : h[1,0])*u[1,0])/Delta;
foreach()
h[] += dt*dh[];
boundary ({h,hu});
}
Run
qcc -g -O3 -o disperse disperse.c -lm
./disperse
or with make
make disperse.tst;make disperse/plots
make disperse.c.html ; open disperse.c.html
Results
The abscissa is x the ordonate is time t.
First we verify the propagation of a pulse which does not change during the propagation
set output 'dalembert.png'
set pm3d map
set palette gray negative
unset colorbox
set tmargin at screen 0.95
set bmargin at screen 0.15
set rmargin at screen 0.95
set lmargin at screen 0.15
set xlabel "x"
set ylabel "t"
unset key
splot 'out-0' u 1:4:2
Second we observe that dispersion distroys the signal and breaks it in small waves traveling faster. Note the initial wave generates a small one running left and reflecting on the border.
set output 'disperse.png'
set pm3d map
set palette gray negative
unset colorbox
set tmargin at screen 0.95
set bmargin at screen 0.15
set rmargin at screen 0.95
set lmargin at screen 0.15
set xlabel "x"
set ylabel "t"
unset key
splot 'out-0.1' u 1:4:2
Next steps
- couple with Shallow Water solvers
- compare to the case with \sigma negative (see mascaret)
ready for new site 09/05/19