Télécharger rayoh-2D.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : rayoh-2D.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 ou avec un milieu infini. *
  11. * Comparaison avec une solution analytique en régime permanent. *
  12. * *
  13. * Calcul 2D-plan. *
  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. ************************************************************************
  36.  
  37. graph = faux ;
  38.  
  39. xquad = vrai ;
  40.  
  41. si xquad ;
  42. option echo 1 dime 2 elem qua8 ;
  43. sinon ;
  44. option echo 1 dime 2 elem qua4 ;
  45. finsi;
  46.  
  47. *-------------------------------
  48. * Données physiques du problème
  49. *-------------------------------
  50.  
  51. T1 = 1000. ; T3 = 2000. ;
  52. e1 = 0.1 ; e2inf = 0.6 ; e2sup = 0.3 ; e3 = 0.9 ;
  53. lamb = 2.e3 ;
  54.  
  55. *-------------------------------
  56. * Solution analytique
  57. *-------------------------------
  58.  
  59. r1 = 1. - e1 ; r2i = 1. - e2inf ;
  60. r3 = 1. - e3 ; r2s = 1. - e2sup ;
  61. alpha1= (e1*e2inf) / (1.-(r1*r2i)) ;
  62. alpha3= (e3*e2sup) / (1.-(r3*r2s)) ;
  63. M1 = T1**4. ; M3 = T3**4. ;
  64. M2 = ((alpha1*M1) + (alpha3*M3)) / (alpha1+alpha3) ;
  65. T2_ana = M2**(0.25) ;
  66.  
  67. *-------------------------------
  68. * Maillage
  69. *-------------------------------
  70.  
  71. p = 0. 0. ; q = 1. 0. ; s = droi 2 p q ;
  72. *s1inf = s tour (0.5 0.) 30. ;
  73. s1inf = s ;
  74. plaq1 = s1inf tran 1 (0. 0.001) ;
  75. s1sup = inve ( cote 3 plaq1) ;
  76.  
  77. s2inf = s1inf plus (0. 0.1) ;
  78. plaq2 = s2inf tran 1 (0. 0.001) ;
  79. s2sup = inve ( cote 3 plaq2) ;
  80.  
  81. rac12 = racc 0.1 s1sup s2inf ;
  82.  
  83. tout = plaq1 et plaq2 ;
  84.  
  85. si graph ; trac tout ; finsi ;
  86.  
  87. *-------------------------------
  88. * Conduction
  89. *-------------------------------
  90.  
  91. mcd = mode tout thermique ;
  92. k = mate mcd 'K' lamb ;
  93. cnd = cond mcd k ;
  94.  
  95.  
  96. *-------------------------------
  97. * Linéarisation du rayonnement :
  98. * - définition de modèles de convection
  99. * - creation des SEG2 définissant les relations entre les supports
  100. * des champs définis sur s1sup et s2inf (pour le face à face)
  101. *-------------------------------
  102.  
  103. * face à face
  104.  
  105. mcv12 = mode sac1² tlermique convection ;
  106.  
  107. si xquad ;
  108. x_s1sup = chan s1sup SEG2 ;
  109. x_s2inf = chan s2inf SEG2 ;
  110. x_rac12 = racc 0.1 x_s1sup x_s2inf ;
  111.  
  112. rel12 = diff (x_rac12 chan ligne) (x_s1sup et x_s2inf) ;
  113. sinon;
  114. rel12 = diff (rac12 chan ligne) (s1sup et s2inf) ;
  115. finsi;
  116.  
  117. * milieu infini
  118.  
  119. mcv2s = mode s2sup thermique convection ;
  120.  
  121.  
  122.  
  123. *-------------------------------
  124. * Mod�les de rayonnement et �missivit�s
  125. *-------------------------------
  126.  
  127. * face � face
  128. mr1 = mode s1sup thermique rayonnement FAC_A_FAC
  129. s1sup s2inf mcv12 rel12 cons 'FAC1';
  130. emi1 = mate mr1 'EMIS' e1 ;
  131.  
  132. mr2i = mode s2inf thermique rayonnement FAC_A_FAC
  133. s1sup s2inf mcv12 rel12 cons 'FAC1';
  134.  
  135. emi2i= mate mr2i 'EMIS' e2inf ;
  136.  
  137. * milieu infini
  138.  
  139. mr2s = mode s2sup thermique rayonnement INFINI CONS 'INF1';
  140. emi2s= mate mr2s 'EMIS' e2sup 'E_IN' e3 ;
  141.  
  142. t3inf = manu chpo s2sup 1 'T' T3 ;
  143.  
  144.  
  145. *-------------------------------
  146. * Conditions aux limites
  147. *-------------------------------
  148.  
  149. c1 = bloq plaq1 'T' ;
  150. tim1 = depi c1 T1 ;
  151.  
  152. *-------------------------------
  153. * Algorithme de résolution
  154. *-------------------------------
  155.  
  156. *** Initialisation de la température ...
  157.  
  158. tp = manu chpo tout 1 'T' T3 natu 'DIFFUS';
  159. nbpsup = nbno tout ;
  160.  
  161. *** Résolution (par itérations) ...
  162.  
  163. * Coeff. de relaxation ...
  164. alfa = 0.3 ;
  165.  
  166. listemp = prog ;
  167. listres = prog ;
  168.  
  169. maxiter = 100 ;
  170. critconv = 1.e-5 ;
  171. * opti echo 1 ;
  172. REPE bloc1 ;
  173.  
  174. nbiter = &bloc1 ;
  175.  
  176. * traitement du rayonnement:
  177. * face à face (1,2)
  178.  
  179. t1sup = redu (exco 'T' tp 'T') s1sup;
  180. t2inf = redu (exco 'T' tp 'T') s2inf;
  181. h12 = HRAYO mcv12 mr1 emi1 t1sup mr2i emi2i t2inf rel12 ;
  182.  
  183. cr12 = cond mcv12 h12 ;
  184.  
  185. * milieu infini (2,3)
  186.  
  187. t2sup = redu (exco 'T' tp 'T') s2sup;
  188. h23 = HRAYO mcv2s mr2s emi2s t2sup t3inf ;
  189.  
  190.  
  191. * mr23 = MATE mcv2s 'H' h23 ;
  192. cr23 = cond mcv2s h23 ;
  193. f3 = conv mcv2s h23 t3inf ;
  194.  
  195. crtot= cr12 et cr23 ;
  196.  
  197. cndtot = crtot et cnd et c1 ;
  198. residu = (cndtot * tp) - f3 ;
  199. normres = ((xtx (exco 'Q' residu)) / nbpsup) ** 0.5 ;
  200. * mess ' La norme du flux résiduel = ' normres ;
  201. si((nbiter > 1) et (normres < critconv)) ;
  202. mess 'Convergence à l itération ' (nbiter-1) ;
  203. quitter bloc1 ;
  204. finsi ;
  205. si(nbiter > maxiter) ;
  206. mess 'Erreur ! Pas de convergence après ' (nbiter-1)
  207. ' itérations !' ;
  208. erre 2 ;
  209. quitter bloc1 ;
  210. finsi ;
  211. listres = listres et (prog normres) ;
  212.  
  213. * mess '---------------------------------------' ;
  214. * mess 'Itération N° ' &bloc1 ;
  215.  
  216. tt = resou cndtot (tim1 et f3) ;
  217. dt = exco 'T' (tt - tp) 'T' ;
  218. normdt = ((xtx dt) / nbpsup) ** 0.5 ;
  219. * mess ' La norme de delta t = ' normdt ;
  220.  
  221. tn = (alfa * tt) + ((1.-alfa) * tp) ;
  222. tp = tn ;
  223.  
  224.  
  225. FIN bloc1 ;
  226. * opti echo 1 ;
  227.  
  228. *-------------------------------
  229. * Post-traitement ...
  230. *-------------------------------
  231.  
  232. titre 'Champ de température final' ;
  233. si(graph) ; trac plaq2 (exco 'T' tn) ; finsi ;
  234.  
  235. Temp = exco 'T' tt 'T' ;
  236. T2 = redu plaq2 Temp ;
  237.  
  238. RESI=maxi( ABS((T2 - T2_ana)/T2_ana));
  239. mess 'Solution analytique ' T2_ana ;
  240. mess 'Erreur relative ' (100*RESI) '%' ;
  241. SI (RESI &lt;EG 1E-4);
  242. ERRE 0;
  243. SINO;
  244. ERRE 5;
  245. FINSI;
  246.  
  247. fin;
  248.  
  249.  
  250.  
  251.  
  252.  
  253.  
  254.  
  255.  
  256.  
  257.  
  258.  
  259.  
  260.  
  261.  
  262.  
  263.  
  264.  

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