Gerris par un nul

But du jeu: calculer l'écoulement dans un tuyau avec une constriction. Comparer à FreeFem++ et aux calculs de l'article Lagrée Chouly

Rappel on résout NS, on adimensionne avec la largeur du canal

...
Le Reynolds est (la hauteur du canal est notée H0 dans le code gerris par la suite):
Re = U0 H0/nu
La vitesse avec dimension est u=U0(y/H0) (1-y/H0)
le gradient de pression sera -dp/dx = 2 mu U0^2/H0^2
la pression sans dimension (rhoU0^2) est donc 2 /Re avec Re= U0 H0/nu

Ce choix de Reynolds permet d'avoir une valeur modérée pour Re =750 et de comparer avec les modèle asymptotiques.
Notre choix donne une pente de vitesse unité. Construit avec la vitesse max le Reynolds est Rm = Re/4.
Construit avec la vitesse moyenne le Reynolds n'est que Ra = Re/6.

1) Première étape

compiler gerris!

2) Premier Exemple

En premier exemple, on calcule Poiseuille sur un carré de taille 1 en x et en y. En effet gerris travaille par défaut dans la boîte unité.

Le code le plus simple serait:

######################################################################
# Poiseuille par PYL, sauver dans "pois0.gfs"
# 01/10

# definition de 1 boite avec 0 connections
1 0 GfsSimulation GfsBox GfsGEdge{} {

# valeur du Reynolds 250
SourceViscosity {} 1./250
# precision 2**(-4.) = 1/16=0.03
Refine 4
# valeurs qui vont sortir pour entrer dans gfsview
# tous les 10 pas de calcul
OutputSimulation { istep = 10 } stdout
}
#conditions aux limites
GfsBox {
# a gauche profil
          left = GfsBoundary {
           GfsBcDirichlet U   (1 - (y+0.5))*(y+0.5)
                      }
# a droite pression imposee, invariance de U en translation
          right = GfsBoundary {
            GfsBcDirichlet P 0
            GfsBcNeumann U 0
          }
# en haut vitesse nulle
          top = GfsBoundary {
            GfsBcDirichlet U 0
            GfsBcDirichlet V 0
          }
# en bas vitesse nulle
        bottom = GfsBoundary {
            GfsBcDirichlet U 0
            GfsBcDirichlet V 0
          }
}
# fin de fichier
######################################################################

on le sauve dans pois0.gfs
pour le lancer:
gerris2D pois0.gfs | gfsview2D

la fenêtre de gfsview s'ouvre
En cliquant sur "Linear " on a la pression, en cliquant sur "Vectors" on peut obtenir l'image suivante que l'on voit à droite

oh le beau Poiseuille

Ce programme ne s'arrète jamais, il n'affiche rien non plus dasn le terminal.
On en profite pour sauver la visualisation dans un fichier "vue0.gfv" qui contient la position de la fenêtre de vision et les champs que l'on veut tracer

La seconde version consiste à rajouter une initialisation non nulle (par défaut les champs semblent nuls), avec Init{}
puis une sortie texte donnant le temps en cours,
puis un test d'arrêt à convergence

#####################################################################
# Poiseuille par PYL, sauver dans "pois0.gfs"
# 01/10

# definition de 1 boite avec 0 connections
1 0 GfsSimulation GfsBox GfsGEdge{} {

# valeur du Reynolds 250
SourceViscosity {} 1./250
# precision 2**(-4.) = 1/16=0.03
Refine 4
# Poiseuille entre -0.5 et 0.5 donne au temps initial
Init {} { U = 0 }
# sortie tous les 10 pas de calculs du temps en cours
OutputTime { istep = 10 } stderr
# valeurs qui vont sortir pour entrer dans gfsview
# tous les 10 pas de calcul
OutputSimulation { istep = 10 } stdout
# arret lorsque la variation de U devient "petite"
EventStop { istep = 10 } U 0.5e-8 DU
}
#conditions aux limites
GfsBox {
# a gauche profil
          left = GfsBoundary {
           GfsBcDirichlet U   (1 - (y+0.5))*(y+0.5)
                      }
# a droite pression imposee, invariance de U en translation
          right = GfsBoundary {
            GfsBcDirichlet P 0
            GfsBcNeumann U 0
          }
# en haut vitesse nulle
          top = GfsBoundary {
            GfsBcDirichlet U 0
            GfsBcDirichlet V 0
          }
# en bas vitesse nulle
        bottom = GfsBoundary {
            GfsBcDirichlet U 0
            GfsBcDirichlet V 0
          }
}
# fin de fichier
######################################################################

