sandbox/yonghui/vortex/rpi.c
Rayleigh-Plateau instability
we want to study the growthrate of a axisymmetric Rayleigh-Plateau instability in the linear regime. The dimensionless NS equation can be written in the form \displaystyle \bar \rho \frac{D \bar u}{D \bar t } = - \nabla \bar P + Oh \sqrt{ \frac{R_0}{\lambda} } \bar \mu \Delta \bar u + \bar \gamma \kappa \delta_s n with he Ohnesorge number Oh = \frac{\mu_1}{\sqrt {\rho_1 \gamma R_0}} as the reverse Re number, here we proposed the caracteristic length L_0 = \lambda.
The values f_{sim} of this simulation can be scaled to f_d as u_{sim}= \sqrt{\frac{\lambda}{R_0}} u_d, l_{sim} = \frac{R_0}{\lambda} l_d, t_{sim} = \sqrt{\frac{R_0^3}{\lambda^3}} t_d if we propose L_0 = R_0
#include <gsl/gsl_fit.h> //need gsl lib
#pragma autolink -lgsl -lgslcblas
//#include "grid/multigrid.h"
#include "axi.h"
#include "navier-stokes/centered.h"
#include "two-phase.h"
#include "tension.h"
#include "view.h"
// AMR switch
#define STATIC 1
#define ADAPT 0
#define LVMAX 9
#define LL 9. // $\lambda$ / $R_0$
#define Oh 1e-3
#define mur 1. // ratio of mu2 / mu1
#define rhor 1. // ratio of rho2 / rho
#define tend 0.55 // pinchoff moment ~ tend 1.5
FILE *fp, *fp2, *fp3;main event
The simulation is periodic at left and right.
We set the initial perturbation amplitude as \epsilon = 0.01 \frac{\Delta x}{R_0} as a optimal choice of error study.
double hh[99999],tt[99999],epsilon;
int kkk,jjj;
int main() {
fp=fopen("fposition","w");
fp2=fopen("lineardata","w");
fp3=fopen("slope","w");
fprintf (fp3,"#alpha slope N oh\n");
#if STATIC
N = 32;
#else
N=pow(2,LVMAX);
#endif
periodic(right);
epsilon = 0.01*LL/(pow(2,LVMAX));
rho1 = 1., mu1 = Oh*sqrt(1. / LL);
rho2 = rhor, mu2 = mu1*mur;
f.sigma=1.;
kkk=0;
jjj =0;
fprintf(stderr,"now N= %d,test kkk %d jjj%d\n" ,N,kkk, jjj);
run();
}boundary conditions
classic condition of no penetrating wall at top
u.n[top] = dirichlet(0.);
p[top] = neumann(0.);init event
The two phases are seperated by the initial interface r = r_0 (1+\epsilon cos(ks)).
With the fraction field f, the inner liquid noted liquid 1 is presented by f = 1, and the external fluid noted liquid 2 is presented by f = 0.
We assume the error is related to the Taylor expension of the radial velocity at order \alpha, as presented by equation: \displaystyle L = L_0 - int\left[\frac{1}{\alpha}\,log_2\left(\frac{{\frac{\partial^{\alpha}u}{\partial r^{\alpha}}}(r=r_0)}{\frac{\partial^{\alpha}u}{\partial r^{\alpha}}(r)}\right)\right] The theorytic profil of the radial velocity perturbation of RPI in the linear regime is under the form of Bessel function, such that u_{r1}' = C_1 I_1(kr) and u_{r2}' = C_2 K_1(kr).
The grid redistribution with \alpha = 2 is implemented with the switch “STATIC”.
event init (t = 0) {
#if STATIC
double position[] = {0.7,1.13,1.78,2.66,3.70,4.85,6.05};
for(int GG= 1 ; GG <= 6 ; GG++){
refine( (y >0.03 ) && (y < position[GG]/(2.*M_PI) )
&& (level <(LVMAX+ 1 -GG)) ); }
refine ((y < 0.03 ) && (level <(LVMAX-1)));
#endif
fraction (f, 1./LL*(1.+ epsilon*cos(2.*M_PI*x)) -y );
//general fraction field view
view(fov=10.,tx=-0.5,ty=-0.2,width=1280,height=640);
draw_vof("f", lw = 1.5);
cells();
save ("fim.png");
fprintf(stderr,"now cell = %ld\n" , grid->tn);
}output
The movies of interface and pressure variation.
event movies (t += tend/100; t <= tend) {
char legend[2000]; //time
sprintf(legend, "t = %0.2g, \n LL = %g", t, LL);
//general view
view(fov=10.,tx=-0.5,ty=-0.2,width=1280,height=640);
squares("p");
draw_vof("f", lw = 1.5);
draw_string(legend, 1, size = 30., lw = 2.);
save ("pressure.mp4");
//movie 2 for detatch details
//view(fov=3.,tx=-0.6,ty=-0.05,width=1280,height=640);
//box();
//squares("p",min=0,max=1);
//draw_vof("f", lw = 1.5);
//draw_string(legend, 1, size = 30., lw = 2.);
//save ("detail.mp4");
}Here we save the position of maximum radius. A control parameter “scale” is used to verify if the simulation is still in the linear regime, that \frac{\eta_1 -\eta_2}{r_0} = \frac{r_{max} + r_{min}}{r_0} -2 < 0.02, with \eta the radius perturbation at max/min radius. The simulation will stoped if the linear condition is not satisfied.
If the small perturbation theory is satisfied, the perturbation growth with time \eta(r,t) = \epsilon(r) e^{st}, with s the growthrate.
We use a LLS as a linear regression method to calculate the growthrate \displaystyle ln \, (\frac{\lambda}{r_0} \eta) = C_{st} + \bar{s} ( t (\frac{\lambda}{r_0})^{\frac{3}{2}})
event slope (i += 1; t <= tend) {
scalar pos2[];
position(f,pos2,{0,1});
double scale = (statsf(pos2).max + statsf(pos2).min)*LL-2.;
double hh2 =statsf(pos2).max*LL-1.;
double tt2 =t*pow(LL, 1.5);
fprintf(fp," %g %g %g %g %g %g %d\n", tt2, hh2, scale,
(statsf(pos2).max-1./LL)*N, epsilon, L0, N);
double endeta;
if ((t*pow(LL, 1.5) > 8.) && (fabs(scqle) < 1e-2) ){
hh[kkk] = log(hh2);
tt[kkk] = tt2;
fprintf(fp2,"%d %d %g %g %g %g\n", N, kkk, t, hh2, tt2,scale);
kkk += 1;
endeta = (statsf(pos2).max-1./LL)*N;
}
if( ((t*pow(LL, 1.5) > 8.) && (fabs(scale) > 1e-2)) || (t*pow(LL, 1.5) > 14.)){
double c0, c1, cov00, cov01, cov11, sumsq;
gsl_fit_linear (tt, 1, hh, 1, kkk-2,
&c0, &c1, &cov00, &cov01, &cov11,
&sumsq);
fprintf (stderr,"#best fit: ln(h) = %g + %g t\n", c0, c1);
fprintf (stderr,"#covariance matrix:\n");
fprintf (stderr,"#[ %g, %g\n #%g, %g]\n",
cov00, cov01, cov01, cov11);
fprintf (stderr,"#sumsq = %g\n", sumsq);
fprintf (fp3,"%g %g %d %g %g\n", epsilon, L0, N, c1, endeta );
return 1;
}
}Adapt wavelet function.
As for now, the non-static AMR seems worse than the static AMR in the linear regime.
