Sources Esope | Cast3M

bsur2
C BSUR2     SOURCE    CB215821  16/04/21    21:15:25     8920
      SUBROUTINE BSUR2 (X,DX,XL,RUG,XW,XN,TN,EN,BN,KIMP,PSLIM,REL,
     &      RINDEX,P1,T1,QAE,QEE1,PHI1,P2,T2,U2,QEE2,PHI2,QW2,RE,H,
     &      PSQ,NPP,ITP,PF,PP,DPF,DPP,RECU,XKUL,XKUT1,XKUT2,XKUT3,XKUT4)
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION XN(NPP),TN(NPP),EN(NPP),BN(NPP)
 
C     operateur FUITE
C>>>  superheated vapour
C     version with wall condensation
C     version gaz reel
C     QAE:  QA.E.B  constant air flowrate
C     QEEi: QEi.E.B water flowrate
C     indice 1,2 : entree et sortie de troncon
C     RECU,XKUL,XKUT1,XKUT2,XKUT3,XKUT4 : coef lois de frot utilisateur
C
C     H correspond au flux d'energie total cede a la paroi
C     HM (convection) est utilise pour le calcul de T
C
      IF(KIMP.GE.1) THEN
        write(6,*) 'entree bsur2'
      ENDIF
 
      HW=0.D0
      RINDEX=1.0
      DQEE = 0.D0
 
C     IMPLEMENTATION DES PROPRIETES PHYSIQUES
      CALL BPHYS(T0,P0,RA,RV,CA,XCV,XCL,XLAT0,XROL,XKL,XKT,XREL)
 
C     TEMPERATURE PAROI A ENTREE (KELVIN)
      CALL BTPAR(XN,TN,X,TP,NPP,ITP)
      TP=TP+T0
 
C     EPAISSEUR DE FISSURE EN ENTREE
      CALL BTPAR(XN,EN,X,E1,NPP,ITP)
C     etendue a l'entree
      CALL BTPAR(XN,BN,X,B1,NPP,ITP)
 
C--------------------------------------------------------
C     DEFINITION DES DIFFERENTS DEBITS
C--------------------------------------------------------
C      debit air
       QA=QAE/E1/B1
 
C      debit eau (vapeur+liquide)
       QE1=QEE1/E1/B1
 
C      debit total
       Q=QA+QE1
C
C     DEFINITION DES MASSES VOLUMIQUES
C     CONSTITUANT:
C                 : air
C                 : vapeur
 
      PV1=PHI1*P1
 
C     masse volumique de air
      ROA1=(P1-PV1)/RA/T1
C     WRITE(6,*) 'T1,P1,PV1,RA',T1,P1,PV1,RA
 
C     masse volumique vapeur
      ROV1=ROVAP0(PV1,T1)
 
C     masse volumique du melange
      RO1 = ROA1 + ROV1
C     write(6,*) ' ROA ROV RO ',ROA1,ROV1,RO1
 
C--------------------------------------------------------
C     CALCUL DE LA VITESSE DE ECOULEMENT EN ENTREE
C--------------------------------------------------------
 
C     calcul de la vitesse
      U1=QA/ROA1
C     write(6,*) ' U1 ',U1
 
 
C--------------------------------------------------------
C     CALCUL DU COEFFICIENT DE COMPRESSIBILITE Z de vapeur
C     CONSTITUANT:
C                : vapeur
C--------------------------------------------------------
C     calcul du ceofficient de compressibilite Z vapeur
      ZV   = ZVAP0(ROV1,T1)
C     WRITE(6,*) 'Z',ZV
 
 
C--------------------------------------------------------
C     CALCUL DES CONDUCTIVITES THERMIQUES
C     CONSTITUANT:
C                : air
C--------------------------------------------------------
 
C        calcul la conductivite de air
         XLA = BLA(T1,T0)
 
 
 
C--------------------------------------------------------
C     CALCUL DES CHALEURS SPECIFIQUES
C     CONSTITUANT:
C                : vapeur
C----------------------------------------------------
 
         CV1 = DHVDT0(PV1,T1)
C        WRITE(6,*) 'cv1',CV1
 
 
 
