* fichier :  rayo-axi-1.dgibi
************************************************************************
* Section : Thermique Rayonnement
* Section : Thermique Convection
* Section : Thermique Conduction
************************************************************************

************************************************************************
*
*  CONDUCTION-RAYONNEMENT
*  2D- axisymetrique      permanent
*  facteurs de forme : option convexe
*
*  Reference : Engelman Int.J.for Num.Fluids 1991 Vol.13
*
*  Le domaine est un cylindre droit R=0.0053m H=0.063m
*  On determine la temperature sur tout le domaine en prenant en 
*  compte la conduction a l interieur du domaine et le rayonnement 
*  entre les faces delimitant le domaine.
*   (conductivité : 0.026 W/m/K)
*  Les surfaces superieure et inferieure sont a temperature imposee.
*
*              306 K
*             Ec=0.065
*          -----<------
*         '            '
*         '            ' Ev=0.45
*        axe           '
*         '            '
*          ----->------
*             Ef=0.88
*             298 K
*
* RESULTAT: temperature a mi-hauteur du cylindre sur l axe
*                           TRIO-EF
*           maillage 10*50  301.6 K
*           maillage  5*20  301.6 K
************************************************************************

*** Options ...

    opti echo 0 dime 2 mode axis elem qua4 ;

    graph = faux ;

*** Paramètres ...

    R = (0.0106/2.) ; L = 0.0633 ;

    NR = 5  ; NZ = 20 ; NZ2 = NZ/2 ;
*   NR = 10 ; NZ = 50 ; NZ2 = NZ/2 ;

*   proprietes physiques
    lamb = 0.026 ;
    e_sup=0.065 ;e_inf=0.88  ; e_late = 0.45 ;
    T_sup =  306. ; T_inf = 298.  ; T_ref = 300. ;

*** Points ...

    dens (R / NR) ;

    P4 = 0. L  ; P3 = R  L  ;
    P42= 0. (0.5*L) ;
    P1 = 0. 0. ; P2 = R  0. ;

*** Lignes ...

    L41 = D NZ2 P4 P42 ; L42 = D NZ2 P42 P1 ; L4 = L41 et L42 ;

    L1 = D NR P1 P2 ; L3 = D NZ P2 P3 ;
    L2 = D NR P3 P4 ;
    L5 = INVE (L4) ;

    cavite = L1 ET L3 ET L2 ;
    tout = L1 L3 L2 L4 DALLE PLAN ;
    nbpsup = nbno tout ;

    titr 'Le maillage du cylindre' ;
    si(graph) ; trac tout ; finsi ;

*** Modélisation ...

    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'  ;
    mr3 = modeli l3   thermique rayonnement 'CAVITE' 'CONVEXE' 
       CONS 'CAV1'   ;
    mrt  = mr1 et mr2 et mr3  ;

    e1 = mate mr1 'EMIS' e_inf ;
    e3 = mate mr3 'EMIS' e_late ;
    e2 = mate mr2 'EMIS' e_sup ;
    e = e1 et e2 et e3 ;

*** Matrice de rayonnement ...

*    opti 'IMPI' 3 ;
    fft = ffor mrt e;
*   fft = ffor mrt        ;
*    opti 'IMPI' 3 ;

    chamr = raye mrt  fft e ;
*    opti 'IMPI' 0 ;

*** Conditions aux limites ...

    c1 = bloq l1   'T' ;
    tim1 = depi c1 T_inf ;
    c2 = bloq l2   'T' ;
    tim2 = depi c2 T_sup ;

*** Initialisation de la température ...

    tp = manu chpo tout 1 'T' T_sup natu 'DIFFUS' ;

*    tmoye  = (T_sup+T_inf) / 2. ;
*    tecart = (T_sup-T_inf) / 2. ;
*    tp = (manu chpo tout 1 'T' tmoye natu 'DIFFUS') +
*         (manu chpo l2 1 'T' tecart natu 'DIFFUS') +
*         (manu chpo l1 1 'T' (-1*tecart) natu 'DIFFUS') ;

*** Résolution ...

*   Coeff. de relaxation ...
    alfa = 0.4 ;

    listres = prog ;

    maxiter = 100 ;
    critconv = 1.e-8 ;
    opti echo 0 ;
    REPE bloc1 ;

       nbiter = &bloc1 ;

       t_cavi   = redu tp cavite ;
       te_cavi   =  chan 'CHAM' t_cavi   mrt 'GRAVITE' ;
       cr = rayn mrt  chamr te_cavi    ;

       cndtot = cr et cnd et c1 et c2 ;

       residu = cndtot * tp ;
       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 ;

       tt = resou cndtot (tim1 et tim2) ;

       dt = exco 'T' (tt - tp) 'T' ;
       normdt = ((xtx dt) / nbpsup) ** 0.5 ;
       mess '   La norme de delta t = ' normdt ;

       tn = (alfa * tt) + ((1.-alfa) * tp) ;
       tp =tn ;

       mess '   Température obtenue = ' (extr tt 'T' p42) ;

    FIN bloc1 ;
*   opti echo 1 ;

*** Post-traitement ...

    t = exco 'T' tt 'T';
    ta = t/T_ref ;

    t5= redu t l5   ;
*    list t5 ;
    ta5= redu ta l5  ;

    titr 'T/T0 sur l axe vertical';
    ev5 = evol 'CHPO' ta5   l5 ;
    si(graph) ; dess ev5 ; finsi ;

    opti isov lign ;
    titr 'Le champ de temperature final ';
    si(graph) ; trace t tout (cont tout) ; finsi ;

    tana = 301.6 ;
    tsol = extr tt 'T' p42;
    mess 'La solution = ' tsol ;

    RESI=ABS((tsol -tana)/tana);
    mess 'Erreur relative ' (100*RESI) '%' ;
    SI(RESI <EG 1E-4);
      ERRE  0;
    SINO;
      ERRE  5;
    FINSI;

*** Bye ...

    FIN;