complet=faux;
****************************************************************
* fichier : rayo_abs-2D-1.dgibi
****************************************************************
*
*
* Calcul de la température d'une cavité carrée contenant
* un milieu absorbant (Tg,k_abs)
*
* DONNEES
* cavité cylindrique de cote 1.
*
* 1000K e=1.0
* -----
* e=0.5 | | e=0.5
* | |
* -----
* 2000K e=0.1
*
*
* RESULTATS: flux rayonné total sur la cavité
* solution numérique de référence
*
* On teste 2 méthodes:
* 1- la méthode basée sur le calcul de la matrice de rayonnement
* (opérateur RAYE) R telle que Phi = R. T**4 - R. Tg**4
* 2- la méthode basée sur le calcul de la température de
* rayonnement Trad telle que Phi = emis*sigma*(T**4-Trad**4)
*
****************************************************************
*** Options ...
* option echo 0 dime 2 elem qua8 ;
option echo 1 dime 2 elem qua4 ;
graph = faux ;
*** Solutions numeriques de référence à la creation du test
* flux total rayonné sur la cavité évalué avec les méthodes 1 et 2
si (complet); n = 100 ;
fref1 = -2.13324E6;
fref2 = -2.13291E6;
sinon; n= 20 ;
fref1 = -2.16680E6;
fref2 = -2.16711E6;
finsi;
*** Paramètres ...
* l = 1.e-2 ;
l = 1. ;
* epaisseur des parois : dz = 10 mm
*<
* dz = 0.01 ;
dz = 0.05*l ;
*>
dzp = dz ;
dzn = -1. * dz ;
* proprietes physiques
lam = 10.;
*>
e_sup=1.0 ; e_inf=0.1 ; e_late = 0.5 ;
* e_sup=1.0 ; e_inf=1.0 ; e_late = 1.0 ;
T_sup = 1000. ; T_inf = 2000. ;
*<
sig = 5.67e-8;
Nr = (sig* (T_sup**3.) * L) / lam ;
* mess 'Nr ' Nr;
* milieu absorbant : coeff et température (K)
* epaisseur optique
tau = 1.0 ;
k_abs = -1. * tau / l ;
* k_abs = -1. * k_abs;
* mess ' k_abs: ' k_abs;
Tg = 2500. ;
*** Points ...
dens l ;
* n = 4 ;
p4 = 0. l ; p3 = l l ;
p1 = 0. 0. ; p2 = l 0. ;
q4 = 0. l ; q3 = l l ;
q1 = 0. 0. ; q2 = l 0. ;
l1 = d n p1 p2 ; l2 = d n p3 p4 ;
l3 = d n q2 q3 ; l4 = d n q4 q1 ;
late = l3 et l4 ;
* la surface laterale est disjointe des deux autres
cavite = l1 et l2 et late ;
*<
elim (1e-5*l) cavite ;
*>
ne = 5;
z1 = l1 trans ne (0. dzn) ;
z2 = l2 trans ne (0. dzp) ;
z3 = l3 trans ne (dzp 0.) ;
z4 = l4 trans ne (dzn 0.) ;
tout= z1 et z2 et z3 et z4;
titr 'Le maillage de la cavite' ;
si(graph) ; trac tout ; finsi ;
*** Modélisation ...
* conduction
lamb = lam/dz ;
mcd = modl tout thermique ;
k = mate mcd 'K' lamb ;
cnd = cond mcd k ;
*
mr1= modl l1 thermique rayonnement 'CAVITE' 'CONS' 'TROU1';
mr2= modl l2 thermique rayonnement 'CAVITE' 'CONS' 'TROU1';
mrl= modl late thermique rayonnement 'CAVITE' 'CONS' 'TROU1';
mrt = mr1 et mr2 et mrl ;
*** Emissivités
e1 = mate mr1 'EMIS' e_inf CABS k_abs 'TABS' tg;
e2 = mate mr2 'EMIS' e_sup CABS k_abs 'TABS' tg ;
el = mate mrl 'EMIS' e_late CABS k_abs 'TABS' tg;
e = e1 et e2 et el ;
*** Calcul des facteurs de forme
* facteurs de forme généralisés
* opti impi 1;
fft = ffor mrt e;
* opti impi 0;
*opti donn 5;
*** Conditions aux limites ...
c1 = bloq z1 'T' ;
tim1 = depi c1 T_inf ;
c2 = bloq z2 'T' ;
tim2 = depi c2 T_sup ;
*** méthode 1
* ---------
* calcul de la matrice de rayonnement
chamr = raye mrt fft e ;
*** Initialisation de la température ...
tp = manu chpo tout 1 'T' T_inf natu 'DIFFUS';
* calcul du flux associé au milieu absorbant
*<
Tg_abs = manu chpo tout 1 'T' Tg natu 'DIFFUS';
T_abs = redu (exco 'T' Tg_abs 'T') cavite ;
Te_abs = chan 'CHAM' T_abs mrt 'GRAVITE' ;
cr_abs = rayn mrt chamr Te_abs ;
Fg_abs = cr_abs * T_abs ;
* list Fg_abs;
*>
*** Résolution (par itérations) ...
* Coeff. de relaxation ...
alfa = 0.42 ;
nbpsup = nbno tout ;
listemp = prog ;
listres = prog ;
maxiter = 100 ;
critconv = 1.e-5 ;
REPE bloc1 maxiter ;
nbiter = &bloc1 ;
t_cavi = redu (exco 'T' tp 'T') cavite ;
te_cavi = chan 'CHAM' t_cavi mrt 'GRAVITE' ;
cr = rayn mrt chamr te_cavi ;
cndtot = cnd et cr et c1 et c2 ;
residu = (cndtot * tp) - Fg_abs ;
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 ;
tt1 = resou cndtot (tim1 et tim2 et Fg_abs) ;
dt = exco 'T' (tt1- tp) 'T' ;
normdt = ((xtx dt) / nbpsup) ** 0.5 ;
* mess ' La norme de delta t = ' normdt ;
tn = (alfa * tt1) + ((1.-alfa) * tp) ;
tp =tn ;
FIN bloc1 ;
*
*** Post-traitement ...
mess '1- calcul avec matrice de rayonnement';
* tsol = extr tt1 'T' q2 ;
* mess 'La solution obtenue = ' tsol ;
* flux rayonne
fray = (cr * tt1) - Fg_abs ;
fray_1 = maxi (resu fray);
mess ' bilan global cavite: ' fray_1 ;
* mess ' flux surface inf: ' (maxi (resu (redu fray l1))) ;
* mess ' flux surface sup: ' (maxi (resu (redu fray l2))) ;
* mess ' flux surface late: ' (maxi (resu (redu fray late))) ;
* flux associé à la condition de température imposée
* mess ' Timp surf. inf.: ' (maxi (resu (reac c1 tt1))) ;
* mess ' Timp surf. sup.: ' (maxi (resu (reac c2 tt1))) ;
t = exco 'T' tt1 'T' ;
titr 'Le champ de temperature final' ;
si(graph) ; trac tout t ; finsi ;
tlate = redu t late ; comm list tlate ;
*! trac (redu tt1 z3 ) z3 ;
ev1 = evol 'CHPO' (redu t l3) l3;
ev1 = ev1 coul 'VERT' ;
RESI1 =ABS((fray_1 - fref1)/fref1);
mess 'Erreur relative - methode 1' (100*RESI1) '%' ;
* opti donn 5;
* méthode 2
* ---------
*** Initialisation de la température ...
tp = manu chpo tout 1 'T' T_inf natu 'DIFFUS';
*** Résolution (par itérations) ...
* Coeff. de relaxation ...
alfa = 1.0 ;
nbpsup = nbno tout ;
listemp = prog ;
listres = prog ;
maxiter = 100 ;
critconv = 1.e-5 ;
REPE bloc2 ;
nbiter = &bloc2 ;
t_cavi = redu (exco 'T' tp 'T') cavite ;
te_cavi = chan 'CHAM' t_cavi mrt 'GRAVITE' ;
* calcul de la temperature de rayonnement associée: trad
trad = raye mrt fft e te_cavi 1e-7 ;
* calcul du coefficient d'echange
hrad = HRCAV mrt e te_cavi trad ;
trad_n1= chan 'CHPO' mrt trad ;
trad_n = nomc trad_n1 'T' 'NATU' 'DIFFUS';
* pour la condition de convection
cr = cond mrt hrad;
f = conv mrt hrad trad_n;
cndtot = cnd et c1 et c2 et cr ;
residu = (cndtot * tp) -f ;
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 bloc2 ;
finsi ;
si(nbiter > maxiter) ;
mess 'Erreur ! Pas de convergence après ' (nbiter-1)
' itérations !' ;
erre 2 ;
quitter bloc2 ;
finsi ;
listres = listres et (prog normres) ;
* mess '---------------------------------------' ;
* mess 'Itération N° ' &bloc2 ;
tt2 = resou cndtot (tim1 et tim2 et f) ;
dt = exco 'T' (tt2- tp) 'T' ;
normdt = ((xtx dt) / nbpsup) ** 0.5 ;
* mess ' La norme de delta t = ' normdt ;
tn = (alfa * tt2) + ((1.-alfa) * tp) ;
tp =tn ;
mess ' Température obtenue = ' (extr tt2 'T' q2) ;
FIN bloc2 ;
*opti donn 5;
*temps impr place sgac ;
*** Post-traitement ...
mess '2- calcul avec température de rayonnement';
t = exco 'T' tt2 'T' ;
titr 'Le champ de temperature final' ;
si(graph) ; trac tout t ; finsi ;
tlate = redu t late ;
* mess ' T laterale: min et max ' (mini tlate) (maxi tlate);
ev2 = evol 'CHPO' (redu t l3) l3;
ev2 = ev2 coul 'BLEU' ;
* flux rayonne
fray = (cr * tt2)- f;
fray_2 = maxi (resu fray);
mess ' bilan global cavite: ' fray_2 ;
* mess ' flux surface inf: ' (maxi (resu (redu fray l1))) ;
* mess ' flux surface sup: ' (maxi (resu (redu fray l2))) ;
* mess ' flux surface late: ' (maxi (resu (redu fray late))) ;
* flux associé à la condition de température imposée
* mess ' Timp surf. inf.: ' (maxi (resu (reac c1 tt2))) ;
* mess ' Timp surf. sup.: ' (maxi (resu (reac c2 tt2))) ;
si graph;
TAB1 = table;
TAB1. 'TITRE' = table;
TAB1. 'TITRE' . 1 = 'MOT' 'methode 1' ;
TAB1. 'TITRE' . 2 = 'MOT' 'methode 2' ;
dess (ev1 et ev2) mima lege TAB1 ;
finsi;
RESI2 =ABS((fray_2 - fref2)/fref2);
mess 'Erreur relative - methode 2' (100*RESI2) '%' ;
SI ((RESI1 <EG 1E-4) et (RESI2 <EG 1E-4));
ERRE 0;
SINO;
ERRE 5;
FINSI;
fin;
© Cast3M 2003 - All rights reserved.
Disclaimer