Dans le fichier "vue0.gfv " il y a:

# GfsView 2D
View {
tx = 0 ty = 0
sx = 1 sy = 1 sz = 1
q0 = 0 q1 = 0 q2 = 0 q3 = 1
fov = 30
r = 0.3 g = 0.4 b = 0.6
res = 1
lc = 0.001
reactivity = 0.1
}
Vectors {
r = 0 g = 0 b = 0
shading = Constant
maxlevel = -1
} {
n.x = 0 n.y = 0 n.z = 1
pos = 0
} P {
amin = 1
amax = 1
cmap = Jet
} U V {
scale = 1.219
use_scalar = 0
}
Linear {
r = 1 g = 1 b = 1
shading = Constant
maxlevel = -1
} {
n.x = 0 n.y = 0 n.z = 1
pos = 0
} P {
amin = 1
amax = 1
cmap = Jet
} 0 {
reversed = 0
use_scalar = 1
}

On sauve dans pois0.gfs et l'autre dans vue0.gfv pour le lancer le calcul:
gerris2D pois0.gfs | gfsview2D vue0.gfv
la fenêtre de gfsview s'ouvre avec la vue précédente. Le calcul s'arrête tout seul au bout de 160 itérations.

3) Troisième Exemple

On continue sur Poiseuille sur un carré de taille différente de 1 en x et en y. En effet gerris travaille par défaut dans la boîte unité. On peut changer la taille en utilisant
PhysicalParams { L = 10 }
le carré fait alors 10x10 et non plus 1x1. Attention cela induit des bugs avec Gfs2oogl
On peut mettre la longueur en variable définie par "Define", de même pour la hauteur du canal que l'on note: H0.
on met une paroi en utilisant Solid
si on se donne une fonction f(x), on veut que la paroi supérieure soit f(x) alors on code
Solid(-y+f(x))
de même si on se donne une fonction g(x), on veut que la paroi inférieur soit g(x) alors on code
Solid( y -g(x))
le solide est là où la fonction est négative.

En pratique on met deux bosses exponentielles en H0-a*H0*exp(-pow((x)/H0,2)))
le canal est fermé de 2aH0, on prend a petit....

On enlève donc les conditions en top et bottom car l'adhérence est assurée par "Solid"
on a une nouvelle instruction RefineSolid

on rafine toujours, mais on a augmenté la longueur, donc le pas d'espace sera 2**(-6.)*10=0.15625 au plus
et par l' adaptation automatique de maillage sur la variation de vorticite dx=L0/2^7=10/128=0.08 au moins
AdaptVorticity { istep = 1 } { maxlevel = 7 cmax = 1e-2 }

Le code le plus simple serait:

# #

######################################################################
# Poiseuille par PYL, sauver dans "pois1.gfs"
# 01/10
# valeurs de dimension
Define L0 10
Define H0 1
Define a 0.01
Define Re 750.

1 0 GfsSimulation GfsBox GfsGEdge{} {
# attention la boite est L0xL0
PhysicalParams { L = L0 }
# on s'arrete au bout du temps 500 quoiqu'il arrive
Time {end = 500   }
# viscosite Re
SourceViscosity {} 1./Re
# le domaine est de longuer L0, les deux plans sont entre y=0 et y=H0
# paroi sup
Solid (-y+H0-a*H0*exp(-pow((x)/H0,2)))
#paroi inf
Solid ( y+0.-a*H0*exp(-pow((x)/H0,2)))
# au moins 64 cases, soir dx=0.16
Refine 6
# raffinement pour la paroi...
RefineSolid 9
# on a Poiseuille entre y=0 et y=H0
   Init {} { U = (1 - (y/H0))*(y/H0) }
# adaptation automatique de maillage sur la vorticite
# dx=L0/2^7=10/128=0.08
   AdaptVorticity { istep = 1 } { maxlevel = 7 cmax = 1e-2 }
OutputSimulation { istep = 10 } stdout
OutputTime { istep = 10 } stderr
# sortie des valeurs dans un fichier
OutputSimulation { step= 10 } VAL/knd-%f.gfs
OutputSimulation { start = end } VAL/end.gfs
#test d'arret
EventStop { istart = 20 istep = 1} U 5e-3 DU
}
#condition de vitesse en entree donnee
# et pression nulle en sortie ainsi que dx(U)=0
GfsBox { left = GfsBoundary {
           GfsBcDirichlet U   (1 - (y/H0))*(y/H0)}
         right = GfsBoundary {
           GfsBcDirichlet P 0
           GfsBcNeumann U 0 }
}

