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