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

IMPORT "johnpye/extpy/extpy";
IMPORT "johnpye/roots";

MODEL dsgsat3;
	n IS_A integer_constant;
	n :== 30;(* with L = 10m: 5,6,7,8,9,10,11 *)
	         (* with L =  5m: 2,3,4,5,7,9,11,12,13,1415,16 *)

	dz IS_A real_constant;
	L IS_A real_constant;
	L :== 11 {m};
	dz :== L / (n-1);

	nodes,butfirst1,upwind4,central IS_A set OF integer_constant;
	nodes :== [1..n];
	butfirst1 :== nodes - [1];
	upwind4 :== nodes - [1,2,n];
	central :== nodes - [1,n];

	(* temporal derivatives *)
	drho_dt[butfirst1] IS_A density_rate;
	(* dmdot_dt[butfirst1] IS_A mass_rate_rate; *)
	du_dt[butfirst1] IS_A specific_energy_rate;
	dTw_dt[butfirst1] 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;

	(* fluid properties *)
	node[nodes] IS_A satsteamstream;
	
	(* flow properties *)
	vel[nodes] IS_A speed;
	T_w[butfirst1] IS_A temperature;
	T[nodes] IS_A temperature;

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

	(* system dynamics *)	
	qdot_t[butfirst1], qdot_l[butfirst1] IS_A power_per_length;
	qdot_s IS_A power_per_length;

	(* some aliases just for easier review of the state of the model *)
	x[nodes] IS_A fraction;
	mdot[nodes] IS_A mass_rate;
	p[nodes] IS_A pressure;
	rho[nodes] IS_A mass_density;
	u[nodes] IS_A specific_energy;
	v[nodes] IS_A specific_volume;

(*	Qdot_t IS_A energy_rate;
	Qdot_t = SUM[qdot_t[i] | i IN butfirst1]; *)

	FOR i IN nodes CREATE
		z_vel[i]: vel[i] = v[i]*mdot[i]/A;
	END FOR;

	FOR i IN nodes 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;
		rho[i],  node[i].rho ARE_THE_SAME;
		u[i],    node[i].u ARE_THE_SAME;
		v[i],    node[i].v ARE_THE_SAME;
	END FOR;

	en_upwind4,en_central,mom_upwind4,mom_central,mass_upwind4,mass_central IS_A set OF integer_constant;

	(* mass conservation *)
	mass_upwind4 :== upwind4;
	mass_central :== central - mass_upwind4;
	FOR i IN mass_upwind4 CREATE (* 4-pt upwind biased *)
		z_massbal1[i]: A * drho_dt[i] * dz =
				- (mdot[i+1] + 6.*mdot[i] - 3.*mdot[i-1] - 2.*mdot[i-2]) / 6.;
	END FOR;
	FOR i IN mass_central CREATE
		z_massbal2[i]: A * drho_dt[i] * dz =
				- (mdot[i+1] - mdot[i-1]) / 2.;
	END FOR;
	FOR i IN butfirst1 - mass_upwind4 - mass_central CREATE
		z_massbal[i]: A * drho_dt[i] * dz = - (mdot[i] - mdot[i-1]);
	END FOR;

	(* energy conservation *)
	en_upwind4 :== [];
	en_central :== central - en_upwind4;
	FOR i IN en_upwind4 CREATE
		z_enbal2[i]: dz * (qdot_t[i] - rho[i] * A * du_dt[i]) =
			 + mdot[i] * (node[i+1].u + 6.*u[i] - 3.*u[i-1] - 2.*u[i-1]) / 6.
			 + (p[i+1]*node[i+1].v*mdot[i+1] - p[i-1]*v[i-1]*mdot[i-1]) / 2.;
	END FOR;
	FOR i IN en_central CREATE
		z_enbal1[i]: dz * (qdot_t[i] - rho[i] * A * du_dt[i]) =
			 + mdot[i] * (u[i] - u[i-1]) (* NOTE: not central *)
			 + (p[i+1]*v[i+1]*mdot[i+1] - p[i-1]*v[i-1]*mdot[i-1]) / 2.;
	END FOR;
	FOR i IN butfirst1 - en_upwind4 - en_central CREATE
		z_enbal[i]: dz * (qdot_t[i] - rho[i] * A * du_dt[i]) =
			 + mdot[i] * (u[i] - u[i-1])
			 + (p[i]*v[i]*mdot[i] - p[i-1]*v[i-1]*mdot[i-1]);
	END FOR;

	(* stationary momentum *)
	FOR i IN butfirst1 CREATE
		z_mombal[i]: p[i] = p[i-1] - dz  * f/D/2 * rho[i] * vel[i]^2;
	END FOR;
	mom_upwind4 :== [];
	mom_central :== [];