C--------------------------------------------------------
C     CALCUL DES VISCOSITES DYNAMIQUES
C     CONSTITUANT:
C                : air + vapeur
C--------------------------------------------------------
 
       XMU=BMUG(T1,PHI1,T0)
C      WRITE(6,*) 'xmu',XMU
 
 
C--------------------------------------------------------
C     CALCUL DU NOMBRE DE PRANDT DE ECOULEMENT
C     CONSTITUANT:
C                : air + vapeur
C--------------------------------------------------------
 
        PR=XMU*CV1/BLA(T1,T0)
C       WRITE(6,*) 'PR',PR
C--------------------------------------------------------
C     CALCUL DES ENTHALPIES SPECIFIQUES
C     CONSTITUANT:
C                : eau liquide
C                : vapeur
C--------------------------------------------------------
         HVS01 = HVS0(PV1,T1)
         HLS01 = HLS0(P1,T1)
C        WRITE(6,*) 'HV',HVS01
C        WRITE(6,*) 'HL',HLS01
 
C--------------------------------------------------------
C     CALCUL CHALEUR LATENTE DU MELANGE
C            (utile pour la condensation en film)
C     CONSTITUANT:
C                : eau+vapeur
C--------------------------------------------------------
 
         XLAT = (HVS01-HLS01)
C        WRITE(6,*) 'L',XLAT
 
C--------------------------------------------------------
C     ACTUALISATION DES VARIABLES
C     CONSTITUANT:
C                : vitesse U
C                : pression totale P
C                : temperature T
C----------------------------------------------------
 
      E=E1
      B=B1
      U=U1
      T=T1
      P=P1
      RO=RO1
 
 
      IF ((X+DX).GT.1.) DX=1-X
      X=X+DX
 
      DU=0
      T2=0
      P2=0
 
      CALL BTPAR(XN,EN,X,E2,NPP,ITP)
      CALL BTPAR(XN,BN,X,B2,NPP,ITP)
 
      NITER=0
 
10    CONTINUE
 
      NITER=NITER+1
      IF(KIMP.GE.1) THEN
        write(*,*) 'NITER= ',NITER
      ENDIF
 
      TT2=T2
      PP2=P2
 
C--------------------------------------------------------
C     CALCUL DU NOMBRE DE REYNOLDS
C     CONSTITUANT:
C                : melange
C-------------------------------------------------------
       RE=Q*2*E/XMU
      IF(KIMP.GE.1) THEN
        WRITE(6,*) 'bsur2: RE',RE
        WRITE(6,*) 'bsur2: Q,E,B,XMU',Q,E,B,XMU
      ENDIF
C--------------------------------------------------------
C     CALCUL DE LA PERTE DE CHARGE DP DUE AU FROTTEMENT
C     CONSTITUANT:
C                : Melange
C--------------------------------------------------------
 
      DPP = XL*Q*Q/4/E
      BK = BKFRO(RE,REL,XKL,XKT,2*E,RUG,RECU,XKUL,XKUT1,XKUT2,XKUT3,
     *XKUT4)
      DP=DX*BK*DPP/RO
C     WRITE(6,*) 'XL Q E ', XL,Q,E
C     WRITE(6,*) 'DX RO DP', DX,RO,DP
 
      P2=P1-DP
         IF (P2.LT.PSLIM) THEN
            PSQ = -1.
            RETURN
         ENDIF
 
 
C--------------------------------------------------------
C     CALCUL DU TERME SOURCE S POUR LE CALCUL DE T
C     CONSTITUANT:
C                : melange
C----------------------------------------------------
 
      S=0.D0
 
C--------------------------------------------------------
C     CALCUL DU COEFFICIENT ECHANGE THERMIQUE H
C     CONSTITUANT:
C                : Melange
C----------------------------------------------------
 
      HM = BHECH(RE,REL,PR,XLA,E)
C     write(6,*) RE,REL,PR
C     write(6,*) XLA,E
C     write(6,*) ' HM ',HM
 
C if vapour at the inlet
 
      IF(PHI1.GT.1.D-2) THEN
C
C wall condensation
         D=2.5D-5
         A=2.3D-5
 
C--------------------------------------------------------
C     CALCUL PRESSION DE SATURATION
C     CONSTITUANT:
C                : vapeur
C----------------------------------------------------
 
         PP=PSATT0(TP)
