Télécharger kepsilon.procedur

Retour à la liste

Numérotation des lignes :

  1. * KEPSILON PROCEDUR CB215821 20/11/04 21:18:38 10766
  2. 'DEBPROC' KEPSILON ;
  3. ARGU RX*TABLE ;
  4. ************************************************************************
  5. * Ro UN Mu DT (GB T)
  6. *
  7. * EN
  8. * KN
  9. * -> CHPOINT générés dans la table inco
  10. * TKTE teta=k/epsilon
  11. * TETK i =epsilon/k alias TKTE
  12. * NUTI valeur intermédiaire de NUT
  13. * FI inconnue Fi
  14. * PRODT Prodution turbulente
  15. * TKTI teta=k/epsilon intermédiaire
  16. * Ksi facteur de déséquilibre : nut P / epsilon
  17. * MUF viscosité dynamique effective
  18. * (tourbillonnaire+moléculaire)
  19. ************************************************************************
  20. rv=rx.'EQEX' ;
  21. iarg=rx.'IARG' ;
  22. *NASTOK = rv.'NAVISTOK' ;
  23. $mod=rx.'DOMZ' ;
  24. Dg=doma $mod 'XXDIAGSI' ;
  25. *mess ' DEBUT KEPSILON';
  26.  
  27. * Lecture du 1er Argument la densité
  28. Si(ega ('TYPE' rx.'ARG1') 'MOT ');
  29. Ro = rv.inco.(rx.'ARG1');
  30. Sinon ;
  31. Ro = rx.'ARG1';
  32. Finsi ;
  33.  
  34. * Lecture du 3ème Argument la viscosité cinématique
  35. Si(ega ('TYPE' rx.'ARG3') 'MOT ');
  36. Mu = rv.inco.(rx.'ARG3');
  37. Sinon ;
  38. Mu = rx.'ARG3';
  39. Finsi ;
  40. Mus2= Mu*(0.5);
  41.  
  42. Si (EGA (TYPE Ro) 'FLOTTANT');
  43. iRo = 1. / Ro ;
  44. Sinon ;
  45. iRo = inve Ro;
  46. Finsi ;
  47.  
  48. Si (EGA (TYPE Mu) 'FLOTTANT');
  49. Mum = Mu ;
  50. iMu = 1./MU ;
  51. Sinon ;
  52. Mum = (Maxi Mu) + (Mini Mu) * 0.5;
  53. iMu = inve MU ;
  54. Finsi ;
  55.  
  56. Nu = Mu * iRo;
  57. iNu = Ro * iMu;
  58.  
  59. Si (EGA (TYPE Nu) 'FLOTTANT');
  60. Num = Nu*0.5 ;
  61. Sinon ;
  62. Num = (Maxi Nu) + (Mini Nu) * 0.25;
  63. Finsi ;
  64.  
  65. Si( (non(ega rx.'KOPT'.'KIMPL' 1)) et
  66. (non(ega rx.'KOPT'.'KFORM' 1)) );
  67. mess ' KEPSILON KIMPL' (rx.'KOPT'.'KIMPL');
  68. mess ' KEPSILON KFORM' (rx.'KOPT'.'KFORM');
  69. mess 'Options non prevues IMPL obligatoire ' ;
  70. ERRE 21 ;
  71. Finsi ;
  72.  
  73. ******* Options du K-epsilon *******************************************
  74. EDPNUT=VRAI;Tmin=FAUX;
  75. EDPFI =FAUX;
  76. Kbw=FAUX ;KCnu=FAUX;RNG=FAUX;Filtre=FAUX;CSTE=FAUX;M2M=FAUX;Kimpr=FAUX;
  77. KRet=FAUX;KLbr=FAUX;Kchn=FAUX;Brls=FAUX;Kkl=FAUX;Kbrey=FAUX;
  78. Brjl=FAUX;Brlb=FAUX;
  79. V0=FAUX;TETA=VRAI;PERI=FAUX;
  80.  
  81. Si ('EXIST' rv 'ALGO_KEPSILON');
  82. Si('EGA' (TYPE rv.'ALGO_KEPSILON') LISTMOTS);
  83. ko=0;
  84. lko=dime rv.'ALGO_KEPSILON';
  85. Si ('EXIST' rv.'ALGO_KEPSILON' 'IMPR') ; Kimpr =VRAI ;ko=ko+1; Finsi;
  86. Si ('EXIST' rv.'ALGO_KEPSILON' 'RNG') ; RNG =VRAI ;ko=ko+1; Finsi;
  87. Si ('EXIST' rv.'ALGO_KEPSILON' 'Filtre') ; Filtre =VRAI ;ko=ko+1; Finsi;
  88. Si ('EXIST' rv.'ALGO_KEPSILON' 'Bw') ; Kbw =VRAI ;ko=ko+1; Finsi;
  89. Si ('EXIST' rv.'ALGO_KEPSILON' 'Cnu') ; KCnu =VRAI ;ko=ko+1; Finsi;
  90. Si ('EXIST' rv.'ALGO_KEPSILON' 'Nut') ; EDPNUT =VRAI ;EDPFI=FAUX ;ko=ko+1; Finsi;
  91. Si ('EXIST' rv.'ALGO_KEPSILON' 'Fi') ; EDPNUT =FAUX ;EDPFI=VRAI ;ko=ko+1; Finsi;
  92. Si ('EXIST' rv.'ALGO_KEPSILON' 'M2M') ; M2M =VRAI ;ko=ko+1; Finsi;
  93. Si ('EXIST' rv.'ALGO_KEPSILON' 'CSTE') ; CSTE =VRAI ;ko=ko+1; Finsi;
  94. Si ('EXIST' rv.'ALGO_KEPSILON' 'KL') ; EDPNUT =FAUX ;EDPFI=FAUX; Kkl=VRAI;ko=ko+1; Finsi;
  95. Si ('EXIST' rv.'ALGO_KEPSILON' 'KLbr') ; KLbr =VRAI ;ko=ko+1; Finsi;
  96. Si ('EXIST' rv.'ALGO_KEPSILON' 'Chien') ; Kchn =VRAI ;ko=ko+1; Finsi;
  97. Si ('EXIST' rv.'ALGO_KEPSILON' 'Sharma') ; Brls =VRAI ;ko=ko+1; Finsi;
  98. Si ('EXIST' rv.'ALGO_KEPSILON' 'Jones') ; Brjl =VRAI ;ko=ko+1; Finsi;
  99. Si ('EXIST' rv.'ALGO_KEPSILON' 'Lam' ) ; Brlb =VRAI ;ko=ko+1; Finsi;
  100.  
  101. COMM 'Les trois suivants ne sont pas dans la notice de KEPSILON !!! ';
  102. Si ('EXIST' rv.'ALGO_KEPSILON' 'V0') ; V0 =VRAI ;ko=ko+1; Finsi;
  103. Si ('EXIST' rv.'ALGO_KEPSILON' 'Ret') ; KRet =VRAI ;ko=ko+1; Finsi;
  104. Si ('EXIST' rv.'ALGO_KEPSILON' 'Perio') ; PERI =VRAI ;ko=ko+1; Finsi;
  105.  
  106. Si (NON (EGA ko lko));
  107. Mess ' ' ;Mess ' ' ;
  108. Mess '********************************************************************';
  109. Mess '******* Il y a une option invalide pour le modele K-Epsilon *******';
  110. Mess '******* ou bien elle apparait plusieurs fois *******';
  111. Mess '******* Liste de rv. ALGO_KEPSILON pour controle *******';
  112. Mess ' ';
  113. list rv.'ALGO_KEPSILON';
  114. Mess '********************************************************************';
  115. ERRE 21 ;
  116. Finsi ;
  117. Finsi ;
  118.  
  119. Sinon ;
  120. Si Kimpr;
  121. mess '*********** Utilisation Standard du modele K-Epsilon **********';
  122. Finsi ;
  123. Finsi ;
  124.  
  125.  
  126. KFORM=rx.'KOPT'.'KFORM' ;
  127. IDCEN=rx.'KOPT'.'IDCEN' ;
  128. CMD =rx.'KOPT'.'CMD' ;
  129. IDIV =rx.'KOPT'.'IDIV' ;
  130.  
  131. rx.'KOPT'.'IDIV' =0 ;
  132. rx.'KOPT'.'IDCEN' =2 ;
  133. * rx.'KOPT'.'CMD' =0.5;
  134. rx.'KOPT'.'CMD' =0.2;
  135. NBTETA=10;
  136.  
  137. Si V0 ;
  138. rx.'KOPT'.'CMD' =0.2 ; EDPNUT=VRAI;EDPFI=FAUX;
  139. Tmin=FAUX;Kbw=FAUX;KCnu=FAUX;M2M=FAUX;CSTE=FAUX;
  140. Filtre=FAUX;Kret=FAUX;KLbr=FAUX;Kchn=FAUX;Brls=FAUX;Brjl=FAUX;
  141. Brlb=FAUX;PERI =FAUX;
  142. TETA=VRAI
  143. Finsi ;
  144.  
  145. Si(Kkl ou KLbr ou Kchn ou Brls ou Brjl ou Brlb);
  146. Kbrey= VRAI;
  147. Tmin = FAUX;
  148. Finsi;
  149.  
  150. Si Kimpr ;
  151. lopt = chai ' ' ;
  152. Si V0 ; lopt = chai lopt ' V0 ' ; Finsi ;
  153. Si EDPFI ; lopt = chai lopt ' FI ' ; Finsi ;
  154. Si EDPNUT ; lopt = chai lopt ' NUT ' ; Finsi ;
  155. Si RNG ; lopt = chai lopt ' RNG ' ; Finsi ;
  156. Si CSTE ; lopt = chai lopt ' CSTE ' ; Finsi ;
  157. Si M2M ; lopt = chai lopt ' M2M ' ; Finsi ;
  158. Si Kbw ; lopt = chai lopt ' Bw ' ; Finsi ;
  159. Si KCnu ; lopt = chai lopt ' Cnu ' ; Finsi ;
  160. Si KRet ; lopt = chai lopt ' Ret ' ; Finsi ;
  161. Si Kkl ; lopt = chai lopt ' KL ' ; Finsi ;
  162. Si KLbr ; lopt = chai lopt ' KLbr' ; Finsi ;
  163. Si Kchn ; lopt = chai lopt ' Chien ' ; Finsi ;
  164. Si Brls ; lopt = chai lopt ' Sharma' ; Finsi ;
  165. Si Brjl ; lopt = chai lopt ' Jones ' ; Finsi ;
  166. Si Brlb ; lopt = chai lopt ' Lam ' ; Finsi ;
  167. Si Peri ; lopt = chai lopt ' Peri ' ; Finsi ;
  168. Si Filtre ; lopt = chai lopt ' Filtre (URANS) ' ; Finsi ;
  169. Si ('EXIST' rv 'ALGO_KEPSILON');
  170. mess '* K-epsilon : Options : ' lopt ;
  171. Finsi ;
  172. Finsi ;
  173.  
  174. IPT=rv.'PASDETPS'.'NUPASDT';
  175.  
  176. ******* Options du K-epsilon * Fin *************************************
  177.  
  178. ************** Constantes diverses **********
  179. * Constante de Bradshaw Bw de Smagorinski Cs
  180. Bw=0.3 ;
  181. Cs=0.1 ;
  182. Ka=0.41;
  183. * Paramètre de Menter Nut P < Ksim * epsilon (ksim=10)
  184. ksim=10.;
  185.  
  186. ************** Constantes du K epsilon ****************
  187. cnuo=0.09;
  188. c1=1.44;
  189. c2=1.92;
  190. sgk=1. ;
  191. sge=1.3;
  192.  
  193. Si (Kkl ou KLbr);
  194. ****** Constantes du modèle K-L ou K-L Bas Reynolds ***
  195. cnuo=0.09;
  196. c1=1.44;
  197. cl=1./Ka;
  198. cd=0.9;
  199. Aetmu=63.;
  200. Aeteps=3.8;
  201. sgk=0.9;
  202. finsi;
  203.  
  204. Si Kchn;
  205. ****** Constantes K epsilon Bas Reynolds Chien ********
  206. cnuo=0.09;
  207. c1=1.35;
  208. c2=1.8 ;
  209. sgk=1. ;
  210. sge=1.3;
  211. finsi;
  212.  
  213. Si RNG ;
  214. ************** Constantes du RNG K epsilon ************
  215. cnuo=0.0845;
  216. c1=1.42;
  217. c2=1.68;
  218. sgk=0.7179;
  219. sge=0.7179;
  220. etao=4.377;
  221. beta=0.012;
  222. etai=1./etao;
  223. finsi;
  224.  
  225. Si M2M ;
  226. ************** Constantes de Mohammadi Medic **********
  227. cnuo=0.09;
  228. c1=0.1296;
  229. c2=11./6.;
  230. sgk=1.;
  231. sge=Ka*Ka*(cnuo**0.5)/(c2*cnuo - c1);
  232. finsi;
  233.  
  234. Si CSTE;
  235. ************** Constantes lues dans rv.inco ***********
  236. cnuo=rv.inco.'cnu';
  237. c2=rv.inco.'c2';
  238. sgk=rv.inco.'sgk';
  239. sge=rv.inco.'sge';
  240. c1=c2 - (Ka*Ka/sge/(cnuo**0.5)) ;
  241. finsi;
  242.  
  243. Si Kimpr;
  244. mess '* K-epsilon, Constantes: cnuo=' cnuo ' c1=' c1 ' c2=' c2;
  245. mess '* K-epsilon, Constantes: sgk=' sgk ' sge=' sge;
  246. Finsi;
  247.  
  248. rcnu = cnuo**0.5 ;
  249.  
  250. * : mp=-9. ; mq = 5. ;
  251. * : mp= 3. ; mq =-2. ;
  252. * : mp=-17.; mq = 9. ;
  253. * : mp=(-1.)*c1; mq = 1. ;
  254. * : mp=2. ; mq =-1. ;
  255. * : mp=-2. ; mq = 1. ;
  256. * : mp=(-1.)*c2; mq = 1. ;
  257. * Fi : mp=-3. ; mq = 2. ;
  258.  
  259. **************** INITIALISATIONs ***********************
  260. *TYPINV = 4 ; cette option est prematureeé
  261. TYPINV = 3 ;
  262.  
  263. Si(non (exist rx 'rxk'));
  264.  
  265. rxt1=table 'KIZX';
  266. rxt2=table 'KIZX';
  267. rxk=table 'KIZX';
  268. rxe=table 'KIZX';
  269. rx.'rxt1'=rxt1;
  270. rx.'rxti'=rxti;
  271. rx.'rxt2'=rxt2;
  272. rx.'rxk'=rxk;
  273. rx.'rxe'=rxe;
  274.  
  275. mtinv = (eqex) . 'METHINV' ;
  276. rv.'METH_KEP'=mtinv;
  277. mtinv.'TYPINV'=TYPINV;
  278. mtinv.'IMPINV'=1;
  279. mtinv.'TYRENU'='SLOA';
  280. mtinv.'PCMLAG'='APR2';
  281. mtinv.'OUBMAT' = 0 ;
  282. mtinv.'SCALING' = 0 ;
  283. mtinv.'NITMAX'=500;
  284. mtinv.'RESID'=1.e-12;
  285. mtinv.'PRECOND'=3;
  286. mtinv.'FCPRECT'=1;
  287. mtinv.'FCPRECI'=1;
  288. mtinv.'BCGSBTOL'=1.D-40;
  289.  
  290. *********** étape de diffusion de K *********************
  291. rxk.'EQEX'=rx.'EQEX';
  292. rxk.'NOMZONE'=' ';
  293. rxk.'DOMZ'=rx.'DOMZ';
  294. rxk.'TDOMZ'=0;
  295. rxk.'NOMOPER'=mot 'LAPN';
  296. rxk.'KOPT'=rx.'KOPT';
  297. rxk.'IARG'=1 ;
  298. rxk.'ARG1'=1.e-10;
  299.  
  300. *********** étape de diffusion de Epsilon ***************
  301. rxe.'EQEX'=rx.'EQEX';
  302. rxe.'NOMZONE'=' ';
  303. rxe.'DOMZ'=rx.'DOMZ';
  304. rxe.'TDOMZ'=0;
  305. rxe.'NOMOPER'=mot 'LAPN';
  306. rxe.'KOPT'=rx.'KOPT';
  307. rxe.'IARG'=1 ;
  308. rxe.'ARG1'=1.e-10;
  309.  
  310. Si (NON (EGA (dime rx.'LISTINCO') 2));
  311. mess ' Erreur dans la procedure Kepsilon !';
  312. mess ' Il doit y avoir deux inconnues !';
  313. ERRE 21 ;
  314. Finsi;
  315. MI1=EXTR rx.'LISTINCO' 1 ;
  316. MI2=EXTR rx.'LISTINCO' 2 ;
  317. rxk.'LISTINCO'=MOTS MI1 ;
  318. rxe.'LISTINCO'=MOTS MI2 ;
  319. Ko=(abs rv.inco.MI1) + 1.e-10;
  320. Eo=(abs rv.inco.MI2) + 1.e-10;
  321. To=Ko*(inve Eo);
  322. rv.inco.'TKTE' =To ;
  323.  
  324. Fio=inve (To*To*Ko);
  325. rv.inco.'Fio' =Fio;
  326.  
  327. ************************ dT/dt **************************
  328.  
  329. rxt1.'EQEX'=rx.'EQEX';
  330. rxt1.'NOMZONE'=' ';
  331. rxt1.'DOMZ'=rx.'DOMZ';
  332. rxt1.'TDOMZ'=0;
  333. rxt1.'NOMOPER'=mot 'TSCA_TK';
  334. rxt1.'KOPT'=rx.'KOPT';
  335. rxt1.'LISTINCO'=MOTS 'TKTE' ;
  336. rxt1.'IARG'=3;
  337. *rxt1.'ARG1'=Nu ;
  338. rxt1.'ARG2'=rx.'ARG2';
  339. rxt1.'ARG3'=0. ;
  340.  
  341. rxt2.'EQEX'=rx.'EQEX';
  342. rxt2.'NOMZONE'=' ';
  343. rxt2.'DOMZ'=rx.'DOMZ';
  344. rxt2.'TDOMZ'=0;
  345. rxt2.'NOMOPER'=mot 'TSCA_TK';
  346. rxt2.'KOPT'=rx.'KOPT';
  347. *rxt2.'LISTINCO'=MOTS 'NUTI' ;
  348. rxt2.'LISTINCO'=MOTS 'NUTI' ;
  349. rxt2.'IARG'=3;
  350. *rxt2.'ARG1'=Nu;
  351. rxt2.'ARG2'=rx.'ARG2';
  352. rxt2.'ARG3'=0. ;
  353. Si (NON ('EXIST' (rv.'INCO') 'NUTI'));
  354. rv.inco.'NUTI'=kcht $mod scal sommet Nu;
  355. Finsi ;
  356. Si (NON ('EXIST' (rv.'INCO') 'FI'));
  357. rv.inco.'FI'=kcht $mod scal sommet 0.;
  358. Finsi ;
  359. Sinon;
  360.  
  361. rxk=rx.'rxk';
  362. mtinv=rv.'METH_KEP';
  363. rxe=rx.'rxe';
  364. rxt1=rx.'rxt1';
  365. rxt2=rx.'rxt2';
  366.  
  367. ******** On crée les matrices de périodicité ************
  368. Si Peri;
  369. *st1 mat1 = KOPS 'MATRIK' ;
  370. *----------------------------------------
  371. mr1=rv.'INCO'.'M1PERIODIC';
  372. mr2=rv.'INCO'.'M2PERIODIC';
  373. mati=rela 'UX' mr1 - 'UX' mr2;
  374. sti =depi mati 0.;
  375. mati=kops 'RIMA' mati;
  376. rxk.'PERIODIC_mat'=mati;
  377. rxk.'PERIODIC_st' =sti;
  378. *----------------------------------------
  379. Finsi;
  380. **** Fin On crée les matrices de périodicité ************
  381. Finsi ;
  382.  
  383. ********************** Fin Initialisations *******************
  384.  
  385. Si('EXIST' rxk 'PERIODIC_mat');
  386. matipc=rxk.'PERIODIC_mat';
  387. stipc =rxk.'PERIODIC_st' ;
  388. Finsi;
  389.  
  390. MI1=EXTR rxk.'LISTINCO' 1 ;
  391. MI2=EXTR rxe.'LISTINCO' 1 ;
  392.  
  393. MEN=FAUX ;
  394. Si (exist rv.'CLIM' MI1);
  395. climk=nomc MI1 (exco MI1 rv.clim) ;
  396. climk=((abs climk) + 1.e-10);
  397. Si(exist rv 'climk');
  398. xo=nomc MI1 (exco MI1 rv.'CLIM');
  399. x1=redu rv.'climk' (extr xo MAILLAGE) ;
  400. x2=rv.'climk' - x1;
  401. climk=x2 + xo;
  402. Finsi ;
  403. MEN=VRAI;
  404. Sinon ;
  405. Si(exist rv 'climk');
  406. climk=abs (rv.'climk');
  407. Sinon;
  408. mess ' Erreur dans la procedure Kepsilon !';
  409. mess ' Une condition limite sur K est necessaire !';
  410. ERRE 21 ;
  411. Finsi ;
  412. Finsi ;
  413.  
  414. Si (exist rv.'CLIM' MI2);
  415. clime=nomc MI2 (exco MI2 rv.clim) ;
  416. clime=((abs clime) + 1.e-10);
  417. MEN=VRAI;
  418. Si(exist rv 'clime');
  419. xo=nomc MI2 (exco MI2 rv.'CLIM');
  420. x1=redu rv.'clime' (extr xo MAILLAGE) ;
  421. x2=rv.'clime' - x1;
  422. clime=x2 + xo;
  423. Finsi ;
  424. Sinon ;
  425. Si(exist rv 'clime');
  426. clime=abs (rv.'clime');
  427. Sinon;
  428. mess ' Erreur dans la procedure Kepsilon !';
  429. mess ' Une condition limite sur Epsilon est necessaire !';
  430. ERRE 21 ;
  431. Finsi ;
  432. Finsi ;
  433.  
  434. Si MEN ;
  435. lic=extr rv.clim 'COMP' ;
  436. nbic=dime lic ;
  437. climr a = kops 'MATRIK' ;
  438.  
  439. repeter bloc1 nbic ;
  440. mic=extr &bloc1 lic ;
  441. Si((non (ega mic MI1)) et (non (ega mic MI2)));
  442. climr = climr et (nomc mic (exco mic (rv.clim)));
  443. Finsi ;
  444. Fin Bloc1 ;
  445. rv.'CLIM'=climr ;
  446. lic=extr climr 'COMP' ;
  447. Finsi;
  448.  
  449. rv.'climk'=climk;
  450. rv.'clime'=clime;
  451.  
  452. dtn=rx.'ARG4';
  453. tdt=type dtn;
  454. alfa=rv.'ALFA';
  455. si (ega tdt 'MOT ');
  456. * dtn = 'MOT' ('TEXTE' ('CHAINE' dtn)) ;
  457. si(ega ('MOT' dtn) ('MOT' 'DELTAT'));
  458. dtn=rv.'PASDETPS'.'DELTAT';
  459. sinon ;
  460. dtn=rv.'INCO'.dtn;
  461. finsi;
  462. finsi;
  463. dtn=dtn*alfa;
  464. * mess ' Kepsilon : Pas de temps ' dtn ;
  465.  
  466. dt=dtn ;
  467. dti=1./dt;
  468. mdt=(-1.)*dt;
  469.  
  470. *mess 'On calcule les grandeurs initiales Ko Eo To ';
  471. Ko=(abs rv.inco.MI1) + 1.e-10 ;
  472. Eo=(abs rv.inco.MI2) + 1.e-10 ;
  473.  
  474. mtinv.'XINIT'=Ko;
  475. To=kcht $mod scal sommet (Ko*(inve Eo));
  476. rv.inco.'TKTE'=To ;
  477. Fio=inve (To*To*Ko);
  478.  
  479. Io=kcht $mod scal sommet (Eo*(inve Ko));
  480. rv.inco.'TETK'=Io ;
  481.  
  482. * Muto est repris dans inco
  483. Si(exist rv.inco 'Muto');
  484. Muto=rv.inco.'Muto';
  485. Muto=Muto + 1.e-10 ;
  486. Sinon;
  487. Muto=Ro * cnuo * Ko * To ;
  488. Finsi;
  489. Nuto = Muto*iRo;
  490.  
  491. ******************** Calcul de P ***********************
  492. * mess 'AVANT CALCUL P';
  493. MIU=rx.'ARG2' ;
  494. UN=rv.inco.miu;
  495.  
  496. Si(ega iarg 4);
  497. P = PRODT UN $mod ;
  498. Sinon ;
  499. Si(ega iarg 6);
  500. si(ega ('TYPE' rx.'ARG5') 'MOT ');
  501. m4=rv.inco.(rx.'ARG5');
  502. sinon;
  503. m4=rx.'ARG5';
  504. finsi;
  505. si(ega ('TYPE' rx.'ARG6') 'MOT ');
  506. m5=rv.inco.(rx.'ARG6');
  507. sinon;
  508. m5=rx.'ARG6';
  509. finsi;
  510. P= PRODT UN $mod m4 m5 ;
  511. Sinon;
  512. mess ' Nombre d argument incorrect dans KEPSILON ';
  513. erreur KEPSILON ;
  514. Finsi;
  515. Finsi;
  516.  
  517. lc = extr UN 'COMP';
  518. mdu2= PSCA UN lc UN lc ;
  519. mdu2= mdu2 + 1.e-10;
  520.  
  521. Si (non Kbrey);
  522. * Pmax
  523. * Réalisabilité sur P dés/activée [A] (P < 10 epsilon/nut (Menter))
  524. Pmax =(ksim*Eo)*(inve Muto);
  525. a =Pmax ; al=0.7 ; ala=al*a ;
  526. * b=ala*al*((1.-al)**(-1.));
  527. b2=ala*(1.+ (ala*al) - ala)*((1.-al)**(-1.));
  528. *b3 = ala*(1.+ (ala*ala*al) - (ala*ala))*((1.-al)**(-1.));
  529. ik=masq P 'INFERIEUR' ala ;
  530. *P =(ik*P)+((1.-ik)*a*(P+b)*(inve(a+P+b)));
  531. P =(ik*P)+((1.-ik)*a*((P*P+b2)*(inve(a+(P*P+b2)))));
  532. *P =(ik*P)+((1.-ik)*a*((P*P*P+b3)*(inve(a+(P*P*P+b3)))));
  533. * Pmax
  534. Finsi;
  535.  
  536. rv.inco.'PRODT'= P ;
  537. *mess ' P = ' (maxi P ) (mini P );
  538.  
  539. rP=(ABS(P+1.e-20))**0.5 ;
  540.  
  541. ****************** FIN CALCUL P **********************
  542. ********************************************************
  543.  
  544.  
  545. Cnu = kcht $mod scal sommet (cnuo);
  546. Ak = Cnu * To * To * P ;
  547.  
  548. Si KCnu ;
  549. * On recalcule Cnu *************************************
  550. * alternative hyperbolique
  551. mess ' desequilibre Ak ' (maxi Ak) (mini Ak) ;
  552. alfa = 0.103 ; beta = 0.15 ;
  553. Cnu = alfa*(inve (Ak + beta));
  554. * alternative exponentielle Alfa = Log (cnu(0)/cnu(1))
  555. * Ak =0 -> cnu=0.7 Ak = 1 cnu=0.09 -> alfa = 2.0512
  556. * Cnu = cnuo*( (exp (2.0512*(1.-Ak))) + 1.e-10);
  557. * On recalcule Cnu Fin *********************************
  558. ********************************************************
  559. Finsi ;
  560.  
  561. **************** Calcul constante c2 *******************
  562.  
  563. *f3=1.;
  564. ATETA=(c1-1.); BTETA=c2-1.;
  565.  
  566. Si RNG;
  567. eta=rP*To;
  568. eta3=eta*eta*eta;
  569. cp2=(cnu*eta3*(1.-(etai*eta))) *
  570. (inve(1.+(beta*eta3))) ;
  571. *cc2=c2 + cp2;
  572. c2=c2 + cp2;
  573. oubli cp2; oubli eta; oubli eta3;
  574. finsi;
  575.  
  576. cnuo=0.09; c1=1.44; c2=1.92; sgk=1.; sge=1.3; alf=1.;
  577. Ret = Ro * Ko * To * iMu ;
  578. fmu = 1. ;
  579. Rocnu=Ro*Cnu;
  580. Epsok=0.;
  581. Epsoe=0.;
  582.  
  583. **************** Modèle K-L ****************************
  584.  
  585. Si Kkl;
  586. cnuo=0.09; c1=1.44; c2=1.92; sgk=1.; sge=1.3; alf=1.;
  587. fmu = 1. ;
  588. Cnu=Cnuo ;
  589. Rocnu=Ro*Cnuo ;
  590. Ret = Muto * iMu * (1./Cnu);
  591. Finsi ;
  592.  
  593. **************** Modèle K-L Bas-Reynolds ***************
  594.  
  595. Si KLbr;
  596. cnuo=0.09; c1=1.44; c2=1.92; sgk=1.; sge=1.3; alf=1.;
  597. fmu = 1. ;
  598. Cnu=Cnuo ;
  599. Rocnu=Ro*Cnu ;
  600. dparoi=rv.'INCO'.'dparoi';
  601. yet=dparoi*(Ko**0.5)*(iNu);
  602. lmu=cl*dparoi*(1. - (exp ( (-1./Aetmu)*yet)));
  603. leps=cl*dparoi*(1. - (exp ( (-1./Aeteps)*yet)));
  604. leps=leps+1.e-10;
  605. Eo=(Ko**1.5)*(inve leps);
  606. ATETA=lmu*(inve leps)/2.; BTETA=0.5;
  607. *f3=lmu*(inve leps);
  608. Finsi ;
  609.  
  610. **************** Modèle Bas-Reynolds de Chien **********
  611.  
  612. Si Kchn;
  613. cnuo=0.09; c1=1.35; c2=1.8; sgk=1.; sge=1.3;
  614. yplus=rv.'INCO'.'yplus';
  615. dparoi=rv.'INCO'.'dparoi';
  616. fmu = 1. - (exp (-0.0115*yplus));
  617. Cnu=Cnuo*fmu;
  618. Rocnu= Ro*cnu;
  619. Ret = Ro*Ko*Ko*iMu*(inve Eo);
  620. f2 = 1. - (0.22 * (exp((-1./36.)*Ret*Ret)));
  621. *E=(-2.)*Mu*Eo*(dparoi**(-2.))*(exp ((-0.5)*yplus));
  622. *epso=2.*Mu*Ko*(dparoi**(-2.));
  623. *alf=1. + ((inve Eo)*(E*To + epso));
  624.  
  625. Epsoe=2.*Mu*(dparoi**(-2.))*(exp ((-0.5)*yplus));
  626. epsok=2.*Mu*(dparoi**(-2.));
  627. alf=1.;
  628.  
  629. ATETA=(c1-1.);
  630. BTETA = (c2*f2) - alf;
  631.  
  632. Finsi ;
  633.  
  634. **************** Modèle Bas-Reynolds de Launder-Sharma *
  635.  
  636. Si Brls;
  637. lux=mots 'UX' 'UY';
  638. si(EGA (VALE DIME) 3); lux=lux et (mots 'UZ');Finsi;
  639. cnuo=0.09; c1=1.44; c2=1.92; sgk=1.; sge=1.3;
  640. Ret = Ro*Ko*Ko*iMu*(inve Eo);
  641. fmu = exp(-3.5*(inve((1.+ (Ret*(1./50.)))**2.)));
  642. Cnu = Cnuo*fmu;
  643. Rocnu=Ro*Cnu ;
  644. rKo=Ko**0.5;
  645. grk=kops 'GRADS' rKo $mod;
  646. epsok=2.*Mu*(psca grk lux grk lux)*(inve Ko);
  647. Grp=Kops 'GRADS' P $mod;
  648. D2u2=1./4.*(psca Grp lux Grp lux)*(inve P);
  649.  
  650. E=2.*mu*Muto*D2u2;
  651. f2 = 1. - (0.3 * (exp((-1.)*Ret*Ret)));
  652. alf=1. + ((inve Eo)*E*To);
  653.  
  654. ATETA=(c1-1.);
  655. BTETA = (c2*f2) - alf;
  656.  
  657. Finsi ;
  658.  
  659. **************** Modèle Bas-Reynolds de Jones-Launder **
  660.  
  661. Si Brjl;
  662. lux=mots 'UX' 'UY';
  663. si(EGA (VALE DIME) 3); lux=lux et (mots 'UZ');Finsi;
  664. cnuo=0.09; c1=1.55; c2=2.; sgk=1.; sge=1.3;
  665. Ret = Ro*Ko*Ko*iMu*(inve Eo);
  666. fmu = exp(-2.5*(inve(1.+ (Ret*(1./50.)))));
  667. Cnu = Cnuo*fmu;
  668. Rocnu=Ro*Cnu ;
  669. rKo=Ko**0.5;
  670. grk=kops 'GRADS' rKo $mod;
  671. epso=2.*Mu*(psca grk lux grk lux);
  672. Grp=Kops 'GRADS' P $mod;
  673. D2u2=1./4.*(psca Grp lux Grp lux)*(inve P);
  674. *Grp=Kops 'GRADS' (exco un 'UX') $mod;
  675. *Grp=Kops 'GRADS' (exco grp 'UY') $mod;
  676. *Grp=exco grp 'UY';
  677. *mess 'Grp ' (mini grp) (maxi grp);
  678. *D2u2=Grp*Grp;
  679.  
  680. E=2.*mu*Muto*D2u2;
  681. f2 = 1. - (0.3 * (exp((-1.)*Ret*Ret)));
  682. alf=1. + ((inve Eo)*(E*To + epso));
  683.  
  684. ATETA=(c1-1.);
  685. BTETA = (c2*f2) - alf;
  686.  
  687. Finsi ;
  688.  
  689. **************** Modèle Bas-Reynolds de Lam-Bremhorst **
  690.  
  691. Si Brlb;
  692. cnuo=0.09; c1=1.44; c2=1.92; sgk=1.; sge=1.3;
  693. Ry = (Ko**0.5)*dparoi*iMu*Ro;
  694. Ret = Muto * iMu * (1./cnuo);
  695. fmu = ((1.-(exp(-0.0165*Ry)))**2.)*(1.+(20.5*(inve Ret)));
  696. f1=1.+((0.05*(inve fmu))**3.);
  697. Cnu = Cnuo*fmu;
  698. Rocnu=Ro*Cnu ;
  699. f2 = 1. - (exp((-1.)*Ret*Ret));
  700. c1 = c1 * f1 ;
  701. c2 = c2 * f2 ;
  702. E=0.;
  703. epso=0.;
  704. alf=1. + ((inve Eo)*(E*To + epso));
  705.  
  706. ATETA=(c1-1.);
  707. BTETA = (c2*f2) - alf;
  708.  
  709. Finsi ;
  710.  
  711. ***** Solution Teta(t) et Equation sur Nuti ************
  712.  
  713. Si TETA;
  714. *................ EDP TETA .............................................
  715.  
  716. *TSCA
  717. *mess ' TSCA sur TETA' ;
  718. rxt1.'NOMOPER'=mot 'TSCA_T';
  719. rxt1.'LISTINCO'=MOTS 'TKTE' ;
  720. rxt1.'IARG'=3;
  721. rxt1.'ARG1'=Num;
  722. rxt1.'ARG2'=UN ;
  723. rxt1.'ARG3'=0. ;
  724. st mat= TSCA rxt1 ;
  725.  
  726. climt= nomc 'TKTE'
  727. ((exco MI1 (rv.'climk'))*(inve (exco MI2 (rv.'clime'))));
  728.  
  729. T=To ;
  730. alt = 1.e-1;
  731. alt = 0. ;
  732. REPETER BCLTETA NBTETA;
  733.  
  734. Stp = BTETA + (To/dt) + (cnu*ATETA*P*T*T);
  735. Stn = (2.*cnu*ATETA*P*T) + (1./dt) ;
  736.  
  737. *MDIA
  738. *mess ' MDIA sur TETA' ;
  739. rxt1.'NOMOPER'=mot 'MDIA_T';
  740. rxt1.'LISTINCO'=MOTS 'TKTE' ;
  741. rxt1.'IARG'=1;
  742. rxt1.'ARG1'=kcht $mod scal sommet (Stn+alt) ;
  743. st1 mat1 = MDIA rxt1 ;
  744. St1 = kcht $mod scal sommet comp 'TKTE'
  745. (nomc 'TKTE' (Dg*(Stp+(alt*T))));
  746. St1 = St et st1 ;
  747. mat1 = mat et mat1 ;
  748.  
  749. mtinv.'MATASS'=mat1;
  750. mtinv.'MAPREC'=mat1;
  751. mtinv.'XINIT'=nomc 'TKTE' T;
  752. mtinv.'IMPINV'=0;
  753. mtinv.'TYPINV'=TYPINV;
  754.  
  755. Si PERI ;
  756. MPI1='TKTE';
  757. mat1=mat1 et (kops 'CHANINCO' mati
  758. ('MOTS' 'LX' 'UX')
  759. ('MOTS' (chai 'LX' MPI1) MPI1)
  760. ('MOTS' 'FLX' 'FX')
  761. ('MOTS' (chai 'LX' MPI1) MPI1) );
  762. st1 =st1 et (exco sti 'FLX' (chai 'LX' MPI1));
  763. Finsi;
  764.  
  765. *mess 'Equation sur TKTE';
  766. T1='KRES' mat1 'TYPI' mtinv 'CLIM' climt 'SMBR' st1 'IMPR' 0 ;
  767. T1= nomc 'SCAL' T1;
  768.  
  769. Si Tmin;
  770. * Réalisabilité sur TETA dés/activée [A]
  771. iTetamin = (1./Num)*mdu2;
  772. *iTetamin = (1./(Num*0.01))*mdu2;
  773. a =iTetamin ; al=0.3; ala=al*a ;
  774. b=al*ala*((1.-al)**(-1.));
  775. * b2=ala*(1.+ (ala*al) - ala)*((1.-al)**(-1.));
  776. *b3 = ala*(1.+ (ala*ala*al) - (ala*ala))*((1.-al)**(-1.));
  777. iT1 = inve T1;
  778. ik=masq iT1 'INFERIEUR' ala ;
  779. iT1=(ik*iT1)+((1.-ik)*a*(iT1+b)*(inve(a+iT1+b)));
  780. *iT1=(ik*iT1)+((1.-ik)*a*((iT1*iT1+b2)*(inve(a+(iT1*iT1+b2)))));
  781. *iT1=(ik*iT1)+((1.-ik)*a*((iT1*iT1*iT1+b3)*(inve(a+(iT1*iT1*iT1+b3)))));
  782. T1 = inve iT1 ;
  783. Finsi ;
  784.  
  785. Si (non Kbrey);
  786. * Tmax
  787. * Réalisabilité sur T1 dés/activée [A] (nut P < 10 epsilon (Menter))
  788. *Tmax =(((ABS P)**0.5)*cnuo) inve ;
  789. Tmax =(inve(cnuo/ksim*(ABS P)))**0.5;
  790. a =Tmax ; al=0.7 ; ala=al*a ;
  791. * b=ala*al*((1.-al)**(-1.));
  792. b2=ala*(1.+ (ala*al) - ala)*((1.-al)**(-1.));
  793. *b3 = ala*(1.+ (ala*ala*al) - (ala*ala))*((1.-al)**(-1.));
  794. ik=masq T1 'INFERIEUR' ala ;
  795. *iE1=(ik*iE1)+((1.-ik)*a*(iE1+b)*(inve(a+iE1+b)));
  796. T1=(ik*T1)+((1.-ik)*a*((T1*T1+b2)*(inve(a+(T1*T1+b2)))));
  797. *iE1=(ik*iE1)+((1.-ik)*a*((iE1*iE1*iE1+b3)*(inve(a+(iE1*iE1*iE1+b3)))));
  798. * Tmax
  799. Finsi;
  800.  
  801. T1= (ABS T1) + 1.e-10 ;
  802. ERRT= (maxi (ABS (T1 - T)))/(maxi T1);
  803. *mess ' iter =' &BCLTETA ' ERRT=' errt ;
  804. T = T1 ;
  805.  
  806. Si (ERRT < 1.e-12) ; QUITTER BCLTETA ; Finsi ;
  807. Si ('EGA' &BCLTETA 12);
  808. mess ' Non convergence sur TETA en 10 iterations: ERRT=' errt ;
  809. Finsi ;
  810. FIN BCLTETA ;
  811. *mess ' Iteration ' &BCLTETA ' ERRT=' errt ;
  812. *............ Fin EDP TETA .............................................
  813. Finsi;
  814.  
  815.  
  816. *................ EDP Nut ..............................................
  817. Si EDPNUT ;
  818. Nutj = Nuto;
  819.  
  820. Snp= cnu*T*P*(2.-c1) ;
  821. Snn= (inve T)*((2.*alf)-c2) ;
  822.  
  823. *MDIA
  824. *mess ' MDIA sur NUTI' ;
  825. rxt2.'NOMOPER'=mot 'MDIA_N';
  826. rxt2.'LISTINCO'=MOTS 'NUTI' ;
  827. rxt2.'IARG'=1;
  828. rxt2.'ARG1'=kcht $mod scal sommet ((1./dt) + Snn);
  829. st2 mat2 = MDIA rxt2 ;
  830.  
  831. St2 = kcht $mod scal sommet comp 'NUTI'
  832. (nomc 'NUTI' (Dg*Nuto*((1./dt) + Snp)));
  833. * (nomc 'NUTI' (Dg*((Nuto*(1./dt)) + (Nutj*Snp))));
  834.  
  835. *TSCA
  836. *mess ' TSCA sur NUTI' ;
  837. rxt2.'NOMOPER'=mot 'TSCA_N';
  838. rxt2.'LISTINCO'=MOTS 'NUTI' ;
  839. rxt2.'IARG'=4;
  840. rxt2.'ARG1'=1.;
  841. rxt2.'ARG2'=UN ;
  842. rxt2.'ARG3'=(Mus2*iRo);
  843. rxt2.'ARG4'=0. ;
  844. st mat = TSCA rxt2 ;
  845. St2= St et st2 ;
  846. mat2= mat et mat2 ;
  847.  
  848. climn = nomc 'NUTI'
  849. (cnu*((abs(nomc climk 'SCAL'))**2.)*(inve(abs (nomc clime 'SCAL'))));
  850. Nutv=Nutj;
  851.  
  852. mtinv.'MATASS'=mat2;
  853. mtinv.'MAPREC'=mat2;
  854. mtinv.'XINIT'=nomc 'NUTI' Nutv;
  855. mtinv.'IMPINV'=0;
  856. mtinv.'TYPINV'=TYPINV;
  857.  
  858. Si PERI ;
  859. MPI1='NUTI';
  860. mat2=mat2 et (kops 'CHANINCO' mati
  861. ('MOTS' 'LX' 'UX')
  862. ('MOTS' (chai 'LX' MPI1) MPI1)
  863. ('MOTS' 'FLX' 'FX')
  864. ('MOTS' (chai 'LX' MPI1) MPI1) );
  865. st2 =st2 et (exco sti 'FLX' (chai 'LX' MPI1));
  866. Finsi;
  867.  
  868. *mess 'Equation sur NUTI';
  869. Nuti='KRES' mat2 'TYPI' mtinv 'CLIM' climn 'SMBR' st2 'IMPR' 0 ;
  870. Nuti=nomc 'SCAL' Nuti;
  871. Nuti = (masq Nuti 'SUPERIEUR' 0.)*Nuti + 1.e-10;
  872.  
  873. Ki = kcht $mod scal sommet (Nuti * (inve(cnu * T)));
  874.  
  875. Si(NON Kbrey);
  876. * Réalisabilité sur Ki dés/activée [A]
  877. a = 0.5*mdu2; al=0.9 ; ala=al*a ;
  878. * b=ala*al*((1.-al)**(-1.));
  879. * b2=ala*(1.+ (ala*al) - ala)*((1.-al)**(-1.));
  880. b3 = ala*(1.+ (ala*ala*al) - (ala*ala))*((1.-al)**(-1.));
  881. ik=masq Ki 'INFERIEUR' ala ;
  882. *Ki=(ik*Ki)+((1.-ik)*a*(Ki+b)*(inve(a+Ki+b)));
  883. *Ki=(ik*Ki)+((1.-ik)*a*((Ki*Ki+b2)*(inve(a+(Ki*Ki+b2)))));
  884. Ki=(ik*Ki)+((1.-ik)*a*((Ki*Ki*Ki+b3)*(inve(a+(Ki*Ki*Ki+b3)))));
  885. Finsi;
  886.  
  887. Ki = kcht $mod scal sommet (Nuti * (inve(cnu * T)));
  888. Ei = kcht $mod scal sommet (Nuti * (inve(cnu * T * T)));
  889.  
  890. Finsi ;
  891. *............ Fin EDP Nut ..............................................
  892.  
  893. *................ EDP Fi ...............................................
  894. Si EDPFI ;
  895.  
  896. mp=-3. ; mq = 2. ;
  897.  
  898. Sfn= ((inve To)*((mq*c2) + mp)) - (cnu*To*P*(mp + (mq*c1)));
  899. Sfp= 0.;
  900.  
  901. *Sfn= (cnu*To*P*(mp + (mq*c1)));
  902. *Sfp= ((inve To)*((mq*c2) + mp)) ;
  903.  
  904. *MDIA
  905. *mess ' MDIA sur FI' ;
  906. rxt2.'NOMOPER'=mot 'MDIA_F';
  907. rxt2.'LISTINCO'=MOTS 'FI' ;
  908. rxt2.'IARG'=1;
  909. rxt2.'ARG1'=kcht $mod scal sommet Sfn ;
  910. st2 mat2 = MDIA rxt2 ;
  911. St2 = kcht $mod scal sommet comp 'FI'
  912. (nomc 'FI' (Dg*Sfp));
  913.  
  914. *TSCA
  915. *mess ' TSCA sur FI' ;
  916. rxt2.'NOMOPER'=mot 'TSCA_F';
  917. rxt2.'LISTINCO'=MOTS 'FI' ;
  918. rxt2.'IARG'=3;
  919. rxt2.'ARG1'=Num;
  920. rxt2.'ARG2'=UN ;
  921. rxt2.'ARG3'=0. ;
  922. st mat= TSCA rxt2 ;
  923. St2 = St et st2 ;
  924. mat2 = mat et mat2 ;
  925.  
  926. climf= nomc 'FI'
  927. ((exco MI1 (rv.'climk'))**mp)*((exco MI2 (rv.'clime'))**mq);
  928.  
  929. mtinv.'MATASS'=mat2;
  930. mtinv.'MAPREC'=mat2;
  931. mtinv.'XINIT'=nomc 'FI' Fio;
  932. mtinv.'IMPINV'=0;
  933. mtinv.'TYPINV'=TYPINV;
  934. *mess 'Equation sur FI';
  935. Fi='KRES' mat2 'TYPI' mtinv 'CLIM' climf 'SMBR' st2 'IMPR' 0 ;
  936. Fi= nomc 'SCAL' Fi ;
  937. Fi= (ABS Fi) + 1.e-10 ;
  938.  
  939. Ki =kcht $mod scal sommet ((Fi**(1./(mp+mq)))*(T**(mq/(mp+mq))));
  940.  
  941. Si(NON Kbrey);
  942. * Réalisabilité sur Ki dés/activée [A]
  943. a = 0.5*mdu2; al=0.9 ; ala=al*a ;
  944. * b=ala*al*((1.-al)**(-1.));
  945. * b2=ala*(1.+ (ala*al) - ala)*((1.-al)**(-1.));
  946. b3 = ala*(1.+ (ala*ala*al) - (ala*ala))*((1.-al)**(-1.));
  947. ik=masq Ki 'INFERIEUR' ala ;
  948. *Ki=(ik*Ki)+((1.-ik)*a*(Ki+b)*(inve(a+Ki+b)));
  949. *Ki=(ik*Ki)+((1.-ik)*a*((Ki*Ki+b2)*(inve(a+(Ki*Ki+b2)))));
  950. Ki=(ik*Ki)+((1.-ik)*a*((Ki*Ki*Ki+b3)*(inve(a+(Ki*Ki*Ki+b3)))));
  951. Finsi;
  952.  
  953. Ei =kcht $mod scal sommet (Ki*(inve T ));
  954. Nuti=kcht $mod scal sommet (cnu*Ki*Ki*(inve Ei));
  955.  
  956. Finsi ;
  957. *............ Fin EDP Fi ...............................................
  958.  
  959. *................ modèle K-L ...........................................
  960. Si Kkl;
  961. Echl= rv.inco.'Echl';
  962. T2=T**(-2.);
  963. Ki =kcht $mod scal sommet ((Echl**2.)*(T2));
  964. ** Réalisabilité sur Ki dés/activée [A]
  965. * a = 0.5*mdu2; al=0.9 ; ala=al*a ;
  966. ** b=ala*al*((1.-al)**(-1.));
  967. ** b2=ala*(1.+ (ala*al) - ala)*((1.-al)**(-1.));
  968. * b3 = ala*(1.+ (ala*ala*al) - (ala*ala))*((1.-al)**(-1.));
  969. * ik=masq Ki 'INFERIEUR' ala ;
  970. **Ki=(ik*Ki)+((1.-ik)*a*(Ki+b)*(inve(a+Ki+b)));
  971. **Ki=(ik*Ki)+((1.-ik)*a*((Ki*Ki+b2)*(inve(a+(Ki*Ki+b2)))));
  972. * Ki=(ik*Ki)+((1.-ik)*a*((Ki*Ki*Ki+b3)*(inve(a+(Ki*Ki*Ki+b3)))));
  973. Ei = cd*(Ki**1.5)*(Echl**(-1.));
  974.  
  975. Nuti=kcht $mod scal sommet (cnu*(Ki**0.5)*Echl);
  976. Finsi;
  977. *............ Fin modèle K-L ...........................................
  978.  
  979. *................ modèle KLbr ..........................................
  980. Si KLbr;
  981. T2=T**(-2.);
  982. Ki =kcht $mod scal sommet ((leps**2.)*(T2));
  983.  
  984. Si(NON Kbrey);
  985. * Réalisabilité sur Ki dés/activée [A]
  986. a = 0.5*mdu2; al=0.9 ; ala=al*a ;
  987. * b=ala*al*((1.-al)**(-1.));
  988. * b2=ala*(1.+ (ala*al) - ala)*((1.-al)**(-1.));
  989. b3 = ala*(1.+ (ala*ala*al) - (ala*ala))*((1.-al)**(-1.));
  990. ik=masq Ki 'INFERIEUR' ala ;
  991. *Ki=(ik*Ki)+((1.-ik)*a*(Ki+b)*(inve(a+Ki+b)));
  992. *Ki=(ik*Ki)+((1.-ik)*a*((Ki*Ki+b2)*(inve(a+(Ki*Ki+b2)))));
  993. Ki=(ik*Ki)+((1.-ik)*a*((Ki*Ki*Ki+b3)*(inve(a+(Ki*Ki*Ki+b3)))));
  994. Finsi;
  995.  
  996. Ei = (Ki**1.5)*(inve leps);
  997. Nuti=kcht $mod scal sommet (cnu*(Ki**0.5)*lmu);
  998.  
  999. Finsi;
  1000. *............ Fin modèle KLbr ..........................................
  1001.  
  1002. *********************** Fin Equation sur Nuti et Teta(t) ************
  1003.  
  1004.  
  1005. ********************** On conserve les résultats intermédiaires *****
  1006. rv.inco.'TKTI'=kcht $mod scal sommet T ;
  1007. rv.inco.'NUTI' = Nuti ;
  1008.  
  1009. ****** Etape de diffusion equation sur K ***************************
  1010. *DFDT
  1011. * mess ' Etape de convection/diffusion + DFDT sur K' ;
  1012. rxk.'KOPT'.'KFORM'=0 ;
  1013. rxk.'KOPT'.'IDCEN'=1 ;
  1014. rxk.'LISTINCO'=MOTS MI1 ;
  1015. rxk.'IARG'=3;
  1016. rxk.'ARG1'=Ro;
  1017. rxk.'ARG2'=kcht $mod scal sommet Ki;
  1018. rxk.'ARG3'=rx.'ARG4';
  1019. sk1 mak1= DFDT rxk ;
  1020.  
  1021. rxk.'KOPT'.'KFORM'=KFORM;
  1022. *TSCAL
  1023. rxk.'IARG' = 1 ;
  1024. rxk.'ARG1' = kcht $mod scal sommet (Mus2 + (Muto * (1./sgk)));
  1025. sko mak = LAPN rxk ;
  1026.  
  1027. *MDIA
  1028. *mess ' MDIA pour les termes Bas-Reynolds' ;
  1029. rxk.'NOMOPER'=mot 'MDIA_K';
  1030. rxk.'IARG'=1;
  1031. rxk.'ARG1'=kcht $mod scal sommet Epsok ;
  1032. skm matk = MDIA rxk ;
  1033. *St2 = kcht $mod scal sommet comp MI1
  1034. * (nomc MI1 (Dg*Sfp));
  1035.  
  1036. sk = sk1 et sko ;
  1037. mak = mak1 et mak et matk;
  1038.  
  1039. climk= nomc MI1 (rv.'climk') ;
  1040.  
  1041. Si PERI ;
  1042. MPI1=MI1;
  1043. mak=mak et (kops 'CHANINCO' mati
  1044. ('MOTS' 'LX' 'UX')
  1045. ('MOTS' (chai 'LX' MPI1) MPI1)
  1046. ('MOTS' 'FLX' 'FX')
  1047. ('MOTS' (chai 'LX' MPI1) MPI1) );
  1048. sk =sk et (exco sti 'FLX' (chai 'LX' MPI1));
  1049. Finsi;
  1050.  
  1051. mtinv.'MATASS'=mak ;
  1052. mtinv.'MAPREC'=mak ;
  1053. mtinv.'XINIT'=nomc MI1 Ki ;
  1054. mtinv.'IMPINV'=0;
  1055. mtinv.'TYPINV'=TYPINV;
  1056.  
  1057. * mess 'Equation de diffusion sur K ';
  1058. K1='KRES' mak 'TYPI' mtinv 'CLIM' climk 'SMBR' sk 'IMPR' 0 ;
  1059. K1=(abs K1) + 1.e-10 ;
  1060. K1=kcht $mod scal sommet (nomc 'SCAL' K1);
  1061. *mess 'maxi mini K1 apres Diff ' (maxi K1) (mini K1) ;
  1062.  
  1063. ****** Etape de diffusion equation sur Epsilon *********************
  1064. Si ((NON Kkl) et (NON KLbr));
  1065. *DFDT
  1066. *mess ' Etape de diffusion DFDT sur EPSILON' ;
  1067. rxe.'KOPT'.'KFORM'=0 ;
  1068. rxe.'KOPT'.'IDCEN'=1 ;
  1069. rxe.'LISTINCO'=MOTS MI2 ;
  1070. rxe.'IARG'=3;
  1071. rxe.'ARG1'=Ro;
  1072. rxe.'ARG2'=kcht $mod scal sommet Ei;
  1073. rxe.'ARG3'=rx.'ARG4';
  1074. se1 mae1= DFDT rxe ;
  1075.  
  1076. rxe.'KOPT'.'KFORM'=KFORM;
  1077. *TSCA
  1078. rxe.'IARG'=1;
  1079. rxe.'ARG1' = kcht $mod scal sommet (Mus2 + (Muto * (1./sge)));
  1080. seo mae = LAPN rxe ;
  1081.  
  1082. *MDIA
  1083. *mess ' MDIA pour les termes Bas-Reynolds' ;
  1084. rxe.'NOMOPER'=mot 'MDIA_E';
  1085. rxe.'IARG'=1;
  1086. rxe.'ARG1'=kcht $mod scal sommet Epsoe ;
  1087. ske maed = MDIA rxe ;
  1088. *St2 = kcht $mod scal sommet comp MI2
  1089. * (nomc MI2 (Dg*Sfp));
  1090.  
  1091. se = se1 et seo ;
  1092. mae = mae1 et mae et maed;
  1093.  
  1094. clime= rv.'clime' ;
  1095.  
  1096. Si PERI ;
  1097. MPI1=MI2;
  1098. mae=mae et (kops 'CHANINCO' mati
  1099. ('MOTS' 'LX' 'UX')
  1100. ('MOTS' (chai 'LX' MPI1) MPI1)
  1101. ('MOTS' 'FLX' 'FX')
  1102. ('MOTS' (chai 'LX' MPI1) MPI1) );
  1103. se =se et (exco sti 'FLX' (chai 'LX' MPI1));
  1104. Finsi;
  1105.  
  1106. mtinv.'MATASS'=mae ;
  1107. mtinv.'MAPREC'=mae ;
  1108. mtinv.'XINIT'=nomc MI2 Eo ;
  1109. mtinv.'IMPINV'=0;
  1110. mtinv.'TYPINV'=TYPINV;
  1111.  
  1112. * mess 'Equation de diffusion sur Epsilon ';
  1113. E1='KRES' mae 'TYPI' mtinv 'CLIM' clime 'SMBR' se 'IMPR' 0 ;
  1114. E1=(Abs E1) + 1.e-10 ;
  1115. E1=kcht $mod scal sommet (nomc 'SCAL' E1);
  1116. * mess 'maxi mini E1 apres Diff ' (maxi E1) (mini E1) ;
  1117. Finsi;
  1118.  
  1119. ****** Cas K-L *****************************************************
  1120. Si Kkl ;
  1121. * mess 'Modèle K-L ';
  1122. E1=cd*(K1**1.5)*(rv.inco.'Echl'**(-1.));
  1123. E1=kcht $mod scal sommet (nomc 'SCAL' E1);
  1124. * mess 'maxi mini E1 apres Diff ' (maxi E1) (mini E1) ;
  1125. Finsi;
  1126.  
  1127. ****** Avancement en temps ******************************
  1128.  
  1129. Si(NON Kbrey);
  1130. Mut = Rocnu * K1 * K1 * (INVE E1) ;
  1131.  
  1132. * Réalisabilité sur K1 dés/activée [A]
  1133. a = 0.5*mdu2; al=0.9 ; ala=al*a ;
  1134. * b=ala*al*((1.-al)**(-1.));
  1135. * b2=ala*(1.+ (ala*al) - ala)*((1.-al)**(-1.));
  1136. b3 = ala*(1.+ (ala*ala*al) - (ala*ala))*((1.-al)**(-1.));
  1137. ik=masq K1 'INFERIEUR' ala ;
  1138. *K1=(ik*K1)+((1.-ik)*a*(K1+b)*(inve(a+K1+b)));
  1139. *K1=(ik*K1)+((1.-ik)*a*((K1*K1+b2)*(inve(a+(K1*K1+b2)))));
  1140. K1=(ik*K1)+((1.-ik)*a*((K1*K1*K1+b3)*(inve(a+(K1*K1*K1+b3)))));
  1141.  
  1142. * Réalisabilité sur E1 dés/activée [A] (nut P < 10 epsilon (Menter))
  1143. * iEmin =10*((K1*((ABS P)**0.5)) inve) ;
  1144. iEmin =(((cnu*K1*K1*(ABS P)) inve)*ksim)**0.5 ;
  1145. a =iEmin ; al=0.7 ; ala=al*a ;
  1146. * b=ala*al*((1.-al)**(-1.));
  1147. b2=ala*(1.+ (ala*al) - ala)*((1.-al)**(-1.));
  1148. *b3 = ala*(1.+ (ala*ala*al) - (ala*ala))*((1.-al)**(-1.));
  1149. iE1 = inve E1 ;
  1150. ik=masq iE1 'INFERIEUR' ala ;
  1151. *iE1=(ik*iE1)+((1.-ik)*a*(iE1+b)*(inve(a+iE1+b)));
  1152. iE1=(ik*iE1)+((1.-ik)*a*((iE1*iE1+b2)*(inve(a+(iE1*iE1+b2)))));
  1153. *iE1=(ik*iE1)+((1.-ik)*a*((iE1*iE1*iE1+b3)*(inve(a+(iE1*iE1*iE1+b3)))));
  1154. E1 = inve iE1 ;
  1155.  
  1156. finsi;
  1157.  
  1158.  
  1159. Si Kbw ;
  1160. *mess ' Réalisabilité Bw ';
  1161. * Réalisabilité sur Nut dés/activée [A]
  1162. Mutmax = Ro * Bw * K1 * (inve (P**0.5));
  1163. a =Mutmax ; al=0.9 ; ala=al*a ;
  1164. * b=ala*al*((1.-al)**(-1.));
  1165. * b2=ala*(1.+ (ala*al) - ala)*((1.-al)**(-1.));
  1166. b3 = ala*(1.+ (ala*ala*al) - (ala*ala))*((1.-al)**(-1.));
  1167. ik=masq Mut 'INFERIEUR' ala ;
  1168. *Mut=(ik*Mut)+((1.-ik)*a*(Mut+b)*(inve(a+Mut+b)));
  1169. *Mut=(ik*Mut)+((1.-ik)*a*((Mut*Mut+b2)*(inve(a+(Mut*Nut+b2)))));
  1170. Mut=(ik*Mut)+((1.-ik)*a*((Mut*Mut*Mut+b3)*(inve(a+(Mut*Mut*Mut+b3)))));
  1171. *Teta1 = K1 * (inve E1);
  1172. mess ' MUT1?';
  1173. E1 = Rocnu * K1 * K1 * (inve Mut) ;
  1174. mess ' MUT1? OK';
  1175. Finsi ;
  1176.  
  1177. Si Filtre;
  1178. Ech = (K1**1.5)*(inve E1);
  1179. Echmax=rv.inco.'Echl';
  1180. mess ' ** Modele K-Epsilon filtre **';
  1181. a =Echmax ; al=0.9; ala=al*a ;
  1182. * b=ala*al*((1.-al)**(-1.));
  1183. * b2=ala*(1.+ (ala*al) - ala)*((1.-al)**(-1.));
  1184. b3 = ala*(1.+ (ala*ala*al) - (ala*ala))*((1.-al)**(-1.));
  1185. ik=masq Ech 'INFERIEUR' ala ;
  1186. *Ech=(ik*Ech)+((1.-ik)*a*(Ech+b)*(inve(a+Ech+b)));
  1187. *Ech=(ik*Ech)+((1.-ik)*a*((Ech*Ech+b2)*(inve(a+(Ech*Ech+b2)))));
  1188. Ech=(ik*Ech)+((1.-ik)*a*((Ech*Ech*Ech+b3)*(inve(a+(Ech*Ech*Ech+b3)))));
  1189. E1 = (K1**1.5)*(inve Ech) ;
  1190. Mut = Rocnu * K1 * K1 * (inve E1) ;
  1191. Finsi ;
  1192.  
  1193. Si (NON KLbr);
  1194. Mut = Rocnu * K1 * K1 * (INVE E1) ;
  1195. Finsi;
  1196.  
  1197. ****** Cas KLbr ****************************************************
  1198. Si KLbr ;
  1199. E1=(K1**1.5)*(leps**(-1.));
  1200. E1=kcht $mod scal sommet E1;
  1201. Mut = Rocnu * (K1**0.5) * lmu;
  1202. Finsi;
  1203. ****** Cas KLbr Fin ************************************************
  1204.  
  1205.  
  1206. dk= maxi (abs ( K1 - Ko))/(maxi Ko) ;
  1207. Ko=K1 ;
  1208. de= maxi (abs ( E1 - Eo))/(maxi Eo) ;
  1209. Eo=E1 ;
  1210. Ksi = Mut * P * (inve E1) * iRo ;
  1211. rv.inco.'Ksi'=Ksi ;
  1212. rv.inco.MI1=kcht $mod scal sommet K1;
  1213. rv.inco.MI2=kcht $mod scal sommet E1;
  1214.  
  1215.  
  1216. rv.inco.'MUF' = kcht $mod scal sommet (Mut + Mu);
  1217. rv.inco.'Muto'= kcht $mod scal sommet Mut;
  1218. rv.inco.MI2=kcht $mod scal sommet E1;
  1219.  
  1220.  
  1221. FIDT=RV.'FIDT';
  1222. INDF = IPT - (IPT/FIDT * FIDT);
  1223. Si (EGA INDF 0);
  1224. Mess ' Ecart relatif entre deux pas de temps : sur K ' dk
  1225. ' sur Epsilon ' de ;
  1226. Finsi ;
  1227. ****** Historiques **************************************
  1228. Si (exist rv 'HIST') ;
  1229. H=rv.'HIST';
  1230. Si(non(exist H 'KFIH'));
  1231. H.'KFIH'=20;
  1232. Finsi ;
  1233. KFIH=H.'KFIH';
  1234. INDH = IPT - (IPT/KFIH * KFIH);
  1235. Si (EGA INDH 0);
  1236.  
  1237. Si (exist (rv.'HIST') MI1);
  1238.  
  1239. M1=H.(CHAI '$' MI1);
  1240. EV1=H.MI1;
  1241. ltps=extr EV1 'ABSC' 1 ;
  1242. nlt=dime ltps;
  1243. Si(ega nlt 0);
  1244. tps= 0. ;
  1245.  
  1246. Sinon;
  1247.  
  1248. tps=(extr ltps nlt) + dt ;
  1249. finsi ;
  1250. ltps=ltps et (prog tps) ;
  1251. nb=nbel M1 ;
  1252. repeter bloc1 nb;
  1253. pi=poin M1 &bloc1;
  1254. lki =(extr EV1 'ORDO' &bloc1) et (prog (extr Ko 'SCAL' pi));
  1255. Si (EGA &bloc1 1);
  1256. evi = (evol 'MANU' 'TEMPS' ltps MI1 lki) ;
  1257. Sinon;
  1258. evi = evi et (evol 'MANU' 'TEMPS' ltps MI1 lki) ;
  1259. Finsi ;
  1260. Fin bloc1 ;
  1261. H.MI1=evi;
  1262.  
  1263. Finsi ;
  1264.  
  1265. Si (exist (rv.'HIST') MI2);
  1266.  
  1267. M2=H.(CHAI '$' MI2);
  1268. EV2=H.MI2;
  1269. nb=nbel M2 ;
  1270. repeter bloc1 nb;
  1271. pi=poin M2 &bloc1;
  1272. lki =(extr EV2 'ORDO' &bloc1) et (prog (extr Eo 'SCAL' pi));
  1273. Si (EGA &bloc1 1);
  1274. evi = (evol 'MANU' 'TEMPS' ltps MI2 lki) ;
  1275. Sinon;
  1276. evi = evi et (evol 'MANU' 'TEMPS' ltps MI2 lki) ;
  1277. Finsi ;
  1278. Fin bloc1 ;
  1279. H.MI2=evi;
  1280. Finsi ;
  1281. Finsi ;
  1282. Finsi ;
  1283. ****** Petit menage avant retour ************************
  1284. rx.'KOPT'.'IDIV'=IDIV ;
  1285.  
  1286. as2 ama1 = 'KOPS' 'MATRIK' ;
  1287.  
  1288. *Mess ' FIN PROC KEPSILON' ;
  1289. RESPRO as2 ama1 ;
  1290.  
  1291. FINPROC ;
  1292.  
  1293.  
  1294.  
  1295.  

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