* fichier : rayoh-3D.dgibi
************************************************************************
* *
* Rayonnement thermique en milieu transparent: *
* *
* Test de bon fonctionnement de la procédure HRAYO permettant de *
* traiter le rayonnement face à face par comparaison à une solution * *
* analytique en régime permanent. *
* *
* Calcul 3D. *
* *
* 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 *
* Remarque: dans le cas quadratique, et pour le face à face chaque *
* frontière doit etre un maillage simple (TRI6 uniquement *
* ou QUA8 uniquement) *
************************************************************************
graph = faux ;
xquad = faux ;
si (xquad) ;
option echo 0 dime 3 elem pr15 ;
sinon ;
option echo 0 dime 3 elem cub8 ;
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
*-------------------------------
oeil=10 10 10;
p1=0 0 0 ; p2=1 0 0 ; p3=1 1 0 ; p4=0 1 0;
l1=p1 d 1 p2;
l2=p2 d 1 p3;
l3=p3 d 1 p4;
l4=p4 d 2 p1;
si xquad ;
* on recherche un objet simple
opti elem tri6 ;
s1= (l1 et l2 et l3 et l4) surf plan ;
option echo 1 dime 3 elem pr15 ;
sinon ;
s1= (l1 et l2 et l3 et l4) surf plan ;
finsi ;
n=1 ;
v = 0 0 0.1 ; v1 = (-1)*v ;
w = 0 0 0.001 ; w1 = (-1)*w ;
s2inf = s1 ;
s2sup = s2inf plus w ;
plaq2 = s2inf volu n s2sup ;
s1sup = s2inf plus v1;
s1inf = s1sup plus w1;
plaq1 = s1inf volu n s1sup ;
rac12 = liai 0.1 s1sup s2inf ;
tout = plaq1 et plaq2 ;
si graph ; trac oeil tout ; finsi ;
*-------------------------------
* Conduction
*-------------------------------
mcd = mode tout thermique ;
k = mate mcd 'K' lamb ;
cnd = cond mcd k ;
*-------------------------------
* Linéarisation du rayonnement :
* - 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_s1 = (s1sup chan TRI6) chan TRI3 ;
x_s2 = (s2inf chan TRI6) chan TRI3 ;
xx_s1 = x_s1 chan ligne;
xx_s2 = x_s2 chan ligne;
x_s12 = liai 0.1 x_s1 x_s2;
xx_s12 = x_s12 chan ligne ;
rel12 = diff xx_s12 (xx_s1 et xx_s2) ;
sinon ;
sr12 = rac12 chan ligne ;
sr1 = s1sup chan ligne ;
sr2 = s2inf chan ligne ;
rel12 = diff sr12 (sr1 et sr2) ;
finsi ;
si graph ; trac oeil rel12 ; 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 S2inf mcv12 rel12 CONS 'FAC1';
emi1 = mate mr1 'EMIS' e1 ;
mr2i = mode s2inf thermique rayonnement FAC_A_FAC
s1sup S2inf 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-2 ;
* 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 ;
* mr12 = matr mcv12 'H' h12 ;
* mr12 = h12 ;
* cr12 = cond mcv12 mr12 ;
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