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
./cerfComme 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> quitOn peut ainsi comparer la solution numérique à la solution exacte.
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 :
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
./cerf2tecOn 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).
On peut également visualiser la solution numérique avec Visit ou Paraview:
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
doneComme précédemment, on peut comparer avec la solution exacte.
On peut également visualiser la solution numérique avec Visit ou Paraview :
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).
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
./cerfComme 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> quitOn peut ainsi comparer la solution numérique à la solution exacte.
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
doneSur 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.
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\).
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
doneSur 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.
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.
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.
