diff options
Diffstat (limited to 'Donjon/src/THMH2O.f')
| -rw-r--r-- | Donjon/src/THMH2O.f | 389 |
1 files changed, 389 insertions, 0 deletions
diff --git a/Donjon/src/THMH2O.f b/Donjon/src/THMH2O.f new file mode 100644 index 0000000..eb22012 --- /dev/null +++ b/Donjon/src/THMH2O.f @@ -0,0 +1,389 @@ +*DECK THMH2O + SUBROUTINE THMH2O(ITIME,I,J,K,K0,PINLET,MFLOW,HMAVG,ENT,HD,IFLUID, + > IHCONV,KHCONV,ISUBM,RADCL,ZF,VCOOL,IDFM,PHI,XFL,EPS,SLIP, + > ACOOL,PCH,DZ,TCALO,RHO,RHOLAV,RHOG,TSCLAD,KWA,VGJprime,HLV) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Nucleate boiling correlations along a single coolant channel. +* +*Copyright: +* Copyright (C) 2012 Ecole Polytechnique de Montreal. +* +*Author(s): +* A. Hebert, P. Gallet +* +*Parameters: input +* ITIME type of calculation (0=steady-state; 1=transient). +* I position of channel alon X-axis +* J position of channel alon Y-axis +* K position of channel alon Z-axis +* K0 onser of nuclear boiling point +* PINLET pressure in Pascal +* MFLOW massic coolant flow rate in Kg/m^2/s +* HMAVG averaged enthalpy +* ENT four values of enthalpy in J/Kg to be used in Gaussian +* integration +* HD hydraulic diameter in m +* IFLUID type of fluid (0=H2O; 1=D2O). +* IHCONV flag indicating HCONV chosen (0=default/1=user-provided). +* KHCONV fixed user-provided HCONV value in W/m^2/K. +* ISUBM subcooling model (0: one-phase; 1: Jens-Lottes model; +* 2: Saha-Zuber model). +* RADCL outer clad radius in m +* ZF parameters used to compute heat flux on clad surface in +* transient cases. +* PHI heat flow exchanged between clad and fluid in W/m^2. +* Given in steady-state cases. +* XFL input coolant flow quality +* EPS input coolant void fraction +* SLIP input slip ratio of vapor phase speed to liquid phase speed. +* ACOOL coolant cross section area in m^2. +* PCH heating perimeter in m. +* DZ axial mesh width in m. +* VCOOL local coolant velocity +* IDFM flag indicating if the drift flux model is to be used +* (0=HEM1(no drift velocity)/1=EPRI/2=MODEBSTION/3=GERAMP/4=CHEXAL) +* +*Parameters: output +* PHI heat flow exchanged between clad and fluid in W/m^2. +* Computed in transient cases. +* XFL output coolant flow quality +* EPS output coolant void fraction +* SLIP output slip ratio of vapor phase speed to liquid phase speed. +* TCALO coolant temperature in K +* RHO coolant density in Kg/m^3 +* RHOLAV liquid density in kg/m^3 +* RHOG vapour density in kg/m^3 +* TSCLAD clad temperature in K +* KWA flow regime (=0: single-phase; =1: subcooled; =2: nucleate +* boiling; =3 superheated steam) +* VGJ drift velocity in m/s +* VGJprime +* HLV delta between liquid and vapour enthalpies +* +*----------------------------------------------------------------------- +* +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER I,J,K,K0,IFLUID,IHCONV,ISUBM,KWA + REAL PINLET,MFLOW,HMAVG,ENT(4),HD,KHCONV,RADCL,ZF(2),PHI,TCALO, + > RHO,RHOLAV,TSCLAD,XFL,EPS,SLIP,ACOOL,PCH,DZ,VCOOL,VGJprime +*---- +* LOCAL VARIABLES +*---- + REAL W(4),HL(4),JL,JG,REL,PRL,VGJ,UL + CHARACTER HSMG*131 + LOGICAL LFIRST +*---- +* SAVE VARIABLES +*---- + SAVE DHSUB,DSAT,W + DATA W /0.347855,0.652145,0.652145,0.347855/ +*---- +* COMPUTE THE PROPERTIES OF THE SATURATED STEAM +*---- + IF(HMAVG.LT.0.0) CALL XABORT('THMH2O: NEGATIVE INPUT ENTHALPY.') + IF(IFLUID.EQ.0) THEN + CALL THMSAT(PINLET,TSAT) + CALL THMTX(TSAT,0.0,RHOL,HLSAT,ZKL,ZMUL,CPL) + CALL THMTX(TSAT,1.0,RHOG,HGSAT,ZKG,ZMUG,CPG) + ELSE IF(IFLUID.EQ.1) THEN + CALL THMHST(PINLET,TSAT) + CALL THMHTX(TSAT,0.0,RHOL,HLSAT,ZKL,ZMUL,CPL) + CALL THMHTX(TSAT,1.0,RHOG,HGSAT,ZKG,ZMUG,CPG) + ENDIF +*---- +* COMPUTE THE DENSITY AND TEMPERATURE OF THE LIQUID +*---- + HL(1)=MIN1(ENT(1),HLSAT) + HL(2)=MIN1(ENT(2),HLSAT) + HL(3)=MIN1(ENT(3),HLSAT) + HL(4)=MIN1(ENT(4),HLSAT) + CALL THMPH(IFLUID,PINLET,HL(1),R11,TL1) + CALL THMPH(IFLUID,PINLET,HL(2),R11,TL2) + CALL THMPH(IFLUID,PINLET,HL(3),R11,TL3) + CALL THMPH(IFLUID,PINLET,HL(4),R11,TL4) + IF(IFLUID.EQ.0) THEN + CALL THMPT(PINLET,TL1,RHO1,R2,R3,R4,CP1) + CALL THMPT(PINLET,TL2,RHO2,R2,R3,R4,CP2) + CALL THMPT(PINLET,TL3,RHO3,R2,R3,R4,CP3) + CALL THMPT(PINLET,TL4,RHO4,R2,R3,R4,CP4) + IF(ABS(TSAT-TL1).LT.0.1) CALL THMTX(TSAT,0.0,RHO1,R2,R3,R4,CP1) + IF(ABS(TSAT-TL2).LT.0.1) CALL THMTX(TSAT,0.0,RHO2,R2,R3,R4,CP2) + IF(ABS(TSAT-TL3).LT.0.1) CALL THMTX(TSAT,0.0,RHO3,R2,R3,R4,CP3) + IF(ABS(TSAT-TL4).LT.0.1) CALL THMTX(TSAT,0.0,RHO4,R2,R3,R4,CP4) + ELSE IF(IFLUID.EQ.1) THEN + CALL THMHPT(PINLET,TL1,RHO1,R2,R3,R4,CP1) + CALL THMHPT(PINLET,TL2,RHO2,R2,R3,R4,CP2) + CALL THMHPT(PINLET,TL3,RHO3,R2,R3,R4,CP3) + CALL THMHPT(PINLET,TL4,RHO4,R2,R3,R4,CP4) + IF(ABS(TSAT-TL1).LT.0.1) CALL THMHTX(TSAT,0.0,RHO1,R2,R3,R4,CP1) + IF(ABS(TSAT-TL2).LT.0.1) CALL THMHTX(TSAT,0.0,RHO2,R2,R3,R4,CP2) + IF(ABS(TSAT-TL3).LT.0.1) CALL THMHTX(TSAT,0.0,RHO3,R2,R3,R4,CP3) + IF(ABS(TSAT-TL4).LT.0.1) CALL THMHTX(TSAT,0.0,RHO4,R2,R3,R4,CP4) + ENDIF + TL=0.5*(W(1)*TL1+W(2)*TL2+W(3)*TL3+W(4)*TL4) + RHOLAV=0.5*(W(1)*RHO1+W(2)*RHO2+W(3)*RHO3+W(4)*RHO4) + CPLAV=0.5*(W(1)*CP1+W(2)*CP2+W(3)*CP3+W(4)*CP4) +*---- +* COMPUTE THE STEAM FLOW QUALITY AND LIQUID ENTHALPY +* Reference: R. T. Lahey Jr. and F. J. Moody, "The thermal hydraulics +* of a Boiling water nuclear reactor," American Nuclear Society, 1977. +* Equation (5.177), page 224 +* F2: Thermodynamic quality +*---- + TSCLAD=600.0 + IF(K0.GT.0) TSCLAD=TSAT+DSAT + XFL0=XFL + EPS0=EPS + SLIP0=SLIP + LFIRST=.TRUE. + HLAVG=HMAVG + F2=0.0 + F3=0.0 + IF(K0.GT.0) THEN + HLV=HGSAT-HLSAT + IF((HLV.GT.0.0).AND.(DHSUB.GT.0.0)) THEN + F2=(HMAVG-HLSAT)/HLV + F3=(DHSUB/HLV)*EXP(-(HMAVG-HLSAT)/DHSUB-1.0) + ENDIF + IF(HMAVG.GE.HGSAT) THEN + XFL=1.0 + EPS=1.0 + SLIP=1.0 + HLAVG=0.0 + ELSE + IF(ISUBM.EQ.1) THEN +* Use the Paul Gallet thesis model. + PI=RHOLAV*CPLAV*(TSCLAD-TL)/(RHOG*HLV) + XFL=XFL0+PCH*PHI*DZ/(MFLOW*ACOOL*HLV)/(1.0+PI) + ELSE IF(ISUBM.EQ.2) THEN +* Use a profile fit model. + XFL=MAX(XFL0,(F2+F3)/(1.0+F3)) + ENDIF + HLAVG=MIN(HLSAT,(HMAVG-XFL*HGSAT)/(1.0-XFL)) + ENDIF +*---- +* RECOMPUTE THE LIQUID PROPERTIES +*---- + IF(HLAVG.GT.0.0) THEN + CALL THMPH(IFLUID,PINLET,HLAVG,RHOL,TL) + IF(IFLUID.EQ.0) THEN + CALL THMPT(PINLET,TL,R1,R2,ZKL,ZMUL,CPL) + ELSE IF(IFLUID.EQ.1) THEN + CALL THMHPT(PINLET,TL,R1,R2,ZKL,ZMUL,CPL) + ENDIF +*---- +* COMPUTE THE COOLANT VOID FRACTION AND SLIP RATIO +* A drift-flux model is proposed by means of the concentration +* parameter CO and the drift velocity VGJ (their correspondent +* correlations are supposed to work fine under different flow regimes. +*---- + IF(HGSAT.GT.HLSAT) THEN + CO=1.13 + PR=PINLET/10**6 + SIGM=-7.2391E-6*PR**3+2.8345E-4*PR**2-5.1566E-3*PR+4.2324E-2 + VGJ=1.18*((SIGM*9.81*(RHOL-RHOG))/RHOL**2)**0.25 + F4=CO*((XFL*RHOL)+((1.0-XFL)*RHOG))+(RHOL*RHOG*VGJ/MFLOW) + EPS=(XFL*RHOL)/F4 + JL=(1.0-XFL)*MFLOW/RHOL + JG=XFL*MFLOW/RHOG + IF(EPS.NE.0) SLIP=JG*(1.0-EPS)/(JL*EPS) + ENDIF + ELSE +* superheated steam + CALL THMPH(IFLUID,PINLET,HMAVG,RHOG,TCALO) + IF(IFLUID.EQ.0) THEN + CALL THMPT(PINLET,TCALO,R1,R2,ZKG,ZMUG,CPG) + ELSE IF(IFLUID.EQ.1) THEN + CALL THMHPT(PINLET,TCALO,R1,R2,ZKG,ZMUG,CPG) + ENDIF + ENDIF + ENDIF +*---- +* COMPUTE THE FLUID PROPERTIES +* RHO: fluid density +* REL: Reynolds number of liquid phase +* PRL: Prandtl number of liquid phase +*---- + IF(XFL.EQ.0.0) THEN +* One phase liquid + TB=TSAT-0.1 + IF(TL.LT.TB) THEN + TCALO=TL + ELSE + TCALO=TB + ENDIF + IF(IFLUID.EQ.0) THEN + CALL THMPT(PINLET,TCALO,R1,R2,ZKONE,ZMUONE,CPONE) + ELSE IF(IFLUID.EQ.1) THEN + CALL THMHPT(PINLET,TCALO,R1,R2,ZKONE,ZMUONE,CPONE) + ENDIF + RHO=RHOLAV + REL=MFLOW*HD/ZMUONE + PRL=ZMUONE*CPONE/ZKONE + ELSE IF(HMAVG.LT.HGSAT) THEN +* Two-phase flow + IF((IFLUID.EQ.0).AND.(IDFM.GT.0)) THEN + CALL THMDFM(PINLET,VCOOL,HMAVG,HD,TL,TSAT,IDFM,EPS,XFL, + > RHO,RHOL,RHOG,VGJ,VGJprime,C0,HLV) + CALL THMTX(TCALO, 0.0, RHO11, H11, ZK11, ZMUL, CPL) + UL = VCOOL - (EPS / (1.0 - EPS))*RHOG/RHO * VGJprime + REL = ABS(UL*RHOL) * HD / ZMUL + PRL = ZMUL*CPL/ZKL + ELSE + REL=MFLOW*(1.0-XFL)*HD/ZMUL + PRL=ZMUL*CPL/ZKL + ENDIF + TCALO=EPS*TSAT+(1.0-EPS)*TL + ZKONE=ZKL + CPONE=CPL + RHO=EPS*RHOG+(1.0-EPS)*RHOL + JL=(1.0-XFL)*MFLOW/RHOL + JG=XFL*MFLOW/RHOG + IF(EPS.NE.0) THEN + SLIP=JG*(1.0-EPS)/(JL*EPS) + ENDIF + ELSE +* superheated steam + RHO=RHOG + REL=MFLOW*HD/ZMUG + PRL=ZMUG*CPG/ZKG + ENDIF +*---- +* THERMAL EXCHANGE BETWEEN CLAD AND FLUID USING THE DITTUS AND BOELTER +* CORRELATION (SINGLE PHASE) OR CHEN CORRELATION (SATURATED BOILING) +*---- + IF(IHCONV.EQ.0) THEN + ITER=0 + KWA=99 + DO + ITER=ITER+1 + IF(ITER.GT.500) THEN + WRITE(HSMG,'(30HTHMH2O: HCONV FAILURE IN SLICE,I5,1H.)') K + CALL XABORT(HSMG) + ENDIF + HA=0.023*(ZKONE/HD)*REL**0.8*PRL**0.4 + F=1.0 + S=1.0 + IF((XFL.EQ.XFL0).OR.(TSCLAD.LE.TSAT).OR.(KWA.EQ.0)) THEN +* Single-phase convection. Use Dittus-Boelter correlation + KWA=0 + HB=0.0 + K0=0 + XFL=XFL0 + EPS=EPS0 + SLIP=SLIP0 + ELSE IF(HMAVG.LT.HGSAT) THEN +* Subcooled convection. Use Dittus-Boelter and Forster-Zuber +* correlations +* XM: Martinelli parameter +* F: Reynolds number factor +* S: nucleate boiling suppression factor +* SIGM: surface tension in N/m +* HA: Dittus-Boelter coefficient +* HB: Forster-Zuber coefficient +* + IF(HMAVG.LT.HLSAT) THEN + KWA=1 + ELSE + KWA=2 + ENDIF + XM=(XFL/(1.0-XFL))**0.9*(RHOL/RHOG)**0.5*(ZMUG/ZMUL)**0.1 + IF(XM.LE.0.100207) THEN + F=1.0 + ELSE + F=2.35*(0.213+XM)**0.736 + ENDIF + RE=REL*F**1.25 + S=1.0/(1.0+2.53E-6*RE**1.17) + PR=PINLET/10**6 + SIGM=-7.2391E-6*PR**3+2.8345E-4*PR**2-5.1566E-3*PR+4.2324E-2 + HA=0.023*(ZKL/HD)*REL**0.8*PRL**0.4 + DTSAT=TSCLAD-TSAT + IF(IFLUID.EQ.0) THEN + CALL THMSAP(PW, TSCLAD) + ELSE + CALL THMHSP(PW, TSCLAD) + ENDIF + DP=PW-PINLET +* Forster-Zuber equation + HLV=HGSAT-HLSAT + HB=0.00122*ZKL**0.79*CPL**0.45*RHOL**0.49/(ZMUL**0.29* + > SIGM**0.5*HLV**0.24*RHOG**0.24)*DTSAT**0.24*DP**0.75 + ELSE +* Superheated steam. Use Mokry correlation + KWA=3 + IF(IFLUID.EQ.0) THEN + CALL THMPT(PINLET,TSCLAD,RHOW,R2,R3,R4,R5) + ELSE IF(IFLUID.EQ.1) THEN + CALL THMHPT(PINLET,TSCLAD,RHOW,R2,R3,R4,R5) + ENDIF + HA=0.0061*(ZKG/HD)*REL**0.904*PRL**0.684*(RHOW/RHO)**0.564 + HB=0.0 + ENDIF +* Chen correlation + HCONV=F*HA+S*HB + IF(HCONV.LE.0.0) THEN + WRITE(HSMG,'(34HTHMH2O: DRY OUT REACHED IN CHANNEL,3I5)') + > I,J,K + CALL XABORT(HSMG) + ENDIF + IF(ITIME.EQ.0) THEN + TWAL=(PHI+S*HB*TSAT+F*HA*TCALO)/(S*HB+F*HA) + ELSE + ZNUM=ZF(1)+RADCL*S*HB*TSAT+RADCL*F*HA*TCALO + ZDEN=ZF(2)+RADCL*S*HB+RADCL*F*HA + TWAL=MAX(273.15,ZNUM/ZDEN) + PHI=MAX(0.0,(ZF(1)-TWAL*ZF(2))/RADCL) + ENDIF + IF(ABS(TSCLAD-TWAL).GT.1.0E-5*TSCLAD) THEN + TSCLAD=TWAL + ELSE + EXIT + ENDIF + ENDDO + ELSE IF(IHCONV.EQ.1) THEN + IF(ITIME.EQ.0) THEN + TSCLAD=TCALO+PHI/KHCONV + ELSE + RCHC=RADCL*KHCONV + TSCLAD=MAX(273.15,(ZF(1)+RCHC*TCALO)/(ZF(2)+RCHC)) + PHI=(ZF(1)-TSCLAD*ZF(2))/RADCL + ENDIF + ENDIF +*---- +* COMPUTE INITIAL BULK LIQUID ENTHALPY SUBCOOLING DHSUB +*---- + IF((ISUBM.GT.0).AND.(K0.EQ.0).AND.LFIRST) THEN + DTSUB=0.0 + IF(ISUBM.EQ.1) THEN +* Bowring correlation +* Reference: R. W. Bowring, "Physical Model, Based on Bubble +* Detachment, and Calculation of Steam Voidage in the Subcooled +* Region of a Heated Channel," OECD Report HPR-10, 1962. +* Equation3 (3) and (17) + VC=MFLOW/RHOL + ETA=14.0+0.1*PINLET/1.01325E+05 + DTSUB=ETA*PHI/VC*1.0E-6 + ELSE IF(ISUBM.EQ.2) THEN +* Saha-Zuber subcooling model +* PE: Peclet number + PE=MFLOW*CPL*HD/ZKL + IF(PE.LE.70000.0) THEN + DTSUB=PHI*HD/(455.0*ZKL) + ELSE +* reactor conditions + DTSUB=154.0*PHI/(MFLOW*CPL) + ENDIF + ENDIF + IF(TCALO.GE.TSAT-DTSUB) K0=K + DSAT=TSCLAD-TCALO-DTSUB + DHSUB=CPL*DTSUB + LFIRST=.FALSE. + ENDIF + RETURN + END |