#####################################################################

on le sauve dans pois1.gfs
on crée un directory VAL
mkdir VAL

on lance avec l'option "-m" qui permet de tenir compte des "Define"
gerris2D -m pois1.gfs | gfsview2D vue1.gfv

voila la pression avec gnuplot

on l'obtient grâce à l'utilitaire qui permet de récupérer des valeurs:

gfs2oogl2D -c P -o -i < VAL/end.gfs > VAL/a.dat

on peut sortir le fichier de pression et le tracer dans gnuplot:
gnuplot> p'gerris_par_un_nul/VAL/val.dat' u 1:4,-(x-.5)/750*2*10
attention, il y a un bug, car la sortie de gfs2oogl2D donne une sortie entre -0.5 et 0.5. mais il faut en plus multiplier par 10

arrêt au t=70. Tous les fichiers intermédiaires sont dans VAL. On peut les relire avec gfsview, par exemple le dernier crée:
gfsview2D VAL/end.gfs
dans linear on peut remplacer le champ P par défaut par
P+(x+0.5)/750*2
on voit donc l'erreur faible induite par la bosse sur le poiseuille

4) Exemple quatrième

Le code est encore légèrement modifié, on rajoute une sauvegarde le long d'un "fil chaud" défini par une suite de points et on crée ce fichier dans un fichier de command qui va lancer le programme. On propose un autre fichier de commande pour aider à tracer avec gnuplot.

# #

######################################################################
# Poiseuille par PYL, sauver dans "pois1.gfs"
# 01/10
# valeurs de dimension
Define L0 10
Define H0 1
Define a 0.01
Define Re 750.

1 0 GfsSimulation GfsBox GfsGEdge{} {
# attention la boite est L0xL0
PhysicalParams { L = L0 }
# on s'arrete au bout du temps 500 quoiqu'il arrive
Time {end = 500   }
# viscosite Re
SourceViscosity {} 1./Re
# le domaine est de longuer L0, les deux plans sont entre y=0 et y=H0
# paroi sup
Solid (-y+H0-a*H0*exp(-pow((x)/H0,2)))
#paroi inf
Solid ( y+0.-a*H0*exp(-pow((x)/H0,2)))
# au moins 64 cases, soir dx=0.16
Refine 6
# raffinement pour la paroi...
RefineSolid 9
# on a Poiseuille entre y=0 et y=H0
   Init {} { U = (1 - (y/H0))*(y/H0) }
# adaptation automatique de maillage sur la vorticite
# dx=L0/2^7=10/128=0.08
   AdaptVorticity { istep = 1 } { maxlevel = 7 cmax = 1e-2 }
OutputSimulation { istep = 10 } stdout
OutputTime { istep = 10 } stderr
# sortie des valeurs dans un fichier
OutputSimulation { step= 10 } VAL/knd-%f.gfs
OutputSimulation { start = end } VAL/end.gfs
# sortie sur une ligne au milieu definie par le fichier "Paroi_bas.dat"
  OutputLocation { istep = 10 } VAL/Centre.data Paroi_bas.dat
#test d'arret
EventStop { istart = 20 istep = 1} U 5e-3 DU
}
#condition de vitesse en entree donnee
# et pression nulle en sortie ainsi que dx(U)=0
GfsBox { left = GfsBoundary {
           GfsBcDirichlet U   (1 - (y/H0))*(y/H0)}
         right = GfsBoundary {
           GfsBcDirichlet P 0
           GfsBcNeumann U 0 }
}

#####################################################################

