Télécharger la plaquette Cast3M


fr en

Télécharger wsgg.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  wsgg.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4.  
  5. ******************  CONVECTION-RAYONNEMENT + VAPEUR D'EAU  ********** *
  6. *                                                                     *
  7. * Couplage convection naturelle laminaire/rayonnement milieu absorbant*
  8. *                                                                     *
  9. * Convection naturelle horizontale dans une cavite carree contenant   *
  10. * de la vapeur d'eau a la pression de 1 atm.                          *
  11. * La temperature des parois varie entre 1000K et 2000K.               *
  12. *                                                                     *
  13. * Par rapport au test de base cvry-2D-1.dgibi le milieu est un gas    *
  14. * reel modelise par le modele WSGG (Weighted Sum of Gray Gases)       *
  15. * dans lequel le gas est decompose en la somme de 3 gaz gris affectés *
  16. * chacun d'un facteur de ponderation dependant polynomialement de     *
  17. * la temperature et d'un coefficient d'absorption constant.           *
  18. *                                                                     *
  19. * Données: test de Tan avec Ra = 1.e4 et Nr = 1.                      *
  20. * Validation numerique par comparaison avec Trio-EF                   *
  21. *                                                                     *
  22. * Formulation EF                                                      *
  23. * Algorithme: implicite en régime permanent avec 6 inconnues          *
  24. *             convection traitée par pénalisation                     *
  25. *             rayonnement traité par la methode des harmoniques       *
  26. *             sphériques (P1) + modele WSGG a 3 gaz                   *
  27. *             non linearite des proprietes radiatives                 *
  28. *                                                                     *
  29. * opérateurs utilisés : DUDW, KONV, LAPN, MDIA                        *
  30. *                                                                     *
  31. * References:                                                         *
  32. * - pour le couplage convection rayonnement et l'approximation P1     *
  33. *    test cvry-2D-1.dgibi                                             *
  34. *    Z.Tan J.R.Howell, I.J.H.M.T Vol.34 1991                          *
  35. *    rapport DMT/96.030 et DMT/94.505                                 *
  36. * - pour les proprietes radiatives de la vapeur d'eau:                *
  37. *    T.H.Smith J.H.T ASME Vol.104 1982                                *
  38. *    rapport DMT/94.505                                               *
  39. *                                                                     *
  40. ***********************************************************************
  41.  
  42. complet = faux ;
  43. post = vrai ; graph = faux ;
  44.  
  45. option echo 0 dime 2 elem qua4 ;
  46.  
  47. *----------------------------------
  48. *   Parametres physiques
  49. *----------------------------------
  50.  
  51. Ra  = 1.e4 ; Nr = 1.0 ;
  52.  
  53. om = 0. ;
  54. Tref = 1000. ; L = 1. ;
  55. phi0 = 1. ; th  = 2.  ; tc = 1. ;
  56. Pr = 0.71 ; Pe = 1. ; Re = Pe/Pr ;
  57. a = 1.  ;
  58. emis = 1.  ;
  59.  
  60. * valeurs de reference pour les Nusselt - version du 8 jan. 98
  61. * test de Tan
  62.   Nu_Tref = 2.15 ; Nu_Iref = 13.1 ;
  63.  
  64. *----------------------------------
  65. *   Parametres numeriques
  66. *----------------------------------
  67.  
  68. OMEGA = 0.6   ;
  69. EPSS  = 1.E-10;
  70. maxiter = 60 ;
  71. si complet; EPS = 1.E-4 ; sinon ; EPS = 1.E-2 ; finsi;
  72.  
  73. * ---------------------------------
  74. * Re    : nombre de Reynolds
  75. * Ra    : nombre de Rayleigh
  76. * Pr    : nombre de Prandtl
  77. * Pe    : nombre de Peclet = Pr * Re
  78. * Nr    : nombre caracteristique du rayonnement
  79.  
  80. * Tref  : temperature de reference en degre K
  81. * L     : longueur caracteristique
  82. *         (pour le calcul des proprietes radiatives)
  83. * t0    : temperature de reference pour le pb adimensionne
  84. * th    : hot temperature 
  85. * tc    : cold temperature
  86. * emis  : emissivite des parois
  87. * a     : rapport d'allongement (largeur/hauteur)
  88. * ki    : coefficient d'absorption du gas i
  89. * taui  : epaisseur optique (ki * L)
  90. * om    : albedo  ( coeff.diffus /(coeff.diffus + coeff.absorb) )
  91. * phi0  : paramètre de forme
  92. *
  93. * ---------------------------------
  94. * OMEGA   : facteur de relaxation
  95. * EPS     : precision sur les itérations internes
  96. * EPSS    : facteur de penalisation
  97. * maxiter : nombre maximum d'iterations
  98.  
  99. *----------------------------------
  100. *   Proprietes radiatives
  101. *----------------------------------
  102. *     H2O Ptot=1atm
  103. * decomposition de la vapeur d'eau en 3 gas gris
  104. * definition des coefficient d'absorption ki et des coefficients
  105. * pour le calcul des facteurs de ponderation wi
  106. * taui designe l'epaisseur optique correspondante
  107.  
  108.   k1 = 0.4496 ; k2 = 7.113 ; k3 = 119.7 ; k4 = 1.e-5 ;
  109.   tau1 = k1*L ; tau2 = k2*L ; tau3 = k3*L ; tau4 = k4*L ;
  110.  
  111.  a0= 6.324e-1 ; a1=-8.358E-4 ; a2= 6.135e-7 ; a3= -13.03e-11 ;
  112.  b0=-2.016e-2 ; b1= 7.145E-4 ; b2=-5.212e-7 ; b3=  9.868e-11 ;
  113.  c0= 3.500e-1 ; c1=-5.040E-4 ; c2= 2.425e-7 ; c3= -3.888e-11 ;
  114.  
  115. * utilisation de la procedure CWSGG pour le calcul des wi
  116. *  w1 = a0+a1*T+a2*T*T+a3*T*T*T ;
  117. *  w2 = b0+b1*T+b2*T*T+b3*T*T*T ;
  118. *  w3 = c0+c1*T+c2*T*T+c3*T*T*T ;
  119. *  w4= 1.0 -w1-w2-w3 ;
  120.  
  121. *---------------------------------------------------------------
  122.  'DEBPROC' CWSGG T*'CHPOINT'
  123.            a0*'FLOTTANT' a1*'FLOTTANT' a2*'FLOTTANT' a3*'FLOTTANT' ;
  124. * calcul de la loi polynomiale donnant les facteurs de ponderation
  125. * en fonction de la temperature en degre K
  126. * w  = a0+a1*T+a2*T*T+a3*T*T*T ;
  127.   b2 = 'KOPS' a2 + ('KOPS' T * a3) ;
  128.   b1 = 'KOPS' a1 + ('KOPS' T * b2) ;
  129.   w  = 'KOPS' a0 + ('KOPS' T * b1) ;
  130. 'FINPROC' w    ;
  131. *---------------------------------------------------------------
  132.  
  133.  
  134.  
  135. *----------------------------------
  136. *   Maillage
  137. *----------------------------------
  138.  
  139. p1= 0.  0.  ; 
  140. p2= 0.5 0.  ;
  141. p3= 0.5 0.5 ; 
  142. p4= 0.  0.5 ;
  143. p5= 1.  0.  ; 
  144. p6= 0.  1.  ;
  145.  
  146. * maill. fin
  147. *di = 0.005 ; df = 0.055 ;
  148. *vi = 0.011 ; vf = df    ;
  149.  
  150. si complet ;
  151. * maill. moyen
  152.  di = 0.005 ; df = 0.1 ;
  153.  vi = 0.02  ; vf = df  ;
  154. sinon ;
  155. * maill. grossier
  156.  di = 0.01 ; df = 0.2 ;
  157.  vi = 0.1  ; vf = df  ;
  158. finsi;
  159.  
  160. * ligne horiz.
  161.  c1a=  droit p1 p2 dini di dfin df ;
  162.  c1b = droit p2 p5 dini df dfin di ;
  163.  
  164. * ligne vert.
  165.  c4a=  droit p6 p4 dini vi dfin vf ;
  166.  c4b = droit p4 p1 dini vf dfin vi ;
  167.  
  168. l1 = c1a et c1b   ;
  169. l4 = c4a et c4b   ;
  170. l3inv= l1 plus (0. 1.) ;  l3 = inve ( l3inv ) ;
  171. l2inv= l4 plus (1. 0.) ;  l2 = inve ( l2inv ) ;
  172. fron = l1 et l2 et l3 et l4 ;
  173. l13  = l1 et l3 ;
  174. elim 1e-5 (l1 et l2 et l3 et l4 ) ;
  175.  
  176. mt = l1 l2 l3 l4 dalle plan ;
  177. elim 1.e-5  mt ;
  178. tass mt ;
  179. orient mt ;
  180.  si graph ; trac mt ; finsi ;
  181.  
  182. *----------------------------------
  183. *   Tables domaine
  184. *----------------------------------
  185.  
  186. $mt   = doma mt                      impr ;
  187. $l13  = doma l13  'INCL' $mt  0.0001 impr ;
  188. $fron = doma fron 'INCL' $mt  0.0001 impr ;
  189.  
  190.  
  191. *----------------------------------
  192. *   Definition de coefficients divers
  193. *----------------------------------
  194.  
  195. t0  = (th+tc)/2. ;
  196. Gr = Ra/Pr ; gg = Gr/(Re*Re) ;
  197.  
  198. ALFV = 1./Re ;
  199. ALFM = 1   ;
  200.  
  201. E = emis/(4.-2.*emis) ;
  202. 3Nr = 3. * Nr ;
  203.  
  204. 3Ntau1 = (-1.0)/(3Nr*tau1) ;
  205. 3Ntau2 = (-1.0)/(3Nr*tau2) ;
  206. 3Ntau3 = (-1.0)/(3Nr*tau3) ;
  207. 3Ntau4 = (-1.0)/(3Nr*tau4) ;
  208.  
  209. ATE1 = (3.*E)*(a*tau1)   ;
  210. ATT1 = (3.*((a*tau1)**2))*(1.-om) ;
  211.  
  212. ATE2 = (3.*E)*(a*tau2)   ;
  213. ATT2 = (3.*((a*tau2)**2))*(1.-om) ;
  214.  
  215. ATE3 = (3.*E)*(a*tau3)   ;
  216. ATT3 = (3.*((a*tau3)**2))*(1.-om) ;
  217.  
  218. ATE4 = (3.*E)*(a*tau4)   ;
  219. ATT4 = (3.*((a*tau4)**2))*(1.-om) ;
  220.  
  221.  
  222. *----------------------------------
  223. *   Initialisation du champ de temperature
  224. *----------------------------------
  225.  
  226.  RT = EQEX  $MT            'OPTI'  'EF'  'IMPL'
  227.       ZONE  $MT   'OPER'  LAPN  1.     INCO 'TT'  
  228. ;
  229.  RT = EQEX  RT
  230.       CLIM
  231.       'TT'  TIMP   l4    th   
  232.       'TT'  TIMP   l2    tc 
  233. ;
  234. rt.INCO       = TABLE 'INCO'             ;
  235. rt.INCO.'TT'  = kcht $mt scal sommet t0  ;
  236. rt.'NITER' = 1    ;
  237. rt.'OMEGA' = 1.   ;
  238.  
  239. exec rt ;
  240.  
  241.  
  242. *----------------------------------
  243. *   Definition du systeme lineaire
  244. *----------------------------------
  245.  
  246. *  bilan masse et qdm (vitesse UN)
  247. *  bilan energie      (temperature)
  248. *  bilan sur les 4 intensites radiatives de la decomposition
  249. *      (I1, I2, I3, I4)
  250.  
  251.  RV = EQEX  $MT           'OPTI'  'EF'  'IMPL'
  252.       ZONE  $MT   'OPER'  DUDW  EPSS            INCO 'UN'
  253.       ZONE  $MT   'OPER'  KONV  1. 'UN'  ALFV   INCO 'UN'
  254.       ZONE  $MT   'OPER'  LAPN  ALFV            INCO 'UN'
  255.       ZONE  $MT   'OPER'  MDIA  'GB'            INCO 'TN'  'UN'
  256. ;
  257. *  inconnue T
  258.  RV = EQEX  RV            'OPTI'  'EF'  'IMPL'
  259.       ZONE  $MT   'OPER ' KONV  Pe  'UN'  ALFM  INCO 'TN'
  260.       ZONE  $MT   'OPER'  LAPN   ALFM           INCO 'TN'
  261. ;
  262.  RV = EQEX  RV            'OPTI'  'EF'  'IMPL'
  263.       ZONE  $MT   'OPER'  MDIA  'CT1'           INCO 'TN'
  264.       ZONE  $MT   'OPER'  MDIA  'CI1T'          INCO 'I1' 'TN'
  265.       ZONE  $l13  'OPER'  MDIA  'CdT1'          INCO 'TN'
  266.       ZONE  $l13  'OPER'  MDIA  'CdI1T'         INCO 'I1' 'TN'
  267. ;
  268.  RV = EQEX  RV            'OPTI'  'EF'  'IMPL'
  269.       ZONE  $MT   'OPER'  MDIA  'CT2'           INCO 'TN'
  270.       ZONE  $MT   'OPER'  MDIA  'CI2T'          INCO 'I2' 'TN'
  271.       ZONE  $l13  'OPER'  MDIA  'CdT2'          INCO 'TN'
  272.       ZONE  $l13  'OPER'  MDIA  'CdI2T'         INCO 'I2' 'TN'
  273. ;
  274.  RV = EQEX  RV            'OPTI'  'EF'  'IMPL'
  275.       ZONE  $MT   'OPER'  MDIA  'CT3'           INCO 'TN'
  276.       ZONE  $MT   'OPER'  MDIA  'CI3T'          INCO 'I3' 'TN'
  277.       ZONE  $l13  'OPER'  MDIA  'CdT3'          INCO 'TN'
  278.       ZONE  $l13  'OPER'  MDIA  'CdI3T'         INCO 'I3' 'TN'
  279. ;
  280.  RV = EQEX  RV            'OPTI'  'EF'  'IMPL'
  281.       ZONE  $MT   'OPER'  MDIA  'CT4'           INCO 'TN'
  282.       ZONE  $MT   'OPER'  MDIA  'CI4T'          INCO 'I4' 'TN'
  283.       ZONE  $l13  'OPER'  MDIA  'CdT4'          INCO 'TN'
  284.       ZONE  $l13  'OPER'  MDIA  'CdI4T'         INCO 'I4' 'TN'
  285. ;
  286. * inconnue I1
  287.  RV = EQEX  RV            'OPTI'  'EF'  'IMPL'
  288.       ZONE  $MT   'OPER'  LAPN   1.             INCO 'I1'
  289.       ZONE  $MT   'OPER'  MDIA  'CI1'           INCO 'I1'
  290.       ZONE  $MT   'OPER'  MDIA  'CTI1'          INCO 'TN' 'I1'
  291.       ZONE  $fron 'OPER'  MDIA  'CdI1'          INCO 'I1'
  292.       ZONE  $fron 'OPER'  MDIA  'CdTI1'         INCO 'TN' 'I1'
  293. ;
  294. * inconnue I2
  295.  RV = EQEX  RV            'OPTI'  'EF'  'IMPL'
  296.       ZONE  $MT   'OPER'  LAPN   1.             INCO 'I2'
  297.       ZONE  $MT   'OPER'  MDIA  'CI2'           INCO 'I2'
  298.       ZONE  $MT   'OPER'  MDIA  'CTI2'          INCO 'TN' 'I2'
  299.       ZONE  $fron 'OPER'  MDIA  'CdI2'          INCO 'I2'
  300.       ZONE  $fron 'OPER'  MDIA  'CdTI2'         INCO 'TN' 'I2'
  301. ;
  302. * inconnue I3
  303.  RV = EQEX  RV            'OPTI'  'EF'  'IMPL'
  304.       ZONE  $MT   'OPER'  LAPN   1.             INCO 'I3'
  305.       ZONE  $MT   'OPER'  MDIA  'CI3'           INCO 'I3'
  306.       ZONE  $MT   'OPER'  MDIA  'CTI3'          INCO 'TN' 'I3'
  307.       ZONE  $fron 'OPER'  MDIA  'CdI3'          INCO 'I3'
  308.       ZONE  $fron 'OPER'  MDIA  'CdTI3'         INCO 'TN' 'I3'
  309. ;
  310. * inconnue I4
  311.  RV = EQEX  RV            'OPTI'  'EF'  'IMPL'
  312.       ZONE  $MT   'OPER'  LAPN   1.             INCO 'I4'
  313.       ZONE  $MT   'OPER'  MDIA  'CI4'           INCO 'I4'
  314.       ZONE  $MT   'OPER'  MDIA  'CTI4'          INCO 'TN' 'I4'
  315.       ZONE  $fron 'OPER'  MDIA  'CdI4'          INCO 'I4'
  316.       ZONE  $fron 'OPER'  MDIA  'CdTI4'         INCO 'TN' 'I4'
  317. ;
  318.  RV = EQEX RV
  319.       CLIM
  320.       'UN'  UIMP   fron   0.0 
  321.       'UN'  VIMP   fron   0.0 
  322.       'TN'  TIMP   l4     th
  323.       'TN'  TIMP   l2     tc
  324. ;
  325. rv.INCO        = TABLE 'INCO'                          ;
  326. rv.INCO.'TN'   = kcht $mt scal sommet rt.INCO.'TT'     ;
  327. rv.INCO.'I1'   = kcht $mt scal sommet 1.               ;
  328. rv.INCO.'I2'   = kcht $mt scal sommet 1.               ;
  329. rv.INCO.'I3'   = kcht $mt scal sommet 1.               ;
  330. rv.INCO.'I4'   = kcht $mt scal sommet 1.               ;
  331. rv.INCO.'GB'   = kcht $mt vect sommet (0. (-1.*phi0*gg)) ;
  332. rv.INCO.'UN'   = kcht $mt vect sommet (1.e-5 1.e-5)    ;
  333.  
  334. rv.'NITER'  = 1      ;
  335. rv.'OMEGA'  = OMEGA  ;
  336. *rv.'IMPR'   = 1      ;
  337. rv.'ITMA'   = 1      ;
  338.  
  339.  
  340. *----------------------------------
  341. *   Algorithme
  342. *----------------------------------
  343.  
  344.  repeter bloc2   ;
  345.  
  346.  nbiter = &bloc2   ;
  347.  
  348. * calcul des facteurs de ponderation pour les 3 gas
  349.  
  350.  T = rv.INCO.'TN' ;
  351.  4T3 = 4.0 * (T ** 3) ;
  352.  
  353.  TK = kops Tref * T ;
  354.  w1 = CWSGG  TK a0 a1 a2 a3 ;
  355.  w2 = CWSGG  TK b0 b1 b2 b3 ;
  356.  w3 = CWSGG  TK c0 c1 c2 c3 ;
  357.    cc1 = kops w1 + w2 ;
  358.    cc2 = kops w3 + cc1    ;
  359.  w4 = kops 1. - cc2  ;
  360.  
  361.  w4T3_1 = kops w1 * 4T3 ;
  362.  w4T3_2 = kops w2 * 4T3 ;
  363.  w4T3_3 = kops w3 * 4T3 ;
  364.  w4T3_4 = kops w4 * 4T3 ;
  365.  
  366. *----- actualisation des coefficients -----------------------------
  367.  
  368. * gaz 1
  369. rv.INCO.'CI1'   = ATT1;
  370. rv.INCO.'CTI1'  = kops (-1.*ATT1) * w4T3_1;
  371. rv.INCO.'CdI1'  = kcht $fron scal sommet ATE1 ;
  372. rv.INCO.'CdTI1' = kcht $fron scal sommet (kops (-1.*ATE1) * w4T3_1) ;
  373.  
  374. rv.INCO.'CI1T'  =  3Ntau1 * rv.INCO.'CI1' ;
  375. rv.INCO.'CT1'   = kops 3Ntau1 * rv.INCO.'CTI1' ;
  376.             CF1 = kops 3Ntau1 * rv.INCO.'CdI1' ;
  377. rv.INCO.'CdI1T' = kcht $l13 scal sommet CF1  ;
  378.             CF2 = kops 3Ntau1 * rv.INCO.'CdTI1';
  379. rv.INCO.'CdT1'  = kcht $l13 scal sommet CF2 ;
  380.  
  381. * gaz 2
  382. rv.INCO.'CI2'   = ATT2 ;
  383. rv.INCO.'CTI2'  = kops (-1.*ATT2) * w4T3_2;
  384. rv.INCO.'CdI2'  = kcht $fron scal sommet ATE2 ;
  385. rv.INCO.'CdTI2' = kcht $fron scal sommet (kops (-1.*ATE2) * w4T3_2) ;
  386.  
  387. rv.INCO.'CI2T'  =  3Ntau2 * rv.INCO.'CI2' ;
  388. rv.INCO.'CT2'   = kops 3Ntau2 * rv.INCO.'CTI2' ;
  389.             CF1 = kops 3Ntau2 * rv.INCO.'CdI2' ;
  390. rv.INCO.'CdI2T' = kcht $l13 scal sommet CF1  ;
  391.             CF2 = kops 3Ntau2 * rv.INCO.'CdTI2';
  392. rv.INCO.'CdT2'  = kcht $l13 scal sommet CF2 ;
  393.  
  394. * gaz 3
  395. rv.INCO.'CI3'   = ATT3 ;
  396. rv.INCO.'CTI3'  = kops (-1.*ATT3) * w4T3_3;
  397. rv.INCO.'CdI3'  = kcht $fron scal sommet ATE3 ;
  398. rv.INCO.'CdTI3' = kcht $fron scal sommet (kops (-1.*ATE3) * w4T3_3) ;
  399.  
  400. rv.INCO.'CI3T'  =  3Ntau3 * rv.INCO.'CI3' ;
  401. rv.INCO.'CT3'   = kops 3Ntau3 * rv.INCO.'CTI3' ;
  402.             CF1 = kops 3Ntau3 * rv.INCO.'CdI3' ;
  403. rv.INCO.'CdI3T' = kcht $l13 scal sommet CF1  ;
  404.             CF2 = kops 3Ntau3 * rv.INCO.'CdTI3';
  405. rv.INCO.'CdT3'  = kcht $l13 scal sommet CF2 ;
  406.  
  407. * gaz 4
  408. rv.INCO.'CI4'   = ATT4 ;
  409. rv.INCO.'CTI4'  = kops (-1.*ATT4) * w4T3_4;
  410. rv.INCO.'CdI4'  = kcht $fron scal sommet ATE4 ;
  411. rv.INCO.'CdTI4' = kcht $fron scal sommet (kops (-1.*ATE4) * w4T3_4) ;
  412.  
  413. rv.INCO.'CI4T'  =  3Ntau4 * rv.INCO.'CI4' ;
  414. rv.INCO.'CT4'   = kops 3Ntau4 * rv.INCO.'CTI4' ;
  415.             CF1 = kops 3Ntau4 * rv.INCO.'CdI4' ;
  416. rv.INCO.'CdI4T' = kcht $l13 scal sommet CF1  ;
  417.             CF2 = kops 3Ntau4 * rv.INCO.'CdTI4';
  418. rv.INCO.'CdT4'  = kcht $l13 scal sommet CF2 ;
  419.  
  420. * --- fin actualisation des coefficients --------------------------
  421.  
  422. UN = copie rv.INCO.'UN' ;
  423. TN = copie rv.INCO.'TN' ;
  424. I1 = copie rv.INCO.'I1' ;
  425. I2 = copie rv.INCO.'I2' ;
  426. I3 = copie rv.INCO.'I3' ;
  427. I4 = copie rv.INCO.'I4' ;
  428.  
  429.  exec rv ;
  430.  
  431. XT  = (maxi (abs (kops ((rv.INCO.'TN') - TN) / TN))) ;
  432. XI1 = (maxi (abs (kops ((rv.INCO.'I1') - I1) / I1))) ;
  433. XI2 = (maxi (abs (kops ((rv.INCO.'I2') - I2) / I2))) ;
  434. XI3 = (maxi (abs (kops ((rv.INCO.'I3') - I3) / I3))) ;
  435. XI4 = (maxi (abs (kops ((rv.INCO.'I4') - I4) / I4))) ;
  436.  
  437. MESS ' ERREUR RELATIVE TN ' XT 'I1 ' XI1 'I2' XI2 'I3' XI3 'I4' XI4;
  438.  
  439. Tmax = maxi (rv.INCO.'TN') ;
  440. Tmin = mini (rv.INCO.'TN') ;
  441.  
  442. *MESS ' EXTREMUM TN      MAX' Tmax 'MIN ' Tmin ;
  443. * mess ' iteration '  nbiter ;
  444.  
  445. * --- test de convergence -------------------------------------------
  446.  
  447.  si ((XT < EPS) et (XI1 < EPS)) ;
  448.  quit bloc2 ;
  449.  sinon ;
  450.    si (nbiter > maxiter) ;
  451.     mess ' pas de convergence' ; erre 5 ;
  452.    finsi ;
  453.  finsi      ;
  454.  
  455. fin bloc2 ;
  456.  
  457. IN = rv.INCO.'I1' ;
  458.  
  459. *----------------------------------
  460. *   Calcul des Nusselts
  461. *----------------------------------
  462.  
  463. si post ;
  464. *
  465. * calcul des Nusselt dus a la conduction et au rayonnement
  466. * sur la paroi chaude (l4)
  467.  
  468. $l4 = doma l4 ;
  469.  
  470. * Nusselt du a la conduction
  471.  
  472. gradt = KOPS rv.'INCO'.'TN' 'GRAD' $mt ;
  473. dTdXe = KCHT $mt 'SCAL' 'CENTRE' (exco 'UX' gradt) ;
  474. dTdX = ELNO $MT (kops (-1.) * dTdXe) ;
  475.  
  476. Nu_T  = KCHT $l4   'SCAL' 'SOMMET' (kops phi0 * dTdX) ;
  477.  
  478. * Nusselt du au rayonnement
  479.  
  480. I1_l4  = redu rv.INCO.'I1' l4  ;
  481. I2_l4  = redu rv.INCO.'I2' l4  ;
  482. I3_l4  = redu rv.INCO.'I3' l4  ;
  483. I4_l4  = redu rv.INCO.'I4' l4  ;
  484.  
  485. T_l4  = redu rv.INCO.'TN' l4  ;
  486. 4T4 = 4.0 * (T_l4 ** 4) ;
  487.  
  488. dI1dn =  ATE1 * (I1_l4 - (w1 * 4T4)) ;
  489. dI2dn =  ATE2 * (I2_l4 - (w2 * 4T4)) ;
  490. dI3dn =  ATE3 * (I3_l4 - (w3 * 4T4)) ;
  491. dI4dn =  ATE4 * (I4_l4 - (w4 * 4T4)) ;
  492.  
  493. Nu_I1 = KCHT $l4   'SCAL' 'SOMMET' (kops dI1dn * 3Ntau1) ;
  494. Nu_I2 = KCHT $l4   'SCAL' 'SOMMET' (kops dI2dn * 3Ntau2) ;
  495. Nu_I3 = KCHT $l4   'SCAL' 'SOMMET' (kops dI3dn * 3Ntau3) ;
  496. Nu_I4 = KCHT $l4   'SCAL' 'SOMMET' (kops dI4dn * 3Ntau4) ;
  497.  
  498. Nu_I = kops (kops Nu_I1 + Nu_I2) + (kops Nu_i3 + Nu_I4) ;
  499. *Nu_I =  kops Nu_I1 + Nu_I2 ;
  500.  
  501. * Nusselt global pour la paroi
  502.  
  503. dy = 'DOMA' $l4   'VOLUME' ;
  504. dyt = SOMT dy ;
  505. *message 'LONGUEUR DE LA PAROI = ' dyt ;
  506.  
  507. Nu_Te = noel $l4 Nu_T  ;
  508. Nu_Tdy = kops Nu_Te  '*' dy ;
  509. Nu_Tdyt = somt Nu_Tdy ;
  510. Nu_Tm = Nu_Tdyt / dyt ;
  511. message 'Nusselt de conduction moyen: ' Nu_Tm ;
  512.  
  513. Nu_Ie = noel $l4 Nu_I  ;
  514. Nu_Idy = kops Nu_Ie  '*' dy ;
  515. Nu_Idyt = somt Nu_Idy ;
  516. Nu_Im = Nu_Idyt / dyt ;
  517. message 'Nusselt de rayonnement moyen: ' Nu_Im ;
  518.  
  519. * comparaison aux valeurs de reference
  520.  
  521. resT = abs((Nu_Tm-Nu_Tref/Nu_Tref));
  522. resI = abs((Nu_Im-Nu_Iref/Nu_Iref));
  524. si ((resT < 1e-2) et (resI < 1e-2 )) ; erre 0 ;
  525. sinon ; erre 5 ; finsi ;
  526.  
  527. finsi ;
  528.  
  529. *----------------------------------
  530. *   Graphique
  531. *----------------------------------
  532. si graph ;
  533.  
  534.   titre 'CONDUCTION - RAYONNEMENT  Ra = 1.e5 Nr = 1. ';
  535.   f_homo = 0.01  ;
  536.  
  537. * température relative comprise entre -0.5 et 0.5
  538.  
  539. TY = kops rv.INCO.'TN' - T0 ;
  540. TX = kops phi0 * TY ;
  541.  
  542. LTRA = prog -0.5  pas 0.1  0.5 ;
  543. opti isov ligne ;
  544. trac TX mt (cont mt) LTRA ;
  545.  
  546. * intensites radiatives
  547.  
  548. trac rv.INCO.'I1' 10 mt (cont mt) ;
  549. trac rv.INCO.'I2' 10 mt (cont mt) ;
  550. trac rv.INCO.'I3' 10 mt (cont mt) ;
  551. trac rv.INCO.'I4' 10 mt (cont mt) ;
  552.  
  553. * vitesse
  554.  v1 = vect (rv.INCO.'UN') f_homo ux uy jaune ;
  555.  trac v1 mt (cont mt) ;
  556.  
  557. T1 = evol 'CHPO' rv.INCO.'TN' l1    ;
  558. T3 = evol 'CHPO' rv.INCO.'TN' l3inv ;
  559.  
  560. TAB1    = TABLE ;
  561. TAB1. 1 = 'MARQ CROI REGU ' ;
  562. TAB1. 2 = 'MARQ TRIA REGU ' ;
  563. TAB1.'TITRE'   = TABLE ;
  564. TAB1.'TITRE'. 1 = MOT ' PAROI_INFERIEURE' ;
  565. TAB1.'TITRE'. 2 = MOT ' PAROI_SUPERIEURE' ;
  566.  
  567. dess (T1 et T3) MIMA  LEGE TAB1
  568.     'TITR' 'TEMPERATURE'
  569.     'TITX' ' x ' 'TITY' ' T ' ;
  570.  
  571. *opti donn 5 ;
  572.  
  573. * evolution des Nusselt sur la paroi chaude
  574.  
  575. Nu_Tl = evol 'CHPO' Nu_T SCAL l4  ;
  576. Nu_Il = evol 'CHPO' Nu_I SCAL l4 ;
  577.  
  578. Nu_l = Nu_Tl + Nu_Il ;
  579.  
  580.  TAB3   = TABLE ;
  581.  TAB3.1 = 'MARQ CROI REGU ' ;
  582.  TAB3.2 = 'MARQ TRIA REGU ';
  583.  TAB3.3 = 'MARQ LOSA REGU '      ;
  584.  TAB3.'TITRE'   = TABLE ;
  585.  TAB3.'TITRE'. 1 = MOT 'CONDUCTION' ;
  586.  TAB3.'TITRE'. 2 = MOT 'RAYONNEMENT';
  587.  TAB3.'TITRE'. 3 = MOT 'TOTAL'      ;
  588.  
  589.  dess (Nu_Tl et Nu_Il et Nu_l) LEGE TAB3
  590.     'TITR' 'Nombre de Nusselt local sur la paroi chaude '
  591.     'TITX' ' y ' 'TITY' ' Nusselt ' ;
  592.  
  593. *opti donn 5 ;
  594.  
  595. finsi ;
  596.  
  597. fin ;
  598.  
  599.  
  600.  
  601.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales