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
		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;
	t := 0 {s};
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);
	EXTERNAL solve(SELF);
	RUN fixed_states;
END on_load;
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 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;
	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;

		(*
		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 := 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	johnpye	923	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	johnpye	1046	REQUIRE "steam/satsteamstream.a4c";
16	johnpye	923
17			MODEL dsgsat2;
18			n IS_A integer_constant;
19	johnpye	1055	n :== 5;
20	johnpye	923
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	johnpye	1055	t := 0 {s};
118	johnpye	923	END values;
119			METHOD on_load;
120			RUN specify;
121			FOR i IN [1..n] DO
122			RUN node[i].bound_self;
123			END FOR;
124			RUN values;
125	johnpye	1049	RUN ode_init;
126	johnpye	1046	EXTERNAL solvernotes(SELF);
127	johnpye	1055	EXTERNAL solve(SELF);
128			RUN fixed_states;
129	johnpye	923	END on_load;
130	johnpye	1049	METHOD fixed_states;
131	johnpye	1046	qdot_s := 10 {W/m};
132	johnpye	923	FOR i IN [2..n] DO
133	johnpye	1049	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	johnpye	923	END FOR;
138	johnpye	1049	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	johnpye	1050	t.obs_id := 1;
151	johnpye	923
152			FOR i IN [2..n] DO
153			drho_dt[i].ode_id := 4i; node[i].rho.ode_id := 4i;
154			drho_dt[i].ode_type := 2; node[i].rho.ode_type := 1;
155
156			dmdot_dt[i].ode_id := 4i+1; node[i].mdot.ode_id := 4i+1;
157			dmdot_dt[i].ode_type := 2; node[i].mdot.ode_type := 1;
158
159			drhou_dt[i].ode_id := 4i+2; rhou[i].ode_id := 4i+2;
160			drhou_dt[i].ode_type := 2; rhou[i].ode_type := 1;
161
162			dTw_dt[i].ode_id := 4i+3; T_w[i].ode_id := 4i+3;
163			dTw_dt[i].ode_type := 2; T_w[i].ode_type := 1;
164
165			(*
166			p[i].obs_id := 4*i;
167			x[i].obs_id := 4*i+1;
168			qdot_t[i].obs_id := 4*i+2;
169			T_w[i].obs_id := 4*i+3;
170			*)
171			END FOR;
172
173	johnpye	1055	FOR i IN [1..n] DO
174	johnpye	1050	p[i].obs_id := 1 + 4*i;
175			x[i].obs_id := 2 + 4*i+1;
176	johnpye	923	END FOR;
177	johnpye	1055	FOR i IN [] DO
178	johnpye	1050	(* qdot_t[i].obs_id := 3 + 4i; )
179	johnpye	1055	node[i].mdot.obs_id := 3 + 4*i;
180	johnpye	1050	T_w[i].obs_id := 4 + 4*i;
181	johnpye	923	END FOR;
182			END ode_init;
183
184			METHOD fix_outlet_quality;
185			FIX x[n];
186			FREE node[1].mdot;
187			END fix_outlet_quality;
188
189	johnpye	1049	METHOD reinit;
190			RUN on_load;
191			EXTERNAL solve(SELF);
192			RUN fixed_states;
193			END reinit;
194
195	johnpye	923	END dsgsat2;
196			ADD NOTES IN dsgsat2;
197			'QRSlv' iterationlimit {50}
198			END NOTES;
199
200