Télécharger ray.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : ray.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4.  
  5. ***********************************************************************
  6. * CAVITE RECTANGULAIRE - Couplage Convection Naturelle_Rayonnement *
  7. * DEN/DM2S/SFME/LTMF Juillet2002 *
  8. * A.Stietel / G.Forestier *
  9. * (méthode 1: avec calcul de la matride de rayonnement) *
  10. * *
  11. * 01/2006: mise au point d'une nouvelle méthode pour le rayonnement *
  12. * (méthode 2) *
  13. * le flux rayonné est mis sous la forme Phi=emis*sigma*(T**4-Trad**4) *
  14. * - Trad est évalué par l'opérateur RAYE (nouvelle option) *
  15. * en fonction des facteurs de forme, de l'emissivité et de la *
  16. * température à une itération donnée *
  17. * - la linéarisation se fait dans la procédure HRCAVNS *
  18. * - l'opérateur ECHI traite ensuite la condition aux limites *
  19. ***********************************************************************
  20.  
  21. COMPLET = faux ;
  22. SI ( COMPLET ) ;
  23. ni = 300 ; itt = 1. ; ome = 0.3 ;
  24. nn = 12 ;
  25. SINON ;
  26. ni = 30 ; itt = 1. ; ome = 0.3 ;
  27. nn = 6 ;
  28. FINSI ;
  29. GRAPH = faux ;
  30.  
  31.  
  32.  
  33. opti dime 2 elem qua8 ;
  34. opti isov suli ;
  35.  
  36. ******************
  37. * MAILLAGE *
  38. ******************
  39. l = 5.; h = 1.;
  40.  
  41. si ( l >eg h) ;
  42. ny = nn ; nx = ny * l / h ; nx = enti nx ;
  43. sinon ;
  44. nx = nn ; ny = nx * h / l ; ny = enti ny ;
  45. finsi ;
  46.  
  47.  
  48. ii = 3. ;ij = 1. / ii ;
  49.  
  50. d1 = l / nx / 2 * ij ; d2 = l / nx / 2 * ii;
  51. e1 = h / ny / 2 * ij ; e2 = h / ny / 2 * ii ;
  52.  
  53. discr = quaf ;
  54. discr = macro ;
  55.  
  56. epsm = 0.0001 ;
  57.  
  58. p1 = 0. 0. ; p2 = l 0. ; p3 = l h ; p4 = 0 h ;
  59. pb = (l/2) 0. ; ph = (l/2) h ; pd = l (h/2) ; pg = 0 (h/2) ;
  60. pm = (l/2) (h/2) ;
  61.  
  62.  
  63. bas = p1 d dini d1 dfin d2 pb d dini d2 dfin d1 p2 ;
  64. hau = p3 d dini d1 dfin d2 ph d dini d2 dfin d1 p4 ;
  65. gau = p4 d dini e1 dfin e2 pg d dini e2 dfin e1 p1 ;
  66. dro = p2 d dini e1 dfin e2 pd d dini e2 dfin e1 p3 ;
  67. cenv = pb d dini e1 dfin e2 pm d dini e2 dfin e1 ph ;
  68. miv = pg d dini d1 dfin d2 pm d dini d2 dfin d1 pd ;
  69.  
  70. cav = bas et dro et hau et gau ;
  71. dom = dall plan bas dro hau gau ;
  72.  
  73. *trac (cav et miv);
  74.  
  75. domf = chan dom quaf ;
  76. cavf = chan cav quaf ;
  77.  
  78. elim (domf et cavf et cenv et miv ) epsm ;
  79.  
  80. *trac dom ;
  81.  
  82. mdom = mode domf navier_stokes discr ;
  83. mcav = mode cavf navier_stokes discr ;
  84. qav = doma mcav maillage;
  85.  
  86. m_sray = mcav ;
  87.  
  88. ********************
  89. * PARAMETRES *
  90. ********************
  91. t0 = 1000. ; rapt = 0.7 ;
  92. tc = t0 ; tf = t0 * rapt ;
  93. re = 1.43 ;
  94. pr = 0.7 ;
  95. gr = 1.e5 ;
  96. ttm = tc + tf * 0.5 ;
  97. tm = 0.5 ;
  98. gb = 0. (-1.*gr /re/re) ;
  99. nu = 1./re ;
  100. pe = re * pr ;
  101. nrf = 3. ;
  102. sig = 5.67e-8 ;
  103. kko = pe * sig * t0 *t0 *t0 /nrf ;
  104. kla = sig * t0 *t0 *t0 /nrf ;
  105. *mess 'kko ' kko ' kla ' kla ;
  106.  
  107.  
  108. ****************************
  109. * DONNEES NUMERIQUES *
  110. ****************************
  111.  
  112. * création d'une table pour le rayonnement en cavité
  113.  
  114. mrm = mode qav thermique rayonnement cavite;
  115. marm = mate mrm emis 1. ;
  116.  
  117. tabra = table ;
  118.  
  119. tabra . ma_rai = qav ;
  120. tabra . mm_rai = mrm ;
  121.  
  122. tabra . mm_ns = m_sray;
  123.  
  124. tabra . mt_rai = marm ;
  125.  
  126. * calcul des facteurs de forme
  127. tabra . mf_rai = ffor mrm marm ;
  128.  
  129. * calcul de la matrice de rayonnement
  130. * cette étape n'est nécessaire que pour la méthode 1
  131. tabra . mr_rai = raye mrm (tabra . mf_rai) marm ;
  132.  
  133.  
  134.  
  135. *********************
  136. * CONVERGENCE *
  137. *********************
  138. debproc calcul ;
  139.  
  140. argu rvx*table ;
  141. rt = rvx . eqex ;
  142.  
  143. jt = dime (rt . inco . it ) ;
  144. jt = jt + 1 ;
  145.  
  146. unn = rt . inco . un ;
  147. un1 = rt . inco . u1 ;
  148.  
  149. tnn = rt . inco . tn ;
  150. tn1 = rt . inco . t1 ;
  151.  
  152. mdo = rt . 1calcul . domz ;
  153. mma = doma mdo maillage ;
  154.  
  155. unx = kcht mdo scal sommet (exco ux unn) ;
  156. uny = kcht mdo scal sommet (exco uy unn) ;
  157. u1x = kcht mdo scal sommet (exco ux un1) ;
  158. u1y = kcht mdo scal sommet (exco uy un1) ;
  159.  
  160. rt . inco . u1 = kcht mdo vect sommet unn ;
  161. rt . inco . t1 = kcht mdo scal sommet tnn ;
  162.  
  163. rt . inco . tet = kcht mdo scal sommet ((tnn - tf)/(tc - tf)) ;
  164.  
  165.  
  166. erux = kops unx - u1x ;
  167. elix = maxi erux abs ;
  168. elix = (log(elix + 1.e-10))/( log 10.) ;
  169. eruy = kops uny - u1y ;
  170. eliy = maxi eruy abs ;
  171. eliy = (log(eliy + 1.e-10))/( log 10.) ;
  172. ertt = kops tnn - tn1 ;
  173. elit = maxi ertt abs ;
  174. elit = (log(elit + 1.e-10))/( log 10.) ;
  175.  
  176. err = maxi (prog elix eliy) ;
  177. mess 'err' elix eliy elit jt ;
  178.  
  179. rt . inco . it = (rt . inco . it) et (prog jt) ;
  180. rt . inco . erx= (rt . inco . erx) et (prog elix);
  181. rt . inco . ery= (rt . inco . ery) et (prog eliy);
  182. rt . inco . ert= (rt . inco . ert) et (prog elit);
  183.  
  184. as af = kops matrik ;
  185. finproc as af ;
  186.  
  187. *-------------------------------------------------------------
  188.  
  189. 'DEBPROC' HRCAVNS;
  190. *---------------------------PROCEDURE RAYONNEMENT--------------------------*
  191.  
  192. ****************************************************************************
  193. * Procédure permettant de prendre en compte le rayonnement dans *
  194. * les operateurs eqex et exec (couplage convection-rayonnement) *
  195. * dans l'environnement du modèle Navier-Stokes *
  196. * Calcule le champ de température Trad et le coefficient d'echange Hrad *
  197. * correspondant à une linéarisation du rayonnement dans une cavité *
  198. * qui sont utlisés par l'opérateur ECHI s'appuyant sur la cavité *
  199. * *
  200. * (issu de la procédure RAY) *
  201. ****************************************************************************
  202.  
  203. **** création d'une table ****
  204.  
  205. 'ARGU' rvx * 'TABLE' ;
  206. rv = rvx . 'EQEX' ;
  207.  
  208. **** récuperation du champ par points des températures ****
  209.  
  210. line =rvx.'LISTINCO';
  211. noinco = 'EXTR' line 1 ;
  212. t = rv.'INCO'.noinco;
  213. tray = rvx . 'ARG1';
  214.  
  215. **** traitements ****
  216.  
  217. *maillage rayonnement*
  218. cavi = tray . 'MA_RAI';
  219.  
  220. *modele rayonnement*
  221. mrt = tray . 'MM_RAI' ;
  222.  
  223. *modele navier-stokes*
  224. m_ns = tray . 'MM_NS' ;
  225.  
  226. *emissivite
  227. emis = tray . 'MT_RAI' ;
  228.  
  229. *facteurs de forme
  230. ff = tray . 'MF_RAI' ;
  231.  
  232. mrn = rvx . 'DOMZ' ;
  233.  
  234. * mess ' ';
  235.  
  236. t_cavi = 'REDU' t cavi ;
  237. * mess ' t_cavi ' (mini t_cavi) (maxi t_cavi);
  238.  
  239. tc = 'EXCO' t_cavi scal 'T';
  240.  
  241. tre = 'CHAN' 'CHAM' tc mrt 'GRAVITE';
  242. * opti impi -1;
  243. * trad = 'RAYE' mrt ff emis tre 1e-8;
  244. trad = 'RAYE' mrt ff emis tre ;
  245. * opti impi 0;
  246. mess ' TRAD: ' (mini trad) (maxi trad);
  247.  
  248. *
  249. * transformation des CHAMELEM en scal centre pour C_fluide
  250. * 1/ temperature resultat de RAYE
  251. *
  252. trad = 'NOMC' trad 'SCAL';
  253. trad_n = 'CHAN' 'CHPO' mrt trad ;
  254. rv . inco . 'TRAD' = 'KCHT' m_ns scal sommet trad_n;
  255. trad = 'NOEL' m_ns (rv . inco . 'TRAD');
  256.  
  257. * 2/ émissivité
  258.  
  259. em_1 = 'NOMC' emis 'SCAL';
  260. em_n = 'CHAN' 'CHPO' mrt em_1 ;
  261.  
  262. em_s = 'NOEL' m_ns em_n;
  263.  
  264. te_cavi = 'NOEL' m_ns t_cavi;
  265. *
  266. * calcul du coefficient d'echange linéarisé
  267. *
  268. sig = 5.67e-8;
  269. te2 = te_cavi * te_cavi;
  270. trad2 = trad*trad;
  271.  
  272. hrad = em_s * sig * (te2 + trad2) * (te_cavi + trad);
  273. mess ' HRAD: ' (mini hrad) (maxi hrad);
  274.  
  275. rv . inco . 'HRAD' = 'KCHT' m_ns scal centre hrad;
  276.  
  277. as af = 'KOPS' 'MATRIK' ;
  278.  
  279. 'FINP' as af ;
  280.  
  281.  
  282. *---------------------FIN PROCEDURE RAYONNEMENT--------------------------*
  283.  
  284.  
  285. *-------------------------------------------------------------
  286. pip = elem (doma mdom centrep1) (lect 1) ;
  287.  
  288.  
  289. xx = coor 1 dom ;
  290.  
  291.  
  292.  
  293.  
  294.  
  295. **************************
  296. * SCHEMA IMPLICITE *
  297. **************************
  298.  
  299. * méthode 1
  300. * ---------
  301.  
  302. *-------------------------------------------------------------
  303. rw = eqex mdom niter ni omega ome
  304.  
  305. zone mdom oper calcul
  306. opti ef impl centrep1
  307. zone mdom oper kbbt 1. inco un pn
  308. zone mdom oper ns 1. 'UN' nu gb tet tm inco un
  309. zone mdom oper konv kko un 1. inco tn
  310. zone mdom oper lapn kla inco tn
  311. ;
  312. rw = eqex rw
  313. opti ef impl centree
  314. zone m_sray oper ray tabra inco tn
  315. ;
  316. rw = eqex rw
  317. clim un uimp cav 0.
  318. un vimp cav 0.
  319. tn timp gau tc
  320. tn timp dro tf
  321. pn timp pip 0.
  322. ;
  323.  
  324. rw . inco = table inco ;
  325. rw . inco . tabray = tabra ;
  326.  
  327. rw . inco . un = kcht mdom vect sommet (0. 0.) ;
  328. rw . inco . u1 = kcht mdom vect sommet (0. 0.) ;
  329. rw . inco . tn = kcht mdom scal sommet (tf - tc * xx / l + tc ) ;
  330. rw . inco . t1 = kcht mdom scal sommet (tf - tc * xx / 1 + tc ) ;
  331. rw . inco . pn = kcht mdom scal centrep1 0. ;
  332. rw . inco . dt = dtt ;
  333.  
  334. rw . inco . it = prog ;
  335. rw . inco . erx = prog ;
  336. rw . inco . ery = prog ;
  337. rw . inco . ert = prog ;
  338.  
  339. *-------------------------------------------------------------
  340.  
  341. * méthode 2
  342. * ---------
  343.  
  344. rw2= eqex mdom niter ni omega ome
  345.  
  346. zone mdom oper calcul
  347. opti ef impl centrep1
  348. zone mdom oper kbbt 1. inco un pn
  349. zone mdom oper ns 1. 'UN' nu gb tet tm inco un
  350. zone mdom oper konv kko un 1. inco tn
  351. zone mdom oper lapn kla inco tn
  352. ;
  353. rw2 = eqex rw2
  354. opti ef impl centree
  355. zone m_sray oper hrcavns tabra inco tn
  356. zone m_sray oper echi 'HRAD' 'TRAD' inco tn
  357. ;
  358. rw2 = eqex rw2
  359. clim un uimp cav 0.
  360. un vimp cav 0.
  361. tn timp gau tc
  362. tn timp dro tf
  363. pn timp pip 0.
  364. ;
  365.  
  366. rw2 . inco = table inco ;
  367. rw2 . inco . tabray = tabra ;
  368.  
  369. rw2 . inco . un = kcht mdom vect sommet (0. 0.) ;
  370. rw2 . inco . u1 = kcht mdom vect sommet (0. 0.) ;
  371. rw2 . inco . tn = kcht mdom scal sommet (tf - tc * xx / l + tc ) ;
  372. rw2 . inco . t1 = kcht mdom scal sommet (tf - tc * xx / 1 + tc ) ;
  373. rw2 . inco . pn = kcht mdom scal centrep1 0. ;
  374. rw2 . inco . dt = dtt ;
  375.  
  376. rw2 . inco . it = prog ;
  377. rw2 . inco . erx = prog ;
  378. rw2 . inco . ery = prog ;
  379. rw2 . inco . ert = prog ;
  380. *-------------------------------------------------------------
  381.  
  382. *******************
  383. * EXECUTION *
  384. *******************
  385. * titr 'impl macro' ;
  386.  
  387. exec rw ;
  388.  
  389. exec rw2;
  390.  
  391.  
  392. ***************
  393. * TESTS *
  394. ***************
  395.  
  396. * valeurs de référence pour le calcul non-complet
  397.  
  398. Lref = prog
  399. 7.0000E+2 7.1468E+2 7.2605E+2 7.3691E+2 7.4622E+2 7.5505E+2
  400. 7.6300E+2 7.7055E+2 7.7677E+2 7.8233E+2 7.8694E+2 7.9101E+2
  401. 7.9430E+2 7.9721E+2 7.9945E+2 8.0133E+2 8.0291E+2 8.0452E+2
  402. 8.0590E+2 8.0730E+2 8.0860E+2 8.0994E+2 8.1144E+2 8.1316E+2
  403. 8.1474E+2 8.1636E+2 8.1798E+2 8.1977E+2 8.2156E+2 8.2345E+2
  404. 8.2524E+2 8.2727E+2 8.2918E+2 8.3129E+2 8.3340E+2 8.3553E+2
  405. 8.3775E+2 8.4006E+2 8.4233E+2 8.4489E+2 ;
  406.  
  407. * méthode 1
  408. * ---------
  409. tt = rw . inco . tn ;
  410. Ttest = evol chpo hau tt;
  411. Ttest = extr Ttest ordo ;
  412. Ttest = extr Ttest (lect 1 pas 1 40);
  413. list Ttest;
  414. *opti sauv 'FORMAT' 'resu' ;
  415. *sauv 'FORMAT' Ttest;
  416.  
  417.  
  418. ERt=SOMM( abs ( Lref - Ttest)) ;
  419. ERm=maxi( abs ( Lref - Ttest)) ;
  420. mess ' méthode 1- Ecart sur profil de température ' ert erm;
  421.  
  422. * méthode 2
  423. * ---------
  424. tt = rw2. inco . tn ;
  425. Ttest = evol chpo hau tt;
  426. Ttest = extr Ttest ordo ;
  427. Ttest = extr Ttest (lect 1 pas 1 40);
  428. list Ttest;
  429. *opti sauv 'FORMAT' 'resu' ;
  430. *sauv 'FORMAT' Ttest;
  431.  
  432.  
  433. ERt2=SOMM( abs ( Lref - Ttest)) ;
  434. ERm2=maxi( abs ( Lref - Ttest)) ;
  435. mess ' méthode 2 - Ecart sur profil de température ' ert2 erm2;
  436.  
  437. SI ((ERm < 0.1) et (ERm2 < 0.7));
  438. ERREUR 0 ;
  439. SINON;
  440. ERREUR 5 ;
  441. FINSI;
  442.  
  443. ****************************************
  444. * PROCEDURE (TRACEE DES COURBES) *
  445. ****************************************
  446. si graph ;
  447.  
  448. schema= 'conv-rayo';
  449.  
  450. evx = evol bleu manu (rw . inco . it) (rw . inco . erx) ;
  451. evy = evol vert manu (rw . inco . it) (rw . inco . ery) ;
  452. evt = evol roug manu (rw . inco . it) (rw . inco . ery) ;
  453. titr 'Test de convergence, ' schema;
  454. dess (evx et evy et evt) ;
  455.  
  456. tt = rw . inco . tn ;
  457. tr = ( tt - tf ) / ( tc - tf ) ;
  458. titr 'Isothermes, ' schema;
  459. trac 50 tr dom cav ;
  460.  
  461.  
  462. vn = rw . inco . un ;
  463. titr 'Vitesses, ' schema;
  464. vv = vect vn (1.e2/gr)ux uy roug ;
  465. titr 'Vitesses, ' schema;
  466. trac vv dom cav ;
  467.  
  468. evx= EXCO vn ux;
  469. evy= EXCO vn uy;
  470. evux = evol rouge chpo evx cenv ;
  471. evuy = evol bleu chpo evy cenv ;
  472. titr 'Vitesses sur axe médian des y,'schema;
  473. dess (evux et evuy);
  474.  
  475. evux = evol rouge chpo evx miv ;
  476. evuy = evol bleu chpo evy miv ;
  477. titr 'Vitesses sur axe médian des x,'schema;
  478. dess (evux et evuy);
  479.  
  480. grat = kops tt grad mdom ;
  481. gratx = kcht mdom scal centre (exco ux grat) ;
  482. gratx = elno mdom gratx ;
  483. gratx = -1. * gratx ;
  484. gratxd = redu dro gratx ;
  485. evd = evol chpo gratxd scal dro ;
  486. gratxg = redu gau gratx ;
  487. evg = evol chpo gratxg scal gau ;
  488. titr ' gradient de température laterale, ' schema;
  489. dess (evd et evg) ;
  490.  
  491.  
  492. tb = redu bas tt ;
  493. tb = redu bas tr ;
  494. etb = evol vert chpo tb scal bas ;
  495. ho = inve hau ;
  496. th = redu ho tt ;
  497. th = redu ho tr ;
  498. eth = evol rouge chpo th scal ho ;
  499. titr 'Evolution température sur parois supérieure et
  500. inférieure,'schema;
  501. dess (etb et eth) ;
  502.  
  503.  
  504.  
  505. vv = rw . inco . un ;
  506. sw = kops vn rot mdom ;
  507.  
  508. rk = eqex mdom opti ef impl
  509. zone mdom oper lapn 1. inco psi
  510. zone mdom oper fimp sw inco psi
  511.  
  512. clim psi timp cav 0. ;
  513.  
  514. rk . inco . psi = kcht mdom scal sommet 0. ;
  515.  
  516. exec rk ;
  517.  
  518. trac (rk . inco . psi) dom cav 14 ;
  519.  
  520. finsi;
  521.  
  522.  
  523. **********************************
  524. * SAUVEGARDE DES RESULTATS *
  525. **********************************
  526. *opti sauv '/test4/ttmf4/rectangulaire/implrayo_proc.sav' ;
  527. *sauv ;
  528.  
  529.  
  530. fin ;
  531.  
  532.  
  533.  
  534.  
  535.  
  536.  
  537.  
  538.  
  539.  
  540.  
  541.  
  542.  
  543.  
  544.  
  545.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales