Esope source files

   1 : C BCOND     SOURCE    DURA      06/07/11    21:16:54     5481           
   2 :       SUBROUTINE BCOND (X,DX,XL,RUG,XW,XN,TN,EN,BN,KIMP,PSLIM,REL,
   3 :      &         P1,PS1,T1,Y1,QAE,QEE1,PHI1,
   4 :      &         P2,T2,Y2,U2,QEE2,PHI2,QW2,RE,H,PSQ,RINDEX,
   5 :      &         NPP,ITP,PF,PP,DPF,DPP,RECU,XKUL,XKUT1,XKUT2,XKUT3,XKUT4)
   6 : 
   7 :       IMPLICIT INTEGER(I-N)
   8 :       IMPLICIT REAL*8 (A-H,O-Z)
   9 : C
  10 : C     operateur FUITE
  11 : C>>>  gaz + liquide
  12 : C
  13 : C     QAE :  QA*E invariant air flowrate
  14 : C     QEEi : QEEi*E invariant water flowrate
  15 : C     XW=1.(0.) with (without) wall condensation
  16 : C     indice 1,2 : entree et sortie de troncon
  17 : C     RECU,XKUL,XKUT1,XKUT2,XKUT3,XKUT4 : coef lois de frot utilisateur
  18 : C   H correspond a l'energie totale cedee a la paroi
  19 : C     H=HM+RAPH*HW  transmis
  20 : C     H=HM+HW       pour le calcul de T2
  21 : C
  22 :       DIMENSION XN(NPP),TN(NPP),EN(NPP),BN(NPP)
  23 : 
  24 :       RINDEX=1.0
  25 : 
  25 :       CALL BPHYS(T0,P0,RA,RV,CA,CV,CL,XLAT0,ROL,XKL,XKT,XREL)
  26 : 
  27 : C     temperature de paroi a l'entree
  27 :       CALL BTPAR(XN,TN,X,TP,NPP,ITP)
  28 :       TP=TP+T0
  29 : C     epaisseur a l'entree
  29 :       CALL BTPAR(XN,EN,X,E1,NPP,ITP)
  30 : C     etendue a l'entree
  30 :       CALL BTPAR(XN,BN,X,B1,NPP,ITP)
  31 : 
  32 :       IF(KIMP.GE.1) THEN
  33 :         WRITE(6,*) ' entree bcond X,E1,B1= ',X,E1,B1
  34 :       ENDIF
  35 : 
  36 :       QA1=QAE/E1/B1
  37 :       QE1=QEE1/E1/B1
  38 :       Q = QA1 + QE1
  39 : 
  40 :       ROV=PS1/RV/T1
  41 :       ROA=(P1-PS1)/RA/T1
  42 : C*
  43 :       AL1=BALPHA(Y1,ROL,ROV)
  44 :       RO1=AL1*(ROA+ROV)+(1.-AL1)*ROL
  45 : C*
  46 :       U1=QA1/AL1/ROA
  47 : 
  48 :       U=U1
  49 :       Y=Y1
  50 :       AL=AL1
  51 :       RO=RO1
  52 :       PS=PS1
  53 :       P=P1
  54 :       T=T1
  55 :       QE=QE1
  56 :       QA=QA1
  57 :       QX=QA+QE
  58 : 
  59 :       IF ((X+DX).GT.1.D0) DX=1.D0-X
  60 : 
  61 :       X=X+DX
  62 :       DU=0
  63 :       T2=0
  64 :       P2=0
  65 : c     ouverture en sortie de maille
  65 :       CALL BTPAR(XN,EN,X,E2,NPP,ITP)
  66 : c     etendue en sortie de maille
  66 :       CALL BTPAR(XN,BN,X,B2,NPP,ITP)
  67 :       QA2=QAE/E2/B2
  68 :       E = E1
  69 :       B = B1
  70 : 
  71 :       NITER = 0
  72 :   10  CONTINUE
  73 : 
  74 :       NITER = NITER + 1
  75 : c      write(6,*) 'bcond NITER P2 T2 ',NITER,P2,T2
  76 :          TT2=T2
  77 :          PP2=P2
  78 :          YY2=Y2
  79 : 
  80 :          XMU=BMUG(T,PS/P,T0)
  81 :          XMU=AL*XMU+BMUL(T,T0)*(1.-AL)
  82 :          RE=QX*2*E/XMU
  83 :          DPP = XL*QX *QX /4/E
  84 :          BK=BKFRO(RE,REL,XKL,XKT,2*E,RUG,RECU,XKUL,XKUT1,XKUT2,XKUT3,
  85 :      *   XKUT4)
  86 :          DP=DX*BK*DPP/RO
  87 : 
  88 :          P2=P1-DP
  89 :          IF (P2.LT.PSLIM) THEN
  90 :             PSQ = -1.
  91 :             RETURN
  92 :          ENDIF
  93 :          CP = (QA*CA+QE*(Y*CV+(1-Y)*CL))/QX
  94 :          PR=XMU*CP/BLA(T,T0)
  95 :          XLA=BLA(T,T0)
  96 : C*
  97 :          XLAT=XLAT0+(CV-CL)*T
  98 : C*
  99 :          FT=(-5/T+7000/T/T)
 100 :          HM = BHECH(RE,REL,PR,XLA,E)
 101 :          H = HM
 102 :          HW=0.D0
 103 :            IF(T1.GE.TP) THEN
 104 :            HW=XW*BWALL(XLAT,CV,PS/P,H,FT)
 105 :            H = HM+HW
 106 :            ENDIF
 107 : 
 108 :          PI=PS/P
 109 :          BB0 = (RA*QA/RV)*PI/(1.-PI)/(1-PI)
 110 :          BB = XLAT*BB0/QX/CP
 111 : 
 112 : C rappel DP=P1-P2 et non l'inverse
 113 : 
 114 :          SP=-BB*DP/P
 115 :          SU=-U*DU/CP
 116 :          S=SP
 117 : 
 118 :          BFT = BB * FT
 119 :          HBT=H/(1+BFT)
 120 :          S=S/(1+BFT)
 121 : 
 122 :          T2 = BT2(T1,TP,S,HBT,PR,RE,XLA,XL,DX)
 123 : 
 124 :          DT = T2-T1
 125 : 
 126 :          PS2=BPSAT(T2)
 127 :          ROV=PS2/RV/T2
 128 :          ROA=(P2-PS2)/RA/T2
 129 : 
 130 :          P=(P1+P2)/2
 131 :          T=(T1+T2)/2
 132 :          PS=(PS1+PS2)/2
 133 : 
 134 :          PSI = PS2/P2
 135 : 
 136 :          DQEE= -2.*HW*DX*XL*B*(T-TP)/XLAT
 137 :          QEE2=QEE1+DQEE
 138 :          QE2=QEE2/E2/B2
 139 :          QEEX=(QEE1+QEE2)/2.
 140 : 
 141 :          Y2=(RA*QA2/RV)*(PSI/(1.-PSI))/QE2
 142 : C
 143 : C        on raffine le maillage dans 2 cas
 144 : C
 145 :          IF ( ((QEE2/QEE1).LT.1.E-5).OR.(Y2.GE.1.)) THEN
 146 :           IF(KIMP.GE.1) THEN
 147 :            WRITE(6,*) 'bcond RINDEX05 X QEE2/QEE1 Y2',X,QEE2/QEE1,Y2
 148 :           ENDIF
 149 :           RINDEX=0.5
 150 : 
 151 :           RETURN
 152 :          ENDIF
 153 : 
 154 :          AL2=BALPHA(Y2,ROL,ROV)
 155 :          RO2=AL2*(ROA+ROV)+(1.-AL2)*ROL
 156 :          U2=(QAE/E2/B2)/AL2/ROA
 157 : 
 158 :          DU=U2-U1
 159 :          RO=(RO1+RO2)/2
 160 :          E=(E1+E2)/2
 161 :          QA=QAE/E/B
 162 :          QE=QEEX/E/B
 163 :          QX=QA+QE
 164 :          AL=(AL1+AL2)/2.
 165 : 
 166 :          DY=Y2-Y1
 167 :          Y=(Y1+Y2)/2
 168 : 
 169 :          QW2=(QEE2-QEE1)/DX/XL/B
 170 : 
 171 :          IF (KIMP.GE.1) THEN
 172 : c           write(6,*) 'bcond P2 T2 Y2 QA2 QE2',P2,T2,Y2,QA2,QE2
 173 :          ENDIF
 174 : 
 175 :       ERT=DABS((T2-TT2)/T2)
 176 :       ERP=DABS((P2-PP2)/P2)
 177 :       ERY=DABS((Y2-YY2)/Y2)
 178 : 
 179 :       IF (((ERT.GT.1E-4).OR.(ERP.GT.1E-4).OR.(ERY.GT.1E-4))
 180 :      $  .AND.(NITER.LE.10))        GOTO 10
 181 : 
 182 :          IF(KIMP.GE.2) THEN
 183 :          WRITE(6,*) ' HM HW HW/HM ',HM,HW,HW/HM
 184 :          ENDIF
 185 :          PHI2=PHI1
 186 : 
 187 : C*    bilan d'energie
 188 : 
 189 :       HH1=QAE*CA*T1+QEE1*(Y1*(CV*T1+XLAT0)+(1.-Y1)*CL*T1)
 190 :       HH2=QAE*CA*T2+QEE2*(Y2*(CV*T2+XLAT0)+(1.-Y2)*CL*T2)
 191 : C !   PF = HH2-HH1
 192 :       PF = (HH2-HH1) - DQEE*(CV*T1+XLAT0)
 193 : 
 194 :       RAPH=(CV*T+XLAT0)/XLAT
 195 : C pour l'energie cedee a la paroi
 196 :       H  =  HM+HW*RAPH
 197 :       PP = H *(T-TP)*2*DX*XL*B
 198 : 
 199 :       IF(KIMP.EQ.-3) THEN
 200 :       write(6,2100) X,RAPH,PF,PP,HH1,HH2
 201 : 2100  format(1X,'bcond X RH DH EC HH1 HH2 ',6E12.5)
 202 :       ENDIF
 203 :       IF(KIMP.NE.0.AND.X.EQ.1.) THEN
 204 :         write(6,2110) X,RE,BK
 205 : 2110    format(1X,'bcond X RE BKRO ',3E12.5)
 206 :       ENDIF
 207 : 
 208 :       DPF=PF/DX/XL/B
 209 :       DPP=PP/DX/XL/B
 210 : 
 211 :       RETURN
 212 :       END
 213 :  
 214 :  
 215 :  
 216 :  
 217 :  
 218 :  
 219 :