Télécharger la plaquette Cast3M


fr en

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

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  rayo-2D-4.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4.  
  5. *************************************************************
  6. *                                                           *
  7. *   Calcul d'une plaque infinie (de largeur 2L) soumise à   *
  8. *   de la convection et du rayonnement.                     *
  9. *                                                           *
  10. *   Modélisation plane.                                     *
  11. *                                                           *
  12. *   Auteurs : Michel Bulik & Nadia Coulon                   *
  13. *                                                           *
  14. *   Date : Octobre 1995                                     *
  15. *                                                           *
  16. *   Références :                                            *
  17. *   ------------                                            *
  18. *   [1] Klaus-Jürgen Bathe & Mohammad R. Khoshgoftaar,      *
  19. *       Finite element formulation and solution of          *
  20. *       non-linear heat transfer, Nuclear Engineering and   *
  21. *       Design, v. 51 (1979), pp. 389-401                   *
  22. *                                                           *
  23. *   [2] J. Joly, Cas tests non linéaires de validation pour *
  24. *       DELFINE, Note technique EMT.SMTS.TTMF 84/29         *
  25. *                                                           *
  26. *   [3] Haji-Sheikh, A. & Sparrow, A. M., The solution of   * 
  27. *       Heat Conduction Problems by Probability Methods,    *
  28. *       Transactions of ASME, Journal of Heat Transfer, 89  *
  29. *       (1967), pp. 121-130                                 *
  30. *                                                           *
  31. *************************************************************
  32. *
  33. * La solution de référence a été obtenue par une méthode de
  34. * Monte-Carlo décrite dans [3] et que réalise le programme
  35. * suivant en C. Après sa compilation il convient de demander
  36. * l'éditions de liens avec les subroutines quarti, cubic et
  37. * quadra de Castem2000.
  38. *
  39. *       #include<stdio.h>
  40. *       
  41. *       #define NB_TOTAL        100000
  42. *       #define NB_INTERVAL     20
  43. *       #define NB_PAS          8000
  44. *       
  45. *       double  actuel[NB_INTERVAL+1], suivant[NB_INTERVAL+1], phi[NB_PAS+1] ;
  46. *       int     les_nk[NB_PAS+1] ;
  47. *       
  48. *       #define Bi      0.2
  49. *       #define Gamma   1.0
  50. *       
  51. *       double  petit_l, delta_t ;
  52. *       double  A_0, A_1, A_2, A_3, A_4 ;
  53. *       
  54. *       char    nomfichier[128] ;
  55. *       FILE    *f_out ;
  56. *       
  57. *       void main()
  58. *       {
  59. *               int     i, j, k, approx ;
  60. *       
  61. *               /*** Initialisation ***/
  62. *               for(i=0;i<=NB_INTERVAL;i++) actuel[i] = 0 ;
  63. *               actuel[0] = NB_TOTAL ;
  64. *               les_nk[0] = actuel[0] ;
  65. *       
  66. *               /*** Boucle sur les pas - calcul des n_k ***/
  67. *               for(i=1;i<=NB_PAS;i++) {
  68. *       
  69. *                       for(j=0;j<=NB_INTERVAL;j++) suivant[j] = 0 ;
  70. *       
  71. *                       for(j=0;j<=NB_INTERVAL;j++) {
  72. *                               if(j==0)                suivant[j+1] += actuel[j] ;
  73. *                               else if(j==NB_INTERVAL) suivant[j-1] += actuel[j] ;
  74. *                               else {
  75. *                                                       suivant[j-1] += actuel[j] * 0.5 ;
  76. *                                                       suivant[j+1] += actuel[j] * 0.5 ;
  77. *                               }
  78. *                       }
  79. *       
  80. *                       for(j=0;j<=NB_INTERVAL;j++) actuel[j] = suivant[j] ;
  81. *       
  82. *                       les_nk[i] = actuel[0] ;
  83. *               }
  84. *       
  85. *               petit_l = 1.0 / NB_INTERVAL ;
  86. *       
  87. *               A_1 = 2.0 + (Bi * petit_l) ;
  88. *               A_2 = 0.0 ;
  89. *               A_3 = 0.0 ;
  90. *               A_4 = Gamma * petit_l ;
  91. *       
  92. *               phi[0] = 1.0 ;
  93. *       
  94. *               /*** Boucle sur les pas - calcul des phi ***/
  95. *               for(i=1;i<=NB_PAS;i++) {
  96. *       
  97. *                       double  phiact ;
  98. *                       double  q1, q2, q3, q4 ;
  99. *                       int     nroot ;
  100. *       
  101. *                       A_0 = -1.0 - phi[i-1] ;
  102. *                       if(i > 1) {
  103. *                               double  lephi, tototo ;
  104. *       
  105. *                               for(j=1;j<=i-1;j++) {
  106. *                                       lephi = phi[i - j] ;
  107. *                                       tototo = les_nk[j] * lephi * petit_l *
  108. *                                               (Bi + (Gamma * lephi * lephi * lephi)) ;
  109. *                                       tototo /= NB_TOTAL ;
  110. *                                       A_0 += tototo ;
  111. *                               }
  112. *                       }
  113. *       
  114. *                       if( Gamma > 0.0) {
  115. *                               quarti(&A_4,&A_3,&A_2,&A_1,&A_0,&q1,&q2,&q3,&q4,&nroot) ;
  116. *                               phiact = q1 ;
  117. *                       }
  118. *                       else {
  119. *                               phiact = -A_0 / A_1 ;
  120. *                       }
  121. *       
  122. *                       phi[i] = phiact ;
  123. *               }
  124. *       
  125. *               sprintf(nomfichier,"tw.Bi=%3.1f.Gamma=%3.1f.xy\0",Bi,Gamma) ;
  126. *               f_out = fopen(nomfichier,"w") ;
  127. *               if(f_out == NULL) exit(5) ;
  128. *               delta_t = petit_l * petit_l / 2 ;
  129. *               for(i=1;i<=NB_PAS;i++) fprintf(f_out,"%15g %15g\n",i*delta_t,phi[i]) ;
  130. *               fclose(f_out);
  131. *       }
  132. *
  133. *************************************************************************
  134.  
  135. *** Options ...
  136.  
  137.     opti dime 2 elem seg2 echo 1 ;
  138.  
  139. *** Paramètres ...
  140.  
  141.     L = 1. ;
  142.     ep = L / 10 ;
  143.  
  144.     graph = faux ;
  145.  
  146.     Tempini = 273. ;
  147.     Tempext = 0. * Tempini ;
  148.     cte_sb = 5.673e-8 ;
  149.     Gamma = 1. ;
  150.     si ( neg Gamma 0. ) ;
  151.        lambda = (cte_sb * Tempini * Tempini * Tempini * L) / Gamma ;
  152.        mess 'Lambda = ' lambda ;
  153.        valemis = 1. ;
  154.     sinon ;
  155.        lambda = 1. ;
  156.        mess 'Lambda choisie arbitrairement = ' lambda ;
  157.        valemis = 0. ;
  158.     finsi ;
  159.     Bi = 0.2 ;
  160.  
  161. *** Points ...
  162.  
  163.     dens 0.1 ;
  164.     p1 = 0 0 ;
  165.     p2 = 0 ep ;
  166.  
  167.     vechoriz = L 0 ;
  168.  
  169. *** Lignes ...
  170.  
  171.     li1 = p1 d 1 p2 ;
  172.  
  173. *** Surface ...
  174.  
  175.     opti elem qua4 ;
  176.     su1 = li1 tran vechoriz dini 0.1 dfin 0.01 ;
  177.  
  178.     li2 = cote 3 su1 ;
  179.     p3 = li2 poin proc vechoriz ;
  180.  
  181.     li3 = cote 2 su1 ;
  182.  
  183.     si(graph) ;
  184.        titr 'Le maillage de la plaque' ;
  185.        trac su1 ;  
  186.     finsi ;
  187.  
  188. *** Modèles ...
  189.  
  190.     mocnd = mode su1 thermique ;
  191.     macnd = mate mocnd 'K' lambda 'RHO' 1. 'C' lambda ;
  192.  
  193.     mocnv = mode li2  thermique convection CONS 'CCONV';
  194.     macnv = mate mocnv 'H' (Bi*lambda/L) ;
  195.  
  196.     moray = mode li2  thermique rayonnement 
  197.        infini CONS 'CRAYO' ;
  198.     maray = mate moray 'EMIS' valemis ;
  199.  
  200. *** Préparation de la table pour PASAPAS ...
  201.  
  202.     tabnl = table ;
  203.  
  204.     tabnl . MODELE = mocnd et mocnv  et moray;
  205.     tabnl . CARACTERISTIQUES = macnd et macnv et maray ;
  206.  
  207.     lr1 = prog   0 100 ;
  208.     lr2 = prog   1   1 ;
  209.     ev1 = evol manu 't' lr1 'g(t)' lr2 ;
  210.     chtext = manu chpo li2 1 'T' Tempext ;
  211.     tabnl . CHARGEMENT =( char 'TECO' ev1 chtext) et 
  212.                         ( char 'TERA' ev1 chtext)  ;
  213.  
  214.     tabnl . TEMPS_CALCULES = prog 0. pas 2.e-4 0.001
  215.                                      pas 5.e-4 0.01 
  216.                                      pas 0.005 0.1
  217.                                      pas  0.05  1. 
  218.                                      pas   0.2 10. ;
  219.  
  220.     tabnl . TEMPERATURES = table ;
  221.     tabnl . TEMPERATURES . 0 = manu chpo su1 1 'T' Tempini ;
  222.  
  223. *    tabnl . RAYONNEMENT = table ;
  224. *    tabnl . RAYONNEMENT . 1 = table ;
  225. *    tabnl . RAYONNEMENT . 1 . TYPE = 'INFINI' ;
  226. *    tabnl . RAYONNEMENT . 1 . MODELE = moray ;
  227.  
  228.     tabnl . 'PROCEDURE_THERMIQUE' = 'DUPONT' ;
  229.     tabnl . 'CTE_STEFAN_BOLTZMANN' = cte_sb ;
  230.  
  231.     tabnl . 'RELAXATION_THETA' = 1. ;
  232.  
  233. *** Appel à PASAPAS ...
  234.  
  235.     pasapas tabnl ;
  236.  
  237. *** Petit post-traitement ...
  238.  
  239.     listt = prog ;
  240.     listtp1 = prog ;
  241.     listtp1r = prog ;
  242.     listtp3 = prog ;
  243.     listtp3r = prog ;
  244.     nbpas = dime (tabnl . TEMPS) ;
  245.     repeter surpas nbpas ;
  246.        lindice = &surpas - 1 ;
  247.        listt = listt et (prog (tabnl . TEMPS . lindice)) ;
  248.        valtp1 = extr (tabnl . TEMPERATURES . lindice) T p1 ;
  249.        listtp1 = listtp1 et (prog valtp1) ;
  250.        valtp1r = (valtp1 - Tempext) / (Tempini - Tempext) ;
  251.        listtp1r = listtp1r et (prog valtp1r) ;
  252.        valtp3 = extr (tabnl . TEMPERATURES . lindice) T p3 ;
  253.        listtp3 = listtp3 et (prog valtp3) ;
  254.        valtp3r = (valtp3 - Tempext) / (Tempini - Tempext) ;
  255.        listtp3r = listtp3r et (prog valtp3r) ;
  256.  
  257. *  ... Vérification s'il n'y a pas d'oscillations spatiales ...
  258.        si(faux) ;
  259.           chtit = chai 'Profil de temperature au temps ' 
  260.                   (tabnl . TEMPS . lindice) ;
  261.           titr chtit ;
  262.           profil = evol chpo (tabnl . TEMPERATURES . lindice) 
  263.                    'T' li3 ;
  264.           dess profil ;
  265.        finsi ;
  266.  
  267.     fin surpas ;
  268.  
  269.     titr 'Evolution de la temperature au centre de la plaque' ;
  270.     evtp1 = evol manu 't' listt 'T' listtp1 ;
  271.     titr 'Evolution de la temperature relative au centre de la plaque' ;
  272.     evtp1r = evol manu 't' listt 'T' listtp1r ;
  273.  
  274.     titr 'Evolution de la temperature au bord de la plaque' ;
  275.     evtp3 = evol manu 't' listt 'T' listtp3 ;
  276.     titr 'Evolution de la temperature relative au bord de la plaque' ;
  277.     evtp3r = evol manu 't' listt 'T' listtp3r ;
  278.  
  279.     si(graph) ;
  280.        dess evtp1  ;
  281.        dess evtp1r ;
  282.        dess evtp3  ;
  283.        dess evtp3r ;
  284.     finsi ;
  285.  
  286. *** Vérification si c'est OK ...
  287.  
  288.     listtref = prog 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. ;
  289.     listfref = prog 0.566607 0.437195 0.348971 0.283422 0.232317 
  290.                     0.19139 0.158117 0.130835 0.108357 0.0897874 ;
  291.     titr 'Evolution de la temperature relative au bord de la plaque'
  292.          ' (M-C)' ;
  293.     evref = evol manu 't' listtref 'T' listfref ;
  294.     si(graph) ;
  295.        dess (evref et evtp3r) ;
  296.     finsi ;
  297.  
  298.     nbtests = dime listtref ;
  299.     maxerr = 0. ;
  300.     repeter surtst nbtests ;
  301.        instant = extr listtref &surtst ;
  302.        valref  = extr listfref &surtst ;
  303.        tcalc = peche tabnl TEMPERATURES instant ;
  304.        valcalc = extr tcalc T p3 ;
  305.        valcalcr = (valcalc - Tempext) / (Tempini - Tempext) ;
  306.        errrel = valcalcr - valref ;
  307.        mess 'Erreur relative à t =' instant 'est de' errrel ;
  308.        errrela = abs errrel ;
  309.        si(errrela > maxerr) ;
  310.           maxerr = errrela ;
  311.        finsi ;
  312.     fin surtst ;
  313.     si(maxerr > 0.01) ;
  314.        erre 5 ;
  315.     finsi ;
  316.  
  317. *** Bye ...
  318.  
  319.    fin ;
  320.  
  321.  
  322.  
  323.  
  324.  
  325.  

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