Télécharger la plaquette Cast3M


fr en

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

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  rayo-2D-4-bis.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
  76. *                                                       suivant[j+1] += actuel[j
  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 * l
  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,
  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
  130. *               fclose(f_out);
  131. *       }
  132. *
  133. *************************************************************************
  134.  
  135.  
  136. *** Options ...
  137.  
  138.     opti dime 2 elem seg2 echo 1 ;
  139.  
  140. *** Paramètres ...
  141.  
  142.     L = 1. ;
  143.     ep = L / 10 ;
  144.  
  145.     graph = faux ;
  146.  
  147.     Tempini = 273. ;
  148.     Tempext = 0. * Tempini ;
  149.     cte_sb = 5.673e-8 ;
  150.     Gamma = 1. ;
  151.     si ( neg Gamma 0. ) ;
  152.        lambda = (cte_sb * Tempini * Tempini * Tempini * L) / Gamma ;
  153.        mess 'Lambda = ' lambda ;
  154.        valemis = 1. ;
  155.     sinon ;
  156.        lambda = 1. ;
  157.        mess 'Lambda choisie arbitrairement = ' lambda ;
  158.        valemis = 0. ;
  159.     finsi ;
  160.     Bi = 0.2 ;
  161.  
  162. *** Points ...
  163.  
  164.     dens 0.1 ;
  165.     p1 = 0 0 ;
  166.     p2 = 0 ep ;
  167.  
  168.     vechoriz = L 0 ;
  169.  
  170. *** Lignes ...
  171.  
  172.     li1 = p1 d 1 p2 ;
  173.  
  174. *** Surface ...
  175.  
  176.     opti elem qua4 ;
  177.     su1 = li1 tran vechoriz dini 0.1 dfin 0.01 ;
  178.  
  179.     li2 = cote 3 su1 ;
  180.     p3 = li2 poin proc vechoriz ;
  181.  
  182.     li3 = cote 2 su1 ;
  183.  
  184.     si(graph) ;
  185.        titr 'Le maillage de la plaque' ;
  186.        trac su1 ;  
  187.     finsi ;
  188.  
  189. *** Modèles ...
  190.  
  191.     mocnd = mode su1 thermique ;
  192.     macnd = mate mocnd 'K' lambda 'RHO' 1. 'C' lambda ;
  193.  
  194.     mocnv = mode li2  thermique convection CONS 'CCONV';
  195.     macnv = mate mocnv 'H' (Bi*lambda/L) ;
  196.  
  197.     moray = mode li2 thermique rayonnement INFINI
  198.     CONS 'CRAYO' ;
  199.     maray = mate moray 'EMIS' valemis 'E_IN' 1.;
  200.  
  201. *** Préparation de la table pour PASAPAS ...
  202.  
  203.     tabnl = table ;
  204.  
  205.     tabnl . MODELE = mocnd et mocnv et moray;
  206.     tabnl . CARACTERISTIQUES = macnd et macnv et maray ;
  207.  
  208.     lr1 = prog   0 100 ;
  209.     lr2 = prog   1   1 ;
  210.     ev1 = evol manu 't' lr1 'g(t)' lr2 ;
  211.     chtext = manu chpo li2 1 'T' Tempext ;
  212.     tabnl . CHARGEMENT =( char 'TECO' ev1 chtext) et 
  213.                         ( char 'TERA' ev1 chtext)  ;
  214.  
  215.     tabnl . TEMPS_CALCULES = prog 0. pas 2.e-4 0.001
  216.                                      pas 5.e-4 0.01 
  217.                                      pas 0.005 0.1
  218.                                      pas  0.05  1. 
  219.                                      pas   0.2 10. ;
  220.  
  221.     tabnl . TEMPERATURES = table ;
  222.     tabnl . TEMPERATURES . 0 = manu chpo su1 1 'T' Tempini ;
  223.  
  224.  
  225.     tabnl . 'PROCEDURE_THERMIQUE' = 'DUPONT' ;
  226.     tabnl . 'CTE_STEFAN_BOLTZMANN' = cte_sb ;
  227.  
  228.     tabnl . 'RELAXATION_THETA' = 1. ;
  229.  
  230. *** Appel à PASAPAS ...
  231.  
  232.     pasapas tabnl ;
  233.  
  234. *** Petit post-traitement ...
  235.  
  236.     listt = prog ;
  237.     listtp1 = prog ;
  238.     listtp1r = prog ;
  239.     listtp3 = prog ;
  240.     listtp3r = prog ;
  241.     nbpas = dime (tabnl . TEMPS) ;
  242.     repeter surpas nbpas ;
  243.        lindice = &surpas - 1 ;
  244.        listt = listt et (prog (tabnl . TEMPS . lindice)) ;
  245.        valtp1 = extr (tabnl . TEMPERATURES . lindice) T p1 ;
  246.        listtp1 = listtp1 et (prog valtp1) ;
  247.        valtp1r = (valtp1 - Tempext) / (Tempini - Tempext) ;
  248.        listtp1r = listtp1r et (prog valtp1r) ;
  249.        valtp3 = extr (tabnl . TEMPERATURES . lindice) T p3 ;
  250.        listtp3 = listtp3 et (prog valtp3) ;
  251.        valtp3r = (valtp3 - Tempext) / (Tempini - Tempext) ;
  252.        listtp3r = listtp3r et (prog valtp3r) ;
  253.  
  254. *  ... Vérification s'il n'y a pas d'oscillations spatiales ...
  255.        si(faux) ;
  256.           chtit = chai 'Profil de temperature au temps ' 
  257.                   (tabnl . TEMPS . lindice) ;
  258.           titr chtit ;
  259.           profil = evol chpo (tabnl . TEMPERATURES . lindice) 
  260.                    'T' li3 ;
  261.           dess profil ;
  262.        finsi ;
  263.  
  264.     fin surpas ;
  265.  
  266.     titr 'Evolution de la temperature au centre de la plaque' ;
  267.     evtp1 = evol manu 't' listt 'T' listtp1 ;
  268.     titr 'Evolution de la temperature relative au centre de la plaque' ;
  269.     evtp1r = evol manu 't' listt 'T' listtp1r ;
  270.  
  271.     titr 'Evolution de la temperature au bord de la plaque' ;
  272.     evtp3 = evol manu 't' listt 'T' listtp3 ;
  273.     titr 'Evolution de la temperature relative au bord de la plaque' ;
  274.     evtp3r = evol manu 't' listt 'T' listtp3r ;
  275.  
  276.     si(graph) ;
  277.        dess evtp1  ;
  278.        dess evtp1r ;
  279.        dess evtp3  ;
  280.        dess evtp3r ;
  281.     finsi ;
  282.  
  283. *** Vérification si c'est OK ...
  284.  
  285.     listtref = prog 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. ;
  286.     listfref = prog 0.566607 0.437195 0.348971 0.283422 0.232317 
  287.                     0.19139 0.158117 0.130835 0.108357 0.0897874 ;
  288.     titr 'Evolution de la temperature relative au bord de la plaque'
  289.          ' (M-C)' ;
  290.     evref = evol manu 't' listtref 'T' listfref ;
  291.     si(graph) ;
  292.        dess (evref et evtp3r) ;
  293.     finsi ;
  294.  
  295.     nbtests = dime listtref ;
  296.     maxerr = 0. ;
  297.     repeter surtst nbtests ;
  298.        instant = extr listtref &surtst ;
  299.        valref  = extr listfref &surtst ;
  300.        tcalc = peche tabnl TEMPERATURES instant ;
  301.        valcalc = extr tcalc T p3 ;
  302.        valcalcr = (valcalc - Tempext) / (Tempini - Tempext) ;
  303.        errrel = valcalcr - valref ;
  304.        mess 'Erreur relative à t =' instant 'est de' errrel ;
  305.        errrela = abs errrel ;
  306.        si(errrela > maxerr) ;
  307.           maxerr = errrela ;
  308.        finsi ;
  309.     fin surtst ;
  310.     si(maxerr > 0.01) ;
  311.        erre 5 ;
  312.     finsi ;
  313.  
  314. *** Bye ...
  315.  
  316.    fin ;
  317.  
  318.  
  319.  
  320.  
  321.  
  322.  

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