trunk/models/heatex.a4c

(*	ASCEND modelling environment
	Copyright (C) 1998, 2006, 2007 Carnegie Mellon University

	This program is free software; you can redistribute it and/or modify
	it under the terms of the GNU General Public License as published by
	the Free Software Foundation; either version 2, or (at your option)
	any later version.

	This program is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU General Public License for more details.
	
	You should have received a copy of the GNU General Public License
	along with this program.  If not, see <http://www.gnu.org/licenses/>.
*)
REQUIRE "atoms.a4l";
PROVIDE "heatex.a4c";
IMPORT "conopt";
IMPORT "lrslv";
IMPORT "cmslv";
(*
	The following the conditional model is discussed in the PhD thesis of
	Vicente Rico-Ramirez,
	https://pse.cheme.cmu.edu/ascend/ftp/pdfThesis/victhesis.pdf

	The original problem was presented by Joe Zaher

	Zaher, Conditional Modeling. Ph.D. Thesis, Carnegie Mellon University,
	Pittsburgh, PA. 1995.

	The problem represents the condensation of vapor because the heat
	transfer between the vapor and cooling water in heat exchanger. There 
	are two disjunctive statement which represent the possible condensation
	of vapor in the different sections of the heat exchanger. This example
	shows how we can represent a conditional model containing differential 
	equations. It represent a problem that we can solve with the ASCEND
	conditional solver, CMSlv.
*)

