Sources Esope | Cast3M

fcolcp
C FCOLCP    SOURCE    CHAT      05/01/12    23:57:19     5004
      SUBROUTINE FCOLCP(NESP,
     &                  GAMG,ROG,PG,UNG,UTG,
     &                  GAMD,ROD,PD,UND,UTD,
     &                  YG,YD,FLU,
     &                  CELLT,
     &                  LOGNC,LOGAN,MESERR)
C
C*** N.B. Cela est une modification non optimisée de la subroutine FCOLTE:
C
C    OU
C
C    NSCA = 0
C    NORDP1 = 1
C    XST = 0.0D0
C    RTOTG = RTOTD = 1.0D0
C    ACVTOG(1) = RTOTG / (GAMG - 1.0D0)
C    ACVTOD(1) = RTOTG / (GAMD - 1.0D0)
C    TG = PG / ROG
C    TD = PD / ROD
C    ETHERG =  ACVTOG(1) * RTOTG
C    ETHERD =  ACVTOD(1) * RTOTD
C
C
C      SUBROUTINE FCOLTE(NESP,NSCA,NORDP1,
C     &                  XST,
C     &                  GAMG,RTOTG,ACVTOG,ROG,PG,TG,UNG,UTG,ETHERG,
C     &                  GAMD,RTOTD,ACVTOD,ROD,PD,TD,UND,UTD,ETHERD,
C     &                  YG,YD,SCAG,SCAD,
C     &                  FLU,CELLT,
C     &                  LOGNC,MESERR,LOGAN)
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  FCOLTE
C
C DESCRIPTION       :  Formulation Volumes Finis  pour les Equations
C                      d'Euler Multi-Especes relatives à un melange
C                      de gaz parfaits.
C
C                      Calcul du flux aux interfaces avec la methode
C                      de Colella-Glaz, avec les chocs qui remplacent
C                      les ondes de rarefaction + entropy fix.
C
C                      (voir:
C                       1)  BECCANTINI, Thèse de doctorat)
C
C LANGUAGE          :  FORTRAN 77
C
C AUTEUR            :  A. BECCANTINI DRN/DMT/SEMT/LTMF
C
C************************************************************************
C
C APPELES           :  RACCOL
C
C************************************************************************
C
C**** Entrées:
C
C     NESP            =  nombre d'especes considérées dans les Equations
C                        d'Euler
C
C     NSCA            =  nombre de scalaires passifs a transporter
C
C     NORDP1          =  (ordre des polynoms cv(T)) + 1
C
C     XST             =  droite characteristique ou je cherche la solution
C
C     GAMG, GAMD      =  les "gamma" du gaz (gauche et droite)
C
C     RTOTG, RTOTD    =  les constantes des gaz
C
C     ACVTOG, ACVTOD  =  les sommes des coefficients des cv,i
C                        ACVTOT(j) = \sum_{i=1,nesp} Y_i*ACV_{i,j}
C                        j = 1 , NORDP1
C
C     ROG, ROD        =  les densités
C
C     PG, PD          =  les pressions
C
C     TG, TD          =  les temperatures
C
C     UNG, UND        =  vitesses normales
C
C     UTG, UTD        =  vitesses tangentielles
C
C     ETHERG, ETHERD  =  les energies
C
C     YG, YD          =  tables des fractiones massiques
C
C     SCAG, SCAD      =  tables des scalaires passifs
C
C**** Sorties:
C
C     FLU             =  table du flux a l'interface dans le repaire
C                        (n,t), i.e.
C                        (rho*un, rho*un*un + p, rho*un*ut, rho*un*ht,
C                         rho*un*y1, ...)
C
C     CELLT           =  condition de stabilité, i.e.
C
C                        dT/diamax &lt; cellt
C
C     LOGNC           =  si TRUE, la methode de newton pour le calcul
C                        des etats intermediares ne converges pas
C
C     MESERR          =  message d'erreuer
C
C     LOGAN           =  si TRUE, une anomalie a été detectée
C
C************************************************************************
C
C HISTORIQUE (Anomalies et modifications éventuelles)
C
C HISTORIQUE :  Créé le 08.02.00
C
C               21.02.00 transport de scalaires passifs
C
C************************************************************************
C
C N.B.: Toutes les variables sont DECLAREES
C
C
      IMPLICIT INTEGER(I-N)
      INTEGER NESP, NORDP1, I1, NSCA
      REAL*8  GAMG,RTOTG,ACVTOG(1),ROG,PG,TG,UNG,UTG,ETHERG
     &       ,GAMD,RTOTD,ACVTOD(1),ROD,PD,TD,UND,UTD,ETHERD
     &       ,YG(*),YD(*),SCAG(1),SCAD(1)
     &       ,FLU(*),CELLT
     &       ,PSROG,PSROD,ASONG,ASOND,PVAL,UVAL
     &       ,HTHEG,HTHED
     &       ,ROGS,PGS,TGS,ETHEGS,UGS,GAMGS
     &       ,RODS,PDS,TDS,ETHEDS,UDS,GAMDS
     &       ,XST,UNGXST,UNDXST,UGSXST,UDSXST,DGSXST,DDSXST,DGS,DDS
     &       ,ROU,ECFLUX,XLIMG,XLIMD,ULIMG,ULIMD
     &       ,UNGMCG,UNDPCD,PFLUX,ROFLUX,UNFLUX,UTFLUX,HTFLUX
     &       ,FC1,FC2,FPOIS1,FPOIS2,FPOIS3,FPOIS4,SIG12,SIG45,SIG3
     &       ,FPOISG,FPOISD,EPSU
     &       ,UGMCGS,UDPCDS,VITG,VITGS,VITD,VITDS,SIGCHO,SIG123,SIG456
     &       ,SIG56
     &       ,FUNT, TFLUX, RTOTFL, ETHFLU
      CHARACTER*(40) MESERR
      LOGICAL LOGAN, LOGNC
