Accueil | Cast3M

* fichier :  rayo-2D-3.dgibi
************************************************************************
************************************************************************
 
************************************************************************
*                                                                      *
* Rayonnement thermique en milieu transparent:                         *
*                                                                      *
* vérification du bon fonctionnement de l'opérateur FFOR dans le cas   *
* général (traitement des parties cachées) sur un cas thermiquement    *
* simple avec axe de symétrie.                                         *
*                                                                      *
* Test 2D-plan sur un problème à symétrie cylindrique 1D:              *
*                                --------------------                  *
* on considère 3 cylindres coaxiaux de rayons 25mm, 50mm et 100mm      *
* conducteurs et d'emissivité 1. La température des cylindres intérieur*
* et extérieur est uniforme (respectivement 1000K et 2000K) et on      *
* détermine la température d'équilibre du cylindre intermédiaire.      *
* (conventions: indices 1,2,3 pour les cylindres de rayon 25,50,100mm) *
*                                                                      *
* Jeu de données identique au cas 2D-axis :rayo-axi-2.dgibi            *
* hormis la valeur des 2 paramètres:                                   *
*            'xrev'  vrai : géométrie axisymétrique                    *
*                    faux :           2D plane                         *
*            'x2cav' vrai : on considère les 2 cavités séparément      *
*                    faux :              l'ensemble des 2 cavités      *
*                                                                      *
************************************************************************
 
xrev = faux ; x2cav = vrai ;
 
graph = faux ;
 
si xrev ;
 option echo 1 dime 2 elem qua4 mode axis ;
sinon;
 option echo 1 dime 2 elem qua4  ;
finsi;
 
*-------------------------------
* Données physiques du problème
*-------------------------------
 
r1 = 25e-3 ; r2 = 50e-3 ; r3 = 100e-3 ;
T1 = 1000. ; T3 = 2000. ;
emis = 1. ; lamb = 20. ;
 
*-------------------------------                                        
* Solution analytique
*-------------------------------
 
si xrev ;
 f21 = (r1*r1)/(r2*r2) ;
sinon;
 f21 = r1/r2 ;
finsi;
f22 = 1. - f21 ;
M1 = T1**4. ; M3 = T3**4. ;
M2 = ( (f21*M1) + M3 )/(2.-f22) ;
T2_ana =  M2**(0.25) ;
 
*-------------------------------
*   Maillage 
*------------------------------- 
 
C  =  0. 0. ;
 
* sphère 1
 
p1_sud =  0. ( -1. * r1 )  ;
p1_equa =  r1  0. ;
p1_nord =  0. r1 ;
d1 = r1/3 ;
c1_inf = cerc p1_sud C p1_equa DINI d1 DFIN d1;
c1_sup = cerc p1_equa C p1_nord DINI d1 DFIN d1;
c1_tot = c1_inf et c1_sup ;
c1_int = c1_tot homo C 0.999 ;
c1_ext = c1_tot homo C 1.001 ;
sphere1 = regl c1_int c1_ext ;
 
* sphère 2
 
p2_sud =  0. ( -1. * r2 )  ;
p2_equa =  r2  0. ;
p2_nord =  0. r2 ;
d2 = r2/3 ;
c2_inf = cerc p2_sud C p2_equa DINI d2 DFIN d2;
c2_sup = cerc p2_equa C p2_nord DINI d2 DFIN d2;
c2_tot = c2_inf et c2_sup ;
c2_int = c2_tot homo C 0.999 ;
c2_ext = c2_tot homo C 1.001 ;
sphere2 = regl c2_int c2_ext ;
 
* sphère 3
 
p3_sud =  0. ( -1. * r3 )  ;
p3_equa =  r3  0. ;
p3_nord =  0. r3 ;
d3 = r3/3 ;
c3_inf = cerc p3_sud C p3_equa DINI d3 DFIN d3;
c3_sup = cerc p3_equa C p3_nord DINI d3 DFIN d3;
c3_tot = c3_inf et c3_sup ;
c3_int = c3_tot homo C 0.999 ;
c3_ext = c3_tot homo C 1.001 ;
sphere3 = regl c3_int c3_ext ;
 
