droplet bagmode breakup under hight density ratio .

[张祥]

Dear all,

I am a student want to simulate 2 dimension droplet bagmode breakup under high density ratio. and now I have completed low density ration simulation .such as rhol/rhog = 80

But when i use high density ratio,such as rhol/rhog = 100/1.2,mu suface become unstable . what should I do to solve this problem ? I guess this is because the momentum unconserved .But I dont know how to solve .

and my total program is below 

#define dimension 2

#include "axi.h"

#include "navier-stokes/centered.h"

#include "two-phase.h"

#include "navier-stokes/conserving.h"

#include "tension.h"

#include "view.h"

// ============para===================

#define RHOL 800.

#define MUL 50.e-3

#define SIGMA 0.03

/**Air properties*/

#define RHOG 10.

#define MUG 1.e-3

/**Drop diameter and stream velocity*/

#define DIAMETER 6.7e-3

#define USTREAM 3.

int minlevel = 6;  

int MAXLEVEL = 9;  

double maxruntime = HUGE;

u.n[left]  = dirichlet(USTREAM);

u.t[left]  = dirichlet(0);

p[left]    = neumann(0);

u.n[right] = neumann(0);

p[right]   = dirichlet(0);

// =========================================

int main (int argc, char * argv[]) {

  L0 = 10*DIAMETER;

  X0 = 0.;

  Y0 = 0.;

  // nondimensional para

  double Re = RHOG * USTREAM * DIAMETER / MUG;

  double We = RHOG * USTREAM * USTREAM * DIAMETER / SIGMA;

  double Oh = MUL / sqrt(RHOL * SIGMA * DIAMETER);

  double rho_ratio = RHOL / RHOG;

  double mu_ratio  = MUL / MUG;

  fprintf(stderr, "\n==== Dimensionless Numbers ====\n");

  fprintf(stderr, "Reynolds number (Re):       %.2f\n", Re);

  fprintf(stderr, "Weber number (We):          %.2f\n", We);

  fprintf(stderr, "Ohnesorge number (Oh):      %.4f\n", Oh);

  fprintf(stderr, "Density ratio (rho_l/rho_g): %.2f\n", rho_ratio);

  fprintf(stderr, "Viscosity ratio (mu_l/mu_g): %.2f\n", mu_ratio);

  fprintf(stderr, "===============================\n\n");

 

  double minDelta = L0/pow(2.0,MAXLEVEL);

  double sigmaTSLimit = sqrt(min(RHOL, RHOG)*pow(minDelta,3.0)/(SIGMA));

  double muTSLimit = min(MUL, MUG)*minDelta/SIGMA;

  DT = min(sigmaTSLimit,muTSLimit);

  CFL = 0.4;

  f.sigma = SIGMA;

  rho1 = RHOL;

  rho2 = RHOG;

  mu1 = MUL;

  mu2 = MUG;

  TOLERANCE = 1e-4;

  run();

 

}

event init (t = 0) {

  refine (sq(x-2*DIAMETER) + sq(y) - sq(0.6*DIAMETER) < 0 && sq(x - 2*DIAMETER) + sq(y) - sq(0.4*DIAMETER) > 0 && level < MAXLEVEL);

  fraction (f, sq(0.5*DIAMETER) - (sq(x - 2*DIAMETER) + sq(y)));

}

event adapt (i++) {

  double uemax = 5e-2;

  adapt_wavelet ({f,u}, (double[]){1e-3,uemax,uemax,uemax}, MAXLEVEL, minlevel);

}

event end (t= 0.06) {

}

//==================log=================

event logfile (i += 100 ,first) {

  if (i == 0){

    fprintf (ferr,

       "t dt mgp.i mgpf.i mgu.i perf.t perf.speed grid->tn\n");

  }

  fprintf (ferr,

     "%g %g %d %d %d"

     "%.2e %.2e %ld \n",

     t, dt, mgp.i, mgpf.i, mgu.i,

     perf.t, perf.speed, grid->tn);

  fflush (ferr);

}

//===================ppm=========================

//====================================

int frame = 0;

event movies ( t += 1.e-3 )

{

  view (fov = 12, width = 1920, height = 1080,tx = -0.6 , ty = -0.2);

  clear();

  draw_vof ("f", filled = 1, lw = 1, lc = {1, 1, 1}, fc = {0.6, 0.8, 1.0});

  cells();

  char name[256];

 

  sprintf(name, "ppm/snapshot-%d.ppm", frame++);

  save(name);

}

//===================field=========================

//====================================

// int num2 = 0;

// event output_slice (t += 1.e-3) {

//   char name[80];

//   sprintf(name, "field/%d.csv", num2++);  

//   FILE * fp = fopen(name, "w");

//   fprintf(fp, "x,y,f,u_x,u_y,p\n");

//   foreach() {

//     fprintf(fp, "%g,%g,%g,%g,%g,%g\n",

//             x, y, f[], u.x[], u.y[], p[]);

//   }

//   fclose(fp);

// }

//===================dump=================

//====================================

// event dumps ( t += 1.e-2) {

//   char fn[99];

//   sprintf (fn, "dump/dump-%g", t);

//   dump (file = fn);  // 储存当前网格、变量等

// }