Re: Bubble cannot burst with two-phase-compressible

[MANDEEP SAINI]

Hi,

I think you are very under resolved. You need many more points to accurately resolve the thin gap. Perhaps you can first try the testcase from the sandbox in your local machine to be sure you have the correct headers.

Best

Mandeep

On Sunday, September 3, 2023 at 2:01:23 PM UTC+2 wengeh...@gmail.com wrote:

Hello everyone,

     I'm trying to use the two-phase-compressible solver listed in Dr. Mandeep Saini's sandbox (http://basilisk.fr/sandbox/msaini/) to simulate a bubble burst problem.

     The simulation setup is simple: a high pressure single air bubble inside a water droplet; I want to simulate the bubble expanding and burst process and the following droplet dynamics.

      However, the simulation always stops when the bubble is about to burst (see the moves attached) no matter how large the initial pressure I imposed inside the bubble.

     Specifically, I find that when the bubble is about to burst, dtmax, and dt is approaching 0. Near the burst point position, the velocity uf.x[], uf.y[] and pressure p[] are becoming very large and the projection Poisson solver fails.

Now, I have no idea about this issue and any suggestion about this issue is really appreciated.

Thank you all very much.

Here are the codes:

//#include "axi.h"
#include "tpccontact.h"
#include "compressible-tension.h"
#include "view.h"

int MINLEVEL = 7;
int MAXLEVEL = 8;
double CFLac = 0.4;

double RhoL = 1., RhoG = 0.001;
double pg0;
double p0 = 1.;
double pinf = 0.5;
double tend = 1.;

double Rbub = 1.;
double lambda = 20.;
double Reynolds = 1000.*0.1;
double Web = 100./72.*1000.*0.1;

scalar keliq[];

uf.n[right] = neumann(0.);
p[right]    = neumann(0.);//dirichlet(pinf);
q.n[right]  = neumann(0.);
//f[right]    = dirichlet(1);
fE1[right] = neumann(0.);
fE2[right] = neumann(0.);

uf.n[top] = neumann(0.);
p[top]    = neumann(0.);//dirichlet(pinf);
q.n[top]  = neumann(0.);
//f[top]    = dirichlet(1.);
fE1[top] = neumann(0.);
fE2[top] = neumann(0.);

uf.n[left] = neumann(0.);
p[left]    = neumann(0.);//dirichlet(pinf);
q.n[left]  = neumann(0.);
//f[top]    = dirichlet(1.);
fE1[left] = neumann(0.);
fE2[left] = neumann(0.);


f[left] = 0;
uf.t[bottom] = neumann(0.);

vector h[];
double theta0 = 130.;
h.t[left] = contact_angle (theta0*pi/180.);


event stability(i++)
{
  double cspeed;
  foreach()
    {
      double fc = clamp (f[],0.,1.);
      double invgammaavg = fc/(gamma1 - 1.) + (1. - fc)/(gamma2 - 1.);
      double PIGAMMAavg = (fc*PI1*gamma1/(gamma1 - 1.) +(1. - fc)*PI2*gamma2/(gamma2 - 1.));
      double cspeedsq = (p[]*(invgammaavg + 1.) + PIGAMMAavg)/invgammaavg/(frho1[]+frho2[]);
      if (cspeedsq > 0.)
cspeed = sqrt(cspeedsq);
      else
cspeed = sqrt(gamma1*(pinf + PI1));
      double dtmaxac = CFLac*Delta/cspeed;
      dtmax = min(dtmax,dtmaxac);
    }
}

int main() {

  tend *= 0.915;
  pg0 = p0 + 1.*f.sigma/0.5;
  f.gradient = frho1.gradient = frho2.gradient = p.gradient = q.x.gradient = q.y.gradient = minmod2;
  f.sigma= 1.;
  mu1 = 1.;
  mu2 = mu1*0.01;
  gamma2 = 1.4;
  gamma1 = 7.14;
  PI1 = 1./sq(0.007)/7.14 - pinf;
  L0 = 10;;  
  //TOLERANCE = 1e-4;
  X0 = 0.;
  init_grid(1<<MINLEVEL);
  run();
  free_grid();
}



event init (i = 0) {
     
  if (!restore (file = "restart"))
    {
      vertex scalar phi_o[];
      vertex scalar phi_i[];
      //scalar f_trail_o[];
      //scalar f_trail_i[];
      #if 1
      foreach_vertex()
      phi_o[] = -(sq(x-7) + sq(y) - Rbub*Rbub*4);;
      fractions(phi_o, f_trail_o);

      foreach_vertex()
      phi_i[] = (sq(x-6.2) + sq(y) - Rbub*Rbub);
      fractions(phi_i, f_trail_i);


      //serve as a marker
      foreach()
         f[] = 2.;

      foreach(){
      if(f_trail_i[] < 1. || f_trail_i[] == 0.)
      f[] = f_trail_i[];
      if(f_trail_o[] < 1. || f_trail_o[] == 0.)
      f[] = f_trail_o[];
      }

      foreach(){
      if(f[] == 2.)
      f[] = 1.;
      }
      boundary({f});
      #endif
     
      #if 0
      vertex scalar phi[];
      foreach_vertex() {
phi[] = sq(x) + sq(y) - Rbub*Rbub;
      }
      boundary ({phi});
      fractions (phi, f);
      boundary({f});
      #endif
      #if 0
      foreach() {
frho1[]  = f[]*rhoL;
frho2[]  = (1. - f[])*rhoG;

double pL = pinf;//*(1. - Rbub/sqrt(sq(x) + sq(y))) + (pg0 - 2*f.sigma/Rbub/Rbub)*Rbub/sqrt(sq(x) + sq(y));
double pg = pg0;

fE1[]   = f[]*(pL/(gamma1 - 1.) + PI1*gamma1/(gamma1 - 1.));
fE2[]   = (1.-f_trail_i[])*10*pg/(gamma2 - 1.)+ (1.-f_trail_o[])*pL/(gamma2 - 1.);


p[] = (fE1[] + fE2[] - PIGAMMAavg)/invgammaavg;
q.x[] = 0.;
q.y[] = 0.;
      }
      #endif
     double u0 = 0.;
     foreach(){
          frho1[] = RhoL*f[];
          frho2[] = RhoG*(1. - f_trail_i[])*10 + RhoG*(1. - f_trail_o[]);
          p[] = 10*pg0*f[] + p0*(1. - f_trail_o[]) + 200*p0*(1. - f_trail_i[]);
          q.x[] = 0.;
          q.y[] = 0.;
          fE1[] = f[]*(p[] + gamma1*PI1)/(gamma1 - 1.) + 0.5*frho1[]*pow(u0,2);
          fE2[] = (1. - f[])*(p[] + gamma2*PI2)/(gamma2 - 1.) + 0.5*frho2[]*pow(u0,2);  
    }
     
      boundary ((scalar *){q,frho1,frho2,p,fE1,fE2});
      dump("init");
    }
}


#if 0
scalar unorm[],fscalar[],pscalar[];
event adapt (i++) {  
  foreach(){
    unorm[] = norm(q)/(frho1[]+frho2[]);
    fscalar[] = f[];
    pscalar[] = p[];
  }
  boundary ((scalar *){fscalar,unorm});
 
  adapt_wavelet((scalar *){fscalar,pscalar,unorm},(double[]){0.001,0.001,0.002},maxlevel = MAXLEVEL);
  //adapt_wavelet((scalar *){fscalar},(double[]){0.001},maxlevel = MAXLEVEL);
}
#endif

event movies (t += 0.001; t<= 0.1) {
//event movies (i++; t<= tend) {
  //view(fov = 1., ty = -0.02, quat = {0,0,-cos(pi/4.),cos(pi/4.)}, width = 3960, height = 3960);
  view(tx = -0.5);
  char name2[80];
  sprintf (name2,"T0.52.mp4");
  foreach(){
    f[] = clamp(f[],0,1);
  }
  boundary ({f});
  clear();
  cells();
  //box (notics = true);
  draw_vof("f", lw = 2);
  //draw_vof("f", filled = 1, fc = {0.9, 0.4, 0.2});
  scalar Ek[];
  foreach(){
  foreach_dimension(){
     Ek[] += sq(q.x[]);
  }
  }
  draw_vof("f", filled = 1, fc = {0.9, 0.4, 0.2});
  squares ("Ek", min = -10000., max = 10000., linear = false, map = cool_warm);
  mirror ({0,1}) {
  //cells();
  //box (notics = true);
  draw_vof("f", lw = 2);
  //draw_vof("f", filled = 1, fc = {0.9, 0.4, 0.2});
  squares ("p", min = -100., max = 1000., linear = false, map = cool_warm);
  }
  save(name2);
}