Télécharger rayo_abs-axi-2.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : rayo_abs-axi-2.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4.  
  5. complet=faux;
  6. ****************************************************************
  7. * Calcul de la température d'une cavité contenant
  8. * un milieu absorbant (Tg,k_abs)
  9. *
  10. * DONNEES
  11. * cavité cylindrique de cote 1.
  12. *
  13. * 1000K e=1.0
  14. * -----
  15. * axe | | e=0.5
  16. * | |
  17. * -----
  18. * 2000K e=0.1
  19. *
  20. *
  21. * RESULTATS: flux rayonné total sur la cavité
  22. * solution numérique de référence
  23. *
  24. * On teste 2 méthodes:
  25. * 1- la méthode basée sur le calcul de la matrice de rayonnement
  26. * (opérateur RAYE) R telle que Phi = R. T**4 - R. Tg**4
  27. * 2- la méthode basée sur le calcul de la température de
  28. * rayonnement Trad telle que Phi = emis*sigma*(T**4-Trad**4)
  29. *
  30. ****************************************************************
  31.  
  32. *** Options ...
  33.  
  34. * option echo 0 dime 2 elem qua8 mode axis;
  35. option echo 0 dime 2 elem qua4 mode axis;
  36.  
  37. graph = faux ;
  38.  
  39. *** Solutions numeriques de référence à la creation du test
  40. * flux total rayonné sur la cavité évalué avec les méthodes 1 et 2
  41.  
  42. si (complet); n = 100 ;
  43. fref1 = -6.83843E6;
  44. fref2 = -6.83870E6;
  45. sinon; n= 20 ;
  46. fref1 = -6.95186E6;
  47. fref2 = -6.95813E6;
  48. finsi;
  49.  
  50. *** Paramètres ...
  51.  
  52. l = 1.0 ;
  53.  
  54. * epaisseur des parois : dz = 10 mm
  55. *<
  56. * dz = 0.01 ;
  57. dz = 0.05 ;
  58. *>
  59. dzp = dz ;
  60. dzn = -1. * dz ;
  61.  
  62. * proprietes physiques
  63. lam = 10.;
  64. *>
  65. e_sup=1.0 ; e_inf=0.1 ; e_late = 0.5 ;
  66. * e_sup=1.0 ; e_inf=1.0 ; e_late = 1.0 ;
  67. T_sup = 1000. ; T_inf = 2000. ;
  68. *<
  69. sig = 5.67e-8;
  70. Nr = sig* (T_sup**3.) * L /lam ; mess 'Nr= ' Nr;
  71.  
  72. * milieu absorbant : coeff et température (K)
  73. k_abs = -1.e-0 ;
  74. Tg = 2500. ;
  75.  
  76. *** Points ...
  77.  
  78. dens l ;
  79. * n = 4 ;
  80.  
  81. p4 = 0. l ; p3 = l l ;
  82. p1 = 0. 0. ; p2 = l 0. ;
  83. q4 = 0. l ; q3 = l l ;
  84. q1 = 0. 0. ; q2 = l 0. ;
  85. l1 = d n p1 p2 ; l2 = d n p3 p4 ;
  86. l3 = d n q2 q3 ;
  87. late = l3 ;
  88.  
  89. * la surface laterale est disjointe des deux autres
  90. cavite = l1 et l2 et late ;
  91. *<
  92. elim 1e-5 cavite ;
  93. *>
  94.  
  95. ne = 5;
  96. z1 = l1 trans ne (0. dzn) ;
  97. z2 = l2 trans ne (0. dzp) ;
  98. z3 = l3 trans ne (dzp 0.) ;
  99.  
  100. tout= z1 et z2 et z3 ;
  101. titr 'Le maillage de la cavite' ;
  102. si(graph) ; trac tout ; finsi ;
  103.  
  104. *** Modélisation ...
  105.  
  106. * conduction
  107. lamb = lam/dz ;
  108.  
  109. mcd = modeli tout thermique ;
  110. k = mate mcd 'K' lamb ;
  111. cnd = cond mcd k ;
  112.  
  113. *
  114. mr1= modeli l1 thermique rayonnement 'CAVITE' 'CONVEXE'
  115. 'CONS' 'CAV1';
  116. mr2= modeli l2 thermique rayonnement 'CAVITE' 'CONVEXE'
  117. 'CONS' 'CAV1';
  118. mrl= modeli late thermique rayonnement'CAVITE' 'CONVEXE'
  119. 'CONS' 'CAV1' ;
  120. mrt = mr1 et mr2 et mrl ;
  121.  
  122. *** Emissivités
  123.  
  124. e1 = mate mr1 'EMIS' e_inf 'CABS' k_abs 'TABS' Tg ;
  125. e2 = mate mr2 'EMIS' e_sup 'CABS' k_abs 'TABS' Tg ;
  126. el = mate mrl 'EMIS' e_late 'CABS' k_abs 'TABS' Tg ;
  127. e = e1 et e2 et el ;
  128.  
  129. *** Calcul des facteurs de forme
  130.  
  131. * facteurs de forme "classiques" pour memoire
  132. * ffc = ffor mrt e ;
  133. * facteurs de forme généralisés
  134. fft = ffor mrt e ;
  135.  
  136. *** Conditions aux limites ...
  137.  
  138. c1 = bloq z1 'T' ;
  139. tim1 = depi c1 T_inf ;
  140. c2 = bloq z2 'T' ;
  141. tim2 = depi c2 T_sup ;
  142.  
  143. *** méthode 1
  144. * ---------
  145.  
  146. * calcul de la matrice de rayonnement
  147.  
  148. chamr = raye mrt fft e ;
  149.  
  150.  
  151. *** Initialisation de la température ...
  152.  
  153. tp = manu chpo tout 1 'T' T_inf natu 'DIFFUS';
  154.  
  155. * calcul du flux associé au milieu absorbant
  156. *<
  157. Tg_abs = manu chpo tout 1 'T' Tg natu 'DIFFUS';
  158. T_abs = redu (exco 'T' Tg_abs 'T') cavite ;
  159. Te_abs = chan 'CHAM' T_abs mrt 'GRAVITE' ;
  160.  
  161. cr_abs = rayn mrt chamr Te_abs ;
  162.  
  163. Fg_abs = cr_abs * T_abs ;
  164. * list Fg_abs;
  165. *>
  166.  
  167.  
  168. *** Résolution (par itérations) ...
  169.  
  170. * Coeff. de relaxation ...
  171. alfa = 0.42 ;
  172.  
  173. nbpsup = nbno tout ;
  174.  
  175. listemp = prog ;
  176. listres = prog ;
  177.  
  178. maxiter = 100 ;
  179. critconv = 1.e-5 ;
  180.  
  181. REPE bloc1 maxiter ;
  182.  
  183. nbiter = &bloc1 ;
  184.  
  185. t_cavi = redu (exco 'T' tp 'T') cavite ;
  186. te_cavi = chan 'CHAM' t_cavi mrt 'GRAVITE' ;
  187. cr = rayn mrt chamr te_cavi ;
  188.  
  189. cndtot = cnd et cr et c1 et c2 ;
  190. residu = (cndtot * tp) - Fg_abs ;
  191.  
  192. normres = ((xtx (exco 'Q' residu)) / nbpsup) ** 0.5 ;
  193. * mess ' La norme du flux résiduel = ' normres ;
  194. si((nbiter > 1) et (normres < critconv)) ;
  195. mess 'Convergence à l itération ' (nbiter-1) ;
  196. quitter bloc1 ;
  197. finsi ;
  198. si(nbiter > maxiter) ;
  199. mess 'Erreur ! Pas de convergence après ' (nbiter-1)
  200. ' itérations !' ;
  201. erre 2 ;
  202. quitter bloc1 ;
  203. finsi ;
  204. listres = listres et (prog normres) ;
  205.  
  206. * mess '---------------------------------------' ;
  207. * mess 'Itération N° ' &bloc1 ;
  208.  
  209. tt1 = resou cndtot (tim1 et tim2 et Fg_abs) ;
  210.  
  211. dt = exco 'T' (tt1- tp) 'T' ;
  212. normdt = ((xtx dt) / nbpsup) ** 0.5 ;
  213. * mess ' La norme de delta t = ' normdt ;
  214.  
  215. tn = (alfa * tt1) + ((1.-alfa) * tp) ;
  216. tp =tn ;
  217.  
  218. FIN bloc1 ;
  219. *
  220.  
  221.  
  222. *** Post-traitement ...
  223.  
  224. mess '1- calcul avec matrice de rayonnement';
  225.  
  226. * tsol = extr tt1 'T' q2 ;
  227. * mess 'La solution obtenue = ' tsol ;
  228.  
  229. * flux rayonne
  230.  
  231. fray = (cr * tt1) - Fg_abs ;
  232. fray_1 = maxi (resu fray);
  233. mess ' bilan global cavite: ' fray_1 ;
  234.  
  235. mess ' flux surface inf: ' (maxi (resu (redu fray l1))) ;
  236. mess ' flux surface sup: ' (maxi (resu (redu fray l2))) ;
  237. mess ' flux surface late: ' (maxi (resu (redu fray late))) ;
  238.  
  239. * flux associé à la condition de température imposée
  240.  
  241. mess ' Timp surf. inf.: ' (maxi (resu (reac c1 tt1))) ;
  242. mess ' Timp surf. sup.: ' (maxi (resu (reac c2 tt1))) ;
  243.  
  244. t = exco 'T' tt1 'T' ;
  245. titr 'Le champ de temperature final' ;
  246. si(graph) ; trac tout t ; finsi ;
  247. tlate = redu t late ; comm list tlate ;
  248.  
  249. *! trac (redu tt1 z3 ) z3 ;
  250.  
  251. ev1 = evol 'CHPO' (redu t l3) l3;
  252. ev1 = ev1 coul 'VERT' ;
  253.  
  254. RESI1 =ABS((fray_1 - fref1)/fref1);
  255. mess 'Erreur relative - methode 1' (100*RESI1) '%' ;
  256.  
  257. * opti donn 5;
  258.  
  259. * méthode 2
  260. * ---------
  261.  
  262. *** Initialisation de la température ...
  263.  
  264. tp = manu chpo tout 1 'T' T_inf natu 'DIFFUS';
  265.  
  266. *** Résolution (par itérations) ...
  267.  
  268. * Coeff. de relaxation ...
  269.  
  270. alfa = 1.0 ;
  271.  
  272. nbpsup = nbno tout ;
  273.  
  274. listemp = prog ;
  275. listres = prog ;
  276.  
  277. maxiter = 100 ;
  278. critconv = 1.e-5 ;
  279.  
  280. REPE bloc2 ;
  281.  
  282. nbiter = &bloc2 ;
  283.  
  284. t_cavi = redu (exco 'T' tp 'T') cavite ;
  285. te_cavi = chan 'CHAM' t_cavi mrt 'GRAVITE' ;
  286.  
  287. * calcul de la temperature de rayonnement associée: trad
  288. trad = raye mrt fft e te_cavi 1e-7 ;
  289.  
  290. * calcul du coefficient d'echange
  291. hrad = HRCAV mrt e te_cavi trad ;
  292.  
  293. trad_n1= chan 'CHPO' mrt trad ;
  294. trad_n = nomc trad_n1 'T' 'NATU' 'DIFFUS';
  295.  
  296. * pour la condition de convection
  297. cr = cond mrt hrad;
  298. f = conv mrt hrad trad_n;
  299.  
  300. cndtot = cnd et c1 et c2 et cr ;
  301.  
  302. residu = (cndtot * tp) -f ;
  303. normres = ((xtx (exco 'Q' residu)) / nbpsup) ** 0.5 ;
  304. * mess ' La norme du flux résiduel = ' normres ;
  305. si((nbiter > 1) et (normres < critconv)) ;
  306. mess 'Convergence à l itération ' (nbiter-1) ;
  307. quitter bloc2 ;
  308. finsi ;
  309. si(nbiter > maxiter) ;
  310. mess 'Erreur ! Pas de convergence après ' (nbiter-1)
  311. ' itérations !' ;
  312. erre 2 ;
  313. quitter bloc2 ;
  314. finsi ;
  315. listres = listres et (prog normres) ;
  316.  
  317. * mess '---------------------------------------' ;
  318. * mess 'Itération N° ' &bloc2 ;
  319.  
  320. tt2 = resou cndtot (tim1 et tim2 et f) ;
  321.  
  322. dt = exco 'T' (tt2- tp) 'T' ;
  323. normdt = ((xtx dt) / nbpsup) ** 0.5 ;
  324. * mess ' La norme de delta t = ' normdt ;
  325.  
  326. tn = (alfa * tt2) + ((1.-alfa) * tp) ;
  327. tp =tn ;
  328.  
  329. * mess ' Température obtenue = ' (extr tt2 'T' q2) ;
  330.  
  331. FIN bloc2 ;
  332.  
  333.  
  334. *temps impr place sgac ;
  335.  
  336.  
  337. *** Post-traitement ...
  338.  
  339. mess '2- calcul avec température de rayonnement';
  340.  
  341. t = exco 'T' tt2 'T' ;
  342. titr 'Le champ de temperature final' ;
  343. si(graph) ; trac tout t ; finsi ;
  344.  
  345. tlate = redu t late ;
  346.  
  347. * mess ' T laterale: min et max ' (mini tlate) (maxi tlate);
  348.  
  349. ev2 = evol 'CHPO' (redu t l3) l3;
  350. ev2 = ev2 coul 'BLEU' ;
  351.  
  352. * flux rayonne
  353.  
  354. fray = (cr * tt2)- f;
  355. fray_2 = maxi (resu fray);
  356. mess ' bilan global cavite: ' fray_2 ;
  357. mess ' flux surface inf: ' (maxi (resu (redu fray l1))) ;
  358. mess ' flux surface sup: ' (maxi (resu (redu fray l2))) ;
  359. mess ' flux surface late: ' (maxi (resu (redu fray late))) ;
  360.  
  361. * flux associé à la condition de température imposée
  362.  
  363. mess ' Timp surf. inf.: ' (maxi (resu (reac c1 tt2))) ;
  364. mess ' Timp surf. sup.: ' (maxi (resu (reac c2 tt2))) ;
  365.  
  366.  
  367. RESI2 =ABS((fray_2 - fref2)/fref2);
  368. mess 'Erreur relative - methode 2' (100*RESI2) '%' ;
  369.  
  370. SI ((RESI1 &lt;EG 1E-4) et (RESI2 &lt;EG 1E-4));
  371. ERRE 0;
  372. SINO;
  373. ERRE 5;
  374. FINSI;
  375.  
  376. si graph;
  377. TAB1 = table;
  378. TAB1. 'TITRE' = table;
  379. TAB1. 'TITRE' . 1 = 'MOT' 'methode 1' ;
  380. TAB1. 'TITRE' . 2 = 'MOT' 'methode 2' ;
  381.  
  382. dess (ev1 et ev2) mima lege TAB1 ;
  383. finsi;
  384.  
  385. fin;
  386.  
  387.  
  388.  
  389.  
  390.  
  391.  
  392.  
  393.  
  394.  
  395.  
  396.  

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