on le sauve dans pois1.gfs
on crée un directory VAL
mkdir VAL

on lance avec l'option "-m" qui permet de tenir compte des "Define"
on a un autre fichier de vue vue1.gfv

gerris2D -m pois1.gfs | gfsview2D vue1.gfv

On met ce ceci dans run.sh

#!/bin/bash
#
#
export LANG=C
echo "lancement du calcul"
rm VAL/knd*

file="Paroi_bas.dat"
awk 'BEGIN{
      L0 = 10
      H0 = 1
      r = 0.*H0
      wd = L0/2.
      xc = -0.*L0/10.
      print -0.501 * L0   " " 0.0 " " 0.0
      for (x = -wd ; x <= wd; x += 0.01)
        print x " " (H0/2. + r*exp(-(x-xc)/H0*(x-xc)/H0)) " 0 ";
      print 0.501 * L0   " " 0.0 " " 0.0
    }' > $file

gerris2D -m pois1.gfs | gfsview2D vue1.gfv

echo "fin normale ?"

gfs2oogl2D -c P -o -i < VAL/end.gfs > VAL/a.dat

On peut mettre ensuite ceci dans le fichier fin.sh
#!/bin/bash
#
#
export LANG=C
echo "liigne de test dans z.txt"
echo "awk '{if (\$1>=tt) print \$0 }' VAL/Centre.data > VAL/Cf.data" > z.txt
echo "extrait le temps final mis dans tfinal.txt"
tfinal=`tail -n 1 VAL/Centre.data`
echo $tfinal > tfinal.txt
awk '{print $1}' tfinal.txt > dump.txt
echo `cat dump.txt`
rm tfinal.txt
echo "substit"
echo "s/tt/"`cat dump.txt`"/g">b.txt
sed `cat b.txt` z.txt>zz.txt
echo `head -n 1 VAL/Cf.data`
echo "transfer des dernieres donnes dans Cf.data"
source zz.txt
rm b.txt
rm z.txt
rm zz.txt

gfs2oogl2D -c P -o -i < `ls VAL/* | tail -n 1 ` > VAL/a.dat

le fichier run.sh crée un fichier Parois_bas.dat (le centre du canal en fait). Les champs sont sauvés le long des points définis dans ce fichier et mis dans le fichier
VAL/Centre.data
suivant:
# 1:T 2:X 3:Y 4:Z 5:P 6:Pmac 7:U 8:V 9:DU
Ils sont tous sauvés. Mais on ne veut que le temps final. Avec fin.sh on extrait le temps final et on met le dernier paquet dans le fichier Cf.data que l'on trace.
on trace la pression calculée et la pression théorique
p'VAL/Cf.data' u 2:5,-(x-5)/750*2
v reste petit
U est quasi égal à 0.25 et varie sur la bosse.

5) Exemple POISEUILLE

ET POISEUILLE dans tout ça, et oui, on a mis une bosse, en fait si on met pas de bosse. Ca marche pas, la pression oscille....

consulter le cas test http://gfs.sourceforge.net/tests/tests/poiseuille.html

6) Exemple FINAL, comparaison FreeFem++

On augmente maintenant la taille de la bosse, on modifie run.sh et pois2.gfs

# #

######################################################################
# Poiseuille par PYL, sauver dans "pois2.gfs"
# 01/10
# valeurs de dimension
Define L0 10
Define H0 1
Define a 0.1
Define Re 750.

