* fichier : rayo_abs-2D-1.dgibi ************************************************************************ * Section : Thermique Diffusion * Section : Thermique Rayonnement ************************************************************************ complet=faux; **************************************************************** * * * 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 = modeLI tout thermique ; k = mate mcd 'K' lamb ; cnd = cond mcd k ; * mr1= modeli l1 thermique rayonnement 'CAVITE' 'CONS' 'TROU1'; mr2= modeli l2 thermique rayonnement 'CAVITE' 'CONS' 'TROU1'; mrl= modeli 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