models/steam/dsgsat2.a4c

REQUIRE "ivpsystem.a4l";
REQUIRE "atoms.a4l";
REQUIRE "johnpye/thermo_types.a4c";

IMPORT "johnpye/extpy/extpy";
IMPORT "johnpye/solve";
IMPORT "johnpye/solvernotes";

(*
	This model uses some ASCEND models from the freesteam library. See
	http://freesteam.sf.net/ for more information. This model doesn't actually
	require compiled binaries of freesteam, so you can just download the .a4c
	files if you wish.
*)
REQUIRE "steam/satsteamstream.a4c";

MODEL dsgsat2;
	n IS_A integer_constant;
	n :== 5;

	(* temporal derivatives *)
	drho_dt[2..n] IS_A density_rate;
	dmdot_dt[2..n] IS_A mass_rate_rate;
	drhou_dt[2..n] IS_A power_per_volume;
	dTw_dt[2..n] IS_A temperature_rate;

	(* wall properties *)
	rho_w IS_A mass_density;
	D, D_2 IS_A distance;
	c_w IS_A specific_heat_capacity;
	A, A_w IS_A area;
	h_int IS_A heat_transfer_coefficient; (* internal *)
	h_ext IS_A heat_transfer_coefficient; (* external *)
	A = 1{PI}*D^2/4;
	A_w = 1{PI}*(D_2^2 - D^2)/4;
	dz IS_A distance;
	L IS_A distance;
	dz = L / n;

	(* fluid properties *)
	node[1..n] IS_A satsteamstream;
	
	(* flow properties *)
	vel[1..n] IS_A speed;
	T_w[1..n] IS_A temperature;

	(* constants, for the moment: *)
	f IS_A positive_factor;
	mu_f IS_A viscosity;
	T_amb IS_A temperature;

	(* system dynamics *)	
	qdot_t[2..n], qdot_l[2..n] IS_A power_per_length;
	qdot_s IS_A power_per_length;
	rhou[1..n] IS_A	energy_per_volume;

	FOR i IN [1..n] CREATE
		vel[i] = node[i].v*node[i].mdot/A;
		rhou[i] = node[i].rho * node[i].u;
	END FOR;

	(* some aliases just for easier review of the state of the model *)
	x[1..n] IS_A fraction;
	mdot[1..n] IS_A mass_rate;
	p[1..n] IS_A pressure;
	FOR i IN [1..n] CREATE
		x[i], node[i].x ARE_THE_SAME;
		mdot[i], node[i].mdot ARE_THE_SAME;
		p[i], node[i].p ARE_THE_SAME;
	END FOR;

	(* differential equations *)
	FOR i IN [2..n] CREATE
		z_massbal[i]: A * drho_dt[i] = - (node[i].mdot - node[i-1].mdot)/dz;
		z_mombal[i]:  1/A*dmdot_dt[i] = -(node[i].p-node[i-1].p)/dz 
						- f/D/2*node[i].rho*node[i].v^2*(
							node[i].rho*vel[i]^2 - node[i-1].rho*vel[i-1]^2
						)/dz;
		z_enbal[i]: A * drhou_dt[i] = qdot_t[i] - (node[i].Hdot - node[i-1].Hdot)/dz;
		z_wall[i]:  rho_w*A_w*c_w*dTw_dt[i] = qdot_s - qdot_l[i] - qdot_t[i];
		z_loss[i]:  qdot_l[i] = h_ext*(1{PI}*D_2)*(T_w[i] - T_amb);
		z_trans[i]: qdot_t[i] = h_int*(1{PI}*D)  *(T_w[i] - node[i].T);
	END FOR;

	t IS_A time;
METHODS
METHOD on_load;
	RUN configure_onload;
	EXTERNAL solve(SELF);
	RUN configure_steady;
END on_load;
(*----------------------- a reliable on-load configuration -------------------*)
METHOD configure_onload;
	RUN default_self;
	RUN reset;
	RUN bound_self;
	RUN values;
	RUN ode_init;
	EXTERNAL solve(SELF);
END configure_onload;	
METHOD default_self;
	A := 5 {cm^2};
	FOR i IN [1..n] DO
		RUN node[i].default_self;
		mdot[i] := 0.001 {kg/s};
		vel[i] := 0.01 {m/s};
		x[i] := 0.1;
	END FOR;
