1 : * SPAL PROCEDUR GOUNAND 11/05/24 21:16:18 6976 2 : ************************************************************************ 3 : * NOM : SPAL 4 : * DESCRIPTION : Calcule la viscosité turbulente grâce au modèle de 5 : * Spalart-Allmaras 6 : ************************************************************************ 7 : * HISTORIQUE : 6/10/2010 : JCARDO : création de la procédure 8 : * HISTORIQUE : 28/02/2011 : JCARDO : - ajout de 'DUMP' et 'METHINV' 9 : * - branchement sur le nouveau PRODT 10 : * - gestion du cas où DFDT est absent 11 : * - légère amélioration du Newton 12 : * - correction de quelques MESS 13 : * HISTORIQUE : 14 : * HISTORIQUE : 15 : ************************************************************************ 16 : * Prière de PRENDRE LE TEMPS DE COMPLÉTER LES COMMENTAIRES 17 : * en cas de modification de ce sous-programme afin de faciliter 18 : * la maintenance ! 19 : ************************************************************************ 20 : * SYNTAXE (cf. opérateur EQEX) 21 : * 22 : * 'ZONE' $MD 'OPER' 'SPAL' 'RHO' 'UN' 'MU' ('DT') 23 : * ('PERIODIC' GEOM1 GEOM2) 24 : * 'INCO' 'NU0' 25 : * 26 : ************************************************************************ 27 : 28 : 29 : DEBP SPAL ; 30 : ARGU RX*'TABLE' ; 31 : 32 : RV = RX . 'EQEX' ; 33 : DG = DOMA $MD 'XXDIAGSI' ; 34 : 35 : NDIM = VALE 'DIME' ; 36 : 37 : * Recherche-t-on directement un régime permanent (pas de DFDT) ? 38 : SI (NON (EXIS RX 'RINSTAT')) ; 39 : RX . 'RINSTAT' = FAUX ; 40 : REPE BLOCO (DIME (RV . 'LISTOPER')) ; 41 : MOPER = EXTR &BLOCO (RV . 'LISTOPER') ; 42 : SI (EGA MOPER 'DFDT ') ; 43 : RX . 'RINSTAT' = VRAI ; 44 : FINS ; 45 : FIN BLOCO ; 46 : FINS ; 47 : 48 : 49 : 50 : * +====================================================================+ 51 : * | | 52 : * | L E C T U R E D E S A R G U M E N T S | 53 : * | | 54 : * +====================================================================+ 55 : 56 : 57 : IARG = RX . 'IARG' ; 58 : 59 : SI (RX . 'RINSTAT') ; 60 : IREQ = 4 ; 61 : SINON ; 62 : IREQ = 3 ; 63 : FINS ; 64 : 65 : SI (IARG < IREQ) ; 66 : MESS '******************************************************' ; 67 : MESS '/!\ ERREUR dans SPAL :' ; 68 : MESS (CHAI ' La procédure SPAL requiert au moins ' 69 : IREQ ' arguments') ; 70 : MESS '******************************************************' ; 71 : QUIT SPAL ; 72 : FINS ; 73 : 74 : 75 : 76 : * +=================+ 77 : * | MASSE VOLUMIQUE | 78 : * +=================+ 79 : 80 : RHO = RX . 'ARG1' ; 81 : TRHO = TYPE RHO ; 82 : 83 : SI (EGA TRHO 'MOT') ; 84 : RHO = RV . 'INCO' . RHO ; 85 : TRHO = TYPE RHO ; 86 : FINS ; 87 : 88 : SI ((NEG TRHO 'FLOTTANT') ET (NEG TRHO 'CHPOINT')) ; 89 : MESS '******************************************************' ; 90 : MESS '/!\ ERREUR dans SPAL :' ; 91 : MESS ' INCO."RHO" renvoie à un objet de type incorrect' ; 92 : MESS '******************************************************' ; 93 : QUIT SPAL ; 94 : FINS ; 95 : 96 : SI (EGA TRHO 'FLOTTANT') ; 97 : USRHO = 1. / RHO ; 98 : SINON ; 99 : USRHO = INVE RHO ; 100 : FINS ; 101 : 102 : 103 : 104 : * +=====================+ 105 : * | VITESSE D'ADVECTION | 106 : * +=====================+ 107 : 108 : UN = RX . 'ARG2' ; 109 : TUN = TYPE UN ; 110 : 111 : SI (EGA TUN 'MOT') ; 112 : UN = RV . 'INCO' . UN ; 113 : TUN = TYPE UN ; 114 : FINS ; 115 : 116 : SI (NEG TUN 'CHPOINT') ; 117 : MESS '******************************************************' ; 118 : MESS '/!\ ERREUR dans SPAL :' ; 119 : MESS ' INCO."UN" doit renvoyer à un objet de type CHPOINT' ; 120 : MESS '******************************************************' ; 121 : QUIT SPAL ; 122 : FINS ; 123 : 124 : LCO = EXTR UN 'COMP' ; 125 : NCO = DIME LCO ; 126 : 127 : 128 : 129 : 130 : * +=================================+ 131 : * | VISCOSITÉ DYNAMIQUE MOLÉCULAIRE | 132 : * +=================================+ 133 : 134 : MU = RX . 'ARG3' ; 135 : TMU = TYPE MU ; 136 : 137 : SI (EGA TMU 'MOT') ; 138 : MU = RV . 'INCO' . MU ; 139 : TMU = TYPE MU ; 140 : FINS ; 141 : 142 : SI ((NEG TMU 'FLOTTANT') ET (NEG TMU 'CHPOINT')) ; 143 : MESS '******************************************************' ; 144 : MESS '/!\ ERREUR dans SPAL :' ; 145 : MESS ' INCO."MU" renvoie à un objet de type incorrect' ; 146 : MESS '******************************************************' ; 147 : QUIT SPAL ; 148 : FINS ; 149 : 150 : * Viscosité cinématique moléculaire 151 : NU = MU * USRHO ; 152 : 153 : SI (EGA TMU 'FLOTTANT') ; 154 : USNU = 1. / NU ; 155 : SINON ; 156 : USNU = INVE NU ; 157 : FINS ; 158 : 159 : 160 : 161 : 162 : * +===============================+ 163 : * | DURÉE DU PAS DE TEMPS COURANT | 164 : * +===============================+ 165 : 166 : SI (RX . 'RINSTAT') ; 167 : DT = RX . 'ARG4' ; 168 : TDT = TYPE DT ; 169 : 170 : SI (EGA TDT 'MOT') ; 171 : DT = RV . 'INCO' . DT ; 172 : TDT = TYPE DT ; 173 : FINS ; 174 : 175 : SI (NEG TDT 'FLOTTANT') ; 176 : MESS '******************************************************' ; 177 : MESS '/!\ ERREUR dans SPAL :' ; 178 : MESS ' INCO."DT" renvoie à un objet de type incorrect' ; 179 : MESS '******************************************************' ; 180 : QUIT SPAL ; 181 : FINS ; 182 : FINS ; 183 : 184 : 185 : 186 : 187 : * +===========================+ 188 : * | CONDITIONS DE PERIODICITÉ | 189 : * +===========================+ 190 : 191 : B_CYCL = FAUX ; 192 : SI (IARG > IREQ) ; 193 : MCLE = RX . (CHAI 'ARG' (IREQ+1)) ; 194 : 195 : SI (EGA MCLE 'PERIODIC') ; 196 : B_CYCL = VRAI ; 197 : 198 : GEOM1 = RX . (CHAI 'ARG' (IREQ+2)) ; 199 : GEOM2 = RX . (CHAI 'ARG' (IREQ+3)) ; 200 : 201 : TG1 = TYPE GEOM1 ; 202 : TG2 = TYPE GEOM2 ; 203 : 204 : SI ((NEG TG1 'MAILLAGE') OU (NEG TG2 'MAILLAGE')) ; 205 : MESS '******************************************************' ; 206 : MESS '/!\ ERREUR dans SPAL :' ; 207 : MESS ' On attendait 2 objets de type MAILLAGE' ; 208 : MESS '******************************************************' ; 209 : QUIT SPAL ; 210 : FINS ; 211 : SINON ; 212 : MESS '******************************************************' ; 213 : MESS '/!\ ERREUR dans SPAL :' ; 214 : MESS ' On attendait le mot-clé "PERIODIC"' ; 215 : MESS '******************************************************' ; 216 : QUIT SPAL ; 217 : FINS ; 218 : FINS ; 219 : 220 : 221 : 222 : 223 : * +=========================================+ 224 : * | NOM DE L'INCONNUE DE VISCOSITÉ MODIFIÉE | 225 : * +=========================================+ 226 : 227 : B_INCO = FAUX ; 228 : SI (EXIS RX 'LISTINCO') ; 229 : SI (EGA (DIME (RX . 'LISTINCO')) 1) ; 230 : B_INCO = VRAI ; 231 : FINS ; 232 : FINS ; 233 : 234 : SI (NON B_INCO) ; 235 : MESS '******************************************************' ; 236 : MESS '/!\ ERREUR dans SPAL :' ; 237 : MESS ' On attendait le nom donné à la seule inconnue du' ; 238 : MESS ' modèle (la viscosité modifiée)' ; 239 : MESS '******************************************************' ; 240 : QUIT SPAL ; 241 : FINS ; 242 : 243 : _NU0 = EXTR (RX . 'LISTINCO') 1 ; 244 : 245 : * Valeurs minimales admises pour NU0 et S0 (si B_POSI est VRAI) 246 : NU0_MIN = NU*1.E-5 ; 247 : SI (EGA TMU 'CHPOINT') ; 248 : NU0_MIN = MINI NU0_MIN ; 249 : FINS ; 250 : S0_MIN = 1.E-10 ; 251 : 252 : 253 : 254 : 255 : 256 : 257 : 258 : * +====================================================================+ 259 : * | | 260 : * | I N I T I A L I S A T I O N D U M O D È L E | 261 : * | | 262 : * +====================================================================+ 263 : * 264 : * Ces opérations sont effectuées uniquement lors du premier appel à 265 : * l'opérateur SPAL (ou après que RX.'INITOK' a été détruit) : 266 : * 267 : * - Création/mise à jour de la table RV.'SPALART_ALLMARAS' contenant 268 : * les options de calcul personnalisables du modèle 269 : * 270 : * - Création des tables requises par les opérateurs TSCA et DFDT 271 : * 272 : * - Création des matrices pour les conditions de périodicité 273 : * 274 : * - Initialisation éventuelle de la variable interne du modèle 275 : * (viscosité modifiée NU0) 276 : * 277 : * - Création d'une sous-table dans 'INCO' qui contient les variables 278 : * intermédiaires du modèle si 'DUMP'=VRAI 279 : * 280 : 281 : 282 : SI (NON (EXIS RX 'INITOK')) ; 283 : 284 : 285 : * +===================================================================+ 286 : * | VALEURS PAR DÉFAUT DES PARAMÈTRES AVANCÉS PERSONNALISABLES | 287 : * +===================================================================+ 288 : 289 : * Version du modèle à utiliser 290 : * ---------------------------- 291 : KVERS = 'ORIG' ; 292 : 293 : 294 : * Nom de l'inconnue contenant la viscosité totale 295 : * ----------------------------------------------- 296 : _MUF = 'MUFN' ; 297 : 298 : 299 : * Constantes du modèle de Spalart-Allmaras 300 : * ---------------------------------------- 301 : KCONST = TABL ; 302 : KCONST . 'SIGMA' = 2./3. ; 303 : KCONST . 'CB1' = 0.1355 ; 304 : KCONST . 'CB2' = 0.622 ; 305 : KCONST . 'KAPPA' = 0.41 ; 306 : KCONST . 'CW1' = ((KCONST.'CB1')/((KCONST.'KAPPA')**2)) 307 : + ((1.+(KCONST.'CB2'))/(KCONST.'SIGMA')) ; 308 : KCONST . 'CW2' = 0.3 ; 309 : KCONST . 'CW3' = 2. ; 310 : KCONST . 'CV1' = 7.1 ; 311 : 312 : 313 : * Mesure scalaire du tenseur gradient des vitesses 314 : * ------------------------------------------------ 315 : KTGRAD = 'TOROT' ; 316 : 317 : 318 : * Instant auquel est renvoyé le résultat dans MUFN 319 : * ------------------------------------------------ 320 : KMUFN = 'APRES' ; 321 : 322 : 323 : * Algorithme de traitement des termes sources 324 : * ------------------------------------------- 325 : KSRC = 'ALGO1' ; 326 : 327 : 328 : * Configuration de l'algorithme de Newton 329 : * --------------------------------------- 330 : NEWTON = TABL ; 331 : NEWTON . 'CRIT' = 1.E-10 ; 332 : NEWTON . 'IMAX' = 10 ; 333 : NEWTON . 'OMEGA' = 1. ; 334 : 335 : 336 : * Options de la méthode d'inversion 337 : * --------------------------------- 338 : METINV = COPI (RV . 'METHINV') ; 339 : 340 : 341 : * Options de la discrétisation temporelle DFDT 342 : * -------------------------------------------- 343 : KOPT2 = COPI (RX . 'KOPT') ; 344 : KOPT2 . 'IDCEN' = 1 ; 345 : 346 : 347 : * Etat des différents verrous numériques 348 : * -------------------------------------- 349 : VERROU = TABL ; 350 : VERROU . 'POSITIF' = VRAI ; 351 : VERROU . 'DURBIN' = FAUX ; 352 : 353 : 354 : * Sauvegarder les variables internes? 355 : * ----------------------------------- 356 : KDUMP = FAUX ; 357 : 358 : 359 : 360 : * +===================================================================+ 361 : * | SURCHARGEMENT DES VALEURS PAR DÉFAUT PAR LES VALEURS UTILISATEUR | 362 : * +===================================================================+ 363 : * 364 : * Les valeurs par défaut définies ci-dessus sont placées dans une table 365 : * temporaire TABSA. Ensuite, la table RV.'SPALART_ALLMARAS' créée par 366 : * l'utilisateur est parcourue (si elle existe) et toutes ses valeurs 367 : * viennent surcharger celles de TABSA. Enfin, on met à jour la table 368 : * RV.'SPALART_ALLMARAS' en la remplaçant par la table TABSA. 369 : * 370 : 371 : TABSA = TABL ; 372 : TABSA . 'KVERS' = KVERS ; 373 : TABSA . 'NOMMUF' = _MUF ; 374 : TABSA . 'KCONST' = KCONST ; 375 : TABSA . 'KTGRAD' = KTGRAD ; 376 : TABSA . 'KMUFN' = KMUFN ; 377 : TABSA . 'KSRC' = KSRC ; 378 : TABSA . 'NEWTON' = NEWTON ; 379 : TABSA . 'METHINV' = METINV ; 380 : TABSA . 'KOPT2' = KOPT2 ; 381 : TABSA . 'VERROU' = VERROU ; 382 : TABSA . 'DUMP' = KDUMP ; 383 : 384 : SI (EXIS RV 'SPALART_ALLMARAS') ; 385 : RVSA = RV . 'SPALART_ALLMARAS' ; 386 : 387 : SI ((DIME RVSA) > 0) ; 388 : TIDX0 = INDE RVSA ; 389 : REPE BLOC1 (DIME TIDX0) ; 390 : MIDX1 = CHAI TIDX0 . &BLOC1 ; 391 : OBJ1 = RVSA . MIDX1 ; 392 : SI (EGA (TYPE OBJ1) 'TABLE') ; 393 : SI ((DIME OBJ1) > 0) ; 394 : TIDX1 = INDE OBJ1 ; 395 : REPE BLOC2 (DIME TIDX1) ; 396 : MIDX2 = CHAI TIDX1 . &BLOC2 ; 397 : TABSA . MIDX1 . MIDX2 = OBJ1 . MIDX2 ; 398 : FIN BLOC2 ; 399 : FINS ; 400 : SINON ; 401 : TABSA . MIDX1 = OBJ1 ; 402 : FINS ; 403 : FIN BLOC1 ; 404 : FINS ; 405 : 406 : FINS ; 407 : 408 : RV . 'SPALART_ALLMARAS' = TABSA ; 409 : 410 : 411 : 412 : 413 : 414 : * +===================================================================+ 415 : * | CRÉATION DES TABLES DE L'ÉQUATION DE TRANSPORT | 416 : * +===================================================================+ 417 : * 418 : * Création des tables de type 'KIZX' qui sont fournies à chaque pas de 419 : * temps aux opérateurs TSCA et DFDT. 420 : 421 : 422 : * Partie commune aux tables de TSCA et de DFDT 423 : RX0 = TABL 'KIZX' ; 424 : RX0 . 'EQEX' = RV ; 425 : RX0 . 'NOMZONE' = ' ' ; 426 : RX0 . 'DOMZ' = $MD ; 427 : RX0 . 'TDOMZ' = 0 ; 428 : RX0 . 'LISTINCO' = MOTS _NU0 ; 429 : 430 : 431 : * Table de l'opérateur TSCA 432 : RX1 = COPI RX0 ; 433 : RX1 . 'NOMOPER' = MOT 'TSCA_N' ; 434 : RX1 . 'KOPT' = (RX . 'KOPT') ; 435 : RX1 . 'IARG' = 3 ; 436 : RX1 . 'ARG1' = 'RHO' ; 437 : RX1 . 'ARG2' = 'UN' ; 438 : RX1 . 'ARG3' = 'SADIFF' ; 439 : 440 : RX . 'RX1' = RX1 ; 441 : 442 : 443 : * Table de l'opérateur DFDT (si résolution instationnaire) 444 : SI (RX . 'RINSTAT') ; 445 : RX2 = COPI RX0 ; 446 : RX2 . 'NOMOPER' = MOT 'DFDT_N' ; 447 : RX2 . 'KOPT' = (RV . 'SPALART_ALLMARAS' . 'KOPT2') ; 448 : RX2 . 'IARG' = 3 ; 449 : RX2 . 'ARG1' = 'RHO' ; 450 : RX2 . 'ARG2' = _NU0 ; 451 : RX2 . 'ARG3' = 'DT' ; 452 : 453 : * Gestion du décentrement pour DFDT (par défaut on l'a désactivé) 454 : IDCEN = (RX2 . 'KOPT' . 'IDCEN') ; 455 : SI ((EGA IDCEN 2) OU (EGA IDCEN 3)) ; 456 : RX2 . 'IARG' = 5 ; 457 : RX2 . 'ARG4' = 'UN' ; 458 : RX2 . 'ARG5' = _MUF ; 459 : FINS ; 460 : 461 : RX . 'RX2' = RX2 ; 462 : FINS ; 463 : 464 : 465 : 466 : 467 : 468 : 469 : 470 : * +===================================================================+ 471 : * | CRÉATION DES MATRICES POUR LES CONDITIONS PÉRIODIQUES | 472 : * +===================================================================+ 473 : 474 : SI (B_CYCL) ; 475 : MATC = RELA 'UX' GEOM1 - 'UX' GEOM2 ; 476 : STC = DEPI MATC 0. ; 477 : MATC = KOPS 'RIMA' MATC ; 478 : 479 : MATC = KOPS 'CHANINCO' MATC 480 : (MOTS 'LX' 'UX') (MOTS 'LX' _NU0) 481 : (MOTS 'FLX' 'FX') (MOTS 'LX' _NU0) ; 482 : STC = EXCO STC 'FLX' 'LX' ; 483 : 484 : RX . 'MATC' = MATC ; 485 : RX . 'STC' = STC ; 486 : FINS ; 487 : 488 : 489 : 490 : 491 : * +===================================================================+ 492 : * | CRÉATION DE L'INCONNUE INTERNE (VISCOSITÉ MODIFIÉE) | 493 : * +===================================================================+ 494 : 495 : SI (NON (EXIS (RV . 'INCO') _NU0)) ; 496 : NU0 = KCHT $MD 'SCAL' 'SOMMET' NU0_MIN ; 497 : RV . 'INCO' . _NU0 = NU0 ; 498 : FINS ; 499 : 500 : 501 : 502 : 503 : 504 : * +===================================================================+ 505 : * | CRÉATION DE LA TABLE CONTENANT LES VARIABLES INTERNES | 506 : * +===================================================================+ 507 : * Cette table ne sera utilisée que si 'DUMP'=VRAI 508 : 509 : SI (NON (EXIS (RV . 'INCO') 'SPAL')) ; 510 : RV . 'INCO' . 'SPAL' = TABL ; 511 : FINS ; 512 : 513 : 514 : 515 : 516 : * INITIALISATION TERMINÉE! 517 : RX . 'INITOK' = 1 ; 518 : FINS ; 519 : 520 : 521 : 522 : 523 : 524 : 525 : 526 : 527 : 528 : * +====================================================================+ 529 : * | | 530 : * | P R É P A R A T I O N D E S C A L C U L S | 531 : * | | 532 : * +====================================================================+ 533 : * 534 : 535 : 536 : * +===================================================================+ 537 : * | RÉCUPÉRATION DES PARAMÈTRES PERSONNALISABLES | 538 : * +===================================================================+ 539 : 540 : RVSA = RV . 'SPALART_ALLMARAS' ; 541 : 542 : 543 : * Version du modèle à utiliser 544 : * ---------------------------- 545 : KVERS = RVSA . 'KVERS' ; 546 : 547 : 548 : * Nom de l'inconnue contenant la viscosité totale 549 : * ----------------------------------------------- 550 : _MUF = RVSA . 'NOMMUF' ; 551 : 552 : 553 : * Constantes du modèle de Spalart-Allmaras 554 : * ---------------------------------------- 555 : SIGMA = RVSA . 'KCONST' . 'SIGMA' ; 556 : CB1 = RVSA . 'KCONST' . 'CB1' ; 557 : CB2 = RVSA . 'KCONST' . 'CB2' ; 558 : KAPPA = RVSA . 'KCONST' . 'KAPPA' ; 559 : CW1 = RVSA . 'KCONST' . 'CW1' ; 560 : CW2 = RVSA . 'KCONST' . 'CW2' ; 561 : CW3 = RVSA . 'KCONST' . 'CW3' ; 562 : CV1 = RVSA . 'KCONST' . 'CV1' ; 563 : 564 : 565 : * Mesure scalaire du tenseur gradient des vitesses 566 : * ------------------------------------------------ 567 : KTGRAD = RVSA . 'KTGRAD' ; 568 : 569 : 570 : * Instant auquel est renvoyé le résultat dans MUFN 571 : * ------------------------------------------------ 572 : KMUFN = RVSA . 'KMUFN' ; 573 : 574 : * Incompatibilité, car les conditions périodiques ne sont vérifiées 575 : * qu'après l'équation de transport, soit KMUFN='AVANT' ou KMUFN='APRES' 576 : SI (B_CYCL ET (EGA KMUFN 'DEMI')) ; 577 : MESS '******************************************************' ; 578 : MESS '/!\ ERREUR dans SPAL :' ; 579 : MESS ' options KMUFN="DEMI" et "PERIODIC" incompatibles' ; 580 : MESS '******************************************************' ; 581 : QUIT SPAL ; 582 : FINS ; 583 : 584 : 585 : * Algorithme de traitement des termes sources 586 : * ------------------------------------------- 587 : KSRC = RVSA . 'KSRC' ; 588 : 589 : * Incompatibilité, car l'algorithme 1 est nécessairement instationnaire 590 : SI ((EGA KSRC 'ALGO1') ET (NON (RX . 'RINSTAT'))) ; 591 : MESS '******************************************************' ; 592 : MESS '/!\ ERREUR dans SPAL :' ; 593 : MESS ' KSRC="ALGO1" requiert la présence de DFDT' ; 594 : MESS '******************************************************' ; 595 : QUIT SPAL ; 596 : FINS ; 597 : 598 : * Incompatibilité, car il n'y a pas de demi pas de temps avec 'ALGO2' 599 : SI ((EGA KSRC 'ALGO2') ET (EGA KMUFN 'DEMI')) ; 600 : MESS '******************************************************' ; 601 : MESS '/!\ ERREUR dans SPAL :' ; 602 : MESS ' options KSRC="ALGO2" et KMUFN="DEMI" incompatibles' ; 603 : MESS '******************************************************' ; 604 : QUIT SPAL ; 605 : FINS ; 606 : 607 : 608 : * Configuration de l'algorithme de Newton 609 : * --------------------------------------- 610 : EPS1 = RVSA . 'NEWTON' . 'CRIT' ; 611 : IMAX = RVSA . 'NEWTON' . 'IMAX' ; 612 : OMEGA = RVSA . 'NEWTON' . 'OMEGA' ; 613 : 614 : 615 : * Options de la méthode d'inversion 616 : * --------------------------------- 617 : METINV = RVSA . 'METHINV' ; 618 : 619 : 620 : * Etat des différents verrous numériques 621 : * -------------------------------------- 622 : B_POSI = RVSA . 'VERROU' . 'POSITIF' ; 623 : *B_DURB = RVSA . 'VERROU' . 'DURBIN' ; 624 : 625 : 626 : * Sauvegarder les variables internes? 627 : * ----------------------------------- 628 : KDUMP = RVSA . 'DUMP' ; 629 : 630 : 631 : 632 : 633 : * +===================================================================+ 634 : * | RÉCUPÉRATION DES VARIABLES AU PAS DE TEMPS COURANT | 635 : * +===================================================================+ 636 : 637 : * Viscosité modifiée au début du pas de temps 638 : NU0 = RV . 'INCO' . _NU0 ; 639 : 640 : 641 : 642 : * Mesure scalaire du tenseur gradient des vitesses 643 : S = (PRODT $MD KTGRAD UN) ** 0.5 ; 644 : 645 : 646 : 647 : * Distance à la paroi (peut changer à chaque pas de temps) 648 : DPAROI = RV . 'PAROIS' . 'DIST' ; 649 : KPAR = DPAROI MASQ 'EGAL' 0. ; 650 : 651 : * /!\ La distance à la paroi intervient uniquement au dénominateur, et 652 : * doit donc être strictement positive : 653 : DPAROI = DPAROI + (KPAR * 1.E-10) ; 654 : 655 : * /!\ USLM2, inverse du carré de DPAROI, est donc singulier à la paroi. 656 : * Or ce terme apparait toujours multiplié par NU0, qui vaut 0 à la 657 : * paroi (condition limite recommandée!). 658 : * => C'est ainsi le produit (USLM2*NU0) vaut bien 0 à la paroi. 659 : * 660 : * Mais dans la pratique: 661 : * - un limiteur peut empêcher les valeurs strict. nulles dans NU0 662 : * - la condition limite à la paroi peut être définie non nulle 663 : * Pour éviter que le produit (USLM2*NU0) ne fasse apparaître de 664 : * très grandes valeurs, on force donc USLM2=0 sur les parois. 665 : USLM2 = (1. - KPAR) * (INVE ((KAPPA*DPAROI)**2)) ; 666 : 667 : 668 : 669 : 670 : 671 : * +===================================================================+ 672 : * | RÉCUPÉRATION DES CONDITIONS AUX LIMITES | 673 : * +===================================================================+ 674 : 675 : * /!\ Les conditions aux limites portent sur la viscosité modifiée ! 676 : * ****************** 677 : SI (EXIS (RV . 'CLIM') _NU0) ; 678 : CLVAL = EXCO (RV . 'CLIM') _NU0 ; 679 : RV . 'SPALART_ALLMARAS' . 'CLIM' = COPI CLVAL ; 680 : RV . 'CLIM' = ENLE (RV . 'CLIM') _NU0 ; 681 : SINON ; 682 : SI (EXIS (RV . 'SPALART_ALLMARAS') 'CLIM') ; 683 : CLVAL = RV . 'SPALART_ALLMARAS' . 'CLIM' ; 684 : SINON ; 685 : SI (NON B_CYCL) ; 686 : MESS '******************************************************' ; 687 : MESS '/!\ ERREUR dans SPAL :' ; 688 : MESS ' Il faut définir les conditions limites pour ' _NU0 ; 689 : MESS '******************************************************' ; 690 : QUIT SPAL ; 691 : FINS ; 692 : FINS ; 693 : FINS ; 694 : CLSPG = EXTR CLVAL 'MAILLAGE' ; 695 : 696 : * Méthode de pénalisation pour imposer les conditions limites 697 : SI (EGA KSRC 'ALGO1') ; 698 : PENAVAL = 1.E30 ; 699 : PENAMAT = DG + (REDU (PENAVAL - DG) CLSPG) ; 700 : PENAST = (CLVAL * PENAVAL) ; 701 : FINS ; 702 : 703 : 704 : 705 : 706 : 707 : 708 : 709 : 710 : 711 : 712 : * +====================================================================+ 713 : * | | 714 : * | T R A I T E M E N T D E S T E R M E S S O U R C E S | 715 : * | | 716 : * +====================================================================+ 717 : * 718 : * Si 'KSRC'='ALGO1', on procède à l'ensemble des itérations internes 719 : * Si 'KSRC'='ALGO2', on quitte le bloc avant la fin du premier passage 720 : * 721 : 722 : * Initialisation de la variable interne NUD1 723 : * ========================================== 724 : 725 : NUD1 = NU0 ; 726 : * -------------------- /!\ VERROU NUMÉRIQUE /!\ -------------------- 727 : SI (B_POSI) ; 728 : NUD1 = KOPS NUD1 '|<' NU0_MIN ; 729 : FINS ; 730 : 731 : SI (EGA KMUFN 'AVANT') ; 732 : NU1 = NUD1 ; 733 : FINS ; 734 : 735 : 736 : 737 : * Itérations internes pour résoudre F1(NUD1)=0 738 : * ============================================ 739 : 740 : * Dénominateur de la norme infinie utilisée pour le critère de sortie 741 : DENERR = (MAXI NU0 'ABS') + 1.E-30 ; 742 : 743 : REPE BLOC0 IMAX ; 744 : 745 : * Calcul du terme source complet pour l'itération interne courante 746 : * ---------------------------------------------------------------- 747 : 748 : KSI = NUD1 * USNU ; 749 : FV1 = (KSI**3) * ( INVE ((KSI**3) + (CV1**3)) ) ; 750 : DENFV2 = INVE (1. + (FV1*KSI)) ; 751 : FV2 = 1. - (KSI*DENFV2) ; 752 : FV3 = 1. ; 753 : S0 = (S*FV3) + (FV2*NUD1*USLM2) ; 754 : * ------------------- /!\ VERROU NUMÉRIQUE /!\ ------------------- 755 : SI (B_POSI) ; 756 : S0 = KOPS S0 '|<' S0_MIN ; 757 : FINS ; 758 : USS0 = INVE S0 ; 759 : 760 : R = NUD1 * USS0 * USLM2 ; 761 : * ------------------- /!\ VERROU NUMÉRIQUE /!\ ------------------- 762 : R = KOPS R '>|' 10. ; 763 : G = R + (CW2 * ((R**6) - R)) ; 764 : DENFW6 = INVE ((G**6) + (CW3**6)) ; 765 : GLIM = ((1.+(CW3**6))*DENFW6) ** (1./6.) ; 766 : FW = G * GLIM ; 767 : 768 : * (B contient la somme des 3 termes non linéaires du second membre) 769 : GDNUD1 = KOPS NU0 'GRADS' $MD ; 770 : * GDNUD1 = KOPS NUD1 'GRADS' $MD ; 771 : LCO1 = EXTR GDNUD1 'COMP' ; 772 : B1 = (CB2 / SIGMA) * (PSCA GDNUD1 GDNUD1 LCO1 LCO1) ; 773 : B2 = CB1 * S0 * NUD1 ; 774 : B3 = CW1 * FW * (NUD1*NUD1) * USLM2*(KAPPA*KAPPA) ; 775 : B = B1 + B2 - B3 ; 776 : 777 : 778 : 779 : * ALGO2 780 : *********************************************************************** 781 : *********************************************************************** 782 : * On n'a pas besoin d'aller plus loin pour 'ALGO2': on sépare juste 783 : * les parties positive et négative de B avant de sortir 784 : SI (EGA KSRC 'ALGO2') ; 785 : 786 : * REMARQUE: Pour le modèle de base, si B_POSI=VRAI et si CW2<1, alors 787 : * PROD est du signe de CB1(>0) et DEST de celui de CW1(>0) 788 : * => plus besoin de séparer S+ et S-, c'est déjà fait! 789 : * 790 : SI ((EGA KVERS 'ORIG') ET B_POSI ET (CW2 < 1)) ; 791 : BPLUS = B1 + B2 ; 792 : BMOINS = B3 ; 793 : SINON ; 794 : BPLUS = KOPS (B1 + B2) '|<' 0. ; 795 : BMOINS = KOPS B3 '|<' 0. ; 796 : FINS ; 797 : 798 : * /!\ SEMBLE TRÈS INSTABLE => NON RECOMMANDÉ 799 : * MASK1 = B MASQ 'SUPERIEUR' 0. ; 800 : * BPLUS = B * MASK1 ; 801 : * BMOINS = B * (1. - MASK1) ; 802 : 803 : QUIT BLOC0 ; 804 : 805 : FINS ; 806 : *********************************************************************** 807 : *********************************************************************** 808 : 809 : 810 : 811 : * Calcul des dérivées exactes des variables précédentes 812 : * ----------------------------------------------------- 813 : 814 : FV2PRIME = ((3.*KSI*FV1*(1.-FV1)) - 1.) * (DENFV2*DENFV2) ; 815 : S0PRIME = (FV2 + (KSI*FV2PRIME)) * USLM2 ; 816 : RPRIME = (S0 - (NUD1*S0PRIME)) * (USS0*USS0) * USLM2 ; 817 : GPRIME = RPRIME * (1. + (CW2 * ((6.*(R**5)) - 1.))) ; 818 : FWPRIME = GPRIME * GLIM * (CW3**6) * DENFW6 ; 819 : 820 : B2PRIME = CB1 * (S0 + (NUD1*S0PRIME)) ; 821 : B3PRIME = CW1 * NUD1 * USLM2*(KAPPA*KAPPA) 822 : * ((2.*FW) + (NUD1*FWPRIME)) ; 823 : BPRIME = B2PRIME - B3PRIME ; 824 : 825 : 826 : 827 : 828 : * Calcul de la nouvelle valeur de NUD1 829 : * ------------------------------------ 830 : * (F1 est la fonction de NUD1 à annuler: elle prend en compte la 831 : * discrétisation temporelle et les conditions limites) 832 : 833 : F1 = (PENAMAT*NUD1) - PENAST - (DG * (NU0 + (DT*B))) ; 834 : F1PRIME = PENAMAT - (DG * DT * BPRIME) ; 835 : * ------------------- /!\ VERROU NUMÉRIQUE /!\ ------------------- 836 : F1PRIME = F1PRIME + ((F1PRIME MASQ 'EGAL' 0.) * 1.E-30) ; 837 : 838 : NUD2 = NUD1 - (F1 * (INVE F1PRIME)) ; 839 : * -------------------- /!\ VERROU NUMÉRIQUE /!\ ------------------ 840 : SI (B_POSI) ; 841 : NUD2 = KOPS NUD2 '|<' NU0_MIN ; 842 : FINS ; 843 : 844 : 845 : 846 : 847 : * Critère d'arrêt: fin des itérations internes 848 : * -------------------------------------------- 849 : XERR = (MAXI (NUD2-NUD1) 'ABS') / DENERR ; 850 : 851 : * Mise à jour de la variable interne (avec relaxation éventuelle) 852 : NUD1 = (OMEGA * NUD2) + ((1.-OMEGA) * NUD1) ; 853 : 854 : 855 : SI ((XERR < EPS1) OU (EGA &BLOC0 IMAX)) ; 856 : MESS (CHAI 'SPAL [sources]: iter ' &BLOC0 ' / ' IMAX 857 : ' :: norm inf = ' XERR) ; 858 : QUIT BLOC0 ; 859 : FINS ; 860 : 861 : 862 : FIN BLOC0 ; 863 : 864 : 865 : * Mise à jour de la table 'INCO' pour l'étape suivante 866 : * ==================================================== 867 : 868 : RV . 'INCO' . 'SADIFF' = ((NU + NUD1) * (1./SIGMA)) ; 869 : RV . 'INCO' . _NU0 = NUD1 ; 870 : SI (EGA KMUFN 'DEMI') ; 871 : NU1 = NUD1 ; 872 : FINS ; 873 : 874 : 875 : 876 : 877 : 878 : 879 : * +====================================================================+ 880 : * | | 881 : * | É Q U A T I O N D E T R A N S P O R T | 882 : * | | 883 : * +====================================================================+ 884 : * 885 : 886 : * Calcul des matrices de l'équation de transport 887 : * ============================================== 888 : * /!\ Tout changement dans les options de discrétisation 'KOPT' ne sera 889 : * effectif qu'après destruction de l'indice RX.'INITOK' 890 : 891 : ST1 MAT1 = TSCA (RX . 'RX1') ; 892 : SI (RX . 'RINSTAT') ; 893 : ST2 MAT2 = DFDT (RX . 'RX2') ; 894 : MAT1 = MAT1 ET MAT2 ; 895 : ST1 = ST1 ET ST2 ; 896 : FINS ; 897 : 898 : * ALGO2 899 : *********************************************************************** 900 : *********************************************************************** 901 : SI (EGA KSRC 'ALGO2') ; 902 : ST2 = NOMC (BPLUS*DG) _NU0 'NATURE' 'DISCRET' ; 903 : BDIA = BMOINS * DG * (INVE NUD1) ; 904 : MAT2 = KOPS 'MATDIAGO' (NOMC _NU0 BDIA) 'MATRIK' ; 905 : 906 : MAT1 = MAT1 ET MAT2 ; 907 : ST1 = ST1 ET ST2 ; 908 : FINS ; 909 : *********************************************************************** 910 : *********************************************************************** 911 : 912 : 913 : 914 : * Récupération des matrices de périodicité 915 : * ======================================== 916 : SI (B_CYCL) ; 917 : MAT1 = MAT1 ET (RX . 'MATC') ; 918 : ST1 = ST1 ET (RX . 'STC') ; 919 : FINS ; 920 : 921 : 922 : 923 : * Initialisation de l'algorithme d'inversion 924 : * ========================================== 925 : METINV . 'MATASS' = MAT1 ; 926 : METINV . 'MAPREC' = MAT1 ; 927 : METINV . 'XINIT' = NOMC NUD1 _NU0 ; 928 : METINV . 'IMPINV' = 0 ; 929 : 930 : 931 : 932 : * Résolution de l'équation de transport 933 : * ===================================== 934 : NUD3 = KRES MAT1 'SMBR' ST1 935 : 'CLIM' (NOMC CLVAL _NU0) 936 : 'TYPI' METINV 937 : 'IMPR' 0 ; 938 : 939 : SI (EXIS METINV 'CONVINV') ; 940 : LCONV1 = METINV . 'CONVINV' ; 941 : NCONV1 = DIME LCONV1 ; 942 : MESS (CHAI 'SPAL [transport]: iter ' (NCONV1 / 2) ' / ' 943 : ((METINV.'NITMAX') / 2) 944 : ' :: norm inf = ' (EXTR LCONV1 NCONV1)) ; 945 : FINS ; 946 : 947 : NUD3 = EXCO NUD3 _NU0 ; 948 : * -------------------- /!\ VERROU NUMÉRIQUE /!\ ------------------ 949 : SI (B_POSI) ; 950 : NUD3 = KOPS NUD3 '|<' NU0_MIN ; 951 : FINS ; 952 : 953 : 954 : 955 : * Mise à jour de la table 'INCO' 956 : * ============================== 957 : RV . 'INCO' . _NU0 = NUD3 ; 958 : SI (EGA KMUFN 'APRES') ; 959 : NU1 = NUD3 ; 960 : FINS ; 961 : 962 : * Sauvegarde des variables internes 963 : SI (KDUMP) ; 964 : RV . 'INCO' . 'SPAL' . 'KSI' = KSI ; 965 : RV . 'INCO' . 'SPAL' . 'FV1' = FV1 ; 966 : RV . 'INCO' . 'SPAL' . 'FV2' = FV2 ; 967 : RV . 'INCO' . 'SPAL' . 'FV3' = FV3 ; 968 : RV . 'INCO' . 'SPAL' . 'S' = S ; 969 : RV . 'INCO' . 'SPAL' . 'S0' = S0 ; 970 : RV . 'INCO' . 'SPAL' . 'R' = R ; 971 : RV . 'INCO' . 'SPAL' . 'G' = G ; 972 : RV . 'INCO' . 'SPAL' . 'FW' = FW ; 973 : RV . 'INCO' . 'SPAL' . 'B1' = B1 ; 974 : RV . 'INCO' . 'SPAL' . 'B2' = B2 ; 975 : RV . 'INCO' . 'SPAL' . 'B3' = B3 ; 976 : SI (EGA KSRC 'ALGO1') ; 977 : RV . 'INCO' . 'SPAL' . 'FV2PRIME' = FV2PRIME ; 978 : RV . 'INCO' . 'SPAL' . 'S0PRIME' = S0PRIME ; 979 : RV . 'INCO' . 'SPAL' . 'RPRIME' = RPRIME ; 980 : RV . 'INCO' . 'SPAL' . 'GPRIME' = GPRIME ; 981 : RV . 'INCO' . 'SPAL' . 'FWPRIME' = FWPRIME ; 982 : RV . 'INCO' . 'SPAL' . 'B2PRIME' = B2PRIME ; 983 : RV . 'INCO' . 'SPAL' . 'B3PRIME' = B3PRIME ; 984 : RV . 'INCO' . 'SPAL' . 'F1' = F1 ; 985 : RV . 'INCO' . 'SPAL' . 'F1PRIME' = F1PRIME ; 986 : FINS ; 987 : FINS ; 988 : 989 : * Calcul de la viscosité turbulente effective à l'instant défini par 990 : * le paramètre KMUFN puis mise à jour de l'indice défini par le 991 : * paramètre NOMMUF. 992 : KSI = NU1 * USNU ; 993 : FV1 = (KSI**3) * ( INVE ((KSI**3) + (CV1**3)) ) ; 994 : MUF = RHO * (NU + (FV1*NU1)) ; 995 : MUF = KCHT $MD 'SCAL' 'SOMMET' MUF ; 996 : RV . 'INCO' . _MUF = MUF ; 997 : 998 : 999 : 1000 : 1001 : * Fin de la procédure 1002 : ST0 MAT0 = KOPS 'MATRIK' ; 1003 : FINP ST0 MAT0 ; 1004 : 1005 : 1006 :
© Cast3M 2003 - All rights reserved.
Disclaimer