Télécharger jetkei.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : jetkei.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ************************************************************************
  5. * jetkei.dgibi
  6. *
  7. * jet 2D axi 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. graph = faux;
  17. *graph = vrai;
  18. *opti trace 'PSC' ;
  19. COMPLET = FAUX;
  20.  
  21. *--------------------------- Numérique ---------------------------------
  22. Si COMPLET;
  23.  
  24. KMAIL=1;
  25. zr=1./3.;
  26. zr=1.;
  27. NBIT=2000;
  28. NBIT=10000;
  29. NBIT=3000;
  30. NBIT=10000;
  31. CFL=5.;
  32. CFL= 1;
  33. NBIT=400;
  34. NBIT=3000;
  35. *NBIT=800;
  36. *NBIT=3 ;
  37. CFL=5.;
  38. DISCR = MACRO;
  39. KPRESS= CENTREP1;
  40. DISCR = LINE;
  41. KPRESS= CENTRE;
  42. sinon;
  43.  
  44. KMAIL=1;
  45. zr=2.;
  46. NBIT=80;
  47. CFL=15.;
  48. DISCR = LINE;
  49. KPRESS= CENTRE;
  50. finsi;
  51.  
  52. *--------------------------- procédure ---------------------------------
  53. *$$$$ OUTFLOW
  54. 'DEBPROC' OUTFLOW ;
  55. ARGU RX*TABLE ;
  56. ************************************************************************
  57. * ZONE $Mt OPER 'OUTFLOW' $out Ro Un Muf INCO 'UN'
  58. *
  59. * UN
  60. * MUF viscosité dynamique effective
  61. * (tourbillonnaire+moléculaire)
  62. ************************************************************************
  63. rv=rx.'EQEX' ;
  64. iarg=rx.'IARG' ;
  65. *mess 'iarg='iarg;
  66. $mt=rx.'DOMZ' ;
  67.  
  68. * Lecture du 1er Argument objet mmodel
  69. Si(ega ('TYPE' rx.'ARG1') 'MMODEL ');
  70. $mfront=rx.'ARG1';
  71. Dg=doma $mfront 'XXDIAGSI' ;
  72. Sinon ;
  73. erreur 2;
  74. quitter outflow;
  75. Finsi ;
  76.  
  77. * Lecture du 2ème Argument la densité
  78. Si(ega ('TYPE' rx.'ARG2') 'MOT ');
  79. Ro = rv.inco.(rx.'ARG2');
  80. Sinon ;
  81. Ro = rx.'ARG2';
  82. Finsi ;
  83.  
  84. * Lecture du 3ème Argument la vitesse
  85. Si(ega ('TYPE' rx.'ARG3') 'MOT ');
  86. un = rv.inco.(rx.'ARG3');
  87. Sinon ;
  88. un = rx.'ARG3';
  89. Finsi ;
  90.  
  91. * Lecture du 4ème Argument la viscosité
  92. Si(ega ('TYPE' rx.'ARG4') 'MOT ');
  93. Muf = rv.inco.(rx.'ARG4');
  94. Sinon ;
  95. Muf = rx.'ARG4';
  96. Finsi ;
  97.  
  98. **************** INITIALISATIONs ***********************
  99. Si (Exist rx 'rxtm');
  100. rxtm=rx.'rxtm';
  101. lmx=rxtm.'lmx';
  102. lm0=rxtm.'lm0';
  103. lm=rxtm.'lm';
  104. nc=rxtm.'nc';
  105. MSV = chai rxtm.'msv' ;
  106. MI1 = rxtm.'MI1' ;
  107. KPRES=rxtm.'KPRES';
  108. Sinon;
  109.  
  110. rxtm=table 'KIZX';
  111. rx.'rxtm'=rxtm;
  112. Si (NON (EGA (dime rx.'LISTINCO') 1));
  113. mess ' Erreur dans la procédure Outflow !';
  114. mess ' Il doit y avoir une inconnue !';
  115. quitter outflow ;
  116. Finsi;
  117. MI1=EXTR rx.'LISTINCO' 1 ;
  118. lmx=(extr (rv.inco.MI1) 'COMP');
  119. lm0=(extr un 'COMP');
  120. nc=dime lm0;
  121.  
  122. Si(nc > 1);
  123. *mess 'Outflow Cas vectoriel ' MI1;
  124. MSV=chai 'VECT';
  125. lm=mots (chai 1 MI1) (chai 2 MI1);
  126. Si (EGA (vale dime) 3);
  127. lm=lm et (mots (chai 3 MI1));
  128. Finsi;
  129. kpr=rx.'KOPT'.'KPOIN';
  130. Si(EGA KPR 2);KPRES='CENTRE';Finsi;
  131. Si(EGA KPR 4);KPRES='CENTREP1';Finsi;
  132. Si(EGA KPR 5);KPRES='MSOMMET';Finsi;
  133. Sinon;
  134. *mess 'Outflow Cas scalaire ' MI1;
  135. MSV=chai 'SCAL';
  136. lm=mots MI1;
  137. Finsi;
  138. rxtm.'lmx'=lmx;
  139. rxtm.'lm0'=lm0;
  140. rxtm.'lm'=lm;
  141. rxtm.'nc'=nc;
  142. rxtm.'msv'=chai MSV;
  143. rxtm.'MI1'=MI1;
  144. rxtm.'KPRES'=KPRES;
  145.  
  146. rxtm.'LISTINCO'=MOTS MI1 ;
  147. rxtm.'DOMZ'=$mfront;
  148. rxtm.'KOPT'=rx.'KOPT';
  149. rxtm.'IARG'=1;
  150. rxtm.'EQEX'=rx.'EQEX';
  151. rxtm.'NOMZONE'=' ';
  152. rxtm.'TDOMZ'=0;
  153. rxtm.'NOMOPER'=mot 'MDIA_T';
  154. Finsi;
  155. **************** Fin Initialisations *******************
  156.  
  157. mfront=doma $mfront maillage ;
  158. nj=doma $mt 'NORMALEV';
  159. *unj= vect nj 0.1 ux uy jaune ;
  160. *trace unj mt;
  161. njf = redu nj mfront ;
  162. *unj= vect njf 0.1 ux uy rouge ;
  163. *trace unj mt TITR ' FRONTIERE';
  164. xn = rv.inco.MI1;
  165. ujf = redu un mfront ;
  166. us=psca ujf lm0 njf lm0;
  167. ius=masq us 'INFERIEUR' 0.;
  168. rv.inco.'us'=us;
  169. rv.inco.'ius'=ius;
  170. Si(Ega nc 1);
  171. us=2.*(redu xn mfront);
  172. Finsi;
  173. rxtm.'ARG1'=kcht $mfront scal sommet ((-1.)*Ro*us*ius);
  174.  
  175. st1 mat1 = MDIA rxtm ;
  176. *st1 mat1 = KOPS 'MATRIK' ;
  177.  
  178. Dgi=Dg*ius;
  179. Si(nc > 1);
  180. puj=0.;
  181. Si(EXIST (rv.inco) 'PRESSION');
  182. pn=redu (elno $mt rv.inco.'PRESSION' KPRES) mfront;
  183. puj=Ro*pn*us;
  184. Finsi;
  185. gru = redu (Muf*(kops 'GRADS' (exco un 'UX') $mt)) mfront;
  186. grv = redu (Muf*(kops 'GRADS' (exco un 'UY') $mt)) mfront;
  187. mgunj = (nomc (Dgi*((psca gru lm0 njf lm0)-puj)) (extr lm 1))
  188. + (nomc (Dgi*((psca grv lm0 njf lm0)-puj)) (extr lm 2)) ;
  189. Si(Ega (vale dime) 3);
  190. mgunj = mgunj + (nomc (Dgi*((psca grv lm0 njf lm0)-puj)) (extr lm 3));
  191. Finsi;
  192. Sinon;
  193. grt = redu (Muf*(kops 'GRADS' xn $mt)) mfront ;
  194. mgunj = (nomc (Dgi*(psca grt lm0 njf lm0)) (extr lm 1)) ;
  195. Finsi;
  196. St1 = kcht $mfront MSV sommet comp lm ((-1.)*mgunj);
  197.  
  198. RESPRO St1 mat1 ;
  199. FINPROC ;
  200. *----------------------- fin procédure ---------------------------------
  201.  
  202.  
  203.  
  204.  
  205. *--------------------------- maillage ----------------------------------
  206.  
  207. TITRE 'JET' ;
  208. OPTI MODE AXIS ;
  209. OPTI DIME 2 ELEM QUA8 ;
  210.  
  211. DJ = 2.e-2 ; RJ = DJ/2. ; RM = 50.*RJ ;
  212. epsi=1.e-7;
  213.  
  214. * points :
  215. DJJ=1.*DJ;
  216. P00=0. 0.; PJ0=RJ 0.; PR0=RM 0.; PJ5=RJ (50.*DJJ);
  217. P02=0. (20.*DJJ); P03=0. (30.*DJJ); P04=0. (40.*DJJ); P05=0. (50.*DJJ);
  218. PR2=RM (20.*DJJ); PR3=RM (30.*DJJ); PR4=RM (40.*DJJ); PR5=RM (50.*DJJ);
  219.  
  220. * segments verticaux :
  221. Si (EGA KMAIL 1);
  222. A02 = DROI P00 P02 dini (0.5*zr*DJ) dfin (1.*zr*DJ) ;
  223. A23 = DROI P02 P03 dini (1.*zr*DJ) dfin (1.2*zr*DJ) ;
  224. A34 = DROI P03 P04 dini (1.2*zr*DJ) dfin (1.3*zr*DJ) ;
  225. A45 = DROI P04 P05 dini (1.3*zr*DJ) dfin (1.4*zr*DJ) ;
  226. Finsi;
  227.  
  228. Si (EGA KMAIL 2);
  229. A02 = DROI -40 P00 P02 dini (0.25*DJ) dfin (0.6*DJ) ;
  230. A23 = DROI -12 P02 P03 dini (0.45*DJ) dfin (0.6*DJ) ;
  231. A34 = DROI 8 P03 P04 ;
  232. A45 = DROI -6 P04 P05 dini (0.6*DJ) dfin (1.*DJ) ;
  233. Finsi;
  234.  
  235. Si (EGA KMAIL 3);
  236. A02 = DROI -80 P00 P02 dini (0.1*DJ) dfin (0.6*DJ) ;
  237. A23 = DROI -12 P02 P03 dini (0.45*DJ) dfin (0.6*DJ) ;
  238. A34 = DROI 8 P03 P04 ;
  239. A45 = DROI -6 P04 P05 dini (0.6*DJ) dfin (1.*DJ) ;
  240. Finsi;
  241.  
  242. AXE = A02 ET A23 ET A34 ET A45 ;
  243. BORD= AXE plus PR0 ;
  244.  
  245. * segments horizontaux :
  246. Si (EGA KMAIL 1);
  247. JET = DROI 2 P00 PJ0 ;
  248. BAS2 = DROI PJ0 PR0 dini (zr*RJ/2.) dfin (5.*zr*RJ) ;
  249. BAS = JET ET BAS2 ;
  250. Finsi;
  251.  
  252. Si (EGA KMAIL 2);
  253. JET = DROI 5 P00 PJ0 ;
  254. BAS2 = DROI PJ0 PR0 dini (RJ/5.) dfin (10.*RJ) ;
  255. BAS = JET ET BAS2 ;
  256. Finsi;
  257.  
  258. Si (EGA KMAIL 3);
  259. JET = DROI 10 P00 PJ0 ;
  260. BAS2 = DROI PJ0 PR0 dini (RJ/10.) dfin (10.*RJ) ;
  261. BAS = JET ET BAS2 ;
  262. Finsi;
  263.  
  264. HAUT = BAS plus P05;
  265. elim (haut et axe et bord et bas) epsi;
  266.  
  267. * domaine total :
  268. MT = DALL BAS BORD (inve HAUT) (inve AXE) 'PLAN' ;
  269. CMT = cont MT ;
  270.  
  271. Mmt = chan QUAF MT ;
  272. Mjet = chan QUAF jet ;
  273. Mbord= chan QUAF bord;
  274. Mbas2= chan QUAF bas2;
  275. Maxe = chan QUAF axe ;
  276. Mhaut= chan QUAF haut;
  277. Mout = chan QUAF (bord et haut);
  278.  
  279. elim (Mmt et Mjet et Mbord et Mbas2 et Maxe et Mhaut) epsi;
  280. elim (Mmt et Mout) epsi;
  281.  
  282. $MT = mode Mmt 'NAVIER_STOKES' DISCR;
  283. $JET = mode Mjet 'NAVIER_STOKES' DISCR;
  284. $BORD = mode Mbord 'NAVIER_STOKES' DISCR;
  285. $BAS2 = mode Mbas2 'NAVIER_STOKES' DISCR;
  286. $AXE = mode Maxe 'NAVIER_STOKES' DISCR;
  287. $HAUT = mode Mhaut 'NAVIER_STOKES' DISCR;
  288. $out = mode Mout 'NAVIER_STOKES' DISCR;
  289.  
  290. Ml20d = Mhaut moins (P05 moins P02) coul verte;
  291. Ml30d = Mhaut moins (P05 moins P03) coul verte;
  292. Ml40d = Mhaut moins (P05 moins P04) coul verte;
  293.  
  294. elim (Mmt et Ml20d et Ml30d et Ml40d) epsi;
  295.  
  296. MT = doma $MT maillage ;
  297. axe = doma $axe maillage ;
  298. jet = doma $jet maillage ;
  299. bas2= doma $bas2 maillage ;
  300. bord= doma $bord maillage ;
  301. out = doma $out maillage ;
  302. haut= doma $haut maillage;
  303.  
  304. doma $mt 'IMPR';
  305.  
  306. Si Graph;
  307. * On prend soin de vérifier que l'allongement n'est pas supérieur à 20
  308. diamin=doma $mt 'DIAMIN';
  309. diamax=doma $mt 'DIAMAX';
  310. ksi=diamax*(inve diamin);
  311. ksi=elno $mt ksi centre;
  312. trace ksi mt (cont mt) TITR ' Facteur d allongemet du maillage';
  313. trace MT;
  314. Finsi;
  315.  
  316. *--------------------------- maillage ----------------------------------
  317.  
  318.  
  319.  
  320. *------------------------ donnees physiques ----------------------------
  321. *
  322. NUF = 1.5E-5 ;
  323. REJ = 2.e5 ;
  324. REJ = 100. ;
  325. *REJ = 30. ;
  326. REJ = 2.e4 ;
  327. REJ = 3000. ;
  328. UJ0 = REJ*NUF/DJ ; mess 'vitesse d injection (m/s) =' UJ0;
  329. DT=CFL*RJ/UJ0 ;
  330. mess ' CFL inj = ' CFL ' DT=' DT;
  331. tasser jet;
  332. xjet=coor 1 jet;
  333. Dgj=doma $jet 'XXDIAGSI';
  334. a=(SOMT (Dgj*((RJ - xjet)**(1./7.))))*(1./PI/RJ/RJ);
  335. UJ=UJ0/a*((RJ - xjet)**(1./7.));
  336. *dess (evol chpo uj jet);
  337.  
  338. *opti donn 5;
  339. KJ = (0.02*UJ0)**2. ; mess 'KJ =' KJ ;
  340. EJ =(KJ**1.5)/DJ ; mess 'EJ =' EJ ;
  341. NUTj = 0.09*KJ*KJ/EJ ; mess 'NUTJ =' nutj ;
  342. KA = KJ*1.E-5 ;
  343. NUTAsnu=1.e-4;
  344. cnu=0.09;
  345. EA=cnu*ka*ka/nuf;
  346.  
  347. *-------------------------- equations ----------------------------------
  348. *$out = $haut;
  349. RV = EQEX $MT 'DUMP' 'ITMA' NBIT
  350. 'OPTI' 'EF' 'SUPG' 'IMPL' KPRESS
  351. * 'ZONE' $MT OPER OUTFLOW $out 1. 'UN' 'MUF' 'INCO' 'KN'
  352. * 'ZONE' $MT OPER OUTFLOW $out 1. 'UN' 'MUF' 'INCO' 'EN'
  353. 'ZONE' $MT OPER OUTFLOW $out 1. 'UN' 'MUF' 'INCO' 'UN'
  354. 'ZONE' $MT 'OPER' 'KEPSILON' 1. 'UN' NUF 'DT' 'INCO' 'KN' 'EN'
  355. 'ZONE' $MT 'OPER' 'NS' 1. 'UN' 'MUF' 'INCO' 'UN';
  356. RV = EQEX RV
  357. 'OPTI' 'EFM1' 'CENTREE'
  358. 'ZONE' $MT 'OPER' DFDT 1. 'UN' 'DT' 'INCO' 'UN'
  359. ;
  360.  
  361. RV = EQEX RV
  362. 'CLIM' 'UN' UIMP JET 0. 'UN' VIMP JET UJ
  363. 'UN' VIMP BAS2 0. 'UN' UIMP AXE 0.
  364. 'KN' TIMP JET KJ 'EN' TIMP JET EJ
  365. * 'UN' VIMP BORD 0.
  366. * 'KN' TIMP Bord 0. 'EN' TIMP Bord 1.e-10
  367. 'KN' TIMP BAS2 0. 'EN' TIMP BAS2 1.e-10;
  368. ;
  369. *RV.'ALGO_KEPSILON'= MOTS 'RNG';
  370. *RV.'ALGO_KEPSILON'= MOTS 'Ret';
  371.  
  372. RVP = EQEX 'OPTI' 'EF' KPRESS
  373. 'ZONE' $MT 'OPER' 'KBBT' -1. 'INCO' 'UN' 'PRES' ;
  374.  
  375. RV.'PROJ'= RVP ;
  376.  
  377. *------------------------ initialisations ------------------------------
  378.  
  379. RV.INCO = TABLE 'INCO' ;
  380. RV.'INCO'.'UN' = KCHT $MT VECT SOMMET (1.E-7 1.E-7) ;
  381. RV.'INCO'.'PRES' = 'KCHT' $MT 'SCAL' KPRESS 0.;
  382. RV.'INCO'.'KN' = KCHT $MT SCAL SOMMET 1.E-7 ;
  383. RV.'INCO'.'EN' = KCHT $MT SCAL SOMMET 1.E-5 ;
  384. RV.'INCO'.'NUT' = KCHT $MT SCAL CENTRE 1.E-6 ;
  385. RV.'INCO'.'MUF' = KCHT $MT SCAL sommet nuf ;
  386. RV.'INCO'.'DT' = DT ;
  387.  
  388. *------------------------ historiques ----------------------------------
  389.  
  390. P11 = MT POIN 'PROC' ((25.*RJ) (10.*DJ));
  391. P12 = MT POIN 'PROC' ((25.*RJ) (20.*DJ));
  392. P14 = MT POIN 'PROC' ((25.*RJ) (40.*DJ));
  393. P15 = MT POIN 'PROC' ((25.*RJ) (50.*DJ));
  394. LH = P02 et P03 et P04 et P05 et P11 et P12 et P14 et P15 ;
  395. HIS = KHIS 'UN' 1 LH 'UN' 2 LH 'KN' LH 'EN' LH ;
  396. RV.'HIST' = HIS ;
  397.  
  398. *------------------------ resolution -----------------------------------
  399.  
  400. lh= (POIN MT PROC (0.05 0.01)) et (POIN MT PROC (0.05 0.05))
  401. et (POIN MT PROC (0.05 0.1)) et (POIN MT PROC (0.05 0.15))
  402. et (POIN MT PROC (0.05 0.2)) et (POIN MT PROC (0.05 0.7))
  403. et (POIN MT PROC (0.05 0.3)) et (POIN MT PROC (0.05 0.6))
  404. et (POIN MT PROC (0.05 0.4)) et (POIN MT PROC (0.05 0.5));
  405.  
  406. his=khis 'UN' 1 lh
  407. 'UN' 2 lh
  408. 'KN' lh
  409. 'EN' lh;
  410. his.'KFIH'=1;
  411. rv.'HIST'=his;
  412.  
  413. EXEC rv;
  414.  
  415. *opti sauv 'JETKEI1';
  416. *sauv;
  417. *------------------------ post-traitement ------------------------------
  418. UN = (RV.'INCO'.'UN');
  419. uaxe = ((exco UN 'UY') redu axe)*(1./UJ0);
  420. evuaxe= evol chpo uaxe axe ;
  421.  
  422. Si graph
  423. trace ((cont mt) et Ml20d et Ml30d et Ml40d);
  424.  
  425. dessin his.'TABD' his.'KN' TITR ' historiques: k';
  426. dessin his.'TABD' his.'EN' TITR ' historiques: epsilon';
  427. dessin his.'TABD' his.'1UN' TITR ' historiques: ux';
  428. dessin his.'TABD' his.'2UN' TITR ' historiques: uy';
  429.  
  430. UNV = VECT UN (0.1/UJ0) UX UY JAUNE;
  431. TRAC UNV CMT TITRE 'VITESSES AIR' ;
  432.  
  433. dess evuaxe TITR 'Vitesse sur l axe';
  434.  
  435. z= (coor 2 axe) + 1.e-3;
  436. vaxe= 5.8 * Dj * (inve z);
  437. evax = evol chpo vaxe axe ;
  438. dess (evax et evuaxe) ybor 0.1 1.2 TITR ' Vitesses sur l axe';
  439.  
  440. u20d = (exco un 'UY') redu Ml20d;
  441. um20d = maxi u20d;
  442. u20d = evol chpo (u20d*(1./um20d)) ml20d ;
  443.  
  444. u30d = (exco un 'UY') redu Ml30d;
  445. um30d = maxi u30d;
  446. u30d = evol chpo (u30d*(1./um30d)) ml30d ;
  447.  
  448. u40d = (exco un 'UY') redu Ml40d;
  449. um40d = maxi u40d;
  450. u40d = evol chpo (u40d*(1./um40d)) ml40d ;
  451.  
  452. z20d = maxi (coor 2 ml20d);
  453. z30d = maxi (coor 2 ml30d);
  454. z40d = maxi (coor 2 ml40d);
  455. r20d = (extr u20d 'ABSC')*(1./z20d);
  456. r30d = (extr u30d 'ABSC')*(1./z30d);
  457. r40d = (extr u40d 'ABSC')*(1./z40d);
  458.  
  459. u20d = evol manu r20d (extr u20d 'ORDO');
  460. u30d = evol manu r30d (extr u30d 'ORDO');
  461. u40d = evol manu r40d (extr u40d 'ORDO');
  462.  
  463. dess (u20d et u30d et u40d) XBOR 0. 0.4
  464. TITR ' Profil radiaux de vitesse';
  465.  
  466.  
  467. nutsnu = rv.inco.'MUF' * (1./NUF) ;
  468. trace nutsnu mt (cont mt) TITR 'Muf / Mu' ;
  469. trace (log nutsnu) mt (cont mt) TITR ' Log (Muf / Mu)' ;
  470.  
  471. trace rv.inco.'KN' mt (cont mt) TITR ' Kn ';
  472. trace (log rv.inco.'KN') mt (cont mt) TITR ' Log Kn ';
  473. trace rv.inco.'EN' mt (cont mt) TITR ' En ';
  474. trace (log rv.inco.'EN') mt (cont mt) TITR ' Log En ';
  475. Finsi;
  476. *-------------------- fin post-traitement ------------------------------
  477.  
  478. Si (NON COMPLET);
  479. uaxref=prog
  480. 1.8650 1.8896 1.9826 1.8759 1.4946
  481. 1.1368 0.91643 0.74427 0.60911 0.51210
  482. 0.43502 0.37462 0.32726 0.28785 0.25592
  483. 0.23163 0.21113 0.19329 0.17719 0.16177
  484. 0.14722 0.13685 0.13161 0.12988 0.12893
  485. 0.12679 0.12180 0.11305;
  486. zaxe = extr evuaxe 'ABSC';
  487. list (extr evuaxe 'ORDO');
  488. evuaxer=evol manu zaxe uaxref ;
  489. m=(INTG evuaxer 'ABSO') SOMM ;
  490. delt2= SOMM (INTG (evuaxer - evuaxe) 'ABSO')/m;
  491. mess ' deltaL2 : m=' m ' delt2=' delt2;
  492. Si (delt2 > 1.e-3);ERREUR 5;Finsi;
  493. FINSI;
  494.  
  495. FIN ;
  496.  
  497.  
  498.  
  499.  
  500.  
  501.  
  502.  

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