Télécharger ylap12.eso

Retour à la liste

ylap12
C YLAP12    SOURCE    OF166741  24/12/13    21:17:40     12097          
      SUBROUTINE YLAP12(MU,KAPPA,CV,IROC,IVITC,IGRVC,
     &     IGRTC,IVIMP,ITAUIM,IQIMP,
     &     MELEMC,MELEMF,MELEFL,ISURF,INORM,IDIAM,
     &     ICHFLU,DT)
C
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  YLAP12
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 faces
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
     &       , ICOORX, IGEOM
 
      REAL*8 MU,KAPPA,CV,LAMBRO,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
      LAMBRO = KAPPA / CV
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(IROC,MPROC,TYPE,IGEOM)
      CALL LICHT(IVITC,MPVITC,TYPE,IGEOM)
      CALL LICHT(IGRVC,MPGRVF,TYPE,IGEOM)
      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(IVIMP .NE. 0)THEN
         CALL LICHT(IVIMP,MPVIMP,TYPE,IGEOM)
C        SEGACT*MOD MPVIMP
         CALL KRIPAD(IGEOM,MLEVIM)
C        SEGACT IGEOM
C        SEGACT MLEVIM
         MELEME = IGEOM
         SEGDES MELEME
      ENDIF
      IF(ITAUIM .NE. 0)THEN
         CALL LICHT(ITAUIM,MPTAUI,TYPE,IGEOM)
C        SEGACT*MOD MPTAUI
         CALL KRIPAD(IGEOM,MLETAI)
C        SEGACT IGEOM
C        SEGACT MLETAI
         MELEME = IGEOM
         SEGDES MELEME
      ENDIF
      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
         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(IVIMP .NE. 0)THEN
            NLFVI = MLEVIM.LECT(NGCF)
         ELSE
            NLFVI = 0
         ENDIF
C
         IF(ITAUIM .NE. 0)THEN
            NLFTI = MLETAI.LECT(NGCF)
         ELSE
            NLFTI = 0
         ENDIF
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
            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
C
C********** Tenseur de contraintes visqueux
C
            DUXXF = MPGRVF.VPOCHA(NLCF,1)
            DUXYF = MPGRVF.VPOCHA(NLCF,2)
            DUYXF = MPGRVF.VPOCHA(NLCF,3)
            DUYYF = MPGRVF.VPOCHA(NLCF,4)
C
            IF (NLFTI .GT. 0) THEN
               TAUXX = MPTAUI.VPOCHA(NLFTI,1)
               TAUYY = MPTAUI.VPOCHA(NLFTI,2)
               TAUXY = MPTAUI.VPOCHA(NLFTI,3)
            ELSE
               DSTDU = 2.0D0/3.0D0 * (DUXXF + DUYYF)
               TAUXX = MU * (2.0D0 * DUXXF - DSTDU)
               TAUXY = MU * (DUXYF + DUYXF)
               TAUYY = MU * (2.0D0 * DUYYF - DSTDU)
            ENDIF
C
C********** Vitesse
C
            IF( NLFVI .GT. 0) THEN
               UXF = MPVIMP.VPOCHA(NLFVI,1)
               UYF = MPVIMP.VPOCHA(NLFVI,2)
            ELSE
C************* Gauche
               UXG = MPVITC.VPOCHA(NLCEG,1)
               UYG = MPVITC.VPOCHA(NLCEG,2)
C************* Droite
               UXD = MPVITC.VPOCHA(NLCED,1)
               UYD = MPVITC.VPOCHA(NLCED,2)
C************* Face
               UXF = UMALPH * UXG + ALPHA * UXD
               UYF = UMALPH * UYG + ALPHA * UYD
C************* Correction de la vitesse lineaire exacte
               UXF = UXF +
     &              (DUXXF * ((XFMXG * UMALPH)+ (XFMXD * ALPHA))) +
     &              (DUXYF * ((YFMYG * UMALPH)+ (YFMYD * ALPHA)))
               UYF = UYF +
     &              (DUYXF * ((XFMXG * UMALPH)+ (XFMXD * ALPHA))) +
     &              (DUYYF * ((YFMYG * UMALPH)+ (YFMYD * ALPHA)))
            ENDIF
C
C********** Flux de chaleur
C
            IF(NLFQI .GT. 0)THEN
               QX = MPQIMP.VPOCHA(NLFQI,1)
               QY = MPQIMP.VPOCHA(NLFQI,2)
            ELSE
C************* Gauche
               DTXF = MPGRTF.VPOCHA(NLCF,1)
               DTYF = MPGRTF.VPOCHA(NLCF,2)
