Télécharger rayoh-3D.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : rayoh-3D.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4.  
  5. ************************************************************************
  6. * *
  7. * Rayonnement thermique en milieu transparent: *
  8. * *
  9. * Test de bon fonctionnement de la procédure HRAYO permettant de *
  10. * traiter le rayonnement face à face par comparaison à une solution *
  11. * analytique en régime permanent. *
  12. * *
  13. * Calcul 3D. *
  14. * *
  15. * on considère 3 plaques parallèles rayonnantes : *
  16. * *
  17. * -------- T3,e3 milieu infini *
  18. * *
  19. * e2sup *
  20. * -------- T2 *
  21. * e2inf face à face *
  22. * *
  23. * -------- T1,e1 *
  24. * *
  25. * La température des plaques supérieure et inférieure est imposée et *
  26. * uniforme (respectivement 1000K et 2000K) et on détermine la tempéra- *
  27. * ture d'équilibre de la plaque intermédiaire. *
  28. * Le rayonnement entre les plaques 1 et 2 est modélisé par un rayonne- *
  29. * ment face à face, celui entre les plaques 2 et 3 par un rayonnement *
  30. * avec un milieu infini, le milieu infini étant representé par la *
  31. * plaque 3 non maillée. *
  32. * *
  33. * paramètre xquad vrai : éléments quadratiques *
  34. * faux : éléments linéaires *
  35. * Remarque: dans le cas quadratique, et pour le face à face chaque *
  36. * frontière doit etre un maillage simple (TRI6 uniquement *
  37. * ou QUA8 uniquement) *
  38. ************************************************************************
  39.  
  40. graph = faux ;
  41.  
  42. xquad = faux ;
  43.  
  44. si (xquad) ;
  45. option echo 0 dime 3 elem pr15 ;
  46. sinon ;
  47. option echo 0 dime 3 elem cub8 ;
  48. finsi ;
  49.  
  50. *-------------------------------
  51. * Données physiques du problème
  52. *-------------------------------
  53.  
  54. T1 = 1000. ; T3 = 2000. ;
  55. e1 = 0.1 ; e2inf = 0.6 ; e2sup = 0.3 ; e3 = 0.9 ;
  56. lamb = 2.e3 ;
  57.  
  58. *-------------------------------
  59. * Solution analytique
  60. *-------------------------------
  61.  
  62. r1 = 1. - e1 ; r2i = 1. - e2inf ;
  63. r3 = 1. - e3 ; r2s = 1. - e2sup ;
  64. alpha1= (e1*e2inf) / (1.-(r1*r2i)) ;
  65. alpha3= (e3*e2sup) / (1.-(r3*r2s)) ;
  66. M1 = T1**4. ; M3 = T3**4. ;
  67. M2 = ((alpha1*M1) + (alpha3*M3)) / (alpha1+alpha3) ;
  68. T2_ana = M2**(0.25) ;
  69.  
  70. *-------------------------------
  71. * Maillage
  72. *-------------------------------
  73.  
  74. oeil=10 10 10;
  75. p1=0 0 0 ; p2=1 0 0 ; p3=1 1 0 ; p4=0 1 0;
  76. l1=p1 d 1 p2;
  77. l2=p2 d 1 p3;
  78. l3=p3 d 1 p4;
  79. l4=p4 d 2 p1;
  80. si xquad ;
  81. * on recherche un objet simple
  82. opti elem tri6 ;
  83. s1= (l1 et l2 et l3 et l4) surf plan ;
  84. option echo 1 dime 3 elem pr15 ;
  85. sinon ;
  86. s1= (l1 et l2 et l3 et l4) surf plan ;
  87. finsi ;
  88.  
  89. n=1 ;
  90. v = 0 0 0.1 ; v1 = (-1)*v ;
  91. w = 0 0 0.001 ; w1 = (-1)*w ;
  92.  
  93.  
  94. s2inf = s1 ;
  95. s2sup = s2inf plus w ;
  96. plaq2 = s2inf volu n s2sup ;
  97.  
  98. s1sup = s2inf plus v1;
  99. s1inf = s1sup plus w1;
  100. plaq1 = s1inf volu n s1sup ;
  101.  
  102. rac12 = liai 0.1 s1sup s2inf ;
  103.  
  104. tout = plaq1 et plaq2 ;
  105.  
  106. si graph ; trac oeil tout ; finsi ;
  107.  
  108. *-------------------------------
  109. * Conduction
  110. *-------------------------------
  111.  
  112. mcd = mode tout thermique ;
  113. k = mate mcd 'K' lamb ;
  114. cnd = cond mcd k ;
  115.  
  116. *-------------------------------
  117. * Linéarisation du rayonnement :
  118. * - modèles de convection
  119. * - creation des SEG2 définissant les relations entre les supports
  120. * des champs définis sur s1sup et s2inf (pour le face à face)
  121. *-------------------------------
  122.  
  123.  
  124. * face à face
  125.  
  126. mcv12 = mode rac12 thermique convection ;
  127.  
  128. si xquad ;
  129.  
  130. x_s1 = (s1sup chan TRI6) chan TRI3 ;
  131. x_s2 = (s2inf chan TRI6) chan TRI3 ;
  132. xx_s1 = x_s1 chan ligne;
  133. xx_s2 = x_s2 chan ligne;
  134. x_s12 = liai 0.1 x_s1 x_s2;
  135. xx_s12 = x_s12 chan ligne ;
  136.  
  137. rel12 = diff xx_s12 (xx_s1 et xx_s2) ;
  138.  
  139. sinon ;
  140. sr12 = rac12 chan ligne ;
  141. sr1 = s1sup chan ligne ;
  142. sr2 = s2inf chan ligne ;
  143. rel12 = diff sr12 (sr1 et sr2) ;
  144.  
  145. finsi ;
  146.  
  147. si graph ; trac oeil rel12 ; finsi ;
  148.  
  149. * milieu infini
  150.  
  151. mcv2s = mode s2sup thermique convection ;
  152.  
  153.  
  154.  
  155. *-------------------------------
  156. * Mod�les de rayonnement et �missivit�s
  157. *-------------------------------
  158.  
  159. * face � face
  160. mr1 = mode s1sup thermique rayonnement FAC_A_FAC
  161. s1sup S2inf mcv12 rel12 CONS 'FAC1';
  162. emi1 = mate mr1 'EMIS' e1 ;
  163.  
  164. mr2i = mode s2inf thermique rayonnement FAC_A_FAC
  165. s1sup S2inf mcv12 rel12 CONS 'FAC1' ;
  166. emi2i= mate mr2i 'EMIS' e2inf ;
  167.  
  168. * milieu infini
  169. mr2s = mode s2sup thermique rayonnement 'INFINI' CONS 'INF1' ;
  170. emi2s= mate mr2s 'EMIS' e2sup 'E_IN' e3 ;
  171.  
  172. t3inf = manu chpo s2sup 1 'T' T3 ;
  173.  
  174. *-------------------------------
  175. * Conditions aux limites
  176. *-------------------------------
  177.  
  178. c1 = bloq plaq1 'T' ;
  179. tim1 = depi c1 T1 ;
  180.  
  181. *-------------------------------
  182. * Algorithme de résolution
  183. *-------------------------------
  184.  
  185. *** Initialisation de la température ...
  186.  
  187. tp = manu chpo tout 1 'T' T3 natu 'DIFFUS';
  188. nbpsup = nbno tout ;
  189.  
  190. *** Résolution (par itérations) ...
  191.  
  192. * Coeff. de relaxation ...
  193. alfa = 0.3 ;
  194.  
  195. listemp = prog ;
  196. listres = prog ;
  197.  
  198. maxiter = 100 ;
  199. critconv = 1.e-2 ;
  200. * opti echo 1 ;
  201. REPE bloc1 ;
  202.  
  203. nbiter = &bloc1 ;
  204.  
  205. * traitement du rayonnement
  206.  
  207. * face à face (1,2)
  208.  
  209. t1sup = redu (exco 'T' tp 'T') s1sup;
  210. t2inf = redu (exco 'T' tp 'T') s2inf;
  211. h12 = HRAYO mcv12 mr1 emi1 t1sup mr2i emi2i t2inf rel12 ;
  212.  
  213. * mr12 = MATE mcv12 'H' h12 ;
  214. * mr12 = h12 ;
  215. * cr12 = cond mcv12 mr12 ;
  216. cr12 = cond mcv12 h12 ;
  217.  
  218. * milieu infini (2,3)
  219.  
  220. t2sup = redu (exco 'T' tp 'T') s2sup;
  221. h23 = HRAYO mcv2s mr2s emi2s t2sup t3inf ;
  222.  
  223. * mr23 = MATE mcv2s 'H' h23 ;
  224. cr23 = cond mcv2s h23 ;
  225. f3 = conv mcv2s h23 t3inf ;
  226.  
  227. crtot= cr12 et cr23 ;
  228.  
  229. cndtot = crtot et cnd et c1 ;
  230. residu = (cndtot * tp) - f3 ;
  231. normres = ((xtx (exco 'Q' residu)) / nbpsup) ** 0.5 ;
  232. * mess ' La norme du flux résiduel = ' normres ;
  233. si((nbiter > 1) et (normres < critconv)) ;
  234. mess 'Convergence à l itération ' (nbiter-1) ;
  235. quitter bloc1 ;
  236. finsi ;
  237. si(nbiter > maxiter) ;
  238. mess 'Erreur ! Pas de convergence après ' (nbiter-1)
  239. ' itérations !' ;
  240. erre 2 ;
  241. quitter bloc1 ;
  242. finsi ;
  243. listres = listres et (prog normres) ;
  244.  
  245. * mess '---------------------------------------' ;
  246. * mess 'Itération N° ' &bloc1 ;
  247.  
  248. tt = resou cndtot (tim1 et f3) ;
  249. dt = exco 'T' (tt - tp) 'T' ;
  250. normdt = ((xtx dt) / nbpsup) ** 0.5 ;
  251. * mess ' La norme de delta t = ' normdt ;
  252.  
  253. tn = (alfa * tt) + ((1.-alfa) * tp) ;
  254. tp = tn ;
  255.  
  256.  
  257. FIN bloc1 ;
  258. * opti echo 1 ;
  259.  
  260. *-------------------------------
  261. * Post-traitement ...
  262. *-------------------------------
  263.  
  264. titre 'Champ de température final' ;
  265. si(graph) ; trac plaq2 (exco 'T' tn) ; finsi ;
  266.  
  267. Temp = exco 'T' tt 'T' ;
  268. T2 = redu plaq2 Temp ;
  269. * list T2 ;
  270.  
  271. RESI=maxi( ABS((T2 - T2_ana)/T2_ana));
  272. mess 'Solution analytique ' T2_ana ;
  273. mess 'Erreur relative ' (100*RESI) '%' ;
  274. SI (RESI &lt;EG 1E-4);
  275. ERRE 0;
  276. SINO;
  277. ERRE 5;
  278. FINSI;
  279.  
  280. fin;
  281.  
  282.  
  283.  
  284.  
  285.  
  286.  

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