* 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 <EG 2E-3);
ERRE 0;
SINO;
ERRE 5;
FINSI;
fin;
© Cast3M 2003 - All rights reserved.
Disclaimer