Télécharger jetplankei.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : jetplankei.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ************************************************************************
  5. * jetplankei.dgibi
  6. *
  7. * jet 2D Plan monophasique incompressible
  8. *
  9. * pour fiche de validation du K-epsilon
  10. *
  11. * Pierre Cornet , sept 97
  12. * jpm , mars 06 : adaptation pour K-epsilon implicite (kei)
  13. *
  14. ************************************************************************
  15. *
  16. COMPLET = FAUX;
  17. *--------------------------- Numérique ---------------------------------
  18. Si COMPLET;
  19. NBIT=400;
  20. DT = 0.002 ;
  21. DISCR = MACRO;
  22. KPRESS= CENTREP1;
  23. sinon;
  24. NBIT=80;
  25. DT = 0.01 ;
  26. DISCR = LINE;
  27. KPRESS= CENTRE;
  28. finsi;
  29.  
  30. graph = faux;
  31. *graph = vrai;
  32. *opti trace 'PSC' ;
  33. *----------------------- fin Numérique ---------------------------------
  34.  
  35. *--------------------------- procédure ---------------------------------
  36. *$$$$ OUTFLOW
  37. 'DEBPROC' OUTFLOW ;
  38. ARGU RX*TABLE ;
  39. ************************************************************************
  40. * ZONE $Mt OPER 'OUTFLOW' $out Ro Un Muf INCO 'UN'
  41. *
  42. * UN
  43. * MUF viscosité dynamique effective
  44. * (tourbillonnaire+moléculaire)
  45. ************************************************************************
  46. rv=rx.'EQEX' ;
  47. iarg=rx.'IARG' ;
  48. *mess 'iarg='iarg;
  49. $mt=rx.'DOMZ' ;
  50.  
  51. * Lecture du 1er Argument objet mmodel
  52. Si(ega ('TYPE' rx.'ARG1') 'MMODEL ');
  53. $mfront=rx.'ARG1';
  54. Dg=doma $mfront 'XXDIAGSI' ;
  55. Sinon ;
  56. erreur 2;
  57. quitter outflow;
  58. Finsi ;
  59.  
  60. * Lecture du 2ème Argument la densité
  61. Si(ega ('TYPE' rx.'ARG2') 'MOT ');
  62. Ro = rv.inco.(rx.'ARG2');
  63. Sinon ;
  64. Ro = rx.'ARG2';
  65. Finsi ;
  66.  
  67. * Lecture du 3ème Argument la vitesse
  68. Si(ega ('TYPE' rx.'ARG3') 'MOT ');
  69. un = rv.inco.(rx.'ARG3');
  70. Sinon ;
  71. un = rx.'ARG3';
  72. Finsi ;
  73.  
  74. * Lecture du 4ème Argument la viscosité
  75. Si(ega ('TYPE' rx.'ARG4') 'MOT ');
  76. Muf = rv.inco.(rx.'ARG4');
  77. Sinon ;
  78. Muf = rx.'ARG4';
  79. Finsi ;
  80.  
  81. **************** INITIALISATIONs ***********************
  82. Si (Exist rx 'rxtm');
  83. rxtm=rx.'rxtm';
  84. lmx=rxtm.'lmx';
  85. lm0=rxtm.'lm0';
  86. lm=rxtm.'lm';
  87. nc=rxtm.'nc';
  88. MSV = chai rxtm.'msv' ;
  89. MI1 = rxtm.'MI1' ;
  90. KPRES=rxtm.'KPRES';
  91. Sinon;
  92.  
  93. rxtm=table 'KIZX';
  94. rx.'rxtm'=rxtm;
  95. Si (NON (EGA (dime rx.'LISTINCO') 1));
  96. mess ' Erreur dans la procédure Outflow !';
  97. mess ' Il doit y avoir une inconnue !';
  98. quitter outflow ;
  99. Finsi;
  100. MI1=EXTR rx.'LISTINCO' 1 ;
  101. lmx=(extr (rv.inco.MI1) 'COMP');
  102. lm0=(extr un 'COMP');
  103. nc=dime lm0;
  104.  
  105. Si(nc > 1);
  106. *mess 'Outflow Cas vectoriel ' MI1;
  107. MSV=chai 'VECT';
  108. lm=mots (chai 1 MI1) (chai 2 MI1);
  109. Si (EGA (vale dime) 3);
  110. lm=lm et (mots (chai 3 MI1));
  111. Finsi;
  112. kpr=rx.'KOPT'.'KPOIN';
  113. Si(EGA KPR 2);KPRES='CENTRE';Finsi;
  114. Si(EGA KPR 4);KPRES='CENTREP1';Finsi;
  115. Si(EGA KPR 5);KPRES='MSOMMET';Finsi;
  116. Sinon;
  117. *mess 'Outflow Cas scalaire ' MI1;
  118. MSV=chai 'SCAL';
  119. lm=mots MI1;
  120. Finsi;
  121. rxtm.'lmx'=lmx;
  122. rxtm.'lm0'=lm0;
  123. rxtm.'lm'=lm;
  124. rxtm.'nc'=nc;
  125. rxtm.'msv'=chai MSV;
  126. rxtm.'MI1'=MI1;
  127. rxtm.'KPRES'=KPRES;
  128.  
  129. rxtm.'LISTINCO'=MOTS MI1 ;
  130. rxtm.'DOMZ'=$mfront;
  131. rxtm.'KOPT'=rx.'KOPT';
  132. rxtm.'IARG'=1;
  133. rxtm.'EQEX'=rx.'EQEX';
  134. rxtm.'NOMZONE'=' ';
  135. rxtm.'TDOMZ'=0;
  136. rxtm.'NOMOPER'=mot 'MDIA_T';
  137. Finsi;
  138. **************** Fin Initialisations *******************
  139.  
  140. mfront=doma $mfront maillage ;
  141. nj=doma $mt 'NORMALEV';
  142. *unj= vect nj 0.1 ux uy jaune ;
  143. *trace unj mt;
  144. njf = redu nj mfront ;
  145. *unj= vect njf 0.1 ux uy rouge ;
  146. *trace unj mt TITR ' FRONTIERE';
  147. xn = rv.inco.MI1;
  148. ujf = redu un mfront ;
  149. us=psca ujf lm0 njf lm0;
  150. ius=masq us 'INFERIEUR' 0.;
  151. rv.inco.'us'=us;
  152. rv.inco.'ius'=ius;
  153. Si(Ega nc 1);
  154. us=2.*(redu xn mfront);
  155. Finsi;
  156. rxtm.'ARG1'=kcht $mfront scal sommet ((-1.)*Ro*us*ius);
  157.  
  158. st1 mat1 = MDIA rxtm ;
  159. *st1 mat1 = KOPS 'MATRIK' ;
  160.  
  161. Dgi=Dg*ius;
  162. Si(nc > 1);
  163. puj=0.;
  164. Si(EXIST (rv.inco) 'PRESSION');
  165. pn=redu (elno $mt rv.inco.'PRESSION' KPRES) mfront;
  166. puj=Ro*pn*us;
  167. Finsi;
  168. gru = redu (Muf*(kops 'GRADS' (exco un 'UX') $mt)) mfront;
  169. grv = redu (Muf*(kops 'GRADS' (exco un 'UY') $mt)) mfront;
  170. mgunj = (nomc (Dgi*((psca gru lm0 njf lm0)-puj)) (extr lm 1))
  171. + (nomc (Dgi*((psca grv lm0 njf lm0)-puj)) (extr lm 2)) ;
  172. Si(Ega (vale dime) 3);
  173. mgunj = mgunj + (nomc (Dgi*((psca grv lm0 njf lm0)-puj)) (extr lm 3));
  174. Finsi;
  175. Sinon;
  176. grt = redu (Muf*(kops 'GRADS' xn $mt)) mfront ;
  177. mgunj = (nomc (Dgi*(psca grt lm0 njf lm0)) (extr lm 1)) ;
  178. Finsi;
  179. St1 = kcht $mfront MSV sommet comp lm ((-1.)*mgunj);
  180.  
  181. RESPRO St1 mat1 ;
  182. FINPROC ;
  183. *----------------------- fin procédure ---------------------------------
  184.  
  185.  
  186. *--------------------------- maillage ----------------------------------
  187.  
  188. TITRE 'JET' ;
  189. OPTI DIME 2 ELEM QUA8 ;
  190. *
  191. DJ = 2.e-2 ; RJ = DJ/2. ; RM = 150.*RJ ;
  192. DJ = 2.e-2 ; RJ = DJ/2. ; RM = 50.*RJ ;
  193.  
  194. * points :
  195. P00=0. 0.; PJ0=RJ 0.; PR0=RM 0.; PJ5=RJ (50.*DJ);
  196. P02=0. (20.*DJ); P03=0. (30.*DJ); P04=0. (40.*DJ); P05=0. (50.*DJ);
  197. PR2=RM (20.*DJ); PR3=RM (30.*DJ); PR4=RM (40.*DJ); PR5=RM (50.*DJ);
  198.  
  199. * segments verticaux :
  200. A02 = DROI -20 P00 P02 dini (0.5*DJ) dfin (1.25*DJ) ;
  201. A23 = DROI -6 P02 P03 dini (0.9*DJ) dfin (1.25*DJ) ;
  202. A34 = DROI 4 P03 P04 ;
  203. A45 = DROI -3 P04 P05 dini (1.25*DJ) dfin (2.*DJ) ;
  204.  
  205. B02 = DROI -20 PR0 PR2 dini (0.5*DJ) dfin (1.25*DJ) ;
  206. B23 = DROI -6 PR2 PR3 dini (0.9*DJ) dfin (1.25*DJ) ;
  207. B34 = DROI 4 PR3 PR4 ;
  208. B45 = DROI -3 PR4 PR5 dini (1.25*DJ) dfin (2.*DJ) ;
  209.  
  210. AXE = A02 ET A23 ET A34 ET A45 ;
  211. BORD= B02 ET B23 ET B34 ET B45 ;
  212.  
  213. * segments horizontaux :
  214. JET = DROI 2 P00 PJ0 ;
  215. BAS2 = DROI PJ0 PR0 dini (RJ/2.) dfin (10.*RJ) ; BAS = JET ET BAS2 ;
  216.  
  217. HAU1 = DROI 2 P05 PJ5 ;
  218. HAU2 = DROI PJ5 PR5 dini (RJ/2.) dfin (10.*RJ) ; HAUT=HAU1 ET HAU2 ;
  219.  
  220. * domaine total :
  221. MT = DALL BAS BORD (inve HAUT) (inve AXE) 'PLAN' ;
  222. CMT = cont MT ;
  223.  
  224. Mmt = chan QUAF MT ;
  225. Mjet = chan QUAF jet ;
  226. Mbord= chan QUAF bord;
  227. Mbas2= chan QUAF bas2;
  228. Maxe = chan QUAF axe ;
  229. Mhaut= chan QUAF haut;
  230.  
  231. elim (Mmt et Mjet et Mbord et Mbas2 et Maxe et Mhaut) 1.e-5;
  232.  
  233. $MT = mode Mmt 'NAVIER_STOKES' DISCR;
  234. $JET = mode Mjet 'NAVIER_STOKES' DISCR;
  235. $BORD = mode Mbord 'NAVIER_STOKES' DISCR;
  236. $BAS2 = mode Mbas2 'NAVIER_STOKES' DISCR;
  237. $AXE = mode Maxe 'NAVIER_STOKES' DISCR;
  238. $HAUT = mode Mhaut 'NAVIER_STOKES' DISCR;
  239.  
  240. Ml20d = Mhaut moins (P05 moins P02) coul verte;
  241. Ml30d = Mhaut moins (P05 moins P03) coul verte;
  242. Ml40d = Mhaut moins (P05 moins P04) coul verte;
  243.  
  244. elim (Mmt et Ml20d et Ml30d et Ml40d) 1.e-5;
  245.  
  246. MT = doma $MT maillage ;
  247. axe = doma $axe maillage ;
  248. jet = doma $jet maillage ;
  249. bord= doma $bord maillage ;
  250. bas2= doma $bas2 maillage ;
  251. haut= doma $haut maillage ;
  252.  
  253. Si Graph;trace MT;Finsi;
  254. *------------------------ donnees physiques ----------------------------
  255. *
  256. NUF = 1.5E-5 ;
  257. REJ = 2.e4 ;
  258. UJ = REJ*NUF/DJ ; mess 'vitesse d injection (m/s) =' UJ ;
  259.  
  260. * KJ = 0.05*UJ*UJ ; mess 'KJ =' KJ ;
  261. * EJ = 0.02*(UJ**3.)/DJ ; mess 'EJ =' EJ ;
  262. KJ = 1.E-3 ;
  263. EJ = 6.E-3 ;
  264. NUTj = 0.09*KJ*KJ/EJ ; mess 'NUTJ =' nutj ;
  265. KA = 1.E-7 ;
  266. EA = 1.E-5 ;
  267. L0 = 25.*DJ ;
  268.  
  269. *-------------------------- equations ----------------------------------
  270.  
  271. RV = EQEX $MT 'DUMP' 'ITMA' NBIT
  272. 'OPTI' 'EF' 'SUPG' 'IMPL'
  273. * 'ZONE' $MT OPER OUTFLOW $haut 1. 'UN' 'MUF' 'INCO' 'KN'
  274. * 'ZONE' $MT OPER OUTFLOW $haut 1. 'UN' 'MUF' 'INCO' 'EN'
  275. 'ZONE' $MT OPER OUTFLOW $haut 1. 'UN' 'MUF' 'INCO' 'UN'
  276. 'ZONE' $MT 'OPER' 'KEPSILON' 1. 'UN' NUF 'DT' 'INCO' 'KN' 'EN'
  277. 'ZONE' $MT 'OPER' 'NS' 1. 'UN' 'MUF' 'INCO' 'UN'
  278. 'OPTI' 'EFM1' 'CENTREE'
  279. 'ZONE' $MT 'OPER' DFDT 1. 'UN' 'DT' 'INCO' 'UN';
  280.  
  281. RV = EQEX RV
  282. 'CLIM' 'UN' UIMP JET 0. 'UN' VIMP JET UJ
  283. 'UN' VIMP BAS2 0. 'UN' UIMP AXE 0.
  284. 'UN' VIMP BORD 0.
  285. 'KN' TIMP JET KJ 'EN' TIMP JET EJ
  286. 'KN' TIMP BORD 0. 'EN' TIMP BORD EA ;
  287.  
  288. *RV.'ALGO_KEPSILON'= MOTS 'RNG';
  289.  
  290. RVP = EQEX 'OPTI' 'EF' KPRESS
  291. 'ZONE' $MT 'OPER' 'KBBT' -1. 'INCO' 'UN' 'PRES' ;
  292.  
  293. RV.'PROJ'= RVP ;
  294.  
  295. *------------------------ initialisations ------------------------------
  296.  
  297. RV.INCO = TABLE 'INCO' ;
  298. RV.'INCO'.'UN' = KCHT $MT VECT SOMMET (1.E-7 1.E-7) ;
  299. RV.'INCO'.'PRES' = 'KCHT' $MT 'SCAL' KPRESS 0.;
  300. RV.'INCO'.'KN' = KCHT $MT SCAL SOMMET 1.E-7 ;
  301. RV.'INCO'.'EN' = KCHT $MT SCAL SOMMET 1.E-5 ;
  302. RV.'INCO'.'MUF' = KCHT $MT SCAL SOMMET 1.E-6 ;
  303. RV.'INCO'.'DT' = DT ;
  304.  
  305. *------------------------ historiques ----------------------------------
  306.  
  307. P11 = MT POIN 'PROC' ((25.*RJ) (10.*DJ));
  308. P12 = MT POIN 'PROC' ((25.*RJ) (20.*DJ));
  309. P14 = MT POIN 'PROC' ((25.*RJ) (40.*DJ));
  310. P15 = MT POIN 'PROC' ((25.*RJ) (50.*DJ));
  311. LH = P02 et P03 et P04 et P05 et P11 et P12 et P14 et P15 ;
  312. HIS = KHIS 'UN' 1 LH 'UN' 2 LH 'KN' LH 'EN' LH ;
  313. RV.'HIST' = HIS ;
  314.  
  315. *------------------------ resolution -----------------------------------
  316.  
  317. lh= (POIN MT PROC (0.05 0.01)) et (POIN MT PROC (0.05 0.05))
  318. et (POIN MT PROC (0.05 0.1)) et (POIN MT PROC (0.05 0.15))
  319. et (POIN MT PROC (0.05 0.2)) et (POIN MT PROC (0.05 0.7))
  320. et (POIN MT PROC (0.05 0.3)) et (POIN MT PROC (0.05 0.6))
  321. et (POIN MT PROC (0.05 0.4)) et (POIN MT PROC (0.05 0.5));
  322.  
  323. his=khis 'UN' 1 lh
  324. 'UN' 2 lh
  325. 'KN' lh
  326. 'EN' lh;
  327. his.'KFIH'=1;
  328. rv.'HIST'=his;
  329.  
  330. EXEC rv;
  331.  
  332. UN = (RV.'INCO'.'UN');
  333. uaxe = ((exco UN 'UY') redu axe)*(1./UJ);
  334. evuaxe= evol chpo uaxe axe ;
  335.  
  336. *------------------------ post-traitement ------------------------------
  337. Si Graph;
  338. trace ((cont mt) et Ml20d et Ml30d et Ml40d);
  339.  
  340. dessin his.'TABD' his.'KN' TITR ' historiques: k';
  341. dessin his.'TABD' his.'EN' TITR ' historiques: epsilon';
  342. dessin his.'TABD' his.'1UN' TITR ' historiques: ux';
  343. dessin his.'TABD' his.'2UN' TITR ' historiques: uy';
  344.  
  345. UNV = VECT UN 5.E-3 UX UY ;
  346. TRAC UNV CMT TITRE 'VITESSES AIR' ;
  347.  
  348. dess evuaxe TITR 'Vitesse sur l axe';
  349.  
  350. z= (coor 2 axe) + 1.e-3;
  351. vaxe= 5.8 * Dj * (inve z);
  352. evax = evol chpo vaxe axe ;
  353. dess (evax et evuaxe) ybor 0.1 1.2 TITR ' Vitesses sur l axe';
  354.  
  355. u20d = (exco un 'UY') redu Ml20d;
  356. um20d = maxi u20d;
  357. u20d = evol chpo (u20d*(1./um20d)) ml20d ;
  358.  
  359. u30d = (exco un 'UY') redu Ml30d;
  360. um30d = maxi u30d;
  361. u30d = evol chpo (u30d*(1./um30d)) ml30d ;
  362.  
  363. u40d = (exco un 'UY') redu Ml40d;
  364. um40d = maxi u40d;
  365. u40d = evol chpo (u40d*(1./um40d)) ml40d ;
  366.  
  367. z20d = maxi (coor 2 ml20d);
  368. z30d = maxi (coor 2 ml30d);
  369. z40d = maxi (coor 2 ml40d);
  370. r20d = (extr u20d 'ABSC')*(1./z20d);
  371. r30d = (extr u30d 'ABSC')*(1./z30d);
  372. r40d = (extr u40d 'ABSC')*(1./z40d);
  373.  
  374. u20d = evol manu r20d (extr u20d 'ORDO');
  375. u30d = evol manu r30d (extr u30d 'ORDO');
  376. u40d = evol manu r40d (extr u40d 'ORDO');
  377.  
  378. dess (u20d et u30d et u40d) XBOR 0. 0.4
  379. TITR ' Profil radiaux de vitesse';
  380.  
  381.  
  382. nutsnu = rv.inco.'MUF' * (1./NUF) ;
  383. trace nutsnu mt (cont mt) TITR 'Muf / Nu' ;
  384.  
  385. trace rv.inco.'KN' mt (cont mt) TITR ' Kn ';
  386. trace rv.inco.'EN' mt (cont mt) TITR ' En ';
  387. Finsi;
  388. *-------------------- fin post-traitement ------------------------------
  389.  
  390. Si (NON COMPLET);
  391. uaxref=prog
  392. 1.0000 1.0250 0.97229 1.0009 0.96387
  393. 0.99073 0.97449 1.0046 1.0016 1.0228
  394. 1.0149 1.0024 0.96638 0.93515 0.89487
  395. 0.85870 0.81994 0.78738 0.75198 0.72043
  396. 0.69096 0.66619 0.64104 0.61681 0.59295
  397. 0.57226 0.55194 0.52699 0.50559 0.48760
  398. 0.47019 0.44757 0.42364 0.40301 ;
  399. zaxe = extr evuaxe 'ABSC';
  400. list (extr evuaxe 'ORDO');
  401. evuaxer=evol manu zaxe uaxref ;
  402. m=(INTG evuaxer 'ABSO') SOMM ;
  403. delt2= SOMM (INTG (evuaxer - evuaxe) 'ABSO')/m;
  404. mess ' deltaL2 : m=' m ' delt2=' delt2;
  405. Si (delt2 > 8.e-3);ERREUR 5;Finsi;
  406. FINSI;
  407.  
  408.  
  409.  
  410.  
  411.  
  412. FIN ;
  413.  
  414.  
  415.  
  416.  
  417.  
  418.  
  419.  
  420.  
  421.  
  422.  
  423.  
  424.  
  425.  

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