Télécharger cvry-2D-1.dgibi

Retour à la liste

* 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 - Tous droits réservés.
Mentions légales