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 :
© Cast3M 2003 - All rights reserved.
Disclaimer