Download villers_platten.dgibi
* fichier : villers_platten.dgibi
***************************************************************
* *
* NOM : villers_platten.dgibi *
* DESCRIPTION : fiche de validation CASTEM2000 *
* Mécanique des Fluides *
* Convection naturelle laminaire + *
* convection thermocapillaire (Marangoni) *
* en cavité rectangulaire (2D) *
* FONCTIONS *
* TESTEES : 'NS' algorithme semi-implicite *
* option SUPG *
* approximation de Boussinesq *
* éléments QUA8 'MACRO' *
* 'TSCA' option SUPG *
* 'TOIMP' pour imposer la contrainte *
* tangentielle de Marangoni en surface *
* 'KOPS' 'GRAD' pour le calcul de cette *
* contrainte (éléments SEG2) *
* *
* AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF) *
* e-mail : gounand@semt2.smts.cea.fr *
* *
***************************************************************
* *
* VERSION : v1, 09/10/97, version initiale *
* HISTORIQUE : v1, 09/10/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 :
*
* ---- Convection naturelle laminaire + convection
* thermocapillaire (Marangoni) sur la surface libre
* (supposée plane).
* Algorithme semi implicite
* opérateur NS (option SUPG)
*
* phi = 0
* - - - - - - - - - - - -
* | |
* T1 = 0 | H | T2 = A = L/H
* | L |
* ------------------------
* phi = 0
*
* U = 0 sur les parois sauf sur la surface ou on applique
* la condition de Levich :
* rho nu (dUx/dz) = (dsigma/dT) (dT/dx).
*
* La structure de l'écoulement ne dépend que des 4 nombres
* adimensionnés suivants :
* - le rapport d'allongement A = L / H ;
* - le nombre de Prandtl : Pr = nu / kappa ;
* - le nombre de Rayleigh :
* Ra = (g beta (T2 - T1) h**3) / (kappa nu)
* = Gr * (A Pr) Gr : nombre de Grashof
* - le nombre de Marangoni :
* Ma = -( dsigma/dT (T2 - T1) h) / (rho alpha nu)
* = Ma2 * (A Pr) Ma2 : autre définition du
* nb de Marangoni
*************************************************************
* Références : - Gounand, S. 1997
* Simulation Numérique de cellules
* thermoconvectives.
* Rapport CEA/DMT 97/357.
* - Villers, D. & Platten, J.K. 1991
* Coupled buoyancy and Marangoni convection
* in acetone : experiments and comparison
* with numerical simulations
* J. Fluid Mech. Vol. 234 pp. 487-510.
*
***************************************************************
*
* JEU DE DONNEE :
*
prg = ('CHAINE' 'Villers & Platten') ;
**** Options
* court = VRAI : le cas-test s'effectue en une quinzaine
* de secondes.
* Les paramètres sont les suivants :
* A = 4 ; Pr = 4 ;
* Ra = 0 ; Ma = 1000
* => une cellule stationnaire
* La comparaison s'effectue sur le profil
* Ux(y) le long de la médiane verticale.
* Il doit se rapprocher de la solution
* obtenue pour une cavité infinie, lorsque
* l'effet Maragoni est la seule cause du mvt :
*
* Ux(z) proportionnel à z(z-2/3)
*
* court = FAUX : maillage raffiné ; on utilise les
* paramètres suivants :
* A = 12 ; Pr = 4.24 ;
* Ra = 4300 ; Ma = 7600
* configuration à une cellule + (1-2)
* rouleaux périphériques.
* La comparaison s'effectue sur le profil
* Ux(y) le long de la médiane verticale.
* Il doit se rapprocher des solutions
* numériques et expérimentales obtenue par
* Villers et Platten (voir références)
* environ 2h 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' ou 'SUPGCC'
*
court = VRAI ;
graph = FAUX ;
inta = FAUX ;
typelem = 'QUA8' ;
discconv = 'SUPG' ;
'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
** -> CALCMAR : calcul de la contrainte de Marangoni
** (cisaillement) sur la surface.
** (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 incluant
** la fonction de courant (qui est calculee)
** (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 Prandtl : ' Pr ) ;
'MESSAGE' ('CHAINE' 'Nombre de Rayleigh : ' ('ENTIER' Ra) ) ;
'MESSAGE' ('CHAINE' 'Nombre de Marangoni : ' ('ENTIER' Ma) ) ;
'MESSAGE' ('CHAINE' 'Nombre de Grashof : ' Gr ) ;
'MESSAGE' ('CHAINE' 'Nombre de Marangon2 : ' Ma2) ;
'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' ;
** CALCMAR
'DEBPROC' CALCMAR ;
'ARGUMENT' trx*'TABLE ' ;
trv = (trx . 'EQEX') ;
tinco = (trv . 'INCO') ;
$surf = (trx . 'DOMZ') ;
clevich = (-1.0D0 '*' (tinco . 'CLEV')) ;
tn = (rv . 'INCO' . 'TN') ;
* On restreint le champ de température à la surface
tn2 = 'KCHT' $surf 'SCAL' 'SOMMET' tn ;
tinco . 'TSURF' = 'KCHT' $surf 'VECT' 'CENTRE'
(('KOPS' tn2 'GRAD' $surf)
'*' clevich) ;
'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') ;
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
th . 'ltps' = ((th . 'ltps') 'ET' ('PROG' (tp . 'TPS'))) ;
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 ;
hlerr = 'ENLEVER' (th . 'lerr') 1 ;
thl = (th . 'lpdt') ;
hlpdt = 'ENLEVER' (thl . 'min') 1 ;
hlconv = 'ENLEVER' (thl . 'conv') 1 ;
hldiff = 'ENLEVER' (thl . 'diff') 1 ;
* Et meme les deux premiers pour l'erreur sur les vitesses
hltp2 = 'ENLEVER' hltps 1 ;
hler2 = 'ENLEVER' hlerr 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 |Erreur Linf Vit.| = fn(t) (convergence)
everr = 'EVOL' 'MANU' 'TEMPS' hltp2 'ERRLINF'
(('LOG' hler2) '/' ('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 ;
** On calcule tout d'abord la fonction de courant
vitesse = (tinco . 'UN') ;
sw = (-1.) * ('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 ;
*-*-*-*-
psi = (rk . 'INCO' . 'PSI') ;
** On effectue une affinité si la cavité est trop allongée
augmen = ((tinco . 'ASPECT') '>' 6.0D0) ;
'SI' augmen ;
dy = 'NOMC' 'UY' (('COORDONNEE' 2 ($dom . 'SOMMET'))
'ET' ('COORDONNEE' 2 ($dom . 'CENTRE')))
'NATU' 'DISCRET' ;
dx = 'NOMC' 'UX' (dy '*' 0.0D0) 'NATU' 'DISCRET' ;
dv = 2.0D0 '*' (dx 'ET' dy) ;
svform = 'FORME' ;
'FORME' dv ;
'FINSI' ;
** Puis on trace les différents champs :
tphy = (trv . 'PASDETPS' . 'TPS') ;
* Vitesses
vref = (tinco . 'UREF') ;
titv = ('CHAINE' 'Champ de vitesse (VECT SOMMET) à t = '
tphy) ;
vv = 'VECTEUR' vitesse (2.0D0 '/' vref) 'UX' 'UY' 'JAUNE' ;
'TRACER' vv cdom 'TITR' titv ;
* Fonction de courant
'OPTION' 'ISOV' 'SULI' ;
titpsi = ('CHAINE' 'Fonction de courant (SCAL SOMMET) à t = '
tphy) ;
'TRACER' psi dom cdom 14 'TITR' titpsi ;
* température
tn = (tinco . 'TN') ;
tittn = ('CHAINE' 'Champ de température (SCAL SOMMET) à t = '
tphy) ;
'TRACER' tn dom cdom 14 'TITR' tittn ;
* 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 ;
** On rétablit le rapport d'aspect initial
'SI' augmen ; 'FORME' svform ; 'FINSI' ;
'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 : paramètres du fluide
* tg, td : températures sur les murs gauches et droits
h = 1.0D0 ;
'SI' court ;
nh = 3 ; aspect = 4.0D0 ;
'SINON' ;
nh = 15 ; aspect = 12.0D0 ;
'FINSI' ;
nl = 'ENTIER' (nh '*' 4) ;
l = h '*' aspect ;
egeom = 1.0D-4 ;
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 ;
pAB = lp 0.0D0 ; pCD = lp h ;
bas = pA 'DROIT' nl pB ; haut = pD 'DROIT' nl pC ;
rwall = pB 'DROIT' nh pC ; lwall = pA 'DROIT' nh pD ;
mediany = pAB 'DROIT' nh pCD ;
mt = 'DALLER' haut rwall bas lwall ;
$mt = 'DOMA' mt 'MACRO' 'IMPR' ;
mt = ($mt . 'MAILLAGE') ; 'TASSER' mt ;
haut = 'CHANGER' haut SEG2 ; bas = 'CHANGER' bas SEG2 ;
lwall = 'CHANGER' lwall SEG2 ; rwall = 'CHANGER' rwall SEG2 ;
mediany = 'CHANGER' mediany SEG2 ;
enleelem = haut 'ELEM' 1 ;
haut2 = 'DIFF' haut enleelem ;
cmt = (lwall 'ET' bas 'ET' rwall 'ET' haut) ;
$haut = 'DOMA' haut 'INCL' $mt egeom ;
$haut2 = 'DOMA' haut2 'INCL' $mt egeom ;
$bas = 'DOMA' bas 'INCL' $mt egeom ;
'ELIMINATION' (mt 'ET' bas 'ET' rwall 'ET' haut
'ET' haut2 'ET' lwall 'ET' cmt 'ET' mediany) egeom ;
**** Initialisation de la table des historiques
*
thist = 'TABLE' 'HIST' ;
thist . 'tpscpu' = ('TEMPS' 'NOEC') ;
* Listes
thist . 'ltps' = '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' 'CLPU' 'CBOU' 'TN' 'TREF'
'INCO' 'UN'
'OPTI' discconv
'ZONE' $mt 'OPER' 'TSCA' 'CLPT' 'UN' 0.0D0
'INCO' 'TN' ;
rv = 'EQEX' rv
'ZONE' $haut 'OPER' 'CALCMAR'
'ZONE' $haut 'OPER' 'TOIMP' 'TSURF'
'INCO' 'UN' ;
tbas = 'KCHT' $bas 'SCAL' 'SOMMET'
(tg '+' ((('COORDONNEE' 1 bas)
'/' l)
'*' (td '-' tg)) ) ;
* Conditions aux limites :
* traitement special pour pA et pB communs respectivement
* à lwall-bas et a bas-rwall sinon les temperatures valent
* deux fois la valeur qu'on veut leur donner
* On utilise 'MANUEL' car 'CLIM' ne prend que les objets de
* type 'MAILLAGE'
*
rv = 'EQEX' rv
'CLIM' 'UN' 'UIMP' (lwall 'ET' rwall 'ET' bas) 0.0D0
'CLIM' 'UN' 'VIMP' (lwall 'ET' rwall 'ET' haut 'ET' bas) 0.0D0
'CLIM' 'TN' 'TIMP' lwall tg
'CLIM' 'TN' 'TIMP' ('MANUEL' 'POI1' pA) (-1.0D0 '*' tg)
'CLIM' 'TN' 'TIMP' rwall td
'CLIM' 'TN' 'TIMP' ('MANUEL' 'POI1' pB) (-1.0D0 '*' td)
'CLIM' 'TN' 'TIMP' bas tbas ;
'SI' graph ;
rv = 'EQEX' rv
'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
*
* Pr : Nombre de Prandtl
* Ra : Nombre de Rayleign (Gr : Grashof)
* Ma : Nombre de Marangoni (Ma2 : autre définition)
* tref : température de référence (approx. de Boussinesq)
* uref : échelle de vitesse (cf. adimensionnement)
* ampvit : facteur d'amplification pour le tracé des vitesses
'SI' court ;
Pr = 4.0D0 ; Ra = 0.0D0 ; Ma = 1000.0D0 ;
'SINON' ;
Pr = 4.24D0 ; Ra = 4300.0D0 ; Ma = 7600.0D0 ;
'FINSI' ;
Gr = Ra '/' (Pr '*' aspect) ; Ma2 = Ma '/' (Pr '*' aspect) ;
sMa2 = (Ma2 '**' 0.5D0) ;
clpu = 1.0D0 '/' sMa2 ;
cbou = Gr '/' Ma2 ;
clpt = 1.0D0 '/' (sMa2 '*' Pr) ;
clev = 1.0D0 ;
uref = sMa2 ;
tref = (tg + td) / 2.0D0 ;
ampvit = 2.0D0 '/' uref ;
*
** 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 = 10 ;
alpha = 0.9D0 ; beta = 1.0D0 ;
'SINON' ;
*nbiter = 10 ; nptrc = 500 ; fidt = 100 ;
nbiter = 10000 ; nptrc = 2000 ; fidt = 100 ;
alpha = 0.5D0 ; beta = 1.0D0 ;
'FINSI' ;
*
** Remplissage des tables
*
rv . 'INCO' . 'ASPECT' = aspect ;
rv . 'INCO' . 'CLPU' = clpu ;
rv . 'INCO' . 'CBOU' = g '*' cbou ;
rv . 'INCO' . 'TREF' = tref ;
rv . 'INCO' . 'UREF' = uref ;
rv . 'INCO' . 'CLPT' = clpt ;
rv . 'INCO' . 'CLEV' = clev ;
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' ;
****
***
**
*
vit = (rv . 'INCO' . 'UN') ;
evvm = 'EVOL' 'BLEU' 'CHPO' vit 'UX' mediany ;
** Le cas-test donne-t-il un résultat correct ??
'OPTION' 'ECHO' 0 ;
'SI' court ;
'OPTION' 'ECHO' 1 ;
*
** Vérification :
** La composante horizontale de la vitesse le long de la
** médiane verticale doit etre prop. à z(z-2/3)
*
cfpoly evpoly = @POMI evvm 2 ;
cfprop = (cfpoly . 2) ;
cflin = (cfpoly . 1) '/' cfprop ;
cfcte = (cfpoly . 0) '/' cfprop ;
cflana = -2.0D0 '/' 3.0D0 ;
* Résultat obtenu le 14/10/97 :
cfpref = -.3128910770680292D+01 ;
critere1 = 'ABS' ((cfprop '-' cfpref) '/' cfpref) ;
critere2 = 'ABS' ((cflin '-' cflana) '/' cflana) ;
critere3 = 'ABS' cfcte ;
* Résultats obtenus le 14/10/97 :
* cflref = -.6784204323293079D+00 ;
* cfcref = .5876882831320513D-02 ;
* correspondant à :
* critere1 = 2.83861D-16
* critere2 = 1.76306D-02
* critere3 = 5.87688D-03
limite1 = 1.D-07 ;
limite2 = 2.D-02 ;
limite3 = 1.D-02 ;
'SI' graph ;
tabg = 'TABLE' ; tabg . 'TITRE' = 'TABLE' ;
tabg . 1 = 'NOLI MARQ CROI' ;
tabg . 'TITRE' . 1 = 'CHAINE' 'Résultat Castem2000' ;
tabg . 'TITRE' . 2 = 'CHAINE' 'Parabole théorique' ;
evpoly = ('COULEUR' evpoly 'JAUN') ;
titg = 'CHAINE' 'Ux(z) sur la médiane verticale' ;
'DESSIN' (evvm 'ET' evpoly) 'MIMA' 'TITR' titg
'LEGE' tabg ;
'FINSI' ;
* Tps CPU (le 10/10/97) ; 1797 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' '------------' ;
'MESSAGE' 'Comparaison de ux(z) sur la médiane verticale à : ' ;
'MESSAGE' 'alpha (z**2 - 2/3 z + 0)' ;
'SAUTER' 'LIGNE' ;
'MESSAGE' ('CHAINE'
'Erreur/Limite sur le coefficient alpha :') ;
'MESSAGE' ('CHAINE' 'FORMAT' '(D24.16)' critere1 ' / ' limite1) ;
'MESSAGE' ('CHAINE'
'Erreur/Limite sur le coefficient de z (2/3)') ;
'MESSAGE' ('CHAINE' 'FORMAT' '(D24.16)' critere2 ' / ' limite2) ;
'MESSAGE' ('CHAINE'
'Erreur/Limite sur le coefficient constant (0)') ;
'MESSAGE' ('CHAINE' 'FORMAT' '(D24.16)' critere3 ' / ' limite3) ;
'MESSAGE' ('CHAINE' 'Temps CPU / Temps CPU (10/10/97) : ') ;
'MESSAGE' ('CHAINE' ('ENTIER' (thist . 'tpscpu')) ' / ~1800'
' centième(s) de secondes') ;
'SI' ('OU' ('OU' (critere1 '>' limite1)
(critere2 '>' limite2))
(critere3 '>' limite3)) ;
'ERREUR' 5 ;
'FINSI' ;
'MESSAGE' 'Tout s_est bien passé' ;
'SINON' ;
'OPTION' 'ECHO' 1 ;
**
** Vérification pour le cas long :
** on compare ux(z) sur la médiane avec des résultats
** experimentaux (attention, l'expérience n'est pas 2D =>
** vitesse plus faible a cause des parois latérales)
** et une autre simulation numérique (effectuée avec les
** memes parametres que la notre)..
** Voir l'article de Villers et Platten pour plus de précisions
**
zexp = 'PROG' 0.0D0
0.1D0 0.2D0 0.35D0 0.57D0 0.75D0 0.97D0
1.2D0 1.5D0 1.8D0 2.0D0 2.2D0 2.5D0 ;
uxexp = 'PROG' 0.0D0
0.6D0 0.75D0 0.95D0 1.05D0 0.95D0 0.8D0
0.4D0 -0.4D0 -1.0D0 -2.0D0 -2.5D0 -3.6D0 ;
zpas = (2.5D0 '/' 32) ;
zsvp = 'PROG' 0.0D0 'PAS' zpas 'NPAS' 32 ;
uxsvp = 'PROG' 0.0D0 0.22D0 0.4D0 0.6D0 0.75D0 0.86D0
0.98D0 1.05D0 1.13D0 1.15D0 1.16D0 1.15D0
1.13D0 1.08D0 1.02D0 0.92D0 0.81D0 0.7D0
0.55D0 0.4D0 0.2D0 0.01D0 -0.2D0 -0.43D0
-0.7D0 -0.95D0 -1.2D0 -1.5D0 -1.8D0 -2.15D0
-2.5D0 -2.85D0 -3.25D0 ;
*
* On remet l'évolution obtenu avec C2000 en dimensionné
*
* mu à 20°C densité
nuacet = 0.3265D-3 '/' 870.0D0 ;
lref = 2.5D0 ;
uref = ((nuacet '*' 1.0D6) '/' lref) '*' sMa2 ;
zdim = ('EXTRAIRE' evvm 'ABSC') '*' lref ;
uxdim = ('EXTRAIRE' evvm 'ORDO') '*' uref ;
* Comparaison des profils :
* on les remet sur zdim par interpolation linéaire
uxexpi = 'IPOL' zdim zexp uxexp ;
uxvpi = 'IPOL' zdim zsvp uxsvp ;
lmux = 'PROG' ('MAXIMUM' uxdim 'ABS')
('MAXIMUM' uxexpi 'ABS')
('MAXIMUM' uxvpi 'ABS') ;
echux = 'MINIMUM' lmux ;
critere1 = ('MAXIMUM' (uxdim '-' uxvpi) 'ABS') '/' echux ;
critere2 = ('MAXIMUM' (uxdim '-' uxexpi) 'ABS') '/' echux ;
* Résultats obtenus le 15/10/97 :
* critere1 = .7530548488495176D-01
* critere2 = .3228687516012702D+00
limite1 = 0.08D0 ;
limite2 = 0.35D0 ;
* trace eventuel
'SI' graph ;
evexp = 'EVOL' 'ROSE' 'MANU' 'Z' zexp 'UX' uxexp ;
evsvp = 'EVOL' 'JAUN' 'MANU' 'Z' zsvp 'UX' uxsvp ;
evexpi = 'EVOL' 'ROUG' 'MANU' 'Z' zdim 'UX' uxexpi ;
evvpi = 'EVOL' 'VERT' 'MANU' 'Z' zdim 'UX' uxvpi ;
evdim = 'EVOL' 'TURQ' 'MANU' 'Z' zdim 'UX' uxdim ;
titd = 'CHAINE' 'Profil Ux(z) sur la médiane verticale' ;
tabd = 'TABLE' ; tabd . 'TITRE' = 'TABLE' ;
tabd . 1 = 'MARQ PLUS TIRR TITR Experience' ;
tabd . 2 = 'MARQ CROI TIRC TITR Vp_num' ;
tabd . 3 = 'MARQ PLUS NOLI TITR Exp_ipol' ;
tabd . 4 = 'MARQ CROI NOLI TITR Vp_num_ipol' ;
tabd . 5 = 'MARQ ETOI NOLI TITR K2000' ;
tabd . 'TITRE' . 1 = 'CHAINE' 'Expérience' ;
tabd . 'TITRE' . 2 = 'CHAINE' 'Numérique (V&P)' ;
tabd . 'TITRE' . 3 = 'CHAINE' 'Exp. (interpolé)' ;
tabd . 'TITRE' . 4 = 'CHAINE' 'V&P (interpolé)' ;
tabd . 'TITRE' . 5 = 'CHAINE' 'Castem 2000' ;
'DESSIN' (evexp 'ET' evsvp 'ET' evexpi 'ET' evvpi 'ET' evdim)
'LEGE' tabd 'TITR' titd ;
'FINSI' ;
* Tps CPU (le 10/10/97) ; 7043 secondes
* (court = FAUX ; graph = VRAI ; inta = FAUX ;)
thist . 'tpscpu' = (thist . 'tpscpu') '+' ('TEMPS' 'NOEC') ;
*
** Affichage des messages et des erreurs éventuelles
*
'OPTION' 'ECHO' 0 ;
'SAUTER' 2 'LIGNE' ;
'MESSAGE' 'Debriefing :' ;
'MESSAGE' '------------' ;
'MESSAGE' 'Comparaison de ux(z) sur la médiane verticale' ;
'SAUTER' 'LIGNE' ;
'MESSAGE' 'par rapport à la simulation numérique de V&P' ;
'MESSAGE' 'Erreur / Limite :' ;
'MESSAGE' ('CHAINE' 'FORMAT' '(D24.16)' critere1 ' / ' limite1) ;
'SAUTER' 'LIGNE' ;
'MESSAGE' 'par rapport à l_expérience de V&P' ;
'MESSAGE' 'Erreur / Limite :' ;
'MESSAGE' ('CHAINE' 'FORMAT' '(D24.16)' critere2 ' / ' limite2) ;
'SAUTER' 'LIGNE' ;
'MESSAGE'
('CHAINE' 'Temps CPU / Temps CPU (15/10/97) (moy2 queue) : ') ;
'MESSAGE' ('CHAINE' ('ENTIER' ((thist . 'tpscpu')
'/' 100.0D0)) ' / ~7040'
' secondes') ;
'SI' ('OU' (critere1 '>' limite1)
(critere2 '>' limite2)) ;
'ERREUR' 5 ;
'FINSI' ;
'MESSAGE' 'Tout s_est bien passé' ;
'FINSI' ;
SUMMARY ;
'OPTION' 'ECHO' 1 ;
'SI' inta ; 'OPTION' 'DONN' 5 ; 'FINSI' ;
'FIN' ;
*
** Fin du 'main'
**************************************************************
© Cast3M 2003 - All rights reserved.
Disclaimer