C BCOND2    SOURCE    DURA      06/07/11    21:16:48     5481           
      SUBROUTINE BCOND2 (X,DX,XL,RUG,XW,XN,TN,EN,BN,KIMP,PSLIM,REL,
     &         P1,PS1,T1,Y1,QAE,QEE1,PHI1,
     &         P2,T2,Y2,U2,QEE2,PHI2,QW2,RE,H,PSQ,RINDEX,
     &         NPP,ITP,PF,PP,DPF,DPP,RECU,XKUL,XKUT1,XKUT2,XKUT3,XKUT4)

      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C     operateur FUITE option gaz reel
C>>>  air + vapeur + liquide
C     QAE :  QA*E*B invariant air flowrate
C     QEEi : QEEi*E*B invariant water flowrate
C     XW=1.(0.) with (without) wall condensation
C     indice 1,2 : entree et sortie de troncon
C     RECU,XKUL,XKUT1,XKUT2,XKUT3,XKUT4 : coef lois de frot utilisateur
C     H correspond a l'energie totale cedee a la paroi
C     H=HM+RAPH*HW  transmis
C     H=HM+HW       pour le calcul de T2
C
      DIMENSION XN(NPP),TN(NPP),EN(NPP),BN(NPP)

      IF(KIMP.GE.1) THEN
      write(6,*) 'entree bcond2 X T1 Y1 PHI1= ',X,T1,Y1,PHI1
      ENDIF
      RINDEX=1.0

      ERPS1 = ABS (1.-(PS1/PSATT0(T1)))
      IF (ERPS1.GT.1D-4) THEN
       WRITE(6,*) ' bcond :erreur sur PS1 ',ERPS1
      ENDIF

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 (homogenes a ro.u)
C-------------------------------------
C      debit air
       QA=QAE/E1/B1
      
C      debit liquide
       QE1=QEE1/E1/B1
 
C      debit vapeur
       QV1=QE1*Y1
 
C      debit liquide
       QL1= (1-Y1)*QE1
 
C      debit total
       Q = QA + QE1

C-------------------------------------
C     DEFINITION DES MASSES VOLUMIQUES 
C     CONSTITUANT:
C                 : air
C                 : vapeur 
C                 : eau
C                 : eau + vapeur
C                 : air + eau + vapeur
C--------------------------------------
C     PS1=PHI1*P1

C     masse volumique de air
      ROA1=(P1-PS1)/RA/T1
C     masse volumique vapeur
      ROV1=ROVAP0(PS1,T1)
      
C     masse volumique eau liquide
      ROL1=ROLIQ0(P1,T1)
  
C     calcul la fraction volumique de gaz
      AL =BALPHA(Y1,ROL1,ROV1)
     
C     calcul la masse volumique air + eau + vapeur           
      RO1=AL *(ROA1+ROV1)+(1.-AL )*ROL1
C-------------------------------------------------
C     CALCUL DE LA VITESSE DE ECOULEMENT EN ENTREE
C------------------------------------------------- 
      
      U1=QA/AL /ROA1

C-------------------------------------------------
C     CALCUL DU COEFFICIENT DE COMPRESSIBILITE Z 
C     CONSTITUANT:
C                : vapeur
C-------------------------------------------------
      ZV = ZVAP0(ROV1,T1)
 
C-------------------------------------------------
C     CALCUL DES VISCOSITES DYNAMIQUES 
C     CONSTITUANT:
C                : air + vapeur
C                : air + eau + vapeur
C-------------------------------------------------

C      viscosite dynamique air+vapeur
       XMU=BMUG(T1,PS1/P1,T0)

C      calcul de la viscosite dynamique air+vapeur+eau_liquide
       XMU=AL*XMU+BMUL(T1,T0)*(1.-AL)
C         WRITE(6,*) 'XMU',XMU

C--------------------------------------------------
C     CALCUL DES ENTHALPIES SPECIFIQUES
C     CONSTITUANT:
C                : eau liquide
C                : vapeur
C--------------------------------------------------------         

C        calcul enthalpie specifique eau liquide
         HLS01 = HLS0(P1,T1)
C        WRITE(6,*) 'HL',HLS01
C        calcul enthalpie specifique de la vapeur
         HVS01 = HVS0(PS1,T1) 
C        WRITE(6,*) 'HV',HVS01

