Accueil | Cast3M

* fichier :  rayo_abs-axi-2.dgibi
************************************************************************
************************************************************************
 
complet=faux;
****************************************************************
* Calcul de la température d'une cavité contenant
* un milieu absorbant (Tg,k_abs)
*
*  DONNEES
*  cavité cylindrique de cote 1.
*
*            1000K e=1.0
*            -----
*        axe |    | e=0.5
*            |    |
*            -----
*            2000K e=0.1
*
*
*  RESULTATS: flux rayonné total sur la cavité
*             solution numérique de référence
*
*  On teste 2 méthodes:
*  1- la méthode basée sur le calcul de la matrice de rayonnement
*     (opérateur RAYE) R telle que Phi = R. T**4 - R. Tg**4
*  2- la méthode basée sur le calcul de la température de
*     rayonnement Trad telle que  Phi = emis*sigma*(T**4-Trad**4)
*
****************************************************************
 
*** Options ...
 
*   option echo 0 dime 2 elem qua8 mode axis;
    option echo 0 dime 2 elem qua4 mode axis;
 
    graph = faux ;
 
*** Solutions numeriques de référence à la creation du test
*   flux total rayonné sur la cavité évalué avec les méthodes 1 et 2
 
    si (complet); n = 100 ;
     fref1 = -6.83843E6;
     fref2 = -6.83870E6;
    sinon;   n= 20 ;
     fref1 = -6.95186E6;
     fref2 = -6.95813E6;
    finsi;
 
*** Paramètres ...
 
    l = 1.0 ;
 
*   epaisseur des parois : dz = 10 mm
*<
*   dz = 0.01  ;
    dz = 0.05 ;
*>
    dzp = dz ;
    dzn = -1. * dz ;
 
*   proprietes physiques
    lam = 10.;
*>
    e_sup=1.0 ; e_inf=0.1 ; e_late = 0.5 ;
*   e_sup=1.0 ; e_inf=1.0 ; e_late = 1.0 ;
    T_sup = 1000. ; T_inf = 2000. ;
*<
    sig = 5.67e-8;
    Nr = sig* (T_sup**3.) * L /lam ; mess 'Nr= ' Nr;
 
*   milieu absorbant : coeff et température (K)
    k_abs = -1.e-0 ;
    Tg = 2500. ;
 
*** Points ...
 
    dens l ;
*   n = 4   ;
 
    p4 = 0. l  ; p3 = l  l  ;
    p1 = 0. 0. ; p2 = l  0. ;
    q4 = 0. l  ; q3 = l  l  ;
    q1 = 0. 0. ; q2 = l  0. ;
    l1 = d n p1 p2 ; l2 = d n p3 p4 ;
    l3 = d n q2 q3 ;
    late = l3  ;
 
*   la surface laterale est disjointe des deux autres
    cavite = l1 et l2 et late ;
*<
    elim 1e-5 cavite ;
*>
 
ne = 5;
    z1 = l1 trans ne (0. dzn) ;
    z2 = l2 trans ne (0. dzp) ;
    z3 = l3 trans ne (dzp 0.) ;
 
    tout= z1 et z2 et z3 ;
    titr 'Le maillage de la cavite' ;
    si(graph) ; trac tout ; finsi ;
 
*** Modélisation ...
 
*   conduction
    lamb = lam/dz ;
 
    mcd  = modeli tout thermique  ;
    k = mate mcd 'K' lamb ;
    cnd = cond mcd k ;
 
*
    mr1= modeli l1  thermique    rayonnement 'CAVITE' 'CONVEXE'  
      'CONS'  'CAV1';
    mr2= modeli l2    thermique    rayonnement  'CAVITE' 'CONVEXE'
     'CONS'  'CAV1';
    mrl= modeli late  thermique    rayonnement'CAVITE' 'CONVEXE'  
        'CONS'  'CAV1' ;
    mrt = mr1 et mr2 et mrl ;
 
*** Emissivités
 
    e1 = mate mr1 'EMIS' e_inf  'CABS' k_abs 'TABS' Tg   ;
    e2 = mate mr2 'EMIS' e_sup    'CABS' k_abs 'TABS' Tg    ;
    el = mate mrl 'EMIS' e_late   'CABS' k_abs 'TABS' Tg    ;
    e = e1 et e2 et el ;
 
*** Calcul des facteurs de forme
 