MODEL heatex;
	points                 IS_A set OF integer_constant;
	components             IS_A set OF symbol_constant;
	k_terms                IS_A set OF integer_constant;
	z[components]          IS_A fraction;
	x[points][components]  IS_A fraction;
	y[points][components]  IS_A fraction;
	Tcr[components]        IS_A critical_temperature;
	Pcr[components]        IS_A critical_pressure;
	Hv[components]         IS_A molar_energy_constant;
	B[components][k_terms] IS_A factor_constant;
	C[k_terms]             IS_A factor_constant;
	D[components][k_terms] IS_A molar_energy_constant;
	E[k_terms]             IS_A factor_constant;
	phi[points]            IS_A fraction;
	eta[points]            IS_A fraction;
	A                      IS_A area;
	Fc                     IS_A molar_rate;
	Fh                     IS_A molar_rate;
	P                      IS_A pressure;
	Tc[points]             IS_A temperature;
	Th[points]             IS_A temperature;
	Hc[points]             IS_A molar_energy;
	Hh[points]             IS_A molar_energy;
	dQ[points]             IS_A energy_rate;
	erscale                IS_A energy_rate;
	entscale               IS_A factor;
	bolphaq[1..2]          IS_A boolean_var;

	(* Boundaries *)
	CONDITIONAL
		cond1: SUM[x[0][i] | i IN components] + phi[0] >= 1.0;
		cond2: SUM[x[2][i] | i IN components] + phi[1] >= 1.0;
	END CONDITIONAL;

	bolphaq[1] == SATISFIED(cond1,1e-08);
	bolphaq[2] == SATISFIED(cond2,1e-08);

	(* Variant Equations *)
	sum0: SUM[x[0][i] | i IN components] = 1.0;
	frac0: phi[0] = 0.0;

	sum2: SUM[x[2][i] | i IN components] = 1.0;
	frac1: phi[1] = 0.0;

	p1: eta[1] = 0.5;
	p2: eta[2] = 0.5;

	sum1: SUM[x[1][i] | i IN components] = 1.0;
	p12: eta[1] = 0.0;

	(* Disjunctive Statements *)
	WHEN (bolphaq[1])
	CASE TRUE:
		USE sum0;
		USE sum1;
	CASE FALSE:
		USE frac0;
		USE p12;
	END WHEN;

	WHEN (bolphaq[2])
	CASE TRUE:
		USE sum2;
		USE p1;
	CASE FALSE:
		USE frac1;
		USE p2;
	END WHEN;

	(* Invariant Equations *)
	phi[2] = 0.0;
	phi[3] = 0.0;

	FOR i IN points CREATE
		FOR j IN components CREATE
			phi[i]*x[i][j] + (1 - phi[i])*y[i][j] = z[j];
		END FOR;
	END FOR;

	FOR i IN points CREATE
		FOR j IN components CREATE
			ln((y[i][j]*P)/(x[i][j]*Pcr[j])) = (Tcr[j]/Th[i])*(
				SUM[B[j][k] *(1 - Th[i]/Tcr[j])^C[k]
				| k IN k_terms] );
		END FOR;
	END FOR;

	FOR i IN points CREATE
		dQ[i]/erscale = 20{BTU/hour/ft^2/R}*(1 + 10*phi[i]) *
			A*(Th[i] - Tc[i])/erscale;
	END FOR;

	FOR i IN [0..2] CREATE
		Fh*(Hh[i+1] - Hh[i])/erscale = Fc*(Hc[i+1] - Hc[i])/erscale;

		Fc*(Hc[i+1] - Hc[i])/erscale = 0.5*(dQ[i] + dQ[i+1]) *
			(eta[i+1] - eta[i])/erscale;
	END FOR;

	FOR i IN points CREATE
		Hc[i]/entscale = 9720{BTU/lb_mole}*(Tc[i]/540{R} - 1)/entscale;
	END FOR;

	FOR j IN points CREATE
		Hh[j]/entscale = (SUM[ z[i] *(
			SUM[ D[i][k]*((Th[j]/540{R})^E[k] - 1)|k IN k_terms])
			- phi[j]*x[j][i]*Hv[i] | i IN components])/entscale;
	END FOR;

	(* Definition of sets *)
	components :== ['B','P5','H'];
	points :== [0..3];
	k_terms :== [1..4];

	(* Constants *)
	Tcr['B'] :== 765.4 {R};
	Tcr['P5'] :== 845.3 {R};
	Tcr['H'] :== 913.3 {R};
	Pcr['B'] :== 37.5 {atm};
	Pcr['P5'] :== 33.3 {atm};
	Pcr['H'] :== 29.3 {atm};
	Hv['B'] :== 9634 {BTU/lb_mole};
	Hv['P5'] :== 11088{BTU/lb_mole};
	Hv['H'] :== 12413 {BTU/lb_mole};
	C[1] :== 1.0;
	C[2] :== 1.5;
	C[3] :== 3.0;
	C[4] :== 6.0;
	E[1] :== 1.0;
	E[2] :== 2.0;
	E[3] :== 3.0;
	E[4] :== 4.0;
	B['B'][1]:== -6.88709;
	B['B'][2]:== 1.15157;
	B['B'][3]:== -1.99873;
	B['B'][4]:== -3.13003;
	B['P5'][1]:== -7.28936;
	B['P5'][2]:== 1.53679;
	B['P5'][3]:== -3.08367;
	B['P5'][4]:== -1.02456;
	B['H'][1]:== -7.46765;
	B['H'][2]:== 1.44211;
	B['H'][3]:== -3.28222;
	B['H'][4]:== -2.50941;
	D['B'][1]:== 1224 {BTU/lb_mole};
	D['B'][2]:== 6410 {BTU/lb_mole};
	D['B'][3]:== -429 {BTU/lb_mole};
	D['B'][4]:== -2 {BTU/lb_mole};
	D['P5'][1]:== -468 {BTU/lb_mole};
	D['P5'][2]:== 9428 {BTU/lb_mole};
	D['P5'][3]:== -998 {BTU/lb_mole};
	D['P5'][4]:== 46 {BTU/lb_mole};
	D['H'][1]:== -569 {BTU/lb_mole};
	D['H'][2]:== 11260 {BTU/lb_mole};
	D['H'][3]:== -1207 {BTU/lb_mole};
	D['H'][4]:== 57 {BTU/lb_mole};

METHODS
    METHOD default_self;
    END default_self;

    METHOD specify;
		FOR j IN components DO
		    FIX z[j];
		    END FOR;
		FIX A;
		FIX P;
		FIX Fc;
		FIX Fh;
		    FIX erscale;
		    FIX entscale;
		FIX eta[0];
		FIX eta[3];
		FIX Tc[0];
		FIX Th[3];
    END specify;

    METHOD values;
		(* fixed *)
		P := 1.0{atm};
		z['B'] := 0.2;
		z['P5'] := 0.5;
		z['H'] := 0.3;
		Fh := 250.0{lb_mole/hour};
		Tc[0] := 540.0{R};
		Th[3] := 600.0{R};
		erscale := 1.0e6{BTU/hour};
		entscale := 1000.0 ;

		(* initial guess *)
		A := 379.12{ft^2};
		Fc := 1104.13{lb_mole/hour};
		x[0]['B'] := 0.04183;
		x[0]['P5'] := 0.32735;
		x[0]['H'] := 0.58922;
		x[1]['B'] := 0.04183;
		x[1]['P5'] := 0.32735;
		x[1]['H'] := 0.58922;
		x[2]['B'] := 0.03558;
		x[2]['P5'] := 0.27044;
		x[2]['H'] := 0.47211;
		x[3]['B'] := 0.03153;
		x[3]['P5'] := 0.23457;
		x[3]['H'] := 0.40038;
		y[0]['B'] := 0.2;
		y[0]['P5'] := 0.5;
		y[0]['H'] := 0.3;
		y[1]['B'] := 0.2;
		y[1]['P5'] := 0.5;
		y[1]['H'] := 0.3;
		y[2]['B'] := 0.2;
		y[2]['P5'] := 0.5;
		y[2]['H'] := 0.3;
		y[3]['B'] := 0.2;
		y[3]['P5'] := 0.5;
		y[3]['H'] := 0.3;
		phi[0] := 0.0;
		phi[1] := 0.0;
		phi[2] := 0.0;
		phi[3] := 0.0;
		eta[0] := 0.0;
		eta[1] := 0.0;
		eta[2] := 0.5;
		eta[3] := 1.0;
		Tc[1] := 540.000{R};
		Tc[2] := 559.820{R};
		Tc[3] := 575.386{R};
		Th[0] := 579.932{R};
		Th[1] := 579.932{R};
		Th[2] := 591.242{R};
		Hc[0] := 0000.0{BTU/lb_mole};
		Hc[1] := 0000.0{BTU/lb_mole};
		Hc[2] := 356.768{BTU/lb_mole};
		Hc[3] := 636.948{BTU/lb_mole};
		Hh[0] := 1212.93{BTU/lb_mole};
		Hh[1] := 1212.93{BTU/lb_mole};
		Hh[2] := 1569.7{BTU/lb_mole};
		Hh[3] := 1849.88{BTU/lb_mole};
		dQ[0] := 199658{BTU/hour};
		dQ[1] := 199658{BTU/hour};
		dQ[2] := 157110{BTU/hour};
		dQ[3] := 123070{BTU/hour};
		eta[1].upper_bound := 0.5;
		eta[2].lower_bound := 0.5;
		eta[2].upper_bound := 1.0;

		(* Initial value for booleans *)
		bolphaq[1] := SATISFIED(cond1,1e-08);
		bolphaq[2] := SATISFIED(cond2,1e-08);
    END values;

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

	METHOD self_test;
		(* nothing tested here as yet *)
	END self_test;

