Sources Esope | Cast3M

yla12t
C YLA12T    SOURCE    OF166741  24/12/13    21:17:38     12097          
C YLAP12    SOURCE    GOUNAND   01/09/10    21:16:17     4198
      SUBROUTINE YLA12T(IGRTC,IQIMP,
     &     MELEMC,MELEMF,MELEFL,ISURF,INORM,IDIAM,
     &     ICHFLU,DT)
C
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  YLA12T
C
C DESCRIPTION       :  Subroutine appellée par YLAP11
C
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
C
C AUTEUR            :  A. BECCANTINI, DRN/DMT/SEMT/LTMF
C
C************************************************************************
C
C
C     ENTRÉES:
C     *******
C
C     MU    (REAL*8) : viscosité dynamique (kg/m^3 * m^2/s dans le SI)
C
C     KAPPA (REAL*8) : conductivité thermique (J / (s m K))
C
C     CV    (REAL*8) : chaleur spécifique à volume constant (J / (kg K))
C
C     IROC           : pointeur du CHAMPOINT "CENTRE" densité (kg/m^3)
C
C     IVITC          : pointeur du CHAMPOINT "CENTRE" vitesse
C
C     IGRVC          : pointeur du CHAMPOINT "FACE" gradient de vitesse
C
C     IGRTC          : pointeur du CHAMPOINT "FACE" gradient de
C                      température
C
C     IVIMP          : pointeur de CHAMPOINT vitesse imposé (sur des
C                      points FACE)
C
C     ITAUIM         : pointeur de CHAMPOINT tenseur de contraintes
C                      visqueux imposé (sur des points FACE)
C
C     IQIMP          : pointeur de CHAMPOINT flux de chaleur imposé
C                      (sur des points FACE)
C
C     MELEMC         : pointeur du maillage "CENTRE"
C
C     MELEMF         : pointeur du maillage "FACE"
C
C     MELEFL         : pointeur du maillage "FACEL"
C
C     ISURF          : pointeur du CHAMPOINT "FACE" qui contient les
C                      surfaces de faces
C
C     INORM          : pointeur du CHAMPOINT "FACE" qui contient les
C                      normales aux facesP12
C
C     IDIAM          : pointeur du CHAMPOINT "CENTRE" qui contient le
C                      diamètre des elts
C
C
C     SORTIES
C     *******
C
C     ICHFLU         : pointeur du CHAMPOINT "FACE" qui contient les
C                      flux diffusives aux interface
C
C     DT (REAL*8)    : pas de temps de stabilité donné par le critère
C                      de Fourier
C
C***********************************************************************
C
C**** N.B.: Traitement des conditions aux bords
C
C     'VIMP' : la vitesse imposé n'importe ou!
C
C     'QIMP' : flux de chaleur imposé n'importe ou
C
C     'TAUI' : tenseur de contraintes visqueux imposé n'importe ou
C
C
      IMPLICIT INTEGER(I-N)
 
-INC PPARAM
-INC CCOPTIO
-INC CCREEL
-INC SMCHPOI
-INC SMELEME
-INC SMCOORD
-INC SMLENTI
-INC SMLMOTS
C
      POINTEUR MPROC.MPOVAL, MPVITC.MPOVAL, MPGRVF.MPOVAL,
     &     MPGRTF.MPOVAL,
     &     MPVIMP.MPOVAL, MPTAUI.MPOVAL, MPQIMP.MPOVAL,
     &     MPSURF.MPOVAL, MPNORM.MPOVAL, MPDIAM.MPOVAL,
     &     MPFLUX.MPOVAL
C
      POINTEUR MELEMC.MELEME,MELEMF.MELEME,MELEFL.MELEME
      POINTEUR MLCENT.MLENTI,MLEVIM.MLENTI,MLEQIM.MLENTI,MLETAI.MLENTI
