Sources Esope | Cast3M

bcond
C BCOND     SOURCE    CB215821  16/04/21    21:15:16     8920
      SUBROUTINE BCOND (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
C     operateur FUITE
C>>>  gaz + liquide
C
C     QAE :  QA*E invariant air flowrate
C     QEEi : QEEi*E 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)
 
      RINDEX=1.0
 
      CALL BPHYS(T0,P0,RA,RV,CA,CV,CL,XLAT0,ROL,XKL,XKT,XREL)
 
C     temperature de paroi a l'entree
      CALL BTPAR(XN,TN,X,TP,NPP,ITP)
      TP=TP+T0
C     epaisseur a l'entree
      CALL BTPAR(XN,EN,X,E1,NPP,ITP)
C     etendue a l'entree
      CALL BTPAR(XN,BN,X,B1,NPP,ITP)
 
      IF(KIMP.GE.1) THEN
        WRITE(6,*) ' entree bcond X,E1,B1= ',X,E1,B1
      ENDIF
 
      QA1=QAE/E1/B1
      QE1=QEE1/E1/B1
      Q = QA1 + QE1
 
      ROV=PS1/RV/T1
      ROA=(P1-PS1)/RA/T1
C*
      AL1=BALPHA(Y1,ROL,ROV)
      RO1=AL1*(ROA+ROV)+(1.-AL1)*ROL
C*
      U1=QA1/AL1/ROA
 
      U=U1
      Y=Y1
      AL=AL1
      RO=RO1
      PS=PS1
      P=P1
      T=T1
      QE=QE1
      QA=QA1
      QX=QA+QE
 
      IF ((X+DX).GT.1.D0) DX=1.D0-X
 
      X=X+DX
      DU=0
      T2=0
      P2=0
c     ouverture en sortie de maille
      CALL BTPAR(XN,EN,X,E2,NPP,ITP)
c     etendue en sortie de maille
      CALL BTPAR(XN,BN,X,B2,NPP,ITP)
      QA2=QAE/E2/B2
      E = E1
      B = B1
 
      NITER = 0
  10  CONTINUE
 
      NITER = NITER + 1
c      write(6,*) 'bcond NITER P2 T2 ',NITER,P2,T2
         TT2=T2
         PP2=P2
         YY2=Y2
 
         XMU=BMUG(T,PS/P,T0)
         XMU=AL*XMU+BMUL(T,T0)*(1.-AL)
         RE=QX*2*E/XMU
         DPP = XL*QX *QX /4/E
         BK=BKFRO(RE,REL,XKL,XKT,2*E,RUG,RECU,XKUL,XKUT1,XKUT2,XKUT3,
     *   XKUT4)
         DP=DX*BK*DPP/RO
 
         P2=P1-DP
         IF (P2.LT.PSLIM) THEN
            PSQ = -1.
            RETURN
         ENDIF
         CP = (QA*CA+QE*(Y*CV+(1-Y)*CL))/QX
         PR=XMU*CP/BLA(T,T0)
         XLA=BLA(T,T0)
C*
         XLAT=XLAT0+(CV-CL)*T
C*
         FT=(-5/T+7000/T/T)
         HM = BHECH(RE,REL,PR,XLA,E)
         H = HM
         HW=0.D0
           IF(T1.GE.TP) THEN
           HW=XW*BWALL(XLAT,CV,PS/P,H,FT)
           H = HM+HW
           ENDIF
 
         PI=PS/P
         BB0 = (RA*QA/RV)*PI/(1.-PI)/(1-PI)
         BB = XLAT*BB0/QX/CP
 
C rappel DP=P1-P2 et non l'inverse
 
         SP=-BB*DP/P
         SU=-U*DU/CP
         S=SP
 
         BFT = BB * FT
         HBT=H/(1+BFT)
         S=S/(1+BFT)
 
         T2 = BT2(T1,TP,S,HBT,PR,RE,XLA,XL,DX)
 
         DT = T2-T1
 
         PS2=BPSAT(T2)
         ROV=PS2/RV/T2
         ROA=(P2-PS2)/RA/T2
 
         P=(P1+P2)/2
         T=(T1+T2)/2
         PS=(PS1+PS2)/2
 
         PSI = PS2/P2
 
         DQEE= -2.*HW*DX*XL*B*(T-TP)/XLAT
         QEE2=QEE1+DQEE
         QE2=QEE2/E2/B2
         QEEX=(QEE1+QEE2)/2.
 
         Y2=(RA*QA2/RV)*(PSI/(1.-PSI))/QE2
C
C        on raffine le maillage dans 2 cas
C
         IF ( ((QEE2/QEE1).LT.1.E-5).OR.(Y2.GE.1.)) THEN
          IF(KIMP.GE.1) THEN
           WRITE(6,*) 'bcond RINDEX05 X QEE2/QEE1 Y2',X,QEE2/QEE1,Y2
          ENDIF
          RINDEX=0.5
 
          RETURN
         ENDIF
 
         AL2=BALPHA(Y2,ROL,ROV)
         RO2=AL2*(ROA+ROV)+(1.-AL2)*ROL
         U2=(QAE/E2/B2)/AL2/ROA
 
         DU=U2-U1
         RO=(RO1+RO2)/2
         E=(E1+E2)/2
         QA=QAE/E/B
         QE=QEEX/E/B
         QX=QA+QE
         AL=(AL1+AL2)/2.
 
         DY=Y2-Y1
         Y=(Y1+Y2)/2
 
         QW2=(QEE2-QEE1)/DX/XL/B
 
         IF (KIMP.GE.1) THEN
c           write(6,*) 'bcond P2 T2 Y2 QA2 QE2',P2,T2,Y2,QA2,QE2
         ENDIF
 
      ERT=ABS((T2-TT2)/T2)
      ERP=ABS((P2-PP2)/P2)
      ERY=ABS((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,*) ' HM HW HW/HM ',HM,HW,HW/HM
         ENDIF
         PHI2=PHI1
 
C*    bilan d'energie
 
      HH1=QAE*CA*T1+QEE1*(Y1*(CV*T1+XLAT0)+(1.-Y1)*CL*T1)
      HH2=QAE*CA*T2+QEE2*(Y2*(CV*T2+XLAT0)+(1.-Y2)*CL*T2)
C !   PF = HH2-HH1
      PF = (HH2-HH1) - DQEE*(CV*T1+XLAT0)
 
      RAPH=(CV*T+XLAT0)/XLAT
C pour l'energie cedee a la paroi
      H  =  HM+HW*RAPH
      PP = H *(T-TP)*2*DX*XL*B
 
      IF(KIMP.EQ.-3) THEN
      write(6,2100) X,RAPH,PF,PP,HH1,HH2
2100  format(1X,'bcond X RH DH EC HH1 HH2 ',6E12.5)
      ENDIF
      IF(KIMP.NE.0.AND.X.EQ.1.) THEN
        write(6,2110) X,RE,BK
2110    format(1X,'bcond X RE BKRO ',3E12.5)
      ENDIF
 
      DPF=PF/DX/XL/B
      DPP=PP/DX/XL/B
 
      RETURN
      END