Télécharger la plaquette Cast3M


fr en

Télécharger rayo_abs-2D-1.dgibi

Retour à la liste

Numérotation des lignes : 

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

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