C
C**** Definition pour le passage "thermally-perfect"
C                               -> "calorically perfect"
C
      NSCA = 0
      NORDP1 = 1
      XST = 0.0D0
      RTOTG = 1.0D0
      RTOTD = 1.0D0
      ACVTOG(1) = RTOTG / (GAMG - 1.0D0)
      ACVTOD(1) = RTOTG / (GAMD - 1.0D0)
      TG = PG / (ROG * RTOTG)
      TD = PD / (ROD * RTOTD)
      ETHERG = ACVTOG(1) * TG
      ETHERD = ACVTOD(1) * TD
      SCAG(1) = 0.0D0
      SCAD(1) = 0.0D0
C
C**** YG, YD, FLU
C
C     Dans le cas Euler monoespece, on doit
C     avoir  :
C     YG(1) = YD(1) = 0.0D0
C
      PSROG = PG / ROG
      PSROD = PD / ROD
      ASONG = SQRT(GAMG * PSROG)
      ASOND = SQRT(GAMD * PSROD)
C
C**** Valeurs de référence
C
      PVAL = MAX(PG,PD)
      UVAL = MAX(ASONG,ASOND)
      EPSU = UVAL*1.0D-6
C
C**** États intermédiaires
C
      HTHEG = ETHERG + PSROG
      HTHED = ETHERD + PSROD
      XLIMG = PSROG / (2.0D0 * HTHEG - PSROG)
      XLIMD = PSROD / (2.0D0 * HTHED - PSROD)
      ULIMG = UNG + SQRT( 2.0D0 * HTHEG *
     &       (1.0D0 - XLIMG) / (1.0D0 + XLIMG))
      ULIMD = UND - SQRT( 2.0D0 * HTHED *
     &       (1.0D0 - XLIMD) / (1.0D0 + XLIMD))
      UNGMCG = UNG - ASONG
      UNDPCD = UND + ASOND
      IF(ULIMD .GE. ULIMG)THEN
C
C******* Vide
C
C        N.B.: la surface de contact n'existe pas: on mit le vide
C
         IF((ULIMG .LE. 0.0D0) .AND. (ULIMD .GE. 0.0D0))THEN
            FLU(1) = 0.0D0
            FLU(2) = 0.0D0
            FLU(3) = 0.0D0
            FLU(4) = 0.0D0
            DO I1 = 1, NESP, 1
               FLU(I1+4) = 0.0D0
            ENDDO
            DO I1 = 1, NSCA, 1
               FLU(I1+4+NESP) = 0.0D0
            ENDDO
C
         ELSE
            ROGS = XLIMG * ROG
            RODS = XLIMD * ROD
C
C********** Les functions pois
C
            FC1 = (ULIMG - XST)/(ULIMG - UNGMCG)
            SIG12 = 0.5D0 * (SIGN(1.0D0,FC1)+1.0D0)
            FC2 = 0.5D0*(FC1+1.0D0-ABS(FC1-1.0D0))
            FPOIS1 = SIG12*FC2
C
            FC1 = (XST - ULIMD)/(UNDPCD - ULIMD)
            SIG45 = 0.5D0 * (SIGN(1.0D0,FC1)+1.0D0)
            FC2 = 0.5D0*(FC1+1.0D0-ABS(FC1-1.0D0))
            FPOIS4 = SIG45*FC2
