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