Télécharger SPAL.procedur

Retour à la liste

Numérotation des lignes :

  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.  
  1007.  

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