Télécharger prodt.procedur

Retour à la liste

Numérotation des lignes :

  1. * PRODT PROCEDUR JC220346 10/11/18 21:15:01 6795
  2. ************************************************************************
  3. * NOM : PRODT
  4. * DESCRIPTION : Renvoie une mesure du tenseur gradient des vitesses.
  5. * Cette grandeur intervient notamment dans l'expression
  6. * du terme de production turbulente.
  7. ************************************************************************
  8. * HISTORIQUE : 14/03/2000 : MAGN : création de la procédure
  9. * HISTORIQUE : 28/10/2010 : JCARDO : généralisation de PRODT
  10. * 1) ajout de 3 autres mesures pour le tenseur gradient,
  11. * en plus du taux de déformation: - taux de rotation
  12. * - tenseur complet
  13. * - expression mixte
  14. * 2) prise en compte du cas où NCO est différent de NDIM
  15. * 3) gestion de composantes de noms quelconques dans UN
  16. * 4) modification de la syntaxe et de la notice
  17. * HISTORIQUE :
  18. * HISTORIQUE :
  19. ************************************************************************
  20. * Prière de PRENDRE LE TEMPS DE COMPLÉTER LES COMMENTAIRES
  21. * en cas de modification de ce sous-programme afin de faciliter
  22. * la maintenance !
  23. ************************************************************************
  24. * Syntaxe :
  25. * _________
  26. *
  27. * CHPO1 = PRODT $MD (| 'TODEF' |) UN (GB TN) ;
  28. * | 'TOROT' |
  29. * | 'COMPL' |
  30. * | 'MIXTE' |
  31. *
  32. *
  33. ************************************************************************
  34.  
  35.  
  36. DEBP PRODT ;
  37.  
  38.  
  39. * Dimension de l'espace physique
  40. NDIM = VALE 'DIME' ;
  41.  
  42. SI (EGA NDIM 1) ;
  43. MESS '******************************************************' ;
  44. MESS ' /!\ ERREUR dans PRODT /!\ ' ;
  45. MESS ' ' ;
  46. MESS ' Cette procédure est inutile pour les problèmes 1D.' ;
  47. MESS '******************************************************' ;
  48. QUIT PRODT ;
  49. FINS ;
  50.  
  51.  
  52.  
  53. * +====================================================================+
  54. * | |
  55. * | L E C T U R E D E S A R G U M E N T S |
  56. * | |
  57. * +====================================================================+
  58.  
  59. * +========+
  60. * | MODÈLE |
  61. * +========+
  62.  
  63. ARGU $MD*'MMODEL' ;
  64.  
  65.  
  66. * +=========+
  67. * | MOT-CLÉ |
  68. * +=========+
  69.  
  70. ARGU MCLE/'MOT' ;
  71. SI (NON (EXIS MCLE)) ;
  72. MCLE = 'TODEF' ;
  73. FINS ;
  74.  
  75. B_TODEF = (EGA MCLE 'TODEF') ;
  76. B_TOROT = (EGA MCLE 'TOROT') ;
  77. B_COMPL = (EGA MCLE 'COMPL') ;
  78. B_MIXTE = (EGA MCLE 'MIXTE') ;
  79.  
  80.  
  81. * +==================+
  82. * | CHAMP DE VITESSE |
  83. * +==================+
  84.  
  85. ARGU UN*'CHPOINT' ;
  86. LCO = EXTR UN 'COMP' ;
  87. NCO = DIME LCO ;
  88.  
  89.  
  90. * +==============================================+
  91. * | CHAMP DE TEMPÉRATURE ET VECTEUR FLOTTABILITÉ |
  92. * +==============================================+
  93.  
  94. ARGU TN/'CHPOINT' ;
  95.  
  96. B_TN = FAUX ;
  97. SI (EXIS TN) ;
  98.  
  99. ARGU GB/'POINT' ;
  100. SI (NON (EXIS GB)) ;
  101. MESS '******************************************************' ;
  102. MESS '/!\ ERREUR dans PRODT :' ;
  103. MESS ' Champ de température fourni => on attendait GB' ;
  104. MESS '******************************************************' ;
  105. ERRE 5 ;
  106. QUIT PRODT ;
  107. FINS ;
  108.  
  109. B_TN = VRAI ;
  110. FINS ;
  111.  
  112.  
  113.  
  114.  
  115. * +====================================================================+
  116. * | |
  117. * | C A L C U L D E S R É S U L T A T S |
  118. * | |
  119. * +====================================================================+
  120.  
  121. * Extraction des composantes du tenseur gradient
  122. * ==============================================
  123.  
  124. UNX = EXCO UN (EXTR LCO 1) ;
  125. GDUX = KOPS 'GRADS' UNX $MD ;
  126. DUDX = EXCO GDUX 'UX' ;
  127. DUDY = EXCO GDUX 'UY' ;
  128.  
  129. SI ((NCO > 1) ET (NDIM > 1)) ;
  130. UNY = EXCO UN (EXTR LCO 2) ;
  131. GDUY = KOPS 'GRADS' UNY $MD ;
  132. DVDX = EXCO GDUY 'UX' ;
  133. DVDY = EXCO GDUY 'UY' ;
  134.  
  135. SI ((NCO > 2) ET (NDIM EGA 3)) ;
  136. UNZ = EXCO UN (EXTR LCO 3) ;
  137. GDUZ = KOPS 'GRADS' UNZ $MD ;
  138. DWDX = EXCO GDUZ 'UX' ;
  139. DWDY = EXCO GDUZ 'UY' ;
  140. DWDZ = EXCO GDUZ 'UZ' ;
  141. SINON ;
  142. DWDX = 0. ;
  143. DWDY = 0. ;
  144. DWDZ = 0. ;
  145. FINS ;
  146.  
  147. SI (NDIM EGA 3) ;
  148. DUDZ = EXCO GDUX 'UZ' ;
  149. DVDZ = EXCO GDUY 'UZ' ;
  150. SINON ;
  151. DUDZ = 0. ;
  152. DVDZ = 0. ;
  153. FINS ;
  154.  
  155. SINON ;
  156. SI (NDIM EGA 3) ;
  157. DUDZ = EXCO GDUX 'UZ' ;
  158. SINON ;
  159. DUDZ = 0. ;
  160. FINS ;
  161. DVDX = 0. ;
  162. DVDY = 0. ;
  163. DVDZ = 0. ;
  164. DWDX = 0. ;
  165. DWDY = 0. ;
  166. DWDZ = 0. ;
  167. FINS ;
  168.  
  169.  
  170.  
  171. * Calcul du carré des normes des parties symétrique et anti-symétrique
  172. * ====================================================================
  173.  
  174. A2 = (DUDY*DUDY) + (DUDZ*DUDZ) + (DVDX*DVDX) + (DVDZ*DVDZ)
  175. + (DWDX*DWDX) + (DWDY*DWDY) ;
  176. A3 = (DUDY*DVDX) + (DUDZ*DWDX) + (DVDZ*DWDY) ;
  177. SI (NON B_TOROT) ;
  178. A1 = (DUDX*DUDX) + (DVDY*DVDY) + (DWDZ*DWDZ) ;
  179. FINS ;
  180.  
  181.  
  182. * Mot-clé 'TODEF': taux de déformation => 2.SijSij (par défaut)
  183. SI (B_TODEF) ;
  184. P = (2. * (A1 + A3)) + A2 ;
  185. FINS ;
  186.  
  187. * Mot-clé 'TOROT': taux de rotation => 2.RijRij
  188. SI (B_TOROT) ;
  189. * /!\ Valeur absolue pour éviter les zéros "négatifs"
  190. P = ABS ( A2 + (-2. * A3) ) ;
  191. FINS ;
  192.  
  193. * Mot-clé 'COMPL': tenseur gradient complet => UijUij
  194. SI (B_COMPL) ;
  195. P = A1 + A2 ;
  196. FINS ;
  197.  
  198. * Mot-clé 'MIXTE': Expression proposée par Daclès-Mariani et al.
  199. SI (B_MIXTE) ;
  200. CPROD = 2. ;
  201. DRR = (ABS (A2 + (-2. * A3))) ** 0.5 ;
  202. DSS = ((2. * (A1 + A3)) + A2) ** 0.5 ;
  203. DMX = KOPS (DSS - DRR) '>|' 0. ;
  204. P = (DRR + (CPROD * DMX)) ** 2 ;
  205. FINS ;
  206.  
  207.  
  208.  
  209.  
  210. * Prise en compte de la contribution des forces de flottabilité
  211. * =============================================================
  212. * (Production turbulente moindre si stratification stable)
  213.  
  214. SI (B_TN) ;
  215. * Prandtl turbulent
  216. SGT = 0.7 ;
  217.  
  218. * Calcul du produit scalaire gB.Grad(T)
  219. GDT = KOPS 'GRADS' TN $MD ;
  220. DTDX = EXCO GDT 'UX' ;
  221. DTDY = EXCO GDT 'UY' ;
  222. AX = COOR 1 GB ;
  223. AY = COOR 2 GB ;
  224. G = (DTDX*AX) + (DTDY*AY) ;
  225. SI (EGA NDIM 3) ;
  226. DTDZ = EXCO GDT 'UZ' ;
  227. AZ = COOR 3 GB ;
  228. G = G + (DTDZ*AZ) ;
  229. FINS ;
  230.  
  231. * G > 0 => dT/dz < 0 => Stratification instable
  232. * G < 0 => dT/dz > 0 => Stratification stable
  233. IG = MASQ G 'INFERIEUR' 0. ;
  234. G = G * IG * (1./SGT) ;
  235.  
  236. * Modélisation selon W.Rodi
  237. C3 = 0.8 ;
  238. PG = P + G ;
  239. IPG = MASQ (ABS PG) 'INFERIEUR' 1.E-30 ;
  240. RF = (-1.) * G * (INVE (PG + (IPG*1.E-30))) ;
  241. P = PG * (1. + (C3 * RF)) ;
  242. FINS ;
  243.  
  244.  
  245. * Empêche une éventuelle division par zéro
  246. P = P + 1.E-30 ;
  247.  
  248. RESP P ;
  249. FINP ;
  250.  
  251.  
  252.  
  253.  

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