* fichier : rayo-2D-2.dgibi
****************************************************************
*
* test 2D couplage conduction-rayonnement
* REFERENCE: cas TPNP 01/89 du guide VPCS
*
* DONNEES
* cavité carrée de cote 1.
*
* 1000K e=1.0
* -----
* e=0.5 | | e=0.5
* | |
* -----
* 2000K e=0.1
*
*
* RESULTATS
* temperature de la surface latérale : 1192 K (solution analytique)
* Remarques: option convexe dans FFOR et éléments quadratiques
*
* 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
* 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 ;
graph = faux ;
*** Solution de référence
tana = 1191.9 ;
*** Paramètres ...
l = 1.0 ;
* epaisseur des parois : dz = 10 mm
dz = 0.01 ;
dzp = dz ;
dzn = -1. * dz ;
* proprietes physiques
lam = 10. ;
e_sup=1.0 ; e_inf=0.1 ; e_late = 0.5 ;
T_sup = 1000. ; T_inf = 2000. ;
*** Points ...
dens l ;
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 1 p1 p2 ; l2 = d 1 p3 p4 ;
l3 = d 1 q2 q3 ; l4 = d 1 q4 q1 ;
late = l3 et l4 ;
* la surface laterale est disjointe des deux autres
cavite = l1 et l2 et late ;
z1 = l1 trans 1 (0. dzn) ;
z2 = l2 trans 1 (0. dzp) ;
z3 = l3 trans 1 (dzp 0.) ;
z4 = l4 trans 1 (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' 'CONVEXE' CONS 'CAV1';
mr2= modl l2 thermique rayonnement 'CAVITE' 'CONVEXE' CONS 'CAV1';
mrl= modl late thermique rayonnement 'CAVITE' 'CONVEXE' CONS 'CAV1';
mrt = mr1 et mr2 et mrl ;
*** Emissivités
e1 = mate mr1 'EMIS' e_inf ;
e2 = mate mr2 'EMIS' e_sup ;
el = mate mrl 'EMIS' e_late ;
e = e1 et e2 et el ;
*** Calcul des facteurs de forme
* opti 'IMPI' 1 ;
fft = ffor mrt e ;
* opti 'IMPI' 0 ;
*** 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';
*** Résolution (par itérations) ...
* Coeff. de relaxation ...
alfa = 0.42 ;
nbpsup = nbno tout ;
listemp = prog ;
listres = prog ;
maxiter = 100 ;
critconv = 1.e-5 ;
* opti echo 1 ;
REPE bloc1 ;
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 ;
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) ;
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 ;
*! mess ' Température obtenue = ' (extr tt1 'T' q2) ;
FIN bloc1 ;
* opti echo 1 ;
*** Post-traitement ...
mess '1- calcul avec matrice de rayonnement';
tsol = extr tt1 'T' q2 ;
mess 'La solution obtenue = ' tsol ;
* flux rayonne
fray = (cr * tt1);
fray_tot = resu fray;
mess ' bilan global cavite: ' (maxi fray_tot) ;
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 ; list tlate ;
RESI1 =ABS((tsol -tana)/tana);
mess 'Erreur relative - methode 1' (100*RESI1) '%' ;
mess ;
* 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 ;
* opti echo 1 ;
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
*opti impi 1;
trad = raye mrt fft e te_cavi 1e-7;
*opti impi 0;
*! mess 'trad ' (mini trad) (maxi trad);
* 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 ;
*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 ;
tsol = extr t 'T' q2 ;
mess 'La solution obtenue = ' tsol ;
* ev3 = evol 'CHPO' (redu t l3) l3;
* dess ev3 mima ;
* flux rayonne
fray = (cr * tt2)- f;
fray_tot = resu fray;
mess ' bilan global cavite: ' (maxi fray_tot) ;
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))) ;
RESI2=ABS((tsol -tana)/tana);
mess 'Erreur relative - méthode 2' (100*RESI2) '%' ;
SI ((RESI1 <EG 1E-4) et (RESI2 <EG 1E-4));
ERRE 0;
SINO;
ERRE 5;
FINSI;
*** Bye ...
FIN;
© Cast3M 2003 - All rights reserved.
Disclaimer