Télécharger la plaquette Cast3M


fr en

Télécharger tlap13.eso

Retour à la liste

Numérotation des lignes : 

tlap13

  1. C TLAP13    SOURCE    CB215821  20/11/25    13:40:59     10792          
  2.       SUBROUTINE TLAP13(NESP,IMUC,IMTC,TSIGK,TSIGE,IGRKEP,
  3.      &     MELEMC,MELEMF,MELEFL,ISURF,INORM,IDIAM,
  4.      &     ICHFLU,DT)
  5. C***********************************************************************
  6. C NOM         : TLAP13
  7. C DESCRIPTION : Calcul des flux diffusifs pour k-\eps equations de
  8. C               turbulence (3D)
  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,ZG,XFMXG,YFMYG,ZFMZG,DRG
  81.      &       ,XD,YD,ZD,XFMXD,YFMYD,ZFMZD,DRD,ALPHA,UMALPH
  82.      &       ,XF,YF,ZF
  83.      &       ,CNX, CNY, CNZ, ORIENT, DIAM, DIAM2, CELL, SURF
  84.       INTEGER NESP
  85.       REAL*8 DKADX,DKADY,DKADZ,DEPDX,DEPDY,DEPDZ
  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.             ZF = MCOORD.XCOOR(ICOORX+2)
  156. C
  157.             ICOORX = ((IDIM + 1) * (NGCEG - 1))+1
  158.             XG = MCOORD.XCOOR(ICOORX)
  159.             YG = MCOORD.XCOOR(ICOORX+1)
  160.             ZG = MCOORD.XCOOR(ICOORX+2)
  161.             XFMXG = XF - XG
  162.             YFMYG = YF - YG
  163.             ZFMZG = ZF - ZG
  164.             DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG)+(ZFMZG*ZFMZG))
  165. C
  166.             ICOORX = ((IDIM + 1) * (NGCED - 1))+1
  167.             XD = MCOORD.XCOOR(ICOORX)
  168.             YD = MCOORD.XCOOR(ICOORX+1)
  169.             ZD = MCOORD.XCOOR(ICOORX+2)
  170.             XFMXD = XF - XD
  171.             YFMYD = YF - YD
  172.             ZFMZD = ZF - ZD
  173.             DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD)+(ZFMZD*ZFMZD))
  174. C
  175. C********** F=G -> DRG = 0 -> ALPHA = 0
  176. C
  177.             ALPHA=DRG/(DRG+DRD)
  178.             UMALPH= 1.0D0 - ALPHA
  179. C
  180. C********** Les valeurs à l'interface
  181. C
  182. C           DRG=0 -> F=G
  183. C
  184.          ELSE
  185. C
  186. C********** MURS
  187. C
  188. C           Etat a gauche = Etat droite
  189. C
  190.             ALPHA=0.0D0
  191.             UMALPH=1.0D0
  192. C
  193. C********** Parametres geometriques
  194. C
  195.             ICOORX = ((IDIM + 1) * (NGCF - 1))+1
  196.             XF = MCOORD.XCOOR(ICOORX)
  197.             YF = MCOORD.XCOOR(ICOORX+1)
  198.             ZF = MCOORD.XCOOR(ICOORX+2)
  199. C
  200.             ICOORX = ((IDIM + 1) * (NGCEG - 1))+1
  201.             XG = MCOORD.XCOOR(ICOORX)
  202.             YG = MCOORD.XCOOR(ICOORX+1)
  203.             ZG = MCOORD.XCOOR(ICOORX+2)
  204.             XFMXG = XF - XG
  205.             YFMYG = YF - YG
  206.             ZFMZG = ZF - ZG
  207.          ENDIF
  208. C
  209. C******* On calcule le sign du pruduit scalare
  210. C        (Normales de Castem) * (vecteur "gauche" -> "centre")
  211. C
  212.          CNX = MPNORM.VPOCHA(NLCF,1)
  213.          CNY = MPNORM.VPOCHA(NLCF,2)
  214.          CNZ = MPNORM.VPOCHA(NLCF,3)
  215.          ORIENT=(CNX * XFMXG)+(CNY * YFMYG)+(CNZ*ZFMZG)
  216.          ORIENT = SIGN(1.0D0,ORIENT)
  217.          IF(ORIENT .NE. 1.0D0)THEN
  218.             MOTERR(1:40)=
  219.      &           'LAPN , subroutine zlap12.eso.           '
  220.             WRITE(IOIMP,*) MOTERR(1:40)
  221.             CALL ERREUR(5)
  222.             GOTO 9999
  223.          ENDIF
  224. C
  225. C******* Les flux aux interfaces
  226. C
  227.          SURF = MPSURF.VPOCHA(NLCF,1)
  228.          DIAM = UMALPH*MPDIAM.VPOCHA(NLCEG,1) +
  229.      &        ALPHA*MPDIAM.VPOCHA(NLCED,1)
  230.          DIAM2=DIAM*DIAM
  231.          MU =  UMALPH*MPMUC.VPOCHA(NLCEG,1) +
  232.      &        ALPHA*MPMUC.VPOCHA(NLCED,1)
  233.          MUT = UMALPH*MPMTC.VPOCHA(NLCEG,1) +
  234.      &        ALPHA*MPMTC.VPOCHA(NLCED,1)
  235. c------------------------------------------
  236.             DKADX=MPGRKE.VPOCHA(NLCF,1)
  237.             DKADY=MPGRKE.VPOCHA(NLCF,2)
  238.             DKADZ=MPGRKE.VPOCHA(NLCF,3)
  239.             DEPDX=MPGRKE.VPOCHA(NLCF,4)
  240.             DEPDY=MPGRKE.VPOCHA(NLCF,5)
  241.             DEPDZ=MPGRKE.VPOCHA(NLCF,6)
  242. * Contribution pour l'espèce IESP
  243.             MPFLUX.VPOCHA(NLCF,IDIM+2+NESP)=(MU+(MUT/TSIGK))*
  244.      $      ((DKADX*CNX)+(DKADY*CNY)+(DKADZ*CNZ))*SURF*(-1.D0)
  245. * Contribution pour l'énergie totale
  246.             MPFLUX.VPOCHA(NLCF,IDIM+3+NESP)=(MU+(MUT/TSIGE))*
  247.      $      ((DEPDX*CNX)+(DEPDY*CNY)+(DEPDZ*CNZ))*SURF*(-1.D0)
  248. * Le pas de temps
  249.             CELL=(4.0D0*(MU+MUT))/DIAM2
  250.             IF(CELL .GT. UNSDT)THEN
  251.                UNSDT=CELL
  252.             ENDIF
  253. C---------------------------------------------------
  254.       ENDDO
  255. C
  256. C
  257.       DT = 1.0D0 / (UNSDT+XPETIT)
  258. C
  259.       SEGDES MELEFL
  260.       SEGDES MELEMF
  261.       SEGDES MELEMC
  262.       SEGDES MPSURF
  263.       SEGDES MPNORM
  264.       SEGDES MPDIAM
  265.       SEGDES MLCENT
  266. C
  267.       SEGDES MPMUC
  268.       SEGDES MPMTC
  269.       SEGDES MPGRKE
  270.       SEGDES MPFLUX
  271. C
  272. C
  273.  9999 CONTINUE
  274.       RETURN
  275.       END
  276.  
  277.  
  278.  
  279.  
  280.  
  281.  
  282.  
  283.  
  284.  
  285.  
  286.  

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