models/steam/dsgsat3.a4c

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

(*
	An attempt to model direct steam generation in pipe flow, limited to the
	saturated regime, and with constant-valued friction factor. External heat
	loss is also simplified.
*)
REQUIRE "steam/satsteamstream.a4c";

MODEL dsgsat3;
	n IS_A integer_constant;
	n :== 4;

	(* temporal derivatives *)
	drho_dt[2..n] IS_A density_rate;
	dmdot_dt[2..n] IS_A mass_rate_rate;
	du_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 *)
	z_A: A = 1{PI}*D^2/4;
	z_Aw: A_w = 1{PI}*(D_2^2 - D^2)/4;
	dz IS_A distance;
	L IS_A distance;
	z_dz: dz = L / (n - 1);

	(* fluid properties *)
	node[1..n] IS_A satsteamstream;
	
	(* flow properties *)
	vel[1..n] IS_A speed;
	T_w[2..n] IS_A temperature;
	T[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;

	FOR i IN [1..n] CREATE
		z_vel[i]: vel[i] = node[i].v*node[i].mdot/A;
	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;
		T[i], node[i].T ARE_THE_SAME;
	END FOR;

	(* differential equations *)
	FOR i IN [2..n] CREATE
		z_massbal[i]: A * drho_dt[i] * dz = - (node[i].mdot - node[i-1].mdot);

		z_enbal[i]: dz * (qdot_t[i] - node[i].rho * A * du_dt[i]) =
			 + node[i].mdot * (node[i].u - node[i-1].u)
			 + (node[i].p*node[i].v*node[i].mdot - node[i-1].p*node[i-1].v*node[i-1].mdot);

		z_mombal[i]:  - dz/A*dmdot_dt[i] =
						(node[i].p-node[i-1].p) 
						+ dz * f/D/2 * node[i].rho * vel[i]^2
						+ (node[i].rho*vel[i]^2 - node[i-1].rho*vel[i-1]^2);

		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);

