Télécharger la plaquette Cast3M


fr en

Télécharger rayo-3D-1.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  rayo-3D-1.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4.  
  5. ****************************************************************
  6. *
  7. *  test 3D couplage conduction-rayonnement
  8. *  REFERENCE: cas TPNV 02/89 du guide VPCS
  9. *
  10. *  DONNEES
  11. *  cavité cubique 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 : 1214K (solution analytique)
  23. *  Remarques :
  24. *   2 types d'éléments pour décrire la cavité
  25. *   option convexe pour FFOR ou bien avec plan de symetrie
  26. *   (dans ce cas la cavite est parallelipipedique 1x2x1 d'ou
  27. *    une solution differente) 
  28. *   calcul en permanent et en transitoire (PASAPAS)
  29. ****************************************************************
  30.  
  31. *** Options ...
  32.  
  33.     option echo 0 dime 3 elem cub8 ;
  34.  
  35.     graph = faux ;
  36.  
  37. * options supplementaires: cavite convexe ou plan de symetrie
  38. *       cvxe = vrai ; syme = faux ;
  39. *     cvxe = faux ; syme = faux ;
  40.         cvxe = faux ; syme = vrai ;
  41.  
  42.      2elem = vrai; comm 2 types elements ;
  43.  
  44. *** Paramètres ...
  45.  
  46. *  proprietes physiques
  47.    lam = 20. ;
  48.    e_sup=1.0 ; e_inf=0.1 ; e_late = 0.5 ;
  49.    T_sup = 1000. ; T_inf = 2000. ;
  50.  
  51. *   solution de reference 
  52.      si (syme) ; Tana = 1205.0 ; sinon; Tana = 1214.0 ; finsi;
  53.  
  54. *** Points ...
  55.  
  56.     dens 1. ;
  57.  
  58.     p1 = -0.5 -0.5  0.  ;  p2 =  0.5 -0.5 0. ;
  59.     p3 =  0.5  0.5  0.  ;  p4 = -0.5  0.5 0. ;
  60.     p1p2 = droi 1 p1 p2 ; p2p3 = droi 1 p2 p3 ;
  61.     p3p4 = droi 1 p3 p4 ; p4p1 = droi 1 p4 p1 ;
  62.  
  63. *   pour le plan de symetrie
  64.     Q4 = P4 plus (0. 0. 1.);
  65.  
  66. *** Surfaces ...
  67.  
  68. *   plan inferieur
  69.     sbase  = p1p2 p2p3 p3p4 p4p1 dall plan ;
  70.  
  71. *   plan supérieur
  72.      si 2elem;
  73.     opti elem tri3 ; comm pour avoir 2 types d elements;
  74.      sinon;
  75.     opti elem qua4 ; 
  76.     finsi ;
  77.     cbase = p1p2 et p2p3 et p3p4 et p4p1 ;
  78.     sb_tri3 = surf cbase plane ;
  79.     shaut = inve( sb_tri3 plus (0. 0. 1. ));
  80.  
  81.     opti elem qua4 ;
  82. *   surface laterale
  83. sfron = inve( sbase tour (0.5 -0.5 0.) (0.5  0.5  0.)    90. ) ;
  84. sarri = inve( sbase tour (-0.5 -0.5 0.) (-0.5  0.5  0.) -90. ) ;
  85. sdroi = inve( sbase tour ( 0.5  0.5 0.) (-0.5  0.5  0.)  90. ) ;
  86. sgauc = inve( sbase tour ( 0.5 -0.5 0.) (-0.5 -0.5  0.) -90. ) ;
  87.  
  88. si syme;
  89. *   sarri est le plan de symetrie: points P1 P4 Q4  
  90.     slate =  sfron           et sdroi et sgauc ;
  91. sinon;
  92.      slate =  sfron et sarri et sdroi et sgauc ;
  93. finsi;
  94.  
  95.     elim 0.01 slate ;
  96.  
  97.     cavite =  sbase et shaut et slate ;
  98.  
  99. *   centre du cube ;
  100.     c = 0. 0. 0.5 ;
  101.  
  102. *   epaisseur des parois : dz=10mmm
  103.     dz = 0.01 ; r = 1.01 ;
  104.  
  105.     sbase_e = sbase homo r c ;
  106.     shaut_e = shaut homo r c ;
  107.     slate_e = slate homo r c ;
  108.  
  109. *** Volumes ...
  110.  
  111.    opti elem cub8 ;
  112.     z1 = volu 1  sbase sbase_e ;
  113.     z2 = volu 1  shaut shaut_e ;
  114.     z3 = volu 1  slate slate_e ;
  115.     tout = z1 et z2 et z3 ;
  116.  
  117.     titr 'Le maillage de la cavite' ;
  118.     si(graph) ; trac tout ; finsi ;
  119.  
  120. *** Modélisation ...
  121.  
  122. *   conduction
  123.     lamb = lam/dz ;
  124.  
  125.     mcd  = modeli tout thermique  ;
  126.     k = mate mcd 'K' lamb ;
  127.     cnd = cond mcd k ;
  128.  
  129. si cvxe;
  130.     mr1= modeli sbase thermique  rayonnement  'CAVITE' 'CONVEXE'
  131.         CONS 'CAV1';
  132.     mr2= modeli shaut thermique rayonnement 'CAVITE' 'CONVEXE'
  133.         CONS 'CAV1';
  134.     mr3= modeli slate thermique  rayonnement 'CAVITE' 'CONVEXE'
  135.         CONS 'CAV1'; 
  136. sinon;
  137.   si syme ;
  138.     mr1= modeli sbase thermique  rayonnement 'CAVITE'
  139.           'SYMETRIE' P4 P1 Q4  CONS 'CAV1'  ;
  140.     mr2= modeli shaut thermique rayonnement 'CAVITE'
  141.            'SYMETRIE' P4 P1 Q4  CONS 'CAV1'  ;
  142.     mr3= modeli slate thermique  rayonnement 'CAVITE'
  143.            'SYMETRIE' P4 P1 Q4  CONS 'CAV1'  ;    
  144.   sinon;
  145.      mr1= modeli sbase thermique  rayonnement 'CAVITE' CONS 'CAV1' ;
  146.      mr2= modeli shaut thermique rayonnement 'CAVITE'  CONS 'CAV1' ;
  147.      mr3= modeli slate thermique  rayonnement 'CAVITE' CONS 'CAV1' ;
  148.   finsi;
  149. finsi;
  150.  
  151.     mrt = mr1 et mr2 et mr3 ;
  152.     e1 = mate mr1 'EMIS' e_inf ;
  153.     e2 = mate mr2 'EMIS' e_sup ;
  154.     e3 = mate mr3 'EMIS' e_late ;
  155.     e = e1 et e2 et e3 ;
  156.  
  157. *** calcul des facteurs de forme
  158.  
  159. *    opti impi 1; comm pour bilans et condition de reciprocite ;
  160.     fft = ffor  mrt  e ;
  161. *    opti 'IMPI' 0 ;
  162.  
  163. *** Quelques verifications ...
  164.  
  165. *   list mrt ; list e ;
  166.  
  167. *   list fft ;
  168.     surf = extr fft SURF 1 1 1;
  169. *   mess ' surfaces: ' ; list surf ;
  170.  
  171.     nface = nbel (extr fft 'MAIL');
  172. *    mess ' nface ' nface ;
  173. *
  174. *   impression des facteurs de forme
  175. *
  176. *    si (non 2elem);
  177.  
  178. *    repe iface nface;
  179. *        fi = extr fft MIDL 1 &iface 1;
  180. *         mess ' ';
  181. *         mess ' elem: ' &iface ' Fij: ' ;
  182. *         list fi ;
  183. *    fin iface ;
  184.  
  185. *    finsi;
  186. *
  187. *  
  188.     tta= prray mrt;
  189. *    list tta ;
  190.     tta_int = tta . 1 ;
  191. *    list tta_int;
  192.  
  193. *    opti donn 5;
  194.  
  195. *** Matrice de  rayonnement ...
  196.  
  197.     chamr = raye mrt  fft e ;
  198.  
  199. *** Conditions aux limites ...
  200.  
  201.     c1 = bloq z1   'T' ;
  202.     tim1 = depi c1 T_inf ;
  203.     c2 = bloq z2   'T' ;
  204.     tim2 = depi c2 T_sup ;
  205.  
  206. *** Initialisation de la température ...
  207.  
  208.     tp = manu chpo tout 1 'T' T_inf natu 'DIFFUS' ;
  209.  
  210. *   calcul en permanent
  211. *   -------------------
  212.  
  213. *** Résolution (par itérations) ...
  214.  
  215. *   Coeff. de relaxation ...
  216.     alfa = 0.4 ;
  217.  
  218.     nbpsup = nbno tout ;
  219.  
  220.     listemp = prog ;
  221.     listres = prog ;
  222.  
  223.     maxiter =  100 ;
  224.     critconv = 1.e-5 ;
  225. *    opti echo 1 ;
  226.     REPE bloc1 ;
  227.  
  228.        nbiter = &bloc1 ;
  229.  
  230.        t_cavi   = redu tp cavite ;
  231.        te_cavi   =  chan 'CHAM' t_cavi   mrt 'GRAVITE' ;
  232.        cr = rayn mrt  chamr te_cavi    ;
  233.  
  234.        cndtot = cnd et cr et c1 et c2 ;
  235.  
  236.        residu = cndtot * tp ;
  237.        normres = ((xtx (exco 'Q' residu)) / nbpsup) ** 0.5 ;
  238. *      mess '   La norme du flux résiduel = ' normres ;
  239.        si((nbiter > 1) et (normres < critconv)) ;
  240.           mess 'Convergence à l itération ' (nbiter-1) ;
  241.           quitter bloc1 ;
  242.        finsi ;
  243.        si(nbiter > maxiter) ;
  244.           mess 'Erreur ! Pas de convergence après ' (nbiter-1)
  245.                ' itérations !' ;
  246.           erre 2 ;
  247.           quitter bloc1 ;
  248.        finsi ;
  249.        listres = listres et (prog normres) ;
  250.  
  251. *      mess '-----------------------------------------' ;
  252. *      mess 'Itération N° ' &bloc1 ;
  253.  
  254.        tt = resou cndtot (tim1 et tim2) ;
  255.  
  256.        dt = exco 'T' (tt - tp) 'T' ;
  257.        normdt = ((xtx dt) / nbpsup) ** 0.5 ;
  258. *      mess '   La norme de delta t = ' normdt ;
  259.  
  260.        tn = (alfa * tt) + ((1.-alfa) * tp) ;
  261.        tp =tn ;
  262.  
  263. *      mess '   Température obtenue = ' (maxi(redu tt slate)) ;
  264.  
  265.     FIN bloc1 ;
  266. *   opti echo 1 ;
  267.  
  268. *** Post-traitement ...
  269.  
  270.     t = exco 'T' tt 'T' ;
  271.     si(graph) ;
  272.        titr 'Le champ de temperature final' ;
  273.        trac tout t ;
  274.     finsi ;
  275.     tlate = redu t slate ;
  276.  
  277.     tsol = maxi(tlate) ;
  278.     mess 'La solution obtenue = ' tsol ;
  279.  
  280.     err_s =  ABS((tsol - tana)/tana) ;
  281.     err_s = err_s * 100. ;
  282.     mess 'Erreur relative en permanent' err_s '%' ;
  283.  
  284.  
  285.  
  286. * calcul en transitoire
  287. * ---------------------
  288.  
  289.     Tinit= manu chpo tout 1 'T' T_inf natu 'DIFFUS';
  290.  
  291.     tabnl = table ; 
  292.     tabnl . 'MODELE' =  mcd  et mrt; 
  293.     tabnl . 'CARACTERISTIQUES' = k et e ; 
  294.  
  295.     listtemp = prog 0  0.5 10 ; 
  296.     listval  = prog 1  1   1 ;
  297.  
  298.     evchar   = evol manu 't' listtemp 'f(t)' listval ; 
  299.     tabnl . 'CHARGEMENT' = char 'TIMP' evchar (tim1 et tim2) ;
  300.     tabnl . 'BLOCAGES_THERMIQUES' = c1 et c2 ; 
  301.  
  302.     tabnl . 'TEMPERATURES' = tabl; 
  303.     tabnl . 'TEMPERATURES' . 0 = Tinit;
  304.  
  305.     tabnl . 'TEMPS_CALCULES' = prog 0 pas 0.5 10 ; 
  306.  
  307.     tabnl . 'PROCEDURE_THERMIQUE' = 'DUPONT' ;
  308. *   tabnl . 'PROCEDURE_THERMIQUE' = 'NONLINEAIRE';
  309.  
  310. *** Résolution (avec PASAPAS) ... 
  311.  
  312.      pasapas tabnl ;
  313.  
  314. *** Petit post-traitement ...
  315.  
  316.     nbpas = dime (tabnl . TEMPS) ; 
  317.     listtime = prog ; 
  318.     listtemp = prog ; 
  319.     repeter surpas nbpas ; 
  320.        lindice = &surpas - 1 ; 
  321.        listtime = listtime et (prog (tabnl . TEMPS . lindice)) ; 
  322.        valT =  redu (tabnl . TEMPERATURES . lindice) slate;
  323. *      mess (mini valT) (maxi valT) ; 
  324.  
  325.        valtemp = mini valT ;
  326.        listtemp = listtemp et (prog valtemp) ; 
  327.     fin surpas ; 
  328.  
  329. *       mess (mini valT) (maxi valT) ; 
  330.  
  331.     titr 'Evolution de temperature laterale' ; 
  332.     evtemp = evol manu 't' listtime 'T' listtemp ; 
  333.     si(graph) ; dess evtemp ; finsi ; 
  334.  
  335.     valobt = extr listtemp (dime listtemp) ; 
  336.     err_t =  abs ((valobt - tana) / tana) ; 
  337.     err_t = 100. * err_t;
  338.  
  339.     mess 'Solution obtenue en transitoire avec PASAPAS: ' valobt ; 
  340.     mess 'soit une erreur de ' err_t  ' %' ;
  341.  
  342. *** Test de bon fonctionnement 
  343.  
  344.     si ( (err_s &lt;EG 2.e-1) et (err_t &lt;EG 2.E-1) );
  345.        erre 0 ; 
  346.     sinon; 
  347.        erre 5; 
  348.     finsi ; 
  349.  
  350.    FIN;
  351.  
  352.  
  353.  
  354.  
  355.  
  356.  
  357.  
  358.  
  359.  
  360.  
  361.  
  362.  
  363.  
  364.  
  365.  

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