Télécharger ray.dgibi

Retour à la liste

* 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=SOMM( 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=SOMM( 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 - Tous droits réservés.
Mentions légales