END heatex;
(* :ex: set ts=4: *)
1	johnpye	1155	(* ASCEND modelling environment
2	johnpye	1259	Copyright (C) 1998, 2006, 2007 Carnegie Mellon University
3	johnpye	1154
4			This program is free software; you can redistribute it and/or modify
5			it under the terms of the GNU General Public License as published by
6			the Free Software Foundation; either version 2, or (at your option)
7			any later version.
8
9			This program is distributed in the hope that it will be useful,
10			but WITHOUT ANY WARRANTY; without even the implied warranty of
11			MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12			GNU General Public License for more details.
13
14			You should have received a copy of the GNU General Public License
15	jpye	2649	along with this program. If not, see <http://www.gnu.org/licenses/>.
16	johnpye	1154	*)
17	aw0a	1	REQUIRE "atoms.a4l";
18			PROVIDE "heatex.a4c";
19	jpye	3308	IMPORT "conopt";
20			IMPORT "lrslv";
21			IMPORT "cmslv";
22	aw0a	1	(*
23	johnpye	1154	The following the conditional model is discussed in the PhD thesis of
24			Vicente Rico-Ramirez,
25			https://pse.cheme.cmu.edu/ascend/ftp/pdfThesis/victhesis.pdf
26	aw0a	1
27	johnpye	1154	The original problem was presented by Joe Zaher
28	aw0a	1
29	johnpye	1154	Zaher, Conditional Modeling. Ph.D. Thesis, Carnegie Mellon University,
30			Pittsburgh, PA. 1995.
31	aw0a	1
32	johnpye	1154	The problem represents the condensation of vapor because the heat
33			transfer between the vapor and cooling water in heat exchanger. There
34			are two disjunctive statement which represent the possible condensation
35			of vapor in the different sections of the heat exchanger. This example
36			shows how we can represent a conditional model containing differential
37			equations. It represent a problem that we can solve with the ASCEND
38			conditional solver, CMSlv.
39	aw0a	1	*)
40
41			MODEL heatex;
42	johnpye	1257	points IS_A set OF integer_constant;
43			components IS_A set OF symbol_constant;
44			k_terms IS_A set OF integer_constant;
45			z[components] IS_A fraction;
46			x[points][components] IS_A fraction;
47			y[points][components] IS_A fraction;
48			Tcr[components] IS_A critical_temperature;
49			Pcr[components] IS_A critical_pressure;
50			Hv[components] IS_A molar_energy_constant;
51			B[components][k_terms] IS_A factor_constant;
52			C[k_terms] IS_A factor_constant;
53			D[components][k_terms] IS_A molar_energy_constant;
54			E[k_terms] IS_A factor_constant;
55			phi[points] IS_A fraction;
56			eta[points] IS_A fraction;
57			A IS_A area;
58			Fc IS_A molar_rate;
59			Fh IS_A molar_rate;
60			P IS_A pressure;
61			Tc[points] IS_A temperature;
62			Th[points] IS_A temperature;
63			Hc[points] IS_A molar_energy;
64			Hh[points] IS_A molar_energy;
65			dQ[points] IS_A energy_rate;
66			erscale IS_A energy_rate;
67			entscale IS_A factor;
68			bolphaq[1..2] IS_A boolean_var;
69	aw0a	1
70	johnpye	1154	(* Boundaries *)
71			CONDITIONAL
72			cond1: SUM[x[0][i] \| i IN components] + phi[0] >= 1.0;
73			cond2: SUM[x[2][i] \| i IN components] + phi[1] >= 1.0;
74			END CONDITIONAL;
75	aw0a	1
76	johnpye	1154	bolphaq[1] == SATISFIED(cond1,1e-08);
77			bolphaq[2] == SATISFIED(cond2,1e-08);
78	aw0a	1
79	johnpye	1154	(* Variant Equations *)
80			sum0: SUM[x[0][i] \| i IN components] = 1.0;
81			frac0: phi[0] = 0.0;
82	aw0a	1
83	johnpye	1154	sum2: SUM[x[2][i] \| i IN components] = 1.0;
84			frac1: phi[1] = 0.0;
85	aw0a	1
86	johnpye	1154	p1: eta[1] = 0.5;
87			p2: eta[2] = 0.5;
88	aw0a	1
89	johnpye	1154	sum1: SUM[x[1][i] \| i IN components] = 1.0;
90			p12: eta[1] = 0.0;
91	aw0a	1
92	johnpye	1154	(* Disjunctive Statements *)
93			WHEN (bolphaq[1])
94	aw0a	1	CASE TRUE:
95	johnpye	1154	USE sum0;
96			USE sum1;
97	aw0a	1	CASE FALSE:
98	johnpye	1154	USE frac0;
99			USE p12;
100			END WHEN;
101	aw0a	1
102	johnpye	1154	WHEN (bolphaq[2])
103	aw0a	1	CASE TRUE:
104	johnpye	1154	USE sum2;
105			USE p1;
106	aw0a	1	CASE FALSE:
107	johnpye	1154	USE frac1;
108			USE p2;
109			END WHEN;
110	aw0a	1
111	johnpye	1154	(* Invariant Equations *)
112			phi[2] = 0.0;
113			phi[3] = 0.0;
114	aw0a	1
115	johnpye	1154	FOR i IN points CREATE
116			FOR j IN components CREATE
117			phi[i]x[i][j] + (1 - phi[i])y[i][j] = z[j];
118			END FOR;
119			END FOR;
120	aw0a	1
121	johnpye	1154	FOR i IN points CREATE
122			FOR j IN components CREATE
123			ln((y[i][j]P)/(x[i][j]Pcr[j])) = (Tcr[j]/Th[i])*(
124			SUM[B[j][k] *(1 - Th[i]/Tcr[j])^C[k]
125			\| k IN k_terms] );
126			END FOR;
127			END FOR;
128	aw0a	1
129	johnpye	1154	FOR i IN points CREATE
130			dQ[i]/erscale = 20{BTU/hour/ft^2/R}(1 + 10phi[i]) *
131			A*(Th[i] - Tc[i])/erscale;
132			END FOR;
133	aw0a	1
134	johnpye	1154	FOR i IN [0..2] CREATE
135			Fh(Hh[i+1] - Hh[i])/erscale = Fc(Hc[i+1] - Hc[i])/erscale;
136	aw0a	1
137	johnpye	1154	Fc(Hc[i+1] - Hc[i])/erscale = 0.5(dQ[i] + dQ[i+1]) *
138			(eta[i+1] - eta[i])/erscale;
139			END FOR;
140	aw0a	1
141	johnpye	1154	FOR i IN points CREATE
142			Hc[i]/entscale = 9720{BTU/lb_mole}*(Tc[i]/540{R} - 1)/entscale;
143			END FOR;
144	aw0a	1
145	johnpye	1154	FOR j IN points CREATE
146			Hh[j]/entscale = (SUM[ z[i] *(
147			SUM[ D[i][k]*((Th[j]/540{R})^E[k] - 1)\|k IN k_terms])
148			- phi[j]x[j][i]Hv[i] \| i IN components])/entscale;
149			END FOR;
150	aw0a	1
151	johnpye	1154	(* Definition of sets *)
152			components :== ['B','P5','H'];
153			points :== [0..3];
154			k_terms :== [1..4];
155	aw0a	1
156	johnpye	1154	(* Constants *)
157			Tcr['B'] :== 765.4 {R};
158			Tcr['P5'] :== 845.3 {R};
159			Tcr['H'] :== 913.3 {R};
160			Pcr['B'] :== 37.5 {atm};
161			Pcr['P5'] :== 33.3 {atm};
162			Pcr['H'] :== 29.3 {atm};
163			Hv['B'] :== 9634 {BTU/lb_mole};
164			Hv['P5'] :== 11088{BTU/lb_mole};
165			Hv['H'] :== 12413 {BTU/lb_mole};
166			C[1] :== 1.0;
167			C[2] :== 1.5;
168			C[3] :== 3.0;
169			C[4] :== 6.0;
170			E[1] :== 1.0;
171			E[2] :== 2.0;
172			E[3] :== 3.0;
173			E[4] :== 4.0;
174			B['B'][1]:== -6.88709;
175			B['B'][2]:== 1.15157;
176			B['B'][3]:== -1.99873;
177			B['B'][4]:== -3.13003;
178			B['P5'][1]:== -7.28936;
179			B['P5'][2]:== 1.53679;
180			B['P5'][3]:== -3.08367;
181			B['P5'][4]:== -1.02456;
182			B['H'][1]:== -7.46765;
183			B['H'][2]:== 1.44211;
184			B['H'][3]:== -3.28222;
185			B['H'][4]:== -2.50941;
186			D['B'][1]:== 1224 {BTU/lb_mole};
187			D['B'][2]:== 6410 {BTU/lb_mole};
188			D['B'][3]:== -429 {BTU/lb_mole};
189			D['B'][4]:== -2 {BTU/lb_mole};
190			D['P5'][1]:== -468 {BTU/lb_mole};
191			D['P5'][2]:== 9428 {BTU/lb_mole};
192			D['P5'][3]:== -998 {BTU/lb_mole};
193			D['P5'][4]:== 46 {BTU/lb_mole};
194			D['H'][1]:== -569 {BTU/lb_mole};
195			D['H'][2]:== 11260 {BTU/lb_mole};
196			D['H'][3]:== -1207 {BTU/lb_mole};
197			D['H'][4]:== 57 {BTU/lb_mole};
198	aw0a	1
199	johnpye	1154	METHODS
200	aw0a	1	METHOD default_self;
201			END default_self;
202
203			METHOD specify;
204	johnpye	1154	FOR j IN components DO
205			FIX z[j];
206			END FOR;
207			FIX A;
208			FIX P;
209			FIX Fc;
210			FIX Fh;
211			FIX erscale;
212			FIX entscale;
213			FIX eta[0];
214			FIX eta[3];
215			FIX Tc[0];
216			FIX Th[3];
217	aw0a	1	END specify;
218
219			METHOD values;
220	johnpye	1154	(* fixed *)
221			P := 1.0{atm};
222			z['B'] := 0.2;
223			z['P5'] := 0.5;
224			z['H'] := 0.3;
225			Fh := 250.0{lb_mole/hour};
226			Tc[0] := 540.0{R};
227			Th[3] := 600.0{R};
228			erscale := 1.0e6{BTU/hour};
229			entscale := 1000.0 ;
230	aw0a	1
231	johnpye	1154	(* initial guess *)
232			A := 379.12{ft^2};
233			Fc := 1104.13{lb_mole/hour};
234			x[0]['B'] := 0.04183;
235			x[0]['P5'] := 0.32735;
236			x[0]['H'] := 0.58922;
237			x[1]['B'] := 0.04183;
238			x[1]['P5'] := 0.32735;
239			x[1]['H'] := 0.58922;
240			x[2]['B'] := 0.03558;
241			x[2]['P5'] := 0.27044;
242			x[2]['H'] := 0.47211;
243			x[3]['B'] := 0.03153;
244			x[3]['P5'] := 0.23457;
245			x[3]['H'] := 0.40038;
246			y[0]['B'] := 0.2;
247			y[0]['P5'] := 0.5;
248			y[0]['H'] := 0.3;
249			y[1]['B'] := 0.2;
250			y[1]['P5'] := 0.5;
251			y[1]['H'] := 0.3;
252			y[2]['B'] := 0.2;
253			y[2]['P5'] := 0.5;
254			y[2]['H'] := 0.3;
255			y[3]['B'] := 0.2;
256			y[3]['P5'] := 0.5;
257			y[3]['H'] := 0.3;
258			phi[0] := 0.0;
259			phi[1] := 0.0;
260			phi[2] := 0.0;
261			phi[3] := 0.0;
262			eta[0] := 0.0;
263			eta[1] := 0.0;
264			eta[2] := 0.5;
265			eta[3] := 1.0;
266			Tc[1] := 540.000{R};
267			Tc[2] := 559.820{R};
268			Tc[3] := 575.386{R};
269			Th[0] := 579.932{R};
270			Th[1] := 579.932{R};
271			Th[2] := 591.242{R};
272			Hc[0] := 0000.0{BTU/lb_mole};
273			Hc[1] := 0000.0{BTU/lb_mole};
274			Hc[2] := 356.768{BTU/lb_mole};
275			Hc[3] := 636.948{BTU/lb_mole};
276			Hh[0] := 1212.93{BTU/lb_mole};
277			Hh[1] := 1212.93{BTU/lb_mole};
278			Hh[2] := 1569.7{BTU/lb_mole};
279			Hh[3] := 1849.88{BTU/lb_mole};
280			dQ[0] := 199658{BTU/hour};
281			dQ[1] := 199658{BTU/hour};
282			dQ[2] := 157110{BTU/hour};
283			dQ[3] := 123070{BTU/hour};
284			eta[1].upper_bound := 0.5;
285			eta[2].lower_bound := 0.5;
286			eta[2].upper_bound := 1.0;
287	aw0a	1
288	johnpye	1154	(* Initial value for booleans *)
289			bolphaq[1] := SATISFIED(cond1,1e-08);
290			bolphaq[2] := SATISFIED(cond2,1e-08);
291	aw0a	1	END values;
292
293	johnpye	1154	METHOD on_load;
294			RUN default_self;
295			RUN reset;
296			RUN values;
297			END on_load;
298
299	johnpye	1257	METHOD self_test;
300			(* nothing tested here as yet *)
301			END self_test;
302
303	aw0a	1	END heatex;
304	johnpye	1154	(* :ex: set ts=4: *)