Accueil | Cast3M

* fichier :  rayoh-3D.dgibi
************************************************************************
************************************************************************
 
************************************************************************
*                                                                      *
* Rayonnement thermique en milieu transparent:                         *
*                                                                      *
* Test de bon fonctionnement de la procédure HRAYO permettant de       *
* traiter le rayonnement face à face par comparaison à une solution    *
* analytique en régime permanent.                                      *
*                                                                      *
* Calcul 3D.                                                           *
*                                                                      *
* 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                            *
* Remarque: dans le cas quadratique, et pour le face à face chaque     *
*           frontière doit etre un maillage simple (TRI6 uniquement    *
*           ou QUA8 uniquement)                                        *
************************************************************************
 
graph = faux ;
 
xquad = faux ;
 
si (xquad) ;
  option echo 0 dime 3 elem pr15 ;
sinon ;
  option echo 0 dime 3 elem cub8 ;
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 
*------------------------------- 
 
oeil=10 10 10;
p1=0 0 0 ; p2=1 0 0 ; p3=1 1 0 ; p4=0 1 0;
l1=p1 d 1 p2;
l2=p2 d 1 p3;
l3=p3 d 1 p4;
l4=p4 d 2 p1;
si xquad ;
* on recherche un objet simple
  opti elem tri6 ;
  s1= (l1 et l2 et l3 et l4) surf plan ;
  option echo 1 dime 3 elem pr15 ;
sinon ;
  s1= (l1 et l2 et l3 et l4) surf plan ;
finsi ;
 
n=1 ;
v = 0 0 0.1 ;  v1 = (-1)*v ;
w = 0 0 0.001 ;  w1 = (-1)*w ;
 
 
 s2inf = s1 ;
 s2sup = s2inf plus w ;
 plaq2 = s2inf volu n s2sup ;
 
 s1sup = s2inf plus v1;
 s1inf = s1sup plus w1;
 plaq1 = s1inf volu n s1sup ;
 
 rac12 = liai 0.1 s1sup s2inf ;
 
 tout = plaq1 et plaq2  ;
 
si graph ; trac oeil tout ; finsi ;
 
*------------------------------- 
*  Conduction
*------------------------------- 
 
  mcd = mode tout thermique ;
  k = mate mcd 'K' lamb ;
  cnd = cond mcd k ;
 
*------------------------------- 
*       Linéarisation du rayonnement :
*    - 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_s1 = (s1sup chan TRI6) chan TRI3 ;
  x_s2 = (s2inf chan TRI6) chan TRI3 ;
  xx_s1 = x_s1 chan ligne;
  xx_s2 = x_s2 chan ligne;
  x_s12 = liai 0.1 x_s1 x_s2;
  xx_s12 = x_s12 chan ligne ;
 
  rel12 = diff xx_s12 (xx_s1 et xx_s2) ;
 
sinon ;
  sr12 = rac12 chan ligne ;
  sr1 = s1sup chan ligne ;
  sr2 = s2inf chan ligne ;
  rel12 = diff sr12 (sr1 et sr2) ;
 
finsi ;
 
  si graph ; trac oeil rel12 ; 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-2 ;
*   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 ;
 
*      mr12 = MATE mcv12 'H' h12 ;
*      mr12 = h12 ;
*      cr12 = cond mcv12 mr12    ;
       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 ;
*   list T2 ;
 
    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;