C
               QX = -1.0D0 * KAPPA * DTXF
               QY = -1.0D0 * KAPPA * DTYF
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
            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
C
C********** Tenseur de contraintes visqueux
C
            DUXXF = MPGRVF.VPOCHA(NLCF,1)
            DUXYF = MPGRVF.VPOCHA(NLCF,2)
            DUYXF = MPGRVF.VPOCHA(NLCF,3)
            DUYYF = MPGRVF.VPOCHA(NLCF,4)
C
            IF (NLFTI .GT. 0) THEN
               TAUXX = MPTAUI.VPOCHA(NLFTI,1)
               TAUYY = MPTAUI.VPOCHA(NLFTI,2)
               TAUXY = MPTAUI.VPOCHA(NLFTI,3)
            ELSE
               DSTDU = 2.0D0/3.0D0 * (DUXXF + DUYYF)
               TAUXX = MU * (2.0D0 * DUXXF - DSTDU)
               TAUXY = MU * (DUXYF + DUYXF)
               TAUYY = MU * (2.0D0 * DUYYF - DSTDU)
            ENDIF
C
C********** Vitesse
C
            IF( NLFVI .GT. 0) THEN
               UXF = MPVIMP.VPOCHA(NLFVI,1)
               UYF = MPVIMP.VPOCHA(NLFVI,2)
            ELSE
               UXF = MPVITC.VPOCHA(NLCEG,1)
               UYF = MPVITC.VPOCHA(NLCEG,2)
C************* Correction de la vitesse lineaire exacte
               UXF = UXF +
     &              (DUXXF * XFMXG ) +
     &              (DUXYF * YFMYG )
               UYF = UYF +
     &              (DUYXF * XFMXG ) +
     &              (DUYYF * YFMYG )
            ENDIF
C
C********** Flux de chaleur
C
            IF(NLFQI .GT. 0)THEN
               QX = MPQIMP.VPOCHA(NLFQI,1)
               QY = MPQIMP.VPOCHA(NLFQI,2)
            ELSE
C************* Gauche
               DTXF = MPGRTF.VPOCHA(NLCF,1)
               DTYF = MPGRTF.VPOCHA(NLCF,2)
C
               QX = -1.0D0 * KAPPA * DTXF
               QY = -1.0D0 * KAPPA * DTYF
C
            ENDIF
         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 ylap12.eso.           '
            WRITE(IOIMP,*) MOTERR(1:40)
            CALL ERREUR(5)
            GOTO 9999
         ENDIF
C
C******* Le flux aux interfaces
C
         SURF = MPSURF.VPOCHA(NLCF,1)
         MPFLUX.VPOCHA(NLCF,1) =  ((TAUXX  * CNX) + (TAUXY  * CNY))
     &        * SURF * (-1.0D0)
         MPFLUX.VPOCHA(NLCF,2) =  ((TAUXY  * CNX) + (TAUYY  * CNY))
     &        * SURF * (-1.0D0)
         MPFLUX.VPOCHA(NLCF,3) = (
     &        ((TAUXX * UXF  + TAUXY * UYF - QX) * CNX) +
     &        ((TAUXY * UXF  + TAUYY * UYF - QY) * CNY))
     &        * SURF * (-1.0D0)
C
C******  Le pas de temps
C
         RO=UMALPH*MPROC.VPOCHA(NLCEG,1) +
     &        ALPHA*MPROC.VPOCHA(NLCED,1)
         DIAM = UMALPH*MPDIAM.VPOCHA(NLCEG,1) +
     &        ALPHA*MPDIAM.VPOCHA(NLCED,1)
         DIAM2=DIAM*DIAM
         CELL  = 4.0D0*MU / (DIAM2*RO)
         CELL  = MAX(CELL, (4.0D0*LAMBRO/(DIAM2*RO)))
C
         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
      SEGSUP MLCENT
C
      SEGDES MPROC
      SEGDES MPVITC
      SEGDES MPGRVF
      SEGDES MPGRTF
      SEGDES MPFLUX
C
      IF(IVIMP .NE. 0) THEN
         SEGDES MPVIMP
         SEGSUP MLEVIM
      ENDIF
      IF(ITAUIM .NE. 0)THEN
         SEGDES MPTAUI
         SEGSUP MLETAI
      ENDIF
      IF(IQIMP .NE. 0)THEN
         SEGDES MPQIMP
         SEGDES MLEQIM
      ENDIF
C
 9999 CONTINUE
      RETURN
      END
 
 
 
 
 
 
 
 
 
 

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