(* -- original formulation --
		z_massbal[i]: A * drho_dt[i] * dz = - (node[i].mdot - node[i-1].mdot);
		z_mombal[i]:  dz/A*dmdot_dt[i] = -(node[i].p-node[i-1].p) 
						- f/D/2*node[i].rho*node[i].v^2*(
							node[i].rho*vel[i]^2 - node[i-1].rho*vel[i-1]^2
						);
		z_enbal[i]: dz * (A * drhou_dt[i] - qdot_t[i]) = - (node[i].Hdot - node[i-1].Hdot);
		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 bound_self;
	vel[1..n].upper_bound := 100 {m/s};
	qdot_l[2..n].lower_bound := 0 {W/m};
	FOR i IN [1..n] DO
		RUN node[i].bound_self;
	END FOR;
END bound_self;
METHOD default_self;
	D := 0.06 {m};
	D_2 := 0.07 {m};
	A_w := 0.25{PI}*D_2^2 -0.25{PI}*D^2;
	FOR i IN [1..n] DO
		RUN node[i].default_self;
	END FOR;
END default_self; 
METHOD values;
	L := 1 {m};
	h_int := 100 {W/m^2/K};
	h_ext := 10 {W/m^2/K};
	f := 0.01;
	node[1].mdot := 0.2 {kg/s};
	node[1].p := 7 {bar};
	node[1].x := 0.2;
	qdot_s := 1000 {W/m^2} * D_2 * 10;
	FOR i IN [2..n] DO
		dmdot_dt[i] := 0.0 {kg/s/s};
		du_dt[i] := 0 {W/m^3};
		node[i].v := 0.2 {L/kg};
		node[i].rho := 6 {kg/L};
		node[i].dp_dT := +0.5 {kPa/K};
		node[i].p := 5 {bar};
	END FOR;
END values;
METHOD on_load;
	RUN configure_steady;
	RUN ode_init;
END on_load;
(*---------------- a physically sensible steady-state configuration-----------*)
METHOD configure_steady;
	RUN default_self;
	RUN ClearAll;
	RUN specify_steady;
	RUN bound_steady;
	RUN values;
END configure_steady;
METHOD bound_steady;
	RUN bound_self;
	T_w[2..n].upper_bound := 1000 {K};
END bound_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];
	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;
	(* FIX 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 du_dt[i]; FREE node[i].T;
		FREE mdot[i]; FIX dmdot_dt[i];
	END FOR;
END specify_steady;
(*------------------------- the dynamic problem ------------------------------*)
METHOD configure_dynamic;
	FOR i IN [2..n] DO
		FREE drho_dt[i];  FIX node[i].rho;
		FREE dmdot_dt[i]; FIX node[i].mdot;
		FREE du_dt[i];    FIX node[i].u;
		FREE dTw_dt[i];   FIX T_w[i];
		FREE node[i].x;
		FREE node[i].T;
	END FOR;
	t := 0 {s};
END configure_dynamic;
METHOD free_states;
	FOR i IN [2..n] DO
		FREE node[i].rho;
		FREE node[i].mdot;
		FREE node[i].u;
		FREE T_w[i];
	END FOR;
END free_states;	
METHOD ode_init;
	(* add the necessary meta data to allow solving with the integrator *)
	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;
		
		du_dt[i].ode_id := 4*i+2;     node[i].u.ode_id := 4*i+2;
		du_dt[i].ode_type := 2;       node[i].u.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 + 10*i;
		x[i].obs_id :=         2 + 10*i;
		node[i].mdot.obs_id := 4 + 10*i;
	END FOR;
	FOR i IN [2..n] DO
		qdot_t[i].obs_id :=    3 + 10*i;
		T_w[i].obs_id :=       5 + 10*i;
		T[i].obs_id :=         6 + 10*i;
	END FOR;
END ode_init;
METHOD fix_outlet_quality;
	FIX x[n];
	FREE node[1].mdot;
END fix_outlet_quality;

END dsgsat3;
ADD NOTES IN dsgsat2;
	'QRSlv' iterationlimit {50}
END NOTES;


1	johnpye	1179	REQUIRE "ivpsystem.a4l";
2			REQUIRE "atoms.a4l";
3			REQUIRE "johnpye/thermo_types.a4c";
4
5			(*
6	johnpye	1275	An attempt to model direct steam generation in pipe flow, limited to the
7			saturated regime, and with constant-valued friction factor. External heat
8			loss is also simplified.
9	johnpye	1179	*)
10			REQUIRE "steam/satsteamstream.a4c";
11
12			MODEL dsgsat3;
13			n IS_A integer_constant;
14	johnpye	1275	n :== 4;
15	johnpye	1179
16			(* temporal derivatives *)
17			drho_dt[2..n] IS_A density_rate;
18			dmdot_dt[2..n] IS_A mass_rate_rate;
19	johnpye	1275	du_dt[2..n] IS_A power_per_volume;
20	johnpye	1179	dTw_dt[2..n] IS_A temperature_rate;
21
22			(* wall properties *)
23			rho_w IS_A mass_density;
24			D, D_2 IS_A distance;
25			c_w IS_A specific_heat_capacity;
26			A, A_w IS_A area;
27			h_int IS_A heat_transfer_coefficient; (* internal *)
28			h_ext IS_A heat_transfer_coefficient; (* external *)
29			z_A: A = 1{PI}*D^2/4;
30			z_Aw: A_w = 1{PI}*(D_2^2 - D^2)/4;
31			dz IS_A distance;
32			L IS_A distance;
33			z_dz: dz = L / (n - 1);
34
35			(* fluid properties *)
36			node[1..n] IS_A satsteamstream;
37
38			(* flow properties *)
39			vel[1..n] IS_A speed;
40			T_w[2..n] IS_A temperature;
41			T[1..n] IS_A temperature;
42
43			(* constants, for the moment: *)
44			f IS_A positive_factor;
45			(* mu_f IS_A viscosity; *)
46			T_amb IS_A temperature;
47
48			(* system dynamics *)
49			qdot_t[2..n], qdot_l[2..n] IS_A power_per_length;
50			qdot_s IS_A power_per_length;
51
52			FOR i IN [1..n] CREATE
53			z_vel[i]: vel[i] = node[i].v*node[i].mdot/A;
54			END FOR;
55
56			(* some aliases just for easier review of the state of the model *)
57			x[1..n] IS_A fraction;
58			mdot[1..n] IS_A mass_rate;
59			p[1..n] IS_A pressure;
60			FOR i IN [1..n] CREATE
61			x[i], node[i].x ARE_THE_SAME;
62			mdot[i], node[i].mdot ARE_THE_SAME;
63			p[i], node[i].p ARE_THE_SAME;
64			T[i], node[i].T ARE_THE_SAME;
65			END FOR;
66
67			(* differential equations *)
68			FOR i IN [2..n] CREATE
69			z_massbal[i]: A * drho_dt[i] * dz = - (node[i].mdot - node[i-1].mdot);
70
71	johnpye	1267	z_enbal[i]: dz * (qdot_t[i] - node[i].rho * A * du_dt[i]) =
72			+ node[i].mdot * (node[i].u - node[i-1].u)
73			+ (node[i].pnode[i].vnode[i].mdot - node[i-1].pnode[i-1].vnode[i-1].mdot);
74	johnpye	1179
75	johnpye	1267	z_mombal[i]: - dz/A*dmdot_dt[i] =
76			(node[i].p-node[i-1].p)
77			+ dz * f/D/2 * node[i].rho * vel[i]^2
78			+ (node[i].rhovel[i]^2 - node[i-1].rhovel[i-1]^2);
79	johnpye	1179
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
84	johnpye	1275	(* -- original formulation --
85			z_massbal[i]: A * drho_dt[i] * dz = - (node[i].mdot - node[i-1].mdot);
86			z_mombal[i]: dz/A*dmdot_dt[i] = -(node[i].p-node[i-1].p)
87			- f/D/2node[i].rhonode[i].v^2*(
88			node[i].rhovel[i]^2 - node[i-1].rhovel[i-1]^2
89			);
90			z_enbal[i]: dz * (A * drhou_dt[i] - qdot_t[i]) = - (node[i].Hdot - node[i-1].Hdot);
91			z_wall[i]: rho_wA_wc_w*dTw_dt[i] = qdot_s - qdot_l[i] - qdot_t[i];
92			z_loss[i]: qdot_l[i] = h_ext(1{PI}D_2)*(T_w[i] - T_amb);
93			z_trans[i]: qdot_t[i] = h_int(1{PI}D) *(T_w[i] - node[i].T);
94			*)
95	johnpye	1179	END FOR;
96
97			t IS_A time;
98			METHODS
99			METHOD bound_self;
100			vel[1..n].upper_bound := 100 {m/s};
101			qdot_l[2..n].lower_bound := 0 {W/m};
102			FOR i IN [1..n] DO
103			RUN node[i].bound_self;
104			END FOR;
105			END bound_self;
106			METHOD default_self;
107			D := 0.06 {m};
108			D_2 := 0.07 {m};
109			A_w := 0.25{PI}D_2^2 -0.25{PI}D^2;
110			FOR i IN [1..n] DO
111			RUN node[i].default_self;
112			END FOR;
113			END default_self;
114			METHOD values;
115	johnpye	1275	L := 1 {m};
116	johnpye	1267	h_int := 100 {W/m^2/K};
117	johnpye	1179	h_ext := 10 {W/m^2/K};
118	johnpye	1275	f := 0.01;
119	johnpye	1267	node[1].mdot := 0.2 {kg/s};
120			node[1].p := 7 {bar};
121	johnpye	1179	node[1].x := 0.2;
122			qdot_s := 1000 {W/m^2} * D_2 * 10;
123			FOR i IN [2..n] DO
124			dmdot_dt[i] := 0.0 {kg/s/s};
125	johnpye	1275	du_dt[i] := 0 {W/m^3};
126	johnpye	1179	node[i].v := 0.2 {L/kg};
127			node[i].rho := 6 {kg/L};
128			node[i].dp_dT := +0.5 {kPa/K};
129	johnpye	1267	node[i].p := 5 {bar};
130	johnpye	1179	END FOR;
131			END values;
132			METHOD on_load;
133			RUN configure_steady;
134			RUN ode_init;
135			END on_load;
136	johnpye	1275	(---------------- a physically sensible steady-state configuration-----------)
137	johnpye	1179	METHOD configure_steady;
138	johnpye	1275	RUN default_self;
139			RUN ClearAll;
140			RUN specify_steady;
141	johnpye	1179	RUN bound_steady;
142			RUN values;
143			END configure_steady;
144			METHOD bound_steady;
145	johnpye	1275	RUN bound_self;
146	johnpye	1179	T_w[2..n].upper_bound := 1000 {K};
147			END bound_steady;
148	johnpye	1275	METHOD specify_steady;
149			(* change to a proper steady-state problem, with fluid properties FREEd *)
150	johnpye	1179	FOR i IN [1..n] DO
151			RUN node[i].specify;
152			FIX dTw_dt[i]; FREE T_w[i];
153			END FOR;
154			FIX node[1].p;
155			FREE node[1].T;
156			FIX qdot_s;
157			FIX D, D_2, L;
158			FIX h_int, c_w, rho_w, h_ext;
159			FIX f;
160			(* FIX mu_f; *)
161			FIX T_amb;
162			(* fix derivatives to zero *)
163			FOR i IN [2..n] DO
164			(* FIX dmdot_dt[i]; FREE node[i].mdot; *)
165			FREE node[i].x; FIX node[i].p;
166			FIX drho_dt[i]; FREE node[i].p;
167			FIX du_dt[i]; FREE node[i].T;
168			FREE mdot[i]; FIX dmdot_dt[i];
169			END FOR;
170	johnpye	1275	END specify_steady;
171			(------------------------- the dynamic problem ------------------------------)
172	johnpye	1179	METHOD configure_dynamic;
173			FOR i IN [2..n] DO
174			FREE drho_dt[i]; FIX node[i].rho;
175			FREE dmdot_dt[i]; FIX node[i].mdot;
176			FREE du_dt[i]; FIX node[i].u;
177			FREE dTw_dt[i]; FIX T_w[i];
178			FREE node[i].x;
179			FREE node[i].T;
180			END FOR;
181			t := 0 {s};
182			END configure_dynamic;
183			METHOD free_states;
184			FOR i IN [2..n] DO
185			FREE node[i].rho;
186			FREE node[i].mdot;
187			FREE node[i].u;
188			FREE T_w[i];
189			END FOR;
190			END free_states;
191			METHOD ode_init;
192	johnpye	1275	(* add the necessary meta data to allow solving with the integrator *)
193	johnpye	1179	t.ode_type := -1;
194
195			FOR i IN [2..n] DO
196			drho_dt[i].ode_id := 4i; node[i].rho.ode_id := 4i;
197			drho_dt[i].ode_type := 2; node[i].rho.ode_type := 1;
198
199			dmdot_dt[i].ode_id := 4i+1; node[i].mdot.ode_id := 4i+1;
200			dmdot_dt[i].ode_type := 2; node[i].mdot.ode_type := 1;
201
202			du_dt[i].ode_id := 4i+2; node[i].u.ode_id := 4i+2;
203			du_dt[i].ode_type := 2; node[i].u.ode_type := 1;
204
205			dTw_dt[i].ode_id := 4i+3; T_w[i].ode_id := 4i+3;
206			dTw_dt[i].ode_type := 2; T_w[i].ode_type := 1;
207			END FOR;
208
209			FOR i IN [1..n] DO
210	johnpye	1267	p[i].obs_id := 1 + 10*i;
211			x[i].obs_id := 2 + 10*i;
212			node[i].mdot.obs_id := 4 + 10*i;
213	johnpye	1179	END FOR;
214			FOR i IN [2..n] DO
215	johnpye	1267	qdot_t[i].obs_id := 3 + 10*i;
216			T_w[i].obs_id := 5 + 10*i;
217			T[i].obs_id := 6 + 10*i;
218	johnpye	1179	END FOR;
219			END ode_init;
220			METHOD fix_outlet_quality;
221			FIX x[n];
222			FREE node[1].mdot;
223			END fix_outlet_quality;
224
225			END dsgsat3;
226			ADD NOTES IN dsgsat2;
227			'QRSlv' iterationlimit {50}
228			END NOTES;
229
230