C        WRITE(6,*) 'TP',TP
C--------------------------------------------------------
C     CALCUL PRESSION PARTIELLE DE VAPEUR
C     CONSTITUANT:
C                : vapeur
C----------------------------------------------------
C        WRITE(6,*) 'ROV1',ROV1
         PV=PVAP0(ROV1,T1)
C        WRITE(6,*) 'PV',PV
 
 
C--------------------------------------------------------
C       CONDENSATION EN FILM
C--------------------------------------------------------
 
 
         RAPP=(P1-PP)/(P1-PV)
 
         IF (RAPP.GT.(1.0001)) THEN
         XJ=XW*(HM/CV1)*(D/A)**(2./3.)*LOG(RAPP)
         ELSE
         XJ=0.D0
         ENDIF
 
         DQEE=-2.*XJ*DX*XL*B
 
         SW=(XLAT/Q/CV1/E)*DQEE
 
         DELTAT = T-TP
         IF (ABS(DELTAT).GT.1D-3) THEN
         HW=XJ*XLAT/(T-TP)
C        write(6,*) ' T TP HW ',T,TP,HW
         ENDIF
 
         QEE2=QEE1+DQEE
         QE2=QEE2/E2/B2
C        write(6,*) ' QE1 QE2  ',QE1,QE2
 
         IF ((QEE2/QEE1).LT.(-1.D-5)) THEN
C          write(6,*) ' relative water flowrate ',QE2/QE1
           RINDEX = 0.5
           RETURN
C          PHI2=0.D0
C          QEE2=0.D0
C          QE2=0.D0
         ELSE
           A=RA*QA/RV/(QE2*ZV)
           PHI2=1./(1.+A)
         ENDIF
 
 
* no vapour at the inlet
       ELSE
         XLAT=(HVS01-HLS01)
         HW=0.D0
         PHI2=0.D0
         QEE2=0.D0
         QE2=0.D0
       ENDIF
 
 
C      write(6,*) ' QE ',QE1,QE2
C--------------------------------------------------------
C     CALCUL DES COEFFICIENTS ALPHAi
C     CONSTITUANT:
C                : on se sert des fonctions de VARI
C--------------------------------------------------------
 
 
C        definition des alphai pour resolution systeme
         alpha1=DHVDT0(PV1,T1)
         alpha2=DHVDP0(PV1,T1)
 
         alpha6=DZVDP0(PV1,T1)
         alpha7=DZVDT0(PV1,T1)
         alpha8=DRVDP0(PV1,T1)
         alpha9=DRVDT0(PV1,T1)
 
C         WRITE(6,*) 'ALPHA',alpha1
C         WRITE(6,*) 'ALPHA',alpha2
C         WRITE(6,*) 'ALPHA',alpha3
C         WRITE(6,*) 'ALPHA',alpha4
C         WRITE(6,*) 'ALPHA',alpha5
C         WRITE(6,*) 'ALPHA',alpha6
C         WRITE(6,*) 'ALPHA',alpha7
C         WRITE(6,*) 'ALPHA',alpha8
C         WRITE(6,*) 'ALPHA',alpha9
 
C--------------------------------------------------------
C     DEFINITION DU Qxx,Qy
C     CONSTITUANT:
C                : on se sert des fonctions de VARI
C--------------------------------------------------------
 
        Qxx= QAE*CA+QEE1*alpha1
        s1=ZV*((alpha9*T1)+ROV1*T1*alpha7)
        s2= RA*T1*(alpha8*ZV+ROV1*alpha6)
        s3= 1-s2
        s4=s1/s3
        Qy= Qxx + QEE1*alpha2*RA*s4
 
C--------------------------------------------------------
C     CALCUL DE LA TEMPERATURE T2 EN SORTIE DU TRONCON
C----------------------------------------------------
 
      HEFF = 2 * HM
C bt3 traite l echange de chaleur par unite de largeur
      Qy = Qy/B
      T2 = BT3(T1,TP,S,HEFF,Qy,DX,XL)
      IF(KIMP.GE.1) THEN
        WRITE(6,*) 'bsur2 : T2= ',T2
      ENDIF
 
