* fichier : rayoh-2D.dgibi
************************************************************************
* *
* Rayonnement thermique en milieu transparent: *
* *
* Test de bon fonctionnement de la procédure HRAYO permettant de *
* traiter le rayonnement face à face ou avec un milieu infini. * *
* Comparaison avec une solution analytique en régime permanent. *
* *
* Calcul 2D-plan. *
* *
* on considère 3 plaques parallèles rayonnantes : *
* *
* -------- T3,e3 milieu infini *
* *
* e2sup *
* -------- T2 *
* e2inf face à face *
* *
* -------- T1,e1 *
* *
* La température des plaques supérieure et inférieure est imposée et *
* uniforme (respectivement 1000K et 2000K) et on détermine la tempéra- *
* ture d'équilibre de la plaque intermédiaire. *
* Le rayonnement entre les plaques 1 et 2 est modélisé par un rayonne- *
* ment face à face, celui entre les plaques 2 et 3 par un rayonnement *
* avec un milieu infini, le milieu infini étant representé par la *
* plaque 3 non maillée. *
* *
* paramètre xquad vrai : éléments quadratiques *
* faux : éléments linéaires *
************************************************************************
graph = faux ;
xquad = vrai ;
si xquad ;
option echo 1 dime 2 elem qua8 ;
sinon ;
option echo 1 dime 2 elem qua4 ;
finsi;
*-------------------------------
* Données physiques du problème
*-------------------------------
T1 = 1000. ; T3 = 2000. ;
e1 = 0.1 ; e2inf = 0.6 ; e2sup = 0.3 ; e3 = 0.9 ;
lamb = 2.e3 ;
*-------------------------------
* Solution analytique
*-------------------------------
r1 = 1. - e1 ; r2i = 1. - e2inf ;
r3 = 1. - e3 ; r2s = 1. - e2sup ;
alpha1= (e1*e2inf) / (1.-(r1*r2i)) ;
alpha3= (e3*e2sup) / (1.-(r3*r2s)) ;
M1 = T1**4. ; M3 = T3**4. ;
M2 = ((alpha1*M1) + (alpha3*M3)) / (alpha1+alpha3) ;
T2_ana = M2**(0.25) ;
*-------------------------------
* Maillage
*-------------------------------
p = 0. 0. ; q = 1. 0. ; s = droi 2 p q ;
*s1inf = s tour (0.5 0.) 30. ;
s1inf = s ;
plaq1 = s1inf tran 1 (0. 0.001) ;
s1sup = inve ( cote 3 plaq1) ;
s2inf = s1inf plus (0. 0.1) ;
plaq2 = s2inf tran 1 (0. 0.001) ;
s2sup = inve ( cote 3 plaq2) ;
rac12 = racc 0.1 s1sup s2inf ;
tout = plaq1 et plaq2 ;
si graph ; trac tout ; finsi ;
*-------------------------------
* Conduction
*-------------------------------
mcd = mode tout thermique ;
k = mate mcd 'K' lamb ;
cnd = cond mcd k ;
*-------------------------------
* Linéarisation du rayonnement :
* - définition de modèles de convection
* - creation des SEG2 définissant les relations entre les supports
* des champs définis sur s1sup et s2inf (pour le face à face)
*-------------------------------
* face à face
mcv12 = mode rac12 thermique convection ;
si xquad ;
x_s1sup = chan s1sup SEG2 ;
x_s2inf = chan s2inf SEG2 ;
x_rac12 = racc 0.1 x_s1sup x_s2inf ;
rel12 = diff (x_rac12 chan ligne) (x_s1sup et x_s2inf) ;
sinon;
rel12 = diff (rac12 chan ligne) (s1sup et s2inf) ;
finsi;
* milieu infini
mcv2s = mode s2sup thermique convection ;
*-------------------------------
* Mod�les de rayonnement et �missivit�s
*-------------------------------
* face � face
mr1 = mode s1sup thermique rayonnement FAC_A_FAC
s1sup s2sup mcv12 rel12 cons 'FAC1';
emi1 = mate mr1 'EMIS' e1 ;
mr2i = mode s2inf thermique rayonnement FAC_A_FAC
s1sup s2sup mcv12 rel12 cons 'FAC1';
emi2i= mate mr2i 'EMIS' e2inf ;
* milieu infini
mr3 = mode s2sup THERMIQUE rayonnement INFINI CONS 'INF1';
emi3 = mate mr3 'EMIS' e3 ;
t3inf = manu chpo s2sup 1 'T' T3 ;
mr2s = mode s2sup thermique rayonnement INFINI CONS 'INF1';
emi2s= mate mr2s 'EMIS' e2sup ;
*-------------------------------
* Conditions aux limites
*-------------------------------
c1 = bloq plaq1 'T' ;
tim1 = depi c1 T1 ;
*-------------------------------
* 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:
* face à face (1,2)
t1sup = redu (exco 'T' tp 'T') s1sup;
t2inf = redu (exco 'T' tp 'T') s2inf;
h12 = HRAYO mcv12 mr1 emi1 t1sup mr2i emi2i t2inf rel12 ;
cr12 = cond mcv12 h12 ;
* milieu infini (2,3)
t2sup = redu (exco 'T' tp 'T') s2sup;
h23 = HRAYO mcv2s mr2s emi2s t2sup mr3 emi3 t3inf ;
* mr23 = matr mcv2s 'H' h23 ;
cr23 = cond mcv2s h23 ;
f3 = conv mcv2s h23 t3inf ;
crtot= cr12 et cr23 ;
cndtot = crtot et cnd et c1 ;
residu = (cndtot * tp) - f3 ;
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 f3) ;
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 plaq2 (exco 'T' tn) ; finsi ;
Temp = exco 'T' tt 'T' ;
T2 = redu plaq2 Temp ;
list T2 ;
RESI=maxi( ABS((T2 - T2_ana)/T2_ana));
mess 'Solution analytique ' T2_ana ;
mess 'Erreur relative ' (100*RESI) '%' ;
SI (RESI <EG 1E-4);
ERRE 0;
SINO;
ERRE 5;
FINSI;
fin;
© Cast3M 2003 - All rights reserved.
Disclaimer