sandbox/yonghui/smalltest/mcvisdrop.c
water drop in viscous environment by momentum conservation methodes
Based on rising.c,
We replace two liquids with each other in order to simulate falling drops of water in a viscous environment.
The main part of the code has not been modified, comparisons will be neglected here. We simply assume the results is good.
#include "axi.h"
#include "navier-stokes/centered.h"
#include "two-phase.h"
#include "navier-stokes/conserving.h"//actually in this simple case it's useless ??
#include "tension.h"
# define LEVEL 8
BC
u.t[right] = dirichlet(0.);
u.t[left] = dirichlet(0.);
u.n[right] = dirichlet(1.45);
u.n[left] = dirichlet(1.45);
p[right] = dirichlet(0);
p[left] = neumann(0);
pf[right] = dirichlet(0);
pf[left] = neumann(0);
uf.n[bottom] = 0.;
uf.n[top] =0.;
Main
We need to set the calculation domain large enough to avoid the drop contact boundaries
int main() {
size (1.5);
init_grid (128);
rho1 = 1000., mu1 = 10.;
rho2 = 100., mu2 = 1., f.sigma = 24.5;
TOLERANCE = 1e-4;
run();
}
initial & acceleration
event init (t = 0) {
mask (y > 0.6 ? top : none);
refine (sq(x - 0.5) + sq(y) - sq(0.25) < 0. && level < LEVEL);
fraction (f, -1.*(sq(x - 0.5) + sq(y) - sq(0.25)));
output_ppm (f, linear = true, n=512, file = "init.png");
}
event acceleration (i++) {
face vector av = a;
foreach_face(x)
av.x[] -= 0.98;
boundary({p});
}
some output event.
I upload a temporary results, as for now , we get a same final shape as the one with level = 12 & solved by [visdrop.c]
For the moment, it’s good, I’ll upload the necessary studies later
event outputvideo (i += 20 ) {
output_ppm (f, linear = true, n=1024,file = "f.mp4");
output_ppm (u.x, linear = true, n=1024, file = "ux.mp4");
output_ppm (p, linear = true, n=1024, file = "p.mp4");
fprintf(stderr,"%d %g\n", i,t);
}
event interface (t = 7) {
output_facets (f, stderr);
}
int k=0;
event velo (t += 0.5) {
if (t > 2.){
k += 1;
char name[40];
char velofield[40];
sprintf (velofield, "ve-%d", k);
sprintf (name, "dump-%d", k);
dump (file = name);
output_field((scalar *){u}, fopen (velofield, "w"));
}
}
event adapt (i++) {
adapt_wavelet ({f,u}, (double[]){5e-4,1e-3,1e-3}, LEVEL,6);
}
I will update the one can be done by ploting code later