C BCOND2 SOURCE CB215821 16/04/21 21:15:17 8920 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=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,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