1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
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
|
