Télécharger rayo-2D-2.dgibi

Retour à la liste

Numérotation des lignes :

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

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