* fichier :  wsgg.dgibi
******************  CONVECTION-RAYONNEMENT + VAPEUR D'EAU  ********** *
*                                                                     *
* Couplage convection naturelle laminaire/rayonnement milieu absorbant*
*                                                                     *
* Convection naturelle horizontale dans une cavite carree contenant   *
* de la vapeur d'eau a la pression de 1 atm.                          *
* La temperature des parois varie entre 1000K et 2000K.               *
*                                                                     *
* Par rapport au test de base cvry-2D-1.dgibi le milieu est un gas    *
* reel modelise par le modele WSGG (Weighted Sum of Gray Gases)       *
* dans lequel le gas est decompose en la somme de 3 gaz gris affectés *
* chacun d'un facteur de ponderation dependant polynomialement de     *
* la temperature et d'un coefficient d'absorption constant.           *
*                                                                     *
* Données: test de Tan avec Ra = 1.e4 et Nr = 1.                      *
* Validation numerique par comparaison avec Trio-EF                   *
*                                                                     *
* Formulation EF                                                      *
* Algorithme: implicite en régime permanent avec 6 inconnues          *
*             convection traitée par pénalisation                     *
*             rayonnement traité par la methode des harmoniques       *
*             sphériques (P1) + modele WSGG a 3 gaz                   *
*             non linearite des proprietes radiatives                 *
*                                                                     *
* opérateurs utilisés : DUDW, KONV, LAPN, MDIA                        *
*                                                                     *
* References:                                                         *
* - pour le couplage convection rayonnement et l'approximation P1     *
*    test cvry-2D-1.dgibi                                             *
*    Z.Tan J.R.Howell, I.J.H.M.T Vol.34 1991                          *
*    rapport DMT/96.030 et DMT/94.505                                 *
* - pour les proprietes radiatives de la vapeur d'eau:                *
*    T.H.Smith J.H.T ASME Vol.104 1982                                *
*    rapport DMT/94.505                                               *
*                                                                     *
***********************************************************************

complet = faux ;
post = vrai ; graph = faux ;

option echo 0 dime 2 elem qua4 ;

*----------------------------------
*   Parametres physiques
*----------------------------------

Ra  = 1.e4 ; Nr = 1.0 ;

om = 0. ;
Tref = 1000. ; L = 1. ;
phi0 = 1. ; th  = 2.  ; tc = 1. ;
Pr = 0.71 ; Pe = 1. ; Re = Pe/Pr ;
a = 1.  ;
emis = 1.  ;

* valeurs de reference pour les Nusselt - version du 8 jan. 98
* test de Tan
  Nu_Tref = 2.15 ; Nu_Iref = 13.1 ;

*----------------------------------
*   Parametres numeriques
*----------------------------------

OMEGA = 0.6   ;
EPSS  = 1.E-10;
maxiter = 60 ;
si complet; EPS = 1.E-4 ; sinon ; EPS = 1.E-2 ; finsi;

* ---------------------------------
* Re    : nombre de Reynolds
* Ra    : nombre de Rayleigh
* Pr    : nombre de Prandtl
* Pe    : nombre de Peclet = Pr * Re
* Nr    : nombre caracteristique du rayonnement

* Tref  : temperature de reference en degre K
* L     : longueur caracteristique
*         (pour le calcul des proprietes radiatives)
* t0    : temperature de reference pour le pb adimensionne
* th    : hot temperature 
* tc    : cold temperature
* emis  : emissivite des parois
* a     : rapport d'allongement (largeur/hauteur)
* ki    : coefficient d'absorption du gas i
* taui  : epaisseur optique (ki * L)
* om    : albedo  ( coeff.diffus /(coeff.diffus + coeff.absorb) )
* phi0  : paramètre de forme
*
* ---------------------------------
* OMEGA   : facteur de relaxation
* EPS     : precision sur les itérations internes
* EPSS    : facteur de penalisation
* maxiter : nombre maximum d'iterations