C--------------------------------------------------------
C     CALCUL CHALEUR LATENTE DU MELANGE
C     CONSTITUANT:
C                : eau+vapeur
C--------------------------------------------------------    

C        calcul la chaleur latente du melange
         XLAT = (HVS01-HLS01)
C        WRITE(6,*) 'XLAT',XLAT
C--------------------------------------------------------
C     CALCUL DES CONDUCTIVITES THERMIQUES
C     CONSTITUANT:
C                : air 
C--------------------------------------------------------  

         XLA = BLA(T1,T0) 
      
C--------------------------------------------------------
C     CALCUL DES CHALEURS SPECIFIQUES
C     CONSTITUANT:
C                : vapeur
C                : eau liquide 
C                : air + vapeur + eau
C--------------------------------------------------------         

C        calcul chaleur specifique de la vapeur
         CV1 = DHVDT0(PS1,T1)

C        calcul chaleur specifique eau liquide
         CL1 = CPLIQ0(P1,HLS01)

         
C        calcul chaleur specifique du melange air vapeur liquide          
         CP1 = (QA*CA+QE1*(Y1*CV1+(1-Y1)*CL1))/(QA+QE1)


C        WRITE(6,*) 'cv1,cl1,cp1',CV1,CL1,CP1
C--------------------------------------------------------
C     ACTUALISATION DES VARIABLES
C     CONSTITUANT:
C                : vitesse U
C                : titre Y
C                : densite du melange RO
C                : pression vapeur PS
C                : pression totale P
C                : temperature T
C                : debit eau (vapeu+liquide) QE
C                : debit total QX
C------------------------------------------------------ 

      E=E1
      B=B1
      U=U1
      Y=Y1
      RO=RO1
      PS=PS1
      P=P1
      T=T1
      QE=QE1
      QX=QA+QE

      
      
      IF ((X+DX).GT.1.D0) DX=1.D0-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
C       write(6,*)
C       write(6,*) ' niter P2 T2 ',NITER,P2,T2
         TT2=T2
         PP2=P2
         YY2=Y2


C--------------------------------------------------------
C     CALCUL DU NOMBRE DE REYNOLDS DE ECOULEMENT 
C     CONSTITUANT:
C                : air + eau + vapeur
C--------------------------------------------------------         
         
C        calcul du nombre de REYNOLDS    
         RE=QX*2*E/XMU
C--------------------------------------------------------
C     CALCUL DE LA PERTE CHARGE DU AU FORCE DE FROTTEMENT 
C     CONSTITUANT:
C                : air + eau + vapeur
C--------------------------------------------------------     
                 
         DPP = XL*QX*QX /4/E
 
C        calcul de la perte de charge DP 
         BK=BKFRO(RE,REL,XKL,XKT,2*E,RUG,RECU,XKUL,XKUT1,XKUT2,XKUT3,
     *   XKUT4)
         DP=DX*BK*DPP/RO
         IF(KIMP.GE.1) THEN
           WRITE(6,*) 'bcond2 ap BKFRO : RE,QX,DPP= ',RE,QX,DPP
           WRITE(6,*) 'bcond2 ap BKFRO : DX,BK,RO= ',DX,BK,RO
         ENDIF        
         P2=P1-DP
         IF (P2.LT.PSLIM) THEN
            PSQ = -1.
            RETURN
         ENDIF

         
C--------------------------------------------------------
C     CALCUL DU NOMBRE DE PRANDT DE ECOULEMENT 
C     CONSTITUANT:
C                : air + eau + vapeur
C--------------------------------------------------------                

C        calcul du nombre de Prandt estime en amont
         PR = XMU*CP1/BLA(T1,T0)
         
         
         
C--------------------------------------------------------
C     CALCUL DU COEFFICIENT ECHANGE THERMIQUE 
C     CONSTITUANT:
C                : air + eau + vapeur
C--------------------------------------------------------                
 
         
C        calcul du coefficient echange thermique         
         HM = BHECH(RE,REL,PR,XLA,E)
         H = HM
         HW=0.D0

C    si condensation en film
           IF(T1.GE.TP) THEN
           HW=XW*BWALL(XLAT,CV1,PS/P,H,(DPSAT0(T1)/PSATT0(T1)))
           H = HM+HW
           ENDIF

         
