/[ascend]/branches/fprops-incomp/models/johnpye/fprops/towerarray.py
ViewVC logotype

Contents of /branches/fprops-incomp/models/johnpye/fprops/towerarray.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3470 - (show annotations) (download) (as text)
Thu Feb 13 01:56:33 2020 UTC (7 months, 1 week ago) by jpye
File MIME type: text/x-python
File size: 12418 byte(s)
add reporting of dp

1 # coding=utf-8
2 import matplotlib
3 matplotlib.use('GtkAgg')
4 import matplotlib.pyplot as pl
5 import numpy as np
6 import math
7
8 class Label:
9 def __init__(self,ix,iy,i):
10 self.ix = ix
11 self.iy = iy
12 self.i = i
13 def __repr__(self):
14 return "{%d,%d}"%(self.ix,self,iy)
15
16 class Node:
17 def __init__(self,x,y,lab):
18 self.x = x
19 self.y = y
20 self.lab = lab
21 self.links_to = []
22 self.links_from = []
23 def distance_from(self,cnode):
24 return math.sqrt((cnode.x-self.x)**2 + (cnode.y-self.y)**2)
25
26 def __repr__(self):
27 return "<Node %s @%.0f,%.0f>" % (repr(self.lab),self.x, self.y)
28
29 class Link:
30 def __init__(self,lab1,lab2,id):
31 self.lab1 = lab1
32 self.lab2 = lab2
33 self.id = id
34 def __repr__(self):
35 return "<Link %d: %s--%s>" % (self.id,repr(self.lab1),repr(self.lab2))
36 def plot(self,arr,xsep,ysep):
37 #print [arr[self.lab1].x,arr[self.lab2].x],[arr[self.lab1].y,arr[self.lab1].y]
38 x = arr[self.lab1].x
39 y = arr[self.lab1].y
40 dx = arr[self.lab2].x - arr[self.lab1].x
41 dy = arr[self.lab2].y - arr[self.lab1].y
42 l = math.sqrt(dx**2+dy**2)
43 hw = 0.2 * xsep
44 hl = 1 * hw
45 l1 = l + hl
46 dx *= l/l1
47 dy *= l/l1
48 pl.arrow(x,y,dx,dy,width=xsep*0.03,head_width=hw,head_length=hl,color='#ff0000')
49 pl.annotate("id=%d"%(self.id,),(x+dx/2,y+dy/2))
50 def node_from(self,arr):
51 return arr[self.lab1]
52 def node_to(self,arr):
53 return arr[self.lab2]
54 def reverse(self):
55 l = self.lab1
56 self.lab1 = self.lab2
57 self.lab2 = l
58 def length(self,arr):
59 dx = arr[self.lab2].x - arr[self.lab1].x
60 dy = arr[self.lab2].y - arr[self.lab1].y
61 l = math.sqrt(dx**2+dy**2)
62 return l
63
64 def create_field_model(nx,ny,Qdot_onecoll_MW=10.,mdot_onecoll=None
65 ,ix_pb=None,iy_pb=None,suffix=None,f=None):
66 xsep = 290.
67 ysep = 240.
68
69 if f is None:
70 import sys
71 f = sys.stdout
72 if ix_pb is None:
73 ix_pb = int(nx)/2
74 if iy_pb is None:
75 iy_pb = int(ny)/2
76 if mdot_onecoll is None:
77 mdot_onecoll = 40 * (Qdot_onecoll_MW / 10)
78 if suffix is None:
79 suffix = "_%dx%d_%fMW"%(nx,ny,Qdot_onecoll_MW)
80
81 ix = np.arange(0,nx)
82 x = ix*xsep;
83 iy = np.arange(0,ny)
84 y = iy*ysep;
85 IX,IY = np.meshgrid(ix,iy)
86 X,Y = np.meshgrid(x,y)
87 lab = zip(IX.ravel(),IY.ravel())
88 loc = zip(X.ravel(),Y.ravel())
89
90 arr = {}
91 for i,l in enumerate(lab):
92 print "node label",l
93 arr[l] = Node(x=loc[i][0],y=loc[i][1],lab=l)
94
95 print "arr =",arr
96
97 centre = arr[(ix_pb,iy_pb)]
98
99 # create the link network...
100 links = []
101 id = 0
102 for ix1 in ix:
103 print "col",ix1
104 for iy1 in np.arange(ny-1):
105 l = Link((ix1,iy1),(ix1,iy1+1),id); id += 1
106 print "row",iy1,"to",iy1+1,":",l
107 links.append(l)
108 if ix1 != 0:
109 print ix1
110 links.append(Link((ix1-1,iy_pb),(ix1,iy_pb),id)); id += 1
111
112 print links
113
114 f.write("""
115 MODEL layout{suffix};
116 n_DI IS_A integer_constant;
117 n_DI :== {n};
118 n_PH IS_A integer_constant;
119 n_PH :== {nl};
120 n_PC ALIASES n_PH;
121 DI[0..n_DI-1] IS_A boiler_simple;
122 TH[0..n_DI-1] IS_A throttle;
123 PH[0..n_PH-1] IS_A pipe_heat_loss_insul;
124 PC[0..n_PC-1] IS_A pipe_heat_loss_insul;
125 JO[0..n_DI-1] IS_A merge_generic;
126 TE[0..n_DI-1] IS_A tee_generic;
127
128 PC[0..n_PH-1].config, PH[0..n_PH-1].config ARE_THE_SAME;
129
130 corr_intconv ALIASES PC[0].config.corr_intconv;
131 k_pipe ALIASES PC[0].config.k_pipe;
132 k_insul ALIASES PC[0].config.k_insul;
133 h_ext ALIASES PC[0].config.h_ext;
134 emiss ALIASES PC[0].config.emiss;
135 eps ALIASES PC[0].config.eps;
136
137 corr_intconv :== 'none';\n\n""".format(
138 n=len(arr)
139 ,nl=len(links)
140 ,suffix=suffix
141 ))
142
143 # reverse mappings
144
145 rarr = {}
146 rlinks = {}
147 for i,lab in enumerate(arr.keys()):
148 rarr[lab] = i
149 for i,l in enumerate(links):
150 rlinks[l] = i
151
152 # create the ASCEND code...
153
154 if 0:
155 f.write("\tNOTES\n");
156 # output the positions of each collector (in comments)
157 for i,lab in enumerate(arr.keys()):
158 print i
159 f.write("\t\t'location' DI[{i}].inlet {{{x},{x}}};\n".format(
160 i=i
161 ,x=arr[lab].x
162 ,y=arr[lab].y
163 ))
164 f.write("\tEND NOTES;\n\n");
165
166 for i,l in enumerate(links):
167 A = l.node_from(arr)
168 B = l.node_to(arr)
169 print "link",l,": A_dist = %f, B_dist = %f" % (A.distance_from(centre), B.distance_from(centre))
170 if A.distance_from(centre) > B.distance_from(centre):
171 l.reverse()
172 B.links_from.append(l)
173 A.links_to.append(l)
174 print "link reverse",l
175 else:
176 A.links_from.append(l)
177 B.links_to.append(l)
178
179 for i,n in enumerate(arr):
180 assert len(arr[n].links_to) <= 1
181 # print "NOTE: Too many links_to node %s: %s"%(n,arr[n].links_to)
182
183 f.write("\t(* links at node {i} {lab} *)\n".format(i=i,lab=n))
184
185 if len(arr[n].links_from) == 0:
186 f.write("\tTE[{i}] IS_REFINED_TO tee_trivial;\n\tJO[{i}] IS_REFINED_TO merge_trivial;\n".format(i=i))
187 else:
188 f.write("\tTE[{i}] IS_REFINED_TO tee_n({nout});\n\tJO[{i}] IS_REFINED_TO merge_n({nout});\n".format(
189 i=i
190 ,lab = n
191 ,nin = len(arr[n].links_to)
192 ,nout = len(arr[n].links_from)+1
193 ))
194
195 for j,lt in enumerate(arr[n].links_to):
196 f.write("\tTE[{i}].inlet, PC[{li}].outlet ARE_THE_SAME;\n".format(
197 i=i
198 ,li=rlinks[lt]
199 ))
200
201 f.write("\tTE[{i}].outlet[0], DI[{i}].inlet ARE_THE_SAME;\n".format(i=i))
202 for j,lt in enumerate(arr[n].links_from):
203 f.write("\tTE[{i}].outlet[{j}], PC[{li}].inlet ARE_THE_SAME;\n".format(
204 i=i
205 ,j=j+1
206 ,li=rlinks[lt]
207 ))
208
209 f.write("\tDI[{i}].outlet, TH[{i}].inlet ARE_THE_SAME;\n".format(i=i));
210 f.write("\tTH[{i}].outlet, JO[{i}].inlet[0] ARE_THE_SAME;\n".format(i=i));
211
212 if len(arr[n].links_from) == 0:
213 f.write("\tTH[{i}] IS_REFINED_TO throttle_trivial;\n".format(i=i));
214
215 for j,lt in enumerate(arr[n].links_from):
216 f.write("\tJO[{i}].inlet[{j}], PH[{li}].outlet ARE_THE_SAME;\n".format(
217 i=i
218 ,j=j+1
219 ,li=rlinks[lt]
220 ))
221
222 for j,lt in enumerate(arr[n].links_to):
223 f.write("\tJO[{i}].outlet, PH[{li}].inlet ARE_THE_SAME;\n".format(
224 i=i
225 ,li=rlinks[lt]
226 ))
227
228
229
230 import sys
231 print "centre.lab",centre.lab
232 print "rarr[centre.lab]",rarr[centre.lab]
233
234
235 f.write("""
236 inlet, outlet IS_A stream_node;
237 inlet, TE[{i}].inlet ARE_THE_SAME;
238 outlet, JO[{i}].outlet ARE_THE_SAME;
239
240 cd ALIASES inlet.cd;
241 cd.component :== 'sodium';
242 cd.type :== 'incomp';
243
244 METHODS
245 """.format(
246 i=rarr[centre.lab]
247 ))
248
249 f.write("METHOD set_pipe_lengths;\n");
250 for i,l in enumerate(links):
251 f.write("""\tFIX PH[{i}].L := {L} {{m}};
252 FIX PC[{i}].L := {L} {{m}};\n""".format(
253 i=i
254 ,L = l.length(arr)
255 ))
256 f.write("END set_pipe_lengths;\n");
257
258
259 f.write("""
260 METHOD default_self;
261 FOR i IN [0..n_DI-1] DO
262 RUN TE[i].default_self;
263 RUN DI[i].default_self;
264 RUN TH[i].default_self;
265 RUN JO[i].default_self;
266 END FOR;
267 FOR i IN [0..n_PH-1] DO
268 RUN PH[i].default_self;
269 RUN PC[i].default_self;
270 END FOR;
271 END default_self;
272 END layout{suffix};
273 """.format(
274 i=rarr[centre.lab]
275 ,suffix=suffix
276 ))
277
278 # note, we don't use str.format() here because lots of units with curly {}:
279 f.write("""
280 UNITS
281 (* currency conversions as of 10 Feb 2020 *)
282 AUD = {USD/1.49};
283 EUR = {1.10*USD};
284 END UNITS;
285 """);
286
287 f.write("""
288 MODEL towerarray{suffix} REFINES layout{suffix};
289 """.format(suffix=suffix))
290
291 f.write("""
292 PC[0..n_PH-1].Vel_out, PH[0..n_PH-1].Vel_out ARE_THE_SAME;
293 Vel ALIASES PC[0].Vel_out;
294
295 PC[0..n_PH-1].t_insul, PH[0..n_PH-1].t_insul ARE_THE_SAME;
296
297 FOR i IN [0..n_PH-1] CREATE
298 (* curve fit to Sch10S pipe sizes *)
299 PC[i].t_pipe = 1{mm} * (PC[i].D/1{mm})^0.3 / 1.222;
300 PH[i].t_pipe = 1{mm} * (PH[i].D/1{mm})^0.3 / 1.222;
301 END FOR;
302 (*PC[0..n_PH-1].t_pipe, PH[0..n_PH-1].t_pipe ARE_THE_SAME;
303 t_pipe ALIASES PC[0].t_pipe;*)
304 t_insul ALIASES PC[0].t_insul;
305
306
307 (* total materials *)
308 m_pipe, m_insul, m_sodium IS_A mass;
309 V_insul IS_A volume;
310 rho_pipe, rho_insul IS_A mass_density;
311 L_pipe_tot IS_A distance;
312 m_pipe = rho_pipe * SUM[PC[i].solid_pipe.V + PH[i].solid_pipe.V | i IN [0..n_PH-1]];
313 V_insul = SUM[PC[i].solid_insul.V + PH[i].solid_insul.V | i IN [0..n_PH-1]];
314 L_pipe_tot = SUM[PC[i].L + PH[i].L | i IN [0..n_PH-1]];
315 m_sodium = SUM[PC[i].solid.V / PC[i].inlet.v + PH[i].solid.V / PH[i].inlet.v | i IN [0..n_PH-1]];
316 m_insul = rho_insul * V_insul;
317
318 (* costs *)
319 C_pipe, C_insul, C_inst_pipe, C_inst_insul, C_supp, C_sodium, C_tot IS_A monetary_unit;
320 c_pipe IS_A cost_per_mass;
321 c_insul IS_A cost_per_volume;
322 C_pipe = c_pipe * m_pipe;
323 C_insul = c_insul * V_insul;
324 C_sodium = c_sodium * m_sodium;
325
326 c_inst_pipe_0, c_inst_insul_0, c_supp_0 IS_A cost_per_length;
327 c_inst_pipe_1, c_inst_insul_1 IS_A cost_per_area;
328 c_sodium IS_A cost_per_mass;
329
330 C_inst_pipe = SUM[PC[i].L * (c_inst_pipe_1*PC[i].D_o + c_inst_pipe_0)
331 + PH[i].L * (c_inst_pipe_1*PH[i].D_o + c_inst_pipe_0) | i IN [0..n_PH-1]];
332
333 C_inst_insul = SUM[PC[i].L * (c_inst_insul_1*PC[i].D_o + c_inst_insul_0)
334 + PH[i].L * (c_inst_insul_1*PH[i].D_o + c_inst_insul_0) | i IN [0..n_PH-1]];
335
336 C_supp = c_supp_0 * SUM[PC[i].L + PH[i].L | i IN [0..n_PH-1]];
337
338 C_tot = C_pipe + C_insul + C_inst_pipe + C_inst_insul + C_supp;
339
340 T_amb IS_A temperature;
341 T_amb, PC[0..n_PH-1].T_amb, PH[0..n_PH-1].T_amb ARE_THE_SAME;
342
343 Qdot_onecoll IS_A energy_rate;
344 Qdot_onecoll, DI[0..n_DI-1].Qdot ARE_THE_SAME;
345 DI[0..n_DI-1].eta ARE_THE_SAME; eta_DI ALIASES DI[1].eta;
346
347 Qdot_coll_tot IS_A energy_rate;
348 Qdot_coll_tot = SUM[DI[i].Qdot | i IN [0..n_DI-1]];
349
350 Qdot_loss_tot IS_A energy_rate;
351 Qdot_loss_tot = SUM[PC[i].Q + PH[i].Q|i IN [0..n_PH-1]];
352
353 Qdot_net IS_A energy_rate;
354 Qdot_net = outlet.mdot*outlet.h - inlet.mdot*inlet.h;
355
356 eta_th_array IS_A fraction;
357 eta_th_array = Qdot_net / Qdot_coll_tot;
358
359 (* try: same mass flow for each dish *)
360 DI[0..n_DI-1].mdot ARE_THE_SAME;
361 mdot_onecoll IS_A mass_rate;
362 mdot_onecoll, DI[0].mdot ARE_THE_SAME;
363
364 T_in ALIASES inlet.T;
365 T_out ALIASES outlet.T;
366
367 p_in ALIASES inlet.p;
368 p_out ALIASES outlet.p;
369 dp IS_A delta_pressure;
370 dp = p_in - p_out;
371
372 METHODS
373 METHOD on_load;
374 RUN default_self;
375
376 RUN set_pipe_lengths;
377
378 FIX Qdot_onecoll := %f {MW}; (* value substituted in script! *)
379 FIX eta_DI := 1;
380
381 FIX Vel := 15 {ft/s}; (* NAK Hbk v3 p6: 4-16" size, 15 *)
382
383 (*FIX t_pipe := 4 {mm}; (* FIXME implement equation here *)*)
384 FIX t_insul := 100 {mm};
385
386 FIX eps := 0.09 {mm};
387 FIX h_ext := 10 {W/m^2/K};
388 FIX emiss := 0.8;
389 FIX k_insul := 0.05 {W/m/K}; (* 'microporous insulation board' @ 700��C, https://is.gd/5j9Gkw *)
390 FIX k_pipe := 22.4 {W/m/K}; (* Haynes 230 @ 700��C: https://is.gd/F0RICh *)
391 FIX rho_pipe := 9 {g/cm^3}; (* Haynes 230 @ 700��C: https://is.gd/F0RICh *)
392 FIX rho_insul := 380 {kg/m^3}; (* 'microporous insulation board', approx value from Felix *)
393
394 FIX c_pipe := 8 {USD/kg} (* 316 stainless steel assumed *);
395
396 FIX c_insul := 10000 {USD/m^3}; (* old quote from brandname product...or maybe 250 USD/kg if alibaba.com *)
397 FIX c_inst_pipe_0 := 22 {AUD/m};
398 FIX c_inst_pipe_1 := 0.3734 {AUD/m/mm};
399 FIX c_inst_insul_0 := 29 {AUD/m};
400 FIX c_inst_insul_1 := 0.6257 {AUD/m/mm};
401 FIX c_supp_0 := 100 {USD/m}; (* this one's in yankee dollars *)
402 FIX c_sodium := 3000 {USD/t}; (* based on a quick look at alibaba https://is.gd/pPKCYt *)
403
404 (*
405 installation cost for pipe, per length installed
406 0.3734 AUD/m/mm * d_o_pipe + 22 AUD/m
407
408 installation cost for insulation, per length installed
409 0.6257 AUD/m/mm * d_o_pipe + 29 AUD/m
410
411 anchors and supports
412 - supports with 2.5 m spacing, 175 AUD/support, plus
413 - 400 AUD/anchor, one per every second support.
414 totals to 100 USD/m.
415
416 insulation cost
417 10000 USD/m^3 (or maybe 250....?)
418 *)
419
420 FIX inlet.p := 50 {bar};
421 FIX inlet.T := 520 {K} + 273.15 {K};
422 FIX mdot_onecoll := %f {kg/s}; (* same for all dishes *)
423 FIX T_amb := 300 {K};
424
425 (* initial guesses to help the solver *)
426 PC[0..n_PH-1].inlet.h := 400 {kJ/kg};
427
428 SOLVER QRSlv;
429 OPTION convopt 'RELNOM_SCALE';
430 END on_load;
431 METHOD correct_dp_and_Tout;
432 FREE Vel;
433 FIX p_in := 5 {bar};
434 FIX p_out := 3 {bar};
435 FREE mdot_onecoll;
436 FIX T_out := 740 {K} + 273.15{K};
437 END correct_dp_and_Tout;
438 METHOD show_temperatures;
439 EXTERNAL disharray_temperature_list(SELF);
440 END show_temperatures;
441 """ % (Qdot_onecoll_MW,mdot_onecoll));
442
443 f.write("END towerarray{suffix};\n".format(suffix=suffix));
444
445 pl.figure()
446 pl.axes().set_aspect('equal', 'datalim')
447 pl.plot(X,Y,'bo')
448 pl.axis([x.min() - xsep/2, x.max() + xsep/2, y.min() - ysep/2, y.max() + ysep/2])
449 for lab in arr:
450 pl.annotate("%s\n%d"%(lab,rarr[lab]),(arr[lab].x,arr[lab].y),ha='right')
451
452 for l in links:
453 print "link",l
454 l.plot(arr,xsep,ysep)
455
456 if __name__=='__main__':
457
458 f = open("layout.a4c","w")
459 f.write("""
460 REQUIRE "johnpye/fprops/pipe.a4c";
461 """);
462
463 create_field_model(1,4,iy_pb=0,suffix="_chain1",f=f)
464
465 create_field_model(5,1,ix_pb=0,Qdot_onecoll_MW=70,suffix="_chain2",f=f)
466
467 f.close()
468
469 pl.show()
470
471

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22