C
            SIG3 = 1.0D0 - (SIG12+SIG45)
            FC1 = (ULIMD - XST)/(ULIMD - ULIMG + EPSU)
            FC1 = 0.5D0*(FC1+1.0D0-ABS(FC1-1.0D0))
            FC2 = 1.0D0 - FC1
            FPOIS2 = SIG12 * (1.0D0 - FPOIS1) + SIG3 * FC1
            FPOIS3 = SIG45 * (1.0D0 - FPOIS4) + SIG3 * FC2
C
C********** La solution pour le calcul du flux
C
            ROFLUX = ROG * FPOIS1 + ROGS * FPOIS2
     &           + RODS * FPOIS3 + ROD * FPOIS4
            TFLUX =  TG * FPOIS1 + TD * FPOIS4
            UNFLUX = UNG * FPOIS1 + ULIMG * FPOIS2
     &           + ULIMD * FPOIS3 + UND * FPOIS4
C
            FPOISG = FPOIS1 + FPOIS2
            FPOISD = FPOIS3 + FPOIS4
            RTOTFL = RTOTG * FPOISG + RTOTD * FPOISD
            PFLUX  = ROFLUX * RTOTFL * TFLUX
            UTFLUX =  UTG * FPOISG + UTD * FPOISD
            ECFLUX = 0.5D0 * (UNFLUX * UNFLUX + UTFLUX * UTFLUX)
C
            FUNT = 1.0D0
            ETHFLU = 0.0D0
            DO I1 = 1, NORDP1, 1
               FUNT = FUNT * TFLUX
               ETHFLU = ETHFLU + (FPOISG * ACVTOG(I1) * FUNT / I1) +
     &              (FPOISD * ACVTOD(I1) * FUNT / I1)
            ENDDO
            HTFLUX = ETHFLU + PFLUX / ROFLUX
C
            ROU = ROFLUX * UNFLUX
            FLU(1) = ROU
            FLU(2) = ROU * UNFLUX + PFLUX
            FLU(3) = ROU * UTFLUX
            FLU(4) = ROU * (HTFLUX + ECFLUX)
            DO I1 = 1, NESP, 1
               FLU(I1+4) = ROU *
     &              (FPOISG * YG(I1) + FPOISD * YD(I1))
            ENDDO
            DO I1 = 1, NSCA, 1
               FLU(I1+4+NESP) = ROU *
     &              (FPOISG * SCAG(I1) + FPOISD * SCAD(I1))
            ENDDO
         ENDIF
         CELLT = 1.0D0/MAX(ABS(UNGMCG),ABS(UNDPCD))
C
C
C******* Test EF
C
C         write(4,'(5D14.6)') XST, FPOIS1, FPOIS2, FPOIS3, FPOIS4
C
      ELSE
         UNGXST=UNG-XST
         UNDXST=UND-XST
         CALL RACCOL(PVAL,UVAL,NORDP1,
     &        RTOTG,ACVTOG,RTOTD,ACVTOD,
     &        ROG,TG,ETHERG,UNGXST,PSROG,GAMG,ASONG,
     &        ROD,TD,ETHERD,UNDXST,PSROD,GAMD,ASOND,
     &        ROGS,PGS,TGS,ETHEGS,UGSXST,DGSXST,GAMGS,
     &        RODS,PDS,TDS,ETHEDS,UDSXST,DDSXST,GAMDS,
     &        LOGNC)
C
         IF(LOGNC)THEN
            MESERR = 'Subroutine RACCOL.ESO, Newton-Raphson.'
C
C******* Message de warning; mais on s'arrete pas
C
         ELSEIF((DGSXST .GT. UGSXST) .OR. (UDSXST .GT. DDSXST))THEN
            LOGAN = .TRUE.
            MESERR = 'Subroutine RACCOL.ESO'
            GOTO 9999
         ENDIF
         UGS=UGSXST+XST
         UDS=UDSXST+XST
         DGS=DGSXST+XST
         DDS=DDSXST+XST
         UGMCGS=UGS-SQRT(GAMGS*PGS/ROGS)
         UDPCDS=UDS+SQRT(GAMDS*PDS/RODS)
 
C
C******* Tests
C
C         write(*,*) ROGS,PGS,TGS,ETHEGS,UGS,GAMGS
C         write(*,*) RODS,PDS,TDS,ETHEDS,UDS,GAMDS
C
C
C******* Les functions pois
C
C        Les vitesses characteristiques des ondes
C
         VITG  = 0.5D0*(DGS+UNGMCG-ABS(DGS-UNGMCG))
         VITGS = 0.5D0*(DGS+UGMCGS+ABS(DGS-UGMCGS))
         VITD  = 0.5D0*(DDS+UNDPCD+ABS(DDS-UNDPCD))
         VITDS = 0.5D0*(DDS+UDPCDS-ABS(DDS-UDPCDS))
