Initial commit from Polytechnique Montreal

1 files changed, 198 insertions, 0 deletions

diff --git a/Donjon/src/THMSAL.f b/Donjon/src/THMSAL.f
new file mode 100644
index 0000000..2121395
--- /dev/null
+++ b/Donjon/src/THMSAL.f
@@ -0,0 +1,198 @@
+*DECK THMSAL
+      SUBROUTINE THMSAL(IMPX,ITIME,I,J,K,K0,MFLOW,HMAVG,ENT,HD,STP,
+     > IHCONV,KHCONV,ISUBM,RADCL,ZF,PHI,XFL,EPS,SLIP,DZ,TCALO,
+     > RHO,RHOLAV,TSCLAD,KWA)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Adaptation of THMH2O.f for convection of Molten Salts using Gnielinski
+* correlation.
+*
+*Copyright:
+* Copyright (C) 2023 Ecole Polytechnique de Montreal.
+*
+*Author(s): 
+* Cristian Garrido Tamm (cristian.garrido@idom.com)
+*
+*Parameters: input
+* IMPX    printing index (=0 for no print).
+* 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
+* 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
+* STP     tpdata object with correlations to obtain properties of molten salt.
+* 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.
+* DZ      axial mesh width in m.
+*
+*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
+* TSCLAD  clad temperature in K
+* KWA     flow regime (=0: single-phase; =1: subcooled; =2: nucleate
+*         boiling; =3 superheated steam)
+*
+*-----------------------------------------------------------------------
+*
+      USE t_saltdata
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(tpdata) STP
+      INTEGER I,J,K,K0,IHCONV,ISUBM,KWA
+      REAL MFLOW,HMAVG,ENT(4),HD,KHCONV,RADCL,ZF(2),PHI,TCALO,RHO,
+     > RHOLAV,TSCLAD,XFL,EPS,SLIP,DZ
+*----
+*  LOCAL VARIABLES
+*----
+      REAL W(4),HL(4)
+      CHARACTER HSMG*131
+*----
+*  SAVE VARIABLES
+*----
+      SAVE W
+      DATA W /0.347855,0.652145,0.652145,0.347855/
+*----
+*  COMPUTE THE DENSITY AND TEMPERATURE OF THE LIQUID
+*----
+      IF(HMAVG.LT.0.0) CALL XABORT('THMSAL: NEGATIVE INPUT ENTHALPY.')
+      IF(ISUBM.NE.0) CALL XABORT('THMSAL: NOT A ONE PHASE FLOW.')
+      CALL THMSST(STP,TSAT,IMPX)
+      HL(1)=ENT(1)
+      HL(2)=ENT(2)
+      HL(3)=ENT(3)
+      HL(4)=ENT(4)
+      CALL THMSH(STP,HL(1),R11,TL1,IMPX)
+      CALL THMSH(STP,HL(2),R11,TL2,IMPX)
+      CALL THMSH(STP,HL(3),R11,TL3,IMPX)
+      CALL THMSH(STP,HL(4),R11,TL4,IMPX)
+      CALL THMSPT(STP,TL1,RHO1,R2,R3,R4,CP1,IMPX)
+      CALL THMSPT(STP,TL2,RHO2,R2,R3,R4,CP2,IMPX)
+      CALL THMSPT(STP,TL3,RHO3,R2,R3,R4,CP3,IMPX)
+      CALL THMSPT(STP,TL4,RHO4,R2,R3,R4,CP4,IMPX)
+      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 FLUID PROPERTIES
+*  RHO: fluid density
+*  REL: Reynolds number of liquid phase
+*  PRL: Prandtl number of liquid phase
+*----
+      IF(XFL.NE.0.0) THEN
+        CALL XABORT('THMSAL: INVALID VALUE OF FLOW QUALITY')
+      ENDIF
+*     One phase liquid
+      TB=TSAT
+      IF(TL.LT.TB) THEN
+        TCALO=TL
+      ELSE
+        TCALO=TB
+      ENDIF
+      CALL THMSPT(STP,TCALO,R1,R2,ZKONE,ZMUONE,CPONE,IMPX)
+      RHO=RHOLAV
+      REL=MFLOW*HD/ZMUONE
+      PRL=ZMUONE*CPONE/ZKONE
+      ZKL=ZKONE
+      XFL0=XFL
+      EPS0=EPS
+      SLIP0=SLIP
+*----
+*  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
+*CGT CHECK IF REYNOLDS AND PRANDTL ARE IN RANGE OF VALIDITY OF
+*    GNIELINSKI CORRELATION
+        TSCLAD=TCALO
+        IF((REL.LT.2300).OR.(REL.GT.1E6)) THEN
+          WRITE(6,*) "  THMSAL: ***WARNING*** REYNOLDS OUT RANGE."
+        ENDIF
+        IF((PRL.LT.0.6).OR.(PRL.GT.1E5)) THEN
+          WRITE(6,*) "  THMSAL: ***WARNING*** PRANDTL OUT RANGE."
+        ENDIF
+        DO
+          ITER=ITER+1
+          IF(ITER.GT.50) THEN
+            WRITE(HSMG,'(30HTHMSAL: HCONV FAILURE IN SLICE,I5,1H.)') K
+            CALL XABORT(HSMG)
+          ENDIF
+*CGT Changed Dittus-Boelter by Gnielinski correlation
+*CGT PRW: Prandtl number of liquid at wall temperature
+          CALL THMSPT(STP,TSCLAD,R1,R2,ZKONE,ZMUONE,CPONE,IMPX)
+          PRW=ZMUONE*CPONE/ZKONE
+          HA=(ZKL/HD)*0.012*(REL**0.87-280)*PRL**0.8*(1+(HD/DZ)
+     >    **(2.0/3.0))*(PRL/PRW)**0.11
+          IF(IMPX.GT.4) THEN
+            WRITE(6,*) 'THMSAL: REL,PRL,PRW,HA=',REL,PRL,PRW,HA
+          ENDIF
+          F=1.0
+          S=1.0
+          IF((XFL.EQ.XFL0).OR.(TSCLAD.LE.TSAT).OR.(KWA.EQ.0)) THEN
+*           Single-phase convection. Use Gnielinski correlation
+            KWA=0
+            HB=0.0
+            K0=0
+            XFL=XFL0
+            EPS=EPS0
+            SLIP=SLIP0
+          ELSE 
+            CALL XABORT('THMSAL: INVALID HEAT TRANSFER REGIME')
+          ENDIF
+*         Chen correlation
+          HCONV=F*HA+S*HB
+          IF(HCONV.LE.0.0) THEN
+            WRITE(HSMG,'(34HTHMSAL: 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
+      RETURN
+      END