sandbox/M1EMN/Exemples/embouch.c

    Propagation d’onde de marée dans une embouchure

    problème

    “La déformation des ondes de marées en pleine mer est assez bien connue. En grand profondeur, cette dynamique peut raisonnablement être modélisée par les équations linéarisées de Saint-Venant. Cependant, lors de son entrée sur le plateau continental et encore plus à l’intérieur des estuaires, l’onde se propage en se déformant.” Il s’agit du projet proposé par wiki.shf-hydro.fr. “A partir de l’ensemble de ces informations, il s’agira de proposer un modèle simplifié d’estuaire. Il s’agit ensuite de mettre clairement en évidence les effets spécifiques de chaque terme dans la déformation du signal temporel de marée.

    Topographie sans dimension et paramètres

    La carte https://www.geoportail.gouv.fr/donnees/carte-littorale nous donne une longueur L=89 000 m, avec profondeur h_0=4 m. A la louche les 4 stations sont réparties en 0 L/4 L/2 3L/4 et L: 1 (Pointe de la Grave), 2 (Richards), 3 (Lamena), 4 (Pauillac), 4 (Bordeaux):

    La vitesse de propagation des ondes est donc
    c_0=\sqrt(9.81*4)= 6.2 m/s, dans cette version simplifiée, on se donne une seule onde de marée, la M2. La pulsation période T=(12*60*60+25*60+12). Sans dimension la hauteur d’entrée est donc, aved A l’amplitude relative de la marée mesurée en h_0: \displaystyle h = 1 + A \cos(2 \pi t). Sur les marégraphes, est environ 0.5 (différence haut moins bas: amplitude de la marée 4 m, donc deux fois l’amplitude d’un cosinus, profondeur 4 m).

    Le domaine fait alors à peu près 0.32 sans dimension car: c_0 T= 380km soit et L=89km. On suppose que tout sort et continue en amont après Bordeaux.

    Au premier ordre en linéarisé on a simplement \displaystyle \frac{\partial u}{\partial t} = -\frac{\partial h}{\partial x},\;\;\;\;\;\;\frac{\partial h}{\partial t} = -\frac{\partial u}{\partial x} la solution de ce problème est h=1 + A \cos (2 \pi (t -x)) et u= A \cos (2 \pi (t -x)).

    Numerics

    An example in simple C (not a Basilisk example). We solve here the problem with a simple scheme with Rusanov flux. \displaystyle \dfrac{U_i^{n+1}-U_i^{n}}{\Delta t}+\dfrac{F^n_{i+1/2}-F^n_{i-1/2}}{\Delta x}=0,

    #include <stdio.h>
    #include <stdlib.h> 
    #include <math.h> 
    #include <string.h>
    
      double*x=NULL,*h=NULL,*u=NULL,*Q=NULL;
      double t,dt,tmax,dx,Z; 
      int nx;

    definition of the fluxes, here a simple Rusanov \displaystyle F_i=f(U_{i-1},U_i)= \dfrac{F(U_{i-1})+F(U_i)}{2}-c\dfrac{U_{i}-U_{i-1}}{2}, which is coded as \displaystyle f(U_G,U_D)= \dfrac{F(U_G)+F(U_D)}{2}-c\dfrac{U_D-U_G}{2} with c=\sup\limits_{U=U_G,U_D}({\sup\limits_{j\in\{1,2\}}} |\lambda_j(U)|), where \lambda_1(U) and \lambda_2(U) are the eigen values.

    double FR1(double ug,double ud,double hg,double hd)
     { double c;
       c=fmax(fabs(ug)+sqrt(hg),fabs(ud)+sqrt(hd));
       return (hg*ug+hd*ud)*0.5-c*(hd-hg)*0.5;
     }
    double FR2(double ug,double ud,double hg,double hd)
     { double c;
       c=fmax(fabs(ug)+sqrt(hg),fabs(ud)+sqrt(hd));
       return (ug*ug*hg + hg*hg/2. + ud*ud*hd + hd*hd/2.)*0.5 - c*(hd*ud-hg*ug)*0.5;
     }
      
    /*     -------------------------------------------------    */        
    int main (int argc, const char *argv[]) {
        int  i,it=0;
        FILE *g;
        double*fp=NULL,*fd=NULL,*un=NULL,*hn=NULL;
        double  q;
        
    // parametres --------------------------------------------------------------  
    dt=0.001;
    tmax=5;
    dx=0.0025;
    nx=128;
    t=0;

    domaine de longueur 2

      fprintf(stderr,"  ---------\\                /        \n"); 
      fprintf(stderr,"            \\             /  -->    \n");   
      fprintf(stderr,"               -----------            \n");   
      fprintf(stderr,"  ____________________________________\n");   
    /* ----------------------------------------------------------------------*/      
          x= (double*)calloc(nx+1,sizeof(double));
          h= (double*)calloc(nx+1,sizeof(double));
          Q= (double*)calloc(nx+1,sizeof(double));
          u= (double*)calloc(nx+1,sizeof(double));
          fp=(double*)calloc(nx+1,sizeof(double));
          fd=(double*)calloc(nx+1,sizeof(double)); 
          un=(double*)calloc(nx+1,sizeof(double));
          hn=(double*)calloc(nx+1,sizeof(double));
          
     fprintf(stderr,"tmax= %lf \n",tmax);
          
    // initialisation cond init ----------------------------      
      for(i=0;i<=nx;i++)
        {  x[i]=(i-nx/2)*dx; 
         Z=0;
         h[i]=1;
         u[i]=0;
         Q[i]=u[i]*h[i];
         }      
    // initialisation du fichier de sortie ----------------------      
          g = fopen("solxhQt.OUT", "w");
          fclose(g); 
                      
        while(t<tmax){   // boucle en temps
          t=t+dt;
          it=it+1;      
     
       for(i=1;i<=nx;i++)
        {    
            fp[i]=FR1(u[i-1],u[i],h[i-1],h[i]);
            fd[i]=FR2(u[i-1],u[i],h[i-1],h[i]);  
        }
       
       for(i=1;i<nx;i++)
        {
          hn[i]=h[i]- dt*(fp[i+1]-fp[i])/dx;   //conservation de la masse
      if(h[i]>0.){                             //conservation quantité de mouvement
          q=h[i]*u[i]-dt*(fd[i+1]-fd[i])/dx ;  
          un[i]=q/hn[i];}
          else{
          un[i]=0.;}
        }

    on ajoute un frottement de Manning donc la vitesse est à corriger de :

    \displaystyle \frac{\partial u}{\partial t} =- n^2\frac{c_0gT}{ h_0^{4/3}}(\frac{ u^2}{ h^{4/3}}) A.N.: cf=sqrt(9.81*4)*9.81*(12*60*60+25*60+12)/(4**(4./3))

    donc pour fond avec n=0.02 (graviers), cf=173 (cf table)

       double cf =173;
       for(i=1;i<nx;i++) { if(h[i]>0.) un[i]= un[i]/(1+cf*dt*fabs(un[i])/pow(hn[i],4./3)); }

    On diminue l’amplitude A car il n’y a aucune dissipation

    on écrit que les ondes entrent et ne sortent pas u +2 c =cst donc u + 2 sqrt(h) est constant.

    en vives eaux 2.5m d’amplitude, (5 m de la basse à la haute)) A=2.5/4 en mortes eaux 1.5 (3 m de bas à haut)

        double A=1.5/4,omega=2*3.141516;
        hn[0]=-hn[1] + (1+ A*cos(omega*t)) *2 ;
        un[0]=un[1]-2*(sqrt(hn[1])-sqrt(hn[0]));  
        hn[nx]=hn[nx-1];
        un[nx]=un[nx-1];  
       
       //swap      
      for(i=0;i<=nx;i++)
        { 
          h[i]=hn[i];
          u[i]=un[i];
          Q[i]=hn[i]*un[i];
        }
    
    
    
      if(it%10==0){
        /* Saving the fields */ 
      g = fopen("solxhQt.OUT", "a");
      for (i=0; i<=nx;i++)
      {      fprintf(g,"%lf %lf %lf %lf \n",x[i],h[i],Q[i],t);}
       fprintf(g,"\n");
       fclose(g);
      }
    #ifdef gnuX
    if(it%20==0){  
      printf("t=%lf\n",t);
      printf("set xlabel \"x\" ; set title \"t=%lf \"\n",t);
         y1=-.25,y2=1.5;   
      printf("h(x)=1+%lf*cos(%lf*(t-x -%lf)) \n",A,omega,nx*dx/2);
    
       
      printf("p[%lf:%lf][%lf:%lf] h(x) t'h exact','-' u 1:2 t'Q'  w l,'-'  t'h' w lp,'-'  t'Z' w l linec -1\n ",
             -nx*dx/2,nx*dx/2,y1,y2); 
        for (i=0; i<=nx;i++)
        {
        printf("%lf %lf \n",x[i],Q[i]);}
        printf("e \n");
        for (i=0; i<=nx;i++)
        {
        printf("%lf %lf \n",x[i],h[i]+Z);}
        printf("e \n");
        for (i=0; i<=nx;i++)
        {
          printf("%lf %lf \n",x[i],Z);}
        printf("e \n");}
    #endif   
          }
          
         free(x);
         free(fp);
         free(fd); 
         free(un);
         free(hn);  
        return 0;
    }

    Run

    Programme en C simple fichier (en cours on a regardé un fichier plus compliqué). Cas de rupture de barrage flux Rusanov (essayer HLC) compilation on le compile et on crée l’exécutable embouch

    cc -O3 -ffast-math -std=c99 -lm embouch.c -o embouch

    pour lancer le programme:

     ./embouch 

    Alternatively, we can compile with gnuX option which allows to pipe in gnuplot

     cc  -O3 -ffast-math -std=c99 -lm  -DgnuX  embouch.c -o embouch    

    to run the code

     ./embouch | gnuplot -persist

    Results

    un fichier appelé solxhQt.OUT avec x h Q et t est créé, pour le tracer, on lance gnuplot:

    Tracé à t=nT+T/5, nT + .1 T, nT + .2T, nT + .3T, de la hauteur d’eau en fonction de x:

     set ylabel "h"; set xlabel "x";  
      p[:][0:]'solxhQt.OUT' u 1:($4==4.?$2:NaN) t 't=0.0',''u 1:($4==4.1?$2:NaN)t't=0.1',''u 1:($4==4.2?$2:NaN)t't=0.2',''u 1:($4==4.3?$2:NaN)t't=0.3'
     # p[ : ][0:]'solxhQt.OUT' u 1:($4==4.2?$2:NaN) not,''u 1:($4==4.4?$2:NaN) not,'' u 1:($4==4.6?$2:NaN) not,'' u 1:($4==4.8?$2:NaN) not,'' u 1:($4==5?$2:NaN) not
    (script)

    (script)

    Ce qui donne pour 1 (Pointe de la Grave), 2 (Richards), 3 (Lamena), 4 (Pauillac), 4 (Bordeaux):

     set ylabel "h"; set xlabel "t";  
     p[2.5:4.5][0:2]'solxhQt.OUT' u 4:($1==-.16?$2:NaN) t'1',''u 4:($1==-.08?$2:NaN) t'2',\
      '' u 4:($1==0?$2:NaN) t'3',''u 4:($1==.08?$2:NaN) t'4',''u 4:($1==.16?$2:NaN) t'5'
    (script)

    (script)

    Go further

    add a slope…

    Do it with Basilsik, use the real … topography

    Links

    Biblio

    version sept 2018 Abidjan PYL