models/johnpye/dsg.a4c

REQUIRE "atoms.a4l";
REQUIRE "johnpye/thermo_types.a4c";
IMPORT "freesteam";
IMPORT "dsg";

(*
	This model requires the 'freesteam' library be installed on your system.
	You must compile it with the 'ASCEND_CONFIG=`which ascend-config`' flag
	to ensure that the required external library is installed where ASCEND
	can find it. Alternatively, you can specify the location by adding to the
	ASCENDLIBRARY path variable. See http://freesteam.sf.net/

	This model also requires the 'dsg' library for Direct Steam Generation
	calculations, available on request from http://pye.dyndns.org/

	For example (assuming you already have the source code for the above):
	export $ASCENDLIBRARY=$HOME/src/ascend/models
	cd $HOME/src/freesteam
	scons -j2 ascend
	export ASCENDLIBRARY=$HOME/src/freesteam/ascend:$ASCENDLIBRARY
	cd $HOME/src/dsg-transient
	scons -j2 ascend
	export ASCENDLIBRARY=$HOME/src/dsg-transient/ascend:$ASCENDLIBRARY
	$HOME/src/ascend/pygtk/ascdev johnpye/dsg.a4c	
*)

MODEL dsg;
	(* temporal derivatives *)
	drho_dt IS_A density_rate;
	dmdot_dt IS_A mass_rate_rate;
	drhou_dt IS_A solver_var;
	dTw_dt IS_A temperature_rate;

	(* spatial derivatives *)
	dmdot_dz IS_A mass_rate_per_length;
	dmdoth_dz IS_A solver_var;
	dekdot_dz IS_A solver_var;
	dp_dz IS_A pressure_per_length;
	drhovel2_dz IS_A pressure_per_length;

	(* wall properties *)
	rhow IS_A mass_density;
	D2 IS_A distance;
	cw IS_A specific_heat_capacity;
	Aw IS_A area;
	hw IS_A heat_transfer_coefficient;

	Aw = 1{PI}*(D2^2 - D^2)/4;

	(* conservation equations *)
	massbal: drho_dt = -1/A * dmdot_dz;
	mombal:  1/A * dmdot_dt = -dp_dz - f/D/2*rho*vel^2 - drhovel2_dz;
	enerbal: drhou_dt = 1/A * ( qdott - dmdoth_dz + mdot * dekdot_dz );
	pipebal: dTw_dt = 1/rhow/Aw/cw * (qdots - qdotl - qdott);

	rho IS_A mass_density;
	mdot IS_A mass_rate;
	p IS_A pressure;
	f IS_A factor; (* pipe friction factor *)
	u IS_A specific_energy;
	h IS_A specific_enthalpy;
	vel IS_A speed;
	qdott, qdotl,qdots IS_A power_per_length;
	mu IS_A viscosity;
	v IS_A specific_volume;
	T, T1, Tw, T2, Tamb IS_A temperature;
	ekdot IS_A specific_power;

	rho * v = 1;
	ekdot = 0.5 * rho * vel^2;

	thermo_props: iapws97_uvTxmu_ph(
		p,h : INPUT;
		u,v,T,x,mu : OUTPUT
	);

	vel = rho*mdot/A;

	Re IS_A factor;
	Re_rel: Re = rho*vel*D/mu;
	eps_on_D IS_A factor;
	eps_on_D_rel: eps_on_D = eps/D;

	D IS_A distance;
	eps IS_A distance;
	A IS_A area;
	A = 1{PI}* D^2 / 4;

	friction: dsg_friction_factor(
		Re, eps_on_D : INPUT;
		f : OUTPUT
	);

	phi2 IS_A factor;
	x IS_A factor;
	twophmult: dsg_phi2_martinelli_nelson(
		p,x : INPUT;
		phi2 : OUTPUT
	);
	(* TODO: compute LIQUID ONLY reynolds number when in the two-phase region...
	that's a bit tricky and needs conditional modelling, so maybe try putting 
	it in the external code. *)

	cavity_losses: dsg_ext_heat_loss(
		T2,Tamb,D2,hw : INPUT;
		qdotl : OUTPUT
	);

	Tw = (T1 + T2)/2;

	mdoth IS_A energy_rate;
	mdoth = mdot * h;

	rhovel2 IS_A solver_var;
	rhovel2 = rho * vel^2;

	(* internal heat transfer *)

	h1 IS_A heat_transfer_coefficient;

	qdott = h1 * 1{PI} * D * (T1 - T);

	twophconv: dsg_conv_coeff_two_phase(
		p, h, mdot, D, f : INPUT;
		h1 : OUTPUT
	);

METHODS

METHOD specify;
	RUN fix_design;
	RUN fix_states;
	RUN fix_spatials;
END specify;

METHOD fix_states;
	FIX T1, T2, p, h, mdot;
END fix_states;

