Télécharger cav_ray_proj.dgibi
* fichier : cav_ray_proj.dgibi ************************************************************************ ************************************************************************ *---------------------------------------------- * Couplage Rayonnement/Convection/Conduction * * Test d'un cavité rayonnante en 2D plan * Le mélange gazeux est transparent * Rayleigh = 1e6 * Résolution par la méthode de Projection * On introduit une paroi pour le rayonnement * afin de récupérer une constante de temps car * l'algorithme de projection est un algorithme * transitoire (DFDT) * Auteur : E. Studer - Mars 2002 * Référence : 'Numerical Analysis of Natural convection * with surface radiation in a square enclosure' * M. Akiyama and Q.P. Chong * Numerical Heat Transfer 31:419-433,1997 *---------------------------------------------- *---------------------------------------------- * Declarations générales *---------------------------------------------- COMPLET = FAUX ; GRAPH = FAUX ; **************** * PARAMETRES : * **************** * nombre d'itérations internes iter = 1 ; omeg = 1.0 ; * nombre d'itérations en temps 'SI' COMPLET ; itma = 400 ; d1 = 0.01 ; d2 = 0.05 ; * elements en vitesse DISCR = 'QUAF' ; 'SINON' ; itma = 20 ; d1 = 0.01 ; d2 = 0.1 ; * elements en vitesse DISCR = 'LINE' ; 'FINSI' ; * pas de temps DT = 1.0 ; * option de décentrement * KSUPG = 'CENTREE' ; KSUPG = 'SUPG' ; * support de la pression KPRES = 'MSOMMET' ; * KPRES = 'CENTRE' ; * KPRES = 'CENTREP1' ; * choix de la méthode de résolution IMPLICITE ou PROJECTION * RESOL = MOT IMPLICITE ; * prise en compte du rayonnement PRAYO = VRAI ; * PRAYO = FAUX ; * Ajout d'une paroi pour conduction et stabilisation PPAROI = VRAI ; * PPAROI = FAUX ; * Parametres du maillage: * Dimension de la cavité L = 1 ; * Parametres physiques: * EPSI = écart de température pour paroi chaude/froide * Température de référence T0 = 293.5 ; * T1 = Température chaude * T2 = Température froide * beta = 1.0 '/' T0 ; gbeta = -9.81 '*' beta ; gb = 0. gbeta ; V0 = ( 'ABS' ( (gbeta '*' (T1 '-' T2)) '*' L ) ) '**' 0.5 ; Pr = 0.71 ; Ra = 1.E6 ; GR = Ra '/' Pr ; Nu = V0 '*' L '/'(Gr '**' 0.5) ; alpha = Nu '/' Pr ; * conductivité de la paroi alphap = 1000.0 ; uref = 1 ; * Rayonnement: emissivité des surfaces 'SI' (PRAYO) ; e_sud = 0.5 ; e_nord = 0.5 ; e_est = 0.5 ; e_ouest = 0.5 ; 'FINSI' ; ************** * MAILLAGE : * ************** p0 = 0. 0. ; p1 = 0.5 0. ; p2 = 1. 0. ; p3 = 1. 0.5 ; p4 = 1. 1. ; p5 = 0.5 1. ; p6 = 0. 1. ; p7 = 0. 0.5 ; sud = d01 'ET' d12 ; est = d23 'ET' d34 ; nord = d45 'ET' d56 ; ouest = d67 'ET' d70 ; parois = (est 'ET' ouest 'ET' nord 'ET' sud) ; 'SI' (PPAROI) ; * ajout d'une paroi pour mettre de la conduction 'SI' (PRAYO) ; mparoi = mest 'ET' mouest 'ET' mnord 'ET' msud ; 'SINON' ; mparoi = mest 'ET' mouest ; 'FINSI' ; 'FINSI' ; * Maillage du carre fluide carre = sud est nord ouest 'DALLER' ; * Maillage total 'SI' (PPAROI) ; mt = carre 'ET' mparoi ; 'SINON' ; mt = carre ; 'FINSI' ; * Tracé des maillages 'SI' GRAPH ; 'FINSI' ; ************************** * MODELE NAVIER-STOKES : * ************************** 'SI' (PPAROI) ; 'SINON' ; 'FINSI' ; 'SI' (PPAROI) ; 'FINSI' ; 'SI' (PPAROI) ; 'FINSI' ; 'SI' (PPAROI) ; 'FINSI' ; 'SI' (PPAROI) ; 'ET' est 'ET' sud 'ET' ouest 'ET' nord) 1.D-4 ; 'ET' est 'ET' sud 'ET' ouest 'ET' nord) 1.D-4 ; 'SINON' ; 'FINSI' ; 'SI' (PPAROI) ; 'FINSI' ; * ***************** * MODELISATION : * ****************** 'SI' (PRAYO) ; * Création des modèles de rayonnement * assemblage des modèles de rayonnement mray = mrsud 'ET' mrest 'ET' mrnord 'ET' mrouest ; * création des matériaux émissivité * assemblage des matériaux e = esud 'ET' eest 'ET' enord 'ET' eouest ; ************************************************* * FACTEURS DE FORME et MATRICE DE RAYONNEMENT : * ************************************************* ********************************************************** * CREATION DE LA TABLE POUR L'OPERATEUR DE RAYONNEMENT : * ********************************************************** TAB_RAY = 'TABLE' ; TAB_RAY.'MAILLAGE' = (sud 'ET' est 'ET' nord 'ET' ouest) ; TAB_RAY.'MODELE' = mray ; TAB_RAY.'MATERIAU' = e ; TAB_RAY.'FFORM' = fft ; TAB_RAY.'MATRICE' = chamr ; * Modèles Navier-Stokes TAB_RAY.'MODELENS' = ($sud 'ET' $est 'ET' $nord 'ET' $ouest) ; 'FINSI' ; ****************************** * OPERATEUR DE RAYONNEMENT : * ****************************** 'DEBP' FRAY ; * récupération du champ de température T = RV.INCO.'TN' ; * récupération de la table TABR associée à l'opérateur TABR = RVX.'ARG1' ; MAILR = TABR.'MAILLAGE' ; * TSU : température de surface * MRS : matrice de rayonnement * MRS : transformation Rigidité/Matrik 'FINP' AS MRS ; ************************* * CREATION DES TABLES : * ************************* 'FIDT' 100 'OPTI' 'EF' KSUPG 'IMPL' KPRES 'ZONE' $carre 'OPER' 'NS' 1. 'UN' Nu gb 'TF' T0 'INCO' 'UN' ; 'OPTI' 'EF' KSUPG 'IMPL' KPRES 'ZONE' $carre 'OPER' 'TSCAL' 1. 'UN' alpha 0. 'INCO' 'TN' ; 'SI' (PPAROI) ; RT = 'EQEX' RT 'OPTI' 'EF' KSUPG 'IMPL' KPRES 'ZONE' $mparoi 'OPER' 'LAPN' alphap 'INCO' 'TN' ; 'FINSI' ; 'SI' ('EGA' RESOL 'PROJECTION') ; RV = 'EQEX' RV 'OPTI' 'EF' 'CENTREE' 'ZONE' $carre 'OPER' 'DFDT' 1.0 'UNM' DT 'INCO' 'UN' ; RT = 'EQEX' RT 'OPTI' 'EF' 'CENTREE' 'ZONE' $mt 'OPER' 'DFDT' 1.0 'TNM' DT 'INCO' 'TN' ; 'FINSI' ; SI ('EGA' RESOL 'IMPLICITE') ; RV = 'EQEX' RV 'OPTI' 'EFM1' 'CENTREE' 'ZONE' $carre 'OPER' 'DFDT' 1.0 'UNM' DT 'INCO' 'UN' ; RT = 'EQEX' RT 'OPTI' 'EFM1' 'CENTREE' 'ZONE' $mt 'OPER' 'DFDT' 1.0 'TNM' DT 'INCO' 'TN' ; 'FINSI' ; 'SI' (PRAYO) ; RT = 'EQEX' RT 'OPTI' 'EF' 'CENTREE' 'ZONE' ($sud 'ET' $est 'ET' $nord 'ET' $ouest) 'OPER' 'FRAY' TAB_RAY 'INCO' 'TN' ; 'FINSI' ; RV = 'EQEX' RV 'CLIM' 'UN' 'UIMP' parois 0. 'UN' 'VIMP' parois 0. ; 'SI' (PPAROI) ; RT = 'EQEX' RT 'CLIM' 'TN' 'TIMP' pouest T1 'TN' 'TIMP' pest T2 ; 'SINON' ; RT = 'EQEX' RT 'CLIM' 'TN' 'TIMP' ouest T1 'TN' 'TIMP' est T2 ; 'FINSI' ; RV.'INCO' = 'TABLE' 'INCO' ; RV.'INCO'.'NITER' = 0 ; RT.'INCO' = RV.'INCO' ; 'SI' ('EGA' RESOL 'PROJECTION') ; RVP = 'EQEX' 'OPTI' 'EF' KPRES 'FINSI' ; 'SI' ('EGA' RESOL 'IMPLICITE') ; RV = 'EQEX' RV 'OPTI' 'EF' 'IMPL' KPRES KSUPG 'FINSI' ; 'SI' ('EGA' RESOL 'PROJECTION') ; 'FINSI' ; rv.'METHINV'.TYPINV=1 ; rv.'METHINV'.IMPINV=0 ; rv.'METHINV'.NITMAX=400; rv.'METHINV'.PRECOND=3 ; rv.'METHINV'.RESID =1.e-10; rv. 'METHINV' . 'FCPRECT'=1 ; rv. 'METHINV' . 'FCPRECI'=1 ; rt.'METHINV'.TYPINV=1 ; rt.'METHINV'.IMPINV=0 ; rt.'METHINV'.NITMAX=400; rt.'METHINV'.PRECOND=3 ; rt.'METHINV'.RESID =1.e-10; rt. 'METHINV' . 'FCPRECT'=1 ; rt. 'METHINV' . 'FCPRECI'=1 ; 'SI' ('EGA' RESOL 'PROJECTION') ; rvp.'METHINV'.TYPINV=1 ; rvp.'METHINV'.IMPINV=0 ; rvp.'METHINV'.NITMAX=400; rvp.'METHINV'.PRECOND=3 ; rvp.'METHINV'.RESID =1.e-8 ; rvp.'METHINV' . 'FCPRECT'=100 ; rvp.'METHINV' . 'FCPRECI'=100 ; 'FINSI' ; **************** * RESOLUTION : * **************** 'REPETER' BOU1 ITMA ; * projection de la température fluide TN = RV.'INCO'.'TN' ; RV.'INCO'.'TF' = 'KCHT' $carre 'SCAL' 'SOMMET' RV.INCO.'NITER' = RV.INCO.'NITER' + 1 ; DD = RV.INCO.'NITER' ; NN = DD '/' 5; LO = (5 '*' NN '-' DD) 'EGA' 0 ; 'SI' (LO) ; UN = RV.INCO.'UN' ; UNM1 = RV.INCO.'UNM1' ; TN = RV.INCO.'TN' ; TNM1 = RV.INCO.'TNM1' ; ELIX = ('LOG' (ELIX '+' 1.e-10)) '/' ('LOG' 10.) ; ELIY = ('LOG' (ELIY '+' 1.e-10)) '/' ('LOG' 10.) ; ELIT = ('LOG' (ELIT '+' 1.e-10)) '/' (LOG 10.) ; RV.INCO.'IT' = (RV.INCO.'IT') 'ET' IT ; RV.INCO.'ERX' = (RV.INCO.'ERX') 'ET' ERX ; RV.INCO.'ERY' = (RV.INCO.'ERY') 'ET' ERY ; RV.INCO.'ERT' = (RV.INCO.'ERT') 'ET' ERT ; 'FINSI' ; EXEC RV ; EXEC RT ; 'FIN' BOU1 ; ****************************** * EXPLOITATION DES RESULTATS * ****************************** 'SI' GRAPH ; 'FINSI' ; * Adimensionnement des températures pour dessin 0.2 0.3 0.4 0.5 ; TN = ((RV.INCO.'TN') '-' T0 ) '/' (T1 '-' T2) ; 'SI' GRAPH ; 'FINSI' ; 'SI' GRAPH ; 'FINSI' ; * Fonction de courant un = RV.INCO.'UN' ; 'ZONE' $carre 'OPER' 'FIMP' sw INCO 'PSI' 'CLIM' 'PSI' 'TIMP' parois 0. ; rk.'INCO' = 'TABLE' 'INCO' ; EXEC rk ; psi = rk.INCO.'PSI' '/' V0 ; 'SI' GRAPH ; 'FINSI' ; (RV.INCO.'ERX') ; 'SI' GRAPH ; 'FINSI' ; 'SI' ('NON' COMPLET) ; * Test de non régression * Erreur de référence sur la température ELITM = -8.30554E-02 ; * Valeur de référence sur le maximum de la vitesse verticale UYMAX = 0.18220; * Valeur de référence sur le maximum de la fonction de courant PSIMAX = 27.652 ; * ERT = 'ABS'( ELIT '-' ELITM) '/' ELITM ; ERU = 'ABS'( UYM '-' UYMAX) '/' UYMAX ; ERP = 'ABS'( PSIM '-' PSIMAX) '/' PSIMAX ; 'SI' (ERT '>' 0.002); 'ERREUR' 5 ; 'FINSI' ; 'SI' (ERU '>' 0.002); 'ERREUR' 5 ; 'FINSI' ; 'SI' (ERP '>' 0.002); 'ERREUR' 5 ; 'FINSI' ; 'FINSI' ; 'FIN' ; * * *
© Cast3M 2003 - Tous droits réservés.
Mentions légales