*DECK THMDRV SUBROUTINE THMDRV(MPTHM,IMPX,IX,IY,NZ,XBURN,VOLXY,HZ,CFLUX,POROS, > FNFU,NFD,NDTOT,IFLUID,SNAME,SCOMP,IGAP,IFUEL,FNAME,FCOMP,FCOOL, > FFUEL,ACOOL,HD,PCH,RAD, > MAXIT1,MAXITL,ERMAXT,SPEED,TINLET,POUTLET, > FRACPU,ICONDF,NCONDF,KCONDF,UCONDF,ICONDC,NCONDC,KCONDC,UCONDC, > IHGAP,KHGAP,IHCONV,KHCONV,WTEFF,IFRCDI,ISUBM,FRO,POW,IPRES,IDFM, > TCOMB, DCOOL,TCOOL,TSURF,HCOOL,PCOOL) * *----------------------------------------------------------------------- * *Purpose: * Driver of the steady-state thermal-hydraulics calculation. * *Copyright: * Copyright (C) 2012 Ecole Polytechnique de Montreal. * *Author(s): * A. Hebert * C. Garrido * 08/2023: Modifications to include Molten Salt heat transfer in coolant * C. Garrido * 07/2024: Modifications to include Molten Salt heat transfer in static * fuel * C. Huet * 02/2025: Modifications to include pressure drop calculation * R. Guasch & M. Bellier * 08/2025: Modifications to include mass+momentum+energy conservation equation * solution using a Drift-Flux Model. * *Parameters: input * MPTHM directory of the THM object containing steady-state * thermohydraulics data. * IMPX printing index (=0 for no print). * IX position of mesh along X direction. * IY position of mesh along Y direction. * NZ number of meshes along Z direction (channel direction). * XBURN burnup distribution in MWday/tonne. * VOLXY mesh area in the radial plane. * HZ Z-directed mesh widths. * CFLUX critical heat flux in W/m^2. * POROS oxyde porosity. * FNFU number of active fuel rods in the fuel bundle. * NFD number of discretization points in fuel region. * NDTOT number of total discretization points in the the fuel * pellet and the cladding. * IFLUID type of fluid (0=H2O; 1=D2O; 2=SALT). * SNAME Name of the molten salt (e.g. "LiF-BeF2") * SCOMP Composition of the molten salt (e.g. "0.66-0.34") * FCOOL power density fraction in coolant. * FFUEL power density fraction in fuel. * ACOOL coolant cross section area in m^2. * HD hydraulic diameter of one assembly in m. * PCH heating perimeter in m. * RAD fuel and clad radii in m. * MAXIT1 maximum number of conduction iterations. * MAXITL maximum number of center-pellet iterations. * ERMAXT convergence criterion. * SPEED inlet flow velocity in m/s. * TINLET inlet temperature in K. * POUTLET outlet pressure in Pa. * FRACPU plutonium fraction in fuel. * ICONDF fuel conductivity flag (0=Stora-Chenebault or COMETHE/ * 1=user-provided polynomial + inverse term). * NCONDF degree of user-provided fuel conductivity polynomial. * KCONDF polynomial coefficients for fuel conductivity in W/m/K^(k+1) * (except for the two last coefficients which belongs to the * inverse term). * UCONDF required unit of temperature in polynomial for fuel * conductivity (KELVIN or CELSIUS). * ICONDC clad conductivity flag (0=default/1=user-provided * polynomial). * NCONDC degree of user-provided clad conductivity polynomial. * KCONDC polynomial coefficients for clad conductivity in W/m/K^(k+1). * UCONDC required unit of temperature in polynomial for clad * conductivity (KELVIN or CELSIUS). * IHGAP flag indicating HGAP chosen (0=default/1=user-provided). * KHGAP fixed user-provided HGAP value in W/m^2/K. * IHCONV flag indicating HCONV chosen (0=default/1=user-provided). * KHCONV fixed user-provided HCONV value in W/m^2/K. * WTEFF surface temperature's weighting factor in effective fuel * temperature. * IFRCDI flag indicating if average approximation is forced during * fuel conductivity evaluation (0=default/1=average * approximation forced). * ISUBM subcooling model (0: one-phase; 1: Jens-Lottes model; * 2: Saha- Zuber model). * FRO radial power form factors. * POW power distribution in W. * IGAP Flag indicating if the gap is considered (0=gap/1=no gap) * IFUEL type of fuel (0=UO2/MOX; 1=SALT). * FNAME Name of the molten salt (e.g. "LiF-BeF2") * FCOMP Composition of the molten salt (e.g. "0.66-0.34") * IPRES flag indicating if pressure is to be computed (0=nonstant/ * 1=variable). * IDFM flag indicating if the drift flux model is to be used * (0=Without modifications(Chexal correlation for epsilon, no drift flux model in the Navier-Stokes equations) * /1=EPRI/2=MODEBSTION/3=GERAMP/4=HEM1(VGJ=0)) * *Parameters: output * TCOMB averaged fuel temperature distribution in K. * DCOOL coolant density distribution in g/cc. * TCOOL coolant temperature distribution in K. * TSURF surface fuel temperature distribution in K. * HCOOL coolant enthalpty distribution in J/kg. * PCOOL coolant pressure distribution in Pa. * *----------------------------------------------------------------------- * USE GANLIB USE t_saltdata *---- * SUBROUTINE ARGUMENTS *---- TYPE(C_PTR) MPTHM INTEGER IMPX,IX,IY,NZ,NFD,NDTOT,IFLUID,MAXIT1,MAXITL,IHGAP,IGAP, > IFUEL,IPRES, IDFM REAL XBURN(NZ),VOLXY,CFLUX,POROS,FRACPU,ERMAXT, > SPEED,TINLET,POUTLET, > FFUEL(NZ),ACOOL(NZ),RAD(NDTOT-1,NZ),FNFU(NZ),FCOOL(NZ),HZ(NZ), > KCONDF(NCONDF+3),KCONDC(NCONDC+1),KHGAP,KHCONV,WTEFF,FRO(NFD-1), > POW(NZ),TCOMB(NZ),DCOOL(NZ),TCOOL(NZ),TSURF(NZ),HCOOL(NZ), > PCOOL(NZ),MUT(NZ), HD(NZ), PCH(NZ) CHARACTER UCONDF*12,UCONDC*12 *---- * LOCAL VARIABLES *---- TYPE(tpdata) STP,FTP PARAMETER (KMAXO=100,MAXNPO=40) REAL TRE11(MAXNPO),RADD(MAXNPO),ENT(4),MFLOW,TLC(NZ) CHARACTER HSMG*131,SNAME*32,SCOMP*32,FNAME*32,FCOMP*32 REAL XS(4),TC1,PC(NZ),TP(NZ),RHOL,XFL(NZ),EPS(NZ),HINLET, > TCLAD(NZ),ENTH(NZ),SLIP(NZ),AGM(NZ),QFUEL(NZ),QCOOL(NZ),K11, > VLIQ(NZ),VVAP(NZ) INTEGER KWA(NZ) REAL XX2(MAXNPO),XX3(MAXNPO),ZF(2) DATA XS/-0.861136,-0.339981,0.339981,0.861136/ REAL TBUL(NZ),VGJprime(NZ),HLV(NZ),DGCOOL(NZ),DLCOOL(NZ) INTEGER I REAL PINLET, ERRV, ERRP, ERRD, NORMV, NORMP, NORMD *---- * ALLOCATABLE ARRAYS *---- REAL, ALLOCATABLE, DIMENSION(:) :: VCOOL,TCENT REAL, ALLOCATABLE, DIMENSION(:,:) :: TEMPT REAL, ALLOCATABLE, DIMENSION(:) :: PTEMP, VTEMP, DTEMP *---- * SCRATCH STORAGE ALLOCATION *---- ALLOCATE(VCOOL(NZ),TEMPT(NDTOT,NZ),TCENT(NZ)) ALLOCATE(PTEMP(NZ), VTEMP(NZ), DTEMP(NZ)) *---- * COMPUTE THE INLET FLOW ENTHALPY AND VELOCITY * INITIALIZE PINLET TO POUTLET, WILL BE UPDATED IF IPRES=1 * ELSE PINLET = POUTLET *---- PINLET = POUTLET IF(NDTOT.GT.MAXNPO) CALL XABORT('THMDRV: MAXNPO OVERFLOW.') IF(IFLUID.EQ.0) THEN CALL THMSAT(PINLET,TSAT) ELSE IF(IFLUID.EQ.1) THEN CALL THMHST(PINLET,TSAT) *CGT TODO: GET ALSO FREEZING?? ELSE IF(IFLUID.EQ.2) THEN CALL THMSGT(SNAME,SCOMP,STP,IMPX) CALL THMSST(STP,TSAT,IMPX) *CGT ENDIF IF (IFUEL.EQ.1) THEN CALL THMSGT(FNAME,FCOMP,FTP,IMPX) ENDIF IF(TINLET.GT.TSAT) THEN WRITE(HSMG,'(27HTHMDRV: INLET TEMPERATURE (,1P,E12.4, > 40H K) GREATER THAN SATURATION TEMPERATURE.)') TINLET CALL XABORT(HSMG) ENDIF IF(IFLUID.EQ.0) THEN CALL THMPT(PINLET,TINLET,RHOIN,HINLET,R3,R4,R5) ELSE IF(IFLUID.EQ.1) THEN CALL THMHPT(PINLET,TINLET,RHOIN,HINLET,R3,R4,R5) ELSE IF(IFLUID.EQ.2) THEN CALL THMSPT(STP,TINLET,RHOIN,HINLET,R3,R4,R5,IMPX) ENDIF MFLOW=SPEED*RHOIN HMSUP=HINLET *---- * INITIALIZE VALUES OF STEAM QUALITIES, VOID FRACTION AND DENSITY * PRESSURE, VELOCITY AND TEMPERATURE OF THE COOLANT ALONG THE CHANNEL. *--- DO K=1,NZ EPS(K)=0.0 XFL(K)=0.0 SLIP(K)=1.0 KWA(K)=0 MUT(K)=0.0 QFUEL(K)=0.0 VGJprime(K)=0.0 HLV(K)=0.0 PCOOL(K)=PINLET VCOOL(K)=MFLOW/RHOIN DCOOL(K)=RHOIN DLCOOL(K)=RHOIN DGCOOL(K)=0.0 TCOOL(K)=TINLET HCOOL(K)=HINLET *---- * COMPUTE THE SATURATION TEMPERATURE AND THE THERMODYNAMIC PROPERTIES * IF THE PRESSURE DROP IS COMPUTED *--- IF (IPRES.EQ.1) THEN IF(POW(K).EQ.0.0) CYCLE IF(IFLUID.EQ.0) THEN CALL THMSAT(PCOOL(K),TSAT) ELSE IF(IFLUID.EQ.1) THEN CALL THMHST(PCOOL(K),TSAT) ENDIF TB=TSAT-0.1 IF(TCOOL(K).LT.TB) THEN IF(IFLUID.EQ.0) THEN CALL THMPT(PCOOL(K),TCOOL(K),RHOIN,H11,K11,MUT(K),C11) ELSE IF(IFLUID.EQ.1) THEN CALL THMHPT(PCOOL(K),TCOOL(K),RHOIN,H11,K11,MUT(K),C11) ELSE IF(IFLUID.EQ.2) THEN CALL THMSPT(STP,TCOOL(K),R11,H11,K11,MUT(K),C11,IMPX) ENDIF ELSE IF(IFLUID.EQ.0) THEN CALL THMPT(PCOOL(K),TB,R11,H11,K11,MUT(K),C11) ELSE IF(IFLUID.EQ.1) THEN CALL THMHPT(PCOOL(K),TB,R11,H11,K11,MUT(K),C11) ELSE IF(IFLUID.EQ.2) THEN CALL THMSPT(STP,TB,R11,H11,K11,MUT(K),C11,IMPX) ENDIF ENDIF ENDIF ENDDO *---- * MAIN LOOP ALONG THE 1D CHANNEL. *---- ERRV = 1.0 ERRP = 1.0 ERRD = 1.0 NORMP = PINLET NORMV = SPEED NORMD = RHOIN I=0 IF (IPRES .EQ. 0) GOTO 30 10 CONTINUE *---- * UPDATE HINLET FUNCTION OF INLET PRESSURE AND TEMPERATURE *---- HMSUP=HINLET SPEED=MFLOW/DCOOL(1) *---- * WHILE LOOP FOR PRESSURE AND VELOCITY CONVERGENCE * CHECK FOR CONVERGENCE *---- IF (I .GT. 1000) GOTO 20 IF ((ERRP.LT.5E-4).AND.(ERRV.LT.5E-4).AND.(IDFM.EQ.0)) GOTO 20 IF ((IDFM.GT.0).AND.(I.GT.3)) THEN IF ((ERRP.LT.5E-4).AND.(ERRV.LT.5E-4).AND.(ERRD.LT.5E-4)) THEN GOTO 20 ENDIF ENDIF I = I + 1 PTEMP = PCOOL VTEMP = VCOOL DTEMP = DCOOL SPEED = MFLOW/DCOOL(1) CALL THMPV(SPEED, PCOOL(NZ), VCOOL, DCOOL, > PCOOL, TCOOL, MUT, XFL, HD, NZ, > HZ, EPS, DLCOOL,DGCOOL, VGJprime, IDFM, ACOOL) * Extrapolate from first two values of PCOOL to get PINLET at first face. * This ensures that computed HINLET is not HCOOL(1) PINLET = (3.0/2.0)*PCOOL(1) - (1.0/2.0)*PCOOL(2) IF (IFLUID.EQ.0) THEN CALL THMPT(PINLET, TINLET, RHOIN, HINLET, R3, R4, R5) ELSE IF(IFLUID.EQ.1) THEN CALL THMHPT(PINLET,TINLET,RHOIN,HINLET,R3,R4,R5) ELSE IF(IFLUID.EQ.2) THEN CALL THMSPT(STP,TINLET,RHOIN,HINLET,R3,R4,R5,IMPX) ENDIF * Update inlet enthalpy based on computed inlet pressure. HMSUP = HINLET 30 CONTINUE *---- * MAIN LOOP ALONG THE 1D CHANNEL. *---- K0=0 ! onset of nucleate boiling point DO K=1,NZ IF(POW(K).EQ.0.0) CYCLE IF(IFLUID.EQ.0) THEN CALL THMSAT(PCOOL(K),TSAT) ELSE IF(IFLUID.EQ.1) THEN CALL THMHST(PCOOL(K),TSAT) ENDIF TBUL(K)=TSAT *---- * COMPUTE THE LINEAR POWER, THE VOLUMIC POWER AND THE THERMAL EXCHANGE * COEFFICIENT OF THE GAP *---- DV=VOLXY*HZ(K) * linear power in W/m POWLIN=(POW(K)/DV)*VOLXY/FNFU(K) * volumic power in W/m^3 QFUEL(K)=POW(K)*FFUEL(K)/DV QCOOL(K)=POW(K)*FCOOL(K)/DV *---- * INITIALIZATION OF PINCELL TEMPERATURES *---- IF(POW(K).EQ.0.0) CYCLE IF(IMPX.GT.4) WRITE(6,190) K DO L=1,NDTOT TRE11(L)=TCOMB(K) ENDDO DO L=1,NDTOT-1 RADD(L)=RAD(L,K) ENDDO *---- * COMPUTE THE POWER DENSITY AND HEAT FLOW ALONG THE CHANNEL *---- * out-of-clad heat flow in W/m2 IF(IMPX.GT.5) WRITE(6,'(15H THMDRV: QFUEL(,I5,2H)=,1P,E12.4, > 6H W/m2.)') K,QFUEL(K) PHI2=0.5*QFUEL(K)*RAD(NFD,K)**2/RAD(NDTOT-1,K) IF(PHI2.GT.CFLUX) THEN WRITE(HSMG,'(23HTHMDRV: THE HEAT FLUX (,1P,E12.4,5H) IS , > 37HGREATER THAN THE CRITICAL HEAT FLUX (,E12.4,2H).)') > PHI2,CFLUX CALL XABORT(HSMG) ENDIF *---- * COMPUTE FOUR VALUES OF ENTHALPY IN J/KG TO BE USED IN GAUSSIAN * INTEGRATION. DELTH1 IS THE ENTHALPY INCREASE IN EACH AXIAL MESH. *---- IF (IDFM.EQ.0) THEN DELTH1=(PCH(K)/ACOOL(K)*PHI2+QCOOL(K))*HZ(K)/MFLOW ELSE DELTH1= (PCH(K)/ACOOL(K)*PHI2+QCOOL(K))*HZ(K)*ACOOL(K) ENDIF IF ((K.GT.1).AND.(IDFM.GT.0)) THEN DELTH1= (PCH(K)/ACOOL(K)*PHI2+QCOOL(K))*HZ(K)*ACOOL(K) DELTH1 = DELTH1 + ((VCOOL(K-1) + EPS(K-1)*(DLCOOL(K-1)- > DGCOOL(K-1))/DCOOL(K-1)*VGJprime(K-1)) > + (VCOOL(K) + EPS(K)*(DLCOOL(K)-DGCOOL(K))/ > DCOOL(K)*VGJprime(K)))/2*(PCOOL(K-1)*ACOOL(K-1)-PCOOL(K) > *ACOOL(K)) DELTH1 = DELTH1 +(EPS(K-1)*DGCOOL(K-1)*(DLCOOL(K-1)/ > DCOOL(K-1))*HLV(K-1)*VGJprime(K-1)*ACOOL(K-1))-(EPS(K)* > DGCOOL(K)*(DLCOOL(K)/DCOOL(K))*HLV(K)*VGJprime(K)*ACOOL(K)) DELTH1 = DELTH1/MFLOW/ACOOL(K) ENDIF DO I1=1,4 POINT=(1.0+XS(I1))/2.0 ENT(I1)=HMSUP+POINT*DELTH1 ENDDO HMSUP=HMSUP+DELTH1 *---- * COMPUTE THE VALUE OF THE DENSITY AND THE CLAD-COOLANT HEAT TRANSFER * COEFFICIENT *---- IF(K.GT.1) THEN XFL(K)=XFL(K-1) EPS(K)=EPS(K-1) SLIP(K)=SLIP(K-1) ENDIF *CGT IF ((IFLUID.EQ.0).OR.(IFLUID.EQ.1)) THEN CALL THMH2O(0,IX,IY,K,K0,PCOOL(K),MFLOW,HMSUP,ENT,HD(K), > IFLUID,IHCONV,KHCONV,ISUBM,RAD(NDTOT-1,K),ZF,VCOOL(K), > IDFM,PHI2,XFL(K),EPS(K),SLIP(K),ACOOL(K),PCH(K),HZ(K),TCALO, > RHO,RHOL,RHOG,TRE11(NDTOT),KWA(K),VGJprime(K),HLV(K)) ELSEIF (IFLUID.EQ.2) THEN CALL THMSAL(IMPX,0,IX,IY,K,K0,MFLOW,HMSUP,ENT,HD(K),STP, > IHCONV,KHCONV,ISUBM,RAD(NDTOT-1,K),ZF,PHI2,XFL(K), > EPS(K),SLIP(K),HZ(K),TCALO,RHO,RHOL,TRE11(NDTOT), > KWA(K)) ENDIF *CGT *---- * STEADY-STATE SOLUTION OF THE CONDUCTION EQUATIONS IN A FUEL PIN. *---- DTINV=0.0 IF(IGAP.EQ.0) THEN CALL THMROD(IMPX,NFD,NDTOT-1,MAXIT1,MAXITL,ERMAXT,DTINV,RADD, > TRE11,TRE11,QFUEL(K),FRO,TRE11(NDTOT),POWLIN,XBURN(K), > POROS,FRACPU,ICONDF,NCONDF,KCONDF,UCONDF,ICONDC,NCONDC, > KCONDC,UCONDC,IHGAP,KHGAP,IFRCDI,TC1,XX2,XX3,ZF) ELSE CALL THMRNG(IMPX,NFD,NDTOT-1,MAXIT1,MAXITL,ERMAXT,DTINV,RADD, > TRE11,TRE11,QFUEL(K),FRO,TRE11(NDTOT),XBURN(K), > POROS,FRACPU,ICONDF,NCONDF,KCONDF,UCONDF,ICONDC,NCONDC, > KCONDC,UCONDC,IFRCDI,IFUEL,FTP,TC1,XX2,XX3,ZF) ENDIF * DO K1=1,NDTOT-1 TRE11(K1)=XX2(K1)+TRE11(NDTOT)*XX3(K1) ENDDO *---- * RECOVER MESHWISE TEMPERATURES AND FLUID DENSITY. BY DEFAULT, USE THE * ROWLANDS FORMULA TO COMPUTE THE EFFECTIVE FUEL TEMPERATURE, OTHERWISE * USE USER-SPECIFIED WEIGHTING FACTOR. *---- TCOMB(K)=(1.0-WTEFF)*TC1+WTEFF*TRE11(NFD) TCENT(K)=TC1 TSURF(K)=TRE11(NFD) TCLAD(K)=TRE11(NDTOT) TCOOL(K)=TCALO DCOOL(K)=RHO DLCOOL(K)=RHOL HCOOL(K)=HMSUP PC(K)=PINLET TP(K)=TCLAD(K) TLC(K)=TCOOL(K) ENTH(K)=HCOOL(K) AGM(K)=MFLOW ! constant flow rate DO K2=1,NDTOT TEMPT(K2,K)=TRE11(K2) ENDDO IF (IPRES .EQ. 0) THEN PCOOL(K)=PINLET VCOOL(K)=MFLOW/DCOOL(K) ENDIF *---- * COMPUTE THE SATURATION TEMPERATURE AND THE THERMODYNAMIC PROPERTIES * IF THE PRESSURE DROP IS COMPUTED *--- IF (IPRES.EQ.1) THEN IF(POW(K).EQ.0.0) CYCLE IF(IFLUID.EQ.0) THEN CALL THMSAT(PCOOL(K),TSAT) ELSE IF(IFLUID.EQ.1) THEN CALL THMHST(PCOOL(K),TSAT) ENDIF TB=TSAT-0.1 IF(TCOOL(K).LT.TB) THEN IF(IFLUID.EQ.0) THEN CALL THMPT(PCOOL(K),TCOOL(K),RHOIN,H11,K11,MUT(K),C11) ELSE IF(IFLUID.EQ.1) THEN CALL THMHPT(PCOOL(K),TCOOL(K),RHOIN,H11,K11,MUT(K),C11) ELSE IF(IFLUID.EQ.2) THEN CALL THMSPT(STP,TCOOL(K),R11,H11,K11,MUT(K),C11,IMPX) ENDIF ELSE IF(IFLUID.EQ.0) THEN CALL THMPT(PCOOL(K),TB,R11,H11,K11,MUT(K),C11) ELSE IF(IFLUID.EQ.1) THEN CALL THMHPT(PCOOL(K),TB,R11,H11,K11,MUT(K),C11) ELSE IF(IFLUID.EQ.2) THEN CALL THMSPT(STP,TB,R11,H11,K11,MUT(K),C11,IMPX) ENDIF ENDIF ENDIF ENDDO *---- * IF THE PRESSURE DROP IS COMPUTED, COMPUTE THE * THE PRESSURE AND VELOCITY RESIDUALS * IF DFM IS ACTIVATED, COMPUTE DCOOL RESIDUALS *---- IF (IPRES .EQ. 0) GOTO 20 ERRV = 0.0 ERRP = 0.0 ERRD = 0.0 NORMV = 0.0 NORMP = 0.0 NORMD = 0.0 DO K=1,NZ * Under relaxation of coolant pressure and velocity. VCOOL(K) = 0.40*VCOOL(K) + (1.00-0.40)*VTEMP(K) PCOOL(K) = 0.40*PCOOL(K) + (1.00-0.40)*PTEMP(K) ERRV = ERRV + (VCOOL(K)-VTEMP(K))**2 NORMV = NORMV + VCOOL(K)**2 ERRP = ERRP + (PCOOL(K)-PTEMP(K))**2 NORMP = NORMP + PCOOL(K)**2 IF (IDFM.GT.0) THEN * Under relaxation of coolant density. DCOOL(K) = 0.40*DCOOL(K) + (1.00-0.40)*DTEMP(K) ERRD = ERRD + (DCOOL(K) - DTEMP(K))**2 NORMD = NORMD + DCOOL(K)**2 ENDIF ENDDO NORMV = SQRT(NORMV) NORMP = SQRT(NORMP) ERRV = SQRT(ERRV) / NORMV ERRP = SQRT(ERRP) / NORMP IF (IDFM.GT.0) THEN NORMD = SQRT(NORMD) ERRD = SQRT(ERRD) / NORMD ENDIF GO TO 10 20 CONTINUE IF (I.GE.1000) THEN PRINT *, 'ERRV =', ERRV PRINT *, 'ERRP =', ERRP PRINT *, 'ERRD =', ERRD CALL XABORT('THMDRV: MAXIMUM NB OF ITERATIONS REACHED.') ELSE IF(IMPX.GT.0) THEN WRITE(6,'(37H THMDRV: CONVERGENCE REACHED AT ITER=,I5,1H.)') I ENDIF *---- * RECONSTRUCT THE PHASE VELOCITIES FROM VCOOL, EPS and VGJ *---- DO K=1,NZ IF (IDFM .GT. 0) THEN VLIQ(K) = VCOOL(K) - (1.0/(1.0- EPS(K)) - DLCOOL(K)/DCOOL(K)) > * VGJprime(K) VVAP(K) = VCOOL(K) + DLCOOL(K)/DCOOL(K) *VGJprime(K) ELSE VLIQ(K) = VCOOL(K) VVAP(K) = VCOOL(K) ENDIF ENDDO *---- * PRINT THE THERMOHYDRAULICAL PARAMETERS *---- IF(IMPX.GT.4) THEN WRITE(6,250) 'POW', POW(:NZ) WRITE(6,250) 'PCOOL', PCOOL(:NZ) WRITE(6,250) 'VCOOL', VCOOL(:NZ) WRITE(6,250) 'DCOOL', DCOOL(:NZ) WRITE(6,250) 'TCOOL', TCOOL(:NZ) WRITE(6,250) 'EPS', EPS(:NZ) WRITE(6,250) 'XFL', XFL(:NZ) WRITE(6,250) 'TSAT', TBUL(:NZ) WRITE(6,250) 'MUT', MUT(:NZ) ENDIF *---- * PRINT THE OUTLET THERMOHYDRAULICAL PARAMETERS *---- IF(IMPX.GT.3) THEN WRITE(6,'(/16H THMDRV: CHANNEL,2I6/1X,27(1H-))') IX,IY WRITE(6,210) ' ____________________________________________', > '_____________________________________________________', > '_____________________________________________________', > '______________' WRITE(6,210) '| | TCOMB | TSURF | DCOOL ', > ' | TCOOL | PCOOL | HCOOL | ', > 'QFUEL | QCOOL | VOID | QUAL |', > ' SLIP | FLOW |', > '| | K | K | Kg/m3 | ', > ' K | Pa | J/Kg | W/m3 ', > ' | W/m3 | | | ', > ' | REGIME |' WRITE(6,210) '|_____|____________|____________|____________', > '_|_____________|_____________|_____________|_________', > '____|_____________|___________|_____________|________', > '_____|________|' DO L=NZ,1,-1 IF(L.EQ.1) THEN WRITE(6,220) '| BOT |',TCOMB(L),' |',TSURF(L), > ' |',DCOOL(L),' |',TCOOL(L),' |',PCOOL(L), > ' |',HCOOL(L),' |',QFUEL(L),' |',QCOOL(L),' |', > EPS(L),' |',XFL(L),' |',SLIP(L),' |',KWA(L),' |' ELSEIF(L.EQ.NZ) THEN WRITE(6,220) '| TOP |',TCOMB(L),' |',TSURF(L), > ' |',DCOOL(L),' |',TCOOL(L),' |',PCOOL(L), > ' |',HCOOL(L),' |',QFUEL(L),' |',QCOOL(L),' |', > EPS(L),' |',XFL(L),' |',SLIP(L),' |',KWA(L),' |' ELSE WRITE(6,230) '| ',L,' |',TCOMB(L),' |',TSURF(L), > ' |',DCOOL(L),' |',TCOOL(L),' |',PCOOL(L), > ' |',HCOOL(L),' |',QFUEL(L),' |',QCOOL(L),' |', > EPS(L),' |',XFL(L),' |',SLIP(L),' |',KWA(L),' |' ENDIF ENDDO WRITE(6,210) '|_____|____________|____________|____________', > '_|_____________|_____________|_____________|_________', > '____|_____________|___________|_____________|________', > '_____|________|' WRITE(6,240) MFLOW ENDIF *---- * MODIFICATION OF THE VECTORS TO FIT THE GEOMETRY OF THE CHANNELS AND * THE BUNDLES AND WRITE THE DATA IN LCM OBJECT THM *---- CALL LCMPUT(MPTHM,'PRESSURE',NZ,2,PCOOL) CALL LCMPUT(MPTHM,'DENSITY',NZ,2,DCOOL) CALL LCMPUT(MPTHM,'LIQUID-DENS',NZ,2,DLCOOL) CALL LCMPUT(MPTHM,'ENTHALPY',NZ,2,HCOOL) CALL LCMPUT(MPTHM,'VELOCITIES',NZ,2,VCOOL) CALL LCMPUT(MPTHM,'V-LIQ',NZ,2,VLIQ) CALL LCMPUT(MPTHM,'V-VAP',NZ,2,VVAP) CALL LCMPUT(MPTHM,'EPSILON',NZ,2,EPS) CALL LCMPUT(MPTHM,'EPSOUT',1,2,EPS(NZ)) CALL LCMPUT(MPTHM,'XFL',NZ,2,XFL) CALL LCMPUT(MPTHM,'CENTER-TEMPS',NZ,2,TCENT) CALL LCMPUT(MPTHM,'COOLANT-TEMP',NZ,2,TCOOL) CALL LCMPUT(MPTHM,'POWER',NZ,2,POW) CALL LCMPUT(MPTHM,'TEMPERATURES',NDTOT*NZ,2,TEMPT) CALL LCMPUT(MPTHM,'PINLET',1,2,PINLET) CALL LCMPUT(MPTHM,'TINLET',1,2,TINLET) CALL LCMPUT(MPTHM,'VINLET',1,2,SPEED) CALL LCMPUT(MPTHM,'POULET',1,2,POUTLET) CALL LCMPUT(MPTHM,'RADII',(NDTOT-1)*NZ,2,RAD) *---- * SCRATCH STORAGE DEALLOCATION *---- DEALLOCATE(TCENT,TEMPT,VCOOL) DEALLOCATE(PTEMP, VTEMP, DTEMP) RETURN * 190 FORMAT(/21H THMDRV: AXIAL SLICE=,I5) 210 FORMAT(1X,A,A,A,A) 220 FORMAT(1X,A,F11.2,A,F11.2,A,F12.4,A,F12.2,A,3P,E12.4, > A,1P,E12.4,A,1P,E12.4,A,1P,E12.4,A,0P,F10.4,A,E12.4,A, > E12.4,A,I5,2X,A) 230 FORMAT(1X,A,I3,A,F11.2,A,F11.2,A,F12.4,A,F12.2,A,3P,E12.4, > A,1P,E12.4,A,1P,E12.4,A,1P,E12.4,A,0P,F10.4,A,E12.4,A, > E12.4,A,I5,2X,A) 240 FORMAT(7H MFLXT=,1P,E12.4,8H Kg/m2/s) 250 FORMAT(9H THMDRV: ,A6,1H:,1P,11E12.4/(4X,12E12.4)) END