C        calcul de la fraction molaire de vapeur                
         PHI=PS/P
         
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(PS ,T )
         alpha2=DHVDP0(PS ,T )
         alpha3=DHLDP0(P ,T )
         HLIQ=HLS0(P ,T )
         alpha4=CPLIQ0(P ,HLIQ )
         alpha5=DPSAT0(T )
         alpha6=DZVDP0(PS ,T )
         alpha7=DZVDT0(PS ,T )
         
         
C--------------------------------------------------------
C     INTEGRATION DES PARAMETRES DE EQUATION THERMIQUE 
C           Qy dT/dx= (H)(Tw-T)+S  resolution BT3
C--------------------------------------------------------    


C        definition de Qc
         Qc=QAE*CA+Y *QEE1*(alpha1+(alpha2*alpha5))+(1-Y )*alpha4*QEE1
C        WRITE(6,*) 'QA QE1 Qc ',QA,QE1,Qc
         
         
C        definition de a         
         a=(RA/RV)*QAE
C        write(6,*) 'a ',a
         
         
C        DEFINITION DE ALPHA9    
         x1=(PHI *alpha5)/(ZV*PSATT0(T )*(1-PHI )*(1-PHI ))
         x2=(PHI *(alpha5*alpha6+alpha7))/(1-PHI )
         x2=x2/ZV/ZV
         alpha9=x1-x2

C        WRITE(6,*) 'ALPHA9 x1 x2 ',alpha9 ,x1,x2
        
C        DEFINITION DE ALPHA10   
         x3=(PHI *XLAT*a)/(ZV*P *(1-PHI )*(1-PHI ))
         x4=alpha3*(1-Y )*QEE1
         alpha10=x3-x4
C        WRITE(6,*) 'alpha10 x3 x4 ',alpha10,x3,x4
                
          
C        calcul du terme source de la thermique S
C        rappel DP=P1-P2 et non l'inverse
         S=-(alpha10)*(DP/(XL*DX))

         Qy=Qc+a*XLAT*alpha9
         
         
C----------------------------------------------------
C       Calcul la temperature T2 de sortie du troncon     
C----------------------------------------------------

         HEFF = 2 * H
C bt3 traite l echange de chaleur par unite d etendue
         Qy = Qy/B
         S = S/B
         T2 = BT3(T1,TP,S,HEFF,Qy,DX,XL)
 
         IF(KIMP.GE.1) THEN
         WRITE(*,*) 'bcond2 ap BT3: niter Qy B HEFF S',niter,Qy,B,HEFF,S
         WRITE(6,*) 'bcond2 ap BT3: T1 T2 TP= ',T1,T2,TP
         ENDIF
         
C        calcul de increment de temperature DT
         DT = T2-T1

C----------------------------------------------------
C        FIN DU TRONCON CALCULS DES NOUVELLES VALEURS    
C----------------------------------------------------

C        calcul de nouvelle pression vapeur saturante
         PS2=PSATT0(T2) 
         
C--------------------------------------------------------
C     DEFINITION DES  NOUVELLES MASSES VOLUMIQUES 
C     CONSTITUANT:
C                 : air
C                 : vapeur 
C                 : eau
C                 : eau + vapeur
C                 : air + eau + vapeur
C--------------------------------------------------------
 
C       masse volumique de air
        ROA2=(P2-PS2)/RA/T2

C       masse volumique vapeur
        ROV2=ROVAP0(PS2,T2)

C--------------------------------------------------------
C     CALCUL DU COEFFICIENT DE COMPRESSIBILITE Z de vapeur
C     CONSTITUANT:
C                : vapeur
C-------------------------------------------------------- 
        Z2 = ZVAP0(ROV2,T2)
        
C       masse volumique eau liquide
        ROL2=ROLIQ0(P2,T2)
      

C--------------------------------------------------------
C     Calcul du dQe 
C--------------------------------------------------------

C        calcul increment de debit eau
         DQEE= -2.*HW*DX*XL*B*(T-TP)/XLAT
         
C        calcul nouveau debit eau (liq+vap)
         QEE2=QEE1+DQEE
         
C        calcul debit eau
         QE2=QEE2/E2/B2

C        calcul du titre de vapeur sortie troncon          
         Y2=(ROV2/ROA2)*(QA/QE2)

         IF(KIMP.GE.3) THEN
         WRITE(6,*) ' ROA2 ROV2 Y2',ROA2,ROV2,Y2
         write(6,*) ' QE1 QE2  ',QE1,QE2
         ENDIF


