* fichier : ray.dgibi
***********************************************************************
* CAVITE RECTANGULAIRE - Couplage Convection Naturelle_Rayonnement *
* DEN/DM2S/SFME/LTMF Juillet2002 *
* A.Stietel / G.Forestier *
* (méthode 1: avec calcul de la matride de rayonnement) *
* *
* 01/2006: mise au point d'une nouvelle méthode pour le rayonnement *
* (méthode 2) *
* le flux rayonné est mis sous la forme Phi=emis*sigma*(T**4-Trad**4) *
* - Trad est évalué par l'opérateur RAYE (nouvelle option) *
* en fonction des facteurs de forme, de l'emissivité et de la *
* température à une itération donnée *
* - la linéarisation se fait dans la procédure HRCAVNS *
* - l'opérateur ECHI traite ensuite la condition aux limites *
***********************************************************************
COMPLET = faux ;
SI ( COMPLET ) ;
ni = 300 ; itt = 1. ; ome = 0.3 ;
nn = 12 ;
SINON ;
ni = 30 ; itt = 1. ; ome = 0.3 ;
nn = 6 ;
FINSI ;
GRAPH = faux ;
opti dime 2 elem qua8 ;
opti isov suli ;
******************
* MAILLAGE *
******************
l = 5.; h = 1.;
si ( l >eg h) ;
ny = nn ; nx = ny * l / h ; nx = enti nx ;
sinon ;
nx = nn ; ny = nx * h / l ; ny = enti ny ;
finsi ;
ii = 3. ;ij = 1. / ii ;
d1 = l / nx / 2 * ij ; d2 = l / nx / 2 * ii;
e1 = h / ny / 2 * ij ; e2 = h / ny / 2 * ii ;
discr = quaf ;
discr = macro ;
epsm = 0.0001 ;
p1 = 0. 0. ; p2 = l 0. ; p3 = l h ; p4 = 0 h ;
pb = (l/2) 0. ; ph = (l/2) h ; pd = l (h/2) ; pg = 0 (h/2) ;
pm = (l/2) (h/2) ;
bas = p1 d dini d1 dfin d2 pb d dini d2 dfin d1 p2 ;
hau = p3 d dini d1 dfin d2 ph d dini d2 dfin d1 p4 ;
gau = p4 d dini e1 dfin e2 pg d dini e2 dfin e1 p1 ;
dro = p2 d dini e1 dfin e2 pd d dini e2 dfin e1 p3 ;
cenv = pb d dini e1 dfin e2 pm d dini e2 dfin e1 ph ;
miv = pg d dini d1 dfin d2 pm d dini d2 dfin d1 pd ;
cav = bas et dro et hau et gau ;
dom = dall plan bas dro hau gau ;
*trac (cav et miv);
domf = chan dom quaf ;
cavf = chan cav quaf ;
elim (domf et cavf et cenv et miv ) epsm ;
*trac dom ;
mdom = mode domf navier_stokes discr ;
mcav = mode cavf navier_stokes discr ;
qav = doma mcav maillage;
m_sray = mcav ;
********************
* PARAMETRES *
********************
t0 = 1000. ; rapt = 0.7 ;
tc = t0 ; tf = t0 * rapt ;
re = 1.43 ;
pr = 0.7 ;
gr = 1.e5 ;
ttm = tc + tf * 0.5 ;
tm = 0.5 ;
gb = 0. (-1.*gr /re/re) ;
nu = 1./re ;
pe = re * pr ;
nrf = 3. ;
sig = 5.67e-8 ;
kko = pe * sig * t0 *t0 *t0 /nrf ;
kla = sig * t0 *t0 *t0 /nrf ;
*mess 'kko ' kko ' kla ' kla ;
****************************
* DONNEES NUMERIQUES *
****************************
* création d'une table pour le rayonnement en cavité
mrm = mode qav thermique rayonnement cavite;
marm = mate mrm emis 1. ;
tabra = table ;
tabra . ma_rai = qav ;
tabra . mm_rai = mrm ;
tabra . mm_ns = m_sray;
tabra . mt_rai = marm ;
* calcul des facteurs de forme
tabra . mf_rai = ffor mrm marm ;
* calcul de la matrice de rayonnement
* cette étape n'est nécessaire que pour la méthode 1
tabra . mr_rai = raye mrm (tabra . mf_rai) marm ;
*********************
* CONVERGENCE *
*********************
debproc calcul ;
argu rvx*table ;
rt = rvx . eqex ;
jt = dime (rt . inco . it ) ;
jt = jt + 1 ;
unn = rt . inco . un ;
un1 = rt . inco . u1 ;
tnn = rt . inco . tn ;
tn1 = rt . inco . t1 ;
mdo = rt . 1calcul . domz ;
mma = doma mdo maillage ;
unx = kcht mdo scal sommet (exco ux unn) ;
uny = kcht mdo scal sommet (exco uy unn) ;
u1x = kcht mdo scal sommet (exco ux un1) ;
u1y = kcht mdo scal sommet (exco uy un1) ;
rt . inco . u1 = kcht mdo vect sommet unn ;
rt . inco . t1 = kcht mdo scal sommet tnn ;
rt . inco . tet = kcht mdo scal sommet ((tnn - tf)/(tc - tf)) ;
erux = kops unx - u1x ;
elix = maxi erux abs ;
elix = (log(elix + 1.e-10))/( log 10.) ;
eruy = kops uny - u1y ;
eliy = maxi eruy abs ;
eliy = (log(eliy + 1.e-10))/( log 10.) ;
ertt = kops tnn - tn1 ;
elit = maxi ertt abs ;
elit = (log(elit + 1.e-10))/( log 10.) ;
err = maxi (prog elix eliy) ;
mess 'err' elix eliy elit jt ;
rt . inco . it = (rt . inco . it) et (prog jt) ;
rt . inco . erx= (rt . inco . erx) et (prog elix);
rt . inco . ery= (rt . inco . ery) et (prog eliy);
rt . inco . ert= (rt . inco . ert) et (prog elit);
as af = kops matrik ;
finproc as af ;
*-------------------------------------------------------------
'DEBPROC' HRCAVNS;
*---------------------------PROCEDURE RAYONNEMENT--------------------------*
****************************************************************************
* Procédure permettant de prendre en compte le rayonnement dans *
* les operateurs eqex et exec (couplage convection-rayonnement) *
* dans l'environnement du modèle Navier-Stokes *
* Calcule le champ de température Trad et le coefficient d'echange Hrad *
* correspondant à une linéarisation du rayonnement dans une cavité *
* qui sont utlisés par l'opérateur ECHI s'appuyant sur la cavité *
* *
* (issu de la procédure RAY) *
****************************************************************************
**** création d'une table ****
'ARGU' rvx * 'TABLE' ;
rv = rvx . 'EQEX' ;
**** récuperation du champ par points des températures ****
line =rvx.'LISTINCO';
noinco = 'EXTR' line 1 ;
t = rv.'INCO'.noinco;
tray = rvx . 'ARG1';
**** traitements ****
*maillage rayonnement*
cavi = tray . 'MA_RAI';
*modele rayonnement*
mrt = tray . 'MM_RAI' ;
*modele navier-stokes*
m_ns = tray . 'MM_NS' ;
*emissivite
emis = tray . 'MT_RAI' ;
*facteurs de forme
ff = tray . 'MF_RAI' ;
mrn = rvx . 'DOMZ' ;
* mess ' ';
t_cavi = 'REDU' t cavi ;
* mess ' t_cavi ' (mini t_cavi) (maxi t_cavi);
tc = 'EXCO' t_cavi scal 'T';
tre = 'CHAN' 'CHAM' tc mrt 'GRAVITE';
* opti impi -1;
* trad = 'RAYE' mrt ff emis tre 1e-8;
trad = 'RAYE' mrt ff emis tre ;
* opti impi 0;
mess ' TRAD: ' (mini trad) (maxi trad);
*
* transformation des CHAMELEM en scal centre pour C_fluide
* 1/ temperature resultat de RAYE
*
trad = 'NOMC' trad 'SCAL';
trad_n = 'CHAN' 'CHPO' mrt trad ;
rv . inco . 'TRAD' = 'KCHT' m_ns scal sommet trad_n;
trad = 'NOEL' m_ns (rv . inco . 'TRAD');
* 2/ émissivité
em_1 = 'NOMC' emis 'SCAL';
em_n = 'CHAN' 'CHPO' mrt em_1 ;
em_s = 'NOEL' m_ns em_n;
te_cavi = 'NOEL' m_ns t_cavi;
*
* calcul du coefficient d'echange linéarisé
*
sig = 5.67e-8;
te2 = te_cavi * te_cavi;
trad2 = trad*trad;
hrad = em_s * sig * (te2 + trad2) * (te_cavi + trad);
mess ' HRAD: ' (mini hrad) (maxi hrad);
rv . inco . 'HRAD' = 'KCHT' m_ns scal centre hrad;
as af = 'KOPS' 'MATRIK' ;
'FINP' as af ;
*---------------------FIN PROCEDURE RAYONNEMENT--------------------------*
*-------------------------------------------------------------
pip = elem (doma mdom centrep1) (lect 1) ;
xx = coor 1 dom ;
**************************
* SCHEMA IMPLICITE *
**************************
* méthode 1
* ---------
*-------------------------------------------------------------
rw = eqex mdom niter ni omega ome
zone mdom oper calcul
opti ef impl centrep1
zone mdom oper kbbt 1. inco un pn
zone mdom oper ns 1. 'UN' nu gb tet tm inco un
zone mdom oper konv kko un 1. inco tn
zone mdom oper lapn kla inco tn
;
rw = eqex rw
opti ef impl centree
zone m_sray oper ray tabra inco tn
;
rw = eqex rw
clim un uimp cav 0.
un vimp cav 0.
tn timp gau tc
tn timp dro tf
pn timp pip 0.
;
rw . inco = table inco ;
rw . inco . tabray = tabra ;
rw . inco . un = kcht mdom vect sommet (0. 0.) ;
rw . inco . u1 = kcht mdom vect sommet (0. 0.) ;
rw . inco . tn = kcht mdom scal sommet (tf - tc * xx / l + tc ) ;
rw . inco . t1 = kcht mdom scal sommet (tf - tc * xx / 1 + tc ) ;
rw . inco . pn = kcht mdom scal centrep1 0. ;
rw . inco . dt = dtt ;
rw . inco . it = prog ;
rw . inco . erx = prog ;
rw . inco . ery = prog ;
rw . inco . ert = prog ;
*-------------------------------------------------------------
* méthode 2
* ---------
rw2= eqex mdom niter ni omega ome
zone mdom oper calcul
opti ef impl centrep1
zone mdom oper kbbt 1. inco un pn
zone mdom oper ns 1. 'UN' nu gb tet tm inco un
zone mdom oper konv kko un 1. inco tn
zone mdom oper lapn kla inco tn
;
rw2 = eqex rw2
opti ef impl centree
zone m_sray oper hrcavns tabra inco tn
zone m_sray oper echi 'HRAD' 'TRAD' inco tn
;
rw2 = eqex rw2
clim un uimp cav 0.
un vimp cav 0.
tn timp gau tc
tn timp dro tf
pn timp pip 0.
;
rw2 . inco = table inco ;
rw2 . inco . tabray = tabra ;
rw2 . inco . un = kcht mdom vect sommet (0. 0.) ;
rw2 . inco . u1 = kcht mdom vect sommet (0. 0.) ;
rw2 . inco . tn = kcht mdom scal sommet (tf - tc * xx / l + tc ) ;
rw2 . inco . t1 = kcht mdom scal sommet (tf - tc * xx / 1 + tc ) ;
rw2 . inco . pn = kcht mdom scal centrep1 0. ;
rw2 . inco . dt = dtt ;
rw2 . inco . it = prog ;
rw2 . inco . erx = prog ;
rw2 . inco . ery = prog ;
rw2 . inco . ert = prog ;
*-------------------------------------------------------------
*******************
* EXECUTION *
*******************
* titr 'impl macro' ;
exec rw ;
exec rw2;
***************
* TESTS *
***************
* valeurs de référence pour le calcul non-complet
Lref = prog
7.0000E+2 7.1468E+2 7.2605E+2 7.3691E+2 7.4622E+2 7.5505E+2
7.6300E+2 7.7055E+2 7.7677E+2 7.8233E+2 7.8694E+2 7.9101E+2
7.9430E+2 7.9721E+2 7.9945E+2 8.0133E+2 8.0291E+2 8.0452E+2
8.0590E+2 8.0730E+2 8.0860E+2 8.0994E+2 8.1144E+2 8.1316E+2
8.1474E+2 8.1636E+2 8.1798E+2 8.1977E+2 8.2156E+2 8.2345E+2
8.2524E+2 8.2727E+2 8.2918E+2 8.3129E+2 8.3340E+2 8.3553E+2
8.3775E+2 8.4006E+2 8.4233E+2 8.4489E+2 ;
* méthode 1
* ---------
tt = rw . inco . tn ;
Ttest = evol chpo hau tt;
Ttest = extr Ttest ordo ;
Ttest = extr Ttest (lect 1 pas 1 40);
list Ttest;
*opti sauv 'FORMAT' 'resu' ;
*sauv 'FORMAT' Ttest;
ERt=resu( abs ( Lref - Ttest)) ;
ERm=maxi( abs ( Lref - Ttest)) ;
mess ' méthode 1- Ecart sur profil de température ' ert erm;
* méthode 2
* ---------
tt = rw2. inco . tn ;
Ttest = evol chpo hau tt;
Ttest = extr Ttest ordo ;
Ttest = extr Ttest (lect 1 pas 1 40);
list Ttest;
*opti sauv 'FORMAT' 'resu' ;
*sauv 'FORMAT' Ttest;
ERt2=resu( abs ( Lref - Ttest)) ;
ERm2=maxi( abs ( Lref - Ttest)) ;
mess ' méthode 2 - Ecart sur profil de température ' ert2 erm2;
SI ((ERm < 0.1) et (ERm2 < 0.7));
ERREUR 0 ;
SINON;
ERREUR 5 ;
FINSI;
****************************************
* PROCEDURE (TRACEE DES COURBES) *
****************************************
si graph ;
schema= 'conv-rayo';
evx = evol bleu manu (rw . inco . it) (rw . inco . erx) ;
evy = evol vert manu (rw . inco . it) (rw . inco . ery) ;
evt = evol roug manu (rw . inco . it) (rw . inco . ery) ;
titr 'Test de convergence, ' schema;
dess (evx et evy et evt) ;
tt = rw . inco . tn ;
tr = ( tt - tf ) / ( tc - tf ) ;
titr 'Isothermes, ' schema;
trac 50 tr dom cav ;
vn = rw . inco . un ;
titr 'Vitesses, ' schema;
vv = vect vn (1.e2/gr)ux uy roug ;
titr 'Vitesses, ' schema;
trac vv dom cav ;
evx= EXCO vn ux;
evy= EXCO vn uy;
evux = evol rouge chpo evx cenv ;
evuy = evol bleu chpo evy cenv ;
titr 'Vitesses sur axe médian des y,'schema;
dess (evux et evuy);
evux = evol rouge chpo evx miv ;
evuy = evol bleu chpo evy miv ;
titr 'Vitesses sur axe médian des x,'schema;
dess (evux et evuy);
grat = kops tt grad mdom ;
gratx = kcht mdom scal centre (exco ux grat) ;
gratx = elno mdom gratx ;
gratx = -1. * gratx ;
gratxd = redu dro gratx ;
evd = evol chpo gratxd scal dro ;
gratxg = redu gau gratx ;
evg = evol chpo gratxg scal gau ;
titr ' gradient de température laterale, ' schema;
dess (evd et evg) ;
tb = redu bas tt ;
tb = redu bas tr ;
etb = evol vert chpo tb scal bas ;
ho = inve hau ;
th = redu ho tt ;
th = redu ho tr ;
eth = evol rouge chpo th scal ho ;
titr 'Evolution température sur parois supérieure et
inférieure,'schema;
dess (etb et eth) ;
vv = rw . inco . un ;
sw = kops vn rot mdom ;
rk = eqex mdom opti ef impl
zone mdom oper lapn 1. inco psi
zone mdom oper fimp sw inco psi
clim psi timp cav 0. ;
rk . inco . psi = kcht mdom scal sommet 0. ;
exec rk ;
trac (rk . inco . psi) dom cav 14 ;
finsi;
**********************************
* SAUVEGARDE DES RESULTATS *
**********************************
*opti sauv '/test4/ttmf4/rectangulaire/implrayo_proc.sav' ;
*sauv ;
fin ;
© Cast3M 2003 - All rights reserved.
Disclaimer