* facteurs de forme "classiques" pour memoire
*   ffc   = ffor  mrt  e  ;
* facteurs de forme généralisés
    fft   = ffor  mrt  e ;
 
*** Conditions aux limites ...
 
    c1 = bloq z1   'T' ;
    tim1 = depi c1 T_inf ;
    c2 = bloq z2   'T' ;
    tim2 = depi c2 T_sup ;
 
*** méthode 1
*   ---------
 
*   calcul de la matrice de rayonnement
 
    chamr = raye mrt  fft e ;
 
 
*** Initialisation de la température ...
 
    tp = manu chpo tout 1 'T' T_inf natu 'DIFFUS';
 
*   calcul du flux associé au milieu absorbant
*<
    Tg_abs = manu chpo tout   1 'T' Tg natu 'DIFFUS';
    T_abs   = redu (exco 'T' Tg_abs 'T') cavite ;
    Te_abs   =  chan 'CHAM' T_abs mrt 'GRAVITE' ;
 
       cr_abs = rayn mrt chamr Te_abs ;
 
       Fg_abs = cr_abs * T_abs ;
*      list Fg_abs;
*>
 
 
*** Résolution (par itérations) ...
 
*   Coeff. de relaxation ...
    alfa = 0.42 ;
 
    nbpsup = nbno tout ;
 
    listemp = prog ;
    listres = prog ;
 
    maxiter = 100 ;
    critconv = 1.e-5 ;
 
    REPE bloc1 maxiter ;
 
       nbiter = &bloc1 ;
 
       t_cavi   = redu (exco 'T' tp 'T') cavite ;
       te_cavi   =  chan 'CHAM' t_cavi mrt 'GRAVITE' ;
       cr = rayn mrt chamr te_cavi ;
 
       cndtot = cnd et cr et c1 et c2 ;
       residu = (cndtot * tp) - Fg_abs ;
 
       normres = ((xtx (exco 'Q' residu)) / nbpsup) ** 0.5 ;
*      mess '   La norme du flux résiduel = ' normres ;
       si((nbiter > 1) et (normres < critconv)) ;
          mess 'Convergence à l itération ' (nbiter-1) ;
          quitter bloc1 ;
       finsi ;
       si(nbiter > maxiter) ;
          mess 'Erreur ! Pas de convergence après ' (nbiter-1)
               ' itérations !' ;
          erre 2 ;
          quitter bloc1 ;
       finsi ;
       listres = listres et (prog normres) ;
 
*      mess '---------------------------------------' ;
*      mess 'Itération N° ' &bloc1 ;
 
       tt1 = resou cndtot (tim1 et tim2 et Fg_abs) ;
 
       dt = exco 'T' (tt1- tp) 'T' ;
        normdt = ((xtx dt) / nbpsup) ** 0.5 ;
*       mess '   La norme de delta t = ' normdt ;
 
       tn = (alfa * tt1) + ((1.-alfa) * tp) ;
       tp =tn ;
 
    FIN bloc1 ;
*
 
 
*** Post-traitement ...
 
    mess '1- calcul avec matrice de rayonnement';
 
*   tsol = extr tt1 'T' q2 ;
*   mess 'La solution obtenue = ' tsol ;
 
*      flux rayonne
 
       fray  =  (cr * tt1) - Fg_abs ;
       fray_1 = maxi (resu fray);
       mess ' bilan global cavite: ' fray_1 ;
 
       mess ' flux  surface inf: ' (maxi (resu (redu fray l1))) ;
       mess ' flux  surface sup: ' (maxi (resu (redu fray l2))) ;
       mess ' flux  surface late: ' (maxi (resu (redu fray late))) ;
 
*      flux associé à la condition de température imposée
 
       mess ' Timp surf. inf.: ' (maxi (resu (reac c1 tt1))) ;
       mess ' Timp surf. sup.: ' (maxi (resu (reac c2 tt1))) ;
 
    t = exco 'T' tt1 'T' ;
    titr 'Le champ de temperature final' ;
    si(graph) ; trac tout t ; finsi ;
    tlate = redu t late ; comm list tlate ;
 
*!  trac (redu tt1 z3 ) z3 ;
 
    ev1 = evol 'CHPO' (redu t  l3) l3;
    ev1 = ev1 coul 'VERT' ;
 
    RESI1 =ABS((fray_1 - fref1)/fref1);
    mess 'Erreur relative - methode 1' (100*RESI1) '%' ;
 
*   opti donn 5;
 