(*
	* momentum conservation *
	FOR i IN mom_upwind4 CREATE
		z_mombal2[i]:  - dz/A * dmdot_dt[i]
			 =	(p[i]-p[i-1])  * backdiff for pressure *
				+ dz * f/D/2 * rho[i] * vel[i]^2
				+ (rho[i+1]*vel[i+1]^2 + 6.*rho[i]*vel[i]^2 - 3.*rho[i-1]*vel[i-1]^2 - 2.*rho[i-2]*vel[i-2]^2) / 6.;
	END FOR;
	FOR i IN mom_central CREATE
		z_mombal1[i]:  - dz/A * dmdot_dt[i]
			 =	(p[i+1]-p[i-1]) / 2.
				+ dz * f/D/2 * rho[i] * vel[i]^2
				+ (rho[i+1]*vel[i+1]^2 - rho[i-1]*vel[i-1]^2) / 2.;
	END FOR;
	FOR i IN butfirst1 - mom_upwind4 - mom_central CREATE
		z_mombal[i]:  - dz/A * dmdot_dt[i]
			 =	(p[i]-p[i-1]) 
				+ dz * f/D/2 * rho[i] * vel[i]^2
				+ (rho[i]*vel[i]^2 - rho[i-1]*vel[i-1]^2);
	END FOR;
*)

	(* internal/external convection, and thermal mass of wall -- no spatial derivs here *)
	FOR i IN butfirst1 CREATE
		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] - T[i]);
	END FOR;

	t IS_A time;

METHODS

	METHOD bound_self;
		vel[nodes].upper_bound := 100 {m/s};
		qdot_l[butfirst1].lower_bound := 0 {W/m};
		FOR i IN nodes DO
			RUN node[i].bound_self;
		END FOR;
	END bound_self;
	METHOD default;
		(* these are initial guesses only; fixed parameters are overwritten by 'values' below *)
		t := 0 {s};
		FOR i IN nodes DO
			T[i] := 298 {K};
			vel[i] := 1 {m/s};
			RUN node[i].default_self;
		END FOR;
		FOR i IN butfirst1 DO
			drho_dt[i] := 0 {kg/m^3/s};
(* 			dmdot_dt[i] := 0 {kg/s/s}; *)
			du_dt[i] := 0 {kJ/kg/s};
			dTw_dt[i] := 0 {K/s};
			qdot_t[i] := 0 {W/m};
			qdot_l[i] := 0 {W/m};
			x[i] := x[1];
		END FOR;
	END default; 
	METHOD specify;
		(* change to a proper steady-state problem, with fluid properties FREEd *)
		FOR i IN nodes DO
			RUN node[i].specify;
			FIX dTw_dt[i];   FREE T_w[i];
		END FOR;
		FIX p[1];
		FREE T[1];
		FIX qdot_s;
		FIX D, D_2;
		FIX h_int, c_w, rho_w, h_ext;
		FIX f;
		(* FIX mu_f; *)
		FIX T_amb;
		(* fix derivatives to zero *)
		FOR i IN butfirst1 DO
			FREE x[i]; FIX p[i]; FREE node[i].mdot;
			FIX drho_dt[i];  FREE p[i];
			FIX du_dt[i]; FREE T[i];
