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
|
!DECK THMDFM
SUBROUTINE THMDFM(PCOOL,VCOOL,HMAVG,HD,TL,TSAT,IDFM,EPS,XFL,RHO,RHOL,RHOG, VGJ, VGJprime, C0, HLV)
!
!-----------------------------------------------------------------------
!
! Purpose:
! Drift-flux Model for the computation of thermohydraulics parameters in two-phase flow
!
!Copyright:
! Copyright (C) 2025 Ecole Polytechnique de Montreal.
!
!Author(s):
! M. Bellier
!
!Parameters: input
! PCOOL pressure in Pascal
! VCOOL coolant velocity in m/s
! HMAVG averaged enthalpy
! HD hydraulic diameter in m
! TL liquid temperature in K
! TSAT saturation temperature in K
! 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)
! EPS input coolant void fraction
!
!
!Parameters: output
! XFL coolant flow quality
! RHO coolant density in Kg/m^3
! RHOL liquid density in kg/m^3
! RHOG vapour density in kg/m^3
! VGJ drift velocity
! C0 concentration parameter
! VGJprime
! HLV delta between liquid and vapour enthaply
!
!-----------------------------------------------------------------------
!
!----
! SUBROUTINE ARGUMENTS
!----
REAL PCOOL,VCOOL,HMAVG,HD,TL,TSAT,EPS,XFL,RHO,RHOL,RHOG, VGJ, VGJprime, C0, HLV
INTEGER IDFM
!----
! LOCAL VARIABLES
!----
REAL EPSold, ERREPS, VLIQ, VVAP, TCALO, HLSAT, HGSAT, ZMUL, ZMUG, CPL, CPG, ZKL, ZKG, ZMU, REY
INTEGER NITER
!----
! INITIALIZE VARIABLES
!----
VGJ = 0
C0 = 1
VGJprime = 0
!----
! MAIN LOOP
!----
NITER=0
ERREPS=1
10 CONTINUE
!----
! SAVE THE OLD EPSILON VALUE
!----
EPSold = EPS
NITER = NITER+1
!----
! TEST ON ERR EPS
!----
IF (NITER .GT. 150) GOTO 20
IF (ERREPS .LT. 1E-8) GOTO 20
!----
! COMPUTE DENSITIES
!----
TCALO=EPS*TSAT+(1.0-EPS)*TL
CALL THMTX(TCALO,0.0,RHOL,HLSAT,ZKL,ZMUL,CPL)
CALL THMTX(TCALO,1.0,RHOG,HGSAT,ZKG,ZMUG,CPG)
RHO = RHOL*(1 - EPS)+ EPS*RHOG
!----
! COMPUTE PHASES VELOCITIES AND REYNOLDS
!----
VLIQ = VCOOL - (1.0/(1.0- EPS) - RHOL/RHO) *VGJprime
VVAP = VCOOL + RHOL/RHO *VGJprime
ZMU = (ZMUL*ZMUG/ (ZMUL*(1.0-EPS) + ZMUG*EPS))
REY = RHO * ABS(VCOOL) * HD / ZMU
!----
! COMPUTE FLOW QUALITY
!----
IF (HLSAT .GT. HMAVG) THEN
XFL = 0
ELSE IF (HMAVG .GT. HGSAT) THEN
XFL = 1
ELSE
XFL = (HMAVG - HLSAT)/(HGSAT - HLSAT)
ENDIF
!----
! COMPUTE VGJ, VGJprime AND C0 AFTER CHOSEN CORRELATION
!----
CALL THMVGJ(VCOOL, RHO, PCOOL, ZMU, XFL, HD, RHOG, RHOL, EPS, IDFM, VGJ, C0)
VGJprime = VGJ + (C0-1)*VCOOL
!----
! COMPUTE HLV
!----
HLV=HGSAT-HLSAT
!----
! COMPUTE NEW EPS VALUE
!----
IF (XFL.EQ.0) THEN
EPS = 0
ELSE IF (XFL.EQ.1) THEN
EPS = 1
ELSE
EPS = XFL / (C0 * (XFL + (RHOG/RHOL) * (1 - XFL)) + (RHOG * VGJ) / (RHOL * VCOOL))
ENDIF
!----
! COMPUTE DELTA BETWEEN EPSold AND EPS
!----
ERREPS = ABS(EPSold - EPS)
GOTO 10
!----
! EXIT LOOP
!----
20 CONTINUE
IF (NITER.GT.150) THEN
PRINT *, 'THMDFM: Maximum number of iterations reached (150)'
ELSE
PRINT *, 'THMDFM: Convergence reached in I = ', NITER, 'iterations'
ENDIF
END
|