1 0 GfsSimulation GfsBox GfsGEdge{} {
# attention la boite est L0xL0
PhysicalParams { L = L0 }
# on s'arrete au bout du temps 500 quoiqu'il arrive
Time {end = 500   }
# viscosite Re
SourceViscosity {} 1./Re
# le domaine est de longuer L0, les deux plans sont entre y=0 et y=H0
# paroi sup
Solid (-y+H0-a*H0*exp(-pow((x)/H0,2)))
#paroi inf
Solid ( y+0.-a*H0*exp(-pow((x)/H0,2)))
# au moins 64 cases, soir dx=0.16
Refine 6
# raffinement pour la paroi...
RefineSolid 9
# on a Poiseuille entre y=0 et y=H0
   Init {} { U = (1 - (y/H0))*(y/H0) }
# adaptation automatique de maillage sur la vorticite
# dx=L0/2^7=10/128=0.08
   AdaptVorticity { istep = 1 } { maxlevel = 7 cmax = 1e-2 }
OutputSimulation { istep = 10 } stdout
OutputTime { istep = 10 } stderr
# sortie des valeurs dans un fichier
OutputSimulation { step= 10 } VAL/knd-%f.gfs
OutputSimulation { start = end } VAL/end.gfs
# sortie sur une ligne au milieu definie par le fichier "Paroi_bas.dat"
  OutputLocation { istep = 10 } VAL/Centre.data Paroi_bas.dat
#test d'arret
EventStop { istart = 20 istep = 1} U 5e-3 DU
}
#condition de vitesse en entree donnee
# et pression nulle en sortie ainsi que dx(U)=0
GfsBox { left = GfsBoundary {
           GfsBcDirichlet U   (1 - (y/H0))*(y/H0)}
         right = GfsBoundary {
           GfsBcDirichlet P 0
           GfsBcNeumann U 0 }
}

#####################################################################

on le sauve dans pois2.gfs
on vérifie que le directory VAL existe tjrs, sinon
mkdir VAL

on lance avec l'option "-m" qui permet de tenir compte des "Define"
gerris2D -m pois2.gfs | gfsview2D vue1.gfv

On met ce ceci dans run.sh

#!/bin/bash
#
#
export LANG=C
echo "lancement du calcul"
rm VAL/knd*

file="Paroi_bas.dat"
awk 'BEGIN{
      L0 = 10
      H0 = 1
      r = 0.1*H0
      wd = L0/2.
      xc = -0.*L0/10.
      print -0.501 * L0   " " 0.0 " " 0.0
      for (x = -wd ; x <= wd; x += 0.01)
        print x " " (0.001 + r*exp(-(x-xc)/H0*(x-xc)/H0)) " 0 ";
      print 0.501 * L0   " " 0.0 " " 0.0
    }' > $file

gerris2D -m pois2.gfs | gfsview2D vue1.gfv

echo "fin normale ?"

gfs2oogl2D -c P -o -i < VAL/end.gfs > VAL/a.dat

On peut mettre ensuite ceci dans le fichier fin.sh
#!/bin/bash
#
#
export LANG=C
echo "liigne de test dans z.txt"
echo "awk '{if (\$1>=tt) print \$0 }' VAL/Centre.data > VAL/Cf.data" > z.txt
echo "extrait le temps final mis dans tfinal.txt"
tfinal=`tail -n 1 VAL/Centre.data`
echo $tfinal > tfinal.txt
awk '{print $1}' tfinal.txt > dump.txt
echo `cat dump.txt`
rm tfinal.txt
echo "substit"
echo "s/tt/"`cat dump.txt`"/g">b.txt
sed `cat b.txt` z.txt>zz.txt
echo `head -n 1 VAL/Cf.data`
echo "transfer des dernieres donnes dans Cf.data"
source zz.txt
rm b.txt
rm z.txt
rm zz.txt

gfs2oogl2D -c P -o -i < `ls VAL/* | tail -n 1 ` > VAL/a.dat

Le fichier run.sh crée un fichier Parois_bas.dat (le bas du canal en fait). Les champs sont sauvés le long des points définis dans ce fichier et mis dans le fichier
VAL/Centre.data
suivant:
# 1:T 2:X 3:Y 4:Z 5:P 6:Pmac 7:U 8:V 9:DU
Ils sont tous sauvés. Mais on ne veut que le temps final. Avec fin.sh on extrait le temps final et on met le dernier paquet dans le fichier Cf.data que l'on trace.
on trace la pression calculée et la pression poiseuille
p'VAL/Cf.data' u 2:5,-(x-5)/750*2
On trace la pression théorique comparée avec FreeFem++
gnuplot> p[-5:5][0:0.035]'VAL/Cf.data'u 2:5,'VAL/a.dat' u ($1*10):2 ,'../../ARTIKLE/fig_nsRNSP_Fig_6/p0_NS_a0p1Re750.txt' u ($1-5):($2) w l,-(x-5)/750*2

On voit que ça colle pas mal

page suivante on augmente la hauteur de la bosse...

Cas suivant