*   méthode 2
*   ---------
 
*** Initialisation de la température ...
 
    tp = manu chpo tout 1 'T' T_inf natu 'DIFFUS';
 
*** Résolution (par itérations) ...
 
*   Coeff. de relaxation ...
 
    alfa = 1.0  ;
 
    nbpsup = nbno tout ;
 
    listemp = prog ;
    listres = prog ;
 
    maxiter =  100 ;
    critconv = 1.e-5 ;
 
    REPE bloc2 ;
 
       nbiter = &bloc2 ;
 
       t_cavi   = redu (exco 'T' tp 'T') cavite ;
       te_cavi   =  chan 'CHAM' t_cavi mrt 'GRAVITE' ;
 
* calcul de la temperature de rayonnement associée: trad
       trad = raye mrt fft e te_cavi 1e-7    ;
 
* calcul du coefficient d'echange
       hrad = HRCAV  mrt e te_cavi trad ;
 
       trad_n1= chan 'CHPO' mrt trad       ;
       trad_n = nomc trad_n1 'T' 'NATU' 'DIFFUS';
 
* pour la condition de convection
       cr = cond mrt hrad;
       f = conv mrt hrad trad_n;
 
       cndtot = cnd et c1 et c2 et cr ;
 
       residu = (cndtot * tp) -f ;
       normres = ((xtx (exco 'Q' residu)) / nbpsup) ** 0.5 ;
*      mess '   La norme du flux résiduel = ' normres ;
       si((nbiter > 1) et (normres < critconv)) ;
          mess 'Convergence à l itération ' (nbiter-1) ;
          quitter bloc2 ;
       finsi ;
       si(nbiter > maxiter) ;
          mess 'Erreur ! Pas de convergence après ' (nbiter-1)
               ' itérations !' ;
          erre 2 ;
          quitter bloc2 ;
       finsi ;
       listres = listres et (prog normres) ;
 
*      mess '---------------------------------------' ;
*      mess 'Itération N° ' &bloc2 ;
 
       tt2 = resou cndtot (tim1 et tim2 et f) ;
 
       dt = exco 'T' (tt2- tp) 'T' ;
        normdt = ((xtx dt) / nbpsup) ** 0.5 ;
*       mess '   La norme de delta t = ' normdt ;
 
       tn = (alfa * tt2) + ((1.-alfa) * tp) ;
       tp =tn ;
 
*      mess '   Température obtenue = ' (extr tt2 'T' q2) ;
 
    FIN bloc2 ;
 
 
*temps impr place sgac ;
 
 
*** Post-traitement ...
 
    mess '2- calcul avec température de rayonnement';
 
    t = exco 'T' tt2 'T' ;
    titr 'Le champ de temperature final' ;
    si(graph) ; trac tout t ; finsi ;
 
    tlate = redu t late ;
 
*   mess ' T laterale: min et max ' (mini tlate) (maxi tlate);
 
    ev2 = evol 'CHPO' (redu t  l3) l3;
    ev2 = ev2 coul 'BLEU' ;
 
*      flux rayonne
 
       fray  =  (cr * tt2)- f;
       fray_2   = maxi (resu fray);
       mess ' bilan global cavite: '  fray_2 ;
       mess ' flux  surface inf: ' (maxi (resu (redu fray l1))) ;
       mess ' flux  surface sup: ' (maxi (resu (redu fray l2))) ;
       mess ' flux  surface late: ' (maxi (resu (redu fray late))) ;
 
*      flux associé à la condition de température imposée
 
       mess ' Timp surf. inf.: ' (maxi (resu (reac c1 tt2))) ;
       mess ' Timp surf. sup.: ' (maxi (resu (reac c2 tt2))) ;
 
 
    RESI2 =ABS((fray_2 - fref2)/fref2);
    mess 'Erreur relative - methode 2' (100*RESI2) '%' ;
 
    SI ((RESI1 &lt;EG 1E-4) et (RESI2 &lt;EG 1E-4));
      ERRE  0;
    SINO;
      ERRE  5;
    FINSI;
 
si graph;
    TAB1 = table;
    TAB1. 'TITRE' = table;
    TAB1. 'TITRE' . 1 = 'MOT' 'methode 1' ;
    TAB1. 'TITRE' . 2 = 'MOT' 'methode 2' ;
 
    dess (ev1 et ev2) mima  lege TAB1 ;
finsi;
 
    fin;