METHOD fix_design;
	(* design parameters, geometry, materials of construction *)
	FIX rhow,D2,cw,hw;
	FIX D, eps;
	FIX qdots;
	FIX Tamb;
END fix_design;

METHOD fix_spatials;
	(* for use when tinkering about with a single node only *)
	FIX dmdot_dz, dmdoth_dz, dekdot_dz, dp_dz, drhovel2_dz;
END fix_spatials;

METHOD fix_temporals;
	(* for steady-state solution *)
	FIX drho_dt;
	FIX dmdot_dt;
	FIX drhou_dt;
	FIX dTw_dt;
END fix_temporals;

METHOD values;
	(* design parameters *)
	D := 60 {mm};
	D2 := 70 {mm};
	eps := 0.05 {mm};
	rhow := 7.8 {g/cm^3};
	cw := 0.47 {J/g/K};
	hw := 10 {W/m^2/K};
	Tamb := 300 {K};

	(* states *)
	T1 := 600 {K};
	T2 := 500 {K};
	p := 10 {bar};
	h := 4000 {kJ/kg};
	mdot := 0.001 {kg/s};

	(* spatial derivs *)
	dp_dz := -500 {Pa/m};
	dmdot_dz := 0 {kg/s/m};
	drhovel2_dz := (0 {kg/m^3}) * (0 {m/s})^2 / (1 {m});
	dekdot_dz := 0 {W/m};
	dmdoth_dz := dmdot_dz * 0 {kJ/kg};

	(* derivative variables *)
	drho_dt := 0 {kg/m^3/s};
	dmdot_dt := 0 {kg/s/s};
	drhou_dt := 0 {kg/m^3*kJ/kg/s};
	dTw_dt := 0 {K/s};

	(* bounds *)
	f.lower_bound := 0.008;
	v.lower_bound := 0.00999 {m^3/kg};
	v.upper_bound := 1/(0.0202 {kg/m^3});
	u.lower_bound := 0 {kJ/kg};
	u.upper_bound := 3663 {kJ/kg};

	(* starting guesses *)
	Re := 10000;
END values;

METHOD on_load;
	RUN reset;
	RUN values;
END on_load;