*----------------------------------
*   Proprietes radiatives
*----------------------------------
*     H2O Ptot=1atm
* decomposition de la vapeur d'eau en 3 gas gris
* definition des coefficient d'absorption ki et des coefficients
* pour le calcul des facteurs de ponderation wi
* taui designe l'epaisseur optique correspondante

  k1 = 0.4496 ; k2 = 7.113 ; k3 = 119.7 ; k4 = 1.e-5 ;
  tau1 = k1*L ; tau2 = k2*L ; tau3 = k3*L ; tau4 = k4*L ;

 a0= 6.324e-1 ; a1=-8.358E-4 ; a2= 6.135e-7 ; a3= -13.03e-11 ;
 b0=-2.016e-2 ; b1= 7.145E-4 ; b2=-5.212e-7 ; b3=  9.868e-11 ;
 c0= 3.500e-1 ; c1=-5.040E-4 ; c2= 2.425e-7 ; c3= -3.888e-11 ;

* utilisation de la procedure CWSGG pour le calcul des wi
*  w1 = a0+a1*T+a2*T*T+a3*T*T*T ;
*  w2 = b0+b1*T+b2*T*T+b3*T*T*T ;
*  w3 = c0+c1*T+c2*T*T+c3*T*T*T ;
*  w4= 1.0 -w1-w2-w3 ;

*---------------------------------------------------------------
 'DEBPROC' CWSGG T*'CHPOINT'
           a0*'FLOTTANT' a1*'FLOTTANT' a2*'FLOTTANT' a3*'FLOTTANT' ;
* calcul de la loi polynomiale donnant les facteurs de ponderation
* en fonction de la temperature en degre K
* w  = a0+a1*T+a2*T*T+a3*T*T*T ;
  b2 = 'KOPS' a2 + ('KOPS' T * a3) ;
  b1 = 'KOPS' a1 + ('KOPS' T * b2) ;
  w  = 'KOPS' a0 + ('KOPS' T * b1) ;
'FINPROC' w    ;
*---------------------------------------------------------------



*----------------------------------
*   Maillage
*----------------------------------

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.
 c1a=  droit p1 p2 dini di dfin df ;
 c1b = droit p2 p5 dini df dfin di ;

* ligne vert.
 c4a=  droit p6 p4 dini vi dfin vf ;
 c4b = droit p4 p1 dini vf dfin vi ;

l1 = c1a et c1b   ;
l4 = c4a 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 ;


*----------------------------------
*   Definition de coefficients divers
*----------------------------------

t0  = (th+tc)/2. ;
Gr = Ra/Pr ; gg = Gr/(Re*Re) ;

ALFV = 1./Re ;
ALFM = 1   ;

E = emis/(4.-2.*emis) ;
3Nr = 3. * Nr ;

3Ntau1 = (-1.0)/(3Nr*tau1) ;
3Ntau2 = (-1.0)/(3Nr*tau2) ;
3Ntau3 = (-1.0)/(3Nr*tau3) ;
3Ntau4 = (-1.0)/(3Nr*tau4) ;

ATE1 = (3.*E)*(a*tau1)   ;
ATT1 = (3.*((a*tau1)**2))*(1.-om) ;

ATE2 = (3.*E)*(a*tau2)   ;
ATT2 = (3.*((a*tau2)**2))*(1.-om) ;

ATE3 = (3.*E)*(a*tau3)   ;
ATT3 = (3.*((a*tau3)**2))*(1.-om) ;

ATE4 = (3.*E)*(a*tau4)   ;
ATT4 = (3.*((a*tau4)**2))*(1.-om) ;


*----------------------------------
*   Initialisation du champ de temperature
*----------------------------------

 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 ;


*----------------------------------
*   Definition du systeme lineaire
*----------------------------------

*  bilan masse et qdm (vitesse UN)
*  bilan energie      (temperature)
*  bilan sur les 4 intensites radiatives de la decomposition
*      (I1, I2, I3, I4)

 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'
;
*  inconnue T
 RV = EQEX  RV            'OPTI'  'EF'  'IMPL'
      ZONE  $MT   'OPER ' KONV  Pe  'UN'  ALFM  INCO 'TN'
      ZONE  $MT   'OPER'  LAPN   ALFM           INCO 'TN'
;
 RV = EQEX  RV            'OPTI'  'EF'  'IMPL'
      ZONE  $MT   'OPER'  MDIA  'CT1'           INCO 'TN'
      ZONE  $MT   'OPER'  MDIA  'CI1T'          INCO 'I1' 'TN'
      ZONE  $l13  'OPER'  MDIA  'CdT1'          INCO 'TN'
      ZONE  $l13  'OPER'  MDIA  'CdI1T'         INCO 'I1' 'TN'
