* fichier : cvry-2D-1.dgibi
****************** CONVECTION-RAYONNEMENT ************************* *
* *
* Couplage convection naturelle laminaire/rayonnement milieu absorbant* *
* *
* Il s'agit du test de De Vahl&Davis -cavité carrée- dans lequel le *
* milieu contenu dans la cavité est radiativement absorbant. * *
* Les parois de la cavité sont supposées grises et diffusantes. * *
* *
* Les données correspondent à celles de 2 cas tests de la littérature * *
* pour le test de validation: Rayleigh=1e5 et Nombre_rayonnement=1. * *
* *
* Formulation EF *
* Algorithme: implicite en régime permanent *
* convection traitée par pénalisation *
* rayonnement traité par la methode des harmoniques *
* sphériques (P1) *
* *
* Cet algorithme est limité aux nombres de Rayleigh faibles *
* (inférieurs à 1.e5 pour le test de Tan). Au dela d'autres *
* algorithnes sont nécessaires. *
* *
* opérateurs utilisés : DUDW, KONV, LAPN, MDIA *
* * *
* References: * *
* Z.Tan J.R.Howell, I.J.H.M.T Vol.34 1991 * *
* A.Yucel S.Acharya M.L.Williams, Num.Heat.Transfer Vol.15 1989 * *
* Rapports DMT/96.030 et DMT/96.059 *
* *
***********************************************************************
tan = vrai ;
complet = faux ;
graph = faux ;
option echo 0 dime 2 elem qua4 ;
* test de Tan
si tan ;
mess ' test de Tan ' ;
Ra = 1.e5 ;
Nr = 1.0 ; om = 0. ;
phi0 = 1. ; th = 2. ; tc = 1. ;
Pr = 0.71 ; Pe = 1. ; Re = Pe/Pr ;
* test de Yucel
sinon ;
mess ' test de Yucel';
Ra = 5.e6 ;
Nr = 0.04 ; om = 0. ;
phi0 = 1.5 ; th = 4./3. ; tc = 2./3. ;
Pr = 0.72 ; Re = 1. ; Pe = Re * Pr ;
finsi ;
* valeurs de reference pour les Nusselt - version du 7 nov. 97
* test de Tan
si complet ;
Nu_Tref = 3.10 ; Nu_Iref = 9.50 ;
sinon;
Nu_Tref = 3.00 ; Nu_Iref = 9.41 ;
finsi ;
OMEGA = 0.6 ;
EPS = 1.E-4 ;
EPSS = 1.E-10;
maxiter = 30 ;
*------------------ geometrie -------------------*
p1= 0. 0. ;
p2= 0.5 0. ;
p3= 0.5 0.5 ;
p4= 0. 0.5 ;
p5= 1. 0. ;
p6= 0. 1. ;
* maill. fin
*di = 0.005 ; df = 0.055 ;
*vi = 0.011 ; vf = df ;
si complet ;
* maill. moyen
di = 0.005 ; df = 0.1 ;
vi = 0.02 ; vf = df ;
sinon ;
* maill. grossier
di = 0.01 ; df = 0.2 ;
vi = 0.1 ; vf = df ;
finsi;
* ligne horiz.
c1 = droit p1 p2 dini di dfin df ;
c1b = droit p2 p5 dini df dfin di ;
* ligne vert.
c4 = droit p6 p4 dini vi dfin vf ;
c4b = droit p4 p1 dini vf dfin vi ;
l1 = c1 et c1b ;
l4 = c4 et c4b ;
l3inv= l1 plus (0. 1.) ; l3 = inve ( l3inv ) ;
l2inv= l4 plus (1. 0.) ; l2 = inve ( l2inv ) ;
fron = l1 et l2 et l3 et l4 ;
l13 = l1 et l3 ;
elim 1e-5 (l1 et l2 et l3 et l4 ) ;
mt = l1 l2 l3 l4 dalle plan ;
elim 1.e-5 mt ;
tass mt ;
orient mt ;
si graph ; trac mt ; finsi ;
* ------------- tables domaine ----------------- *
$mt = doma mt impr ;
$l13 = doma l13 'INCL' $mt 0.0001 impr ;
$fron = doma fron 'INCL' $mt 0.0001 impr ;
* ------------- parametres physiques ----------- *
* Re : nombre de Reynolds
* Ra : nombre de Rayleigh
* Pr : nombre de Prandtl
* Pe : nombre de Peclet = Pr * Re
* Nr : nombre caracteristique du rayonnement
* lambda: conductivite
* h : coefficient d'echange
* t0 : temperature de reference
* th : hot temperature
* tc : cold temperature
* emis : emissivite des parois
* a : rapport d'allongement (largeur/hauteur)
* beta : coefficient d'extinction (coeff.diffus + coeff.absorb)
* tau : epaisseur optique (beta * l)
* om : albedo ( coeff.diffus /(coeff.diffus + coeff.absorb) )
* phi0 : paramètre de forme
*
* ------------- parametres numeriques ---------- *
* OMEGA : facteur de relaxation
* EPS : precision sur les itérations internes
* EPSS : facteur de penalisation
a = 1. ; tau = 1. ;
emis = 1. ;
t0 = (th+tc)/2. ;
Gr = Ra/Pr ; gg = Gr/(Re*Re) ;
ALFV = 1./Re ;
ALFM = 1 ;
* ---------- coefficients algorithme ----------- *
E = emis/(4.-2.*emis) ;
H = (3.*(a*tau)) ;
H1 = (3.*((a*tau)**2))*(1.-om) ;
HE = H*E ;
3Nr = 3.*Nr*phi0 ;
C1 = (-1.*H1)/3Nr ;
C2 = (-1.*HE)/3Nr ;
* ---------- initialisation de la température--- *
RT = EQEX $MT 'OPTI' 'EF' 'IMPL'
ZONE $MT 'OPER' LAPN 1. INCO 'TT'
;
RT = EQEX RT
CLIM
'TT' TIMP l4 th
'TT' TIMP l2 tc
;
rt.INCO = TABLE 'INCO' ;
rt.INCO.'TT' = kcht $mt scal sommet t0 ;
rt.'NITER' = 1 ;
rt.'OMEGA' = 1. ;
exec rt ;
* ----------------- algorithme------------------ *
RV = EQEX $MT 'OPTI' 'EF' 'IMPL'
ZONE $MT 'OPER' DUDW EPSS INCO 'UN'
ZONE $MT 'OPER' KONV 1. 'UN' ALFV INCO 'UN'
ZONE $MT 'OPER' LAPN ALFV INCO 'UN'
ZONE $MT 'OPER' MDIA 'GB' INCO 'TN' 'UN'
;
RV = EQEX RV 'OPTI' 'EF' 'IMPL'
ZONE $MT 'OPER ' KONV Pe 'UN' ALFM INCO 'TN'
ZONE $MT 'OPER' LAPN ALFM INCO 'TN'
ZONE $MT 'OPER' MDIA 'C1' INCO 'IN' 'TN'
ZONE $MT 'OPER' MDIA 'KN' INCO 'TN'
;
RV = EQEX RV 'OPTI' 'EF' 'IMPL'
ZONE $MT 'OPER' LAPN 1. INCO 'IN'
ZONE $MT 'OPER' MDIA 'H1' INCO 'IN'
ZONE $MT 'OPER' MDIA 'K' INCO 'TN' 'IN'
;
RV = EQEX RV 'OPTI' 'EF' 'IMPL'
ZONE $fron 'OPER' MDIA HE INCO 'IN'
ZONE $fron 'OPER' MDIA 'KI' INCO 'TN' 'IN'
ZONE $l13 'OPER' MDIA 'KT' INCO 'TN'
ZONE $l13 'OPER' MDIA 'C2' INCO 'IN' 'TN'
;
RV = EQEX RV
CLIM
'UN' UIMP fron 0.0
'UN' VIMP fron 0.0
'TN' TIMP l4 th
'TN' TIMP l2 tc
;
rv.INCO = TABLE 'INCO' ;
rv.INCO.'TN' = kcht $mt scal sommet rt.INCO.'TT' ;
rv.INCO.'IN' = kcht $mt scal sommet 1. ;
rv.INCO.'GB' = kcht $mt vect sommet (0. (-1.*phi0*gg)) ;
rv.INCO.'UN' = kcht $mt vect sommet (1.e-5 1.e-5) ;
rv.'NITER' = 1 ;
rv.'OMEGA' = OMEGA ;
rv.'IMPR' = 1 ;
rv.'ITMA' = 1 ;
rv.INCO.'C1'= C1;
rv.INCO.'C2'= C2;
rv.INCO.'H1'= H1;
repeter bloc2 ;
nbiter = &bloc2 ;
* -----------actualisation des coefficients----- *
TN3 = kops (rv.INCO.'TN') ** 3 ;
4H1 = -4.*H1 ;
rv.INCO.'K' = kops 4H1 * TN3 ;
rv.INCO.'KN'= kops ((-1. * 4H1) / 3Nr ) * TN3 ;
4HE = -4.*HE ;
4HEN = (4.*HE)/3Nr ;
rv.INCO.'KI' = kcht $fron scal sommet (kops 4HE * TN3) ;
rv.INCO.'KT' = kcht $l13 scal sommet (kops 4HEN * TN3) ;
detr TN3 ;
UN = copie rv.INCO.'UN' ;
TN = copie rv.INCO.'TN' ;
IN = copie rv.INCO.'IN' ;
exec rv ;
*XU = (maxi (abs (kops ((rv.INCO.'UN') - UN) / UN))) ;
XT = (maxi (abs (kops ((rv.INCO.'TN') - TN) / TN))) ;
XI = (maxi (abs (kops ((rv.INCO.'IN') - IN) / IN))) ;
MESS ' ERREUR RELATIVE TN ' XT 'IN ' XI ;
Tmax = maxi (rv.INCO.'TN') ;
Tmin = mini (rv.INCO.'TN') ;
*MESS ' EXTREMUM TN MAX' Tmax 'MIN ' Tmin ;
* mess ' iteration ' nbiter ;
si ((XT < EPS) et (XI < EPS)) ;
quit bloc2 ;
sinon ;
si (nbiter > maxiter) ;
mess ' pas de convergence' ; erre 5 ;
finsi ;
finsi ;
fin bloc2 ;
*-------------post-traitement------------------------------------------
* Nusselt sur la paroi chaude
* ---------------------------
*
* calcul des Nusselt dus a la conduction et au rayonnement
* sur la paroi chaude (l4)
$l4 = doma l4 ;
* Nusselt du a la conduction
gradt = KOPS rv.'INCO'.'TN' 'GRAD' $mt ;
dTdXe = KCHT $mt 'SCAL' 'CENTRE' (exco 'UX' gradt) ;
dTdX = ELNO $MT (kops (-1.) * dTdXe) ;
Nu_T = KCHT $l4 'SCAL' 'SOMMET' (kops phi0 * dTdX) ;
* Nusselt du au rayonnement
I_l4 = redu rv.INCO.'IN' l4 ;
T_l4 = redu rv.INCO.'TN' l4 ;
T4 = kops T_l4 ** 4 ;
4T4 = kops 4 * T4 ;
dIdn = kops HE * (4T4 - I_l4) ;
Nu_I = KCHT $l4 'SCAL' 'SOMMET' (kops dIdn / (3.*Nr)) ;
* Nusselt global pour la paroi
dy = 'DOMA' $l4 'VOLUME' ;
dyt = SOMT dy ;
*message 'LONGUEUR DE LA PAROI = ' dyt ;
Nu_Te = noel $l4 Nu_T ;
Nu_Tdy = kops Nu_Te '*' dy ;
Nu_Tdyt = somt Nu_Tdy ;
Nu_Tm = Nu_Tdyt / dyt ;
message 'Nusselt de conduction moyen: ' Nu_Tm ;
Nu_Ie = noel $l4 Nu_I ;
Nu_Idy = kops Nu_Ie '*' dy ;
Nu_Idyt = somt Nu_Idy ;
Nu_Im = Nu_Idyt / dyt ;
message 'Nusselt de rayonnement moyen: ' Nu_Im ;
* comparaison aux valeurs de reference (version du 6 nov.97)
resT = abs((Nu_Tm-Nu_Tref/Nu_Tref));
resI = abs((Nu_Im-Nu_Iref/Nu_Iref));
mess 'resT' resT 'resI' resI ;
si ((resT < 1e-2) et (resI < 1e-2 )) ; erre 0 ;
sinon ; erre 5 ; finsi ;
*-------------graphique------------------------------------------------
si graph ;
si tan ;
* test de Tan
* opti titre 'CONDUCTION - RAYONNEMENT Ra = 1.e5 Nr = 1. ';
f_homo = 0.001 ;
sinon ;
* test de Yucel
* opti titre 'CONVECTION - RAYONNEMENT Ra = 5.e6 Nr = 0.04 ';
f_homo = 0.0003 ;
finsi ;
* température relative comprise entre -0.5 et 0.5
TY = kops rv.INCO.'TN' - T0 ;
TX = kops phi0 * TY ;
LTRA = prog -0.5 pas 0.1 0.5 ;
*trac TX mt (cont mt) LTRA ;
opti isov ligne ;
trac TX mt (cont mt) LTRA ;
* intensite radiative
trac rv.INCO.'IN' 10 mt (cont mt) ;
* vitesse
v1 = vect (rv.INCO.'UN') f_homo ux uy jaune ;
trac v1 mt (cont mt) ;
T1 = evol 'CHPO' rv.INCO.'TN' l1 ;
T3 = evol 'CHPO' rv.INCO.'TN' l3inv ;
*TAB1=TABLE ;
*TAB1.1= 'MARQ CROI REGU TITR PAROI_INFERIEURE' ;
*TAB1.2= 'MARQ TRIA REGU TITR PAROI_SUPERIEURE' ;
dess (T1 et T3) MIMA ;
* evolution des Nusselt sur la paroi chaude
Nu_Tl = EVOL 'CHPO' Nu_T SCAL l4 ;
*TAB2 = TABLE ;
*TAB2.1 = 'MARQ LOSA' ;
DESS Nu_Tl 'MIMA' ;
*DESS Nu_Tl 'MIMA' TAB2
* 'TITR' ' Ra =1.e5 Nombre de Nusselt local sur la paroi chaude '
* 'TITX' ' y ' 'TITY' ' Nusselt ' ;
Nu_Il = evol 'CHPO' Nu_I l4 ;
dess Nu_Il 'MIMA' ;
* Nusselt total
Nu_l = Nu_Tl + Nu_Il ;
*TAB3 = TABLE ;
*TAB3.1 = 'MARQ CROI REGU ' ;
*TAB3.2 = 'MARQ TRIA REGU ';
*TAB3.3 = 'MARQ LOSA REGU ' ;
*TAB3.'TITRE' = TABLE ;
*TAB3.'TITRE'. 1 = MOT 'CONDUCTION' ;
*TAB3.'TITRE'. 2 = MOT 'RAYONNEMENT';
*TAB3.'TITRE'. 3 = MOT 'TOTAL' ;
dess (Nu_Tl et Nu_Il et Nu_l) ;
*dess (Nu_Tl et Nu_Il et Nu_l) LEGE TAB3
* 'TITR' 'Nombre de Nusselt local sur la paroi chaude '
* 'TITX' ' y ' 'TITY' ' Nusselt ' ;
finsi ;
fin ;
© Cast3M 2003 - All rights reserved.
Disclaimer