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 :== 7;(* 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 :== 10 {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;

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

	(* 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;
	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 :== [];
	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;
	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;

	(* momentum conservation *)
	mom_upwind4 :== [];
	mom_central :== central;
	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 *)
		D := 0.06 {m};
		D_2 := 0.07 {m};
		A_w := 0.25{PI}*D_2^2 -0.25{PI}*D^2;
		A := 1 {m^2};
		T_amb := 298 {K};
		f := 0.0;
		h_ext := 10 {W/m^2/K};
		h_int := 100 {W/m^2/K};
		qdot_s := 1000 {W/m};
		rho_w := 1000 {kg/m^3};
		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
			T_w[i] := 298 {K};
			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];
			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;
		h_int := 100 {W/m^2/K};
		h_ext := 20 {W/m^2/K};
		f := 0.03;
		mdot[1] := 0.26 {kg/s};
		p[1]    := 10 {bar};
		x[1]    := 0.23;
		qdot_s := 1000 {W/m^2} * D_2 * 10;
		FOR i IN butfirst1 DO
			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;
		RUN default_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;

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

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