;
 RV = EQEX  RV            'OPTI'  'EF'  'IMPL'
      ZONE  $MT   'OPER'  MDIA  'CT2'           INCO 'TN'
      ZONE  $MT   'OPER'  MDIA  'CI2T'          INCO 'I2' 'TN'
      ZONE  $l13  'OPER'  MDIA  'CdT2'          INCO 'TN'
      ZONE  $l13  'OPER'  MDIA  'CdI2T'         INCO 'I2' 'TN'
;
 RV = EQEX  RV            'OPTI'  'EF'  'IMPL'
      ZONE  $MT   'OPER'  MDIA  'CT3'           INCO 'TN'
      ZONE  $MT   'OPER'  MDIA  'CI3T'          INCO 'I3' 'TN'
      ZONE  $l13  'OPER'  MDIA  'CdT3'          INCO 'TN'
      ZONE  $l13  'OPER'  MDIA  'CdI3T'         INCO 'I3' 'TN'
;
 RV = EQEX  RV            'OPTI'  'EF'  'IMPL'
      ZONE  $MT   'OPER'  MDIA  'CT4'           INCO 'TN'
      ZONE  $MT   'OPER'  MDIA  'CI4T'          INCO 'I4' 'TN'
      ZONE  $l13  'OPER'  MDIA  'CdT4'          INCO 'TN'
      ZONE  $l13  'OPER'  MDIA  'CdI4T'         INCO 'I4' 'TN'
;
* inconnue I1
 RV = EQEX  RV            'OPTI'  'EF'  'IMPL'
      ZONE  $MT   'OPER'  LAPN   1.             INCO 'I1'
      ZONE  $MT   'OPER'  MDIA  'CI1'           INCO 'I1'
      ZONE  $MT   'OPER'  MDIA  'CTI1'          INCO 'TN' 'I1'
      ZONE  $fron 'OPER'  MDIA  'CdI1'          INCO 'I1'
      ZONE  $fron 'OPER'  MDIA  'CdTI1'         INCO 'TN' 'I1'
;
* inconnue I2
 RV = EQEX  RV            'OPTI'  'EF'  'IMPL'
      ZONE  $MT   'OPER'  LAPN   1.             INCO 'I2'
      ZONE  $MT   'OPER'  MDIA  'CI2'           INCO 'I2'
      ZONE  $MT   'OPER'  MDIA  'CTI2'          INCO 'TN' 'I2'
      ZONE  $fron 'OPER'  MDIA  'CdI2'          INCO 'I2'
      ZONE  $fron 'OPER'  MDIA  'CdTI2'         INCO 'TN' 'I2'
;
* inconnue I3
 RV = EQEX  RV            'OPTI'  'EF'  'IMPL'
      ZONE  $MT   'OPER'  LAPN   1.             INCO 'I3'
      ZONE  $MT   'OPER'  MDIA  'CI3'           INCO 'I3'
      ZONE  $MT   'OPER'  MDIA  'CTI3'          INCO 'TN' 'I3'
      ZONE  $fron 'OPER'  MDIA  'CdI3'          INCO 'I3'
      ZONE  $fron 'OPER'  MDIA  'CdTI3'         INCO 'TN' 'I3'
;
* inconnue I4
 RV = EQEX  RV            'OPTI'  'EF'  'IMPL'
      ZONE  $MT   'OPER'  LAPN   1.             INCO 'I4'
      ZONE  $MT   'OPER'  MDIA  'CI4'           INCO 'I4'
      ZONE  $MT   'OPER'  MDIA  'CTI4'          INCO 'TN' 'I4'
      ZONE  $fron 'OPER'  MDIA  'CdI4'          INCO 'I4'
      ZONE  $fron 'OPER'  MDIA  'CdTI4'         INCO 'TN' 'I4'
;
 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.'I1'   = kcht $mt scal sommet 1.               ;
rv.INCO.'I2'   = kcht $mt scal sommet 1.               ;
rv.INCO.'I3'   = kcht $mt scal sommet 1.               ;
rv.INCO.'I4'   = 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      ;


*----------------------------------
*   Algorithme
*----------------------------------

 repeter bloc2   ;

 nbiter = &bloc2   ;