END dsg;
1	REQUIRE "atoms.a4l";
2	REQUIRE "johnpye/thermo_types.a4c";
3	IMPORT "freesteam";
4	IMPORT "dsg";
5
6	(*
7	This model requires the 'freesteam' library be installed on your system.
8	You must compile it with the 'ASCEND_CONFIG=`which ascend-config`' flag
9	to ensure that the required external library is installed where ASCEND
10	can find it. Alternatively, you can specify the location by adding to the
11	ASCENDLIBRARY path variable. See http://freesteam.sf.net/
12
13	This model also requires the 'dsg' library for Direct Steam Generation
14	calculations, available on request from http://pye.dyndns.org/
15
16	For example (assuming you already have the source code for the above):
17	export $ASCENDLIBRARY=$HOME/src/ascend/models
18	cd $HOME/src/freesteam
19	scons -j2 ascend
20	export ASCENDLIBRARY=$HOME/src/freesteam/ascend:$ASCENDLIBRARY
21	cd $HOME/src/dsg-transient
22	scons -j2 ascend
23	export ASCENDLIBRARY=$HOME/src/dsg-transient/ascend:$ASCENDLIBRARY
24	$HOME/src/ascend/pygtk/ascdev johnpye/dsg.a4c
25	*)
26
27	MODEL dsg;
28	(* temporal derivatives *)
29	drho_dt IS_A density_rate;
30	dmdot_dt IS_A mass_rate_rate;
31	drhou_dt IS_A solver_var;
32	dTw_dt IS_A temperature_rate;
33
34	(* spatial derivatives *)
35	dmdot_dz IS_A mass_rate_per_length;
36	dmdoth_dz IS_A solver_var;
37	dekdot_dz IS_A solver_var;
38	dp_dz IS_A pressure_per_length;
39	drhovel2_dz IS_A pressure_per_length;
40
41	(* wall properties *)
42	rhow IS_A mass_density;
43	D2 IS_A distance;
44	cw IS_A specific_heat_capacity;
45	Aw IS_A area;
46	hw IS_A heat_transfer_coefficient;
47
48	Aw = 1{PI}*(D2^2 - D^2)/4;
49
50	(* conservation equations *)
51	massbal: drho_dt = -1/A * dmdot_dz;
52	mombal: 1/A * dmdot_dt = -dp_dz - f/D/2rhovel^2 - drhovel2_dz;
53	enerbal: drhou_dt = 1/A * ( qdott - dmdoth_dz + mdot * dekdot_dz );
54	pipebal: dTw_dt = 1/rhow/Aw/cw * (qdots - qdotl - qdott);
55
56	rho IS_A mass_density;
57	mdot IS_A mass_rate;
58	p IS_A pressure;
59	f IS_A factor; (* pipe friction factor *)
60	u IS_A specific_energy;
61	h IS_A specific_enthalpy;
62	vel IS_A speed;
63	qdott, qdotl,qdots IS_A power_per_length;
64	mu IS_A viscosity;
65	v IS_A specific_volume;
66	T, T1, Tw, T2, Tamb IS_A temperature;
67	ekdot IS_A specific_power;
68
69	rho * v = 1;
70	ekdot = 0.5 * rho * vel^2;
71
72	thermo_props: iapws97_uvTxmu_ph(
73	p,h : INPUT;
74	u,v,T,x,mu : OUTPUT
75	);
76
77	vel = rho*mdot/A;
78
79	Re IS_A factor;
80	Re_rel: Re = rhovelD/mu;
81	eps_on_D IS_A factor;
82	eps_on_D_rel: eps_on_D = eps/D;
83
84	D IS_A distance;
85	eps IS_A distance;
86	A IS_A area;
87	A = 1{PI}* D^2 / 4;
88
89	friction: dsg_friction_factor(
90	Re, eps_on_D : INPUT;
91	f : OUTPUT
92	);
93
94	phi2 IS_A factor;
95	x IS_A factor;
96	twophmult: dsg_phi2_martinelli_nelson(
97	p,x : INPUT;
98	phi2 : OUTPUT
99	);
100	(* TODO: compute LIQUID ONLY reynolds number when in the two-phase region...
101	that's a bit tricky and needs conditional modelling, so maybe try putting
102	it in the external code. *)
103
104	cavity_losses: dsg_ext_heat_loss(
105	T2,Tamb,D2,hw : INPUT;
106	qdotl : OUTPUT
107	);
108
109	Tw = (T1 + T2)/2;
110
111	mdoth IS_A energy_rate;
112	mdoth = mdot * h;
113
114	rhovel2 IS_A solver_var;
115	rhovel2 = rho * vel^2;
116
117	(* internal heat transfer *)
118
119	h1 IS_A heat_transfer_coefficient;
120
121	qdott = h1 * 1{PI} * D * (T1 - T);
122
123	twophconv: dsg_conv_coeff_two_phase(
124	p, h, mdot, D, f : INPUT;
125	h1 : OUTPUT
126	);
127
128	METHODS
129
130	METHOD specify;
131	RUN fix_design;
132	RUN fix_states;
133	RUN fix_spatials;
134	END specify;
135
136	METHOD fix_states;
137	FIX T1, T2, p, h, mdot;
138	END fix_states;
139
140	METHOD fix_design;
141	(* design parameters, geometry, materials of construction *)
142	FIX rhow,D2,cw,hw;
143	FIX D, eps;
144	FIX qdots;
145	FIX Tamb;
146	END fix_design;
147
148	METHOD fix_spatials;
149	(* for use when tinkering about with a single node only *)
150	FIX dmdot_dz, dmdoth_dz, dekdot_dz, dp_dz, drhovel2_dz;
151	END fix_spatials;
152
153	METHOD fix_temporals;
154	(* for steady-state solution *)
155	FIX drho_dt;
156	FIX dmdot_dt;
157	FIX drhou_dt;
158	FIX dTw_dt;
159	END fix_temporals;
160
161	METHOD values;
162	(* design parameters *)
163	D := 60 {mm};
164	D2 := 70 {mm};
165	eps := 0.05 {mm};
166	rhow := 7.8 {g/cm^3};
167	cw := 0.47 {J/g/K};
168	hw := 10 {W/m^2/K};
169	Tamb := 300 {K};
170
171	(* states *)
172	T1 := 600 {K};
173	T2 := 500 {K};
174	p := 10 {bar};
175	h := 4000 {kJ/kg};
176	mdot := 0.001 {kg/s};
177
178	(* spatial derivs *)
179	dp_dz := -500 {Pa/m};
180	dmdot_dz := 0 {kg/s/m};
181	drhovel2_dz := (0 {kg/m^3}) * (0 {m/s})^2 / (1 {m});
182	dekdot_dz := 0 {W/m};
183	dmdoth_dz := dmdot_dz * 0 {kJ/kg};
184
185	(* derivative variables *)
186	drho_dt := 0 {kg/m^3/s};
187	dmdot_dt := 0 {kg/s/s};
188	drhou_dt := 0 {kg/m^3*kJ/kg/s};
189	dTw_dt := 0 {K/s};
190
191	(* bounds *)
192	f.lower_bound := 0.008;
193	v.lower_bound := 0.00999 {m^3/kg};
194	v.upper_bound := 1/(0.0202 {kg/m^3});
195	u.lower_bound := 0 {kJ/kg};
196	u.upper_bound := 3663 {kJ/kg};
197
198	(* starting guesses *)
199	Re := 10000;
200	END values;
201
202	METHOD on_load;
203	RUN reset;
204	RUN values;
205	END on_load;
206
207	END dsg;