END default_self;
METHOD specify;
	(* put the model in a state that we are sure it will solve: no big blocks, etc *)
	RUN ClearAll;
	FOR i IN [1..n] DO
		RUN node[i].specify;
		FIX dTw_dt[i]; FREE T_w[i];
		FREE rhou[i];
	END FOR;
	FIX qdot_s;
	FIX D, D_2, L;
	FIX h_int, c_w, rho_w, h_ext;
	FIX f, mu_f;
	FIX T_amb;
END specify;
METHOD values;
	node[1].T := 424.98 {K};
	node[1].x := 0.1;
	qdot_s := 0 {W/m};
	D := 60 {mm}; D_2 := 70 {mm};
	L := 100 {m};
	A_w := 1{PI}*D_2^2;
	h_int := 10 {W/m^2/K}; c_w := 0.47 {J/g/K}; rho_w := 7.8 {g/cm^3}; h_ext := 10 {W/m^2/K};
	f := 0.005; mu_f := 4.5e-5 {Pa*s};
	T_amb := 300 {K};
	FOR i IN [2..n] DO
		drho_dt[i] := 0 {kg/m^3/s};
		dmdot_dt[i] := 0 {kg/s/s};
		drhou_dt[i] := 0 {kJ/m^3/s};
		dTw_dt[i] := 0 {K/s};
	END FOR;
	t := 0 {s};
END values;	
METHOD bound_self;
	vel[1..n].upper_bound := 100 {m/s};
	FOR i IN [1..n] DO
		RUN node[i].bound_self;
	END FOR;
END bound_self;
(*---------------- a physically sensible steady-state configuration-----------*)
METHOD configure_steady;
	RUN ClearAll;
	RUN specify_steady;
	RUN values_steady;
END configure_steady;
METHOD specify_steady;
	(* change to a proper steady-state problem, with fluid properties FREEd *)
	FOR i IN [1..n] DO
		RUN node[i].specify;
		FIX dTw_dt[i];   FREE T_w[i];
		FREE rhou[i];
	END FOR;
	FIX node[1].p;
	FREE node[1].T;
	FIX qdot_s;
	FIX D, D_2, L;
	FIX h_int, c_w, rho_w, h_ext;
	FIX f, mu_f;
	FIX T_amb;
	(* fix derivatives to zero *)
	FOR i IN [2..n] DO
		(* FIX dmdot_dt[i]; FREE node[i].mdot; *)
		FREE node[i].x; FIX node[i].p;
		FIX drho_dt[i];  FREE node[i].p;
		FIX drhou_dt[i]; FREE node[i].T;
		FREE mdot[i]; FIX dmdot_dt[i];
	END FOR;
END specify_steady;
METHOD values_steady;
	node[1].p := 5 {bar};
	FOR i IN [2..n] DO
		dmdot_dt[i] := 0.0 {kg/s/s};
		node[i].v := 0.2 {L/kg};
		node[i].rho := 6 {kg/L};
	END FOR;
	(* nothing atm *)
END values_steady;
(*------------------------- the dynamic problem ------------------------------*)
METHOD fixed_pressures;
	RUN ClearAll;
	RUN specify;
	FOR i IN [1..n] DO
		FREE node[i].T;
		FIX node[i].p;
		node[i].p := 5 {bar};
		FREE node[i].x;
		FIX dmdot_dt[i];
(*
		FIX dmdot_dt[i];
		FREE node[i].x;
		dmdot_dt[i] := 0 {kg/s/s};*)
	END FOR;
END fixed_pressures;
METHOD fixed_states;
	qdot_s := 10 {W/m};
	FOR i IN [2..n] DO
		FREE drho_dt[i];  FIX node[i].rho;
		FREE dmdot_dt[i]; FIX node[i].mdot;
		FREE drhou_dt[i]; FIX rhou[i];
		FREE dTw_dt[i];   FIX T_w[i];
	END FOR;
END fixed_states;
METHOD fixed_derivs;
	FOR i IN [2..n] DO
		FIX drho_dt[i];  FREE node[i].rho;
		FIX dmdot_dt[i]; FREE node[i].mdot;
		FIX drhou_dt[i]; FREE rhou[i];
		FIX dTw_dt[i];   FREE T_w[i];
	END FOR;
