Télécharger la plaquette Cast3M


fr en

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

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  rayo-2D-5.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4.  
  5. *************************************************************
  6. *                                                           *
  7. *   Calcul d'une plaque infinie (de largeur 2L) avec tempé- *
  8. *   rature imposée au milieu et soumise au rayonnement.     *
  9. *   Il s'agit de trouver la température au bord en régime   *
  10. *   permanent.                                              *
  11. *                                                           *
  12. *   Néanmoins nous cherchons la solution en faisant un      * 
  13. *   calcul transitoire avec PASAPAS pour tester celle-ci    *
  14. *   en configuration température imposée + rayonnement vers *
  15. *   l'infini.                                               *
  16. *                                                           *
  17. *   Modélisation plane.                                     *
  18. *                                                           *
  19. *   Auteurs : Michel Bulik & Nadia Coulon                   *
  20. *                                                           *
  21. *   Date : Avril 1997                                       *
  22. *                                                           *
  23. *   on introduit une emissivite a l infini differente de 1  *    
  24. *      et une temperature a l'infini differente de 0K       *    
  25. *************************************************************
  26. *
  27. *** Options ...
  28.  
  29.     opti dime 2 elem seg2 echo 1 ;
  30.  
  31. *** Paramètres ...
  32.  
  33.     L = 1. ;
  34.     ep = L / 10 ;
  35.  
  36.     graph = faux ;
  37.  
  38.    Tempini = 1000. ;
  39.    Tempext = 2. * Tempini ;
  40.  
  41.     cte_sb = 5.673e-8 ;
  42.     Gamma = 1. ;
  43.     si ( neg Gamma 0. ) ;
  44.        lambda = (cte_sb * Tempini * Tempini * Tempini * L) / Gamma ;
  45.        mess 'Lambda = ' lambda ;
  46.        valemis = 0.8 ;
  47.         e0  = 0.5 ; 
  48.  
  49.     sinon ;
  50.        lambda = 1. ;
  51.        mess 'Lambda choisie arbitrairement = ' lambda ;
  52.        valemis = 0. ;
  53.        e0  = 1. ;
  54.     finsi ;
  55.  
  56. *** Points ...
  57.  
  58.     dens 0.1 ;
  59.     p1 = 0 0 ;
  60.     p2 = 0 ep ;
  61.  
  62.     vechoriz = L 0 ;
  63.  
  64. *** Lignes ...
  65.  
  66.     li1 = p1 d 1 p2 ;
  67.  
  68. *** Surface ...
  69.  
  70.     opti elem qua4 ;
  71.     su1 = li1 tran vechoriz dini 0.1 dfin 0.01 ;
  72.  
  73.     li2 = cote 3 su1 ;
  74.     p3 = li2 poin proc vechoriz ;
  75.  
  76.     li3 = cote 2 su1 ;
  77.  
  78.     si(graph) ;
  79.        titr 'Le maillage de la plaque' ;
  80.        trac su1 ;  
  81.     finsi ;
  82.  
  83. *** Modèles ...
  84.  
  85.     mocnd = mode su1 thermique ;
  86.     macnd = mate mocnd 'K' lambda 'RHO' 1. 'C' lambda ;
  87.  
  88.     moray = mode li2 thermique  rayonnement  INFINI ;
  89.     maray = mate moray 'EMIS' valemis 'E_IN' e0;
  90.  
  91. *** Température imposée ...
  92.  
  93.     blt = bloq 'T' li1 ;
  94.     ti  = depi blt Tempini ;
  95.     lr1 = prog   0 100 ;
  96.     lr2 = prog   1   1 ;
  97.     ev1 = evol manu 't' lr1 'T' lr2 ;
  98.     chti = char 'TIMP' ev1 ti ;
  99.  
  100. *** Température extérieure ...
  101.  
  102.     chtext = manu chpo li2 1 'T' Tempext nature diffus ;
  103.     chte = char 'TERA' ev1 chtext ;
  104.  
  105. *** Préparation de la table pour PASAPAS ...
  106.  
  107.     tabnl = table ;
  108.  
  109.     tabnl . MODELE = mocnd et moray ;
  110.     tabnl . CARACTERISTIQUES = macnd et maray ;
  111.     tabnl . BLOCAGES_THERMIQUES = blt ;
  112.     tabnl . CHARGEMENT = chti et chte ;
  113.  
  114.     tabnl . TEMPS_CALCULES = prog 0. pas  0.1  1. 
  115.                                      pas  0.2  2. ;
  116.  
  117.     tabnl . TEMPERATURES = table ;
  118.     tabnl . TEMPERATURES . 0 = manu chpo su1 1 'T' Tempini ;
  119.  
  120. *    tabnl . RAYONNEMENT = table ;
  121. *    tabnl . RAYONNEMENT . 1 = table ;
  122. *    tabnl . RAYONNEMENT . 1 . TYPE = 'INFINI' ;
  123. *    tabnl . RAYONNEMENT . 1 . MODELE = moray ;
  124.  
  125.     tabnl . 'PROCEDURE_THERMIQUE' = 'DUPONT' ;
  126.     tabnl . 'CTE_STEFAN_BOLTZMANN' = cte_sb ;
  127.  
  128.     tabnl . 'RELAXATION_THETA' = 1. ;
  129.  
  130. *** Appel à PASAPAS ...
  131.  
  132.     pasapas tabnl ;
  133.  
  134.  
  135. *** Petit post-traitement ...
  136.  
  137.     listt = prog ;
  138.     listtp3 = prog ;
  139.     listtp3r = prog ;
  140.     nbpas = dime (tabnl . TEMPS) ;
  141.     repeter surpas nbpas ;
  142.        lindice = &surpas - 1 ;
  143.        listt = listt et (prog (tabnl . TEMPS . lindice)) ;
  144.        valtp3 = extr (tabnl . TEMPERATURES . lindice) T p3 ;
  145.        listtp3 = listtp3 et (prog valtp3) ;
  146.        valtp3r = (valtp3 - Tempext) / (Tempini - Tempext) ;
  147.        listtp3r = listtp3r et (prog valtp3r) ;
  148.  
  149. *  ... Vérification s'il n'y a pas d'oscillations spatiales ...
  150.        si(faux) ;
  151.           chtit = chai 'Profil de temperature au temps ' 
  152.                   (tabnl . TEMPS . lindice) ;
  153.           titr chtit ;
  154.           profil = evol chpo (tabnl . TEMPERATURES . lindice) 
  155.                    'T' li3 ;
  156.           dess profil ;
  157.        finsi ;
  158.  
  159.     fin surpas ;
  160.  
  161.     titr 'Evolution de la temperature au bord de la plaque' ;
  162.     evtp3 = evol manu 't' listt 'T' listtp3 ;
  163.     titr 'Evolution de la temperature relative au bord de la plaque' ;
  164.     evtp3r = evol manu 't' listt 'T' listtp3r ;
  165.  
  166.     si(graph) ;
  167.        dess evtp3  ;
  168.        dess evtp3r ;
  169.     finsi ;
  170.  
  171. *** Vérification si c'est OK ...
  172.  
  173. ************************************************************************
  174. *                                                                      *
  175. *   On va résoudre analytiquement le problème en écrivant l'égalité    *
  176. *   des flux (T = température au bord de la plaque) :                  *
  177. *                                                                      *
  178. *            T_imp - T                                                 *
  179. *   lambda ------------- = epsilon * sigma * (T^4 - T_ext^4)           *
  180. *                L                                                     *
  181. *   avec epsilon(e,e0)                                                 *
  182. *                                                                      *
  183. ************************************************************************
  184.  
  185.     ee = ((1./valemis) + (1./e0) - 1.); ee = 1./ee ;
  186.     mess ' ee: ' ee ;
  187.  
  188.     A_0 = -1 * ((lambda * Tempini / L) + 
  189.                              (ee * cte_sb * (Tempext**4))) ;
  190.     A_1 = lambda / L ;
  191.     A_2 = 0.0 ;
  192.     A_3 = 0.0 ;
  193.     A_4 = ee * cte_sb ;
  194.  
  195.  
  196.     r1  = racp A_0 A_1 A_2 A_3 A_4 ;
  197.  
  198. *    mess 'Zéros du polynôme :' r1 ' ' r2 ' ' r3 ' ' r4 ;
  199.  
  200.     sol_anlt = r1 extr 2 ;
  201.  
  202.     sol_calc = extr listtp3 (dime listtp3) ;
  203.     mess ' anal ' sol_anlt ; mess ' calc ' sol_calc ;
  204.     erreur = 100 * (sol_calc - sol_anlt) / sol_anlt ;
  205.     mess 'Erreur relative = ' erreur '%' ;
  206.  
  207.     si(abs(erreur) > 0.01) ;
  208.        erre 5 ;
  209.     sinon ;
  210.        erre 0 ;
  211.     finsi ;
  212.  
  213. *** Bye ...
  214.  
  215.     fin ;
  216.  
  217.  
  218.  
  219.  
  220.  
  221.  
  222.  
  223.  
  224.  
  225.  
  226.  
  227.  
  228.  
  229.  
  230.  
  231.  
  232.  

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