* calcul des facteurs de ponderation pour les 3 gas

 T = rv.INCO.'TN' ;
 4T3 = 4.0 * (T ** 3) ;

 TK = kops Tref * T ;
 w1 = CWSGG  TK a0 a1 a2 a3 ;
 w2 = CWSGG  TK b0 b1 b2 b3 ;
 w3 = CWSGG  TK c0 c1 c2 c3 ;
   cc1 = kops w1 + w2 ;
   cc2 = kops w3 + cc1    ;
 w4 = kops 1. - cc2  ;

 w4T3_1 = kops w1 * 4T3 ;
 w4T3_2 = kops w2 * 4T3 ;
 w4T3_3 = kops w3 * 4T3 ;
 w4T3_4 = kops w4 * 4T3 ;

*----- actualisation des coefficients -----------------------------

* gaz 1
rv.INCO.'CI1'   = ATT1;
rv.INCO.'CTI1'  = kops (-1.*ATT1) * w4T3_1;
rv.INCO.'CdI1'  = kcht $fron scal sommet ATE1 ;
rv.INCO.'CdTI1' = kcht $fron scal sommet (kops (-1.*ATE1) * w4T3_1) ;

rv.INCO.'CI1T'  =  3Ntau1 * rv.INCO.'CI1' ;
rv.INCO.'CT1'   = kops 3Ntau1 * rv.INCO.'CTI1' ;
            CF1 = kops 3Ntau1 * rv.INCO.'CdI1' ;
rv.INCO.'CdI1T' = kcht $l13 scal sommet CF1  ;
            CF2 = kops 3Ntau1 * rv.INCO.'CdTI1';
rv.INCO.'CdT1'  = kcht $l13 scal sommet CF2 ;

* gaz 2
rv.INCO.'CI2'   = ATT2 ;
rv.INCO.'CTI2'  = kops (-1.*ATT2) * w4T3_2;
rv.INCO.'CdI2'  = kcht $fron scal sommet ATE2 ;
rv.INCO.'CdTI2' = kcht $fron scal sommet (kops (-1.*ATE2) * w4T3_2) ;

rv.INCO.'CI2T'  =  3Ntau2 * rv.INCO.'CI2' ;
rv.INCO.'CT2'   = kops 3Ntau2 * rv.INCO.'CTI2' ;
            CF1 = kops 3Ntau2 * rv.INCO.'CdI2' ;
rv.INCO.'CdI2T' = kcht $l13 scal sommet CF1  ;
            CF2 = kops 3Ntau2 * rv.INCO.'CdTI2';
rv.INCO.'CdT2'  = kcht $l13 scal sommet CF2 ;

* gaz 3
rv.INCO.'CI3'   = ATT3 ;
rv.INCO.'CTI3'  = kops (-1.*ATT3) * w4T3_3;
rv.INCO.'CdI3'  = kcht $fron scal sommet ATE3 ;
rv.INCO.'CdTI3' = kcht $fron scal sommet (kops (-1.*ATE3) * w4T3_3) ;

rv.INCO.'CI3T'  =  3Ntau3 * rv.INCO.'CI3' ;
rv.INCO.'CT3'   = kops 3Ntau3 * rv.INCO.'CTI3' ;
            CF1 = kops 3Ntau3 * rv.INCO.'CdI3' ;
rv.INCO.'CdI3T' = kcht $l13 scal sommet CF1  ;
            CF2 = kops 3Ntau3 * rv.INCO.'CdTI3';
rv.INCO.'CdT3'  = kcht $l13 scal sommet CF2 ;

* gaz 4
rv.INCO.'CI4'   = ATT4 ;
rv.INCO.'CTI4'  = kops (-1.*ATT4) * w4T3_4;
rv.INCO.'CdI4'  = kcht $fron scal sommet ATE4 ;
rv.INCO.'CdTI4' = kcht $fron scal sommet (kops (-1.*ATE4) * w4T3_4) ;

rv.INCO.'CI4T'  =  3Ntau4 * rv.INCO.'CI4' ;
rv.INCO.'CT4'   = kops 3Ntau4 * rv.INCO.'CTI4' ;
            CF1 = kops 3Ntau4 * rv.INCO.'CdI4' ;
rv.INCO.'CdI4T' = kcht $l13 scal sommet CF1  ;
            CF2 = kops 3Ntau4 * rv.INCO.'CdTI4';
rv.INCO.'CdT4'  = kcht $l13 scal sommet CF2 ;

* --- fin actualisation des coefficients --------------------------

UN = copie rv.INCO.'UN' ;
TN = copie rv.INCO.'TN' ;
I1 = copie rv.INCO.'I1' ;
I2 = copie rv.INCO.'I2' ;
I3 = copie rv.INCO.'I3' ;
I4 = copie rv.INCO.'I4' ;

 exec rv ;

