sandbox/jventre/Stenose/MaskJ.c

    Resolution de Poiseuille dans un cylindre

    Equation à résoudre :

    \displaystyle 0 = -\frac{\partial p}{\partial x} + \frac{1}{r} \frac{\partial }{\partial r} r \frac{\partial u}{\partial r}

    #include "axi.h"
    #include "navier-stokes/centered.h"
    #include "two-phase.h"
    #include <sys/stat.h>
    
    int Re ;  // Reynolds
    double dR = 0.; // degrée d'obstruction
    double xst = 2.5 ; // position de la sténose
    scalar un[];
    
    double stenose(double x){
    variable ‘xb’ set but not used [-Wunused-but-set-variable]
        double xb,xl;
        xl=1.;
        xb=xst-.5*xl;
        //return (fabs((x-xb)/xl-.5)<.5?  (1+cos(pi+2*pi*(x-xb)/xl))/2.   : 0 );
        return  exp(-sq((x-xst)/xl));
    }
    
    int main() {
      
      TOLERANCE = 1e-6;
      L0=25;
      Re= 100;
      // two-phase gives f=1 to 1 and f=0 to phase 2  
       rho1 = 1, rho2 = 1.;
       mu1 = 1./Re ;
       mu2 = 0.;
     
      N = 512 ; 
    
      system("rm p_final_Re_100.csv");
      system("rm u_final_Re_100.csv");
      run();  
    
      // Re = 100;
      // L0 = 18;
      // system("rm p_final_Re_100.csv");
      // system("rm u_final_Re_100.csv");
      // mu1 = 1./Re ;
      // run(); 
    
      // Re = 200;
      // L0=25;
      // system("rm p_final_Re_200.csv");
      // system("rm u_final_Re_200.csv");
      // mu1 = 1./Re ;
      // run(); 
    
      // Re = 300;
      // L0=30;
      // N=1024;
      // system("rm p_final_Re_300.csv");
      // system("rm u_final_Re_300.csv");
      // mu1 = 1./Re ;
      // run(); 
    
      //  Re = 500;
      //  L0=35;
      //  N=1024;
      // system("rm p_final_Re_500.csv");
      // system("rm u_final_Re_500.csv");
      // mu1 = 1./Re ;
      // run(); 
    
    }

    Conditions aux limites

    • vitesse nulle dans la phase du haut
    • dirichlet 1 sur la vitesse en entrée
    • dirichlet sur la pression à l’entrée et à la sortie (8*L/Re en entrée, 0 en sortie pour avoir un gradient = 8/Re)
    u.t[top] = dirichlet(0);
    u.n[top] = dirichlet(0);
    
    u.n[left] = ( y < 1 ? dirichlet(1.):dirichlet(0.));
    p[left] = neumann(0.);
    pf[left]   = neumann(0.);
    
    u.n[right] =  ( y < 1 ? neumann(0.) : dirichlet(0.) ) ;
    p[right]   = dirichlet(0.);
    pf[right]  = dirichlet(0.);

    Définition de l’interface

    event init (t = 0) {
     
      // géométrie de l'interface
      mask (y >  1 - dR*stenose(x)? top :    none);
        
      // on sauve l'interface 
      FILE * fp2 = fopen ("interface.csv", "a");
      output_facets (f, fp2);
      fprintf(fp2,"\n\n");
      fclose (fp2);
    
      // initialisation de un pour la convergence
      foreach()
        un[] = u.x[];
    }

    Test de la convergence sur la vitesse

    event conv (t += 1) {
        double du = change (u.x, un);
        fprintf(stdout,"t= %g %g \n",t, du);
        if (i > 0 && du < 1.0e-4)
            return 1; /* stop */
    }
    
    event bound (i++,last){
    
      // on doit redéfinir à chaque itération l'interface
      fraction (f, (1. - dR *exp(-(x-xst)*(x-xst)))- y  );
    
      // a chaque itération on impose la vitesse nulle en haut de l'interface 
      foreach() {
        if (f[]>1e-6 && f[]<1.-1e-6){
         u.x[]=u.y[]=p[]=0.0;
        } else {
        u.x[] = u.x[]*(f[]);
         u.y[] = u.y[]*(f[]);
         p[] = p[]*(f[]);
       }
      }
    
    }

    Enregistrement des vitesses/pressions finales

    event boucle(t+=1;t<=300){
      // printf(" t = %g \n", t);
    }
    
    
    event sauve(t=end)
    {
      char filename[80];
      sprintf (filename, "p_final_Re_%d.csv", Re);
      FILE * fp = fopen (filename, "a");
      output_field ({p}, fp, linear = true);
      fprintf(fp,"\n\n");
      fclose (fp);
    
      char filename2[80];
      sprintf (filename2, "u_final_Re_%d.csv", Re);
      FILE * fp2 = fopen (filename2, "a");
      output_field ({u.x}, fp2, linear = true);
      fprintf(fp2,"\n\n");
      fclose (fp2);
    
    }

    Run

    pour lancer

    
    make entreepoiseuille.tst
    make entreepoiseuille/plots
    make entreepoiseuille.c.html

    Results

    Exemples

    
    stats 'u_final_Re_100.csv' u 1
    
    plot 'u_final_Re_100.csv' i STATS_blocks-2 us 1:($2==0?$3:NaN) w linespoints title 'Re=100'
    
    set xlabel 'y'
    set ylabel 'u'
    
    rep
    velocity (script)

    velocity (script)

    
    stats 'p_final_Re_100.csv' u 1
    
    p'p_final_Re_100.csv' u 1:($2==0?$3:NaN) title 'entry effect'
    rep 8.*25/100.-8./100.*x, title "poiseuille pressure"
    
    set xlabel 'y'
    set ylabel 'p'
    
    rep
    pressure (script)

    pressure (script)

    bibliography