Télécharger kepsilon.procedur

Retour à la liste

Numérotation des lignes :

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

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