sandbox/M1EMN/Exemples/viscolsqrt.c

    Viscous Collapse:

    Effondrement visqueux: il s’agit de l’effondrement d’un tas 1D initialement rectangulaire et disposé sur une pente Z(x). Comme l’écoulement est très visqueux, l’inertie est vite négligeable: on a vite un équilibre entre le terme de pente (de gravité) qui fait tomber et le frottement visqueux qui freine. Cet écoulement pourrait modéliser de la lave s’écoulant le long d’un volcan (il est générique en géophysique, le terme à changer est celui du frottement, nous verrons sur les autres pages le cas de bingham, fluide à seuil).

    Les équations de Saint Venant \displaystyle \frac{\partial h} {\partial t} + \frac{\partial (Q)}{\partial x}=0 \displaystyle \frac{\partial Q } {\partial t} + \frac{\partial }{\partial x}(\Gamma\frac{Q^2}{ h } + \frac{g}{2} (h_{\;}^2))= -gZ' h - Cf \frac{Q}{h_{\;}^2} avec \Gamma=1 (au lieu de (5/4)) et C_f=3\nu sont résolues pour la partie non visqueuse dans:

    #include "grid/cartesian1D.h"
#include "saint-venant.h"

    On se donne une masse initiale A=\int_{0}^{x_1} h(x,0) dx sur une pente \alpha

    double Cf,h0,alpha,tmax;
#define MAXLEVEL 13
#define MINLEVEL 11

    Initialisation d’un tas de masse 1

    event init (i = 0)
{
  foreach() {
    zb[] = -(x-X0)*alpha;
     h[] = h0*(x>0)*(x<=1./h0);
  }
}

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

event friction (i++) {
  // Poiseuille friction (implicit scheme)
  foreach()
  {  	 double ff = h[] < dry ? HUGE : (1. + Cf*dt/h[]/h[]);
    u.x[] /= ff;}
  boundary ({u.x});
}

 event plot (t=0.000001;t < tmax;t+=50) {	
 	 FILE * fp  = fopen ("sol.dat", "w");
//a=.584;p[][]'sol.dat'u ($1/$3**(1./3)):($2*$3**(1./3)),'sol13.dat' u ($1/$3**(1./3)):($2*$3**(1./3)),a*sqrt(x)*(x<(3./2./a)**(2./3))
   fprintf (stderr," %g \n", t);
   printf ("a=.584;p[:][0:.25] '-' u 1:2 t'calc' w l,''u 1:(a*( sqrt($1/$3))*($1<((9/4/a/a*$3)**(1./3))))t'selfs' w l \n");
     for(double x=X0;x<L0+X0;x+=L0/N){
       printf (" %g %g %g \n",x , interpolate (h, x,t),t);
       fprintf (fp,"%g %g %g \n",x , interpolate (h, x,t),t);}
       fclose(fp);
       // foreach()
       // printf (" %g %g %g \n",x ,h[],t);
       printf ("e\n");
}
//event adapt (i++) { adapt_wavelet ({h}, (double[]){1e-3}, MAXLEVEL, MINLEVEL); }
//event adapt (i++) {
//  astats s = adapt_wavelet ({h}, (double[]){1e-5}, MAXLEVEL, MINLEVEL);
 // fprintf (stderr, "# refined %d cells, coarsened %d cells\n", s.nf, s.nc);
//}

int main() 
{
  L0 = 30;
  X0 = -2.;
  G = 1;
  N = 512;
  DT=0.002;
  h0 = .25;
  Cf = .25; // valeur de 3 nu
  alpha=.25;
  tmax = 1500;

  run();
}

    Run

    Ensuite on compile et on lance:

    qcc -fopenmp  -g -O2 -DTRASH=1 -Wall  viscolsqrt.c   -o viscolsqrt
./viscolsqrt | gnuplot

    avec make

     make viscolsqrt.tst
 make viscolsqrt/plots
 make viscolsqrt.c.html
 
 source c2html.sh viscolsqrt

    Si on veut générer la figure

    ./viscolsqrt > v.out
echo "set term png;set output 'slope.png';set multiplot; set nokey ; set xrange[-5:25] ; set yrange [:1];" > dmp
cat v.out >> dmp
mv dmp v.out 
cat v.out | gnuplot

    ce qui produit la comparaison entre la solution semblable et le calcul: collapse

    Les fronts en fonction du temps sortent de out

    set xlabel "x"
set ylabel "h(x,t)"
p[][0:.5]'out' u 1:2 t'h(x,t)' w l
    heap as a function of time (script)

    heap as a function of time (script)

    Comparaison avec la solution analytique autosemblable F(\eta)

    Pour les pentes fortes telles que |Z'|\gg \partial_x h, on peut simplier les équations car on n’a plus que l’équilibre entre la pente et le frottement: alors Q= -gZ' h^3/(3 \nu) et l’équation de la masse devient \displaystyle \frac{\partial h} {\partial t} - \frac{gZ'}{3\nu} \frac{\partial }{\partial x} h_{\;}^3 =0, l’équation est de la forme
    \displaystyle \frac{\partial h} {\partial t} -k h_{\;}^2 \frac{\partial h}{\partial x} =0

    Huppert en trouve une solution par les caractéristiques, ici, nous nous proposons de trouver la soution sous forme auto semblable.

    Par changement d’échelle, l’invariance par dilatation donne H/T=H^3/X, la masse conservée donne HX=1, ce qui permet de trouver H=T_{\;}^{-1/3} et X=T_{\;}^{1/3} donc une solution auto semblable de la forme t_{\;}^{-1/3}F(x/t_{\;}^{1/3}) convient, par subsitution \displaystyle -\eta F'(\eta )-F(\eta )+3 k F(\eta )^2 F'(\eta )=0 avec F(0)=0 on a la forme implicite \eta=k (F(\eta)^2) et donc la solution est F(\eta)=\sqrt{(\eta)/k}, la solution du problème est simplement pour 0<x<x_f \displaystyle h =t_{\;}^{-1/3} \sqrt{(x)t_{\;}^{-1/3})/k}= \sqrt{\frac{x}{kt} }

    On s’est donné une masse initiale A=\int_{0}^{x_1} h(x,0) dx qui se déplace entre x_f t, position du front et 0 le front arrière qui ne bouge pas, la masse est donc ultérieurement \displaystyle \int_0^{x_f} \sqrt{\frac{x}{kt}} dx = (2/3) \sqrt{x_f^3/(kt)} qui est A par conservation de la masse soit x_f = (\frac{9 A_{\;}^2kt}4)^{1/3}

    a=.584;
set key left
p[-0.25:2.25][0:1]'sol.dat'u ($1/$3**(1./3)):($2*$3**(1./3)) w l ,a*sqrt(x)*(x<(3./2./a)**(2./3)) t'anal'
    selfsimilar solution (script)

    selfsimilar solution (script)

    plus on augmente la résolution, moins la partie gauche va bouger: ici on voit que le bloc a glissé, ce qui est une erreur numérique. Ce glissement est dû à l’étape visqueuse

    Links

    Bibliographie

    OK v2