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){
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
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