Télécharger la plaquette Cast3M


fr en

Télécharger faceaface2.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 du rayonnement face à face                                      *
  10. * Comparaison avec une solution analytique en régime permanent.        *
  11. *                                                                      *
  12. * Calcul par 2 méthodes:                                               *
  13. *      1. algorithme iteratif                                          *
  14. *      2. utilisation de PASAPAS en transitoire jusqu'au permanent     *
  15. *                                                                      *
  16. *                                                                      *
  17. * Calcul 2D-plan.                                                      *
  18. *                                                                      *
  19. * on considère 3 plaques parallèles rayonnantes :                      *
  20. *                                                                      *
  21. *    -------- T3,e3                                                    *
  22. *                                                                      *
  23. *      e2sup           face à face                                     *
  24. *    -------- T2                                                       *
  25. *      e2inf           face à face                                     *
  26. *                                                                      *
  27. *    -------- T1,e1                                                    *
  28. *                                                                      *
  29. * La température des plaques supérieure et inférieure est imposée et   *
  30. * uniforme (respectivement 1000K et 2000K) et on détermine la tempéra- *
  31. * ture d'équilibre de la plaque intermédiaire.                         *
  32. * Le rayonnement entre les plaques 1 et 2 est modélisé par un rayonne- *
  33. * ment face à face, et de meme pour celui entre les plaques 2 et 3     *
  34. *                                                                      *
  35. * paramètre xquad  vrai : éléments quadratiques                        *
  36. *                  faux : éléments linéaires                           *
  37. ************************************************************************
  38.  
  39. graph = FAUX ;
  40.  
  41. xquad = vrai ;
  42.  
  43. si xquad ;
  44.   option echo 1 dime 2 elem qua8  ;
  45. sinon ;
  46.   option echo 1 dime 2 elem qua4  ;
  47. finsi;
  48.  
  49. *-------------------------------
  50. * Données physiques du problème
  51. *-------------------------------
  52.  
  53. T1 = 1000. ; T3 = 2000. ;
  54. e1 = 0.1 ; e2inf = 0.6 ; e2sup = 0.3 ; e3 = 0.9;
  55. lamb = 2.e3 ;
  56.  
  57. *-------------------------------
  58. * Solution analytique
  59. *-------------------------------
  60.  
  61. r1 = 1. - e1 ; r2i = 1. - e2inf ;
  62. r3 = 1. - e3 ; r2s = 1. - e2sup ;
  63. alpha1= (e1*e2inf) / (1.-(r1*r2i)) ;
  64. alpha3= (e3*e2sup) / (1.-(r3*r2s)) ;
  65. M1 = T1**4. ; M3 = T3**4. ;
  66. M2 = ((alpha1*M1) + (alpha3*M3)) / (alpha1+alpha3) ;
  67. T2_ana =  M2**(0.25) ;
  68. *LIST T2_ANA ;
  69.  
  70. *-------------------------------
  71. *   Maillage 
  72. *------------------------------- 
  73.  
  74. p  = 0. 0. ; q  = 1. 0. ; s  = droi 2 p q ;
  75. *s1inf = s tour (0.5 0.) 30. ;
  76. s1inf = s ;
  77. plaq1 = s1inf tran 1 (0. 0.001) ;
  78. s1sup = inve ( cote 3 plaq1) ;
  79.  
  80. s2inf = s1inf plus (0. 0.1) ;
  81. plaq2 = s2inf tran 1 (0. 0.001) ;
  82. s2sup = inve ( cote 3 plaq2) ;
  83.  
  84. rac12 = racc 0.1 s1sup s2inf ;
  85.  
  86. s3inf = s2inf plus (0. 0.1) ;
  87. plaq3 = s3inf tran 1 (0. 0.001) ;
  88. s3sup = inve ( cote 3 plaq3) ;
  89.  
  90. rac23 = racc 0.1 s2Sup s3inf ;
  91.  
  92. tout = plaq1 et plaq2 et plaq3 ;
  93.  
  94. si graph ; trac tout ; finsi ;
  95.  
  96. *------------------------------- 
  97. *  Conduction
  98. *------------------------------- 
  99.  
  100.   mcd = mode tout thermique ;
  101.   k   = mate mcd 'K' lamb 'RHO' 1. 'C' LAMB;
  102.   cnd = cond mcd k ;
  103.  
  104.   cte_sb = 5.673e-8 ;
  105.  
  106. *------------------------------- 
  107. *       Linéarisation du rayonnement :
  108. *    - définition de modèles de convection
  109. *    - creation des SEG2 définissant les relations entre les supports
  110. *      des champs définis sur s1sup et s2inf (pour le face à face)
  111. *------------------------------- 
  112.  
  113. * face à face
  114.  
  115.   mcv12 = mode rac12 thermique convection ;
  116.   mcv23 = mode rac23 thermique convection ;
  117.  
  118.   si xquad ;
  119.    x_s1sup = chan s1sup SEG2 ;
  120.    x_s2inf = chan s2inf SEG2 ;
  121.    x_rac12 = racc 0.1 x_s1sup x_s2inf ;
  122.  
  123.    rel12 = diff (x_rac12 chan ligne) (x_s1sup et x_s2inf) ;
  124.  
  125.    x_s2sup = chan s2sup SEG2 ;
  126.    x_s3inf = chan s3inf SEG2 ;
  127.    x_rac23 = racc 0.1 x_s2sup x_s3inf ;
  128.  
  129.    rel23 = diff (x_rac23 chan ligne) (x_s2sup et x_s3inf) ;
  130.  
  131.   sinon;
  132.    rel12 = diff (rac12 chan ligne) (s1sup et s2inf) ;
  133.  
  134.    rel23 = diff (rac23 chan ligne) (s2sup et s3inf) ;
  135.   finsi;
  136.  
  137.  
  138. *------------------------------- 
  139. *  Modeles de rayonnement et emissivites
  140. *------------------------------- 
  141.  
  142. * face a face
  143.   mr1 = mode s1sup thermique rayonnement FAC_A_FAC 
  144.       s1sup s2inf mcv12 rel12 cons 'FAC1';
  145.   emi1 = mate mr1 'EMIS' e1 ;
  146.  
  147.   mr2i = mode s2inf thermique rayonnement FAC_A_FAC
  148.       s1sup s2inf mcv12 rel12 cons 'FAC1';
  149.   emi2i= mate mr2i 'EMIS' e2inf ;
  150.  
  151.  
  152.   mr2S= mode s2sup thermique rayonnement FAC_A_FAC
  153.       s2sup s3inf mcv23 rel23 cons 'FAC2';
  154.   emi2S= mate mr2S 'EMIS' e2SUP ;
  155.  
  156.   mr3  = mode s3inf thermique rayonnement FAC_A_FAC
  157.       s2sup s3inf mcv23 rel23 cons 'FAC2';
  158.   emi3 = mate mr3  'EMIS' e3 ;
  159.  
  160.  
  161.  
  162. *------------------------------- 
  163. *  Conditions aux limites
  164. *------------------------------- 
  165.  
  166.     c1 = bloq plaq1   'T' ;
  167.     tim1 = depi c1 T1 ;
  168.  
  169.     c3 = bloq plaq3   'T' ;
  170.     tim3 = depi c3 T3 ;
  171.  
  172.     lr1 = prog   0 100 ;
  173.     lr2 = prog   1   1 ;
  174.  
  175.     ev1 = evol manu 't' lr1 'T' lr2 ;
  176.     cht1 = char 'TIMP' ev1 tIM1 ;
  177.     ev3 = evol manu 't' lr1 'T' lr2 ;
  178.     cht3 = char 'TIMP' ev3 tIM3 ;
  179.  
  180. *** Températures initiales ...
  181.  
  182.     T2 = 500. ;
  183.     chINI1 = manu chpo PLAQ1 1 'T' T1 nature diffus ;
  184.     chINI2 = manu chpo PLAQ2 1 'T' T2 nature diffus ;
  185.     chINI3 = manu chpo PLAQ3 1 'T' T3 nature diffus ;
  186.  
  187.  
  188. *----------------------------------------------------------
  189. * 1.  Première méthode de résolution : Algorithme itératif 
  190. *----------------------------------------------------------
  191.  
  192. *** Initialisation de la température ...
  193.  
  194.     tp = CHINI1 ET CHINI2 ET CHINI3 ;
  195.     nbpsup = nbno tout ;
  196.  
  197. *** Résolution (par itérations) ...
  198.  
  199. *   Coeff. de relaxation ...
  200.     alfa = 0.3 ;
  201.  
  202.     listemp = prog ;
  203.     listres = prog ;
  204.  
  205.     maxiter = 100 ;
  206.     critconv = 1.e-5 ;
  207. *   opti echo 1 ;
  208.     REPE bloc1 ;
  209.  
  210.        nbiter = &bloc1 ;
  211.  
  212. *   traitement du rayonnement:
  213. *                              face à face (1,2)
  214.  
  215.        t1sup = redu (exco 'T' tp 'T') s1sup;
  216.        t2inf = redu (exco 'T' tp 'T') s2inf;
  217.        h12 = HRAYO mcv12 mr1 emi1 t1sup mr2i emi2i t2inf rel12 ;
  218.  
  219.        cr12 = cond mcv12  h12      ;
  220.  
  221. *                              face à face (2,3)
  222.  
  223.        t2sup = redu (exco 'T' tp 'T') s2sup;
  224.        t3inf = redu (exco 'T' tp 'T') s3inf;
  225.        h23 = HRAYO mcv23 mr2S emi2S t2sup mr3 emi3 t3INF rel23 ;
  226.  
  227.        cr23 = cond mcv23  h23      ;
  228.  
  229.        crtot= cr12 et cr23  ;
  230.  
  231.        cndtot = crtot et cnd et c1 ET C3;
  232.        residu = (cndtot * tp)  ;
  233.        normres = ((xtx (exco 'Q' residu)) / nbpsup) ** 0.5 ;
  234. *      mess '   La norme du flux résiduel = ' normres ;
  235.        si((nbiter > 1) et (normres < critconv)) ;
  236.           mess 'Convergence à l itération ' (nbiter-1) ;
  237.           quitter bloc1 ;
  238.        finsi ;
  239.        si(nbiter > maxiter) ;
  240.           mess 'Erreur ! Pas de convergence après ' (nbiter-1)
  241.                ' itérations !' ;
  242.           erre 2 ;
  243.           quitter bloc1 ;
  244.        finsi ;
  245.        listres = listres et (prog normres) ;
  246.  
  247. *      mess '---------------------------------------' ;
  248. *      mess 'Itération N° ' &bloc1 ;
  249.  
  250.        tt = resou cndtot (tim1 et TIM3) ;
  251.        dt = exco 'T' (tt - tp) 'T' ;
  252.        normdt = ((xtx dt) / nbpsup) ** 0.5 ;
  253. *      mess '   La norme de delta t = ' normdt ;
  254.  
  255.        tn = (alfa * tt) + ((1.-alfa) * tp) ;
  256.        tp = tn ;
  257.  
  258.  
  259.     FIN bloc1 ;
  260. *   opti echo 1 ;
  261.  
  262. *------------------------------- 
  263. *  Post-traitement ...
  264. *------------------------------- 
  265.  
  266.     Temp = exco 'T' tt 'T' ;
  267.     T2A = redu plaq2   Temp ;
  268.     RESIA=maxi( ABS((T2A - T2_ana)/T2_ana));
  269.  
  270.  
  271. *------------------------------------------------------
  272. * 2.  Deuxieme méthode de résolution : appel a PASAPAS
  273. *------------------------------------------------------
  274.  
  275.  
  276. *** Préparation de la table pour PASAPAS ...
  277.  
  278.     tabnl = table ;
  279.  
  280.     tabnl . MODELE = MCD ET  MR1 ET MR2I ET MR2S ET MR3 ;
  281.  
  282.  
  283.     tabnl . CARACTERISTIQUES = K ET EMI1 ET EMI2I ET EMI2S ET EMI3;
  284.     tabnl . BLOCAGES_THERMIQUES = C1 ET C3 ;
  285.     tabnl . CHARGEMENT = cht1 et cht3 ;
  286.  
  287.     tabnl . TEMPS_CALCULES = prog 0. pas  0.1  1. 
  288.                                      pas  0.2  4. ;
  289.  
  290.     tabnl . TEMPERATURES = table ;
  291.     tabnl . TEMPERATURES . 0 = CHINI1 ET CHINI2 ET CHINI3 ;
  292.  
  293.  
  294.     tabnl . 'PROCEDURE_THERMIQUE' = 'DUPONT' ;
  295.     tabnl . 'CTE_STEFAN_BOLTZMANN' = cte_sb ;
  296.  
  297.     tabnl . 'RELAXATION_THETA' = 1. ;
  298.  
  299. *** Appel à PASAPAS ...
  300.  
  301.     pasapas tabnl ;
  302.  
  303.  
  304. *** Petit post-traitement ...
  305.  
  306.     listt = prog ;
  307.     listtp3 = prog ;
  308.     listtp3r = prog ;
  309.     nbpas = dime (tabnl . TEMPS) ;
  310.     repeter surpas nbpas ;
  311.        lindice = &surpas - 1 ;
  312.        listt = listt et (prog (tabnl . TEMPS . lindice)) ;
  313.        valtp3 = MAXI (REDU tabnl . TEMPERATURES . lindice PLAQ2) ;
  314.        listtp3 = listtp3 et (prog valtp3) ;
  315.     fin surpas ;
  316.  
  317.     titr 'Evolution de la temperature de la plaque2';
  318.     evtp3 = evol manu 't' listt 'T' listtp3 ;
  319.  
  320.     si(graph) ;
  321.        dess evtp3  ;
  322.     finsi ;
  323.  
  324.     T2B = redu tabnl . TEMPERATURES . (nbpas - 1) PLAQ2  ;
  325.     RESIB=maxi( ABS((T2B - T2_ana)/T2_ana));
  326.  
  327. *  comparaison a la solution analytique 
  328.  
  329.     mess 'Solution analytique ' T2_ana ;
  330.  
  331.     mess 'Solution Algorithme ' (MAXI T2A) ;
  332.     mess 'Erreur relative    ' (100*RESIA) '%' ;
  333.  
  334.     mess 'Solution PASAPAS    ' (MAXI T2B) ;
  335.     mess 'Erreur relative    ' (100*RESIB) '%' ;
  336.  
  337.     SI ((RESIA &lt;EG 1.E-4)et(RESIB &lt;EG 1.E-4)) ;
  338.        ERRE  0;
  339.     SINO;
  340.        ERRE  5;
  341.     FINSI;
  342.  
  343.  
  344.     fin ;
  345.  
  346.  
  347.  
  348.  
  349.  
  350.  

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