***************** CAS TEST : TP4.DGIBI *************************
*
* Convection naturelle dans un cylindre uniformement chauffé
* (Navier-Stokes incompressible et approximation de Boussinesq)
*
* Ref biblio : Gartling D.K., A finite element analysis of
* volumetrically heated fluids in an axisymmetric enclosure,
* in Finite Elements in Fluids, vol4, pp233-250, 1982.
*
*------------------------------------------------------------------
*
* Algorithme de projection et élément (v,T)/p MACRO/CENTREP1
*
*------------------------------------------------------------------
*
* Auteur : Frédéric DABBENE
* ttmf3@semt2.smts.cea.fr
*
*------------------------------------------------------------------
* Procédures spécifiques à évoluer et à généraliser
*------------------------------------------------------------------
* COURANT : Calcul de la fonction de courant d'un domaine fermé
* CALCNUSS : Calcul du Nusselt (en adimensionné = grad THETA)
* RESIDU : Calcul du residu en température
*------------------------------------------------------------------
*
*
* COMPLET : Booleen mis à FAUX pour les tests de non régresssion
* GRAPH : Booleen réalisation ou non des post-traitements
* POST1 : Booleen indiquant si on affiche le résidu en NCLK
* POST2 : Booleen indiquant si on affiche la température en NCLK
* IRESU : Critère d'arret sur l'incrément de température
*
COMPLET = FAUX ;
GRAPH = FAUX ;
POST1 = FAUX ;
POST2 = FAUX ;
'SI' COMPLET ;
IRESU = 1.D-6 ;
'SINO' ;
IRESU = 1.D-3 ;
'FINSI' ;
*
*
*
*
*==================================================================
* Calcul de la fonction de courant d'un domaine fermé
*==================================================================
* E/ : UN : CHPO : Champ de vitesse
* E/ : $DOMAINE : MMODEL 'NAVIER_STOKES' volumique
* E/ : $CONTOUR : MMODEL 'NAVIER_STOKES' surfacique
* /S : PSI : CHPO : Fonction de courant
*------------------------------------------------------------------
* On vérifie qu'on est en dimension 2 mais pas que div(UN)=0
* Si le domaine est ouvert, modifier les conditions aux limites
*------------------------------------------------------------------
'DEBPROC' courant un*chpoint $domaine*mmodel $contour*mmodel ;
VAL0 = 'VALE' 'DIME' ;
TEST = 'EGA' VAL0 2 ;
'SI' TEST ;
'MESS' 'Remember that Velocity have to be at zero divergence' ;
'SINON' ;
'ERRE' 832 ;
'QUIT' courant ;
'FINSI' ;
VAL1 = 'VALE' 'MODE' ;
TEST = 'EGA' VAL1 'AXIS' ;
'SI' TEST ;
ROTU0 = 'KOPS' un 'ROT' $domaine ;
XC1 YC1 = 'COOR' ('DOMA' $domaine 'CENTRE') ;
VZ1 = 'NOEL' $domaine ('EXCO' un 'UY' 'SCAL') ;
ROTU0 = 2. * VZ1 '-' (ROTU0 * XC1) ;
'SINO' ;
ROTU0 = 'KOPS' un 'ROT' $domaine ;
'FINSI' ;
CONT0 = 'DOMA' $contour 'MAILLAGE' ;
RK = 'EQEX' $domaine 'OPTI' 'EF' 'IMPL'
'ZONE' $domaine 'OPER' 'LAPN' -1. 'INCO' 'PSI'
'ZONE' $domaine 'OPER' 'FIMP' ROTU0 'INCO' 'PSI'
'CLIM' 'PSI' 'TIMP' CONT0 0. ;
RK . 'INCO' = 'TABLE' 'INCO' ;
RK . 'INCO' . 'PSI' = 'KCHT' $domaine 'SCAL' 'SOMMET' 0. ;
EXEC RK ;
psi = 'COPI' RK . 'INCO' . 'PSI' ;
'FINP' psi ;
*
*==================================================================
* Calcul du Nusselt (en adimensionné = grad THETA)
*==================================================================
* E/ : TN : CHPO : Température
* E/ : $DOMAINE : MMODEL 'NAVIER_STOKES' volumique
* E/ : $COTE : MMODEL 'NAVIER_STOKES' surfacique
* /S : NUSS1 : CHPO : Nusselt local aux sommets
* /S : RES1 : FLOTTANT : Nusselt moyen
*------------------------------------------------------------------
*
'DEBPROC' calcnuss tn*chpoint $domaine*mmodel $cote*mmodel ;
GRADC0 = 'KOPS' tn 'GRAD' $domaine ;
GRADS0 = 'ELNO' $domaine GRADC0 ;
GS0 = 'KCHT' $cote 'VECT' 'SOMMET' GRADS0 ;
GC0 = 'NOEL' $cote GS0 ;
NORM1 = 'DOMA' $domaine 'NORMALE' ;
NCOTE = 'KCHT' $cote 'VECT' 'CENTRE' NORM1 ;
MOT1 = 'MOTS' 'UX' 'UY' ;
NUSS1 = 'PSCA' GC0 NCOTE MOT1 MOT1 ;
S1 = 'DOMA' $cote 'VOLUME' ;
NUSS2 = NUSS1 * S1 ;
RES1 = 'MAXI' ('RESU' NUSS2) ;
'FINP' nuss1 res1 ;
*------------------------------------------------------------------
*
*
*==================================================================
* Calcul du résidu en température et arrêt suivant un critère
*==================================================================
* E/ : RVX : TABLE : TABLE des données créées par EQEX
* ARG1 : Fréquence d'impression
* ARG2 : Critère d'arrêt
* /S : MAT1 : MATRIK : Objet vide
* /S : CHP1 : CHPO : Objet vide
*------------------------------------------------------------------
*
'DEBPROC' residu rvx*table ;
RV = rvx . 'EQEX' ;
FREQ = RVX . 'ARG1' ;
EPS0 = RVX . 'ARG2' ;
NITER = RV . 'NITER' ;
DD = RV . 'PASDETPS' . 'NUPASDT' ;
NN = DD '/' FREQ ;
L0 = 'EGA' (DD '-' (FREQ*NN)) 0 ;
'SI' L0 ;
RANG0 = RV . 'PASDETPS' . 'NUPASDT' ;
TIME0 = RV . 'PASDETPS' . 'TPS' ;
TN0 = RV . 'INCO' . 'TN' ;
TNM0 = RV . 'INCO' . 'TN2' ;
ERR0 = ('MAXIMUM' ('ABS' (TN0 '-' TNM0))) '+' 1.D-20 ;
ERR10 = ('LOG' ERR0 ) '/' ('LOG' 10.) ;
'MESSAGE' 'Résidu en Température au pas'
RANG0 '(' TIME0 ') :' ERR0 ':' ERR10 ;
RV . 'INCO' . 'IT' = RV . 'INCO' . 'IT' 'ET' ('PROG' RANG0) ;
RV . 'INCO' . 'TI' = RV . 'INCO' . 'TI' 'ET' ('PROG' TIME0) ;
RV . 'INCO' . 'ER' = RV . 'INCO' . 'ER' 'ET' ('PROG' ERR10) ;
EV1 = 'EVOL' 'MANUEL' (RV . 'INCO' . 'IT') (RV . 'INCO' . 'ER') ;
Y1 = ('LOG' EPS0) '/' ('LOG' 10) ;
'SI' POST1 ;
X1 = 0. ; X2 = RV . 'ITMA' ;
'DESSIN' EV1 'YBOR' Y1 0. 'NCLK' ;
'FINSI' ;
'SI' POST2 ;
L1 = (PROG 10. PAS 5. 90.)* 1.D-3 ;
trace L1 tn0 DOM1 CNT1 'TITR' 'Température' 'NCLK' ;
'FINSI' ;
'SI' ((ERR10 < Y1) 'ET' (DD > 10)) ;
RV . 'TFINAL' = RV . 'PASDETPS' . 'TPS' ;
'FINSI' ;
'FINSI' ;
RV . 'INCO' . 'TN2' = 'COPIER' RV . 'INCO' . 'TN' ;
mat1 chp1 = 'KOPS' 'MATRIK' ;
'FINP' mat1 chp1 ;
*------------------------------------------------------------------
*
*
'OPTION' 'DIME' 2 'ELEM' 'QUA8' 'MODE' 'AXIS' 'ISOV' 'SULI' ;
*
* NLI0 : Nombre d'isovaleurs
*
NLI0 = 14 ;
*
*------------------------------------------------------------------
* Choix de l'algorithme et des discrétisations en vitesse/pression
*------------------------------------------------------------------
* Les couples vitesse/pression (ICHOI/)
* :( 1/ LINE + CENTRE (déconseillé,--couteux, -stable -précis)
* (P1/P0 et Q1/P0)
* :( 2/ MACRO + CENTRE (déconseillé,--couteux, =stable -précis)
* (P1/P0 stabilisé et Q1/P0 stabilisé)
* 3/ MACRO + CENTREP1 (conseillé, -couteux, +stable, =précis)
* :( (isoP2/isoP1nc à éviter (PAS DE TRIANGLE SVP))
* :) (isoQ2/isoP1nc à consommer sans modération)
* :) 4/ QUAF + CENTREP1 (conseillé, +couteux, +stable, +précis)
* :) 5/ QUAF + MSOMMET (conseillé,++couteux, +stable, +précis)
* Les autres combinaisons ne sont pas dans les notices de
* NAVI, EXEC, EQEX, EQPR et DOMA. Donc ...
*------------------------------------------------------------------
* Les algorithmes (IRESO/)
* :( 1/ avec itérations internes et méthode de projection
* (à utiliser avec précaution : TESTS INSUFFISANTS)
* :) 2/ sans itérations internes et méthode de projection
* :) 3/ full implicite
* :) 4/ semi-explicite et RV . 'PRESSSION' = RVP
* (déconseillé, ancienne syntaxe : nouvelle ci-dessous)
* :( 5/ semi-explicite et RV . 'POISSON' = RVP
* Les autres choix ne sont pas dans les notices de
* NAVI, EXEC, EQEX, EQPR et DOMA. Donc ...
*------------------------------------------------------------------
* Le décentrement (IDCEN/)
* 0/ Pas de décentrement (Galerkin)
* 1/ SUPG (Petrov-Galerkin)
* / SUPGDC (En dernier recours)
*------------------------------------------------------------------
IRESO = 2 ;
ICHOI = 3 ;
IDCEN = 0 ;
*
* On force MACRO CENTRE dans le cas semi-implicite
'SI' ('EGA' IRESO 4) ;
'MESS' 'On force MACRO CENTRE dans le cas semi-explicite' ;
ICHOI = 2 ;
IDCEN = 0 ;
'FINSI' ;
'SI' ('EGA' ICHOI 1) ;
DISCR = 'LINE' ;
KPRES = 'CENTRE' ;
'FINSI' ;
'SI' ('EGA' ICHOI 2) ;
DISCR = 'MACRO' ;
KPRES = 'CENTRE' ;
'FINSI' ;
'SI' ('EGA' ICHOI 3) ;
DISCR = 'MACRO' ;
KPRES = 'CENTREP1' ;
'FINSI' ;
'SI' ('EGA' ICHOI 4) ;
DISCR = 'QUAF' ;
KPRES = 'CENTREP1' ;
'FINSI' ;
'SI' ('EGA' ICHOI 5) ;
DISCR = 'QUAF' ;
KPRES = 'MSOMMET' ;
'FINSI' ;
'SI' ('EGA' IDCEN 0) ;
KSUPG = 'CENTREE' ;
'SINON' ;
KSUPG = 'SUPG' ;
* KSUPG = 'SUPGDC' ;
'FINSI' ;
*
*==========================================================
* Maillage
*==========================================================
*
* Dimensions caractéristiques
L = 1.25 ;
H = 2.5 ;
A = H '/' L ;
AS2 = A '/' 2. ;
FLAG1 = ICHOI < 4 ;
*
* Décallage par rapport à l'axe pour les éléments quadratiques
'SI' FLAG1 ;
RMIN = 0.0 ;
'SINO' ;
RMIN = 0.01 ;
'FINSI' ;
*
* Points du maillage
P0 = 0.5 AS2 ;
P1 = RMIN 0.0 ;
P2 = 0.5 0.0 ;
P3 = 1.0 0.0 ;
P4 = 1.0 AS2 ;
P5 = 1.0 A ;
P6 = 0.5 A ;
P7 = RMIN A ;
P8 = RMIN AS2 ;
*
* Données pour le mailleur
NXY = -4 ;
RAF = 2 ;
NX0 = RAF * NXY ;
NX = RAF * NXY ;
NY = RAF * NXY ;
D0 = 0.10 ;
D1 = 0.10 ;
D2 = 0.10 ;
*
* Droites du maillage filaire
P1P2 = P1 'DROI' NX0 P2 'DINI' D0 'DFIN' D2 ;
P2P3 = P2 'DROI' NX P3 'DINI' D2 'DFIN' D1 ;
P1P3 = P1P2 'ET' P2P3 ;
P3P4 = P3 'DROI' NY P4 'DINI' D1 'DFIN' D2 ;
P4P5 = P4 'DROI' NY P5 'DINI' D2 'DFIN' D1 ;
P3P5 = P3P4 'ET' P4P5 ;
P5P6 = P5 'DROI' NX P6 'DINI' D1 'DFIN' D2 ;
P6P7 = P6 'DROI' NX0 P7 'DINI' D2 'DFIN' D0 ;
P5P7 = P5P6 'ET' P6P7 ;
P7P8 = P7 'DROI' NY P8 'DINI' D1 'DFIN' D2 ;
P8P1 = P8 'DROI' NY P1 'DINI' D2 'DFIN' D1 ;
P7P1 = P7P8 'ET' P8P1 ;
*
* Maillages, sous-maillages et modèles associés
CNT1 = P1P3 'ET' P3P5 'ET' P5P7 'ET' P7P1 ;
CNT2 = P1P3 'ET' P3P5 'ET' P5P7 ;
DOM1 = 'DALL' P1P3 P3P5 P5P7 P7P1 ;
DOM1 = 'CHAN' DOM1 'QUAF' ;
*
$DOM1 = 'MODE' DOM1 'NAVIER_STOKES' DISCR ;
$CNT1 = 'MODE' CNT1 'NAVIER_STOKES' DISCR ;
$CNT2 = 'MODE' CNT2 'NAVIER_STOKES' DISCR ;
*
CNT1 = 'DOMA' $CNT1 'CENTRE' ;
CNT2 = 'DOMA' $CNT2 'CENTRE' ;
DOMF = 'DOMA' $DOM1 'FACE' ;
'ELIM' DOMF (CNT1 'ET' CNT2) 1.D-4 ;
DOM1 = 'DOMA' $DOM1 'MAILLAGE' ;
CNT1 = 'DOMA' $CNT1 'MAILLAGE' ;
CNT2 = 'DOMA' $CNT2 'MAILLAGE' ;
MP1 = ('DOMA' $DOM1 KPRES) 'ELEM' 1 ;
******** ????????????????? Menvlp = doma $vc 'ENVELOPPE' ;
******** ????????????????? nentr = doma $entree 'NORMALEV' ;
'DOMA' $DOM1 'IMPR' ;
*
*==========================================================
* Définition des équations vitesse, pression et température
*==========================================================
*
*
* Paramètres physiques et pas de temps
Pr = 0.65 ;
Gr = 4.4E4 ;
**Gr = 4.4E5 ;
GB1 = Gr * Pr * Pr ;
NU = Pr ;
ALFA = 1. ;
GB = 0. (-1. * GB1) ;
DT = 0.05 ;
*
* Equations en vitesse et température
* ITMA : Nombre de pas de temps
* NITER : Nombre d'itérations internes
* OMEGA : Facteur de relaxation des itérations internes
* FIDT : Nombre maximum de pas de temps
*RV = 'EQEX' $DOM1 'ITMA' 100 'NITER' 10 'OMEGA' 0.9 'FIDT' 10000
'SI' ('EGA' IRESO 1) ;
RV = 'EQEX' $DOM1 'ITMA' 500 'NITER' 5 'OMEGA' 0.2 'FIDT' 10000
'OPTI' 'EF' KSUPG 'IMPL' KPRES
'ZONE' $DOM1 'OPER' residu 1 IRESU
'ZONE' $DOM1 'OPER' 'NS' NU GB 'TN' 0. 'INCO' 'UN'
'ZONE' $DOM1 'OPER' 'LAPN' ALFA 'INCO' 'TN'
'ZONE' $DOM1 'OPER' 'KONV' 1. 'UN' ALFA 'INCO' 'TN'
'ZONE' $DOM1 'OPER' 'FIMP' 1. 'INCO' 'TN'
'OPTI' 'EF' 'CENTREE'
'ZONE' $DOM1 'OPER' 'DFDT' 1. 'UNM' DT 'INCO' 'UN'
'ZONE' $DOM1 'OPER' 'DFDT' 1. 'TNM' DT 'INCO' 'TN' ;
'FINSI' ;
'SI' ('EGA' IRESO 2) ;
RV = 'EQEX' $DOM1 'ITMA' 1000 'NITER' 1 'OMEGA' 1. 'FIDT' 10000
'OPTI' 'EF' KSUPG 'IMPL' KPRES
'ZONE' $DOM1 'OPER' residu 1 IRESU
'ZONE' $DOM1 'OPER' 'NS' NU GB 'TN' 0. 'INCO' 'UN'
'ZONE' $DOM1 'OPER' 'LAPN' ALFA 'INCO' 'TN'
'ZONE' $DOM1 'OPER' 'KONV' 1. 'UN' ALFA 'INCO' 'TN'
'ZONE' $DOM1 'OPER' 'FIMP' 1. 'INCO' 'TN'
'OPTI' 'EF' 'CENTREE'
'ZONE' $DOM1 'OPER' 'DFDT' 1. 'UN' DT 'INCO' 'UN'
'ZONE' $DOM1 'OPER' 'DFDT' 1. 'TN' DT 'INCO' 'TN' ;
'FINSI' ;
'SI' ('EGA' IRESO 3) ;
RV = 'EQEX' $DOM1 'ITMA' 1000 'NITER' 1 'OMEGA' 1. 'FIDT' 10000
'OPTI' 'EF' KSUPG 'IMPL' KPRES
'ZONE' $DOM1 'OPER' residu 1 IRESU
'ZONE' $DOM1 'OPER' 'NS' NU GB 'TN' 0. 'INCO' 'UN'
'ZONE' $DOM1 'OPER' 'LAPN' ALFA 'INCO' 'TN'
'ZONE' $DOM1 'OPER' 'KONV' 1. 'UN' ALFA 'INCO' 'TN'
'ZONE' $DOM1 'OPER' 'FIMP' 1. 'INCO' 'TN'
** 'ZONE' $DOM1 'OPER' 'TSCAL' ALFA 'UN' 1. 'INCO' 'TN'
'OPTI' 'EFM1' 'CENTREE'
'ZONE' $DOM1 'OPER' 'DFDT' 1. 'UN' DT 'INCO' 'UN'
'ZONE' $DOM1 'OPER' 'DFDT' 1. 'TN' DT 'INCO' 'TN' ;
'FINSI' ;
'SI' ('EGA' IRESO 4) ;
RV = 'EQEX' $DOM1 'ITMA' 1000 'ALFA' 1.
'ZONE' $DOM1 'OPER' residu 1 IRESU
'OPTI' 'SUPG'
'ZONE' $DOM1 'OPER' 'NS' NU GB 'TN' 0. 'INCO' 'UN'
'OPTI' 'SUPG'
'ZONE' $DOM1 'OPER' 'TSCAL' ALFA 'UN' 1. 'INCO' 'TN'
'OPTI' 'CENTREE'
'ZONE' $DOM1 'OPER' 'DFDT' 1. 'UN' 'DELTAT' 'INCO' 'UN'
'ZONE' $DOM1 'OPER' 'DFDT' 1. 'TN' 'DELTAT' 'INCO' 'TN' ;
'FINSI' ;
'SI' ('EGA' IRESO 5) ;
RV = 'EQEX' $DOM1 'ITMA' 1000 'ALFA' 1.
'ZONE' $DOM1 'OPER' residu 1 IRESU
'OPTI' 'SUPG'
'ZONE' $DOM1 'OPER' 'NS' NU GB 'TN' 0. 'INCO' 'UN'
'OPTI' 'SUPG'
'ZONE' $DOM1 'OPER' 'TSCAL' ALFA 'UN' 1. 'INCO' 'TN'
'OPTI' 'CENTREE'
'ZONE' $DOM1 'OPER' 'DFDT' 1. 'UN' 'DELTAT' 'INCO' 'UN'
'ZONE' $DOM1 'OPER' 'DFDT' 1. 'TN' 'DELTAT' 'INCO' 'TN' ;
'FINSI' ;
*
* Conditions aux limites en vitesse et température
RV = 'EQEX' RV 'CLIM'
'TN' 'TIMP' CNT2 0.
'UN' 'UIMP' CNT1 0. 'UN' 'VIMP' CNT2 0. ;
*
* Equation en pression avec condition de Dirichlet en un point
'SI' (' Début des tracés
'SI' GRAPH ;
trace psi1 DOM1 CNT1 NLI0 vun 'TITR' 'Fonction de courant' ;
*
* Température
tn = RV . 'INCO' . 'TN' ;
trace tn DOM1 CNT1 NLI0 'TITR' 'Température' ;
*
* Vitesse
un = RV . 'INCO' . 'UN' ;
vun = 'VECT' UN 0.005 'UX' 'UY' 'JAUNE' ;
trace VUN DOM1 CNT1 'TITR' 'Vitesse' ;
*
* Pression
'SI' ('EG' IRESO 4) ;
pe = RV . 'INCO' . 'PRESSION' ;
'FINSI' ;
pn = 'ELNO' $DOM1 ('KCHT' $DOM1 'SCAL' kpres pe) kpres;
trace pn dom1 CNT1 NLI0 'TITR' 'Pression' ;
*
'FINSI' ;
*-------------------------> Fin des tracés
*
*==========================================================
* Test de non régression
*==========================================================
*
'SAUT' 5 'LIGN' ;
VTMAX = maxi (RV . 'INCO' . 'TN') ;
mininu = 'MINI' nug ; maxinu = 'MAXI' nug ;
VTIME = RV . 'PASDETPS' . 'TPS' ;
'SI' COMPLET ;
TTMAX = 0.12727 ;
TNMAX =-3.52366E-2 ;
TNMIN =-0.99271 ;
TNAVE =-5.9121 ;
TIMEM = 3.25 ;
'SINON' ;
TTMAX = 0.12813 ;
TNMAX =-3.50596E-2 ;
TNMIN =-0.96950 ;
TNAVE =-5.8951 ;
TIMEM = 0.5 ;
'FINSI' ;
DTMAX = 1.D-5 ;
DNMAX = 1.D-4 ;
DTIME = 1.D-4 ;
ER1 = VTMAX - TTMAX ABS ; TEST1 = ER1 < DTMAX ;
ER2 = MININU - TNMIN ABS ; TEST2 = ER2 < DNMAX ;
ER3 = RESG - TNAVE ABS ; TEST3 = ER3 < DNMAX ;
ER4 = VTIME - TIMEM ABS ; TEST4 = ER4 < DTIME ;
'MESS' 'CHAMP/VALEUR/CIBLE/ERREUR/TOLERANCE/TEST' ;
'MESS' '----------------------------------------' ;
'MESS' 'Temperature ' VTMAX TTMAX ER1 DTMAX ;
list TEST1 ;
'MESS' 'Nusselt min ' mininu TNMIN ER2 DNMAX ;
list TEST2 ;
'MESS' 'Nusselt moy ' resg TNMIN ER3 DNMAX ;
list TEST3 ;
'MESS' 'Temps final ' VTIME TIMEM ER4 DTIME ;
list TEST4 ;
TEST5 = TEST1 ET TEST2 ET TEST3 ET TEST4 ;
SI TEST5 ;
ERRE 0 ;
SINO ;
ERRE 5 ;
FINSI ;
'FIN' ;
|
Fig 19 -  Température
|
|
|
Fig 20 - Vitesse et ligne de courant
|
|
|
Fig 21 - Pression
|
traduction
2003-11-04
