Télécharger la plaquette Cast3M


fr en

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

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  rayo-axi-2.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4.  
  5. *  pour calcul complet mettre complet à : vrai;
  6. complet = faux;
  7.  
  8. ************************************************************************
  9. * Rayonnement thermique en milieu transparent:                         *
  10. *                                                                      *
  11. * vérification du bon fonctionnement de l'opérateur FFOR dans le cas   *
  12. * général (traitement des parties cachées) sur un cas thermiquement    *
  13. * simple.                                                              *
  14. *                                                                      *
  15. * Test 2D-axisymétrique sur un problème à symétrie sphérique 1D:       *
  16. *                                         ------------------           *
  17. * on considère 3 sphères concentriques de rayons 25mm, 50mm et 100mm   *
  18. * conductrices et d'emissivité 1. La température des sphères intérieure*
  19. * et extérieure est uniforme (respectivement 1000K et 2000K) et on     *
  20. * détermine la température d'équilibre de la sphère intermédiaire.     *
  21. * (conventions: indices 1,2,3 pour les sphères de rayon 25,50,100mm)   *
  22. *                                                                      *
  23. * Jeu de données identique au cas 2D-plan :rayo-2D-3.dgibi             *
  24. * hormis la valeur des 2 paramètres:                                   *
  25. *            'xrev'  vrai : géométrie axisymétrique                    *
  26. *                    faux :           2D plane                         *
  27. *            'x2cav' vrai : on considère les 2 cavités séparément      *
  28. *                    faux :              l'ensemble des 2 cavités      *
  29. *                                                                      *
  30. ************************************************************************
  31.  
  32. xrev = vrai ; x2cav = faux ;
  33. *  attention ne fonctionne peut etre pas si  x2cav= vrai et xrev=vrai
  34. * il faudrait refaire le model sans "symetrie .."
  35. graph = faux ;
  36.  
  37. si xrev ;
  38.  option echo 0 dime 2 elem qua4 mode axis ;
  39. sinon;
  40.  option echo 0 dime 2 elem qua4  ;
  41. finsi;
  42.  
  43. *-------------------------------
  44. * Données physiques du problème
  45. *-------------------------------
  46.  
  47. r1 = 25e-3 ; r2 = 50e-3 ; r3 = 100e-3 ;
  48. T1 = 1000. ; T3 = 2000. ;
  49. emis = 1. ; lamb = 20. ;
  50.  
  51. *-------------------------------                                        
  52. * Solution analytique
  53. *-------------------------------
  54.  
  55. si xrev ;
  56.  f21 = (r1*r1)/(r2*r2) ;
  57. sinon;
  58.  f21 = r1/r2 ;
  59. finsi;
  60. f22 = 1. - f21 ;
  61. M1 = T1**4. ; M3 = T3**4. ;
  62. M2 = ( (f21*M1) + M3 )/(2.-f22) ;
  63. T2_ana =  M2**(0.25) ;
  64.  
  65. *-------------------------------
  66. *   Maillage 
  67. *------------------------------- 
  68.  
  69. C  =  0. 0. ;
  70.  
  71. * sphère 1
  72.  
  73. p1_sud =  0. ( -1. * r1 )  ;
  74. p1_equa =  r1  0. ;
  75. p1_nord =  0. r1 ;
  76. d1 = r1/3 ;
  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. c1_tot = c1_inf et c1_sup ;
  80. c1_int = c1_tot homo C 0.999 ;
  81. c1_ext = c1_tot homo C 1.001 ;
  82. sphere1 = regl c1_int c1_ext ;
  83.  
  84. * sphère 2
  85.  
  86. p2_sud =  0. ( -1. * r2 )  ;
  87. p2_equa =  r2  0. ;
  88. p2_nord =  0. r2 ;
  89. d2 = r2/3 ;
  90. c2_inf = cerc p2_sud C p2_equa DINI d2 DFIN d2;
  91. c2_sup = cerc p2_equa C p2_nord DINI d2 DFIN d2;
  92. c2_tot = c2_inf et c2_sup ;
  93. c2_int = c2_tot homo C 0.999 ;
  94. c2_ext = c2_tot homo C 1.001 ;
  95. sphere2 = regl c2_int c2_ext ;
  96.  
  97. * sphère 3
  98.  
  99. p3_sud =  0. ( -1. * r3 )  ;
  100. p3_equa =  r3  0. ;
  101. p3_nord =  0. r3 ;
  102. si complet;
  103. d3 = r3/3 ;
  104. sinon;
  105. d3 = r3 / 1.5;
  106. finsi;
  107. c3_inf = cerc p3_sud C p3_equa DINI d3 DFIN d3;
  108. c3_sup = cerc p3_equa C p3_nord DINI d3 DFIN d3;
  109. c3_tot = c3_inf et c3_sup ;
  110. c3_int = c3_tot homo C 0.999 ;
  111. c3_ext = c3_tot homo C 1.001 ;
  112. sphere3 = regl c3_int c3_ext ;
  113.  
  114. * on controle l'orientation des elements pour les cavités rayonnantes
  115.  
  116. cavi_12 = (inve c1_ext) et c2_int ;
  117. cavi_23 = (inve c2_ext) et c3_int ;
  118. cavi_tot = cavi_12 et cavi_23 ;
  119.  
  120. tout = sphere1 et sphere2 et sphere3 ;
  121.  
  122. si graph ; 
  123. trac tout  ;
  124. finsi ; 
  125.  
  126. *------------------------------- 
  127. *  Conduction
  128. *------------------------------- 
  129.  
  130. mcd = mode tout thermique ;
  131. k = mate mcd 'K' lamb ;
  132. cnd = cond mcd k ;
  133.  
  134. *------------------------------- 
  135. *  Modèles de rayonnement
  136. *------------------------------- 
  137.  
  138. si x2cav;
  139.   mr12 = mode cavi_12  thermique rayonnement 'CAVITE' 
  140.   SYMETRIE C (0. 1.)  'CONS'  'CAV1' ;
  141.   e12 = mate mr12 'EMIS' emis ;
  142.  
  143.   mr23 = mode cavi_23 thermique rayonnement 'CAVITE' 
  144.     SYMETRIE C (0. 1.)  'CONS'  'CAV1' ;
  145.   e23 = mate mr23 'EMIS' emis ;
  146. sinon;
  147.   mrtot= mode cavi_tot  thermique rayonnement  'CAVITE' 
  148.  SYMETRIE  C (0. 1.)  ; 
  149.   etot = mate mrtot 'EMIS' emis ;
  150. finsi;
  151.  
  152. *------------------------------- 
  153. *  Facteurs de forme et matrices de rayonnement
  154. *------------------------------- 
  155.  
  156. si x2cav ;
  157. *opti 'IMPI' 1 ;
  158.  si xrev ;
  159.   ff12 = ffor mr12  e12;
  160.   ff23 = ffor mr23  e23;
  161.  sinon;
  162.   ff12 = ffor mr12 e12 ;
  163.   ff23 = ffor mr23 e23 ;
  164.  finsi;
  165. *opti 'IMPI' 0 ;
  166.  
  167.  cham12 = raye mr12 ff12 e12 ;
  168.  cham23 = raye mr23 ff23 e23 ;
  169.  
  170. sinon;
  171. *opti 'IMPI' 1 ;
  172.  si xrev ;
  173.   fftot= ffor mrtot  etot;
  174.  sinon;
  175.   fftot = ffor mrtot etot ;
  176.  finsi;
  177. *opti 'IMPI' 0 ;
  178.  
  179.  chamtot = raye mrtot fftot etot ;
  180.  
  181. finsi;
  182.  
  183. *------------------------------- 
  184. *  Conditions aux limites
  185. *------------------------------- 
  186.  
  187.     c1 = bloq sphere1 'T' ;
  188.     tim1 = depi c1 T1 ;
  189.     c3 = bloq sphere3 'T' ;
  190.     tim3 = depi c3 T3 ;
  191.  
  192. *------------------------------- 
  193. *  Algorithme de résolution
  194. *------------------------------- 
  195.  
  196. *** Initialisation de la température ...
  197.  
  198.     tp = manu chpo tout 1 'T' T3 natu 'DIFFUS';
  199.     nbpsup = nbno tout ;
  200.  
  201. *** Résolution (par itérations) ...
  202.  
  203. *   Coeff. de relaxation ...
  204.     alfa = 0.3 ;
  205.  
  206.     listemp = prog ;
  207.     listres = prog ;
  208.  
  209.     maxiter = 100 ;
  210. si complet;
  211.     critconv = 1.e-5 ;
  212. sinon;
  213.    critconv = 1.;
  214. finsi;
  215.  
  216. *   opti echo 1 ;
  217.     REPE bloc1 ;
  218.  
  219.        nbiter = &bloc1 ;
  220.  
  221. *   traitement du rayonnement
  222.  
  223.     si x2cav ;
  224.        t_c12 = redu (exco 'T' tp 'T') cavi_12 ;
  225.        te_c12 = chan 'CHAM' t_c12   mr12 'GRAVITE' ;
  226.        cr12 = rayn mr12  cham12 te_c12 ;
  227.  
  228.        t_c23 = redu (exco 'T' tp 'T') cavi_23 ;
  229.        te_c23 = chan 'CHAM' t_c23   mr23 'GRAVITE' ;
  230.        cr23 = rayn mr23  cham23 te_c23 ;
  231.        crtot = cr12 et cr23 ;
  232.  
  233.     sinon;
  234.        t_tot = redu (exco 'T' tp 'T') cavi_tot;
  235.        te_tot = chan 'CHAM' t_tot mrtot  'GRAVITE' ;
  236.        crtot= rayn mrtot chamtot te_tot ;
  237.  
  238.     finsi;
  239.  
  240.        cndtot = crtot et cnd et c1 et c3 ;
  241.        residu = cndtot * tp ;
  242.        normres = ((xtx (exco 'Q' residu)) / nbpsup) ** 0.5 ;
  243. *      mess '   La norme du flux résiduel = ' normres ;
  244.        si((nbiter > 1) et (normres < critconv)) ;
  245.           mess 'Convergence à l itération ' (nbiter-1) ;
  246.           quitter bloc1 ;
  247.        finsi ;
  248.        si(nbiter > maxiter) ;
  249.           mess 'Erreur ! Pas de convergence après ' (nbiter-1)
  250.                ' itérations !' ;
  251.           erre 2 ;
  252.           quitter bloc1 ;
  253.        finsi ;
  254.        listres = listres et (prog normres) ;
  255.  
  256. *      mess '---------------------------------------' ;
  257. *      mess 'Itération N° ' &bloc1 ;
  258.  
  259.        tt = resou cndtot (tim1 et tim3) ;
  260.        dt = exco 'T' (tt - tp) 'T' ;
  261.        normdt = ((xtx dt) / nbpsup) ** 0.5 ;
  262. *      mess '   La norme de delta t = ' normdt ;
  263.  
  264.        tn = (alfa * tt) + ((1.-alfa) * tp) ;
  265.        tp = tn ;
  266.  
  267.  
  268.     FIN bloc1 ;
  269. *   opti echo 1 ;
  270.  
  271. *------------------------------- 
  272. *  Post-traitement ...
  273. *------------------------------- 
  274.  
  275.     titre 'Champ de température final' ;
  276.     si(graph) ; trac sphere2 (exco 'T' tn) ; finsi ;
  277.  
  278.     Temp = exco 'T' tt 'T' ;
  279.     T2 = redu sphere2 Temp ;
  280.  
  281.     RESI=maxi( ABS((T2 - T2_ana)/T2_ana));
  282.     mess 'Problème à symétrie sphérique';
  283.     mess 'Solution analytique ' T2_ana ;
  284.     mess 'Erreur relative ' (100*RESI) '%' ;
  285.     SI (RESI &lt;EG 2E-3);
  286.        ERRE  0;
  287.     SINO;
  288.        ERRE  5;
  289.     FINSI;
  290.  
  291. fin;
  292.  
  293.  
  294.  
  295.  
  296.  
  297.  

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