Sources Esope | Cast3M

tlap12
C TLAP12    SOURCE    CB215821  20/11/25    13:40:58     10792          
      SUBROUTINE TLAP12(NESP,IMUC,IMTC,TSIGK,TSIGE,IGRKEP,
     &     MELEMC,MELEMF,MELEFL,ISURF,INORM,IDIAM,
     &     ICHFLU,DT)
C***********************************************************************
C NOM         : TLAP12
C DESCRIPTION : Calcul des flux diffusifs pour k-\eps equations de
C               turbulence (2D)
C
C
C LANGAGE     : ESOPE
C AUTEUR      : S. Kudriakov (CEA/DEN/DM2S/SFME/LTMF)
C               mél : skudriakov@cea.fr
C***********************************************************************
C APPELES (UTIL)   : KRIPAD : MELEME -> (num. globale->locale)
C                    LICHT  : Lecture des pointeurs (maillages, valeurs)
C                             d'un objet de type MCHPOI.
C                    ERREUR : Gestion des erreurs par GIBI.
C APPELE PAR       : ZLAP11 : Chapeau de l'opérateur Gibiane 'LAPN'
C                             option 'VF'.
C***********************************************************************
C ENTREES            : NESP (type ENTIER) : number of species explicitly
C                                        treated in the NS equations
C                      IMUC (type MCHPOI) : laminar viscosity
C                      IMTC (type MCHPOI) : turbulent viscosity
C                      TSIGK (type REEL) : turbulent constant \sigma_k
C                      TSIGE (type REEL) : turbulent constant
C                                                     \sigma_{\eps}
C                      IGRKEP (type MCHPOI) : gradients of k and \epsilon
C
C                      MELEMC (type MELEME) : maillage des centres des
C                                             éléments.
C                      MELEMF (type MELEME) : maillage des faces des
C                                             éléments.
C                      MELEFL (type MELEME) : connectivités face-(centre
C                                             gauche, centre droit).
C                      ISURF (type MCHPOI) : surface des faces.
C                      INORM (type MCHPOI) : normale aux faces.
C                      IDIAM (type MCHPOI) : diamètre des éléments.
C ENTREES/SORTIES    : -
C SORTIES            : ICHFLU (type MCHPOI) : flux diffusif aux
C                                             interfaces.
C                      DT     (type REAL*8) : pas de temps de stabilité
C                                             (Fourier)
C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
C***********************************************************************
C VERSION    : v1, 3/12/2003, version initiale
C HISTORIQUE :
C***********************************************************************
C Prière de PRENDRE LE TEMPS de compléter les commentaires
C en cas de modification de ce sous-programme afin de faciliter
C la maintenance !
C***********************************************************************
      IMPLICIT INTEGER(I-N)
 
-INC PPARAM
-INC CCOPTIO
-INC CCREEL
-INC SMCHPOI
-INC SMELEME
-INC SMCOORD
-INC SMLENTI
-INC SMLMOTS
C
      POINTEUR IMUC.MCHPOI,IMTC.MCHPOI,IGRKEP.MCHPOI
      POINTEUR MPMUC.MPOVAL,MPMTC.MPOVAL,MPGRKE.MPOVAL,
     &     MPSURF.MPOVAL, MPNORM.MPOVAL, MPDIAM.MPOVAL,
     &     MPFLUX.MPOVAL
C
      POINTEUR MELEMC.MELEME,MELEMF.MELEME,MELEFL.MELEME
      POINTEUR MLCENT.MLENTI,MLETIM.MLENTI,MLERIM.MLENTI
C
      INTEGER ITIMP
     &       ,ISURF,INORM,IDIAM,ICHFLU
     &       ,NFAC, NLCF, NGCF, NGCF1, NGCEG, NGCED
     &       ,NLCEG,NLCED,NLFTI,NLFRI
     &       , ICOORX, IGEOM
 
      REAL*8 DT, UNSDT
     &       ,XG,YG,XFMXG,YFMYG,DRG
     &       ,XD,YD,XFMXD,YFMYD,DRD,ALPHA,UMALPH
     &       ,XF,YF
     &       ,CNX, CNY, ORIENT, DIAM, DIAM2, CELL, SURF
      INTEGER NESP
      REAL*8 DKADX,DKADY,DEPDX,DEPDY
      REAL*8 MU,MUT,TSIGK,TSIGE
C
      CHARACTER*8 TYPE
C
C**** Initialisation de 1/DT
C
      UNSDT = 0.0D0
C
C**** KRIPAD pour la correspondance global/local de centre
C
      CALL KRIPAD(MELEMC,MLCENT)
C
C     EN KRIPAD
C     SEGACT MELEMC
C     SEGACT MLCENT
C
      CALL LICHT(IMUC,MPMUC,TYPE,IGEOM)
      CALL LICHT(IMTC,MPMTC,TYPE,IGEOM)
      CALL LICHT(IGRKEP,MPGRKE,TYPE,IGEOM)
      CALL LICHT(ISURF,MPSURF,TYPE,IGEOM)
      CALL LICHT(INORM,MPNORM,TYPE,IGEOM)
      CALL LICHT(IDIAM,MPDIAM,TYPE,IGEOM)
      CALL LICHT(ICHFLU,MPFLUX,TYPE,IGEOM)
C
C     EN LICHT
C     SEGACT*MOD MPROC
C     SEGACT*MOD MPTEMC
C     SEGACT*MOD MPSURF
C     SEGACT*MOD MPNORM
C     SEGACT*MOD MPDIAM
C     SEGACT*MOD MPFLUX
C
C---------------------------------------------------------
      SEGACT MELEFL
      SEGACT MELEMF
      NFAC = MELEMF.NUM(/2)
