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 :== 3;

	(* 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
		A * drho_dt[i] = - (node[i].mdot - node[i-1].mdot)/dz;
		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;
		A * drhou_dt[i] = qdot_t[i] - (node[i].Hdot - node[i-1].Hdot)/dz;
		rho_w*A_w*c_w*dTw_dt[i] = qdot_s - qdot_l[i] - qdot_t[i];
		qdot_l[i] = h_ext*(1{PI}*D_2)*(T_w[i] - T_amb);
		qdot_t[i] = h_int*(1{PI}*D)  *(T_w[i] - node[i].T);
	END FOR;

	t IS_A time;
METHODS
METHOD specify;
	RUN node[1].specify;
	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 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;
	(* FIX node[3].rho; *)

END specify;
METHOD values;
	node[1].T := 400 {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;
END values;
METHOD on_load;
	RUN specify;
	FOR i IN [1..n] DO
		RUN node[i].bound_self;
	END FOR;
	RUN values;
	RUN ode_init;
    EXTERNAL solvernotes(SELF);
	(* RUN ode_init;
	RUN solve; *) (* after fixing the states and freeing the derivatives *)
END on_load;
METHOD fixed_states;
	t := 0 {s};
	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 free_derivs;
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 fix_states;
METHOD ode_init;
	(* get the model into the required state for solving as ODE *)
	t.ode_type := -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;

		(*
		p[i].obs_id := 4*i;
		x[i].obs_id := 4*i+1;
		qdot_t[i].obs_id := 4*i+2;
		T_w[i].obs_id := 4*i+3;
		*)
	END FOR;

	FOR i IN [1,n] DO
		p[i].obs_id := 4*i;
		x[i].obs_id := 4*i+1;
	END FOR;
	FOR i IN [2,n] DO
		(* qdot_t[i].obs_id := 4*i+2; *)
		T_w[i].obs_id := 4*i+3;
	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 :== 3;
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	A * drho_dt[i] = - (node[i].mdot - node[i-1].mdot)/dz;
75	1/Admdot_dt[i] = -(node[i].p-node[i-1].p)/dz - f/D/2node[i].rhonode[i].v^2 (node[i].rhovel[i]^2 - node[i-1].rhovel[i-1]^2)/dz;
76	A * drhou_dt[i] = qdot_t[i] - (node[i].Hdot - node[i-1].Hdot)/dz;
77	rho_wA_wc_w*dTw_dt[i] = qdot_s - qdot_l[i] - qdot_t[i];
78	qdot_l[i] = h_ext(1{PI}D_2)*(T_w[i] - T_amb);
79	qdot_t[i] = h_int(1{PI}D) *(T_w[i] - node[i].T);
80	END FOR;
81
82	t IS_A time;
83	METHODS
84	METHOD specify;
85	RUN node[1].specify;
86	FIX qdot_s;
87	FIX D, D_2, L;
88	FIX h_int, c_w, rho_w, h_ext;
89	FIX f, mu_f;
90	FIX T_amb;
91	(* fix derivatives to zero *)
92	FOR i IN [2..n] DO
93	FIX drho_dt[i]; FREE node[i].rho;
94	FIX dmdot_dt[i]; FREE node[i].mdot;
95	FIX drhou_dt[i]; FREE rhou[i];
96	FIX dTw_dt[i]; FREE T_w[i];
97	END FOR;
98	(* FIX node[3].rho; *)
99
100	END specify;
101	METHOD values;
102	node[1].T := 400 {K};
103	node[1].x := 0.1;
104	qdot_s := 0 {W/m};
105	D := 60 {mm}; D_2 := 70 {mm};
106	L := 100 {m};
107	A_w := 1{PI}*D_2^2;
108	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};
109	f := 0.005; mu_f := 4.5e-5 {Pa*s};
110	T_amb := 300 {K};
111	FOR i IN [2..n] DO
112	drho_dt[i] := 0 {kg/m^3/s};
113	dmdot_dt[i] := 0 {kg/s/s};
114	drhou_dt[i] := 0 {kJ/m^3/s};
115	dTw_dt[i] := 0 {K/s};
116	END FOR;
117	END values;
118	METHOD on_load;
119	RUN specify;
120	FOR i IN [1..n] DO
121	RUN node[i].bound_self;
122	END FOR;
123	RUN values;
124	RUN ode_init;
125	EXTERNAL solvernotes(SELF);
126	(* RUN ode_init;
127	RUN solve; ) ( after fixing the states and freeing the derivatives *)
128	END on_load;
129	METHOD fixed_states;
130	t := 0 {s};
131	qdot_s := 10 {W/m};
132	FOR i IN [2..n] DO
133	FREE drho_dt[i]; FIX node[i].rho;
134	FREE dmdot_dt[i]; FIX node[i].mdot;
135	FREE drhou_dt[i]; FIX rhou[i];
136	FREE dTw_dt[i]; FIX T_w[i];
137	END FOR;
138	END free_derivs;
139	METHOD fixed_derivs;
140	FOR i IN [2..n] DO
141	FIX drho_dt[i]; FREE node[i].rho;
142	FIX dmdot_dt[i]; FREE node[i].mdot;
143	FIX drhou_dt[i]; FREE rhou[i];
144	FIX dTw_dt[i]; FREE T_w[i];
145	END FOR;
146	END fix_states;
147	METHOD ode_init;
148	(* get the model into the required state for solving as ODE *)
149	t.ode_type := -1;
150
151	FOR i IN [2..n] DO
152	drho_dt[i].ode_id := 4i; node[i].rho.ode_id := 4i;
153	drho_dt[i].ode_type := 2; node[i].rho.ode_type := 1;
154
155	dmdot_dt[i].ode_id := 4i+1; node[i].mdot.ode_id := 4i+1;
156	dmdot_dt[i].ode_type := 2; node[i].mdot.ode_type := 1;
157
158	drhou_dt[i].ode_id := 4i+2; rhou[i].ode_id := 4i+2;
159	drhou_dt[i].ode_type := 2; rhou[i].ode_type := 1;
160
161	dTw_dt[i].ode_id := 4i+3; T_w[i].ode_id := 4i+3;
162	dTw_dt[i].ode_type := 2; T_w[i].ode_type := 1;
163
164	(*
165	p[i].obs_id := 4*i;
166	x[i].obs_id := 4*i+1;
167	qdot_t[i].obs_id := 4*i+2;
168	T_w[i].obs_id := 4*i+3;
169	*)
170	END FOR;
171
172	FOR i IN [1,n] DO
173	p[i].obs_id := 4*i;
174	x[i].obs_id := 4*i+1;
175	END FOR;
176	FOR i IN [2,n] DO
177	(* qdot_t[i].obs_id := 4i+2; )
178	T_w[i].obs_id := 4*i+3;
179	END FOR;
180	END ode_init;
181
182	METHOD fix_outlet_quality;
183	FIX x[n];
184	FREE node[1].mdot;
185	END fix_outlet_quality;
186
187	METHOD reinit;
188	RUN on_load;
189	EXTERNAL solve(SELF);
190	RUN fixed_states;
191	END reinit;
192
193	END dsgsat2;
194	ADD NOTES IN dsgsat2;
195	'QRSlv' iterationlimit {50}
196	END NOTES;