-
Notifications
You must be signed in to change notification settings - Fork 0
/
tanks.py
257 lines (178 loc) · 10.5 KB
/
tanks.py
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
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
from tools import *
import numpy as np
from scipy.optimize import *
from octopus.main import Fluid
class Propellant(object):
"""Object for handling a propellant, initialised without any idea of the properties
Arguments:
temperature
volume
total mass
Attributes:
pressure
dryness fraction
liquid density
liquid enthalpy
vapour density
vapour enthalpy
Functions:
find_instrinsic_properties
changing one of the arguments post-hoc, calculate all the others and then all attributes
"""
def __init__(self, propellant, temperature, volume, prop_mass, pressure = None):
self.temperature = temperature
self.volume = volume
self.mass = prop_mass
self.propellant_lookup = propellants[propellant]
self.propellant = self.propellant_lookup["prop_name"]
self.phase = self.propellant_lookup["phase"]
self.pressurant_name = self.propellant_lookup["pressurant"]
if self.phase == Propellant_Phase.SELF_PRESSURISING and self.pressurant_name == None:
self.pressurant = None
self.target_pressure = None
self.mass_p = 0
if self.propellant == Propellant_Name.NITROUS:
self.thermophys = nitrous_thermophys
self.den_l, self.den_v, self.h_l, self.h_v, self.latent_heat_vapourisation, self.pressure, ldynvis = self.thermophys(self.temperature)
self.fluid = Fluid('nitrous oxide', method='helmholz', T = self.temperature, P = self.pressure)
self.dryness = ((self.volume / self.mass) - (1/self.den_l)) / ((1/self.den_v) - (1/self.den_l))
self.mass_l = (1 - self.dryness) * self.mass
self.mass_v = self.dryness * self.mass
self.enthalpy = self.mass_l * self.h_l + self.mass_v * self.h_v
if self.mass_l/self.den_l > self.volume:
print("Warning: too much mass of "+self.propellant_lookup["propep_name"])
elif self.phase == Propellant_Phase.LIQUID:
if self.pressurant_name == None:
raise RuntimeError
else:
self.pressurant = Pressurant(self.pressurant_name)
if pressure != None:
self.pressure = pressure
self.target_pressure = pressure
else:
self.pressure = 10 * 1e5
if self.propellant == Propellant_Name.IPA:
self.thermophys = ipa_thermophys
self.den_l, den_v, self.h_l, h_v, latent_heat_vapourisation, pressure, ldynvis = self.thermophys(self.temperature)
self.cv_l = 2680
self.mass_l = prop_mass
self.mass_v = 0
if self.mass_l/self.den_l > self.volume:
print("Warning: too much mass of "+self.propellant_lookup["propep_name"])
if self.pressurant.pressurant_type == Pressurant_Name.HELIUM:
self.pressurant_thermophys = helium_thermophys
self.den_p, self.h_p = self.pressurant_thermophys(self.temperature, self.pressure)
self.mass_p = (self.volume - (self.mass_l/self.den_l))*self.den_p
self.enthalpy = self.mass_l*self.h_l + self.mass_p*self.h_p
elif self.phase == Propellant_Phase.NONE:
self.mass_l = self.mass_v = self.mass_p = 0.0
def self_pressurising_equations_of_state_for_solving(self,arguments,):
"""
A wrapper onto the equations of state from thermophys that allows solving by scipy.optimise.fsolve
arguments [float, float]: [dryness fraction, temperature]
"""
if arguments[1]>309.57 or arguments[1]<183.15 or arguments[0]<0 or arguments[0]>1:
#If either parameter is out of bounds, return a random very large value to encourage the solver to move away
H,V = 1e50,1e50
return(H*V)
else:
dryness, temp = arguments
den_l, den_v, h_l, h_v, c, pres, ldynvis = self.thermophys(temp)
#enthalpy_diff and volume_diff are the difference between predicted and true values, squared to produce a smooth odd function
enthalpy_diff = ((self.mass * dryness * h_v ) + (self.mass * (1 - dryness) * h_l) - self.enthalpy)**2
volume_diff = ((dryness * self.mass / den_v) + ((1-dryness) * self.mass / den_l) - self.volume)**2
return((enthalpy_diff*volume_diff)**2)
def liquid_pressurant_equations_of_state_for_solving(self,old_mass_l):
"""Provides a "solver" for the liquid-pressurant system
Assumes the volume of IPA remains constant (note: this is not entirely valid but is handled in the iterative part of remove_liquid)
Gas behaves adiabatically
"""
gamma = self.pressurant.gamma
old_volume_p = self.volume - (old_mass_l/self.den_l)
new_volume_p = self.volume - (self.mass_l/self.den_l)
new_temp = self.temperature * (old_volume_p/new_volume_p)**(gamma-1)
new_pressure = self.pressure * (old_volume_p/new_volume_p)**(gamma)
return(new_temp, new_pressure)
def remove_liquid(self, mass_out):
if self.phase == Propellant_Phase.SELF_PRESSURISING and self.pressurant == None and self.propellant == Propellant_Name.NITROUS:
self.mass -= mass_out
self.enthalpy -= mass_out * self.h_l
guess = np.array([self.dryness, self.temperature])
bounds = [[0,1],[183.15,309.57]]
new_dryness, new_temp = dual_annealing(self.self_pressurising_equations_of_state_for_solving,bounds,x0=guess, accept=-10.0, maxiter=1000).x
self.dryness = new_dryness
self.temperature = new_temp
self.den_l, self.den_v, self.h_l, self.h_v, self.latent_heat_vapourisation, self.pressure, ldynvis = self.thermophys(self.temperature)
self.mass_l = (1 - self.dryness) * self.mass
self.mass_v = self.dryness * self.mass
self.enthalpy = self.mass_l * self.h_l + self.mass_v * self.h_v
elif self.phase == Propellant_Phase.SELF_PRESSURISING and self.pressurant == None and self.propellant == Propellant_Name.NITROUS and False:
self.mass -= mass_out
self.enthalpy -= mass_out * self.h_l
u = ['p', 'chi']
dryness_estimate = self.dryness#self.mass_v / self.mass
y = [self.pressure, self.dryness]
initial = least_squares(self.fluid.fun_ps, [800, 250], bounds=([0, 0], [np.inf, self.fluid.Tc]), args=[u, y])
den, new_temp = initial.x
properties = self.fluid.get_properties(den, new_temp)
self.dryness = properties["chi"]
self.temperature = new_temp
self.den_l, self.den_v, self.h_l, self.h_v, self.latent_heat_vapourisation, self.pressure, ldynvis = self.thermophys(self.temperature)
self.mass_l = (1 - self.dryness) * self.mass
self.mass_v = self.dryness * self.mass
print(self.enthalpy - (self.mass_l * self.h_l + self.mass_v * self.h_v))
self.enthalpy = self.mass_l * self.h_l + self.mass_v * self.h_v
elif self.phase == Propellant_Phase.LIQUID:
self.mass_l -= mass_out
self.enthalpy -= mass_out * self.h_l
guess = [float(self.temperature)]
bounds = [[180.0, 320.0]]
for i in range(0,10):
new_temp, new_pressure = self.liquid_pressurant_equations_of_state_for_solving(self.mass_l + mass_out)
self.den_l = self.thermophys(new_temp)[0]
self.temp, self.pressure = new_temp, new_pressure
self.den_l, den_v, self.h_l, h_v, latent_heat_vapourisation, pressure, ldynvis = self.thermophys(self.temperature)
self.den_p, self.h_p = self.pressurant_thermophys(self.temperature, self.pressure)
self.enthalpy = self.mass_l*self.h_l + self.mass_p*self.h_p
def pressurant_input_equations_of_state_for_solving(self, temperature_in, pressure_in):
""" Returns the added mass and new temperature when regulated pressurant is added to the system
Source: CUED 1A Thermofluids Example paper 3, Question 6
Bet you never thought you'd see an example paper cited in a real application!
Warning: does not handle a lower pressure in regulator than tank. Currently assumes a check valve, ie will prevent mass flow out
Arguments:
temperature_in float input temperature of pressurant
pressure_in float input pressure and target pressure of pressurant
Returns:
mass_in float mass of pressurant added
final_temp float temperature of pressurant at end
"""
gamma = self.pressurant.gamma
coeffs = [gamma - 1 + (temperature_in/self.temperature),
gamma - (pressure_in/self.pressure) + (temperature_in/self.temperature) + ((self.mass_l * self.cv_l) / (self.mass_p * 8314)),
1 - (pressure_in/self.pressure) + ((self.mass_l * self.cv_l) / (self.mass_p * 8314)) * (1 - (pressure_in/self.pressure))]
roots = np.roots(coeffs)
if roots[0].imag !=0:
roots = [0,0]
f = max(roots)
if f != 0:
mass_in = f * self.mass_p
final_temp = self.temperature * (pressure_in/self.pressure) / (1 + f)
else:
mass_in = 0
final_temp = self.temperature
return(mass_in, final_temp)
def top_up_pressurant(self, input_temperature, target_pressure):
"""Uses an iterative solver
"""
for i in range(0,10):
mass_in, new_temp = self.pressurant_input_equations_of_state_for_solving(input_temperature, target_pressure)
#new_temp, mass_in = dual_annealing(self.pressurant_input_equations_of_state_for_solving, bounds=[[180,320],[0,2]],args=[target_presure ,input_temperature], accept=-7.5, maxiter=2000).x
self.temperature = new_temp
self.mass_p += mass_in
self.pressure = target_pressure
self.den_l, den_v, self.h_l, h_v, latent_heat_vapourisation, pressure, ldynvis = self.thermophys(self.temperature)
self.den_p, self.h_p = self.pressurant_thermophys(self.temperature, self.pressure)
self.enthalpy = self.mass_l*self.h_l + self.mass_p*self.h_p
nitrous = Propellant("nitrous-self-pressurised", 295, 0.15, 110)
#helium = Pressurant(Pressurant_Name.HELIUM)
#ipa = Propellant(Propellant_Phase.LIQUID, Propellant_Name.IPA, 295, 0.1, 50, 10*1e5, helium)