* fichier : benchmark_imst.dgibi
***************************************************************
* *
* NOM : benchmark_imst.dgibi *
* DESCRIPTION : fiche de validation CASTEM2000 *
* Mécanique des Fluides *
* Convection naturelle laminaire *
* en cavité rectangulaire (2D) *
* FONCTIONS *
* TESTEES : 'NS' algorithme semi-implicite *
* option CENTREE *
* approximation de Boussinesq *
* éléments QUA8 'MACRO' *
* *
* AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF) *
* e-mail : gounand@semt2.smts.cea.fr *
* *
***************************************************************
* *
* VERSION : v1, 19/9/97, version initiale *
* HISTORIQUE : v1, 19/9/97, création *
* HISTORIQUE : *
* HISTORIQUE : *
* HISTORIQUE : *
* *
***************************************************************
* *
* Priere de completer la partie HISTORIQUE apres modification *
* du jeu de donnees afin de faciliter la maintenance *
* *
***************************************************************
*
* DESCRIPTION DETAILLEE :
*
* -- IMST -- Convection naturelle laminaire à Prandtl (Pr) = 0
* Algorithme semi implicite
* opérateur NS (option CENTREE)
*
* phi = 0
* ------------------------
* | |
* T1 = 0 | H | T2 = A
* | L |
* ------------------------
* phi = 0
*
* A = L/H = 4 l'allongement
* U = 0 sur les parois
* Le profil de température est imposé T = T1 + (T2-T1)*(x/L)
*
* La structure de l'écoulement ne dépend que
* du nombre de Grashof Gr :
* Gr = g beta (T2 - T1) H**3 A / nu**2
*
* ^
* | Psi max (valeur max de la fonction de courant)
* |
* |
* -| o o
* | unsteady -> o v
* | (3 cells) . v
* | ----------- o v . o
* | steady -> .^ v o
* | (3 cells) . ^ v .
* | ----------- o ^ o
* -| ^ .
* | steady -> . ^ . <-- steady (2 cells)
* | o
* | (1 cell) .
* |
* | .
* | o
* -|-------------------------------------------> Gr
* I I I
* 1.E4 2.E4 3.E4
**************************************************************
* Références : - Gounand, S. 1997
* Simulation numérique de cellules
* thermoconvectives.
* Rapport CEA/DMT 97/357.
* - Le Garrec, S. 1988
* Convection naturelle en cavité pour des
* fluides à faible nombre de Prandtl.
* Rapport CEA/DMT 88/216.
* - Le Garrec, S. & Magnaud, J.-P.
* Numerical Simulation of oscillatory
* convection in low Prandtl Fluids with
* TRIO-EF.
* A GAMM Workshop Edited by Bernard Roux
* Notes on Numerical Fluid Mechanics,Volume 27
* pp. 189-199. (Vieweg, Braunschweig 1990)
*
***************************************************************
*
* JEU DE DONNEE :
*
prg = ('CHAINE' 'Benchmark Imst') ;
**** Options
* court = VRAI : le cas-test s'effectue en une dizaine
* de secondes.
* Gr = 5000 => une cellule stationnaire
* la comparaison s'effectue sur le profil
* Uy(x) le long de la médiane horizontale.
* On vérifie également que l'écoulement
* est centrosymétrique.
* court = FAUX : maillage raffiné ;
* Gr = 40000 => bifurcation 3->2 cellules
* environ 1h30 CPU avec les graphiques
*
* graph = VRAI : on trace les graphiques (historiques, champs)
* On a plus de calculs avec cette option
* (fonction de courant, convergence...)
* On mettra donc graph = faux pour la fiche
* de validation
*
* inta = VRAI si on est en interactif ; 0 sinon
*
* discconv : option de discrétisation des termes de
* convection pour les équations de Navier-Stokes
* choix possibles : 'CENTREE', 'SUPG'...
* voir notice de 'NS'
court = VRAI ;
graph = FAUX ;
inta = FAUX ;
typelem = 'QUA8' ;
discconv = 'CENTREE' ;
'OPTION' 'DIME' 2 'ELEM' typelem ;
'SI' ('NON' inta) ;
'OPTION' 'TRAC' 'PS' ;
'OPTION' 'ECHO' 0 ;
'SINON' ;
**od = 'TEXTE' 'OPTION DONN 3' ;
'FINSI' ;
*************************************************************
**** Définition de quelques procédures :
** -> SUMMARY : affiche un résumé des paramètres
** -> CALCPSI : calcul de la fonction de courant
** (appelée par 'EXEC' via la table 'EQEX')
** -> MAJHIST : mise à jour d'historiques
** (appelée par 'EXEC' via la table 'EQEX')
** -> TRHIST : tracé des historiques
** -> TRCHAMP : tracé des différents champs
** (appelée par 'EXEC' via la table 'EQEX')
** NB : pour simplifier, on ne transmet pas les paramètres à
** SUMMARY. C'est de la mauvaise programmation...
** Les autres procédures sont correctes
** sous toutes réserves
*
** SUMMARY
'DEBPROC' SUMMARY ;
'MESSAGE' ('CHAINE' 'Programme : ' prg) ;
'SAUTER' 'LIGNE' ;
'MESSAGE' ('CHAINE' '*** GEOMETRIE :') ;
'MESSAGE' ('CHAINE' 'Boite carrée : ' l ' x' h) ;
'MESSAGE' ('CHAINE' ' soient : ' nl 'x' nh '='
(nl '*' nh) ' éléments ' typelem) ;
'MESSAGE' ('CHAINE' ' soient : ' (nl '*' nh '*' 4)
' éléments linéaires') ;
'SAUTER' 'LIGNE' ;
'MESSAGE' ('CHAINE' '*** PHYSIQUE :') ;
'MESSAGE' ('CHAINE' 'Rapport d_aspect : ' ('ENTIER' aspect)) ;
'MESSAGE' ('CHAINE' 'Nombre de Grashof : '('ENTIER' Gr) ) ;
'MESSAGE' ('CHAINE' 'Température de référence : ' tref
' comprise entre 0 et ' ('ENTIER' aspect)) ;
'SAUTER' 'LIGNE' ;
'MESSAGE' ('CHAINE' '*** NUMERIQUE :') ;
'MESSAGE' ('CHAINE' 'Discrétisation ' discconv ' des termes de
convection.');
'MESSAGE' ('CHAINE' ' Alfa = ' alpha) ;
'MESSAGE' ('CHAINE' ' Béta = ' beta) ;
'SAUTER' 'LIGNE' ;
rvpdt = rv . 'PASDETPS' ;
'MESSAGE' ('CHAINE' '*** TEMPS :') ;
'MESSAGE' ('CHAINE' 'Nb. itérations : ' ((rvpdt . 'NUPASDT')
'-' 1)) ;
'MESSAGE' ('CHAINE' 'Temps physique : ' (rvpdt . 'TPS')) ;
'MESSAGE' ('CHAINE' 'Temps CPU : '
('ENTIER' (thist . 'tpscpu')) ' centième(s) de seconde') ;
'FINPROC' ;
** CALCPSI
'DEBPROC' CALCPSI ;
'ARGUMENT' trx*'TABLE ' ;
trv = (trx . 'EQEX') ;
tinco = (trv . 'INCO') ;
$dom = (trx . 'DOMZ') ;
cdom = 'CONTOUR' ($dom . 'MAILLAGE') ;
vitesse = 'KCHT' $dom 'VECT' 'SOMMET'
(tinco . 'UN') ;
sw = 'KOPS' vitesse 'ROT' $dom ;
rk = 'EQEX' $dom 'OPTI' 'EF' 'IMPL'
'ZONE' $dom 'OPER' 'LAPN' 1.0D0 'INCO' 'PSI'
'ZONE' $dom 'OPER' 'FIMP' sw 'INCO' 'PSI'
'CLIM' 'PSI' 'TIMP' cdom 0.0D0 ;
rk . 'INCO' = 'TABLE' 'INCO' ;
rk . 'INCO' . 'PSI' = 'KCHT' $dom 'SCAL' 'SOMMET' 0.0D0 ;
*-*-*-*- A remplacer par EXEC qd ca marchera...
EXIC rk ;
*-*-*-*-
tinco . 'PSI' = (rk . 'INCO' . 'PSI') ;
'FINPROC' trx ;
** MAJHIST
'DEBPROC' MAJHIST ;
'ARGUMENT' trx*'TABLE ' ;
trv = (trx . 'EQEX') ;
tinco = (trv . 'INCO') ;
th = (trv . 'HISTOIRE') ;
tp = (trv . 'PASDETPS') ;
vit = (tinco . 'UN') ; vitm1 = (tinco . 'UN-1') ;
vitref = (tinco . 'UREF') ;
psi = (tinco . 'PSI') ;
npdt = (tp . 'NUPASDT') ;
* Calcul de l'erreur relative (norme Linfini) entre deux
* champs de vitesse successifs
errrel = ('MAXIMUM' (vit '-' vitm1) 'ABS') '/' vitref ;
'SI' ('MULT' npdt (trv . 'FIDT')) ;
'MESSAGE' ('CHAINE' 'Pdt : ' npdt
' |Un - Un_1| Linf = ' errrel) ;
'FINSI' ;
tinco . 'UN-1' = 'COPIER' vit ;
* Mise à jour des différentes listes
mpsi = 'MAXIMUM' psi 'ABS' ;
th . 'ltps' = ((th . 'ltps') 'ET' ('PROG' (tp . 'TPS'))) ;
th . 'lpsi' = ((th . 'lpsi') 'ET' ('PROG' mpsi)) ;
th . 'lerr' = ((th . 'lerr') 'ET' ('PROG' errrel)) ;
thl = (th . 'lpdt') ;
thl . 'conv' = ((thl . 'conv') 'ET' ('PROG' (tp . 'DTCONV'))) ;
thl . 'diff' = ((thl . 'diff') 'ET' ('PROG' (tp . 'DTDIFU'))) ;
thl . 'min' = ((thl . 'min') 'ET' ('PROG' (tp . 'DELTAT-1'))) ;
'FINPROC' trx ;
** TRHIST
'DEBPROC' TRHIST ;
'ARGUMENT' th*'TABLE ' ;
* On enlève le premier indice des listes
hltps = 'ENLEVER' (th . 'ltps') 1 ;
hlpsi = 'ENLEVER' (th . 'lpsi') 1 ;
hlerr = 'ENLEVER' (th . 'lerr') 1 ;
thl = (th . 'lpdt') ;
hlpdt = 'ENLEVER' (thl . 'min') 1 ;
hlconv = 'ENLEVER' (thl . 'conv') 1 ;
hldiff = 'ENLEVER' (thl . 'diff') 1 ;
* Pas de temps = fn(t)
evpdt = 'EVOL' 'MANU' 'TEMPS' hltps 'PDT' hlpdt ;
'DESSIN' evpdt 'MIMA' 'TITR' 'Pas de temps' ;
* LOG10(pas de temps) = fn(t)
ltco = ('LOG' (hlconv '*' alpha)) '/' ('LOG' 10.0D0) ;
ltdi = ('LOG' (hldiff '*' alpha)) '/' ('LOG' 10.0D0) ;
ltmi = ('LOG' hlpdt) '/' ('LOG' 10.0D0) ;
evco = 'COULEUR'
('EVOL' 'MANU' 'TEMPS' hltps 'dtconv' ltco) 'BLEU' ;
evdi = 'COULEUR'
('EVOL' 'MANU' 'TEMPS' hltps 'dtdiff' ltdi) 'ROUG' ;
evmi = 'COULEUR'
('EVOL' 'MANU' 'TEMPS' hltps 'dt' ltmi) 'BLAN' ;
tabt = 'TABLE' ; tabt . 'TITRE' = 'TABLE' ;
tabt . 1 = 'TIRC' ;
tabt . 2 = 'TIRL' ;
tabt . 3 = 'NOLI' ;
tabt . 'TITRE' . 1 = 'CHAINE' 'Pdt convection' ;
tabt . 'TITRE' . 2 = 'CHAINE' 'Pdt diffusion' ;
tabt . 'TITRE' . 3 = 'CHAINE' 'Pdt min(conv, diff)' ;
'DESSIN' (evco 'ET' evdi 'ET' evmi)
'MIMA' 'TITR' 'Pas de temps'
'LEGE' tabt ;
* Max |Fn de courant| = fn(t)
evpsi = 'EVOL' 'MANU' 'TEMPS' hltps 'PSI' hlpsi ;
'DESSIN' evpsi 'MIMA' 'TITR' 'Max |Fn de courant|' ;
* Max |Erreur Linf Vit.| = fn(t) (convergence)
everr = 'EVOL' 'MANU' 'TEMPS' hltps 'ERRLINF'
(('LOG' hlerr) '/' ('LOG' 10.0D0)) ;
'DESSIN' everr 'MIMA' 'TITR' 'Log Erreur Linf Vit.' ;
'FINPROC' ;
** TRCHAMP
'DEBPROC' TRCHAMP ;
'ARGUMENT' trx*'TABLE ' ;
trv = (trx . 'EQEX') ;
'SI' ('NON' ('MULT' (trv . 'PASDETPS' . 'NUPASDT')
(trv . 'NPTRC'))) ;
'QUITTER' TRCHAMP ;
'FINSI' ;
tinco = (trv . 'INCO') ;
$dom = (trx . 'DOMZ') ;
dom = ($dom . 'MAILLAGE') ;
cdom = 'CONTOUR' dom ;
tphy = (trv . 'PASDETPS' . 'TPS') ;
* Vitesses
vitesse = (tinco . 'UN') ; vref = (tinco . 'UREF') ;
titv = ('CHAINE' 'Champ de vitesse (VECT SOMMET) à t = '
tphy) ;
vv = 'VECTEUR' vitesse (0.2D0 '/' vref) 'UX' 'UY' 'JAUNE' ;
'TRACER' vv cdom 'TITR' titv ;
* Fonction de courant
'OPTION' 'ISOV' 'SULI' ;
psi = (tinco . 'PSI') ;
titpsi = ('CHAINE' 'Fonction de courant (SCAL SOMMET) à t = '
tphy) ;
'TRACER' psi dom cdom 14 'TITR' titpsi ;
* Pression
* NB : On n'utilise pas 'ELNO' : on dessine une pression
* constante sur un élément.
pression = (trv . 'PRESSION' . 'PRESSION') ;
modelp = 'MODELISER' dom 'THERMIQUE' ;
elemp = 'KCHA' $dom 'CHAM' pression ;
titpres = ('CHAINE' 'Pression (SCAL CENTRE) à t = '
tphy) ;
'TRACER' elemp modelp cdom 'TITR' titpres ;
'FINPROC' ;
*
** Fin de la définition des procédures
**************************************************************
**************************************************************
** Début du 'main'
*
**** Définition des constantes
* l, h, nl, nh : données pour la boite.
* egeom : erreur géométrique pour 'ELIMINATION'
* aspect : rapport d'aspect de la boite
* nu, g, rho : paramètres du fluide
* tg, td : températures sur les murs gauches et droits
h = 1.0D0 ;
'SI' court ; nh = 3 ; 'SINON' ; nh = 5 ; 'FINSI' ;
aspect = 4.0D0 ;
l = h '*' aspect ;
'SI' court ;
nl = 'ENTIER' (nh '*' 3) ;
'SINON' ;
nl = 'ENTIER' (nh '*' aspect) ;
'FINSI' ;
egeom = 1.0D-4 ;
nu = 1.0D0 ; g = (0.D0 -1.D0) ;
tg = 0.D0 ; td = aspect ;
**** Création des maillages
lp = l '/' 2.0D0 ; hp = h '/' 2.0D0 ;
mlp = (-1.0D0) '*' lp ; mhp = (-1.0D0 '*' hp) ;
pmil = lp hp ;
pA = 0.0D0 0.0D0 ; pB = l 0.0D0 ;
pC = l h ; pD = 0.0D0 h ;
pDA = 0.0D0 hp ; pBC = l hp ;
bas = pA 'DROIT' nl pB ; haut = pD 'DROIT' nl pC ;
rwall = pB 'DROIT' nh pC ; lwall = pA 'DROIT' nh pD ;
median = pDA 'DROIT' nl pBC ;
mt = 'DALLER' haut rwall bas lwall ;
mt = chan mt QUAF ;
$mt= MODE mt 'NAVIER_STOKES' MACRO ; DOMA $MT 'IMPR' ;
mt = (DOMA $mt 'MAILLAGE') ; 'TASSER' mt ;
haut = 'CHANGER' haut SEG2 ; bas = 'CHANGER' bas SEG2 ;
lwall = 'CHANGER' lwall SEG2 ; rwall = 'CHANGER' rwall SEG2 ;
median = 'CHANGER' median SEG2 ;
cmt = (lwall 'ET' bas 'ET' rwall 'ET' haut) ;
'ELIMINATION' (mt 'ET' bas 'ET' rwall 'ET' haut
'ET' lwall 'ET' cmt 'ET' median) egeom ;
**** Initialisation de la table des historiques
*
thist = 'TABLE' 'HIST' ;
thist . 'tpscpu' = ('TEMPS' 'NOEC') ;
* Listes
thist . 'ltps' = 'PROG' ;
thist . 'lpsi' = 'PROG' ;
thist . 'lerr' = 'PROG' ;
thist . 'lpdt' = 'TABLE' 'LPDT' ; thl = (thist . 'lpdt') ;
thl . 'conv' = 'PROG' ; thl . 'diff' = 'PROG' ;
thl . 'min' = 'PROG' ;
**** Initialisation des tables rv et rvp
rv = 'EQEX' $mt 'OPTI' discconv
'ZONE' $mt 'OPER' 'NS' 'NU' 'GB' 'TN' 'TREF'
'INCO' 'UN'
'OPTI' 'CENTREE'
'ZONE' $mt 'OPER' 'DFDT' 1. 'UN' 'DELTAT' 'INCO' 'UN'
;
* Conditions aux limites
rv = 'EQEX' rv
'CLIM' 'UN' 'UIMP' cmt 0.0D0
'UN' 'VIMP' cmt 0.0D0 ;
'SI' graph ;
rv = 'EQEX' rv
'ZONE' $mt 'OPER' 'CALCPSI'
'ZONE' $mt 'OPER' 'MAJHIST'
'ZONE' $mt 'OPER' 'TRCHAMP' ;
'FINSI' ;
rvp = 'EQPR' $mt
'ZONE' $mt 'OPER' 'PRESSION' 0.D0
'PIMP' 0.D0 ;
rv . 'INCO' = 'TABLE' 'INCO' ;
rv . 'PRESSION' = rvp ;
rv . 'HISTOIRE' = thist ;
**** Initialisation des inconnues
* V -> 0 ; T gradient constant avec tg sur lwall, td sur rwall
rv . 'INCO' . 'UN-1' = 'KCHT' $mt 'VECT' 'SOMMET' (0.D0 0.D0) ;
rv . 'INCO' . 'UN' = 'KCHT' $mt 'VECT' 'SOMMET' (0.D0 0.D0) ;
rv . 'INCO' . 'TN' = 'KCHT' $mt 'SCAL' 'SOMMET'
(tg '+' ((('COORDONNEE' 1 mt)
'/' l)
'*' (td '-' tg)) ) ;
**** Valeurs des paramètres :
*
** Physiques
*
* Gr : Nombre de Grashof
* tref : température de référence (approx. de Boussinesq)
* Il faut la dissymétriser pour observer la transition
* 3 -> 2 cellules
* uref : échelle de vitesse (cf. adimensionnement)
* ampvit : facteur d'amplification pour le tracé des vitesses
*
'SI' court ;
Gr = 5.0D3 ; tref = (tg '+' td) '/' 2.0D0 ;
'SINON' ;
Gr = 40.0D3 ; tref = tg ;
'FINSI' ;
uref = (Gr '**' 0.5D0) ;
*
** Numériques
*
* nbiter : nombre d'itérations
* nptrc : fréquence pour le tracé des champs (.ps ou à l'écran)
* fidt : fréquence d'affichage des messages à l'écran
* alpha : multiplicateur du pas de temps
* beta : paramètre de stabilisation de la pression
* pour les éléments MACRO
'SI' court ;
nbiter = 250 ; nptrc = 125 ; fidt = 50 ;
'SINON' ;
nbiter = 15000 ; nptrc = 2500 ; fidt = 100 ;
'FINSI' ;
alpha = 0.9D0 ; beta = 1.0D0 ;
*
** Remplissage des tables
*
rv . 'INCO' . 'NU' = nu ;
rv . 'INCO' . 'GB' = g '*' Gr ;
rv . 'INCO' . 'TREF' = tref ;
rv . 'INCO' . 'UREF' = uref ;
rv . 'ALFA' = alpha ;
rv . 'PRESSION' . 'BETA' = beta ;
rv . 'FIDT' = fidt ;
rv . 'NPTRC' = nptrc ;
rv . 'ITMA' = nbiter ;
*
**
***
**** Exécution
SUMMARY ;
'SI' inta ; 'OPTION' 'DONN' 5 ; 'FINSI' ;
EXEC rv ;
thist . 'tpscpu' = (thist . 'tpscpu') '+' ('TEMPS' 'NOEC') ;
'SI' graph ; TRHIST thist ; 'FINSI' ;
****
***
**
*
** Le cas-test donne-t-il un résultat correct ??
'SI' court ;
*
** Vérification de la centrosymétrie sur les vitesse
*
vit = (rv . 'INCO' . 'UN') ;
tvit = vit 'TOURNER' 180.0D0 pmil ;
vit = 'KCHT' $mt 'VECT' 'SOMMET' vit ;
'ELIMINATION' (('EXTRAIRE' tvit 'MAIL') 'ET' mt) egeom ;
tvit = 'KCHT' $mt 'VECT' 'SOMMET' tvit ;
dv = 'KOPS' vit '-' tvit ;
critere1 = ('MAXIMUM' dv 'ABS') '/' uref ;
* Valeur de critere1 au 2/10/97 : 1.718764099177966D-15
limite1 = 1.0D-13 ;
*
** Evolution vitesse verticale le long de la médiane
*
evvm = 'EVOL' 'CHPO' vit 'UY' median ;
evym = 'EXTRAIRE' evvm 'ORDO' ;
*
* Résultat obtenu le 2/10/97
* correspondant à critere2 = 1.507288760336422D-16
evvref = 'PROG' 0.0D0 -.2205372973655706D+02
-.1537657415971809D+02 -.9886959157032784D+01 -.4994589224716669D+01
-.1325630316341117D+01 .5065207872228873D+00 .9297558065189220D+00
.5908058903982639D+00 .2945842774751989D-13 -.5908058903982122D+00
-.9297558065189337D+00 -.5065207872229142D+00 .1325630316341034D+01
.4994589224716784D+01 .9886959157032821D+01 .1537657415971808D+02
.2205372973655705D+02 0.0D0 ;
critere2 = ('MAXIMUM' (evym '-' evvref) 'ABS') '/' uref ;
limite2 = 1.0D-10 ;
* Tps CPU (le 2/10/97) ; 1009 centièmes de secondes
* (court = VRAI ; graph = FAUX ; inta = FAUX ;)
thist . 'tpscpu' = (thist . 'tpscpu') '+' ('TEMPS' 'NOEC') ;
*
** Affichage des messages et des erreurs éventuelles
*
'OPTION' 'ECHO' 0 ;
'SAUTER' 2 'LIGNE' ;
'MESSAGE' 'Debriefing :' ;
'MESSAGE' '------------' ;
'SAUTER' 'LIGNE' ;
'MESSAGE' ('CHAINE'
'Erreur/Limite de centrosymétrie du champ de vitesses :') ;
'MESSAGE' ('CHAINE' 'FORMAT' '(D24.16)' critere1 ' / ' limite1) ;
'MESSAGE' ('CHAINE'
'Erreur/Limite Vy sur la médiane horizontale') ;
'MESSAGE' ('CHAINE' 'FORMAT' '(D24.16)' critere2 ' / ' limite2) ;
'MESSAGE' ('CHAINE' 'Temps CPU / Temps CPU (30/9/97) : ') ;
'MESSAGE' ('CHAINE' ('ENTIER' (thist . 'tpscpu')) ' / ~1000'
' centième(s) de secondes') ;
'SI' ('OU' (critere1 '>' limite1) (critere2 '>' limite2)) ;
'ERREUR' 5 ;
'FINSI' ;
'MESSAGE' 'Tout s_est bien passé' ;
'OPTION' 'ECHO' 1 ;
'FINSI' ;
SUMMARY ;
'SI' inta ; 'OPTION' 'DONN' 5 ; 'FINSI' ;
'FIN' ;
*
** Fin du 'main'
**************************************************************
© Cast3M 2003 - All rights reserved.
Disclaimer