C
      INTEGER IROC,IVITC,IGRVC,IGRTC,IVIMP,ITAUIM,IQIMP
     &       ,ISURF,INORM,IDIAM,ICHFLU
     &       ,NFAC, NLCF, NGCF, NGCF1, NGCEG, NGCED
     &       ,NLCEG,NLCED,NLFVI,NLFTI,NLFQI
     &       ,IGEOM,ICORX1,ICORXG,ICORXD
 
      REAL*8 MU, DT, UNSDT
     &       ,UXG,UYG
     &       ,XG,YG,XFMXG,YFMYG,DRG
     &       ,UXD,UYD
     &       ,XD,YD,XFMXD,YFMYD,DRD,ALPHA,UMALPH
     &       ,UXF,UYF,DUXXF,DUXYF,DUYXF,DUYYF,DTXF,DTYF
     &       ,DSTDU,TAUXX,TAUXY,TAUYY,QX,QY,XF,YF
     &       ,CNX, CNY, ORIENT, RO, DIAM, DIAM2, CELL, SURF
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(IGRTC,MPGRTF,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 MPVITC
C     SEGACT*MOD MPGRVF
C     SEGACT*MOD MPGRTF
C     SEGACT*MOD MPSURF
C     SEGACT*MOD MPNORM
C     SEGACT*MOD MPDIAM
C     SEGACT*MOD MPFLUX
C
      IF(IQIMP .NE. 0)THEN
         CALL LICHT(IQIMP,MPQIMP,TYPE,IGEOM)
C        SEGACT*MOD MPQIMP
         CALL KRIPAD(IGEOM,MLEQIM)
C        SEGACT IGEOM
C        SEGACT MLEQIM
         MELEME = IGEOM
C         SEGDES MELEME
      ENDIF
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
C******* On controlle si sur NGCF on impose de CL
C
C        NLFVI  = numero local du centre de face sul le maillage des
C                 "vitesses" "imposées"
C
C        NLFTI  = numero local du centre de face sul le maillage des
C                 "tau" "imposés"
C
C        NLFQI  = numero local du centre de face sul le maillage des
C                 "q" "imposés"
C
         IF(IQIMP .NE. 0)THEN
            NLFQI = MLEQIM.LECT(NGCF)
         ELSE
            NLFQI = 0
         ENDIF
C
         IF(NGCEG .NE. NGCED)THEN
C
C********** Parametres geometriques
C
 
C CB215821 : Modification pour faire le calcul independamment de la dimension...
           ICORX1 = ((IDIM + 1) * (NGCF  - 1))+1
           ICORXG = ((IDIM + 1) * (NGCEG - 1))+1
           ICORXD = ((IDIM + 1) * (NGCED - 1))+1
           DRG    = 0.
           DRD    = 0.
 
           DO ILADIM=1,IDIM
              XF    = MCOORD.XCOOR(ICORX1+(ILADIM-1))
              XG    = MCOORD.XCOOR(ICORXG+(ILADIM-1))
              XD    = MCOORD.XCOOR(ICORXD+(ILADIM-1))
              XFMXG = XF - XG
              XFMXD = XF - XD
              DRG   = DRG + (XFMXG*XFMXG)
              DRD   = DRD + (XFMXD*XFMXD)
           ENDDO
 
           DRG=SQRT(DRG)
           DRD=SQRT(DRD)
C CB215821 Fin de la modification...
 
C
C********** F=G -> DRG = 0 -> ALPHA = 0
           ALPHA=DRG/(DRG+DRD)
           UMALPH= 1.0D0 - ALPHA
C
C
C********** Les valeurs à l'interface
C
C          DRG=0 -> F=G
C
C
C********** Flux de chaleur
C
c           IF(NLFQI .GT. 0)THEN
c              QX = MPQIMP.VPOCHA(NLFQI,1)
c           ELSE
C************* Gauche
           DTXF = MPGRTF.VPOCHA(NLCF,1)
C
C CORRECTION
           QX =  -1.0D0 * DTXF
C
c           ENDIF
         ELSE
C
C********** MURS
C
C           Etat a gauche = Etat droite
C
 
            ALPHA=0.0D0
            UMALPH=1.0D0
C
C********** Parametres geometriques
C
 
C CB215821 : Tout ce qui est la semble inutile...
C            ICORX1 = ((IDIM + 1) * (NGCF - 1))+1
C            XF = MCOORD.XCOOR(ICORX1)
C            YF = MCOORD.XCOOR(ICORX1+1)
C
C            ICORXG = ((IDIM + 1) * (NGCEG - 1))+1
C            XG = MCOORD.XCOOR(ICORXG)
C            YG = MCOORD.XCOOR(ICORXG+1)
C            XFMXG = XF - XG
C            YFMYG = YF - YG
 
C
C********** Flux de chaleur
C
c            IF(NLFQI .GT. 0)THEN
c               QX = MPQIMP.VPOCHA(NLFQI,1)
c            ELSE
C************* Gauche
               DTXF = MPGRTF.VPOCHA(NLCF,1)
c               DTYF = MPGRTF.VPOCHA(NLCF,2)
C
C CORRECTION
               QX = -1.0D0 * DTXF
C
c            ENDIF
         ENDIF
C
C******* On calcule le sign du pruduit scalare
C        (Normales de Castem) * (vecteur "gauche" -> "centre")
C
c         CNX = MPNORM.VPOCHA(NLCF,1)
c         CNY = MPNORM.VPOCHA(NLCF,2)
c         ORIENT = CNX * XFMXG + CNY * YFMYG
c         ORIENT = SIGN(1.0D0,ORIENT)
c         IF(ORIENT .NE. 1.0D0)THEN
c            MOTERR(1:40)=
c     &           'LAPN , subroutine ylap12.eso.           '
c            WRITE(IOIMP,*) MOTERR(1:40)
c            CALL ERREUR(5)
c            GOTO 9999
c         ENDIF
C
C******* Le flux aux interfaces
C
         SURF = MPSURF.VPOCHA(NLCF,1)
         MPFLUX.VPOCHA(NLCF,1) =
     &        (  QX ) * SURF
C
C******  Le pas de temps
C
c         RO=UMALPH +
c     &        ALPHA
c         DIAM = UMALPH*MPDIAM.VPOCHA(NLCEG,1) +
c     &        ALPHA*MPDIAM.VPOCHA(NLCED,1)
c         DIAM2=DIAM*DIAM
c         CELL  = 4.0D0*MU / (DIAM2*RO)
c         CELL  = MAX(CELL, (4.0D0/(DIAM2*RO)))
C
c         IF(CELL .GT. UNSDT)THEN
c            UNSDT  =  CELL
c         ENDIF
C
      ENDDO
C
C
c      DT = 1.0D0 / (UNSDT + XPETIT)
      DT = 1.0D0
C
C      SEGDES MELEFL
C      SEGDES MELEMF
C      SEGDES MELEMC
C      SEGDES MPSURF
C      SEGDES MPNORM
C      SEGDES MPDIAM
      SEGSUP MLCENT
C
C      SEGDES MPGRTF
C      SEGDES MPFLUX
C
C      IF(IQIMP .NE. 0)THEN
C         SEGDES MPQIMP
C         SEGDES MLEQIM
C      ENDIF
C
 9999 CONTINUE
      END