END fixed_derivs;
METHOD ode_init;
	(* add the necessary meta data to allow solving with the integrator *)
	t.ode_type := -1;
	t.obs_id := 1;

	FOR i IN [2..n] DO
		drho_dt[i].ode_id := 4*i;     node[i].rho.ode_id := 4*i;
		drho_dt[i].ode_type := 2;     node[i].rho.ode_type := 1;

		dmdot_dt[i].ode_id := 4*i+1;  node[i].mdot.ode_id := 4*i+1;
		dmdot_dt[i].ode_type := 2;    node[i].mdot.ode_type := 1;
		
		drhou_dt[i].ode_id := 4*i+2;  rhou[i].ode_id := 4*i+2;
		drhou_dt[i].ode_type := 2;    rhou[i].ode_type := 1;

		dTw_dt[i].ode_id := 4*i+3;    T_w[i].ode_id := 4*i+3;
		dTw_dt[i].ode_type := 2;      T_w[i].ode_type := 1;
	END FOR;

	FOR i IN [1..n] DO
		p[i].obs_id := 1 + 4*i;
		x[i].obs_id := 2 + 4*i+1;
	END FOR;
	FOR i IN [] DO
		(* qdot_t[i].obs_id := 3 + 4*i; *)
		node[i].mdot.obs_id := 3 + 4*i;
		T_w[i].obs_id := 4 + 4*i;
	END FOR;
END ode_init;

METHOD fix_outlet_quality;
	FIX x[n];
	FREE node[1].mdot;
END fix_outlet_quality;

METHOD reinit;
	RUN on_load;
	EXTERNAL solve(SELF);
	RUN fixed_states;
END reinit;

END dsgsat2;
ADD NOTES IN dsgsat2;
	'QRSlv' iterationlimit {50}