XT  = (maxi (abs (kops ((rv.INCO.'TN') - TN) / TN))) ;
XI1 = (maxi (abs (kops ((rv.INCO.'I1') - I1) / I1))) ;
XI2 = (maxi (abs (kops ((rv.INCO.'I2') - I2) / I2))) ;
XI3 = (maxi (abs (kops ((rv.INCO.'I3') - I3) / I3))) ;
XI4 = (maxi (abs (kops ((rv.INCO.'I4') - I4) / I4))) ;

MESS ' ERREUR RELATIVE TN ' XT 'I1 ' XI1 'I2' XI2 'I3' XI3 'I4' XI4;

Tmax = maxi (rv.INCO.'TN') ;
Tmin = mini (rv.INCO.'TN') ;

*MESS ' EXTREMUM TN      MAX' Tmax 'MIN ' Tmin ;
* mess ' iteration '  nbiter ;

* --- test de convergence -------------------------------------------

 si ((XT < EPS) et (XI1 < EPS)) ;
 quit bloc2 ;
 sinon ;
   si (nbiter > maxiter) ;
    mess ' pas de convergence' ; erre 5 ;
   finsi ;
 finsi      ;

fin bloc2 ;

IN = rv.INCO.'I1' ;

*----------------------------------
*   Calcul des Nusselts
*----------------------------------

si post ;
*
* 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

I1_l4  = redu rv.INCO.'I1' l4  ;
I2_l4  = redu rv.INCO.'I2' l4  ;
I3_l4  = redu rv.INCO.'I3' l4  ;
I4_l4  = redu rv.INCO.'I4' l4  ;

T_l4  = redu rv.INCO.'TN' l4  ;
4T4 = 4.0 * (T_l4 ** 4) ;

dI1dn =  ATE1 * (I1_l4 - (w1 * 4T4)) ;
dI2dn =  ATE2 * (I2_l4 - (w2 * 4T4)) ;
dI3dn =  ATE3 * (I3_l4 - (w3 * 4T4)) ;
dI4dn =  ATE4 * (I4_l4 - (w4 * 4T4)) ;

Nu_I1 = KCHT $l4   'SCAL' 'SOMMET' (kops dI1dn * 3Ntau1) ;
Nu_I2 = KCHT $l4   'SCAL' 'SOMMET' (kops dI2dn * 3Ntau2) ;
Nu_I3 = KCHT $l4   'SCAL' 'SOMMET' (kops dI3dn * 3Ntau3) ;
Nu_I4 = KCHT $l4   'SCAL' 'SOMMET' (kops dI4dn * 3Ntau4) ;

Nu_I = kops (kops Nu_I1 + Nu_I2) + (kops Nu_i3 + Nu_I4) ;
*Nu_I =  kops Nu_I1 + Nu_I2 ;

* 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

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 ;

finsi ;

*----------------------------------
*   Graphique
*----------------------------------
si graph ;

  opti titre 'CONDUCTION - RAYONNEMENT  Ra = 1.e5 Nr = 1. ';
  f_homo = 0.01  ;

* 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 ;
opti isov ligne ;
trac TX mt (cont mt) LTRA ;

* intensites radiatives

trac rv.INCO.'I1' 10 mt (cont mt) ;
trac rv.INCO.'I2' 10 mt (cont mt) ;
trac rv.INCO.'I3' 10 mt (cont mt) ;
trac rv.INCO.'I4' 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 ' ;
TAB1. 2 = 'MARQ TRIA REGU ' ;
TAB1.'TITRE'   = TABLE ;
TAB1.'TITRE'. 1 = MOT ' PAROI_INFERIEURE' ;
TAB1.'TITRE'. 2 = MOT ' PAROI_SUPERIEURE' ;

dess (T1 et T3) MIMA  LEGE TAB1
    'TITR' 'TEMPERATURE'
    'TITX' ' x ' 'TITY' ' T ' ;

*opti donn 5 ;

* evolution des Nusselt sur la paroi chaude

Nu_Tl = evol 'CHPO' Nu_T SCAL l4  ;
Nu_Il = evol 'CHPO' Nu_I SCAL l4 ;

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) LEGE TAB3
    'TITR' 'Nombre de Nusselt local sur la paroi chaude '
    'TITX' ' y ' 'TITY' ' Nusselt ' ;

opti donn 5 ;

finsi ;

fin ;