Télécharger la plaquette Cast3M


fr en

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

Retour à la liste

Numérotation des lignes : 

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

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