END NOTES;
1	REQUIRE "ivpsystem.a4l";
2	REQUIRE "atoms.a4l";
3	REQUIRE "johnpye/thermo_types.a4c";
4
5	IMPORT "johnpye/extpy/extpy";
6	IMPORT "johnpye/solve";
7	IMPORT "johnpye/solvernotes";
8
9	(*
10	This model uses some ASCEND models from the freesteam library. See
11	http://freesteam.sf.net/ for more information. This model doesn't actually
12	require compiled binaries of freesteam, so you can just download the .a4c
13	files if you wish.
14	*)
15	REQUIRE "steam/satsteamstream.a4c";
16
17	MODEL dsgsat2;
18	n IS_A integer_constant;
19	n :== 5;
20
21	(* temporal derivatives *)
22	drho_dt[2..n] IS_A density_rate;
23	dmdot_dt[2..n] IS_A mass_rate_rate;
24	drhou_dt[2..n] IS_A power_per_volume;
25	dTw_dt[2..n] IS_A temperature_rate;
26
27	(* wall properties *)
28	rho_w IS_A mass_density;
29	D, D_2 IS_A distance;
30	c_w IS_A specific_heat_capacity;
31	A, A_w IS_A area;
32	h_int IS_A heat_transfer_coefficient; (* internal *)
33	h_ext IS_A heat_transfer_coefficient; (* external *)
34	A = 1{PI}*D^2/4;
35	A_w = 1{PI}*(D_2^2 - D^2)/4;
36	dz IS_A distance;
37	L IS_A distance;
38	dz = L / n;
39
40	(* fluid properties *)
41	node[1..n] IS_A satsteamstream;
42
43	(* flow properties *)
44	vel[1..n] IS_A speed;
45	T_w[1..n] IS_A temperature;
46
47	(* constants, for the moment: *)
48	f IS_A positive_factor;
49	mu_f IS_A viscosity;
50	T_amb IS_A temperature;
51
52	(* system dynamics *)
53	qdot_t[2..n], qdot_l[2..n] IS_A power_per_length;
54	qdot_s IS_A power_per_length;
55	rhou[1..n] IS_A energy_per_volume;
56
57	FOR i IN [1..n] CREATE
58	vel[i] = node[i].v*node[i].mdot/A;
59	rhou[i] = node[i].rho * node[i].u;
60	END FOR;
61
62	(* some aliases just for easier review of the state of the model *)
63	x[1..n] IS_A fraction;
64	mdot[1..n] IS_A mass_rate;
65	p[1..n] IS_A pressure;
66	FOR i IN [1..n] CREATE
67	x[i], node[i].x ARE_THE_SAME;
68	mdot[i], node[i].mdot ARE_THE_SAME;
69	p[i], node[i].p ARE_THE_SAME;
70	END FOR;
71
72	(* differential equations *)
73	FOR i IN [2..n] CREATE
74	z_massbal[i]: A * drho_dt[i] = - (node[i].mdot - node[i-1].mdot)/dz;
75	z_mombal[i]: 1/A*dmdot_dt[i] = -(node[i].p-node[i-1].p)/dz
76	- f/D/2node[i].rhonode[i].v^2*(
77	node[i].rhovel[i]^2 - node[i-1].rhovel[i-1]^2
78	)/dz;
79	z_enbal[i]: A * drhou_dt[i] = qdot_t[i] - (node[i].Hdot - node[i-1].Hdot)/dz;
80	z_wall[i]: rho_wA_wc_w*dTw_dt[i] = qdot_s - qdot_l[i] - qdot_t[i];
81	z_loss[i]: qdot_l[i] = h_ext(1{PI}D_2)*(T_w[i] - T_amb);
82	z_trans[i]: qdot_t[i] = h_int(1{PI}D) *(T_w[i] - node[i].T);
83	END FOR;
84
85	t IS_A time;
86	METHODS
87	METHOD on_load;
88	RUN configure_onload;
89	EXTERNAL solve(SELF);
90	RUN configure_steady;
91	END on_load;
92	(----------------------- a reliable on-load configuration -------------------)
93	METHOD configure_onload;
94	RUN default_self;
95	RUN reset;
96	RUN bound_self;
97	RUN values;
98	RUN ode_init;
99	EXTERNAL solve(SELF);
100	END configure_onload;
101	METHOD default_self;
102	A := 5 {cm^2};
103	FOR i IN [1..n] DO
104	RUN node[i].default_self;
105	mdot[i] := 0.001 {kg/s};
106	vel[i] := 0.01 {m/s};
107	x[i] := 0.1;
108	END FOR;
109	END default_self;
110	METHOD specify;
111	(* put the model in a state that we are sure it will solve: no big blocks, etc *)
112	RUN ClearAll;
113	FOR i IN [1..n] DO
114	RUN node[i].specify;
115	FIX dTw_dt[i]; FREE T_w[i];
116	FREE rhou[i];
117	END FOR;
118	FIX qdot_s;
119	FIX D, D_2, L;
120	FIX h_int, c_w, rho_w, h_ext;
121	FIX f, mu_f;
122	FIX T_amb;
123	END specify;
124	METHOD values;
125	node[1].T := 424.98 {K};
126	node[1].x := 0.1;
127	qdot_s := 0 {W/m};
128	D := 60 {mm}; D_2 := 70 {mm};
129	L := 100 {m};
130	A_w := 1{PI}*D_2^2;
131	h_int := 10 {W/m^2/K}; c_w := 0.47 {J/g/K}; rho_w := 7.8 {g/cm^3}; h_ext := 10 {W/m^2/K};
132	f := 0.005; mu_f := 4.5e-5 {Pa*s};
133	T_amb := 300 {K};
134	FOR i IN [2..n] DO
135	drho_dt[i] := 0 {kg/m^3/s};
136	dmdot_dt[i] := 0 {kg/s/s};
137	drhou_dt[i] := 0 {kJ/m^3/s};
138	dTw_dt[i] := 0 {K/s};
139	END FOR;
140	t := 0 {s};
141	END values;
142	METHOD bound_self;
143	vel[1..n].upper_bound := 100 {m/s};
144	FOR i IN [1..n] DO
145	RUN node[i].bound_self;
146	END FOR;
147	END bound_self;
148	(---------------- a physically sensible steady-state configuration-----------)
149	METHOD configure_steady;
150	RUN ClearAll;
151	RUN specify_steady;
152	RUN values_steady;
153	END configure_steady;
154	METHOD specify_steady;
155	(* change to a proper steady-state problem, with fluid properties FREEd *)
156	FOR i IN [1..n] DO
157	RUN node[i].specify;
158	FIX dTw_dt[i]; FREE T_w[i];
159	FREE rhou[i];
160	END FOR;
161	FIX node[1].p;
162	FREE node[1].T;
163	FIX qdot_s;
164	FIX D, D_2, L;
165	FIX h_int, c_w, rho_w, h_ext;
166	FIX f, mu_f;
167	FIX T_amb;
168	(* fix derivatives to zero *)
169	FOR i IN [2..n] DO
170	(* FIX dmdot_dt[i]; FREE node[i].mdot; *)
171	FREE node[i].x; FIX node[i].p;
172	FIX drho_dt[i]; FREE node[i].p;
173	FIX drhou_dt[i]; FREE node[i].T;
174	FREE mdot[i]; FIX dmdot_dt[i];
175	END FOR;
176	END specify_steady;
177	METHOD values_steady;
178	node[1].p := 5 {bar};
179	FOR i IN [2..n] DO
180	dmdot_dt[i] := 0.0 {kg/s/s};
181	node[i].v := 0.2 {L/kg};
182	node[i].rho := 6 {kg/L};
183	END FOR;
184	(* nothing atm *)
185	END values_steady;
186	(------------------------- the dynamic problem ------------------------------)
187	METHOD fixed_pressures;
188	RUN ClearAll;
189	RUN specify;
190	FOR i IN [1..n] DO
191	FREE node[i].T;
192	FIX node[i].p;
193	node[i].p := 5 {bar};
194	FREE node[i].x;
195	FIX dmdot_dt[i];
196	(*
197	FIX dmdot_dt[i];
198	FREE node[i].x;
199	dmdot_dt[i] := 0 {kg/s/s};*)
200	END FOR;
201	END fixed_pressures;
202	METHOD fixed_states;
203	qdot_s := 10 {W/m};
204	FOR i IN [2..n] DO
205	FREE drho_dt[i]; FIX node[i].rho;
206	FREE dmdot_dt[i]; FIX node[i].mdot;
207	FREE drhou_dt[i]; FIX rhou[i];
208	FREE dTw_dt[i]; FIX T_w[i];
209	END FOR;
210	END fixed_states;
211	METHOD fixed_derivs;
212	FOR i IN [2..n] DO
213	FIX drho_dt[i]; FREE node[i].rho;
214	FIX dmdot_dt[i]; FREE node[i].mdot;
215	FIX drhou_dt[i]; FREE rhou[i];
216	FIX dTw_dt[i]; FREE T_w[i];
217	END FOR;
218	END fixed_derivs;
219	METHOD ode_init;
220	(* add the necessary meta data to allow solving with the integrator *)
221	t.ode_type := -1;
222	t.obs_id := 1;
223
224	FOR i IN [2..n] DO
225	drho_dt[i].ode_id := 4i; node[i].rho.ode_id := 4i;
226	drho_dt[i].ode_type := 2; node[i].rho.ode_type := 1;
227
228	dmdot_dt[i].ode_id := 4i+1; node[i].mdot.ode_id := 4i+1;
229	dmdot_dt[i].ode_type := 2; node[i].mdot.ode_type := 1;
230
231	drhou_dt[i].ode_id := 4i+2; rhou[i].ode_id := 4i+2;
232	drhou_dt[i].ode_type := 2; rhou[i].ode_type := 1;
233
234	dTw_dt[i].ode_id := 4i+3; T_w[i].ode_id := 4i+3;
235	dTw_dt[i].ode_type := 2; T_w[i].ode_type := 1;
236	END FOR;
237
238	FOR i IN [1..n] DO
239	p[i].obs_id := 1 + 4*i;
240	x[i].obs_id := 2 + 4*i+1;
241	END FOR;
242	FOR i IN [] DO
243	(* qdot_t[i].obs_id := 3 + 4i; )
244	node[i].mdot.obs_id := 3 + 4*i;
245	T_w[i].obs_id := 4 + 4*i;
246	END FOR;
247	END ode_init;
248
249	METHOD fix_outlet_quality;
250	FIX x[n];
251	FREE node[1].mdot;
252	END fix_outlet_quality;
253
254	METHOD reinit;
255	RUN on_load;
256	EXTERNAL solve(SELF);
257	RUN fixed_states;
258	END reinit;
259
260	END dsgsat2;
261	ADD NOTES IN dsgsat2;
262	'QRSlv' iterationlimit {50}
263	END NOTES;