C--------------------------------------------------------
C     On raffine le maillage dans 2 cas
C--------------------------------------------------------

         IF ( ((QEE2/QEE1).LT.1.E-5).OR.(Y2.GE.1.)) THEN
C         WRITE(6,*) ' WALL FLUX ',QEE2/QEE1
          RINDEX=0.5
          RETURN
         ENDIF


         QV2 = QE2 * Y2
         QL2 = QE2 * (1.-Y2)


C       calcul la fraction volumique de gaz
        AL2=BALPHA(Y2,ROL2,ROV2)
     
     
C       calcul la masse volumique air + eau + vapeur           
        RO2=AL2*(ROA2+ROV2)+(1.-AL2)*ROL2
         



         
C--------------------------------------------------------
C     CALCUL DE LA NOUVELLE VITESSE DE ECOULEMENT 
C-------------------------------------------------------- 
    
         U2= QA/AL2/ROA2

C--------------------------------------------------------
C     DEFINITION DES  NOUVELLES VALEURS MOYENNES 
C
C                 : Pression
C                 : Temperature
C                 : pression vapeur
C--------------------------------------------------------
 
         P=(P1+P2)/2
         T=(T1+T2)/2
         PS=(PS1+PS2)/2

         PHI = PS/P
         
C--------------------------------------------------------
C     CALCUL DES VITESSES MOYENNES SUR TRONCON
C     CONSTITUANT: melange
C                 : masse volumique
C                 : vitesse 
C                 : debit eau
C                 : debit total
C                 : titre
C--------------------------------------------------------
         
         
C        calcul de la masse volumique moyenne melange   
         RO=(RO1+RO2)/2
         
C        calcul de la vitesse moyenne du melange                 
         U=Q/RO
         
         
C        calcul du debit eau moyen du melange     
         QE=(QE1+QE2)/2.
         
C        calcul du debit total moyen du melange          
         QX=QA+QE

C        calcul du titre moyen sur le troncon    
         Y=(Y1+Y2)/2
         
C-----------------------------------------------------------------
C     CALCUL DU TAUX DE CONDENSATION EN FILM 
C-----------------------------------------------------------------

         QW2=(QEE2-QEE1)/DX/XL/B

         
         
         IF(KIMP.GE.3) write(6,*) ' P2 T2 Y2 ',P2,T2,Y2

C        definition des conditions fin de boucle  


C        critere sur P,T,Y

      ERT=DABS((T2-TT2)/T2)
      ERP=DABS((P2-PP2)/P2)
      ERY=DABS((Y2-YY2)/Y2)

      IF (((ERT.GT.1E-4).OR.(ERP.GT.1E-4).OR.(ERY.GT.1E-4))
     $  .AND.(NITER.LE.10))        GOTO 10

         IF(KIMP.GE.2) THEN
         write(6,2300) X, P2/P0, T2-T0, Y2, RO2, QE1, QE2
2300     FORMAT((1X,'bcond2: X P2 T2 Y2 RO2 QE1 QE2 ',7E12.5))
         ENDIF
         IF(KIMP.GE.3) THEN
         WRITE(6,*) ' HM HW HW/HM ',HM,HW,HW/HM
         ENDIF
         PHI2=PHI1


        
C--------------------------------------------------------
C     BILAN ENERGIE SUR TRONCON
C-------------------------------------------------------- 

C     write(6,*) ' E1 E2 ',E1,E2

      HH1=(QA*CA*T1)+(QV1*HVS0(PS1,T1))+(QL1*HLS0(P1,T1))
      HH1=HH1*E1*B1
      HH2=(QA*CA*T2)+(QV2*HVS0(PS2,T2))+(QL2*HLS0(P2,T2))
      HH2=HH2*E2*B2
C     puissance fluide
C !   PF = HH2-HH1
      PF = (HH2-HH1) - (DQEE*HVS0(PS1,T1))

C     puissance paroi
      PP = H *(T-TP)*2*DX*XL*B


      IF(KIMP.GE.2) THEN
      write(6,2200) X,HH1,HH2,PF,PP
2200  format(1X,'bcond2 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,'bcond2 X RE BKRO ',3E12.5)
      ENDIF
      DPF=PF/DX/XL/B
      DPP=PP/DX/XL/B

      RETURN
      END