Accueil | Cast3M

* fichier :  rayoh-2D.dgibi
************************************************************************
************************************************************************
 
************************************************************************
*                                                                      *
* Rayonnement thermique en milieu transparent:                         *
*                                                                      *
* Test de bon fonctionnement de la procédure HRAYO permettant de       *
* traiter le rayonnement face à face ou avec un milieu infini.         *
* Comparaison avec une solution analytique en régime permanent.        *
*                                                                      *
* Calcul 2D-plan.                                                      *
*                                                                      *
* on considère 3 plaques parallèles rayonnantes :                      *
*                                                                      *
*    -------- T3,e3    milieu infini                                   *
*                                                                      *
*      e2sup                                                           *
*    -------- T2                                                       *
*      e2inf           face à face                                     *
*                                                                      *
*    -------- T1,e1                                                    *
*                                                                      *
* La température des plaques supérieure et inférieure est imposée et   *
* uniforme (respectivement 1000K et 2000K) et on détermine la tempéra- *
* ture d'équilibre de la plaque intermédiaire.                         *
* Le rayonnement entre les plaques 1 et 2 est modélisé par un rayonne- *
* ment face à face, celui entre les plaques 2 et 3 par un rayonnement  *
* avec un milieu infini, le milieu infini étant representé par la      *
* plaque 3 non maillée.                                                *
*                                                                      *
* paramètre xquad  vrai : éléments quadratiques                        *
*                  faux : éléments linéaires                           *
************************************************************************
 
graph = faux ;
 
xquad = vrai ;
 
si xquad ;
  option echo 1 dime 2 elem qua8  ;
sinon ;
  option echo 1 dime 2 elem qua4  ;
finsi;
 
*-------------------------------
* Données physiques du problème
*-------------------------------
 
T1 = 1000. ; T3 = 2000. ;
e1 = 0.1 ; e2inf = 0.6 ; e2sup = 0.3 ; e3 = 0.9 ;
lamb = 2.e3 ;
 
*-------------------------------
* Solution analytique
*-------------------------------
 
r1 = 1. - e1 ; r2i = 1. - e2inf ;
r3 = 1. - e3 ; r2s = 1. - e2sup ;
alpha1= (e1*e2inf) / (1.-(r1*r2i)) ;
alpha3= (e3*e2sup) / (1.-(r3*r2s)) ;
M1 = T1**4. ; M3 = T3**4. ;
M2 = ((alpha1*M1) + (alpha3*M3)) / (alpha1+alpha3) ;
T2_ana =  M2**(0.25) ;
 
*-------------------------------
*   Maillage 
*------------------------------- 
 
p  = 0. 0. ; q  = 1. 0. ; s  = droi 2 p q ;
*s1inf = s tour (0.5 0.) 30. ;
s1inf = s ;
plaq1 = s1inf tran 1 (0. 0.001) ;
s1sup = inve ( cote 3 plaq1) ;
 
s2inf = s1inf plus (0. 0.1) ;
plaq2 = s2inf tran 1 (0. 0.001) ;
s2sup = inve ( cote 3 plaq2) ;
 
rac12 = racc 0.1 s1sup s2inf ;
 
tout = plaq1 et plaq2  ;
 
si graph ; trac tout ; finsi ;
 
*------------------------------- 
*  Conduction
*------------------------------- 
 
  mcd = mode tout thermique ;
  k = mate mcd 'K' lamb ;
  cnd = cond mcd k ;
 
 
*------------------------------- 
*       Linéarisation du rayonnement :
*    - définition de modèles de convection
*    - creation des SEG2 définissant les relations entre les supports
*      des champs définis sur s1sup et s2inf (pour le face à face)
*------------------------------- 
 
* face à face
 
  mcv12 = mode rac12 thermique convection ;
 
  si xquad ;
   x_s1sup = chan s1sup SEG2 ;
   x_s2inf = chan s2inf SEG2 ;
   x_rac12 = racc 0.1 x_s1sup x_s2inf ;
 
   rel12 = diff (x_rac12 chan ligne) (x_s1sup et x_s2inf) ;
  sinon;
   rel12 = diff (rac12 chan ligne) (s1sup et s2inf) ;
  finsi;
 
* milieu infini
 
  mcv2s = mode s2sup thermique convection ;
 
 
 
*------------------------------- 
*  Modï¿½les de rayonnement et ï¿½missivitï¿½s
*------------------------------- 
 
* face ï¿½ face
  mr1 = mode s1sup thermique rayonnement FAC_A_FAC 
      s1sup s2inf mcv12 rel12 cons 'FAC1';
  emi1 = mate mr1 'EMIS' e1 ;
 
  mr2i = mode s2inf thermique rayonnement FAC_A_FAC
   s1sup s2inf mcv12 rel12 cons 'FAC1';
 
  emi2i= mate mr2i 'EMIS' e2inf ;
 
* milieu infini
 
  mr2s = mode s2sup thermique rayonnement INFINI CONS 'INF1';
  emi2s= mate mr2s 'EMIS' e2sup 'E_IN' e3 ;
 
  t3inf = manu chpo s2sup 1 'T' T3 ;
 
 
*------------------------------- 
*  Conditions aux limites
*------------------------------- 
 
    c1 = bloq plaq1   'T' ;
    tim1 = depi c1 T1 ;
 
*------------------------------- 
*  Algorithme de résolution
*------------------------------- 
 
*** Initialisation de la température ...
 
    tp = manu chpo tout 1 'T' T3 natu 'DIFFUS';
    nbpsup = nbno tout ;
 
*** Résolution (par itérations) ...
 
*   Coeff. de relaxation ...
    alfa = 0.3 ;
 
    listemp = prog ;
    listres = prog ;
 
    maxiter = 100 ;
    critconv = 1.e-5 ;
*   opti echo 1 ;
    REPE bloc1 ;
 
       nbiter = &bloc1 ;
 
*   traitement du rayonnement:
*                              face à face (1,2)
 
       t1sup = redu (exco 'T' tp 'T') s1sup;
       t2inf = redu (exco 'T' tp 'T') s2inf;
       h12 = HRAYO mcv12 mr1 emi1 t1sup mr2i emi2i t2inf rel12 ;
 
       cr12 = cond mcv12  h12      ;
 
*                              milieu infini (2,3)
 
       t2sup = redu (exco 'T' tp 'T') s2sup;
       h23 = HRAYO mcv2s mr2s emi2s t2sup t3inf ;
 
 
*      mr23 = MATE mcv2s 'H' h23 ;
       cr23 = cond mcv2s h23     ;
       f3 = conv mcv2s  h23 t3inf ;
 
       crtot= cr12 et cr23  ;
 
       cndtot = crtot et cnd et c1 ;
       residu = (cndtot * tp) - f3 ;
       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 f3) ;
       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 ;
 
 
    FIN bloc1 ;
*   opti echo 1 ;
 
*------------------------------- 
*  Post-traitement ...
*------------------------------- 
 
    titre 'Champ de température final' ;
    si(graph) ; trac plaq2   (exco 'T' tn) ; finsi ;
 
    Temp = exco 'T' tt 'T' ;
    T2 = redu plaq2   Temp ;
 
    RESI=maxi( ABS((T2 - T2_ana)/T2_ana));
    mess 'Solution analytique ' T2_ana ;
    mess 'Erreur relative ' (100*RESI) '%' ;
    SI (RESI &lt;EG 1E-4);
       ERRE  0;
    SINO;
       ERRE  5;
    FINSI;
 
fin;