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

Retour à la liste

Numérotation des lignes :

  1. * fichier : cvry-2D-1.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4.  
  5. ****************** CONVECTION-RAYONNEMENT ************************* *
  6. * *
  7. * Couplage convection naturelle laminaire/rayonnement milieu absorbant*
  8. * *
  9. * Il s'agit du test de De Vahl&Davis -cavité carrée- dans lequel le *
  10. * milieu contenu dans la cavité est radiativement absorbant. *
  11. * Les parois de la cavité sont supposées grises et diffusantes. *
  12. * *
  13. * Les données correspondent à celles de 2 cas tests de la littérature *
  14. * pour le test de validation: Rayleigh=1e5 et Nombre_rayonnement=1. *
  15. * *
  16. * Formulation EF *
  17. * Algorithme: implicite en régime permanent *
  18. * convection traitée par pénalisation *
  19. * rayonnement traité par la methode des harmoniques *
  20. * sphériques (P1) *
  21. * *
  22. * Cet algorithme est limité aux nombres de Rayleigh faibles *
  23. * (inférieurs à 1.e5 pour le test de Tan). Au dela d'autres *
  24. * algorithnes sont nécessaires. *
  25. * *
  26. * opérateurs utilisés : DUDW, KONV, LAPN, MDIA *
  27. * *
  28. * References: *
  29. * Z.Tan J.R.Howell, I.J.H.M.T Vol.34 1991 *
  30. * A.Yucel S.Acharya M.L.Williams, Num.Heat.Transfer Vol.15 1989 *
  31. * Rapports DMT/96.030 et DMT/96.059 *
  32. * *
  33. ***********************************************************************
  34.  
  35. tan = vrai ;
  36. complet = faux ;
  37. graph = faux ;
  38.  
  39. option echo 0 dime 2 elem qua4 ;
  40.  
  41. * test de Tan
  42. si tan ;
  43. mess ' test de Tan ' ;
  44. Ra = 1.e5 ;
  45. Nr = 1.0 ; om = 0. ;
  46. phi0 = 1. ; th = 2. ; tc = 1. ;
  47. Pr = 0.71 ; Pe = 1. ; Re = Pe/Pr ;
  48.  
  49. * test de Yucel
  50. sinon ;
  51. mess ' test de Yucel';
  52.  
  53. Ra = 5.e6 ;
  54. Nr = 0.04 ; om = 0. ;
  55. phi0 = 1.5 ; th = 4./3. ; tc = 2./3. ;
  56. Pr = 0.72 ; Re = 1. ; Pe = Re * Pr ;
  57.  
  58. finsi ;
  59.  
  60. * valeurs de reference pour les Nusselt - version du 7 nov. 97
  61. * test de Tan
  62. si complet ;
  63. Nu_Tref = 3.10 ; Nu_Iref = 9.50 ;
  64. sinon;
  65. Nu_Tref = 3.00 ; Nu_Iref = 9.41 ;
  66. finsi ;
  67.  
  68. OMEGA = 0.6 ;
  69. EPS = 1.E-4 ;
  70. EPSS = 1.E-10;
  71. maxiter = 30 ;
  72.  
  73. *------------------ geometrie -------------------*
  74.  
  75.  
  76. p1= 0. 0. ;
  77. p2= 0.5 0. ;
  78. p3= 0.5 0.5 ;
  79. p4= 0. 0.5 ;
  80. p5= 1. 0. ;
  81. p6= 0. 1. ;
  82.  
  83. * maill. fin
  84. *di = 0.005 ; df = 0.055 ;
  85. *vi = 0.011 ; vf = df ;
  86.  
  87. si complet ;
  88. * maill. moyen
  89. di = 0.005 ; df = 0.1 ;
  90. vi = 0.02 ; vf = df ;
  91. sinon ;
  92. * maill. grossier
  93. di = 0.01 ; df = 0.2 ;
  94. vi = 0.1 ; vf = df ;
  95. finsi;
  96.  
  97. * ligne horiz.
  98. c1 = droit p1 p2 dini di dfin df ;
  99. c1b = droit p2 p5 dini df dfin di ;
  100.  
  101. * ligne vert.
  102. c4 = droit p6 p4 dini vi dfin vf ;
  103. c4b = droit p4 p1 dini vf dfin vi ;
  104.  
  105. l1 = c1 et c1b ;
  106. l4 = c4 et c4b ;
  107. l3inv= l1 plus (0. 1.) ; l3 = inve ( l3inv ) ;
  108. l2inv= l4 plus (1. 0.) ; l2 = inve ( l2inv ) ;
  109. fron = l1 et l2 et l3 et l4 ;
  110. l13 = l1 et l3 ;
  111. elim 1e-5 (l1 et l2 et l3 et l4 ) ;
  112.  
  113. mt = l1 l2 l3 l4 dalle plan ;
  114. elim 1.e-5 mt ;
  115. tass mt ;
  116. orient mt ;
  117. si graph ; trac mt ; finsi ;
  118.  
  119. * ------------- tables domaine ----------------- *
  120.  
  121. $mt = doma mt impr ;
  122. $l13 = doma l13 'INCL' $mt 0.0001 impr ;
  123. $fron = doma fron 'INCL' $mt 0.0001 impr ;
  124.  
  125. * ------------- parametres physiques ----------- *
  126.  
  127. * Re : nombre de Reynolds
  128. * Ra : nombre de Rayleigh
  129. * Pr : nombre de Prandtl
  130. * Pe : nombre de Peclet = Pr * Re
  131. * Nr : nombre caracteristique du rayonnement
  132.  
  133. * lambda: conductivite
  134. * h : coefficient d'echange
  135. * t0 : temperature de reference
  136. * th : hot temperature
  137. * tc : cold temperature
  138. * emis : emissivite des parois
  139. * a : rapport d'allongement (largeur/hauteur)
  140. * beta : coefficient d'extinction (coeff.diffus + coeff.absorb)
  141. * tau : epaisseur optique (beta * l)
  142. * om : albedo ( coeff.diffus /(coeff.diffus + coeff.absorb) )
  143. * phi0 : paramètre de forme
  144. *
  145. * ------------- parametres numeriques ---------- *
  146.  
  147. * OMEGA : facteur de relaxation
  148. * EPS : precision sur les itérations internes
  149. * EPSS : facteur de penalisation
  150.  
  151. a = 1. ; tau = 1. ;
  152. emis = 1. ;
  153.  
  154. t0 = (th+tc)/2. ;
  155.  
  156. Gr = Ra/Pr ; gg = Gr/(Re*Re) ;
  157.  
  158. ALFV = 1./Re ;
  159. ALFM = 1 ;
  160.  
  161. * ---------- coefficients algorithme ----------- *
  162.  
  163. E = emis/(4.-2.*emis) ;
  164. H = (3.*(a*tau)) ;
  165. H1 = (3.*((a*tau)**2))*(1.-om) ;
  166. HE = H*E ;
  167.  
  168. 3Nr = 3.*Nr*phi0 ;
  169.  
  170. C1 = (-1.*H1)/3Nr ;
  171. C2 = (-1.*HE)/3Nr ;
  172.  
  173. * ---------- initialisation de la température--- *
  174.  
  175. RT = EQEX $MT 'OPTI' 'EF' 'IMPL'
  176. ZONE $MT 'OPER' LAPN 1. INCO 'TT'
  177. ;
  178. RT = EQEX RT
  179. CLIM
  180. 'TT' TIMP l4 th
  181. 'TT' TIMP l2 tc
  182. ;
  183. rt.INCO = TABLE 'INCO' ;
  184. rt.INCO.'TT' = kcht $mt scal sommet t0 ;
  185.  
  186. rt.'NITER' = 1 ;
  187. rt.'OMEGA' = 1. ;
  188.  
  189. exec rt ;
  190.  
  191. * ----------------- algorithme------------------ *
  192.  
  193. RV = EQEX $MT 'OPTI' 'EF' 'IMPL'
  194. ZONE $MT 'OPER' DUDW EPSS INCO 'UN'
  195. ZONE $MT 'OPER' KONV 1. 'UN' ALFV INCO 'UN'
  196. ZONE $MT 'OPER' LAPN ALFV INCO 'UN'
  197. ZONE $MT 'OPER' MDIA 'GB' INCO 'TN' 'UN'
  198. ;
  199. RV = EQEX RV 'OPTI' 'EF' 'IMPL'
  200. ZONE $MT 'OPER ' KONV Pe 'UN' ALFM INCO 'TN'
  201. ZONE $MT 'OPER' LAPN ALFM INCO 'TN'
  202. ZONE $MT 'OPER' MDIA 'C1' INCO 'IN' 'TN'
  203. ZONE $MT 'OPER' MDIA 'KN' INCO 'TN'
  204. ;
  205. RV = EQEX RV 'OPTI' 'EF' 'IMPL'
  206. ZONE $MT 'OPER' LAPN 1. INCO 'IN'
  207. ZONE $MT 'OPER' MDIA 'H1' INCO 'IN'
  208. ZONE $MT 'OPER' MDIA 'K' INCO 'TN' 'IN'
  209. ;
  210. RV = EQEX RV 'OPTI' 'EF' 'IMPL'
  211. ZONE $fron 'OPER' MDIA HE INCO 'IN'
  212. ZONE $fron 'OPER' MDIA 'KI' INCO 'TN' 'IN'
  213. ZONE $l13 'OPER' MDIA 'KT' INCO 'TN'
  214. ZONE $l13 'OPER' MDIA 'C2' INCO 'IN' 'TN'
  215. ;
  216. RV = EQEX RV
  217. CLIM
  218. 'UN' UIMP fron 0.0
  219. 'UN' VIMP fron 0.0
  220. 'TN' TIMP l4 th
  221. 'TN' TIMP l2 tc
  222. ;
  223. rv.INCO = TABLE 'INCO' ;
  224. rv.INCO.'TN' = kcht $mt scal sommet rt.INCO.'TT' ;
  225. rv.INCO.'IN' = kcht $mt scal sommet 1. ;
  226. rv.INCO.'GB' = kcht $mt vect sommet (0. (-1.*phi0*gg)) ;
  227. rv.INCO.'UN' = kcht $mt vect sommet (1.e-5 1.e-5) ;
  228.  
  229. rv.'NITER' = 1 ;
  230. rv.'OMEGA' = OMEGA ;
  231. rv.'IMPR' = 1 ;
  232. rv.'ITMA' = 1 ;
  233. rv.INCO.'C1'= C1;
  234. rv.INCO.'C2'= C2;
  235. rv.INCO.'H1'= H1;
  236.  
  237. repeter bloc2 ;
  238.  
  239. nbiter = &bloc2 ;
  240.  
  241. * -----------actualisation des coefficients----- *
  242.  
  243. TN3 = kops (rv.INCO.'TN') ** 3 ;
  244. 4H1 = -4.*H1 ;
  245.  
  246. rv.INCO.'K' = kops 4H1 * TN3 ;
  247. rv.INCO.'KN'= kops ((-1. * 4H1) / 3Nr ) * TN3 ;
  248.  
  249. 4HE = -4.*HE ;
  250. 4HEN = (4.*HE)/3Nr ;
  251.  
  252. rv.INCO.'KI' = kcht $fron scal sommet (kops 4HE * TN3) ;
  253. rv.INCO.'KT' = kcht $l13 scal sommet (kops 4HEN * TN3) ;
  254.  
  255. detr TN3 ;
  256.  
  257. UN = copie rv.INCO.'UN' ;
  258. TN = copie rv.INCO.'TN' ;
  259. IN = copie rv.INCO.'IN' ;
  260.  
  261. exec rv ;
  262.  
  263. *XU = (maxi (abs (kops ((rv.INCO.'UN') - UN) / UN))) ;
  264. XT = (maxi (abs (kops ((rv.INCO.'TN') - TN) / TN))) ;
  265. XI = (maxi (abs (kops ((rv.INCO.'IN') - IN) / IN))) ;
  266.  
  267. MESS ' ERREUR RELATIVE TN ' XT 'IN ' XI ;
  268.  
  269. Tmax = maxi (rv.INCO.'TN') ;
  270. Tmin = mini (rv.INCO.'TN') ;
  271.  
  272. *MESS ' EXTREMUM TN MAX' Tmax 'MIN ' Tmin ;
  273. * mess ' iteration ' nbiter ;
  274.  
  275. si ((XT < EPS) et (XI < EPS)) ;
  276. quit bloc2 ;
  277. sinon ;
  278. si (nbiter > maxiter) ;
  279. mess ' pas de convergence' ; erre 5 ;
  280. finsi ;
  281. finsi ;
  282.  
  283. fin bloc2 ;
  284.  
  285. *-------------post-traitement------------------------------------------
  286.  
  287. * Nusselt sur la paroi chaude
  288. * ---------------------------
  289. *
  290. * calcul des Nusselt dus a la conduction et au rayonnement
  291. * sur la paroi chaude (l4)
  292.  
  293. $l4 = doma l4 ;
  294.  
  295. * Nusselt du a la conduction
  296.  
  297. gradt = KOPS rv.'INCO'.'TN' 'GRAD' $mt ;
  298. dTdXe = KCHT $mt 'SCAL' 'CENTRE' (exco 'UX' gradt) ;
  299. dTdX = ELNO $MT (kops (-1.) * dTdXe) ;
  300.  
  301. Nu_T = KCHT $l4 'SCAL' 'SOMMET' (kops phi0 * dTdX) ;
  302.  
  303. * Nusselt du au rayonnement
  304.  
  305. I_l4 = redu rv.INCO.'IN' l4 ;
  306. T_l4 = redu rv.INCO.'TN' l4 ;
  307.  
  308. T4 = kops T_l4 ** 4 ;
  309. 4T4 = kops 4 * T4 ;
  310. dIdn = kops HE * (4T4 - I_l4) ;
  311.  
  312. Nu_I = KCHT $l4 'SCAL' 'SOMMET' (kops dIdn / (3.*Nr)) ;
  313.  
  314. * Nusselt global pour la paroi
  315.  
  316. dy = 'DOMA' $l4 'VOLUME' ;
  317. dyt = SOMT dy ;
  318. *message 'LONGUEUR DE LA PAROI = ' dyt ;
  319.  
  320. Nu_Te = noel $l4 Nu_T ;
  321. Nu_Tdy = kops Nu_Te '*' dy ;
  322. Nu_Tdyt = somt Nu_Tdy ;
  323. Nu_Tm = Nu_Tdyt / dyt ;
  324. message 'Nusselt de conduction moyen: ' Nu_Tm ;
  325.  
  326. Nu_Ie = noel $l4 Nu_I ;
  327. Nu_Idy = kops Nu_Ie '*' dy ;
  328. Nu_Idyt = somt Nu_Idy ;
  329. Nu_Im = Nu_Idyt / dyt ;
  330. message 'Nusselt de rayonnement moyen: ' Nu_Im ;
  331.  
  332. * comparaison aux valeurs de reference (version du 6 nov.97)
  333.  
  334. resT = abs((Nu_Tm-Nu_Tref/Nu_Tref));
  335. resI = abs((Nu_Im-Nu_Iref/Nu_Iref));
  336. si ((resT < 1e-2) et (resI < 1e-2 )) ; erre 0 ;
  337. sinon ; erre 5 ; finsi ;
  338.  
  339. *-------------graphique------------------------------------------------
  340.  
  341. si graph ;
  342.  
  343. si tan ;
  344. * test de Tan
  345. * opti titre 'CONDUCTION - RAYONNEMENT Ra = 1.e5 Nr = 1. ';
  346. f_homo = 0.001 ;
  347.  
  348. sinon ;
  349.  
  350. * test de Yucel
  351. * opti titre 'CONVECTION - RAYONNEMENT Ra = 5.e6 Nr = 0.04 ';
  352. f_homo = 0.0003 ;
  353.  
  354. finsi ;
  355.  
  356. * température relative comprise entre -0.5 et 0.5
  357.  
  358. TY = kops rv.INCO.'TN' - T0 ;
  359. TX = kops phi0 * TY ;
  360.  
  361. LTRA = prog -0.5 pas 0.1 0.5 ;
  362.  
  363. *trac TX mt (cont mt) LTRA ;
  364. opti isov ligne ;
  365. trac TX mt (cont mt) LTRA ;
  366.  
  367. * intensite radiative
  368.  
  369. trac rv.INCO.'IN' 10 mt (cont mt) ;
  370.  
  371. * vitesse
  372. v1 = vect (rv.INCO.'UN') f_homo ux uy jaune ;
  373. trac v1 mt (cont mt) ;
  374.  
  375. T1 = evol 'CHPO' rv.INCO.'TN' l1 ;
  376. T3 = evol 'CHPO' rv.INCO.'TN' l3inv ;
  377.  
  378. *TAB1=TABLE ;
  379. *TAB1.1= 'MARQ CROI REGU TITR PAROI_INFERIEURE' ;
  380. *TAB1.2= 'MARQ TRIA REGU TITR PAROI_SUPERIEURE' ;
  381.  
  382. dess (T1 et T3) MIMA ;
  383.  
  384. * evolution des Nusselt sur la paroi chaude
  385.  
  386. Nu_Tl = EVOL 'CHPO' Nu_T SCAL l4 ;
  387.  
  388. *TAB2 = TABLE ;
  389. *TAB2.1 = 'MARQ LOSA' ;
  390.  
  391. DESS Nu_Tl 'MIMA' ;
  392. *DESS Nu_Tl 'MIMA' TAB2
  393. * 'TITR' ' Ra =1.e5 Nombre de Nusselt local sur la paroi chaude '
  394. * 'TITX' ' y ' 'TITY' ' Nusselt ' ;
  395.  
  396. Nu_Il = evol 'CHPO' Nu_I l4 ;
  397.  
  398. dess Nu_Il 'MIMA' ;
  399.  
  400. * Nusselt total
  401.  
  402. Nu_l = Nu_Tl + Nu_Il ;
  403.  
  404. *TAB3 = TABLE ;
  405. *TAB3.1 = 'MARQ CROI REGU ' ;
  406. *TAB3.2 = 'MARQ TRIA REGU ';
  407. *TAB3.3 = 'MARQ LOSA REGU ' ;
  408. *TAB3.'TITRE' = TABLE ;
  409. *TAB3.'TITRE'. 1 = MOT 'CONDUCTION' ;
  410. *TAB3.'TITRE'. 2 = MOT 'RAYONNEMENT';
  411. *TAB3.'TITRE'. 3 = MOT 'TOTAL' ;
  412.  
  413. dess (Nu_Tl et Nu_Il et Nu_l) ;
  414.  
  415. *dess (Nu_Tl et Nu_Il et Nu_l) LEGE TAB3
  416. * 'TITR' 'Nombre de Nusselt local sur la paroi chaude '
  417. * 'TITX' ' y ' 'TITY' ' Nusselt ' ;
  418.  
  419. finsi ;
  420.  
  421. fin ;
  422.  
  423.  
  424.  
  425.  
  426.  
  427.  

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