-
Notifications
You must be signed in to change notification settings - Fork 21
Expand file tree
/
Copy pathboiling-2012.fee
More file actions
163 lines (133 loc) · 5.81 KB
/
Copy pathboiling-2012.fee
File metadata and controls
163 lines (133 loc) · 5.81 KB
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
##############################
# vertical boiling channel with arbitrary power distribution
# extended clausse & lahey nodalization
# as presented at ENIEF 2012
# Theler G., Clausse A., Bonetto F.,
# A Moving Boiling-Boundary Model of an Arbitrary-Powered Two-Phase Flow Loop
# Mecanica Computacional Volume XXXI, Number 5, Multiphase Flows, pages 695--720, 2012
# https://cimec.org.ar/ojs/index.php/mc/article/view/4091/4017
# updated to work with FeenoX
# jeremy@seamplex.com
##############################
DEFAULT_ARGUMENT_VALUE 1 uniform
DEFAULT_ARGUMENT_VALUE 2 6.5
DEFAULT_ARGUMENT_VALUE 3 11
##############################
# non-dimensional parameters
##############################
Nsub = $2 # subcooling number
Eu = $3 # euler number
Fr = 1 # froude number
Lambda = 3 # distributed friction number
ki = 6 # inlet head loss coefficient
ke = 2 # outlet head loss coefficient
##############################
# phase-space definition
##############################
N1 = 6 # nodes in the one-phase region
VECTOR l SIZE N1
PHASE_SPACE l ui ue m rhoe phi eta hi
# the boiling frontier is equal to the last one-phase node position
# and we refer to it as lambda throughout the file
ALIAS l[N1] AS lambda
##############################
# DAE solver settings
##############################
end_time = 100 # final integration time
# compute the initial derivatives of the differential objects
# from the variables of both differential and algebraic objects
INITIAL_CONDITIONS_MODE FROM_VARIABLES
# forbid implicit declaration of variables from now on to
# detect typos at parse time
IMPLICIT NONE
##############################
# steady state values
##############################
VAR z'
# include the steady-state power profile from a file:
INCLUDE boiling-2012-$(1).fee
# the transient space-dependant power profile in this case
# is constant and equal to the steady-state profile
q(z,t) = qstar(z)
# functions needed for the steady-state computation
lambdastar(Npch) = root(integral(qstar(z'), z', 0, z) - Nsub/Npch, z, 0, 1)
q2phistar(z,Npch) = integral(qstar(z'), z', lambdastar(Npch), z)
F(Npch) = {
(Nsub/Npch + Nsub*q2phistar(1,Npch))^2/(1 + Npch * q2phistar(1,Npch))
- (Nsub/Npch)^2
+ Lambda*(Nsub/Npch)^2*lambdastar(Npch)
+ Lambda*integral((Nsub/Npch + Nsub*q2phistar(z,Npch))^2/(1 + Npch*q2phistar(z,Npch)), z, lambdastar(Npch), 1)
+ ki*(Nsub/Npch)^2
+ ke*(Nsub/Npch + Nsub*q2phistar(1,Npch))^2 / (1 + Npch*q2phistar(1,Npch))
+ 1/Fr * (lambdastar(Npch) + integral(1/(1 + Npch*q2phistar(z,Npch)), z, lambdastar(Npch), 1))
- Eu }
Npch = root(F(Npch), Npch, Nsub+1e-3, 50)
IF Npch<(Nsub+1e-2)
PRINT TEXT "Npch =" Npch TEXT "should be larger than Nsub =" Nsub SEP " "
ABORT
ENDIF
##############################
# initial conditions
##############################
ui_0 = 0.9*Nsub/Npch # disturbance
hi_0 = -Nsub/Npch
ue_0 = Nsub/Npch + Nsub * integral(qstar(z'), z', lambdastar(Npch), 1)
rhoe_0 = 1/(1 + Npch * integral(qstar(z'), z', lambdastar(Npch), 1))
eta_0 = 1
m_0 = lambdastar(Npch) + integral(1/(1 + Npch * eta_0 * integral(qstar(z'),z',lambdastar(Npch),z)), z, lambdastar(Npch), 1)
l_0[i] = root(hi_0 * i/N1 + integral(qstar(z'), z', 0, z), z, 0, 1)
phi_0 = integral((Nsub/Npch + Nsub*integral(qstar(z'), z', lambdastar(Npch), z))/(1 + Npch*eta*integral(qstar(z'), z', lambdastar(Npch), z)), z, lambdastar(Npch), 1)
# stop the integration if certain variables get out of
# the [0:1] interval -> unstable condition
done = done | m>1 | lambda>1 | ui<0 | ui>1
IF done
PRINT TEXT "\# model is out of bounds" Nsub Npch m lambda ui
ABORT
ENDIF
##############################
# the dynamical system equations
##############################
0 .= -1/hi*hi_dot * (N1 - 1-0.5)*(l[1] - 0) + 0.5*(l_dot[1] + 0) - ui - Nsub/Npch * N1/hi * integral(q(z,t), z, 0, l[1])
# TODO: this used to work in wasora
# 0(i)<2:N1> .= -1/hi*hi_dot * (N1 - i-0.5)*(l(i) - l(i-1)) + 0.5*(l_dot(i) + l_dot(i-1)) - ui - Nsub/Npch * N1/hi * integral(q(z,t), z, l(i-1), l(i))
0 .= -1/hi*hi_dot * (N1 - 2-0.5)*(l[2] - l[2-1]) + 0.5*(l_dot[2] + l_dot[2-1]) - ui - Nsub/Npch * N1/hi * integral(q(z,t), z, l[2-1], l[2])
0 .= -1/hi*hi_dot * (N1 - 3-0.5)*(l[3] - l[3-1]) + 0.5*(l_dot[3] + l_dot[3-1]) - ui - Nsub/Npch * N1/hi * integral(q(z,t), z, l[3-1], l[3])
0 .= -1/hi*hi_dot * (N1 - 4-0.5)*(l[4] - l[4-1]) + 0.5*(l_dot[4] + l_dot[4-1]) - ui - Nsub/Npch * N1/hi * integral(q(z,t), z, l[4-1], l[4])
0 .= -1/hi*hi_dot * (N1 - 5-0.5)*(l[5] - l[5-1]) + 0.5*(l_dot[5] + l_dot[5-1]) - ui - Nsub/Npch * N1/hi * integral(q(z,t), z, l[5-1], l[5])
0 .= -1/hi*hi_dot * (N1 - 6-0.5)*(l[6] - l[6-1]) + 0.5*(l_dot[6] + l_dot[6-1]) - ui - Nsub/Npch * N1/hi * integral(q(z,t), z, l[6-1], l[6])
0 .= ui - ue + Nsub*integral(q(z,t), z, lambda, 1)
0 .= rhoe - 1/(1 + Npch * eta * (1 - lambda))
0 .= lambda - m + integral(1/(1 + Npch*eta*integral(q(z',t),z',lambda,z)), z, lambda, 1)
0 .= m_dot + rhoe*ue - ui
0 .= phi - integral((ui + Nsub*integral(q(z',t), z', lambda, z))/(1 + Npch*eta*integral(q(z',t), z', lambda, z)), z, lambda, 1)
0 .= {
+ ui_dot*lambda
+ ui*l_dot(N1)
+ phi_dot
+ Lambda*(ui^2*lambda +
integral((ui + Nsub * integral(q(z',t), z', lambda, z))^2/
( 1 + Npch*eta*integral(q(z',t), z', lambda, z)),
z, lambda, 1) )
+ rhoe * ue^2
- ui^2
+ ki*ui^2
+ ke*rhoe*ue^2
+ m/Fr
- Eu }
# constant inlet enthalpy and pressure drop
0 .= hi_dot
##############################
# output results
##############################
# write information (commented out) in the output header
IF in_static
PRINT TEXT "\# vertical boiling channel with arbitrary power: $(1) (2012)"
PRINT TEXT "\# Npch = " Npch
PRINT TEXT "\# Nsub = " Nsub
PRINT TEXT "\# Fr = " Fr
PRINT TEXT "\# Lambda = " Lambda
PRINT TEXT "\# ki = " ki
PRINT TEXT "\# ke = " ke
PRINT TEXT "\# Eu = " %.10f Eu
ENDIF
PRINT t lambda ui