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

Retour à la liste

Numérotation des lignes :

  1. * fichier : rayo-3D-2.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4.  
  5. complet = faux;
  6. *
  7. *************************************************
  8. * pour calcul complet mettre complet à : vrai;
  9. ***************************************************
  10. ************************************************************************
  11. * *
  12. * Rayonnement thermique en milieu transparent: *
  13. * *
  14. * vérification du bon fonctionnement de l'opérateur FFOR dans le cas *
  15. * général (traitement des parties cachées) avec plan de symétrie sur *
  16. * un cas thermiquement simple. *
  17. * *
  18. * Test 3D sur un problème à symétrie sphérique 1D: *
  19. * ------------------ *
  20. * on considère 3 sphères concentriques de rayons 25mm, 50mm et 100mm *
  21. * conductrices et d'emissivité 1. La température des sphères intérieure*
  22. * et extérieure est uniforme (respectivement 1000K et 2000K) et on *
  23. * détermine la température d'équilibre de la sphère intermédiaire. *
  24. * (conventions: indices 1,2,3 pour les sphères de rayons 25,50,100mm) *
  25. * *
  26. * *
  27. * Remarques: *
  28. * avec la densité choisie, la cavité comprend 96 éléments *
  29. * ce problème est également traité en 2D-axis: cf. rayo-axi-2.dgibi *
  30. ************************************************************************
  31.  
  32. x2cav = faux ;
  33.  
  34. graph = faux ;
  35.  
  36. option echo 1 dime 3 elem tri3 ;
  37.  
  38. *-------------------------------
  39. * Données physiques du problème
  40. *-------------------------------
  41.  
  42. r1 = 25e-3 ; r2 = 50e-3 ; r3 = 100e-3 ;
  43. T1 = 1000. ; T3 = 2000. ;
  44. emis = 1. ; lamb = 20. ;
  45.  
  46. *-------------------------------
  47. * Solution analytique
  48. *-------------------------------
  49.  
  50. f21 = (r1*r1)/(r2*r2) ;
  51. f22 = 1. - f21 ;
  52. M1 = T1**4. ; M3 = T3**4. ;
  53. M2 = ( (f21*M1) + M3 )/(2.-f22) ;
  54. T2_ana = M2**(0.25) ;
  55.  
  56. *-------------------------------
  57. * Maillage
  58. *-------------------------------
  59.  
  60. C = 0. 0. 0. ;
  61. O1 = 0. 1.E-6 0. ;
  62. O2 = 0. -1.E-6 0. ;
  63. oeil = -10. -10. 10. ;
  64.  
  65. * sphère 1
  66.  
  67. p1_sud = 0. 0. (-1.*r1) ;
  68. p1_equa = r1 0. 0. ;
  69. p1_nord = 0. 0. r1 ;
  70.  
  71. * densité
  72. aa = 2.;
  73. *si complet ; aa = 2. ; sinon; aa = 1.;finsi;
  74. d1 = r1/aa ;
  75.  
  76. * cercle
  77. c1_inf = cerc p1_sud C p1_equa DINI d1 DFIN d1;
  78. c1_sup = cerc p1_equa C p1_nord DINI d1 DFIN d1;
  79.  
  80. c1_droi = c1_inf et c1_sup ;
  81. c1_gauc = inve(SYME c1_droi 'DROIT' p1_sud p1_nord) ;
  82. c1_tot = c1_droi et c1_gauc;
  83. elim 1.e-6 c1_tot ;
  84. * sphere de reference
  85. s1_ref = surf c1_tot 'SPHERIQUE' O2 ;
  86.  
  87. s1_int = s1_ref homo C 0.99 ;
  88. s1_ext = s1_ref homo C 1.01 ;
  89.  
  90. opti elem pri6 ;
  91. sphere1 = volu s1_int s1_ext ;
  92.  
  93. * spheres 2 et 3
  94.  
  95. s2_int = s1_ref homo C (0.999*(r2/r1)) ;
  96. s2_ext = s1_ref homo C (1.001*(r2/r1)) ;
  97. sphere2 = volu s2_int s2_ext ;
  98.  
  99. s3_int = s1_ref homo C (0.999*(r3/r1)) ;
  100. s3_ext = s1_ref homo C (1.001*(r3/r1)) ;
  101. sphere3 = volu s3_int s3_ext ;
  102.  
  103. * on controle l'orientation des elements pour les cavités rayonnantes
  104.  
  105. cavi_12 = (inve s1_ext) et s2_int ;
  106. cavi_23 = (inve s2_ext) et s3_int ;
  107. cavi_tot = cavi_12 et cavi_23 ;
  108.  
  109. tout = sphere1 et sphere2 et sphere3 ;
  110.  
  111. si graph ;
  112. trac tout oeil cache ;
  113. finsi ;
  114.  
  115.  
  116. *-------------------------------
  117. * Conduction
  118. *-------------------------------
  119.  
  120. mcd = mode tout thermique ;
  121. k = mate mcd 'K' lamb ;
  122. cnd = cond mcd k ;
  123.  
  124. *-------------------------------
  125. * Modèles de rayonnement
  126. *-------------------------------
  127. si x2cav;
  128. mr12 = mode cavi_12 thermique rayonnement 'CAVITE' 'SYMETRIE'
  129. C (1. 0. 0.) (0. 0. 1.) 'CONS' 'CAV1' ;
  130. e12 = mate mr12 'EMIS' emis ;
  131.  
  132. mr23 = mode cavi_23 thermique rayonnement 'CAVITE' 'SYMETRIE'
  133. C (1. 0. 0.) (0. 0. 1.) 'CONS' 'CAV1' ;
  134. e23 = mate mr23 'EMIS' emis ;
  135. sinon;
  136. mrtot= mode cavi_tot thermique rayonnement 'CAVITE' 'SYMETRIE'
  137. C (1. 0. 0.) (0. 0. 1.) ;
  138. etot = mate mrtot 'EMIS' emis ;
  139. finsi;
  140.  
  141.  
  142. *-------------------------------
  143. * Facteurs de forme et matrices de rayonnement
  144. *-------------------------------
  145.  
  146. * le plan y = 0. est le plan de symétrie
  147. si x2cav ;
  148. *opti 'IMPI' 1 ;
  149. ff12 = ffor mr12 e12 ;
  150. ff23 = ffor mr23 e23 ;
  151. *opti 'IMPI' 0 ;
  152.  
  153. cham12 = raye mr12 ff12 e12 ;
  154. cham23 = raye mr23 ff23 e23 ;
  155.  
  156. sinon;
  157. *opti 'IMPI' 1 ;
  158. fftot = ffor mrtot etot;
  159. *opti 'IMPI' 0 ;
  160.  
  161. chamtot = raye mrtot fftot etot ;
  162.  
  163. finsi;
  164.  
  165.  
  166.  
  167.  
  168. *-------------------------------
  169. * Conditions aux limites
  170. *-------------------------------
  171.  
  172. c1 = bloq sphere1 'T' ;
  173. tim1 = depi c1 T1 ;
  174. c3 = bloq sphere3 'T' ;
  175. tim3 = depi c3 T3 ;
  176.  
  177. *-------------------------------
  178. * Algorithme de résolution
  179. *-------------------------------
  180.  
  181. *** Initialisation de la température ...
  182.  
  183. tp = manu chpo tout 1 'T' T3 natu 'DIFFUS';
  184. nbpsup = nbno tout ;
  185.  
  186. *** Résolution (par itérations) ...
  187.  
  188. * Coeff. de relaxation ...
  189. alfa = 0.3 ;
  190.  
  191. listemp = prog ;
  192. listres = prog ;
  193.  
  194. maxiter = 100 ;
  195. si complet;
  196. critconv = 1.e-5 ;
  197. sinon;
  198. critconv = 10. ;
  199. finsi;
  200. * opti echo 1 ;
  201. REPE bloc1 ;
  202.  
  203. nbiter = &bloc1 ;
  204.  
  205. * traitement du rayonnement
  206.  
  207.  
  208. t_tot = redu (exco 'T' tp 'T') cavi_tot;
  209. te_tot = chan 'CHAM' t_tot mrtot 'GRAVITE' ;
  210. crtot= rayn mrtot chamtot te_tot ;
  211.  
  212.  
  213.  
  214. cndtot = crtot et cnd et c1 et c3 ;
  215. residu = cndtot * tp ;
  216. normres = ((xtx (exco 'Q' residu)) / nbpsup) ** 0.5 ;
  217. * mess ' La norme du flux résiduel = ' normres ;
  218. si((nbiter > 1) et (normres < critconv)) ;
  219. mess 'Convergence à l itération ' (nbiter-1) ;
  220. quitter bloc1 ;
  221. finsi ;
  222. si(nbiter > maxiter) ;
  223. mess 'Erreur ! Pas de convergence après ' (nbiter-1)
  224. ' itérations !' ;
  225. erre 2 ;
  226. quitter bloc1 ;
  227. finsi ;
  228. listres = listres et (prog normres) ;
  229.  
  230. * mess '---------------------------------------' ;
  231. * mess 'Itération N° ' &bloc1 ;
  232.  
  233. tt = resou cndtot (tim1 et tim3) ;
  234. dt = exco 'T' (tt - tp) 'T' ;
  235. normdt = ((xtx dt) / nbpsup) ** 0.5 ;
  236. * mess ' La norme de delta t = ' normdt ;
  237.  
  238. tn = (alfa * tt) + (|spcn style="color: #000000; font-weight: bold;">(1.-alfa) * tp) ;
  239. tp = tn ;
  240.  
  241.  
  242. FIN bloc1 ;
  243. * opti echo 1 ;
  244.  
  245. *-------------------------------
  246. * Post-traitement ...
  247. *-------------------------------
  248.  
  249. titre 'Champ de température final' ;
  250. si(graph) ; trac sphere2 (exco 'T' tn) ; finsi ;
  251.  
  252. Temp = exco 'T' tt 'T' ;
  253. T2 = redu sphere2 Temp ;
  254. mess 'resultat trouve' (maxi t2);
  255. RESI=maxi( ABS((T2 - T2_ana)/T2_ana));
  256. mess 'Problème à symétrie sphérique';
  257. mess 'Solution analytique ' T2_ana ;
  258. mess 'Erreur relative ' (100*RESI) '%' ;
  259. si complet; valres = 1.2e-2; sinon; valres = 3.e-2; finsi;
  260. SI (RESI &lt;EG valres);
  261. ERRE 0;
  262. SINO;
  263. ERRE 5;
  264. FINSI;
  265.  
  266. fin;
  267.  
  268.  
  269.  
  270.  
  271.  
  272.  

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