Télécharger la plaquette Cast3M


fr en

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