Accueil | Cast3M

* fichier :  rayoh-2D.dgibi
************************************************************************
************************************************************************
 
************************************************************************
*                                                                      *
* Rayonnement thermique en milieu transparent:                         *
*                                                                      *
* Test du rayonnement face à face                                      *
* Comparaison avec une solution analytique en régime permanent.        *
*                                                                      *
* Calcul par 2 méthodes:                                               *
*      1. algorithme iteratif                                          *
*      2. utilisation de PASAPAS en transitoire jusqu'au permanent     *
*                                                                      *
*                                                                      *
* Calcul 2D-plan.                                                      *
*                                                                      *
* on considère 3 plaques parallèles rayonnantes :                      *
*                                                                      *
*    -------- T3,e3                                                    *
*                                                                      *
*      e2sup           face à face                                     *
*    -------- 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, et de meme pour celui entre les plaques 2 et 3     *
*                                                                      *
* 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) ;
*LIST T2_ANA ;
 
*-------------------------------
*   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 ;
 
s3inf = s2inf plus (0. 0.1) ;
plaq3 = s3inf tran 1 (0. 0.001) ;
s3sup = inve ( cote 3 plaq3) ;
 
rac23 = racc 0.1 s2Sup s3inf ;
 
tout = plaq1 et plaq2 et plaq3 ;
 
si graph ; trac tout ; finsi ;
 
*------------------------------- 
*  Conduction
*------------------------------- 
 
  mcd = mode tout thermique ;
  k   = mate mcd 'K' lamb 'RHO' 1. 'C' LAMB;
  cnd = cond mcd k ;
 
  cte_sb = 5.673e-8 ;
 
*------------------------------- 
*  Modeles de rayonnement et emissivites
*------------------------------- 
 
* face a face
  mr1 = mode rac12 thermique rayonnement FACE_A_FACE ;
  emi1 = mate mr1 'EMI1' e1 'EMI2' e2inf ;
 
  mr2S= mode rac23 thermique rayonnement FACE_A_FACE ;
  emi2S= mate mr2S 'EMI1' e2SUP 'EMI2' e3 ;
 
*------------------------------- 
*  Conditions aux limites
*------------------------------- 
 
    c1 = bloq plaq1   'T' ;
    tim1 = depi c1 T1 ;
 
    c3 = bloq plaq3   'T' ;
    tim3 = depi c3 T3 ;
 
    lr1 = prog   0 100 ;
    lr2 = prog   1   1 ;
 
    ev1 = evol manu 't' lr1 'T' lr2 ;
    cht1 = char 'TIMP' ev1 tIM1 ;
    ev3 = evol manu 't' lr1 'T' lr2 ;
    cht3 = char 'TIMP' ev3 tIM3 ;
 
*** Températures initiales ...
 
    T2 = 500. ;
    chINI1 = manu chpo PLAQ1 1 'T' T1 nature diffus ;
    chINI2 = manu chpo PLAQ2 1 'T' T2 nature diffus ;
    chINI3 = manu chpo PLAQ3 1 'T' T3 nature diffus ;
 
 
*----------------------------------------------------------
* 1.  Première méthode de résolution : Algorithme itératif 
*----------------------------------------------------------
 
*** Initialisation de la température ...
 
    tp = CHINI1 ET CHINI2 ET CHINI3 ;
    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)
 
       h12  = HRAYO 'FACE_A_FACE' mr1 emi1 tp ;
       cr12 = cond mr1  h12 ;
 
*                              face à face (2,3)
 
       h23  = HRAYO 'FACE_A_FACE' mr2S emi2S tp ;
       cr23 = cond mr2s h23 ;
 
       crtot  = cr12 et cr23  ;
       cndtot = crtot et cnd et c1 ET C3;
       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 TIM3) ;
       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 ...
*------------------------------- 
 
    Temp = exco 'T' tt 'T' ;
    T2A = redu plaq2   Temp ;
    RESIA=maxi( ABS((T2A - T2_ana)/T2_ana));
 
*------------------------------------------------------
* 2.  Deuxieme méthode de résolution : appel a PASAPAS
*------------------------------------------------------
 
 
*** Préparation de la table pour PASAPAS ...
 
    tabnl = table ;
 
    tabnl . MODELE = MCD ET MR1 ET MR2S ;
 
 
    tabnl . CARACTERISTIQUES = K ET EMI1 ET EMI2S ;
    tabnl . BLOCAGES_THERMIQUES = C1 ET C3 ;
    tabnl . CHARGEMENT = cht1 et cht3 ;
 
    tabnl . TEMPS_CALCULES = prog 0. pas  0.1  1. 
                                     pas  0.2  4. ;
 
    tabnl . TEMPERATURES = table ;
    tabnl . TEMPERATURES . 0 = CHINI1 ET CHINI2 ET CHINI3 ;
 
 
    tabnl . 'PROCEDURE_THERMIQUE' = 'DUPONT' ;
    tabnl . 'CTE_STEFAN_BOLTZMANN' = cte_sb ;
 
    tabnl . 'RELAXATION_THETA' = 1. ;
 
*** Appel à PASAPAS ...
 
    pasapas tabnl ;
 
 
*** Petit post-traitement ...
 
    listt = prog ;
    listtp3 = prog ;
    listtp3r = prog ;
    nbpas = dime (tabnl . TEMPS) ;
    repeter surpas nbpas ;
       lindice = &surpas - 1 ;
       listt = listt et (prog (tabnl . TEMPS . lindice)) ;
       valtp3 = MAXI (REDU tabnl . TEMPERATURES . lindice PLAQ2) ;
       listtp3 = listtp3 et (prog valtp3) ;
    fin surpas ;
 
    titr 'Evolution de la temperature de la plaque2';
    evtp3 = evol manu 't' listt 'T' listtp3 ;
 
    si(graph) ;
       dess evtp3  ;
    finsi ;
 
    T2B = redu tabnl . TEMPERATURES . (nbpas - 1) PLAQ2  ;
    RESIB=maxi( ABS((T2B - T2_ana)/T2_ana));
 
*  comparaison a la solution analytique 
 
    mess 'Solution analytique ' T2_ana ;
 
    mess 'Solution Algorithme ' (MAXI T2A) ;
    mess 'Erreur relative    ' (100*RESIA) '%' ;
 
    mess 'Solution PASAPAS    ' (MAXI T2B) ;
    mess 'Erreur relative    ' (100*RESIB) '%' ;
 
    SI ((RESIA &lt;EG 1.E-4)et(RESIB &lt;EG 1.E-4)) ;
       ERRE  0;
    SINO;
       ERRE  5;
    FINSI;
 
 
    fin ;