Télécharger tlap12.eso

Retour à la liste

Numérotation des lignes :

tlap12
  1. C TLAP12 SOURCE CB215821 20/11/25 13:40:58 10792
  2. SUBROUTINE TLAP12(NESP,IMUC,IMTC,TSIGK,TSIGE,IGRKEP,
  3. & MELEMC,MELEMF,MELEFL,ISURF,INORM,IDIAM,
  4. & ICHFLU,DT)
  5. C***********************************************************************
  6. C NOM : TLAP12
  7. C DESCRIPTION : Calcul des flux diffusifs pour k-\eps equations de
  8. C turbulence (2D)
  9. C
  10. C
  11. C LANGAGE : ESOPE
  12. C AUTEUR : S. Kudriakov (CEA/DEN/DM2S/SFME/LTMF)
  13. C mél : skudriakov@cea.fr
  14. C***********************************************************************
  15. C APPELES (UTIL) : KRIPAD : MELEME -> (num. globale->locale)
  16. C LICHT : Lecture des pointeurs (maillages, valeurs)
  17. C d'un objet de type MCHPOI.
  18. C ERREUR : Gestion des erreurs par GIBI.
  19. C APPELE PAR : ZLAP11 : Chapeau de l'opérateur Gibiane 'LAPN'
  20. C option 'VF'.
  21. C***********************************************************************
  22. C ENTREES : NESP (type ENTIER) : number of species explicitly
  23. C treated in the NS equations
  24. C IMUC (type MCHPOI) : laminar viscosity
  25. C IMTC (type MCHPOI) : turbulent viscosity
  26. C TSIGK (type REEL) : turbulent constant \sigma_k
  27. C TSIGE (type REEL) : turbulent constant
  28. C \sigma_{\eps}
  29. C IGRKEP (type MCHPOI) : gradients of k and \epsilon
  30. C
  31. C MELEMC (type MELEME) : maillage des centres des
  32. C éléments.
  33. C MELEMF (type MELEME) : maillage des faces des
  34. C éléments.
  35. C MELEFL (type MELEME) : connectivités face-(centre
  36. C gauche, centre droit).
  37. C ISURF (type MCHPOI) : surface des faces.
  38. C INORM (type MCHPOI) : normale aux faces.
  39. C IDIAM (type MCHPOI) : diamètre des éléments.
  40. C ENTREES/SORTIES : -
  41. C SORTIES : ICHFLU (type MCHPOI) : flux diffusif aux
  42. C interfaces.
  43. C DT (type REAL*8) : pas de temps de stabilité
  44. C (Fourier)
  45. C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  46. C***********************************************************************
  47. C VERSION : v1, 3/12/2003, version initiale
  48. C HISTORIQUE :
  49. C***********************************************************************
  50. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  51. C en cas de modification de ce sous-programme afin de faciliter
  52. C la maintenance !
  53. C***********************************************************************
  54. IMPLICIT INTEGER(I-N)
  55.  
  56. -INC PPARAM
  57. -INC CCOPTIO
  58. -INC CCREEL
  59. -INC SMCHPOI
  60. -INC SMELEME
  61. -INC SMCOORD
  62. -INC SMLENTI
  63. -INC SMLMOTS
  64. C
  65. POINTEUR IMUC.MCHPOI,IMTC.MCHPOI,IGRKEP.MCHPOI
  66. POINTEUR MPMUC.MPOVAL,MPMTC.MPOVAL,MPGRKE.MPOVAL,
  67. & MPSURF.MPOVAL, MPNORM.MPOVAL, MPDIAM.MPOVAL,
  68. & MPFLUX.MPOVAL
  69. C
  70. POINTEUR MELEMC.MELEME,MELEMF.MELEME,MELEFL.MELEME
  71. POINTEUR MLCENT.MLENTI,MLETIM.MLENTI,MLERIM.MLENTI
  72. C
  73. INTEGER ITIMP
  74. & ,ISURF,INORM,IDIAM,ICHFLU
  75. & ,NFAC, NLCF, NGCF, NGCF1, NGCEG, NGCED
  76. & ,NLCEG,NLCED,NLFTI,NLFRI
  77. & , ICOORX, IGEOM
  78.  
  79. REAL*8 DT, UNSDT
  80. & ,XG,YG,XFMXG,YFMYG,DRG
  81. & ,XD,YD,XFMXD,YFMYD,DRD,ALPHA,UMALPH
  82. & ,XF,YF
  83. & ,CNX, CNY, ORIENT, DIAM, DIAM2, CELL, SURF
  84. INTEGER NESP
  85. REAL*8 DKADX,DKADY,DEPDX,DEPDY
  86. REAL*8 MU,MUT,TSIGK,TSIGE
  87. C
  88. CHARACTER*8 TYPE
  89. C
  90. C**** Initialisation de 1/DT
  91. C
  92. UNSDT = 0.0D0
  93. C
  94. C**** KRIPAD pour la correspondance global/local de centre
  95. C
  96. CALL KRIPAD(MELEMC,MLCENT)
  97. C
  98. C EN KRIPAD
  99. C SEGACT MELEMC
  100. C SEGACT MLCENT
  101. C
  102. CALL LICHT(IMUC,MPMUC,TYPE,IGEOM)
  103. CALL LICHT(IMTC,MPMTC,TYPE,IGEOM)
  104. CALL LICHT(IGRKEP,MPGRKE,TYPE,IGEOM)
  105. CALL LICHT(ISURF,MPSURF,TYPE,IGEOM)
  106. CALL LICHT(INORM,MPNORM,TYPE,IGEOM)
  107. CALL LICHT(IDIAM,MPDIAM,TYPE,IGEOM)
  108. CALL LICHT(ICHFLU,MPFLUX,TYPE,IGEOM)
  109. C
  110. C EN LICHT
  111. C SEGACT*MOD MPROC
  112. C SEGACT*MOD MPTEMC
  113. C SEGACT*MOD MPSURF
  114. C SEGACT*MOD MPNORM
  115. C SEGACT*MOD MPDIAM
  116. C SEGACT*MOD MPFLUX
  117. C
  118. C---------------------------------------------------------
  119. SEGACT MELEFL
  120. SEGACT MELEMF
  121. NFAC = MELEMF.NUM(/2)
  122. C
  123. C**** Boucle sur les faces
  124. C
  125. DO NLCF = 1, NFAC, 1
  126. C
  127. C******* NLCF = numero local du centre de facel
  128. C NGCF = numero global du centre de facel
  129. C NLCF1 = numero local du centre de face
  130. C NGCEG = numero global du centre ELT "gauche"
  131. C NLCEG = numero local du centre ELT "gauche"
  132. C NGCED = numero global du centre ELT "droite"
  133. C NLCED = numero local du centre ELT "droite"
  134. C
  135. NGCF = MELEMF.NUM(1,NLCF)
  136. NGCF1 = MELEFL.NUM(2,NLCF)
  137. IF(NGCF .NE. NGCF1)THEN
  138. MOTERR(1:40)= 'FACEL et FACE = ? '
  139. CALL ERREUR(5)
  140. GOTO 9999
  141. ENDIF
  142. C
  143. NGCEG = MELEFL.NUM(1,NLCF)
  144. NGCED = MELEFL.NUM(3,NLCF)
  145. NLCEG = MLCENT.LECT(NGCEG)
  146. NLCED = MLCENT.LECT(NGCED)
  147. C-------------------------------------------------------
  148. IF(NGCEG .NE. NGCED)THEN
  149. C
  150. C********** Parametres geometriques
  151. C
  152. ICOORX = ((IDIM + 1) * (NGCF - 1))+1
  153. XF = MCOORD.XCOOR(ICOORX)
  154. YF = MCOORD.XCOOR(ICOORX+1)
  155. C
  156. ICOORX = ((IDIM + 1) * (NGCEG - 1))+1
  157. XG = MCOORD.XCOOR(ICOORX)
  158. YG = MCOORD.XCOOR(ICOORX+1)
  159. XFMXG = XF - XG
  160. YFMYG = YF - YG
  161. DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG))
  162. C
  163. ICOORX = ((IDIM + 1) * (NGCED - 1))+1
  164. XD = MCOORD.XCOOR(ICOORX)
  165. YD = MCOORD.XCOOR(ICOORX+1)
  166. XFMXD = XF - XD
  167. YFMYD = YF - YD
  168. DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD))
  169. C
  170. C********** F=G -> DRG = 0 -> ALPHA = 0
  171. C
  172. ALPHA=DRG/(DRG+DRD)
  173. UMALPH= 1.0D0 - ALPHA
  174. C
  175. C********** Les valeurs à l'interface
  176. C
  177. C DRG=0 -> F=G
  178. C
  179. ELSE
  180. C
  181. C********** MURS
  182. C
  183. C Etat a gauche = Etat droite
  184. C
  185. ALPHA=0.0D0
  186. UMALPH=1.0D0
  187. C
  188. C********** Parametres geometriques
  189. C
  190. ICOORX = ((IDIM + 1) * (NGCF - 1))+1
  191. XF = MCOORD.XCOOR(ICOORX)
  192. YF = MCOORD.XCOOR(ICOORX+1)
  193. C
  194. ICOORX = ((IDIM + 1) * (NGCEG - 1))+1
  195. XG = MCOORD.XCOOR(ICOORX)
  196. YG = MCOORD.XCOOR(ICOORX+1)
  197. XFMXG = XF - XG
  198. YFMYG = YF - YG
  199. ENDIF
  200. C
  201. C******* On calcule le sign du pruduit scalare
  202. C (Normales de Castem) * (vecteur "gauche" -> "centre")
  203. C
  204. CNX = MPNORM.VPOCHA(NLCF,1)
  205. CNY = MPNORM.VPOCHA(NLCF,2)
  206. ORIENT = CNX * XFMXG + CNY * YFMYG
  207. ORIENT = SIGN(1.0D0,ORIENT)
  208. IF(ORIENT .NE. 1.0D0)THEN
  209. MOTERR(1:40)=
  210. & 'LAPN , subroutine zlap12.eso. '
  211. WRITE(IOIMP,*) MOTERR(1:40)
  212. CALL ERREUR(5)
  213. GOTO 9999
  214. ENDIF
  215. C
  216. C******* Les flux aux interfaces
  217. C
  218. SURF = MPSURF.VPOCHA(NLCF,1)
  219. DIAM = UMALPH*MPDIAM.VPOCHA(NLCEG,1) +
  220. & ALPHA*MPDIAM.VPOCHA(NLCED,1)
  221. DIAM2=DIAM*DIAM
  222. MU = UMALPH*MPMUC.VPOCHA(NLCEG,1) +
  223. & ALPHA*MPMUC.VPOCHA(NLCED,1)
  224. MUT = UMALPH*MPMTC.VPOCHA(NLCEG,1) +
  225. & ALPHA*MPMTC.VPOCHA(NLCED,1)
  226. c------------------------------------------
  227. DKADX=MPGRKE.VPOCHA(NLCF,1)
  228. DKADY=MPGRKE.VPOCHA(NLCF,2)
  229. DEPDX=MPGRKE.VPOCHA(NLCF,3)
  230. DEPDY=MPGRKE.VPOCHA(NLCF,4)
  231. * Contribution pour l'espèce IESP
  232. MPFLUX.VPOCHA(NLCF,IDIM+2+NESP)=
  233. $ (MU+(MUT/TSIGK))*((DKADX*CNX)+(DKADY*CNY))*SURF*(-1.D0)
  234. * Contribution pour l'énergie totale
  235. MPFLUX.VPOCHA(NLCF,IDIM+3+NESP)=
  236. $ (MU+(MUT/TSIGE))*((DEPDX*CNX)+(DEPDY*CNY))*SURF*(-1.D0)
  237. * Le pas de temps (proportional to \rho) - not implemented
  238. CELL=(4.0D0*(MU+MUT))/DIAM2
  239. IF(CELL .GT. UNSDT)THEN
  240. UNSDT=CELL
  241. ENDIF
  242. C---------------------------------------------------
  243. ENDDO
  244. C
  245. C-----------------------------------------------
  246. DT = 1.0D0 / (UNSDT+XPETIT)
  247. C
  248. SEGDES MELEFL
  249. SEGDES MELEMF
  250. SEGDES MELEMC
  251. SEGDES MPSURF
  252. SEGDES MPNORM
  253. SEGDES MPDIAM
  254. SEGDES MLCENT
  255. C
  256. SEGDES MPMUC
  257. SEGDES MPMTC
  258. SEGDES MPGRKE
  259. SEGDES MPFLUX
  260. C
  261. C
  262. 9999 CONTINUE
  263. RETURN
  264. END
  265.  
  266.  
  267.  
  268.  
  269.  
  270.  
  271.  
  272.  
  273.  
  274.  
  275.  

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