(*			FREE mdot[i]; FIX dmdot_dt[i]; *)
		END FOR;
	END specify;
	METHOD values;
		D := 0.06 {m};
		D_2 := 0.07 {m};
		T_amb := 298 {K};
		h_int := 1000 {W/m^2/K};
		h_ext := 5 {W/m^2/K};
		f := 0.03;
		mdot[1] := 0.26 {kg/s};
		p[1]    := 10 {bar};
		x[1]    := 0.23;
		rho_w := 1000 {kg/m^3};
		qdot_s := 100 {W/m^2} * D_2 * 10 * 11{m} / L;
		FOR i IN butfirst1 DO
			T_w[i] := 298 {K};
(*			dmdot_dt[i] := 0.0 {kg/s/s};*)
			du_dt[i] := 0 {kJ/kg/s};
			v[i] := 0.2 {L/kg};
			rho[i] := 6 {kg/L};
			node[i].dp_dT := +0.5 {kPa/K};
			p[i] := 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;
		EXTERNAL defaultself_visit_childatoms(SELF);
		EXTERNAL defaultself_visit_submodels(SELF); 
		RUN ClearAll;
		RUN specify;
		RUN bound_steady;
		RUN values;
	END configure_steady;
	METHOD bound_steady;
		RUN bound_self;
		T_w[butfirst1].upper_bound := 1000 {K};
	END bound_steady;
	(*------------------------- the dynamic problem ------------------------------*)
	METHOD configure_dynamic;
		FOR i IN butfirst1 DO
			FREE drho_dt[i];  FIX rho[i];
(*			FREE dmdot_dt[i]; FIX mdot[i]; *)
			FREE du_dt[i];    FIX u[i];
			FREE dTw_dt[i];   FIX T_w[i];
			FREE x[i];
			FREE T[i];
		END FOR;
		t := 0 {s};
	END configure_dynamic;
	METHOD free_states;
		FOR i IN butfirst1 DO
			FREE rho[i];
(*			FREE mdot[i]; *)
			FREE u[i];
			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 butfirst1 DO
			drho_dt[i].ode_id := 4*i;     rho[i].ode_id := 4*i;
			drho_dt[i].ode_type := 2;     rho[i].ode_type := 1;

