Accueil | Cast3M

* fichier :  rayo-3D-2.dgibi
************************************************************************
************************************************************************
 
complet = faux;
*
*************************************************
*  pour calcul complet mettre complet à : vrai;
***************************************************
************************************************************************
*                                                                      *
* 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) avec plan de symétrie sur   *
* un cas thermiquement simple.                                         *
*                                                                      *
* Test 3D sur un problème à symétrie sphérique 1D:                     *
*                           ------------------                         *
* on considère 3 sphères concentriques de rayons 25mm, 50mm et 100mm   *
* conductrices et d'emissivité 1. La température des sphères intérieure*
* et extérieure est uniforme (respectivement 1000K et 2000K) et on     *
* détermine la température d'équilibre de la sphère intermédiaire.     *
* (conventions: indices 1,2,3 pour les sphères de rayons 25,50,100mm)  *
*                                                                      *
*                                                                      *
* Remarques:                                                           *
*  avec la densité choisie, la cavité comprend 96 éléments             *
*  ce problème est également traité en 2D-axis: cf. rayo-axi-2.dgibi   *
************************************************************************
 
x2cav = faux ;
 
graph = faux ;
 
 option echo 1 dime 3 elem tri3 ;
 
*-------------------------------
* 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
*-------------------------------
 
f21 = (r1*r1)/(r2*r2) ;
f22 = 1. - f21 ;
M1 = T1**4. ; M3 = T3**4. ;
M2 = ( (f21*M1) + M3 )/(2.-f22) ;
T2_ana =  M2**(0.25) ;
 
*-------------------------------
*   Maillage 
*------------------------------- 
 
C  =  0. 0. 0. ;
O1 = 0. 1.E-6 0. ;
O2 = 0. -1.E-6 0. ;
oeil = -10. -10. 10. ;
 
* sphère 1
 
p1_sud =  0. 0. (-1.*r1)  ;
p1_equa =  r1  0. 0. ;
p1_nord =  0. 0. r1 ;
 
* densité
aa = 2.;
*si complet ; aa = 2. ; sinon; aa = 1.;finsi;
d1 = r1/aa  ;
 
* cercle
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_droi = c1_inf et c1_sup ;
c1_gauc = inve(SYME c1_droi 'DROIT' p1_sud p1_nord) ;
c1_tot = c1_droi et c1_gauc;
elim 1.e-6 c1_tot ;
* sphere de reference
s1_ref = surf c1_tot 'SPHERIQUE' O2 ;
 
s1_int = s1_ref homo C 0.99 ;
s1_ext = s1_ref homo C 1.01 ;
 
opti elem pri6 ;
sphere1 = volu s1_int s1_ext ;
 
* spheres 2 et 3
 
s2_int = s1_ref homo C (0.999*(r2/r1)) ;
s2_ext = s1_ref homo C (1.001*(r2/r1)) ;
sphere2 = volu s2_int s2_ext ;
 
s3_int = s1_ref homo C (0.999*(r3/r1)) ;
s3_ext = s1_ref homo C (1.001*(r3/r1)) ;
sphere3 = volu s3_int s3_ext ;
 
* on controle l'orientation des elements pour les cavités rayonnantes
 
cavi_12 = (inve s1_ext) et s2_int ;
cavi_23 = (inve s2_ext) et s3_int ;
cavi_tot = cavi_12 et cavi_23 ;
 
tout = sphere1 et sphere2 et sphere3 ;
 
si graph ; 
trac tout oeil cache ;
finsi ; 
 
 
*------------------------------- 
*  Conduction
*------------------------------- 
 
mcd = mode tout thermique ;
k   = mate mcd 'K' lamb ;
cnd = cond mcd k ;
 
*------------------------------- 
*  Modèles de rayonnement
*------------------------------- 
si x2cav;
  mr12 = mode cavi_12 thermique rayonnement 'CAVITE' 'SYMETRIE'
   C (1. 0. 0.) (0. 0. 1.)  'CONS' 'CAV1' ;
  e12 = mate mr12 'EMIS' emis ;
 
  mr23 = mode cavi_23 thermique rayonnement 'CAVITE' 'SYMETRIE'
   C (1. 0. 0.) (0. 0. 1.)  'CONS' 'CAV1'  ;
  e23 = mate mr23 'EMIS'   emis ;
sinon;
  mrtot= mode cavi_tot thermique rayonnement 'CAVITE' 'SYMETRIE'
  C (1. 0. 0.) (0. 0. 1.)  ;
  etot = mate mrtot 'EMIS' emis ;
finsi;
 
 
*------------------------------- 
*  Facteurs de forme et matrices de rayonnement
*------------------------------- 
 
* le plan y = 0. est le plan de symétrie
si x2cav ;
*opti 'IMPI' 1 ;
  ff12 = ffor mr12 e12 ;
  ff23 = ffor mr23 e23 ;
*opti 'IMPI' 0 ;
 
 cham12 = raye mr12 ff12 e12 ;
 cham23 = raye mr23 ff23 e23 ;
 
sinon;
*opti 'IMPI' 1 ;
  fftot = ffor mrtot etot;
*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 ;
si complet;
    critconv = 1.e-5 ;
sinon; 
    critconv = 10. ;
finsi;
*   opti echo 1 ;
    REPE bloc1 ;
 
       nbiter = &bloc1 ;
 
*   traitement du rayonnement
 
 
       t_tot = redu (exco 'T' tp 'T') cavi_tot;
       te_tot = chan 'CHAM' t_tot mrtot  'GRAVITE' ;
       crtot= rayn mrtot chamtot te_tot ;
 
 
 
       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 ;
   mess 'resultat trouve' (maxi t2);
    RESI=maxi( ABS((T2 - T2_ana)/T2_ana));
    mess 'Problème à symétrie sphérique';
    mess 'Solution analytique ' T2_ana ;
    mess 'Erreur relative ' (100*RESI) '%' ;
si complet; valres = 1.2e-2; sinon; valres = 3.e-2; finsi;
    SI (RESI &lt;EG valres);
       ERRE  0;
    SINO;
       ERRE  5;
    FINSI;
 
fin;