* fichier : rayo-axi-2.dgibi ************************************************************************ * Section : Thermique Rayonnement * Section : Thermique Convection * Section : Thermique Conduction ************************************************************************ * pour calcul complet mettre complet à : vrai; complet = faux; ************************************************************************ * Rayonnement thermique en milieu transparent: * * * * vérification du bon fonctionnement de l'opérateur FFOR dans le cas * * général (traitement des parties cachées) sur un cas thermiquement * * simple. * * * * Test 2D-axisymétrique sur un problème à symétrie sphérique 1D: * * ------------------ * * on considère 3 sphères concentriques de rayons 25mm, 50mm et 100mm * * conductrices et d'emissivité 1. La température des sphères intérieure* * et extérieure est uniforme (respectivement 1000K et 2000K) et on * * détermine la température d'équilibre de la sphère intermédiaire. * * (conventions: indices 1,2,3 pour les sphères de rayon 25,50,100mm) * * * * Jeu de données identique au cas 2D-plan :rayo-2D-3.dgibi * * hormis la valeur des 2 paramètres: * * 'xrev' vrai : géométrie axisymétrique * * faux : 2D plane * * 'x2cav' vrai : on considère les 2 cavités séparément * * faux : l'ensemble des 2 cavités * * * ************************************************************************ xrev = vrai ; x2cav = faux ; * attention ne fonctionne peut etre pas si x2cav= vrai et xrev=vrai * il faudrait refaire le model sans "symetrie .." graph = faux ; si xrev ; option echo 0 dime 2 elem qua4 mode axis ; sinon; option echo 0 dime 2 elem qua4 ; finsi; *------------------------------- * Données physiques du problème *------------------------------- r1 = 25e-3 ; r2 = 50e-3 ; r3 = 100e-3 ; T1 = 1000. ; T3 = 2000. ; emis = 1. ; lamb = 20. ; *------------------------------- * Solution analytique *------------------------------- si xrev ; f21 = (r1*r1)/(r2*r2) ; sinon; f21 = r1/r2 ; finsi; f22 = 1. - f21 ; M1 = T1**4. ; M3 = T3**4. ; M2 = ( (f21*M1) + M3 )/(2.-f22) ; T2_ana = M2**(0.25) ; *------------------------------- * Maillage *------------------------------- C = 0. 0. ; * sphère 1 p1_sud = 0. ( -1. * r1 ) ; p1_equa = r1 0. ; p1_nord = 0. r1 ; d1 = r1/3 ; c1_inf = cerc p1_sud C p1_equa DINI d1 DFIN d1; c1_sup = cerc p1_equa C p1_nord DINI d1 DFIN d1; c1_tot = c1_inf et c1_sup ; c1_int = c1_tot homo C 0.999 ; c1_ext = c1_tot homo C 1.001 ; sphere1 = regl c1_int c1_ext ; * sphère 2 p2_sud = 0. ( -1. * r2 ) ; p2_equa = r2 0. ; p2_nord = 0. r2 ; d2 = r2/3 ; c2_inf = cerc p2_sud C p2_equa DINI d2 DFIN d2; c2_sup = cerc p2_equa C p2_nord DINI d2 DFIN d2; c2_tot = c2_inf et c2_sup ; c2_int = c2_tot homo C 0.999 ; c2_ext = c2_tot homo C 1.001 ; sphere2 = regl c2_int c2_ext ; * sphère 3 p3_sud = 0. ( -1. * r3 ) ; p3_equa = r3 0. ; p3_nord = 0. r3 ; si complet; d3 = r3/3 ; sinon; d3 = r3 / 1.5; finsi; c3_inf = cerc p3_sud C p3_equa DINI d3 DFIN d3; c3_sup = cerc p3_equa C p3_nord DINI d3 DFIN d3; c3_tot = c3_inf et c3_sup ; c3_int = c3_tot homo C 0.999 ; c3_ext = c3_tot homo C 1.001 ; sphere3 = regl c3_int c3_ext ; * on controle l'orientation des elements pour les cavités rayonnantes cavi_12 = (inve c1_ext) et c2_int ; cavi_23 = (inve c2_ext) et c3_int ; cavi_tot = cavi_12 et cavi_23 ; tout = sphere1 et sphere2 et sphere3 ; si graph ; trac tout ; finsi ; *------------------------------- * Conduction *------------------------------- mcd = mode tout thermique ; k = mate mcd 'K' lamb ; cnd = cond mcd k ; *------------------------------- * Modèles de rayonnement *------------------------------- si x2cav; mr12 = mode cavi_12 thermique rayonnement 'CAVITE' SYMETRIE C (0. 1.) 'CONS' 'CAV1' ; e12 = mate mr12 'EMIS' emis ; mr23 = mode cavi_23 thermique rayonnement 'CAVITE' SYMETRIE C (0. 1.) 'CONS' 'CAV1' ; e23 = mate mr23 'EMIS' emis ; sinon; mrtot= mode cavi_tot thermique rayonnement 'CAVITE' SYMETRIE C (0. 1.) ; etot = mate mrtot 'EMIS' emis ; finsi; *------------------------------- * Facteurs de forme et matrices de rayonnement *------------------------------- si x2cav ; *opti 'IMPI' 1 ; si xrev ; ff12 = ffor mr12 e12; ff23 = ffor mr23 e23; sinon; ff12 = ffor mr12 e12 ; ff23 = ffor mr23 e23 ; finsi; *opti 'IMPI' 0 ; cham12 = raye mr12 ff12 e12 ; cham23 = raye mr23 ff23 e23 ; sinon; *opti 'IMPI' 1 ; si xrev ; fftot= ffor mrtot etot; sinon; fftot = ffor mrtot etot ; finsi; *opti 'IMPI' 0 ; chamtot = raye mrtot fftot etot ; finsi; *------------------------------- * Conditions aux limites *------------------------------- c1 = bloq sphere1 'T' ; tim1 = depi c1 T1 ; c3 = bloq sphere3 'T' ; tim3 = depi c3 T3 ; *------------------------------- * 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 ; si complet; critconv = 1.e-5 ; sinon; critconv = 1.; finsi; * opti echo 1 ; REPE bloc1 ; nbiter = &bloc1 ; * traitement du rayonnement si x2cav ; t_c12 = redu (exco 'T' tp 'T') cavi_12 ; te_c12 = chan 'CHAM' t_c12 mr12 'GRAVITE' ; cr12 = rayn mr12 cham12 te_c12 ; t_c23 = redu (exco 'T' tp 'T') cavi_23 ; te_c23 = chan 'CHAM' t_c23 mr23 'GRAVITE' ; cr23 = rayn mr23 cham23 te_c23 ; crtot = cr12 et cr23 ; sinon; t_tot = redu (exco 'T' tp 'T') cavi_tot; te_tot = chan 'CHAM' t_tot mrtot 'GRAVITE' ; crtot= rayn mrtot chamtot te_tot ; finsi; 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 ... *------------------------------- titre 'Champ de température final' ; si(graph) ; trac sphere2 (exco 'T' tn) ; finsi ; Temp = exco 'T' tt 'T' ; T2 = redu sphere2 Temp ; RESI=maxi( ABS((T2 - T2_ana)/T2_ana)); mess 'Problème à symétrie sphérique'; mess 'Solution analytique ' T2_ana ; mess 'Erreur relative ' (100*RESI) '%' ; SI (RESI