(*			dmdot_dt[i].ode_id := 4*i+1;  mdot[i].ode_id := 4*i+1;
			dmdot_dt[i].ode_type := 2;    mdot[i].ode_type := 1;*)
			
			du_dt[i].ode_id := 4*i+2;     u[i].ode_id := 4*i+2;
			du_dt[i].ode_type := 2;       u[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 nodes DO
			(* p[i].obs_id :=         1 + 10*i; *)
			(* x[i].obs_id :=         2 + 10*i; *)
			(* mdot[i].obs_id := 4 + 10*i; *)
			(* node[i].h.obs_id :=    3 + 10*i; *)
		END FOR;
		FOR i IN butfirst1 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;

		mdot[n].obs_id := 1;
		x[n].obs_id := 1;
		p[n].obs_id := 1;
		vel[n].obs_id := 1;
		T_w[n].obs_id := 2;

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

	METHOD roots;
		EXTERNAL roots(SELF);
	END roots;

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	jpye	1331	IMPORT "johnpye/extpy/extpy";
13			IMPORT "johnpye/roots";
14
15	johnpye	1179	MODEL dsgsat3;
16			n IS_A integer_constant;
17	jpye	1383	n :== 30;(* with L = 10m: 5,6,7,8,9,10,11 *)
18	johnpye	1305	(* with L = 5m: 2,3,4,5,7,9,11,12,13,1415,16 *)
19	johnpye	1179
20	johnpye	1305	dz IS_A real_constant;
21			L IS_A real_constant;
22	jpye	1383	L :== 11 {m};
23	johnpye	1305	dz :== L / (n-1);
24
25			nodes,butfirst1,upwind4,central IS_A set OF integer_constant;
26			nodes :== [1..n];
27			butfirst1 :== nodes - [1];
28			upwind4 :== nodes - [1,2,n];
29			central :== nodes - [1,n];
30
31	johnpye	1179	(* temporal derivatives *)
32	johnpye	1305	drho_dt[butfirst1] IS_A density_rate;
33	jpye	1383	(* dmdot_dt[butfirst1] IS_A mass_rate_rate; *)
34	johnpye	1305	du_dt[butfirst1] IS_A specific_energy_rate;
35			dTw_dt[butfirst1] IS_A temperature_rate;
36	johnpye	1179
37			(* wall properties *)
38			rho_w IS_A mass_density;
39			D, D_2 IS_A distance;
40			c_w IS_A specific_heat_capacity;
41			A, A_w IS_A area;
42			h_int IS_A heat_transfer_coefficient; (* internal *)
43			h_ext IS_A heat_transfer_coefficient; (* external *)
44			z_A: A = 1{PI}*D^2/4;
45			z_Aw: A_w = 1{PI}*(D_2^2 - D^2)/4;
46
47			(* fluid properties *)
48	johnpye	1305	node[nodes] IS_A satsteamstream;
49	johnpye	1179
50			(* flow properties *)
51	johnpye	1305	vel[nodes] IS_A speed;
52			T_w[butfirst1] IS_A temperature;
53			T[nodes] IS_A temperature;
54	johnpye	1179
55	johnpye	1305	(* constant, for the moment: *)
56	johnpye	1179	f IS_A positive_factor;
57			(* mu_f IS_A viscosity; *)
58			T_amb IS_A temperature;
59
60			(* system dynamics *)
61	johnpye	1305	qdot_t[butfirst1], qdot_l[butfirst1] IS_A power_per_length;
62	johnpye	1179	qdot_s IS_A power_per_length;
63
64			(* some aliases just for easier review of the state of the model *)
65	johnpye	1305	x[nodes] IS_A fraction;
66			mdot[nodes] IS_A mass_rate;
67			p[nodes] IS_A pressure;
68			rho[nodes] IS_A mass_density;
69			u[nodes] IS_A specific_energy;
70			v[nodes] IS_A specific_volume;
71	jpye	1383
72			(* Qdot_t IS_A energy_rate;
73			Qdot_t = SUM[qdot_t[i] \| i IN butfirst1]; *)
74
75	johnpye	1305	FOR i IN nodes CREATE
76	jpye	1383	z_vel[i]: vel[i] = v[i]*mdot[i]/A;
77			END FOR;
78
79			FOR i IN nodes CREATE
80	johnpye	1305	x[i], node[i].x ARE_THE_SAME;
81	johnpye	1179	mdot[i], node[i].mdot ARE_THE_SAME;
82	johnpye	1305	p[i], node[i].p ARE_THE_SAME;
83			T[i], node[i].T ARE_THE_SAME;
84			rho[i], node[i].rho ARE_THE_SAME;
85			u[i], node[i].u ARE_THE_SAME;
86			v[i], node[i].v ARE_THE_SAME;
87	johnpye	1179	END FOR;
88
89	johnpye	1312	en_upwind4,en_central,mom_upwind4,mom_central,mass_upwind4,mass_central IS_A set OF integer_constant;
90
91	johnpye	1305	(* mass conservation *)
92	johnpye	1312	mass_upwind4 :== upwind4;
93	jpye	1328	mass_central :== central - mass_upwind4;
94	johnpye	1312	FOR i IN mass_upwind4 CREATE (* 4-pt upwind biased *)
95	johnpye	1305	z_massbal1[i]: A * drho_dt[i] * dz =
96			- (mdot[i+1] + 6.mdot[i] - 3.mdot[i-1] - 2.*mdot[i-2]) / 6.;
97			END FOR;
98	johnpye	1312	FOR i IN mass_central CREATE
99	johnpye	1305	z_massbal2[i]: A * drho_dt[i] * dz =
100			- (mdot[i+1] - mdot[i-1]) / 2.;
101			END FOR;
102	johnpye	1312	FOR i IN butfirst1 - mass_upwind4 - mass_central CREATE
103	johnpye	1305	z_massbal[i]: A * drho_dt[i] * dz = - (mdot[i] - mdot[i-1]);
104			END FOR;
105	johnpye	1277
106	johnpye	1305	(* energy conservation *)
107	johnpye	1312	en_upwind4 :== [];
108	jpye	1328	en_central :== central - en_upwind4;
109	johnpye	1312	FOR i IN en_upwind4 CREATE
110	johnpye	1305	z_enbal2[i]: dz * (qdot_t[i] - rho[i] * A * du_dt[i]) =
111			+ mdot[i] * (node[i+1].u + 6.u[i] - 3.u[i-1] - 2.*u[i-1]) / 6.
112			+ (p[i+1]node[i+1].vmdot[i+1] - p[i-1]v[i-1]mdot[i-1]) / 2.;
113	johnpye	1179	END FOR;
114	johnpye	1312	FOR i IN en_central CREATE
115	johnpye	1305	z_enbal1[i]: dz * (qdot_t[i] - rho[i] * A * du_dt[i]) =
116			+ mdot[i] * (u[i] - u[i-1]) (* NOTE: not central *)
117			+ (p[i+1]v[i+1]mdot[i+1] - p[i-1]v[i-1]mdot[i-1]) / 2.;
118			END FOR;
119	johnpye	1312	FOR i IN butfirst1 - en_upwind4 - en_central CREATE
120	johnpye	1305	z_enbal[i]: dz * (qdot_t[i] - rho[i] * A * du_dt[i]) =
121			+ mdot[i] * (u[i] - u[i-1])
122			+ (p[i]v[i]mdot[i] - p[i-1]v[i-1]mdot[i-1]);
123			END FOR;
124	johnpye	1277
125	jpye	1383	(* stationary momentum *)
126			FOR i IN butfirst1 CREATE
127			z_mombal[i]: p[i] = p[i-1] - dz * f/D/2 * rho[i] * vel[i]^2;
128			END FOR;
129	johnpye	1312	mom_upwind4 :== [];
130	jpye	1383	mom_central :== [];
131			(*
132			* momentum conservation *
133	johnpye	1312	FOR i IN mom_upwind4 CREATE
134	johnpye	1305	z_mombal2[i]: - dz/A * dmdot_dt[i]
135	jpye	1383	= (p[i]-p[i-1]) * backdiff for pressure *
136	johnpye	1305	+ dz * f/D/2 * rho[i] * vel[i]^2
137			+ (rho[i+1]vel[i+1]^2 + 6.rho[i]vel[i]^2 - 3.rho[i-1]vel[i-1]^2 - 2.rho[i-2]*vel[i-2]^2) / 6.;
138	johnpye	1278	END FOR;
139	johnpye	1312	FOR i IN mom_central CREATE
140	johnpye	1305	z_mombal1[i]: - dz/A * dmdot_dt[i]
141			= (p[i+1]-p[i-1]) / 2.
142			+ dz * f/D/2 * rho[i] * vel[i]^2
143			+ (rho[i+1]vel[i+1]^2 - rho[i-1]vel[i-1]^2) / 2.;
144			END FOR;
145	johnpye	1312	FOR i IN butfirst1 - mom_upwind4 - mom_central CREATE
146	johnpye	1305	z_mombal[i]: - dz/A * dmdot_dt[i]
147			= (p[i]-p[i-1])
148			+ dz * f/D/2 * rho[i] * vel[i]^2
149			+ (rho[i]vel[i]^2 - rho[i-1]vel[i-1]^2);
150			END FOR;
151	jpye	1383	*)
152	johnpye	1277
153	johnpye	1305	(* internal/external convection, and thermal mass of wall -- no spatial derivs here *)
154			FOR i IN butfirst1 CREATE
155	johnpye	1278	z_wall[i]: rho_wA_wc_w*dTw_dt[i] = qdot_s - qdot_l[i] - qdot_t[i];
156			z_loss[i]: qdot_l[i] = h_ext(1{PI}D_2)*(T_w[i] - T_amb);
157	johnpye	1305	z_trans[i]: qdot_t[i] = h_int(1{PI}D) *(T_w[i] - T[i]);
158	johnpye	1277	END FOR;
159
160	johnpye	1179	t IS_A time;
161	jpye	1383
162	johnpye	1179	METHODS
163	jpye	1383
164	johnpye	1278	METHOD bound_self;
165	johnpye	1305	vel[nodes].upper_bound := 100 {m/s};
166			qdot_l[butfirst1].lower_bound := 0 {W/m};
167			FOR i IN nodes DO
168	johnpye	1278	RUN node[i].bound_self;
169			END FOR;
170			END bound_self;
171	johnpye	1305	METHOD default;
172			(* these are initial guesses only; fixed parameters are overwritten by 'values' below *)
173	johnpye	1278	t := 0 {s};
174	johnpye	1305	FOR i IN nodes DO
175	johnpye	1312	T[i] := 298 {K};
176			vel[i] := 1 {m/s};
177	johnpye	1278	RUN node[i].default_self;
178			END FOR;
179	johnpye	1305	FOR i IN butfirst1 DO
180	johnpye	1278	drho_dt[i] := 0 {kg/m^3/s};
181	jpye	1383	(* dmdot_dt[i] := 0 {kg/s/s}; *)
182	johnpye	1278	du_dt[i] := 0 {kJ/kg/s};
183			dTw_dt[i] := 0 {K/s};
184			qdot_t[i] := 0 {W/m};
185			qdot_l[i] := 0 {W/m};
186	johnpye	1312	x[i] := x[1];
187	johnpye	1278	END FOR;
188	johnpye	1305	END default;
189	johnpye	1278	METHOD specify;
190			(* change to a proper steady-state problem, with fluid properties FREEd *)
191	johnpye	1305	FOR i IN nodes DO
192	johnpye	1278	RUN node[i].specify;
193			FIX dTw_dt[i]; FREE T_w[i];
194			END FOR;
195	johnpye	1305	FIX p[1];
196			FREE T[1];
197	johnpye	1278	FIX qdot_s;
198	johnpye	1305	FIX D, D_2;
199	johnpye	1278	FIX h_int, c_w, rho_w, h_ext;
200			FIX f;
201			(* FIX mu_f; *)
202			FIX T_amb;
203			(* fix derivatives to zero *)
204	johnpye	1311	FOR i IN butfirst1 DO
205	jpye	1383	FREE x[i]; FIX p[i]; FREE node[i].mdot;
206	johnpye	1305	FIX drho_dt[i]; FREE p[i];
207			FIX du_dt[i]; FREE T[i];
208	jpye	1383	(* FREE mdot[i]; FIX dmdot_dt[i]; *)
209	johnpye	1278	END FOR;
210			END specify;
211			METHOD values;
212	johnpye	1313	D := 0.06 {m};
213			D_2 := 0.07 {m};
214			T_amb := 298 {K};
215	jpye	1383	h_int := 1000 {W/m^2/K};
216			h_ext := 5 {W/m^2/K};
217	johnpye	1307	f := 0.03;
218	johnpye	1305	mdot[1] := 0.26 {kg/s};
219			p[1] := 10 {bar};
220	johnpye	1307	x[1] := 0.23;
221	johnpye	1313	rho_w := 1000 {kg/m^3};
222	jpye	1383	qdot_s := 100 {W/m^2} * D_2 * 10 * 11{m} / L;
223	johnpye	1311	FOR i IN butfirst1 DO
224	johnpye	1313	T_w[i] := 298 {K};
225	jpye	1383	(* dmdot_dt[i] := 0.0 {kg/s/s};*)
226	johnpye	1278	du_dt[i] := 0 {kJ/kg/s};
227	johnpye	1305	v[i] := 0.2 {L/kg};
228			rho[i] := 6 {kg/L};
229	johnpye	1278	node[i].dp_dT := +0.5 {kPa/K};
230	johnpye	1305	p[i] := 5 {bar};
231	johnpye	1278	END FOR;
232			END values;
233			METHOD on_load;
234			RUN configure_steady;
235			RUN ode_init;
236			END on_load;
237			(---------------- a physically sensible steady-state configuration-----------)
238			METHOD configure_steady;
239	jpye	1383	EXTERNAL defaultself_visit_childatoms(SELF);
240			EXTERNAL defaultself_visit_submodels(SELF);
241	johnpye	1278	RUN ClearAll;
242			RUN specify;
243			RUN bound_steady;
244			RUN values;
245			END configure_steady;
246			METHOD bound_steady;
247			RUN bound_self;
248	johnpye	1311	T_w[butfirst1].upper_bound := 1000 {K};
249	johnpye	1278	END bound_steady;
250			(------------------------- the dynamic problem ------------------------------)
251			METHOD configure_dynamic;
252	johnpye	1311	FOR i IN butfirst1 DO
253	johnpye	1305	FREE drho_dt[i]; FIX rho[i];
254	jpye	1383	(* FREE dmdot_dt[i]; FIX mdot[i]; *)
255	johnpye	1305	FREE du_dt[i]; FIX u[i];
256	johnpye	1278	FREE dTw_dt[i]; FIX T_w[i];
257	johnpye	1305	FREE x[i];
258			FREE T[i];
259	johnpye	1278	END FOR;
260			t := 0 {s};
261			END configure_dynamic;
262			METHOD free_states;
263	johnpye	1311	FOR i IN butfirst1 DO
264	johnpye	1305	FREE rho[i];
265	jpye	1383	(* FREE mdot[i]; *)
266	johnpye	1305	FREE u[i];
267	johnpye	1278	FREE T_w[i];
268			END FOR;
269			END free_states;
270			METHOD ode_init;
271			(* add the necessary meta data to allow solving with the integrator *)
272			t.ode_type := -1;
273	johnpye	1179
274	johnpye	1311	FOR i IN butfirst1 DO
275	johnpye	1305	drho_dt[i].ode_id := 4i; rho[i].ode_id := 4i;
276			drho_dt[i].ode_type := 2; rho[i].ode_type := 1;
277	johnpye	1179
278	jpye	1383	(* dmdot_dt[i].ode_id := 4i+1; mdot[i].ode_id := 4i+1;
279			dmdot_dt[i].ode_type := 2; mdot[i].ode_type := 1;*)
280	johnpye	1278
281	johnpye	1305	du_dt[i].ode_id := 4i+2; u[i].ode_id := 4i+2;
282			du_dt[i].ode_type := 2; u[i].ode_type := 1;
283	johnpye	1179
284	johnpye	1278	dTw_dt[i].ode_id := 4i+3; T_w[i].ode_id := 4i+3;
285			dTw_dt[i].ode_type := 2; T_w[i].ode_type := 1;
286			END FOR;
287	johnpye	1179
288	johnpye	1305	FOR i IN nodes DO
289	johnpye	1278	(* p[i].obs_id := 1 + 10i; )
290			(* x[i].obs_id := 2 + 10i; )
291	johnpye	1305	(* mdot[i].obs_id := 4 + 10i; )
292	johnpye	1279	(* node[i].h.obs_id := 3 + 10i; )
293	johnpye	1278	END FOR;
294	johnpye	1311	FOR i IN butfirst1 DO
295	johnpye	1278	(* qdot_t[i].obs_id := 3 + 10i; )
296			(* T_w[i].obs_id := 5 + 10i; )
297			(* T[i].obs_id := 6 + 10i;)
298			END FOR;
299	johnpye	1279
300	johnpye	1305	mdot[n].obs_id := 1;
301			x[n].obs_id := 1;
302	johnpye	1279	p[n].obs_id := 1;
303			vel[n].obs_id := 1;
304	jpye	1383	T_w[n].obs_id := 2;
305	johnpye	1279
306	johnpye	1278	END ode_init;
307			METHOD fix_outlet_quality;
308			FIX x[n];
309	johnpye	1305	FREE mdot[1];
310	johnpye	1278	END fix_outlet_quality;
311	johnpye	1179
312	jpye	1331	METHOD roots;
313			EXTERNAL roots(SELF);
314			END roots;
315
316	johnpye	1179	END dsgsat3;
317			ADD NOTES IN dsgsat2;
318			'QRSlv' iterationlimit {50}
319			END NOTES;