C
C******* ZONES 1,2,3
C
         FC1 = (VITGS-XST)/(ABS(VITGS - VITG)+EPSU)
         SIG12 = 0.5D0*(SIGN(1.0D0,FC1)+1.0D0)
         FC2 = 0.5D0*(FC1 + 1.0D0 - ABS(FC1 - 1.0D0))
         SIGCHO = 0.5D0 *(SIGN(1.0D0,(VITG -VITGS+EPSU))+1.0D0)
         SIG123 = 0.5D0 *(SIGN(1.0D0,((0.5D0*(UGS+UDS))-XST))+1.0D0)
C
         FPOIS2 = (1.0D0 - SIGCHO)*(1.0D0 - FC2*SIG12)*SIG123 +
     &        SIGCHO * (SIG123 - SIG12)
         FPOIS1 = (1.0D0 -FPOIS2)*SIG123
C
C******** ZONES 4,5,6
C
         SIG456 = 1.0D0 - SIG123
         FC1 = (XST - VITDS)/(ABS(VITD-VITDS)+EPSU)
         SIG56 = 0.5D0*(SIGN(1.0D0,FC1)+1.0D0)
         FC2 = 0.5D0 * (FC1 + 1.0D0 - ABS(FC1 - 1.0D0))
         SIGCHO=0.5D0*(SIGN(1.0D0,VITDS-VITD+EPSU)+1.0D0)
C
         FPOIS3= (1.0D0 - SIGCHO)*(1.0D0 - FC2*SIG56)*SIG456 +
     &         SIGCHO * (SIG456 - SIG56)
         FPOIS4 = (1.0D0 - FPOIS3) * SIG456
C
C******* La solution pour le calcul du flux
C        (ro, T, w lineaire sur les  detantes
C         p = ro R T
C         h = sum_i ACVTO(i)/i T^i
C
C
         ROFLUX = ROG * FPOIS1 + ROGS * FPOIS2
     &          + RODS * FPOIS3 + ROD * FPOIS4
         TFLUX =  TG * FPOIS1 + TGS * FPOIS2
     &          + TDS * FPOIS3 + TD * FPOIS4
         UNFLUX = UNG * FPOIS1 + UGS * FPOIS2
     &          + UDS * FPOIS3 + UND * FPOIS4
C
         FPOISG = FPOIS1 + FPOIS2
         FPOISD = FPOIS3 + FPOIS4
         RTOTFL = RTOTG * FPOISG + RTOTD * FPOISD
         PFLUX  = ROFLUX * RTOTFL * TFLUX
         UTFLUX =  UTG * FPOISG + UTD * FPOISD
         ECFLUX = 0.5D0 * (UNFLUX * UNFLUX + UTFLUX * UTFLUX)
C
         FUNT = 1.0D0
         ETHFLU = 0.0D0
         DO I1 = 1, NORDP1, 1
            FUNT = FUNT * TFLUX
            ETHFLU = ETHFLU + (FPOISG * ACVTOG(I1) * FUNT / I1) +
     &         (FPOISD * ACVTOD(I1) * FUNT / I1)
         ENDDO
         HTFLUX = ETHFLU + PFLUX / ROFLUX
C
         ROU = ROFLUX * UNFLUX
         FLU(1) = ROU
         FLU(2) = ROU * UNFLUX + PFLUX
         FLU(3) = ROU * UTFLUX
         FLU(4) = ROU * (HTFLUX + ECFLUX)
         DO I1 = 1, NESP, 1
            FLU(I1+4) = ROU *
     &           (FPOISG * YG(I1) + FPOISD * YD(I1))
         ENDDO
         DO I1 = 1, NSCA, 1
            FLU(I1+4+NESP) = ROU *
     &           (FPOISG * SCAG(I1) + FPOISD * SCAD(I1))
         ENDDO
C
         CELLT = 1.0/MAX(ABS(VITG),ABS(VITD))
C
C******* Tests: les vitesses
C
C         write(3,'(A8,D16.8)') 'VITG   =', VITG
C         write(3,'(A8,D16.8)') 'VITGS  =', VITGS
C         write(3,'(A8,2D16.8)') 'VITS   =', UGS, UDS
C         write(3,'(A8,D16.8)') 'VITDS  =', VITDS
C         write(3,'(A8,D16.8)') 'VITD   =', VITD
C
C
C******* Test EF
C
C         write(4,'(5D14.6)') XST, FPOIS1, FPOIS2, FPOIS3, FPOIS4
C
      ENDIF
C
 9999 CONTINUE
      RETURN
      END