3 Exemples documentés – CERF : Computation, ERosion and Flow V1.0

# 3 Exemples documentés

## 3.1 Exemples avec le modèle de Saint-Venant

Dans ce chapitre, nous présentons la mise en oeuvre de quelques cas tests classiques d’écoulements pour le modèle de Saint-Venant. Les cas traités sont les suivants :

- Résolution d’un problème de Riemann
- Ressaut hydraulique
- Friction, cas test de McDonald
- …

### 3.1.1 Problème de Riemann

On considère le cas classique de la rupture de barrage 1D pour le modèle de Saint-Venant sur un fond plat. On reprend le cas test de Faccanoni and Mangeney (2013) (cas test 5). Le domaine est défini par \(x\in[-2.,2.]\) et \(t\in[0,0.5]\). Les conditions initiales sont données par : \[ \overrightarrow{w_p}(x,0) =\left\{ \begin{array}{ll} \overrightarrow{w_p}_L = <h_L,u_L>^T = <0.0046,0.>^T &\mbox{ si } x<0 \\ \overrightarrow{w_p}_R = <h_R,u_R>^T = <0.1446,0.>^T &\mbox{ si } x>0 \end{array} \right. \]
On ainsi cas de rupture de barrage sans zone sèche, avec un choc et une détente, dont la solution exacte peut être calculée par le solveur de Riemann.

Dans un premier temps, on considère un maillage 1D uniforme et fixe de 400 cellules sur un seul domaine. Le fichier de données correspondant au cas traité est fourni ci-après:

Tip 3.1: Fichier de données ex1_sv_Riemann.inp

    !-----------------------------------------------------------------------
    ! Data file in CERF format
    ! Each block of data is characterized in the first line by a header (4 characters)
    ! then a code (1 integer) and an argument (1 integer)
    ! PHYS MAME INIT BATH COND MESH OBST NUME
    !-----------------------------------------------------------------------

    !-----------------------------------------------------------------------
    ! Reading physical parameters
    ! Header PHYS
    ! 1 keyword (4 characters) and a code (1 integer)
    ! for now we have the same values for all domains
    !-----------------------------------------------------------------------
    ! MODE : physical model
    !        1/ Shallow water
    !        2/ Euler isothermal two-phases               (default 0.)
    !---- for isothermal Euler model
    ! GPES : Gravity Orientation  0/none, 1/x, 2/y, 3/z   (default 2.)
    ! PREF : isothermal pressure of reference             (default 1.d5)
    ! CSSM : Mixture sound velocity (isothermal)          (default 20.)
    ! REFA : Isothermal reference density for Air         (default 1.)
    ! REFW : Isothermal reference density for Water       (default 1000.)
    ! SHAR : isothermal interface sharpening parameter
    !        recommended value 0.01                       (defaut 0.)
    !---- for Shallow water model
    ! WLTR : Value to truncate the Water level            (default 0.0001)
    ! FRMO : Friction Model
    ! 0/none, 1/Manning-Strickler, 2/Darcy-Weisbach       (default 0.)
    ! BAGR : Bathymetry gradient 0/ Given , 1/ numerical approximation
    !-----------------------------------------------------------------------
    PHYS 0
    MODE 1.

    !-----------------------------------------------------------------------
    ! Reading the MAster MEsh
    ! Header  MAME
    !-----------------------------------------------------------------------
    ! code 1 -> creation of a mesh
    !
    !           argument: 0 type "1D shock tube"
    ! next line: xmin,xmax
    ! next line: nx
    !
    !           argument: 1 type "2d box"
    ! next line: xmin,xmax,ymin,ymax
    ! next line: nx,ny
    !
    !           argument: 2 type "3d box"
    ! next line: xmin,xmax,ymin,ymax,zmin,zmax
    ! next line: nx,ny,nz
    !
    ! code 2 -> Reading an ascii file format GMSH V2.2
    ! next line: file name
    !-----------------------------------------------------------------------
    MAME 1 0
    -2. 2.
    400

    !-----------------------------------------------------------------------
    ! Reading Initial conditions
    ! Header INIT
    !-----------------------------------------------------------------------
    ! code 0 -> Definition for a shock tube
    !           primitive variables on left
    !           primitive variables on right
    !           argument: 0
    ! code 1 -> user function
    !           argument: Function's Number
    ! code 2 -> Definition using areas
    !           argument: number of areas
    !           and for each area 2 lignes
    !           - xmin,xmax,ymin,ymax,zmin,zmax
    !           - nvar values ( primitive variables)
    ! Shallow water model: water level, x velocity, y velocity
    ! Air-water flow model: density, x velocity, y velocity,z velocity, pressure
    !-----------------------------------------------------------------------
    INIT 0 0
    0.0046 0. 0.
    0.1446 0. 0.

    !-----------------------------------------------------------------------
    ! Reading bathymetry and friction
    ! Header BATH
    !-----------------------------------------------------------------------
    ! code 0 -> Default : flat bottom and no friction

    ! code 1 -> bathymetry and friction defined by user function (defined in sources/PHY/fct_bathy.f90)
    !           argument: Function's Number
    ! code 2 -> bathymetry and friction defined linearly per block
    !           argument: 0  reading the file
    !           argument:-1  creation of a file containing the coordinates of the vertices of the block
    ! next line: file name
    !-----------------------------------------------------------------------
    BATH    0   0

    !-----------------------------------------------------------------------
    ! Reading boundary condition
    ! (maximum 10 for the moment)
    ! Header COND
    !-----------------------------------------------------------------------
    !
    ! code 0 -> Description in the case  "shock tube" : Nothing
    !
    ! code 1 -> Description in the case  "2d box" (at least 1 line per type) :
    !           kind of boundary condition in xmin, xmax, ymin, ymax
    !
    ! code 2 -> Description in the case  "3d box" (at least 1 line per type) :
    !           kind of boundary condition in xmin, xmax, ymin, ymax, zmin,zmax
    !
    ! code 3 -> Definition per zone for  GMSH mesh
    !           argument: number of zones defining a boundary condition
    !           number of zone, kind of boundary condition
    !
    ! kind of boundary condition:
    ! 0 outlet (we copy)
    ! 1 miror
    ! 2 Dirichlet condition (if 1.e20 then copy), then next line nvar values
    ! 3 user's function , then next line function number
    ! 4 Dirichlet condition on conservative variables  (if 1.e20 then copy),
    !   then next line nvar values
    ! 999 Non used
    !-----------------------------------------------------------------------
    COND    0    0

    !-----------------------------------------------------------------------
    ! Reading mesh parameters in order to define the mesh according to
    ! the Master mesh is the one represented by the domain "0".
    ! Header MESH
    !-----------------------------------------------------------------------
    !
    ! code 0 -> default value
    !
    ! 1 keyword (4 characters) and a real
    ! NBDS : NumBer of DomainS                         (default 001.)
    ! MARL : MAximal Refinement Level                  (default 000.)
    ! COPA : Mesh COarsening PArameter          0<..<1 (default 0.002)
    ! REPA : Mesh REfinement PArameter          0<..<1 (default 0.02)
    ! FDRA : Function defining the mesh refinement to be applied
    ! NZRA : Number of zones where initial blocks are refined (default 1.)
    !        then for each zone
    !        nrb,xmin,xmax,ymin,ymax,zmin,zmax
    !-----------------------------------------------------------------------
    MESH    0
    NBDS 1
    MARL 0
    NZRA 1
    0 -2.  2. -100. 100. -100. 100.

    !-----------------------------------------------------------------------
    ! Reading numerical parameter
    ! Header  NUME
    !-----------------------------------------------------------------------
    ! 1 keyword (4 characters) then a real
    ! code 0 -> default value
    !
    ! TINT : Time interval to compute                    (default 0.)
    ! CCFL : cfl condition                               (default 0.9)
    ! TDOR : Time discretization order 1 ou 2         (defaut 1.)
    ! SDOR : Space discretization order 1 ou 2        (defaut 1.)
    ! LIMI : kind of limiter  0/Barth, 1/Cartesian       (default 0.)
    ! SAVE : kind of backup, three-digit number          (default 001.)
    !        decade  : shock tube backup 0/no 1/yes
    !        unit    : binary backup     0/no 1/yes
    ! NPRO : number of probe  < 10                      (defaut 000.)
    !        then coordinates x,y,z of the probe (1 probe per line)
    !-----------------------------------------------------------------------
    NUME 0
    TINT .5
    CCFL .9
    SDOR 2.
    TDOR 2.
    SAVE 11.

Pour lancer le calcul, comme il n’y a qu’un domaine à traiter sur un seul intervalle de temps :

    ./cerf_input ex1_sv_Riemann.inp
    ./cerf

Comme dans le fichier de données ex1_sv_Riemann.inp on a demandé une sauvegarde au format 1D (paramètre SAVE=11 section NUME), des fichiers de sorties 1D ont été générés. Des fichiers pour la solution numérique en hauteur d’eau et vitesse ‘n_h_001’ et ‘n_u_001’, et des fichiers pour la solution exacte en hauteur d’eau et vitesse : ‘e_h_001’ et ‘e_u_001’. Ces fichiers peuvent être visualisés par exemple avec gnuplot :

    gnuplot
    gnuplot> set term x11
    gnuplot> plot 'e_h_001' w l, 'n_h_001' w l
    gnuplot> plot 'e_u_001' w l, 'n_u_001' w l
    gnuplot> quit

On peut ainsi comparer la solution numérique à la solution exacte.

Comparaison numérique/exacte hauteur d’eau à \(t=0.5s\) , maillage uniforme 400 cellules

Comparaison numérique/exacte vitesse à \(t=0.5s\), maillage uniforme 400 cellules

Pour tester le code avec une zone “sèche” (Faccanoni and Mangeney (2013) cas test 1), on peut par exemple reprendre le même problème avec une condition initiale différente en modifiant la section INIT du fichier de données ex1_sv_Riemann.inp:

    INIT    0   0
    0. 0. 0.
    0.1446 0. 0.

On obtiendra alors :

Comparaison numérique/exacte hauteur d’eau à \(t=0.5s\) , maillage uniforme 400 cellules avec zone sèche

Comparaison numérique/exacte vitesse à \(t=0.5s\), maillage uniforme 400 cellules avec zone sèche

### 3.1.2 Ressaut hydraulique

On considère ici un cas tests stationnaire, sans friction, avec un fond non plat. La solution analytique est issue de Delestre et al. (2013). Sur un domaine \(x\in[0.,25.]\), on considère une topographie définie par :
\[z_b(x) = \left\{ \begin{array}{ll} 0.2 - 0.05 \left( x - 10 \right)^2 & \mbox{ si } 8 < x < 12 \\ 0. & \mbox{ sinon } \end{array}\right. \]

On considère le cas de l’écoulement transcritique avec choc, caractérisé par les conditions aux limites suivantes :

- À gauche, on impose \(h = 0.33 m\) et on laisse le débit libre.
- À droite, on impose \(hu = 0.18 m^2/s\) et on laisse la hauteur d’eau.

Comme préconisé par Delestre et al. (2013), les conditions initiales sont données par : \(h + z_b =0\) et \(u=0\).

On effectue le calcul sur un maillage fixe et régulier de 2500 cellules réparties sur 8 domaines. À la recherche d’un état stationnaire, on effectue le calcul jusqu’à \(T=100s\) . Le fichier de données correspondant au cas traité est fourni ci-après Tip 3.2 :

Tip 3.2: Fichier de données ex2_sv_Bump.inp

    !-----------------------------------------------------------------------
    ! Data file in CERF format
    ! Each block of data is characterized in the first line by a header (4 characters)
    ! then a code (1 integer) and an argument (1 integer)
    ! PHYS MAME INIT BATH COND MESH OBST NUME
    !-----------------------------------------------------------------------

    !-----------------------------------------------------------------------
    ! Reading physical parameters
    ! Header PHYS
    ! 1 keyword (4 characters) and a code (1 integer)
    ! for now we have the same values for all domains
    !-----------------------------------------------------------------------
    ! MODE : physical model
    !        1/ Shallow water
    !        2/ Euler isothermal two-phases               (default 0.)
    !---- for isothermal Euler model
    ! GPES : Gravity Orientation  0/none, 1/x, 2/y, 3/z   (default 2.)
    ! PREF : isothermal pressure of reference             (default 1.d5)
    ! CSSM : Mixture sound velocity (isothermal)          (default 20.)
    ! REFA : Isothermal reference density for Air         (default 1.)
    ! REFW : Isothermal reference density for Water       (default 1000.)
    ! SHAR : isothermal interface sharpening parameter
    !        recommended value 0.01                       (defaut 0.)
    !---- for Shallow water model
    ! WLTR : Value to truncate the Water level            (default 0.0001)
    ! FRMO : Friction Model
    ! 0/none, 1/Manning-Strickler, 2/Darcy-Weisbach       (default 0.)
    ! BAGR : Bathymetry gradient 0/ Given , 1/ numerical approximation
    !-----------------------------------------------------------------------
    PHYS 0
    MODE 1.
    FRMO 0.
    BAGR 0.
    !-----------------------------------------------------------------------
    ! Reading the MAster MEsh
    ! Header  MAME
    !-----------------------------------------------------------------------
    ! code 1 -> creation of a mesh
    !
    !           argument: 0 type "1D shock tube"
    ! next line: xmin,xmax
    ! next line: nx
    !
    !           argument: 1 type "2d box"
    ! next line: xmin,xmax,ymin,ymax
    ! next line: nx,ny
    !
    !           argument: 2 type "3d box"
    ! next line: xmin,xmax,ymin,ymax,zmin,zmax
    ! next line: nx,ny,nz
    !
    ! code 2 -> Reading an ascii file format GMSH V2.2
    ! next line: file name
    !-----------------------------------------------------------------------
    MAME    1    0
    0. 25.
    2500
    !-----------------------------------------------------------------------
    ! Reading Initial conditions
    ! Header INIT
    !-----------------------------------------------------------------------
    ! code 0 -> Definition for a shock tube
    !           primitive variables on left
    !           primitive variables on right
    !           argument: 0
    ! code 1 -> user function
    !           argument: Function's Number
    ! code 2 -> Definition using areas
    !           argument: number of areas
    !           and for each area 2 lignes
    !           - xmin,xmax,ymin,ymax,zmin,zmax
    !           - nvar values ( primitive variables)
    ! Shallow water model: water level, x velocity, y velocity
    ! Air-water flow model: density, x velocity, y velocity,z velocity, pressure
    !-----------------------------------------------------------------------
    INIT    1   1

    !-----------------------------------------------------------------------
    ! Reading bathymetry and friction
    ! Header BATH
    !-----------------------------------------------------------------------
    ! code 0 -> Default : flat bottom and no friction

    ! code 1 -> bathymetry and friction defined by user function (defined in sources/PHY/fct_bathy.f90)
    !           argument: Function's Number
    ! code 2 -> bathymetry and friction defined linearly per block
    !           argument: 0  reading the file
    !           argument:-1  creation of a file containing the coordinates of the vertices of the block
    ! next line: file name
    !-----------------------------------------------------------------------
    BATH    1   1

    !-----------------------------------------------------------------------
    ! Reading boundary condition
    ! (maximum 10 for the moment)
    ! Header COND
    !-----------------------------------------------------------------------
    !
    ! code 0 -> Description in the case  "shock tube" : Nothing
    !
    ! code 1 -> Description in the case  "2d box" (at least 1 line per type) :
    !           kind of boundary condition in xmin, xmax, ymin, ymax
    !
    ! code 2 -> Description in the case  "3d box" (at least 1 line per type) :
    !           kind of boundary condition in xmin, xmax, ymin, ymax, zmin,zmax
    !
    ! code 3 -> Definition per zone for  GMSH mesh
    !           argument: number of zones defining a boundary condition
    !           number of zone, kind of boundary condition
    !
    ! kind of boundary condition:
    ! 0 outlet (we copy)
    ! 1 miror
    ! 2 Dirichlet condition (if 1.e20 then copy), then next line nvar values
    ! 3 user's function , then next line function number
    ! 4 Dirichlet condition on conservative variables  (if 1.e20 then copy),
    !   then next line nvar values
    ! 999 Non used
    !-----------------------------------------------------------------------
    COND    2    0
    4
    1.e20 0.18 0.
    2
    0.33 1.e20 0.
    999
    999
    999
    999

    !-----------------------------------------------------------------------
    ! Reading mesh parameters in order to define the mesh according to
    ! the Master mesh is the one represented by the domain "0".
    ! Header MESH
    !-----------------------------------------------------------------------
    !
    ! code 0 -> default value
    !
    ! 1 keyword (4 characters) and a real
    ! NBDS : NumBer of DomainS                         (default 001.)
    ! MARL : MAximal Refinement Level                  (default 000.)
    ! COPA : Mesh COarsening PArameter          0<..<1 (default 0.002)
    ! REPA : Mesh REfinement PArameter          0<..<1 (default 0.02)
    ! FDRA : Function defining the mesh refinement to be applied
    ! NZRA : Number of zones where initial blocks are refined (default 1.)
    !        then for each zone
    !        nrb,xmin,xmax,ymin,ymax,zmin,zmax
    !-----------------------------------------------------------------------
    MESH 0
    NBDS 8
    NZRA 1
    0 0.  25. -100. 100. -100. 100.

    !-----------------------------------------------------------------------
    ! Reading numerical parameter
    ! Header  NUME
    !-----------------------------------------------------------------------
    ! 1 keyword (4 characters) then a real
    ! code 0 -> default value
    !
    ! TINT : Time interval to compute                    (default 0.)
    ! CCFL : cfl condition                               (default 0.9)
    ! TDOR : Time discretization order 1 ou 2         (defaut 1.)
    ! SDOR : Space discretization order 1 ou 2        (defaut 1.)
    ! LIMI : kind of limiter  0/Barth, 1/Cartesian       (default 0.)
    ! SAVE : kind of backup, three-digit number          (default 001.)
    !        decade  : shock tube backup 0/no 1/yes
    !        unit    : binary backup     0/no 1/yes
    ! NPRO : number of probe  < 10                      (defaut 000.)
    !        then coordinates x,y,z of the probe (1 probe per line)
    !-----------------------------------------------------------------------
    NUME 0
    TINT 100.
    CCFL .9
    SDOR 2.
    TDOR 2.
    SAVE 01.

Pour lancer le calcul, comme il n’y a qu’un domaine à traiter sur un seul intervalle de temps, mais sur 8 domaines, on utilise les commandes suivantes :

    ./cerf_input ex2_sv_Bump.inp
    mpirun -np 8 ./cerf
    ./cerf2tec

On peut récupérer la solution analytique de ce problème sur le site de Swashes. Dans le répertoire doc de CERF on trouvera à titre d’information la solution analytique de \(\eta = h + z_b\) sur 2500 points au format Tecplot. On peut alors comparer la solution numérique à la solution exacte sur la figure ci-après (Figure 3.1).

Figure 3.1: Comparaison numérique (Rouge) / exacte (Noire) niveau d’eau à \(t=100.s\) , maillage uniforme 2500 cellules

On peut également visualiser la solution numérique avec Visit ou Paraview:

Figure 3.2: Isovaleur niveau d’eau à \(t=100.s\) , maillage uniforme 2500 cellules

Sur la figure Figure 3.2, on visualise les isovaleurs du niveau d’eau à \(t=100s\). On remarquera que le maillage a été modifié pour le post-traitement. Les coordonnées des nœuds de la base des hexaèdres correspondent à \(z_b\), alors que les nœuds de la partie supérieure des hexaèdres correspondent à \(z_b+h\). On peut ainsi visualiser la topographie et le niveau d’eau en même temps.

On peut également réaliser ce calcul en utilisant le raffinement de maillage automatique. Pour cela, on modifie le fichier de données en conséquence ex2_sv_Bump.inp. On crée alors un maillage constitué de seulement 250 blocs et on autorise un raffinement de maillage seulement jusqu’au niveau 2. On garde tous les autres paramètres du calcul comme pour le maillage fixe.

Tip 3.3: Fichier de données ex2_sv_Bump_amr.inp

    !-----------------------------------------------------------------------
    ! Data file in CERF format
    ! Each block of data is characterized in the first line by a header (4 characters)
    ! then a code (1 integer) and an argument (1 integer)
    ! PHYS MAME INIT BATH COND MESH OBST NUME
    !-----------------------------------------------------------------------

    !-----------------------------------------------------------------------
    ! Reading physical parameters
    ! Header PHYS
    ! 1 keyword (4 characters) and a code (1 integer)
    ! for now we have the same values for all domains
    !-----------------------------------------------------------------------
    ! MODE : physical model
    !        1/ Shallow water
    !        2/ Euler isothermal two-phases               (default 0.)
    !---- for isothermal Euler model
    ! GPES : Gravity Orientation  0/none, 1/x, 2/y, 3/z   (default 2.)
    ! PREF : isothermal pressure of reference             (default 1.d5)
    ! CSSM : Mixture sound velocity (isothermal)          (default 20.)
    ! REFA : Isothermal reference density for Air         (default 1.)
    ! REFW : Isothermal reference density for Water       (default 1000.)
    ! SHAR : isothermal interface sharpening parameter
    !        recommended value 0.01                       (defaut 0.)
    !---- for Shallow water model
    ! WLTR : Value to truncate the Water level            (default 0.0001)
    ! FRMO : Friction Model
    ! 0/none, 1/Manning-Strickler, 2/Darcy-Weisbach       (default 0.)
    ! BAGR : Bathymetry gradient 0/ Given , 1/ numerical approximation
    !-----------------------------------------------------------------------
    PHYS 0
    MODE 1.

    !-----------------------------------------------------------------------
    ! Reading the MAster MEsh
    ! Header  MAME
    !-----------------------------------------------------------------------
    ! code 1 -> creation of a mesh
    !
    !           argument: 0 type "1D shock tube"
    ! next line: xmin,xmax
    ! next line: nx
    !
    !           argument: 1 type "2d box"
    ! next line: xmin,xmax,ymin,ymax
    ! next line: nx,ny
    !
    !           argument: 2 type "3d box"
    ! next line: xmin,xmax,ymin,ymax,zmin,zmax
    ! next line: nx,ny,nz
    !
    ! code 2 -> Reading an ascii file format GMSH V2.2
    ! next line: file name
    !-----------------------------------------------------------------------
    MAME    1    0
    0. 25.
    2500
    !-----------------------------------------------------------------------
    ! Reading Initial conditions
    ! Header INIT
    !-----------------------------------------------------------------------
    ! code 0 -> Definition for a shock tube
    !           primitive variables on left
    !           primitive variables on right
    !           argument: 0
    ! code 1 -> user function
    !           argument: Function's Number
    ! code 2 -> Definition using areas
    !           argument: number of areas
    !           and for each area 2 lignes
    !           - xmin,xmax,ymin,ymax,zmin,zmax
    !           - nvar values ( primitive variables)
    ! Shallow water model: water level, x velocity, y velocity
    ! Air-water flow model: density, x velocity, y velocity,z velocity, pressure
    !-----------------------------------------------------------------------
    INIT    1   1

    !-----------------------------------------------------------------------
    ! Reading bathymetry and friction
    ! Header BATH
    !-----------------------------------------------------------------------
    ! code 0 -> Default : flat bottom and no friction

    ! code 1 -> bathymetry and friction defined by user function (defined in sources/PHY/fct_bathy.f90)
    !           argument: Function's Number
    ! code 2 -> bathymetry and friction defined linearly per block
    !           argument: 0
    ! next line: file name
    !
    ! code 3 -> bathymetry and friction defined  per block at the highest mesh refinement level
    !           argument: 0
    ! next line: file name
    !-----------------------------------------------------------------------
    BATH    1   1

    !-----------------------------------------------------------------------
    ! Reading boundary condition
    ! (maximum 10 for the moment)
    ! Header COND
    !-----------------------------------------------------------------------
    !
    ! code 0 -> Description in the case  "shock tube" : Nothing
    !
    ! code 1 -> Description in the case  "2d box" (at least 1 line per type) :
    !           kind of boundary condition in xmin, xmax, ymin, ymax
    !
    ! code 2 -> Description in the case  "3d box" (at least 1 line per type) :
    !           kind of boundary condition in xmin, xmax, ymin, ymax, zmin,zmax
    !
    ! code 3 -> Definition per zone for  GMSH mesh
    !           argument: number of zones defining a boundary condition
    !           number of zone, kind of boundary condition
    !
    ! kind of boundary condition:
    ! 0 outlet (we copy)
    ! 1 miror
    ! 2 Dirichlet condition (if 1.e20 then copy), then next line nvar values
    ! 3 user's function , then next line function number
    ! 4 Dirichlet condition on conservative variables  (if 1.e20 then copy),
    !   then next line nvar values
    ! 999 Non used
    !-----------------------------------------------------------------------
    COND    2    0
    4
    1.e20 0.18 0. 0. 0. 0.
    2
    0.33 1.e20 0. 0. 0. 0.
    999
    999
    999
    999

    !-----------------------------------------------------------------------
    ! Reading mesh parameters in order to define the mesh according to
    ! the Master mesh is the one represented by the domain "0".
    ! Header MESH
    !-----------------------------------------------------------------------
    !
    ! code 0 -> default value
    !
    ! 1 keyword (4 characters) and a real
    ! NBDS : NumBer of DomainS                         (default 001.)
    ! MARL : MAximal Refinement Level                  (default 000.)
    ! COPA : Mesh COarsening PArameter          0<..<1 (default 0.002)
    ! REPA : Mesh REfinement PArameter          0<..<1 (default 0.02)
    ! FDRA : Function defining the mesh refinement to be applied
    ! NZRA : Number of zones where initial blocks are refined (default 1.)
    !        then for each zone
    !        nrb,xmin,xmax,ymin,ymax,zmin,zmax
    !-----------------------------------------------------------------------
    MESH 0
    NBDS 8
    MARL 2
    NZRA 2
    0 0.  25. -100. 100. -100. 100.
    2 7.8  12.1 -100. 100. -100. 100.

    !-----------------------------------------------------------------------
    ! Reading numerical parameter
    ! Header  NUME
    !-----------------------------------------------------------------------
    ! 1 keyword (4 characters) then a real
    ! code 0 -> default value
    !
    ! TINT : Time interval to compute                    (default 0.)
    ! CCFL : cfl condition                               (default 0.9)
    ! TDOR : Time discretization order 1 ou 2         (defaut 1.)
    ! SDOR : Space discretization order 1 ou 2        (defaut 1.)
    ! LIMI : kind of limiter  0/Barth, 1/Cartesian       (default 0.)
    ! SAVE : kind of backup, three-digit number          (default 001.)
    !        decade  : shock tube backup 0/no 1/yes
    !        unit    : binary backup     0/no 1/yes
    ! NPRO : number of probe  < 10                      (defaut 000.)
    !        then coordinates x,y,z of the probe (1 probe per line)
    !-----------------------------------------------------------------------
    NUME 0
    TINT 10.
    CCFL .9
    SDOR 2.
    TDOR 2.
    SAVE 01.

Afin de lancer le calcul avec raffinement de maillage, on utilise un script :

    #!/bin/bash
    ./cerf_input ex2_sv_Bump2.inp
    ./cerf2tec
    mv cerfout.dat b_0.dat
    mv cerfamr.dat b_mame_0.dat
    for i in $(seq 1 20 )
    do
    mpirun -np 8 ./cerf
    ./cerf2tec
    mv cerfout.dat b_$i.dat
    ./cerf_amr
    mv cerfamr.dat b_mame_$i.dat
    done

Comme précédemment, on peut comparer avec la solution exacte.

Figure 3.3: Comparaison numérique sur maillage fixe (Rouge) / numérique sur maillage dynamique (Vert) / exacte (Noire) niveau d’eau à \(t=100.s\)

On peut également visualiser la solution numérique avec Visit ou Paraview :

Figure 3.4: Isovaleur niveau d’eau, maillage dynamique 344 cellules à \(t=100.s\)

Que soit sur Figure 3.3 ou Figure 3.4, on remarquera que le raffinement de maillage a permis de réduire le nombre de cellules tout en conservant une bonne précision sur la solution numérique.

### 3.1.3 Friction, cas test de McDonald

On considère ici un cas tests stationnaire, avec friction sur un fond non plat. La solution semi-analytique est issue de Delestre et al. (2013) (cas 3.2.1).

Sur un domaine \(x\in[0.,1000.]\), on considère un écoulement sub-critique avec un débit constant à l’entrée de \(q = h u = 2 m^2/s\) et une hauteur d’eau prescrite à la sortie \(h_{e}(1000)\). Le domaine est initialement sec, avec \(h = 0\) et \(u = 0\). On impose une friction de Manning avec \(n = 0.033\).
La hauteur d’eau est donnée par : \[ h_{ex}(x) = \left( \frac{4}{g} \right)^{1/3} \left( 1 + \frac{1}{2} exp \left( - 16 \left( \frac{x}{1000} - \frac{1}{2} \right) \right) \right)\]

En utilisant Equation 1.6, dans le cas stationnaire, on peut alors retrouver la topographie par intégration de : \[\frac{\partial z_b}{\partial x} = \left( \frac{q^2}{g h^3} - 1 \right) \frac{\partial h}{\partial x} - S_f\]

On effectue le calcul sur un maillage fixe et régulier de 1000 cellules. À la recherche d’un état stationnaire, on effectue le calcul jusqu’à \(T=1500s\) . Le fichier de données correspondant au cas traité est fourni ci-après Tip 3.4 :

Tip 3.4: Fichier de données ex3_sv_McDo.inp

    !-----------------------------------------------------------------------
    ! Data file in CERF format
    ! Each block of data is characterized in the first line by a header (4 characters)
    ! then a code (1 integer) and an argument (1 integer)
    ! PHYS MAME INIT BATH COND MESH OBST NUME
    !-----------------------------------------------------------------------

    !-----------------------------------------------------------------------
    ! Reading physical parameters
    ! Header PHYS
    ! 1 keyword (4 characters) and a code (1 integer)
    ! for now we have the same values for all domains
    !-----------------------------------------------------------------------
    ! MODE : physical model
    !        1/ Shallow water
    !        2/ Euler isothermal two-phases               (default 0.)
    !---- for isothermal Euler model
    ! GPES : Gravity Orientation  0/none, 1/x, 2/y, 3/z   (default 2.)
    ! PREF : isothermal pressure of reference             (default 1.d5)
    ! CSSM : Mixture sound velocity (isothermal)          (default 20.)
    ! REFA : Isothermal reference density for Air         (default 1.)
    ! REFW : Isothermal reference density for Water       (default 1000.)
    ! SHAR : isothermal interface sharpening parameter
    !        recommended value 0.01                       (defaut 0.)
    !---- for Shallow water model
    ! WLTR : Value to truncate the Water level            (default 0.0001)
    ! FRMO : Friction Model
    ! 0/none, 1/Manning-Strickler, 2/Darcy-Weisbach       (default 0.)
    ! BAGR : Bathymetry gradient 0/ Given , 1/ numerical approximation
    !-----------------------------------------------------------------------
    PHYS 0
    MODE 1.
    FRMO 1.

    !-----------------------------------------------------------------------
    ! Reading the MAster MEsh
    ! Header  MAME
    !-----------------------------------------------------------------------
    ! code 1 -> creation of a mesh
    !
    !           argument: 0 type "1D shock tube"
    ! next line: xmin,xmax
    ! next line: nx
    !
    !           argument: 1 type "2d box"
    ! next line: xmin,xmax,ymin,ymax
    ! next line: nx,ny
    !
    !           argument: 2 type "3d box"
    ! next line: xmin,xmax,ymin,ymax,zmin,zmax
    ! next line: nx,ny,nz
    !
    ! code 2 -> Reading an ascii file format GMSH V2.2
    ! next line: file name
    !-----------------------------------------------------------------------
    MAME    1    0
    0. 1000.
    1000
    !-----------------------------------------------------------------------
    ! Reading Initial conditions
    ! Header INIT
    !-----------------------------------------------------------------------
    ! code 0 -> Definition for a shock tube
    !           primitive variables on left
    !           primitive variables on right
    !           argument: 0
    ! code 1 -> user function
    !           argument: Function's Number
    ! code 2 -> Definition using areas
    !           argument: number of areas
    !           and for each area 2 lignes
    !           - xmin,xmax,ymin,ymax,zmin,zmax
    !           - nvar values ( primitive variables)
    ! Shallow water model: water level, x velocity, y velocity
    ! Air-water flow model: density, x velocity, y velocity,z velocity, pressure
    !-----------------------------------------------------------------------
    INIT    2   1
    -1. 1001. -100. 100. -100. 100.
    0. 0. 0.

    !-----------------------------------------------------------------------
    ! Reading bathymetry and friction
    ! Header BATH
    !-----------------------------------------------------------------------
    ! code 0 -> Default : flat bottom and no friction

    ! code 1 -> bathymetry and friction defined by user function (defined in sources/PHY/fct_bathy.f90)
    !           argument: Function's Number
    ! code 2 -> bathymetry and friction defined linearly per block
    !           argument: 0  reading the file
    !           argument:-1  creation of a file containing the coordinates of the vertices of the block
    ! next line: file name
    !-----------------------------------------------------------------------

    BATH    1   2
    !-----------------------------------------------------------------------
    ! Reading boundary condition
    ! (maximum 10 for the moment)
    ! Header COND
    !-----------------------------------------------------------------------
    !
    ! code 0 -> Description in the case  "shock tube" : Nothing
    !
    ! code 1 -> Description in the case  "2d box" (at least 1 line per type) :
    !           kind of boundary condition in xmin, xmax, ymin, ymax
    !
    ! code 2 -> Description in the case  "3d box" (at least 1 line per type) :
    !           kind of boundary condition in xmin, xmax, ymin, ymax, zmin,zmax
    !
    ! code 3 -> Definition per zone for  GMSH mesh
    !           argument: number of zones defining a boundary condition
    !           number of zone, kind of boundary condition
    !
    ! kind of boundary condition:
    ! 0 outlet (we copy)
    ! 1 miror
    ! 2 Dirichlet condition (if 1.e20 then copy), then next line nvar values
    ! 3 user's function , then next line function number
    ! 4 Dirichlet condition on conservative variables  (if 1.e20 then copy),
    !   then next line nvar values
    ! 999 Non used
    !-----------------------------------------------------------------------
    COND    2    0
    4
    1.e20 2. 0.
    2
    0.74832355831838937 1.e20 0. 0.
    999
    999
    999
    999

    !-----------------------------------------------------------------------
    ! Reading mesh parameters in order to define the mesh according to
    ! the Master mesh is the one represented by the domain "0".
    ! Header MESH
    !-----------------------------------------------------------------------
    !
    ! code 0 -> default value
    !
    ! 1 keyword (4 characters) and a real
    ! NBDS : NumBer of DomainS                         (default 001.)
    ! MARL : MAximal Refinement Level                  (default 000.)
    ! COPA : Mesh COarsening PArameter          0<..<1 (default 0.002)
    ! REPA : Mesh REfinement PArameter          0<..<1 (default 0.02)
    ! FDRA : Function defining the mesh refinement to be applied
    ! NZRA : Number of zones where initial blocks are refined (default 1.)
    !        then for each zone
    !        nrb,xmin,xmax,ymin,ymax,zmin,zmax
    !-----------------------------------------------------------------------
    MESH 0
    NBDS 1
    MARL 0
    NZRA 1
    0 0.  1000. -100. 100. -100. 100.

    !-----------------------------------------------------------------------
    ! Reading numerical parameter
    ! Header  NUME
    !-----------------------------------------------------------------------
    ! 1 keyword (4 characters) then a real
    ! code 0 -> default value
    !
    ! TINT : Time interval to compute                    (default 0.)
    ! CCFL : cfl condition                               (default 0.9)
    ! TDOR : Time discretization order 1 ou 2         (defaut 1.)
    ! SDOR : Space discretization order 1 ou 2        (defaut 1.)
    ! LIMI : kind of limiter  0/Barth, 1/Cartesian       (default 0.)
    ! SAVE : kind of backup, three-digit number          (default 001.)
    !        decade  : shock tube backup 0/no 1/yes
    !        unit    : binary backup     0/no 1/yes
    ! NPRO : number of probe  < 10                      (defaut 000.)
    !        then coordinates x,y,z of the probe (1 probe per line)
    !-----------------------------------------------------------------------
    NUME 0
    TINT 50.
    CCFL .9
    SDOR 1.
    TDOR 1.
    SAVE 01.

Pour lancer le calcul, afin de visualiser l’évolution de la recherche de l’état stationnaire, on décompose le temps de calcul en \(30\) étapes en utilisant le script ci-après :

TipScript ex3_sv_McDo.sh

    #!/bin/bash
    ./cerf_input ex3_sv_McDo.inp
    ./cerf2tec
    mv cerfout1d.dat md_0.dat
    for i in $(seq 1 30 )
    do
    ./cerf
    ./cerf2tec
    mv cerfout1d.dat md_$i.dat
    done

On peut récupérer la solution analytique de ce problème sur le site de Swashes. Dans le répertoire doc de CERF on trouvera à titre d’information la solution analytique sur 10000 points au format Tecplot. On peut alors comparer la solution numérique à la solution exacte sur la figure ci-après (Figure 3.5).

Figure 3.5: Comparaison numérique (Noire) / exacte (Rouge) sur le niveau d’eau \(h\) (en haut) et la vitesse \(u\) (en bas), avec un maillage uniforme 1000 cellules à l’ordre 1 en espace et en temps

## 3.2 Exemples avec le modèle bifluide eulérien

Dans ce chapitre, nous présentons la mise en oeuvre de quelques cas tests classiques d’écoulements bifluides. Les cas traités sont les suivants :

- Résolution d’un problème de Riemann
- Rupture de barrage
- Chute d’un cylindre sur un plan d’eau
- …

### 3.2.1 Problème de Riemann

Afin de valider la mise en oeuvre du modèle, on considère le cas classique d’un “tube à choc” avec le modèle isotherme bi-fluide. On reprend le cas test de Frédéric Golay and Helluy (2007). Le domaine est défini par \(x\in[-2.,2.]\) et \(t\in[0,0.25]\). Les conditions initiales sont données par : \[ \overrightarrow{w_p}(x,0) =\left\{ \begin{array}{ll} \overrightarrow{w_p}_L = <\rho_L,u_L,p_L>^T = <0.0046,\;0,\;1.5 \;10^5>^T &\mbox{ si } x<0 \\ \overrightarrow{w_p}_R = <\rho_R,u_R,p_R>^T = <0.1446,\;0.,\;1. \;10^5>^T &\mbox{ si } x>0 \end{array} \right. \]

On considère un maillage 1D uniforme et fixe de \(400\) cellules sur un seul domaine. Le fichier de données correspondant au cas traité est fourni ci-après:

Tip 3.5: Fichier de données ex1_bf_Riemann.inp

    !-----------------------------------------------------------------------
    ! Data file in CERF format
    ! Each block of data is characterized in the first line by a header (4 characters)
    ! then a code (1 integer) and an argument (1 integer)
    ! PHYS MAME INIT BATH COND MESH OBST NUME
    !-----------------------------------------------------------------------

    !-----------------------------------------------------------------------
    ! Reading physical parameters
    ! Header PHYS
    ! 1 keyword (4 characters) and a code (1 integer)
    ! for now we have the same values for all domains
    !-----------------------------------------------------------------------
    ! MODE : physical model
    !        1/  Shallow water
    !        2/  Euler isothermal two-phases              (default 0.)
    !---- for isothermal Euler model
    ! GPES : Gravity Orientation  0/none, 1/x, 2/y, 3/z   (default 2.)
    ! PREF : isothermal pressure of reference             (default 1.d5)
    ! CSSM : Mixture sound velocity (isothermal)          (default 20.)
    ! REFA : Isothermal reference density for Air         (default 1.)
    ! REFW : Isothermal reference density for Water       (default 1000.)
    ! SHAR : isothermal interface sharpening parameter
    !        recommended value 0.01                       (defaut 0.)
    !---- for Shallow water model
    ! WLTR : Value to truncate the Water level            (default 0.0001)
    ! FRMO : Friction Model
    ! 0/none, 1/Manning-Strickler, 2/Darcy-Weisbach       (default 0.)
    ! BAGR : Bathymetry gradient 0/ Given , 1/ numerical approximation
    !-----------------------------------------------------------------------
    PHYS 0
    MODE 2.

    !-----------------------------------------------------------------------
    ! Reading the MAster MEsh
    ! Header  MAME
    !-----------------------------------------------------------------------
    ! code 1 -> creation of a mesh
    !
    !           argument: 0 type "1D shock tube"
    ! next line: xmin,xmax
    ! next line: nx
    !
    !           argument: 1 type "2d box"
    ! next line: xmin,xmax,ymin,ymax
    ! next line: nx,ny
    !
    !           argument: 2 type "3d box"
    ! next line: xmin,xmax,ymin,ymax,zmin,zmax
    ! next line: nx,ny,nz
    !
    ! code 2 -> Reading an ascii file format GMSH V2.2
    ! next line: file name
    !-----------------------------------------------------------------------
    MAME    1    0
    -10. 10.
    400
    !-----------------------------------------------------------------------
    ! Reading Initial conditions
    ! Header INIT
    !-----------------------------------------------------------------------
    ! code 0 -> Definition for a shock tube
    !           primitive variables on left
    !           primitive variables on right
    !           argument: 0
    ! code 1 -> user function
    !           argument: Function's Number
    ! code 2 -> Definition using areas
    !           argument: number of areas
    !           and for each area 2 lignes
    !           - xmin,xmax,ymin,ymax,zmin,zmax
    !           - nvar values ( primitive variables)
    ! Shallow water model: water level, x velocity, y velocity
    ! Air-water flow model: density, x velocity, y velocity,z velocity, pressure
    !-----------------------------------------------------------------------
    INIT    0   0
    1125. 0. 1.5e5
    1. 0. 1.e5

    !-----------------------------------------------------------------------
    ! Reading boundary condition
    ! (maximum 10 for the moment)
    ! Header COND
    !-----------------------------------------------------------------------
    !
    ! code 0 -> Description in the case  "shock tube" : Nothing
    !
    ! code 1 -> Description in the case  "2d box" (at least 1 line per type) :
    !           kind of boundary condition in xmin, xmax, ymin, ymax
    !
    ! code 2 -> Description in the case  "3d box" (at least 1 line per type) :
    !           kind of boundary condition in xmin, xmax, ymin, ymax, zmin,zmax
    !
    ! code 3 -> Definition per zone for  GMSH mesh
    !           argument: number of zones defining a boundary condition
    !           number of zone, kind of boundary condition
    !
    ! kind of boundary condition:
    ! 0 outlet (we copy)
    ! 1 miror
    ! 2 Dirichlet condition (if 1.e20 then copy), then next line nvar values
    ! 3 user's function , then next line function number
    ! 4 Dirichlet condition on conservative variables  (if 1.e20 then copy),
    !   then next line nvar values
    ! 999 Non used
    !-----------------------------------------------------------------------
    COND    0    0

    !-----------------------------------------------------------------------
    ! Reading mesh parameters in order to define the mesh according to
    ! the Master mesh is the one represented by the domain "0".
    ! Header MESH
    !-----------------------------------------------------------------------
    !
    ! code 0 -> default value
    !
    ! 1 keyword (4 characters) and a real
    ! NBDS : NumBer of DomainS                         (default 001.)
    ! MARL : MAximal Refinement Level                  (default 000.)
    ! COPA : Mesh COarsening PArameter          0<..<1 (default 0.002)
    ! REPA : Mesh REfinement PArameter          0<..<1 (default 0.02)
    ! FDRA : Function defining the mesh refinement to be applied
    ! NZRA : Number of zones where initial blocks are refined (default 1.)
    !        then for each zone
    !        nrb,xmin,xmax,ymin,ymax,zmin,zmax
    !-----------------------------------------------------------------------
    MESH 0
    NBDS 1
    MARL 0
    NZRA 1
    0 -10.  10. -100. 100. -100. 100.

    !-----------------------------------------------------------------------
    ! Reading numerical parameter
    ! Header  NUME
    !-----------------------------------------------------------------------
    ! 1 keyword (4 characters) then a real
    ! code 0 -> default value
    !
    ! TINT : Time interval to compute                    (default 0.)
    ! CCFL : cfl condition                               (default 0.9)
    ! TDOR : Time discretization order 1 ou 2         (defaut 1.)
    ! SDOR : Space discretization order 1 ou 2        (defaut 1.)
    ! LIMI : kind of limiter  0/Barth, 1/Cartesian       (default 0.)
    ! SAVE : kind of backup, three-digit number          (default 001.)
    !        decade  : shock tube backup 0/no 1/yes
    !        unit    : binary backup     0/no 1/yes
    ! NPRO : number of probe  < 10                      (defaut 000.)
    !        then coordinates x,y,z of the probe (1 probe per line)
    !-----------------------------------------------------------------------
    NUME 0
    TINT 0.25
    CCFL .9
    SDOR 2.
    TDOR 2.
    SAVE 11.

Pour lancer le calcul, comme il n’y a qu’un domaine à traiter sur un seul intervalle de temps :

    ./cerf_input ex1_bf_Riemann.inp
    ./cerf

Comme dans le fichier de données ex1_bf_Riemann.inp on a demandé une sauvegarde au format 1D (paramètre SAVE=11 section NUME), des fichiers de sorties 1D ont été générés. Des fichiers pour la solution numérique en densité, vitesse et pression ‘n_r_001’, ‘n_u_001’ et ‘n_p_001’, et des fichiers pour la solution exacte en densité, vitesse et pression : ‘e_r_001’, ‘e_u_001’ et ‘e_p_001’. Ces fichiers peuvent être visualisés par exemple avec gnuplot :

    gnuplot
    gnuplot> set term x11
    gnuplot> plot 'e_r_001' w l, 'n_r_001' w l
    gnuplot> plot 'e_u_001' w l, 'n_u_001' w l
    gnuplot> plot 'e_p_001' w l, 'n_p_001' w l
    gnuplot> quit

On peut ainsi comparer la solution numérique à la solution exacte.

Comparaison numérique/exacte de la densité à \(t=0.25s\) , maillage uniforme 400 cellules

Comparaison numérique/exacte de la vitesse à \(t=0.25s\), maillage uniforme 400 cellules

### 3.2.2 Rupture de barrage

On considère ici un cas test 2D de rupture de barrage. Il s’agit d’un cas très largement utilisé et documenté. Dans point de vue expérimental, on peut par exemple s’appuyer sur les expériences réalisées par Martin and Moyce (1952). On considère un domaine de calcul \(x\in[0.,4.]\) et \(y\in[0.,3.]\). Les bords du domaine sont caractérisés par des conditions aux limites de type “miroir” ou Neumann nul. À l’instant initial, le domaine est constitué d’air, hormis une zone définie par \(x\in[0.,2.]\) et \(y\in[0.,3.]\) où l’eau est présente. Le maillage maître est constitué de \(40 \times 30\) blocs qui sont raffinés à un niveau 4 autour de l’interface air / eau. Pour améliorer la qualité de l’interface, on pénalise le calcul par un coefficient de raidissement d’interface de \(0.02\).
Le fichier de données correspondant au cas traité est fourni ci-après :

Tip 3.6: Fichier de données ex2_bf_Dam.inp

    !-----------------------------------------------------------------------
    ! Data file in CERF format
    ! Each block of data is characterized in the first line by a header (4 characters)
    ! then a code (1 integer) and an argument (1 integer)
    ! PHYS MAME INIT BATH COND MESH OBST NUME
    !-----------------------------------------------------------------------

    !-----------------------------------------------------------------------
    ! Reading physical parameters
    ! Header PHYS
    ! 1 keyword (4 characters) and a code (1 integer)
    ! for now we have the same values for all domains
    !-----------------------------------------------------------------------
    ! MODE : physical model
    !        1/ Shallow water
    !        2/ Euler isothermal two-phases               (default 0.)
    !---- for isothermal Euler model
    ! GPES : Gravity Orientation  0/none, 1/x, 2/y, 3/z   (default 2.)
    ! PREF : isothermal pressure of reference             (default 1.d5)
    ! CSSM : Mixture sound velocity (isothermal)          (default 20.)
    ! REFA : Isothermal reference density for Air         (default 1.)
    ! REFW : Isothermal reference density for Water       (default 1000.)
    ! SHAR : isothermal interface sharpening parameter
    !        recommended value 0.01                       (defaut 0.)
    !---- for Shallow water model
    ! WLTR : Value to truncate the Water level            (default 0.0001)
    ! FRMO : Friction Model
    ! 0/none, 1/Manning-Strickler, 2/Darcy-Weisbach       (default 0.)
    ! BAGR : Bathymetry gradient 0/ Given , 1/ numerical approximation
    !-----------------------------------------------------------------------
    PHYS 0
    MODE 3.
    GPES 2
    SHAR 0.02
    !-----------------------------------------------------------------------
    ! Reading the MAster MEsh
    ! Header  MAME
    !-----------------------------------------------------------------------
    ! code 1 -> creation of a mesh
    !
    !           argument: 0 type "1D shock tube"
    ! next line: xmin,xmax
    ! next line: nx
    !
    !           argument: 1 type "2d box"
    ! next line: xmin,xmax,ymin,ymax
    ! next line: nx,ny
    !
    !           argument: 2 type "3d box"
    ! next line: xmin,xmax,ymin,ymax,zmin,zmax
    ! next line: nx,ny,nz
    !
    ! code 2 -> Reading an ascii file format GMSH V2.2
    ! next line: file name
    !-----------------------------------------------------------------------
    MAME    1    1
    0. 4.  0. 3.
    40 30
    !-----------------------------------------------------------------------
    ! Reading Initial conditions
    ! Header INIT
    !-----------------------------------------------------------------------
    ! code 0 -> Definition for a shock tube
    !           primitive variables on left
    !           primitive variables on right
    !           argument: 0
    ! code 1 -> user function
    !           argument: Function's Number
    ! code 2 -> Definition using areas
    !           argument: number of areas
    !           and for each area 2 lignes
    !           - xmin,xmax,ymin,ymax,zmin,zmax
    !           - nvar values ( primitive variables)
    ! Shallow water model: water level, x velocity, y velocity
    ! Air-water flow model: density, x velocity, y velocity,z velocity, pressure
    !-----------------------------------------------------------------------
    INIT    2    2
    -1. 4.1 -1. 3.1 -1. 1.
    1. 0. 0. 0.  1.e5 0.
    -1. 1. -1. 2. -1. 1.
    1000. 0. 0. 0.  1.e5 0.

    !-----------------------------------------------------------------------
    ! Reading boundary condition
    ! (maximum 10 for the moment)
    ! Header COND
    !-----------------------------------------------------------------------
    !
    ! code 0 -> Description in the case  "shock tube" : Nothing
    !
    ! code 1 -> Description in the case  "2d box" (at least 1 line per type) :
    !           kind of boundary condition in xmin, xmax, ymin, ymax
    !
    ! code 2 -> Description in the case  "3d box" (at least 1 line per type) :
    !           kind of boundary condition in xmin, xmax, ymin, ymax, zmin,zmax
    !
    ! code 3 -> Definition per zone for  GMSH mesh
    !           argument: number of zones defining a boundary condition
    !           number of zone, kind of boundary condition
    !
    ! kind of boundary condition:
    ! 0 outlet (we copy)
    ! 1 miror
    ! 2 Dirichlet condition (if 1.e20 then copy), then next line nvar values
    ! 3 user's function , then next line function number
    ! 4 Dirichlet condition on conservative variables  (if 1.e20 then copy),
    !   then next line nvar values
    ! 999 Non used
    !-----------------------------------------------------------------------
    COND    1    0
    1
    1
    1
    1

    !-----------------------------------------------------------------------
    ! Reading mesh parameters in order to define the mesh according to
    ! the Master mesh is the one represented by the domain "0".
    ! Header MESH
    !-----------------------------------------------------------------------
    !
    ! code 0 -> default value
    !
    ! 1 keyword (4 characters) and a real
    ! NBDS : NumBer of DomainS                         (default 001.)
    ! MARL : MAximal Refinement Level                  (default 000.)
    ! COPA : Mesh COarsening PArameter          0<..<1 (default 0.002)
    ! REPA : Mesh REfinement PArameter          0<..<1 (default 0.02)
    ! FDRA : Function defining the mesh refinement to be applied
    ! NZRA : Number of zones where initial blocks are refined (default 1.)
    !        then for each zone
    !        nrb,xmin,xmax,ymin,ymax,zmin,zmax
    !-----------------------------------------------------------------------
    MESH 0
    NBDS 8
    MARL 4
    COPA 0.9
    REPA 3
    NZRA 3
    0 -1. 4.1 -1. 3.1 -1. 1.
    4 -1. 1.5 -1. 2.25 -1. 1.
    0 -1. .75 -1. 1.75 -1. 1.

    !-----------------------------------------------------------------------
    ! Reading numerical parameter
    ! Header  NUME
    !-----------------------------------------------------------------------
    ! 1 keyword (4 characters) then a real
    ! code 0 -> default value
    !
    ! TINT : Time interval to compute                    (default 0.)
    ! CCFL : cfl condition                               (default 0.9)
    ! TDOR : Time discretization order 1 ou 2         (defaut 1.)
    ! SDOR : Space discretization order 1 ou 2        (defaut 1.)
    ! LIMI : kind of limiter  0/Barth, 1/Cartesian       (default 0.)
    ! SAVE : kind of backup, three-digit number          (default 001.)
    !        decade  : shock tube backup 0/no 1/yes
    !        unit    : binary backup     0/no 1/yes
    ! NPRO : number of probe  < 10                      (defaut 000.)
    !        then coordinates x,y,z of the probe (1 probe per line)
    !-----------------------------------------------------------------------
    NUME 0
    TINT 0.01
    CCFL .9
    SDOR 2.
    TDOR 2.
    SAVE 1.

Les vitesses attendues sont de l’ordre de quelques mètres par seconde. Donc, comme les blocs ont une taille de \(0.1 m\), on peut estimer que l’intervalle de temps entre chaque remaillage doit être de \(0.01 s\). Pour accélérer le calcul on décompose le maillage en 8 domaines ( ou encore 8 process Mpi). Pour effectuer un calcul jusqu’à \(T=0.8 s\) on peut utiliser le script suivant :

    #!/bin/bash
    ./cerf_input ex2_bf_Dam.inp
    ./cerf2tec
    mv cerfout.dat d_0.dat
    mv cerfamr.dat d_mame_0.dat
    for i in $(seq 1 80 )
    do
    mpirun -np 8 ./cerf
    ./cerf2tec
    mv cerfout.dat d_$i.dat
    ./cerf_amr
    mv cerfamr.dat d_mame_$i.dat
    done

Sur l’animation ci-après Figure 3.6, la densité du mélange est représentée (Bleu pour l’eau et rouge pour l’air). On vérifie que l’interface est “fine”, c’est-à-dire sur 3-4 cellules.

Figure 3.6: Rupture de barrage, densité

### 3.2.3 Chute d’un cylindre sur un plan d’eau

On considère ici un cas test 2D de chute d’un objet solide sur un plan d’eau. Ce problème d’interaction fluide-structure est traité ici avec l’outil de “fluide-rigide” de CERF (Frederic Golay (2018)). Il s’agit d’un cas très largement utilisé et documenté (voir par exemple H. Sun and Faltinsen (2006), P. Sun, Ming, and Zhang (2015) ). D’un point de vue expérimental, on peut par exemple s’appuyer sur les expériences réalisées par Greenhow and Lin (2015). On considère un domaine de calcul \(x\in[0.,2.]\) et \(y\in[0.,2.]\). Les bords du domaine sont caractérisés par des conditions aux limites de type “miroir” ou Neumann nul. À l’instant initial, le domaine est constitué d’air, hormis une zone définie par \(x\in[0.,2.]\) et \(y\in[0.,0.65]\) où l’eau est présente. Le solide rigide est défini par un maillage au format .OBJ. Les cellules à l’intérieur de cette enveloppe auront un comportement de solide rigide libre, soumis à la pesanteur et de densité \(500 kg/m^3\). Le maillage maître est constitué de \(40 \times 40\) blocs qui sont raffinés à un niveau \(4\) autour de l’interface air / eau et autour du cylindre dont le centre est positionné à \(0.5m\) de la surface du plan d’eau comme illustré sur la figure Figure 3.7. Pour améliorer la qualité de l’interface, on pénalise le calcul par un coefficient de raidissement d’interface de \(0.01\).

Figure 3.7: Descriptif du cas test : chute d’un cylindre sur un plan d’eau

Le fichier de données correspondant au cas traité est fourni ci-après :

Tip 3.7: Fichier de données ex3_bf_cyl.inp

    !-----------------------------------------------------------------------
    ! Data file in CERF format
    ! Each block of data is characterized in the first line by a header (4 characters)
    ! then a code (1 integer) and an argument (1 integer)
    ! PHYS MAME INIT BATH COND MESH OBST NUME
    !-----------------------------------------------------------------------

    !-----------------------------------------------------------------------
    ! Reading physical parameters
    ! Header PHYS
    ! 1 keyword (4 characters) and a code (1 integer)
    ! for now we have the same values for all domains
    !-----------------------------------------------------------------------
    ! MODE : physical model
    !        1/ Shallow water
    !        2/ Euler isothermal two-phases               (default 0.)
    !---- for isothermal Euler model
    ! GPES : Gravity Orientation  0/none, 1/x, 2/y, 3/z   (default 2.)
    ! PREF : isothermal pressure of reference             (default 1.d5)
    ! CSSM : Mixture sound velocity (isothermal)          (default 20.)
    ! REFA : Isothermal reference density for Air         (default 1.)
    ! REFW : Isothermal reference density for Water       (default 1000.)
    ! SHAR : isothermal interface sharpening parameter
    !        recommended value 0.01                       (defaut 0.)
    !---- for Shallow water model
    ! WLTR : Value to truncate the Water level            (default 0.0001)
    ! FRMO : Friction Model
    ! 0/none, 1/Manning-Strickler, 2/Darcy-Weisbach       (default 0.)
    ! BAGR : Bathymetry gradient 0/ Given , 1/ numerical approximation
    !-----------------------------------------------------------------------
    PHYS 0
    MODE 2.
    GPES 2
    SHAR 0.01
    !-----------------------------------------------------------------------
    ! Reading the MAster MEsh
    ! Header  MAME
    !-----------------------------------------------------------------------
    ! code 1 -> creation of a mesh
    !
    !           argument: 0 type "1D shock tube"
    ! next line: xmin,xmax
    ! next line: nx
    !
    !           argument: 1 type "2d box"
    ! next line: xmin,xmax,ymin,ymax
    ! next line: nx,ny
    !
    !           argument: 2 type "3d box"
    ! next line: xmin,xmax,ymin,ymax,zmin,zmax
    ! next line: nx,ny,nz
    !
    ! code 2 -> Reading an ascii file format GMSH V2.2
    ! next line: file name
    !-----------------------------------------------------------------------
    MAME    1    1
    0. 2.  0. 2.
    40 40
    !-----------------------------------------------------------------------
    ! Reading Initial conditions
    ! Header INIT
    !-----------------------------------------------------------------------
    ! code 0 -> Definition for a shock tube
    !           primitive variables on left
    !           primitive variables on right
    !           argument: 0
    ! code 1 -> user function
    !           argument: Function's Number
    ! code 2 -> Definition using areas
    !           argument: number of areas
    !           and for each area 2 lignes
    !           - xmin,xmax,ymin,ymax,zmin,zmax
    !           - nvar values ( primitive variables)
    ! Shallow water model: water level, x velocity, y velocity
    ! Air-water flow model: density, x velocity, y velocity,z velocity, pressure
    !-----------------------------------------------------------------------
    INIT    2    2
    -1. 2.1 -1. 2.1 -1. 1.
    1. 0. 0. 0.  1.e5 0.
    -1. 2. -1. 0.65 -1. 1.
    1000. 0. 0. 0.  1.e5 0.

    !-----------------------------------------------------------------------
    ! Reading boundary condition
    ! (maximum 10 for the moment)
    ! Header COND
    !-----------------------------------------------------------------------
    !
    ! code 0 -> Description in the case  "shock tube" : Nothing
    !
    ! code 1 -> Description in the case  "2d box" (at least 1 line per type) :
    !           kind of boundary condition in xmin, xmax, ymin, ymax
    !
    ! code 2 -> Description in the case  "3d box" (at least 1 line per type) :
    !           kind of boundary condition in xmin, xmax, ymin, ymax, zmin,zmax
    !
    ! code 3 -> Definition per zone for  GMSH mesh
    !           argument: number of zones defining a boundary condition
    !           number of zone, kind of boundary condition
    !
    ! kind of boundary condition:
    ! 0 outlet (we copy)
    ! 1 miror
    ! 2 Dirichlet condition (if 1.e20 then copy), then next line nvar values
    ! 3 user's function , then next line function number
    ! 4 Dirichlet condition on conservative variables  (if 1.e20 then copy),
    !   then next line nvar values
    ! 999 Non used
    !-----------------------------------------------------------------------
    COND    1    0
    1
    1
    1
    1

    !-----------------------------------------------------------------------
    ! Reading mesh parameters in order to define the mesh according to
    ! the Master mesh is the one represented by the domain "0".
    ! Header MESH
    !-----------------------------------------------------------------------
    !
    ! code 0 -> default value
    !
    ! 1 keyword (4 characters) and a real
    ! NBDS : NumBer of DomainS                         (default 001.)
    ! MARL : MAximal Refinement Level                  (default 000.)
    ! COPA : Mesh COarsening PArameter          0<..<1 (default 0.002)
    ! REPA : Mesh REfinement PArameter          0<..<1 (default 0.02)
    ! FDRA : Function defining the mesh refinement to be applied
    ! NZRA : Number of zones where initial blocks are refined (default 1.)
    !        then for each zone
    !        nrb,xmin,xmax,ymin,ymax,zmin,zmax
    !-----------------------------------------------------------------------
    MESH 0
    NBDS 10
    MARL 4
    COPA 1
    REPA 2
    NZRA 4
    0 -1. 2.1 0. 2.1 -1. 1
    4 -1. 2.1 0.6 0.7 -1. 1
    0 -1. 2.1 0. 0.6 -1. 1
    4 0.9 1.1  1.0 1.3  -1. 1

    !-----------------------------------------------------------------------
    ! Reading rigid obstacle parameters
    ! Header OBST
    !-----------------------------------------------------------------------
    !
    ! code 0 -> default value
    !           argument: number of obstacles
    ! 1 line per obstacle: file name, density, kind of motion
    ! (  -1: fixed, 0: free, >0 : function number imposing the movement
    !
    !-----------------------------------------------------------------------
    OBST    0    1
    ex3_bf_cyl.obj 500. 0

    !-----------------------------------------------------------------------
    ! Reading numerical parameter
    ! Header  NUME
    !-----------------------------------------------------------------------
    ! 1 keyword (4 characters) then a real
    ! code 0 -> default value
    !
    ! TINT : Time interval to compute                    (default 0.)
    ! CCFL : cfl condition                               (default 0.9)
    ! TDOR : Time discretization order 1 ou 2         (defaut 1.)
    ! SDOR : Space discretization order 1 ou 2        (defaut 1.)
    ! LIMI : kind of limiter  0/Barth, 1/Cartesian       (default 0.)
    ! SAVE : kind of backup, three-digit number          (default 001.)
    !        decade  : shock tube backup 0/no 1/yes
    !        unit    : binary backup     0/no 1/yes
    ! NPRO : number of probe  < 10                      (defaut 000.)
    !        then coordinates x,y,z of the probe (1 probe per line)
    !-----------------------------------------------------------------------
    NUME 0
    TINT 0.01
    CCFL .9
    SDOR 2.
    TDOR 2.
    TIME 1
    SAVE 1.

Les vitesses attendues sont de l’ordre de quelques mètres par seconde. Donc, comme les blocs ont une taille de \(0.05 m\), on peut estimer que l’intervalle de temps entre chaque remaillage doit être de \(0.01 s\). Pour accélérer le calcul on décompose le maillage en 10 domaines ( ou encore 10 process Mpi). Pour effectuer un calcul jusqu’à \(T=2 s\) on peut utiliser le script suivant :

    #!/bin/bash
    ./cerf_input ex3_bf_cyl.inp
    ./cerf2tec
    mv cerfout.dat c_000.dat
    mv cerfamr.dat c_mame_000.dat
    cp ex3_bf_cyl.obj cs_000.obj

    for i in $(seq 1 9 )
    do
    mpirun -np 10 ./cerf
    ./cerf2tec
    mv cerfout.dat c_00$i.dat
    ./cerf_amr
    mv cerfamr.dat c_mame_00$i.dat
    mv solide_001.obj cs_00$i.obj
    done
    for i in $(seq 10 99 )
    do
    mpirun -np 10 ./cerf
    ./cerf2tec
    mv cerfout.dat c_0$i.dat
    ./cerf_amr
    mv cerfamr.dat c_mame_0$i.dat
    mv solide_001.obj cs_0$i.obj
    done
    for i in $(seq 100 200 )
    do
    mpirun -np 10 ./cerf
    ./cerf2tec
    mv cerfout.dat c_$i.dat
    ./cerf_amr
    mv cerfamr.dat c_mame_$i.dat
    mv solide_001.obj cs_$i.obj
    done

Sur les images ci-après, la densité du mélange est représentée (Bleu pour l’eau et rouge pour l’air). On vérifie que l’interface est “fine”, c’est-à-dire sur 3-4 cellules et que le maillage dynamique “suit” bien les zones d’intérêt. La première figure Figure 3.8, à \(t=0.3 s\), correspond à l’instant où le cylindre touche la surface de l’eau à une vitesse de l’ordre de \(3m/s\), ce qui est en accord avec la théorie (\(\sqrt{2gh}\)) et l’expérience. La seconde figure Figure 3.9, à \(t=0.45 s\), illustre la génération de la vague qui se crée lors de l’entrée du cylindre dans l’eau.

Figure 3.8: Chute d’un cylindre, densité et maillage à \(t=0.3s\)
Figure 3.9: Chute d’un cylindre, densité et maillage à \(t=0.45s\)

Sur l’animation ci-après Figure 3.10, la densité du mélange est représentée (Bleu pour l’eau et rouge pour l’air). Le cylindre est blanc.

Figure 3.10: Chute d’un cylindre: densité

Sur la figure ci-après Figure 3.11, on compare la trajectoire verticale du cylindre au cours du temps avec les résultats expérimentaux de Greenhow and Lin (2015). On constate un bon accord entre les deux approches.

Figure 3.11: Trajectoire verticale du cylindre au cours du temps: comparaison de la simulation numérique (trait plein) avec l’expérience (point)
Delestre, Olivier, Carine Lucas, Pierre-Antoine Ksinant, Frédéric Darboux, Christian Laguerre, T.-N.-Tuoi Vo, François James, and Stéphane Cordier. 2013. “SWASHES: A Compilation of Shallow Water Analytic Solutions for Hydraulic and Environmental Studies.” International Journal for Numerical Methods in Fluids 72 (3): 269–300. https://doi.org/https://doi.org/10.1002/fld.3741.
Faccanoni, Gloria, and Anne Mangeney. 2013. “Exact Solution for Granular Flows.” International Journal for Numerical and Analytical Methods in Geomechanics 37 (10): 1408–33. https://doi.org/https://doi.org/10.1002/nag.2124.
Golay, Frederic. 2018. “Monolithic Fluid Structure Interaction Model, Application to Water Entry Problem.” In Topical Problems of Fluid Mechanics 2018. Prague, Czech Republic: Institute of Thermomechanics AS CR v.v.i. https://doi.org/10.14311/TPFM.2018.016.
Golay, Frédéric, and Philippe Helluy. 2007. “Numerical Schemes for Low Mach Wave Breaking.” International Journal of Computational Fluid Dynamics 21 (2): 69–86. https://doi.org/10.1080/10618560701343382.
Greenhow, Martin, and and Lin. 2015. “Greenhow m. And Lin w-m."nonlinear Free Surface Effects: Experiments and Theory" MIT Internal Report 83-19 (1983),” June. https://doi.org/10.13140/RG.2.1.4129.2964.
Martin, J. C., and W. J. Moyce. 1952. “Part IV. An Experimental Study of the Collapse of Liquid Columns on a Rigid Horizontal Plane.” Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 244 (882): 312–24. http://www.jstor.org/stable/91519.
Sun, Hui, and Odd M. Faltinsen. 2006. “Water Impact of Horizontal Circular Cylinders and Cylindrical Shells.” Applied Ocean Research 28 (5): 299–311. https://doi.org/https://doi.org/10.1016/j.apor.2007.02.002.
Sun, Pengnan, Furen Ming, and Aman Zhang. 2015. “Numerical Simulation of Interactions Between Free Surface and Rigid Body Using a Robust SPH Method.” Ocean Engineering 98: 32–49. https://doi.org/https://doi.org/10.1016/j.oceaneng.2015.01.019.