models/johnpye/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 "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 solvernotes;
	RUN solve;
	RUN ode_init;
	qdot_s := 10 {W/m};
	RUN solve; (* after fixing the states and freeing the derivatives *)
	
END on_load;
METHOD ode_init;
	(* get the model into the required state for solving as ODE *)
	t.ode_type := -1;
	t := 0 {s};
	FOR i IN [2..n] DO
		FREE drho_dt[i];  FREE node[i].rho;
		FREE dmdot_dt[i]; FREE node[i].mdot;
		FREE drhou_dt[i]; FREE rhou[i];
		FREE dTw_dt[i];   FREE T_w[i];
	END FOR;

	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 solvernotes;
	EXTERNAL solvernotes(SELF);
END solvernotes;

METHOD solve;
	EXTERNAL solve(SELF);
END solve;

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 "satsteamstream.a4c";
16
17
18	MODEL dsgsat2;
19	n IS_A integer_constant;
20	n :== 3;
21
22	(* temporal derivatives *)
23	drho_dt[2..n] IS_A density_rate;
24	dmdot_dt[2..n] IS_A mass_rate_rate;
25	drhou_dt[2..n] IS_A power_per_volume;
26	dTw_dt[2..n] IS_A temperature_rate;
27
28	(* wall properties *)
29	rho_w IS_A mass_density;
30	D, D_2 IS_A distance;
31	c_w IS_A specific_heat_capacity;
32	A, A_w IS_A area;
33	h_int IS_A heat_transfer_coefficient; (* internal *)
34	h_ext IS_A heat_transfer_coefficient; (* external *)
35	A = 1{PI}*D^2/4;
36	A_w = 1{PI}*(D_2^2 - D^2)/4;
37	dz IS_A distance;
38	L IS_A distance;
39	dz = L / n;
40
41	(* fluid properties *)
42	node[1..n] IS_A satsteamstream;
43
44	(* flow properties *)
45	vel[1..n] IS_A speed;
46	T_w[1..n] IS_A temperature;
47
48	(* constants, for the moment: *)
49	f IS_A positive_factor;
50	mu_f IS_A viscosity;
51	T_amb IS_A temperature;
52
53	(* system dynamics *)
54	qdot_t[2..n], qdot_l[2..n] IS_A power_per_length;
55	qdot_s IS_A power_per_length;
56	rhou[1..n] IS_A energy_per_volume;
57
58	FOR i IN [1..n] CREATE
59	vel[i] = node[i].v*node[i].mdot/A;
60	rhou[i] = node[i].rho * node[i].u;
61	END FOR;
62
63	(* some aliases just for easier review of the state of the model *)
64	x[1..n] IS_A fraction;
65	mdot[1..n] IS_A mass_rate;
66	p[1..n] IS_A pressure;
67	FOR i IN [1..n] CREATE
68	x[i], node[i].x ARE_THE_SAME;
69	mdot[i], node[i].mdot ARE_THE_SAME;
70	p[i], node[i].p ARE_THE_SAME;
71	END FOR;
72
73	(* differential equations *)
74	FOR i IN [2..n] CREATE
75	A * drho_dt[i] = - (node[i].mdot - node[i-1].mdot)/dz;
76	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;
77	A * drhou_dt[i] = qdot_t[i] - (node[i].Hdot - node[i-1].Hdot)/dz;
78	rho_wA_wc_w*dTw_dt[i] = qdot_s - qdot_l[i] - qdot_t[i];
79	qdot_l[i] = h_ext(1{PI}D_2)*(T_w[i] - T_amb);
80	qdot_t[i] = h_int(1{PI}D) *(T_w[i] - node[i].T);
81	END FOR;
82
83	t IS_A time;
84	METHODS
85	METHOD specify;
86	RUN node[1].specify;
87	FIX qdot_s;
88	FIX D, D_2, L;
89	FIX h_int, c_w, rho_w, h_ext;
90	FIX f, mu_f;
91	FIX T_amb;
92	(* fix derivatives to zero *)
93	FOR i IN [2..n] DO
94	FIX drho_dt[i]; FREE node[i].rho;
95	FIX dmdot_dt[i]; FREE node[i].mdot;
96	FIX drhou_dt[i]; FREE rhou[i];
97	FIX dTw_dt[i]; FREE T_w[i];
98	END FOR;
99	(* FIX node[3].rho; *)
100
101	END specify;
102	METHOD values;
103	node[1].T := 400 {K};
104	node[1].x := 0.1;
105	qdot_s := 0 {W/m};
106	D := 60 {mm}; D_2 := 70 {mm};
107	L := 100 {m};
108	A_w := 1{PI}*D_2^2;
109	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};
110	f := 0.005; mu_f := 4.5e-5 {Pa*s};
111	T_amb := 300 {K};
112	FOR i IN [2..n] DO
113	drho_dt[i] := 0 {kg/m^3/s};
114	dmdot_dt[i] := 0 {kg/s/s};
115	drhou_dt[i] := 0 {kJ/m^3/s};
116	dTw_dt[i] := 0 {K/s};
117	END FOR;
118	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	RUN solvernotes;
126	RUN solve;
127	RUN ode_init;
128	qdot_s := 10 {W/m};
129	RUN solve; (* after fixing the states and freeing the derivatives *)
130
131	END on_load;
132	METHOD ode_init;
133	(* get the model into the required state for solving as ODE *)
134	t.ode_type := -1;
135	t := 0 {s};
136	FOR i IN [2..n] DO
137	FREE drho_dt[i]; FREE node[i].rho;
138	FREE dmdot_dt[i]; FREE node[i].mdot;
139	FREE drhou_dt[i]; FREE rhou[i];
140	FREE dTw_dt[i]; FREE T_w[i];
141	END FOR;
142
143	FOR i IN [2..n] DO
144	drho_dt[i].ode_id := 4i; node[i].rho.ode_id := 4i;
145	drho_dt[i].ode_type := 2; node[i].rho.ode_type := 1;
146
147	dmdot_dt[i].ode_id := 4i+1; node[i].mdot.ode_id := 4i+1;
148	dmdot_dt[i].ode_type := 2; node[i].mdot.ode_type := 1;
149
150	drhou_dt[i].ode_id := 4i+2; rhou[i].ode_id := 4i+2;
151	drhou_dt[i].ode_type := 2; rhou[i].ode_type := 1;
152
153	dTw_dt[i].ode_id := 4i+3; T_w[i].ode_id := 4i+3;
154	dTw_dt[i].ode_type := 2; T_w[i].ode_type := 1;
155
156	(*
157	p[i].obs_id := 4*i;
158	x[i].obs_id := 4*i+1;
159	qdot_t[i].obs_id := 4*i+2;
160	T_w[i].obs_id := 4*i+3;
161	*)
162	END FOR;
163
164	FOR i IN [1,n] DO
165	p[i].obs_id := 4*i;
166	x[i].obs_id := 4*i+1;
167	END FOR;
168	FOR i IN [2,n] DO
169	(* qdot_t[i].obs_id := 4i+2; )
170	T_w[i].obs_id := 4*i+3;
171	END FOR;
172	END ode_init;
173
174	METHOD fix_outlet_quality;
175	FIX x[n];
176	FREE node[1].mdot;
177	END fix_outlet_quality;
178
179	METHOD solvernotes;
180	EXTERNAL solvernotes(SELF);
181	END solvernotes;
182
183	METHOD solve;
184	EXTERNAL solve(SELF);
185	END solve;
186
187	END dsgsat2;
188	ADD NOTES IN dsgsat2;
189	'QRSlv' iterationlimit {50}
190	END NOTES;