C BCOND SOURCE DURA 06/07/11 21:16:54 5481 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=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,*) ' 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