3 Documented examples
3.1 Examples with the Saint-Venant model
In this chapter, we present the implementation of a few classic flow test cases for the Saint-Venant model. The cases covered are as follows:
- Solving a Riemann problem
- Hydraulic jump
- Friction, McDonald test case
- …
3.1.1 Riemann problem
We consider the classic 1D dam-break test case for the Saint-Venant model on a flat bottom. We use the test case from Faccanoni and Mangeney (2013) (test case 5). The domain is defined by \(x\in[-2.,2.]\) and \(t\in[0,0.5]\). The initial conditions are given by: \[
\overrightarrow{w_p}(x,0) =\left\{
\begin{array}{ll}
\overrightarrow{w_p}_L = <h_L,u_L>^T = <0.0046,0.>^T &\mbox{ if } x<0 \\
\overrightarrow{w_p}_R = <h_R,u_R>^T = <0.1446,0.>^T &\mbox{ if } x>0
\end{array}
\right.
\]
This is therefore a dam-break case without a dry zone, with a shock and a rarefaction wave, for which the exact solution can be computed using the Riemann solver.
First, we consider a uniform, fixed 1D mesh of 400 cells over a single domain. The data file corresponding to this test case is provided below:
Tip 3.1: Data file 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.
To run the computation, since there is only one domain to process over a single time interval:
./cerf_input ex1_sv_Riemann.inp
./cerfSince the data file ex1_sv_Riemann.inp requests output in 1D format (parameter SAVE=11, section NUME), 1D output files have been generated. Files for the numerical solution in water depth and velocity ‘n_h_001’ and ‘n_u_001’, and files for the exact solution in water depth and velocity: ‘e_h_001’ and ‘e_u_001’. These files can be visualized, for example, with 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> quitThe numerical solution can thus be compared to the exact solution.
To test the code with a “dry” zone (Faccanoni and Mangeney (2013), test case 1), the same problem can for example be revisited with a different initial condition by modifying the INIT section of the data file ex1_sv_Riemann.inp:
INIT 0 0
0. 0. 0.
0.1446 0. 0. This then gives:
3.1.2 Hydraulic jump
We consider here a steady-state test case, without friction, with a non-flat bottom. The analytical solution is taken from Delestre et al. (2013). On a domain \(x\in[0.,25.]\), we consider a topography defined by:
\[z_b(x) = \left\{ \begin{array}{ll} 0.2 - 0.05 \left( x - 10 \right)^2 & \mbox{ if } 8 < x < 12 \\ 0. & \mbox{ otherwise } \end{array}\right. \]
We consider the case of transcritical flow with a shock, characterized by the following boundary conditions:
- On the left, \(h = 0.33 m\) is imposed and the discharge is left free.
- On the right, \(hu = 0.18 m^2/s\) is imposed and the water depth is left free.
As recommended by Delestre et al. (2013), the initial conditions are given by: \(h + z_b =0\) and \(u=0\).
The computation is carried out on a fixed, regular mesh of 2500 cells distributed over 8 domains. In order to reach a steady state, the computation is run up to \(T=100s\). The data file corresponding to this test case is provided below Tip 3.2:
Tip 3.2: Data file 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.
To run the computation, since there is only one domain to process over a single time interval, but distributed over 8 domains, the following commands are used:
./cerf_input ex2_sv_Bump.inp
mpirun -np 8 ./cerf
./cerf2tecThe analytical solution to this problem can be obtained from the Swashes website. In the doc directory of CERF, the analytical solution of \(\eta = h + z_b\) on 2500 points in Tecplot format is provided for reference. The numerical solution can then be compared to the exact solution in the figure below (Figure 3.1).
The numerical solution can also be visualized with Visit or Paraview:
In figure Figure 3.2, the isovalues of the water level are shown at \(t=100s\). Note that the mesh has been modified for post-processing purposes. The coordinates of the nodes at the base of the hexahedra correspond to \(z_b\), while the nodes on the upper part of the hexahedra correspond to \(z_b+h\). This makes it possible to visualize the topography and the water level at the same time.
This computation can also be performed using automatic mesh refinement. To do this, the data file is modified accordingly ex2_sv_Bump.inp. A mesh consisting of only 250 blocks is then created, and mesh refinement is allowed only up to level 2. All other computation parameters are kept the same as for the fixed mesh.
Tip 3.3: Data file 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.
To run the computation with mesh refinement, the following script is used:
#!/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
doneAs before, the results can be compared with the exact solution.
The numerical solution can also be visualized with Visit or Paraview:
In both Figure 3.3 and Figure 3.4, it can be seen that mesh refinement made it possible to reduce the number of cells while maintaining good accuracy in the numerical solution.
3.1.3 Friction, McDonald test case
We consider here a steady-state test case, with friction, on a non-flat bottom. The semi-analytical solution is taken from Delestre et al. (2013) (case 3.2.1).
On a domain \(x\in[0.,1000.]\), we consider a subcritical flow with a constant discharge at the inlet of \(q = h u = 2 m^2/s\) and a prescribed water depth at the outlet \(h_{e}(1000)\). The domain is initially dry, with \(h = 0\) and \(u = 0\). A Manning friction with \(n = 0.033\) is imposed.
The water depth is given by: \[ 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)\]
Using Equation 1.6, in the steady-state case, the topography can then be recovered by integrating: \[\frac{\partial z_b}{\partial x} = \left( \frac{q^2}{g h^3} - 1 \right) \frac{\partial h}{\partial x} - S_f\]
The computation is carried out on a fixed, regular mesh of 1000 cells. In order to reach a steady state, the computation is run up to \(T=1500s\). The data file corresponding to this test case is provided below Tip 3.4:
Tip 3.4: Data file 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.
To run the computation, in order to visualize the evolution towards the steady state, the computation time is split into \(30\) steps using the script below:
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
The analytical solution to this problem can be obtained from the Swashes website. In the doc directory of CERF, the analytical solution on 10000 points in Tecplot format is provided for reference. The numerical solution can then be compared to the exact solution in the figure below (Figure 3.5).
3.2 Examples with the Eulerian two-fluid model
In this chapter, we present the implementation of a few classic test cases for two-fluid flows. The cases covered are as follows:
- Solving a Riemann problem
- Dam break
- A cylinder falling onto a body of water
- …
3.2.1 Riemann problem
In order to validate the implementation of the model, we consider the classic “shock tube” test case with the isothermal two-fluid model. We use the test case from Frédéric Golay and Helluy (2007). The domain is defined by \(x\in[-2.,2.]\) and \(t\in[0,0.25]\). The initial conditions are given by: \[ \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{ if } x<0 \\ \overrightarrow{w_p}_R = <\rho_R,u_R,p_R>^T = <0.1446,\;0.,\;1. \;10^5>^T &\mbox{ if } x>0 \end{array} \right. \]
We consider a uniform, fixed 1D mesh of \(400\) cells over a single domain. The data file corresponding to this test case is provided below:
Tip 3.5: Data file 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.
To run the computation, since there is only one domain to process over a single time interval:
./cerf_input ex1_bf_Riemann.inp
./cerfSince the data file ex1_bf_Riemann.inp requests output in 1D format (parameter SAVE=11, section NUME), 1D output files have been generated. Files for the numerical solution in density, velocity, and pressure ‘n_r_001’, ‘n_u_001’, and ‘n_p_001’, and files for the exact solution in density, velocity, and pressure: ‘e_r_001’, ‘e_u_001’, and ‘e_p_001’. These files can be visualized, for example, with 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> quitThe numerical solution can thus be compared to the exact solution.
3.2.2 Dam break
We consider here a 2D dam-break test case. This is a widely used and well-documented case. From an experimental standpoint, one can rely for example on the experiments carried out by Martin and Moyce (1952). We consider a computational domain \(x\in[0.,4.]\) and \(y\in[0.,3.]\). The domain boundaries are characterized by “mirror”-type or zero-Neumann boundary conditions. At the initial instant, the domain consists of air, except for a region defined by \(x\in[0.,2.]\) and \(y\in[0.,3.]\) where water is present. The master mesh consists of \(40 \times 30\) blocks, which are refined to level 4 around the air/water interface. To improve the quality of the interface, the computation is penalized using an interface stiffening coefficient of \(0.02\).
The data file corresponding to this test case is provided below:
Tip 3.6: Data file 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.
The expected velocities are on the order of a few meters per second. Therefore, since the blocks have a size of \(0.1 m\), the time interval between each remeshing can be estimated to be \(0.01 s\). To speed up the computation, the mesh is decomposed into 8 domains (or equivalently 8 MPI processes). To run a computation up to \(T=0.8 s\), the following script can be used:
#!/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
doneIn the animation below Figure 3.6, the density of the mixture is shown (blue for water and red for air). The interface can be seen to be “sharp”, i.e., spread over 3-4 cells.
3.2.3 A cylinder falling onto a body of water
We consider here a 2D test case of a solid object falling onto a water surface. This fluid-structure interaction problem is treated here using CERF’s “rigid-fluid” tool (Frederic Golay (2018)). This is a widely used and well-documented case (see for example H. Sun and Faltinsen (2006), P. Sun, Ming, and Zhang (2015)). From an experimental standpoint, one can rely for example on the experiments carried out by Greenhow and Lin (2015). We consider a computational domain \(x\in[0.,2.]\) and \(y\in[0.,2.]\). The domain boundaries are characterized by “mirror”-type or zero-Neumann boundary conditions. At the initial instant, the domain consists of air, except for a region defined by \(x\in[0.,2.]\) and \(y\in[0.,0.65]\) where water is present. The rigid solid is defined by a mesh in .OBJ format. The cells inside this envelope behave as a free rigid solid, subject to gravity, with a density of \(500 kg/m^3\). The master mesh consists of \(40 \times 40\) blocks, which are refined to level \(4\) around the air/water interface and around the cylinder, whose center is positioned \(0.5m\) from the surface of the water, as illustrated in figure Figure 3.7. To improve the quality of the interface, the computation is penalized using an interface stiffening coefficient of \(0.01\).
The data file corresponding to this test case is provided below:
Tip 3.7: Data file 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.
The expected velocities are on the order of a few meters per second. Therefore, since the blocks have a size of \(0.05 m\), the time interval between each remeshing can be estimated to be \(0.01 s\). To speed up the computation, the mesh is decomposed into 10 domains (or equivalently 10 MPI processes). To run a computation up to \(T=2 s\), the following script can be used:
#!/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
doneIn the images below, the density of the mixture is shown (blue for water and red for air). The interface can be seen to be “sharp”, i.e., spread over 3-4 cells, and the dynamic mesh can be seen to “follow” the regions of interest closely. The first figure Figure 3.8, at \(t=0.3 s\), corresponds to the instant when the cylinder touches the water surface at a velocity of about \(3m/s\), which is consistent with theory (\(\sqrt{2gh}\)) and with experiment. The second figure Figure 3.9, at \(t=0.45 s\), illustrates the generation of the wave created as the cylinder enters the water.
In the animation below Figure 3.10, the density of the mixture is shown (blue for water and red for air). The cylinder is shown in white.
In the figure below Figure 3.11, the vertical trajectory of the cylinder over time is compared with the experimental results of Greenhow and Lin (2015). Good agreement is observed between the two approaches.