C
C**** Boucle sur les faces
C
      DO NLCF = 1, NFAC, 1
C
C******* NLCF  = numero local du centre de facel
C        NGCF  = numero global du centre de facel
C        NLCF1 = numero local du centre de face
C        NGCEG = numero global du centre ELT "gauche"
C        NLCEG = numero local du centre ELT "gauche"
C        NGCED = numero global du centre ELT "droite"
C        NLCED = numero local du centre ELT "droite"
C
         NGCF = MELEMF.NUM(1,NLCF)
         NGCF1 = MELEFL.NUM(2,NLCF)
         IF(NGCF .NE. NGCF1)THEN
            MOTERR(1:40)= 'FACEL et FACE = ?                       '
            CALL ERREUR(5)
            GOTO 9999
         ENDIF
C
         NGCEG = MELEFL.NUM(1,NLCF)
         NGCED = MELEFL.NUM(3,NLCF)
         NLCEG = MLCENT.LECT(NGCEG)
         NLCED = MLCENT.LECT(NGCED)
C-------------------------------------------------------
         IF(NGCEG .NE. NGCED)THEN
C
C********** Parametres geometriques
C
            ICOORX = ((IDIM + 1) * (NGCF - 1))+1
            XF = MCOORD.XCOOR(ICOORX)
            YF = MCOORD.XCOOR(ICOORX+1)
C
            ICOORX = ((IDIM + 1) * (NGCEG - 1))+1
            XG = MCOORD.XCOOR(ICOORX)
            YG = MCOORD.XCOOR(ICOORX+1)
            XFMXG = XF - XG
            YFMYG = YF - YG
            DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG))
C
            ICOORX = ((IDIM + 1) * (NGCED - 1))+1
            XD = MCOORD.XCOOR(ICOORX)
            YD = MCOORD.XCOOR(ICOORX+1)
            XFMXD = XF - XD
            YFMYD = YF - YD
            DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD))
C
C********** F=G -> DRG = 0 -> ALPHA = 0
C
            ALPHA=DRG/(DRG+DRD)
            UMALPH= 1.0D0 - ALPHA
C
C********** Les valeurs à l'interface
C
C           DRG=0 -> F=G
C
         ELSE
C
C********** MURS
C
C           Etat a gauche = Etat droite
C
            ALPHA=0.0D0
            UMALPH=1.0D0
C
C********** Parametres geometriques
C
            ICOORX = ((IDIM + 1) * (NGCF - 1))+1
            XF = MCOORD.XCOOR(ICOORX)
            YF = MCOORD.XCOOR(ICOORX+1)
C
            ICOORX = ((IDIM + 1) * (NGCEG - 1))+1
            XG = MCOORD.XCOOR(ICOORX)
            YG = MCOORD.XCOOR(ICOORX+1)
            XFMXG = XF - XG
            YFMYG = YF - YG
         ENDIF
C
C******* On calcule le sign du pruduit scalare
C        (Normales de Castem) * (vecteur "gauche" -> "centre")
C
         CNX = MPNORM.VPOCHA(NLCF,1)
         CNY = MPNORM.VPOCHA(NLCF,2)
         ORIENT = CNX * XFMXG + CNY * YFMYG
         ORIENT = SIGN(1.0D0,ORIENT)
         IF(ORIENT .NE. 1.0D0)THEN
            MOTERR(1:40)=
     &           'LAPN , subroutine zlap12.eso.           '
            WRITE(IOIMP,*) MOTERR(1:40)
            CALL ERREUR(5)
            GOTO 9999
         ENDIF
C
C******* Les flux aux interfaces
C
         SURF = MPSURF.VPOCHA(NLCF,1)
         DIAM = UMALPH*MPDIAM.VPOCHA(NLCEG,1) +
     &        ALPHA*MPDIAM.VPOCHA(NLCED,1)
         DIAM2=DIAM*DIAM
         MU =  UMALPH*MPMUC.VPOCHA(NLCEG,1) +
     &        ALPHA*MPMUC.VPOCHA(NLCED,1)
         MUT = UMALPH*MPMTC.VPOCHA(NLCEG,1) +
     &        ALPHA*MPMTC.VPOCHA(NLCED,1)
c------------------------------------------
            DKADX=MPGRKE.VPOCHA(NLCF,1)
            DKADY=MPGRKE.VPOCHA(NLCF,2)
            DEPDX=MPGRKE.VPOCHA(NLCF,3)
            DEPDY=MPGRKE.VPOCHA(NLCF,4)
* Contribution pour l'espèce IESP
            MPFLUX.VPOCHA(NLCF,IDIM+2+NESP)=
     $      (MU+(MUT/TSIGK))*((DKADX*CNX)+(DKADY*CNY))*SURF*(-1.D0)
* Contribution pour l'énergie totale
            MPFLUX.VPOCHA(NLCF,IDIM+3+NESP)=
     $      (MU+(MUT/TSIGE))*((DEPDX*CNX)+(DEPDY*CNY))*SURF*(-1.D0)
* Le pas de temps (proportional to \rho) - not implemented
            CELL=(4.0D0*(MU+MUT))/DIAM2
            IF(CELL .GT. UNSDT)THEN
               UNSDT=CELL
            ENDIF
C---------------------------------------------------
      ENDDO
C
C-----------------------------------------------
      DT = 1.0D0 / (UNSDT+XPETIT)
C
      SEGDES MELEFL
      SEGDES MELEMF
      SEGDES MELEMC
      SEGDES MPSURF
      SEGDES MPNORM
      SEGDES MPDIAM
      SEGDES MLCENT
C
      SEGDES MPMUC
      SEGDES MPMTC
      SEGDES MPGRKE
      SEGDES MPFLUX
C
C
 9999 CONTINUE
      RETURN
      END