* fichier : rayoh-2D.dgibi ************************************************************************ * Section : Thermique Rayonnement * Section : Thermique Convection * Section : Thermique Conduction ************************************************************************ ************************************************************************ * * * 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 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 mr2s = mode s2sup thermique rayonnement INFINI CONS 'INF1'; emi2s= mate mr2s 'EMIS' e2sup 'E_IN' e3 ; t3inf = manu chpo s2sup 1 'T' T3 ; *------------------------------- * 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 t3inf ; * mr23 = MATE 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 ; RESI=maxi( ABS((T2 - T2_ana)/T2_ana)); mess 'Solution analytique ' T2_ana ; mess 'Erreur relative ' (100*RESI) '%' ; SI (RESI