C--------------------------------------------------------
C     CALCUL DE LA VARIATION DE TEMPERATURE SUR TRONCON
C--------------------------------------------------------
 
 
      DT = T2-T1
 
C--------------------------------------------------------
C     CALCUL DE LA VITESSE DE ECOULEMENT U2 SORTIE TRONCON
C---------------------------------------------------------
 
C     WRITE(6,*) 'P2 RA T2 ',P2,RA,T2
 
C     initialisation avec la valeur a l'entree du troncon
 
      ZTRAN = ZV
      NZ = 0
 20   CONTINUE
      NZ = NZ + 1
 
      AA = (QEE2/QAE)* (RV*ZTRAN/RA)
      PHI2 = AA/(1.+AA)
      PV2 = PHI2 * P2
      ROV2 = ROVAP0(PV2,T2)
      Z2  = ZVAP0(ROV2,T2)
C     write(6,*) ' AA PHI2 ROV2 ',AA,PHI2,ROV2
      DELTAZ = ABS(Z2 - ZTRAN)
C     WRITE(6,*) ' NZ ZVAP ', NZ,Z2
      ZTRAN=Z2
      IF ((DELTAZ.GT.1D-5).AND.(NZ.LE.20)) GOTO 20
 
      IF (KIMP.GE.2) WRITE(6,*) ' NZ ZVAP ', NZ,Z2
 
      ROA2 = (1.-PHI2)*P2/RA/T2
      RO2 = ROA2 + ROV2
 
      U2 = QA/ROA2
 
C     WRITE(6,*) ' ROA ROV RO Z PV ',ROA2,ROV2,RO2,Z2,PV2
C     WRITE(6,*) ' PHI U ',PHI2,U2
 
 
      IF(KIMP.GE.2) THEN
      write(6,1200)  X,(P2/P0),(T2-T0),PV2,Z2
1200  FORMAT(1X,' bsur2 X P2 T2 PV2 Z2 ',5E12.5)
      write(6,1300)  QE1,QE2
1300  FORMAT(1X,'       QE1 QE2 ',2E12.5)
      ENDIF
 
 
C---------------------------------------------------------
C     CALCUL DES VALEURS MOYENNES DE P,T,U,RO
C---------------------------------------------------------
 
      T=(T1+T2)/2
      P=(P1+P2)/2
      RO= (RO1+RO2)/2
C     write(6,*) ' RO1 RO2 ',RO1,RO2
      U=Q/RO
 
C     write(6,*) ' T P RO U ',T,P,RO,U
 
      QE = (QE1 + QE2)/2
      Q = QA + QE
 
      QW2=(QEE2-QEE1)/DX/XL/B
 
      ERT=ABS((T2-TT2)/T2)
      ERP=ABS((P2-PP2)/P2)
 
       IF(((ERT.GT.1E-4).OR.(ERP.GT.1E-4))
     $  .AND.(NITER.LE.10))    GOTO 10
 
 
         IF(KIMP.GE.2) THEN
         write(6,*) ' HM HW HW/HM ',HM,HW,HW/HM
         ENDIF
 
C     calcul des puissances echangees
 
C     RAPH=(CV1*T+XLAT0)/XLAT
C     H=HM+HW*RAPH
      H = HM
 
 
C-------------------------------------------
C   BILAN ENERGIE TOTALE
C-------------------------------------------
 
      HH1=(QAE*CA*T1) + (QEE1*HVS0(PV1,T1))
      HH2=(QAE*CA*T2) + (QEE2*HVS0(PV2,T2))
C     puissance fluide
C !   PF = HH2-HH1
      PF = (HH2-HH1) - (DQEE*HVS0(PV1,T1))
 
C     puissance paroi
      PP = H*(T-TP)*2*DX*XL*B
 
 
      IF(KIMP.GE.2) THEN
        write(6,1000)  X,HH1,HH2,PF,PP
1000    FORMAT(1X,' bsur2 X HH1 HH2 DH HDT ',5E12.5)
      ENDIF
      IF(KIMP.NE.0.AND.X.EQ.1.) THEN
        write(6,2110) X,RE,BK
2110    format(1X,'bsur2 X RE BKRO ',3E12.5)
      ENDIF
      DPF=PF/DX/XL/B
      DPP=PP/DX/XL/B
 
      RETURN
      END