***************************************************** ************************************************************************ ************************************************************************ * fichier : tp4.dgibi * ** modifie le 15/06/2014 passage EQPR -> EQEX * ***************************************************** ***************** 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 ; TEST = 'EGA' VAL0 2 ; 'SI' TEST ; 'SINON' ; 'QUIT' courant ; 'FINSI' ; TEST = 'EGA' VAL1 'AXIS' ; 'SI' TEST ; ROTU0 = 2. * VZ1 '-' (ROTU0 * XC1) ; 'SINO' ; 'FINSI' ; 'ZONE' $domaine 'OPER' 'FIMP' ROTU0 'INCO' 'PSI' 'CLIM' 'PSI' 'TIMP' CONT0 0. ; RK . 'INCO' = 'TABLE' 'INCO' ; EXEC RK ; '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 ; NUSS2 = NUSS1 * S1 ; '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 ; 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 ; Y1 = ('LOG' EPS0) '/' ('LOG' 10) ; 'SI' POST1 ; X1 = 0. ; X2 = RV . 'ITMA' ; 'DESSIN' EV1 'YBOR' Y1 0. 'NCLK' ; 'FINSI' ; 'SI' POST2 ; 'FINSI' ; 'SI' ((ERR10 < Y1) 'ET' (DD > 10)) ; RV . 'TFINAL' = RV . 'PASDETPS' . 'TPS' ; 'FINSI' ; 'FINSI' ; RV . 'INCO' . 'TN2' = 'COPIER' RV . 'INCO' . 'TN' ; 'FINP' mat1 chp1 ; *------------------------------------------------------------------ * * * * 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 * 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 ; * '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 P1P3 = P1P2 'ET' P2P3 ; P3P5 = P3P4 'ET' P4P5 ; P5P7 = P5P6 'ET' P6P7 ; 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 ; * * ******** ????????????????? 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) ; 'OPTI' 'EF' KSUPG 'IMPL' KPRES 'ZONE' $DOM1 'OPER' residu 1 IRESU 'ZONE' $DOM1 'OPER' 'LAPN' ALFA 'INCO' 'TN' 'ZONE' $DOM1 'OPER' 'KONV' 1. 'UN' ALFA 'INCO' 'TN' RV = 'EQEX' RV 'OPTI' 'EF' 'CENTREE' 'ZONE' $DOM1 'OPER' 'DFDT' 1. 'UNM' DT 'INCO' 'UN' 'FINSI' ; 'SI' ('EGA' IRESO 2) ; 'OPTI' 'EF' KSUPG 'IMPL' KPRES 'ZONE' $DOM1 'OPER' residu 1 IRESU 'ZONE' $DOM1 'OPER' 'LAPN' ALFA 'INCO' 'TN' 'ZONE' $DOM1 'OPER' 'KONV' 1. 'UN' ALFA 'INCO' 'TN' RV = 'EQEX' RV 'OPTI' 'EF' 'CENTREE' 'ZONE' $DOM1 'OPER' 'DFDT' 1. 'UN' DT 'INCO' 'UN' 'FINSI' ; 'SI' ('EGA' IRESO 3) ; 'OPTI' 'EF' KSUPG 'IMPL' KPRES 'ZONE' $DOM1 'OPER' residu 1 IRESU 'ZONE' $DOM1 'OPER' 'LAPN' ALFA 'INCO' 'TN' 'ZONE' $DOM1 'OPER' 'KONV' 1. 'UN' ALFA 'INCO' 'TN' ** 'ZONE' $DOM1 'OPER' 'TSCAL' 1. 'UN' ALFA 1. 'INCO' 'TN' RV = 'EQEX' RV 'OPTI' 'EFM1' 'CENTREE' 'ZONE' $DOM1 'OPER' 'DFDT' 1. 'UN' DT 'INCO' 'UN' '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' ('<EG' IRESO 2) ; 'FINSI' ; 'SI' ('EGA' IRESO 3) ; 'FINSI' ; * * Initialisation des champs (table INCO) rv . 'INCO' = TABLE 'INCO' ; * * Champs supplémentaires pour la procédure residu * * Méthode d'inversion du problème en vitesse rv . 'METHINV' . TYPINV = 1 ; rv . 'METHINV' . IMPINV = 0 ; rv . 'METHINV' . NITMAX = 100 ; rv . 'METHINV' . PRECOND = 3 ; rv . 'METHINV' . RESID = 1.e-6 ; rv. 'METHINV' . 'FCPRECT'=1 ; rv. 'METHINV' . 'FCPRECI'=1 ; * * Méthode d'inversion du problème en pression 'SI' ('<EG' IRESO 2) ; rvp . 'METHINV' . TYPINV = 1 ; rvp . 'METHINV' . IMPINV = 0 ; rvp . 'METHINV' . NITMAX = 100 ; rvp . 'METHINV' . PRECOND = 3 ; rvp . 'METHINV' . RESID = 1.e-12 ; rvp.'METHINV' . 'FCPRECT'=100 ; rvp.'METHINV' . 'FCPRECI'=100 ; 'FINSI' ; * * Couplage vitesse/pression : Méthode de projection 'SI' ('<EG' IRESO 2) ; 'FINSI' ; 'SI' ('EGA' IRESO 3) ; * En implicite complet, EXEC sans infos en + 'FINSI' ; * *========================================================== * Résolution *========================================================== * 'OPTION' echo 0 ; EXEC RV ; 'OPTION' echo 0 ; * *========================================================== * Post traitement *========================================================== * * Nusselt nug resg = calcnuss (RV . 'INCO' . 'TN') $DOM1 $CNT2 ; * * Fonction de courant psi1 = courant (RV . 'INCO' . 'UN') $DOM1 $CNT1 ; *-------------------------> Début des tracés 'SI' GRAPH ; trace psi1 DOM1 CNT1 NLI0 'TITR' 'Fonction de courant' ; * * Température tn = RV . 'INCO' . 'TN' ; * * Vitesse un = RV . 'INCO' . 'UN' ; trace VUN DOM1 CNT1 'TITR' 'Vitesse' ; * * Pression 'SI' ('<EG' IRESO 2) ; 'FINSI' ; 'SI' ('EGA' IRESO 3) ; 'FINSI' ; trace pn dom1 CNT1 NLI0 'TITR' 'Pression' ; * 'FINSI' ; *-------------------------> Fin des tracés * *========================================================== * Test de non régression *========================================================== * 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' '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 ; SINO ; FINSI ; 'FIN' ;
© Cast3M 2003 - Tous droits réservés.
Mentions légales