* on controle l'orientation des elements pour les cavités rayonnantes
 
cavi_12 = (inve c1_ext) et c2_int ;
cavi_23 = (inve c2_ext) et c3_int ;
cavi_tot = cavi_12 et cavi_23 ;
 
tout = sphere1 et sphere2 et sphere3 ;
 
si graph ; 
trac tout  ;
finsi ; 
 
*------------------------------- 
*  Conduction
*------------------------------- 
 
mcd = mode tout thermique ;
k = mate mcd 'K' lamb ;
cnd = cond mcd k ;
 
*------------------------------- 
*  Modèles de rayonnement
*------------------------------- 
pptt= 0    1.;
'SI' xrev;
  tt = 'TEXTE' '  ';
 'SINON' ;
  tt=  'TEXTE' 'SYMETRIE  C  PPTT';
 'FINSI';
 
si x2cav;
  mr12 = mode cavi_12  thermique rayonnement 'CAVITE' 
  SYMETRIE C pptt CONS 'CAV1';
  e12 = mate mr12 'EMIS' emis ;
 
  mr23 = mode cavi_23 thermique rayonnement 'CAVITE' 
   tt CONS 'CAV1';
  e23 = mate mr23 'EMIS' emis ;
sinon;
  mrtot= mode cavi_tot thermique rayonnement   'CAVITE' 
   tt;
  etot = mate mrtot 'EMIS' emis ;
finsi;
 
*------------------------------- 
*  Facteurs de forme et matrices de rayonnement
*------------------------------- 
 
si x2cav ;
*opti 'IMPI' 1 ;
 si xrev ;
  ff12 = ffor mr12 e12;
  ff23 = ffor mr23 e23;
 sinon;
  ff12 = ffor mr12  e12;
  ff23 = ffor mr23  e23;
 finsi;
*opti 'IMPI' 0 ;
 
 cham12 = raye mr12 ff12 e12 ;
 cham23 = raye mr23 ff23 e23 ;
 
sinon;
*opti 'IMPI' 1 ;
 si xrev ;
  fftot= ffor mrtot etot;
 sinon;
  fftot = ffor mrtot etot;
 finsi;
*opti 'IMPI' 0 ;
 
 chamtot = raye mrtot fftot etot ;
 
finsi;
 
*------------------------------- 
*  Conditions aux limites
*------------------------------- 
 
    c1 = bloq sphere1 'T' ;
    tim1 = depi c1 T1 ;
    c3 = bloq sphere3 'T' ;
    tim3 = depi c3 T3 ;
 
*------------------------------- 
*  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
 
    si x2cav ;
       t_c12 = redu (exco 'T' tp 'T') cavi_12 ;
       te_c12 = chan 'CHAM' t_c12   mr12 'GRAVITE' ;
       cr12 = rayn mr12  cham12 te_c12 ;
 
       t_c23 = redu (exco 'T' tp 'T') cavi_23 ;
       te_c23 = chan 'CHAM' t_c23   mr23 'GRAVITE' ;
       cr23 = rayn mr23  cham23 te_c23 ;
       crtot = cr12 et cr23 ;
 
    sinon;
       t_tot = redu (exco 'T' tp 'T') cavi_tot;
       te_tot = chan 'CHAM' t_tot mrtot  'GRAVITE' ;
       crtot= rayn mrtot chamtot te_tot ;
 
    finsi;
 
       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 ...
*------------------------------- 
 
    titre 'Champ de température final' ;
    si(graph) ; trac sphere2 (exco 'T' tn) ; finsi ;
 
    Temp = exco 'T' tt 'T' ;
    T2 = redu sphere2 Temp ;
 
    RESI=maxi( ABS((T2 - T2_ana)/T2_ana));
    mess 'Problème à symétrie cylindrique';
    mess 'Solution analytique ' T2_ana ;
    mess 'Erreur relative ' (100*RESI) '%' ;
    SI (RESI &lt;EG 2E-3);
       ERRE  0;
    SINO;
       ERRE  5;
    FINSI;
 
fin;