Bubble in a stagnation point flow.

    We first create the axi-symmetric stagnation point flow. This is done by imposing the analytical boundary conditions (inflow on the top and outflow on the right). After a transient, the flow converges to the analytical solution. The stationary flow is then used to study the deformations and breakup of an initially spherical bubble which lies at the stagnation point.

    The velocity field of a stagnation point flow is (using Basilisk notation, y is the radial direction and x the axial direction.): \displaystyle \bm{u} = -Ex \bm{e_x} +\frac{1}{2}E y \bm{e_y}

    In what follows, the flow enters radially and outs along the axis of symmetry, so that E <0.

    //#include "grid/quadtree.h"
    #include "axi.h"
    #include "navier-stokes/centered.h"
    #include "view.h"
    #include "two-phase.h"
    #include "navier-stokes/conserving.h"
    #include "tension.h"
    #include "navier-stokes/perfs.h"
    #include "maxruntime.h"
    //Physical ratios
    #define RHOR 850.0 //density ratio
    #define MUR 25.0 //kinematic viscosity ratio; real value 55; to be consistent with turbubble simulations
    #define HALFBOX 5.
    double amp=-1., R0=1., SIGMA=1.;//velocity amplitude; bubble radius; surface tension
    double MAXTIME = 5.;

    Boundary conditions are chosen so that they match the extract solution for the flow at a stagnation point.

    //boundary conditions
    u.n[left] = dirichlet(0.);
    u.n[right] = neumann(-amp);
    u.n[top]    = dirichlet(1./2.*amp*2.*HALFBOX);
    u.n[bottom] = dirichlet(0.);
    p[top] = neumann(-1./2.*sq(amp)*HALFBOX);
    pf[top] = neumann(-1./2.*sq(amp)*HALFBOX);
    p[bottom] = neumann(0.);
    pf[bottom] = neumann(0.);
    p[left] = neumann(0.);
    pf[left] = neumann(0.);
    p[right] = dirichlet(-1./8.*sq(amp*y));
    pf[right] = dirichlet(-1./8.*sq(amp*y));
    int main(int argc, char * argv[]){
        if (argc>1){
    	amp = atof(argv[1]);//typical shear
    	R0 = atof(argv[2]);//in practice, R0 is kept constant at 1. It is the length reference
    	SIGMA = atof(argv[3]);
    	MINLEVEL = atoi(argv[4]);
    	MAXLEVEL = atoi(argv[5]);
    	MAXTIME = atof(argv[6]);
        //Physical properties
        rho1 = 1.0;//liquid is the reference
        rho2 = rho1/RHOR;
        mu1 = 0.01;
        mu2 = mu1/MUR;
        f.sigma = SIGMA;
    event init(i = 0){

    The same simulation can be used to create the stagnation point flow and to inject a bubble in a previously computed flow. To inject a bubble in a flow, be sure to include in the folder a dump file called ‘restart’.

        if (!restore("restart") && !restore("dump")){
    	    f[] = 1.0;
        else if (restore(file="restart")){//bubble injection
    	fraction(f, sq(x) + sq(y) - sq(R0));
        p.nodump=false;//just added

    We adapt on both the velocity field and the position of the interface.

    event adapt(i++){//grid adaptation
        double uemax = 1.e-3;
        double femax = 1.e-2;
        adapt_wavelet((scalar *){f, u}, (double[]){femax, uemax, uemax}, MAXLEVEL, MINLEVEL);

    To compute the spherical harmonics decomposition, we extract the coordinates of the interface. This is very low level and should be improved in the future.

    event getInterface(i++){
        char fixname[100];//x
        char fiyname[100];//y
        char fikname[100];//curvature
        sprintf(fixname, "fix_%d.dat", pid() );
        sprintf(fiyname, "fiy_%d.dat", pid() );
        sprintf(fikname, "fik_%d.dat", pid() );
        FILE * fix = fopen(fixname, "a");
        FILE * fiy = fopen(fiyname, "a");
        FILE * fik = fopen(fikname, "a");
        fprintf (fix, "%g ", t);
        fprintf (fiy, "%g ", t);
        fprintf (fik, "%g ", t);
        scalar curve[]; //curvature
        curvature(f, curve);
        //position with the heigh function
        scalar xpos[];
        scalar ypos[];
        position(f, xpos, {1,0});
        position(f, ypos, {0,1});
          if (interfacial(point,f)){
    	fprintf (fix, "%g ", xpos[]);
            fprintf (fiy, "%g ", ypos[]);
    	fprintf (fik, "%g ", curve[]);
        fprintf (fix, "\n");
        fprintf (fiy, "\n");
        fprintf (fik, "\n");
        return 0;

    We generate a movie of the simulation. The bubble will be in white, the * surrounding flow in black.

    event movie(t+=0.05){
      draw_vof("f", filled = 1);
      mirror({0, 1}){
        draw_vof("f", filled = 1);

    We output a dump file at the end of the simulation.

    event end(t = MAXTIME){
        printf("%g", t);