* fichier : rayoh-2D.dgibi ************************************************************************ * Section : Thermique Transitoire * Section : Thermique Rayonnement * Section : Thermique Convection ************************************************************************ ************************************************************************ * * * Rayonnement thermique en milieu transparent: * * * * Test du rayonnement face à face * * Comparaison avec une solution analytique en régime permanent. * * * * Calcul par 2 méthodes: * * 1. algorithme iteratif * * 2. utilisation de PASAPAS en transitoire jusqu'au permanent * * * * * * Calcul 2D-plan. * * * * on considère 3 plaques parallèles rayonnantes : * * * * -------- T3,e3 * * * * e2sup face à face * * -------- 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, et de meme pour celui entre les plaques 2 et 3 * * * * 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) ; *LIST T2_ANA ; *------------------------------- * 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 ; s3inf = s2inf plus (0. 0.1) ; plaq3 = s3inf tran 1 (0. 0.001) ; s3sup = inve ( cote 3 plaq3) ; rac23 = racc 0.1 s2Sup s3inf ; tout = plaq1 et plaq2 et plaq3 ; si graph ; trac tout ; finsi ; *------------------------------- * Conduction *------------------------------- mcd = mode tout thermique ; k = mate mcd 'K' lamb 'RHO' 1. 'C' LAMB; cnd = cond mcd k ; cte_sb = 5.673e-8 ; *------------------------------- * 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 ; mcv23 = mode rac23 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) ; x_s2sup = chan s2sup SEG2 ; x_s3inf = chan s3inf SEG2 ; x_rac23 = racc 0.1 x_s2sup x_s3inf ; rel23 = diff (x_rac23 chan ligne) (x_s2sup et x_s3inf) ; sinon; rel12 = diff (rac12 chan ligne) (s1sup et s2inf) ; rel23 = diff (rac23 chan ligne) (s2sup et s3inf) ; finsi; *------------------------------- * Modeles de rayonnement et emissivites *------------------------------- * face a 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 ; mr2S= mode s2sup thermique rayonnement FAC_A_FAC s2sup s3inf mcv23 rel23 cons 'FAC2'; emi2S= mate mr2S 'EMIS' e2SUP ; mr3 = mode s3inf thermique rayonnement FAC_A_FAC s2sup s3inf mcv23 rel23 cons 'FAC2'; emi3 = mate mr3 'EMIS' e3 ; *------------------------------- * Conditions aux limites *------------------------------- c1 = bloq plaq1 'T' ; tim1 = depi c1 T1 ; c3 = bloq plaq3 'T' ; tim3 = depi c3 T3 ; lr1 = prog 0 100 ; lr2 = prog 1 1 ; ev1 = evol manu 't' lr1 'T' lr2 ; cht1 = char 'TIMP' ev1 tIM1 ; ev3 = evol manu 't' lr1 'T' lr2 ; cht3 = char 'TIMP' ev3 tIM3 ; *** Températures initiales ... T2 = 500. ; chINI1 = manu chpo PLAQ1 1 'T' T1 nature diffus ; chINI2 = manu chpo PLAQ2 1 'T' T2 nature diffus ; chINI3 = manu chpo PLAQ3 1 'T' T3 nature diffus ; *---------------------------------------------------------- * 1. Première méthode de résolution : Algorithme itératif *---------------------------------------------------------- *** Initialisation de la température ... tp = CHINI1 ET CHINI2 ET CHINI3 ; 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 ; * face à face (2,3) t2sup = redu (exco 'T' tp 'T') s2sup; t3inf = redu (exco 'T' tp 'T') s3inf; h23 = HRAYO mcv23 mr2S emi2S t2sup mr3 emi3 t3INF rel23 ; cr23 = cond mcv23 h23 ; crtot= cr12 et cr23 ; cndtot = crtot et cnd et c1 ET C3; 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 ; tt = resou cndtot (tim1 et TIM3) ; 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 ... *------------------------------- Temp = exco 'T' tt 'T' ; T2A = redu plaq2 Temp ; RESIA=maxi( ABS((T2A - T2_ana)/T2_ana)); *------------------------------------------------------ * 2. Deuxieme méthode de résolution : appel a PASAPAS *------------------------------------------------------ *** Préparation de la table pour PASAPAS ... tabnl = table ; tabnl . MODELE = MCD ET MR1 ET MR2I ET MR2S ET MR3 ; tabnl . CARACTERISTIQUES = K ET EMI1 ET EMI2I ET EMI2S ET EMI3; tabnl . BLOCAGES_THERMIQUES = C1 ET C3 ; tabnl . CHARGEMENT = cht1 et cht3 ; tabnl . TEMPS_CALCULES = prog 0. pas 0.1 1. pas 0.2 4. ; tabnl . TEMPERATURES = table ; tabnl . TEMPERATURES . 0 = CHINI1 ET CHINI2 ET CHINI3 ; tabnl . 'PROCEDURE_THERMIQUE' = 'DUPONT' ; tabnl . 'CTE_STEFAN_BOLTZMANN' = cte_sb ; tabnl . 'RELAXATION_THETA' = 1. ; *** Appel à PASAPAS ... pasapas tabnl ; *** Petit post-traitement ... listt = prog ; listtp3 = prog ; listtp3r = prog ; nbpas = dime (tabnl . TEMPS) ; repeter surpas nbpas ; lindice = &surpas - 1 ; listt = listt et (prog (tabnl . TEMPS . lindice)) ; valtp3 = MAXI (REDU tabnl . TEMPERATURES . lindice PLAQ2) ; listtp3 = listtp3 et (prog valtp3) ; fin surpas ; titr 'Evolution de la temperature de la plaque2'; evtp3 = evol manu 't' listt 'T' listtp3 ; si(graph) ; dess evtp3 ; finsi ; T2B = redu tabnl . TEMPERATURES . (nbpas - 1) PLAQ2 ; RESIB=maxi( ABS((T2B - T2_ana)/T2_ana)); * comparaison a la solution analytique mess 'Solution analytique ' T2_ana ; mess 'Solution Algorithme ' (MAXI T2A) ; mess 'Erreur relative ' (100*RESIA) '%' ; mess 'Solution PASAPAS ' (MAXI T2B) ; mess 'Erreur relative ' (100*RESIB) '%' ; SI ((RESIA