test/bug567/brayton.a4c

(*	ASCEND modelling environment
	Copyright (C) 2007, 2008, 2009, 2010 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/>.
*)(*
	This file contains some models of Brayton engines and associated cycles,
	following the development of Çengel & Boles 'Thermodynamcs: An Engineering
	Approach, 6th Ed, McGraw-Hill, 2008.

	Author: John Pye
*)

REQUIRE "atoms.a4l";
REQUIRE "johnpye/thermo_types.a4c";
REQUIRE "johnpye/airprops.a4c";
IMPORT "sensitivity/solve";

(* first some models of air as an ideal gas *)

MODEL ideal_gas_base;
	M IS_A molar_weight_constant;	
	c_p IS_A specific_heat_capacity;
	s IS_A specific_entropy;
	h IS_A specific_enthalpy;
	v IS_A specific_volume;
	T IS_A temperature;
	p IS_A pressure;
	R IS_A specific_gas_constant;

	R :== 1{GAS_C} / M;
	p * v = R * T;

METHODS
METHOD bound_self;
	s.lower_bound := -5 {kJ/kg/K};
END bound_self;
END ideal_gas_base;

(*
	Ideal air assuming ideal gas and constant cp.
*)
MODEL simple_ideal_air
REFINES ideal_gas_base;

	M :== 28.958600656 {kg/kmol};

	c_p = 1.005 {kJ/kg/K};

	T_ref IS_A temperature_constant;
	p_ref IS_A pressure_constant;

	s = c_p * ln(T / T_ref) - R * ln(p / p_ref);
	h = c_p * (T - T_ref);

	T_ref :== 273.15 {K};
	p_ref :== 1 {bar};

METHODS
METHOD on_load;
	RUN ClearAll;
	RUN bound_self;
	FIX T, p;
	T := 300 {K};
	p := 1 {bar};
END on_load;	
END simple_ideal_air;

(*
	Ideal air, using a quartic polynomial for c_p as given in Moran & Shapiro
	4th Ed.
*)
MODEL ideal_air_poly REFINES ideal_gas_base;
	M :== 28.958600656 {kg/kmol};

	a[0..4] IS_A real_constant;
	a[0] :== 3.653;
	a[1] :== -1.337e-3{1/K};
	a[2] :== 3.294e-6{1/K^2};
	a[3] :== -1.913e-9{1/K^3};
	a[4] :== 0.2763e-12{1/K^4};

	T_ref IS_A temperature_constant;
	p_ref IS_A pressure_constant;
	h_ref IS_A real_constant;
	s_ref IS_A real_constant;

	T_ref :== 300 {K};
	p_ref :== 1 {bar};
	h_ref :== -4.40119 {kJ/kg};
	s_ref :== 0. {kJ/kg/K};

	c_p * M = 1{GAS_C} * SUM[a[i]*T^i | i IN [0..4]];

	(h - h_ref) * M = 1{GAS_C} * SUM[a[i]/(i+1) * T^(i+1) | i IN[0..4]];
	
	s0 IS_A specific_entropy;
	s0 = R * (a[0]*ln(T/T_ref) + SUM[a[i]/i * (T - T_ref)^i | i IN[1..4]]) + 1.294559 {kJ/kg/K} + 0.38191663487 {kJ/kg/K};

	s = s0 - R * ln(p/p_ref);

METHODS
METHOD bound_self;
	RUN ideal_gas_base::bound_self;
	s0.lower_bound := -1e20 {kJ/kg/K};
END bound_self;
METHOD on_load;
	RUN ClearAll;
	RUN bound_self;
	FIX T, p;
	T := 200 {K};
	p := 1 {bar};
END on_load;	
END ideal_air_poly;

IMPORT "johnpye/datareader/datareader";
MODEL drconf;
	filename IS_A symbol_constant;
	format IS_A symbol_constant;
	parameters IS_A symbol_constant;
END drconf;
	
MODEL ideal_air_table REFINES ideal_gas_base;
	M :== 28.958600656 {kg/kmol};

	dataparams IS_A drconf;
	dataparams.filename :== 'johnpye/idealair.csv';
	dataparams.format :== 'CSV';
	dataparams.parameters :== '1,4,6';
	
	s0 IS_A specific_entropy;	
	u IS_A specific_energy;
	p_ref IS_A pressure;

	data: datareader(
		T : INPUT;
		u, s0 :OUTPUT;
		dataparams : DATA
	);

	h = u + R * T; 

	s = s0 - R * ln(p/p_ref);
END ideal_air_table;
	
(*
	Thermo properties
*)
MODEL air_state;
	air IS_A ideal_air_poly;
	p ALIASES air.p;
	T ALIASES air.T;
	h ALIASES air.h;
	s ALIASES air.s;
	v ALIASES air.v;
METHODS
	METHOD default;
		p := 10{bar};
		p.nominal := 42 {bar};
		h := 2000 {kJ/kg};

		T := 400 {K};
		v.nominal := 10 {L/kg};
		s := 4 {kJ/kg/K};
	END default;
	METHOD solve;
		EXTERNAL do_solve(SELF);
	END solve;
	METHOD on_load;
		RUN default_all;
		FIX p, h;
	END on_load;
END air_state;


(* a simple connector that includes calculation of steam properties *)
MODEL air_node;
	state IS_A air_state;
	p ALIASES state.p;
	h ALIASES state.h;
	v ALIASES state.v;
	T ALIASES state.T;
	s ALIASES state.s;
	mdot IS_A mass_rate;
METHODS
	METHOD default;
		mdot.nominal := 2 {kg/s};
	END default;
	METHOD solve;
		EXTERNAL do_solve(SELF);
	END solve;
	METHOD on_load;
		RUN default; RUN reset; RUN values;
		FIX p,h;
	END on_load;
END air_node;

MODEL air_equipment;
	inlet "in: inlet air stream" IS_A air_node;
	outlet "out: outlet air stream" IS_A air_node;

	inlet.mdot, outlet.mdot ARE_THE_SAME;
	mdot ALIASES inlet.mdot;
END air_equipment;


MODEL compressor REFINES air_equipment;
	NOTES
		'block' SELF {Simple model of a compressor using isentropic efficiency}
	END NOTES;

	dp IS_A delta_pressure;
	inlet.p + dp = outlet.p;

	outlet_is IS_A air_state;
	outlet_is.p, outlet.p ARE_THE_SAME;

	outlet_is.s, inlet.s ARE_THE_SAME;
	eta IS_A fraction;

	r IS_A factor;
	r * inlet.p = outlet.p;
	
	eta_eq:eta * (inlet.h - outlet.h) = (inlet.h - outlet_is.h);

	(* work done on the environment, will be negative *)
	Wdot IS_A energy_rate;
	Wdot_eq:Wdot * eta = mdot * (inlet.h - outlet_is.h);

	w IS_A specific_energy;
	w_eq:w = eta * (outlet.h - inlet.h);

END compressor;

MODEL compressor_test REFINES compressor;
	(* no equations here *)
METHODS
METHOD on_load;
	FIX inlet.T;
	FIX inlet.p;

	inlet.T := 300 {K};
	inlet.p := 1 {bar};

	FIX r;
	FIX eta;
	FIX mdot;

	r := 8;
	eta := 0.8;
	mdot := 1 {kg/s};
END on_load;
END compressor_test;


MODEL gas_turbine REFINES air_equipment;
	NOTES
		'block' SELF {Simple turbine model}
	END NOTES;

	dp IS_A delta_pressure;
	inlet.p + dp = outlet.p;
	
	outlet_is IS_A air_state;
	outlet_is.p, outlet.p ARE_THE_SAME;
	outlet_is.s, inlet.s ARE_THE_SAME;

	eta IS_A fraction;
	eta_eq:eta * (inlet.h - outlet_is.h) = (inlet.h - outlet.h);
	
	(* work done on the environment, will be positive *)
	Wdot IS_A energy_rate;
	Wedot_eq:Wdot = mdot * (inlet.h - outlet.h);

	w IS_A specific_energy;
	w_eq:w = inlet.h - outlet.h;

	r IS_A factor;
	r * outlet.p = inlet.p;

END gas_turbine;

MODEL gas_turbine_test REFINES gas_turbine;
	(* no equations here *)
METHODS
METHOD on_load;
	FIX inlet.p;
	FIX inlet.T;
	FIX outlet.p;
	FIX eta;
	FIX mdot;

	inlet.p := 15 {bar};
	inlet.T := 1200 {K};
	outlet.p := 1 {bar};
	eta := 0.85;
	mdot := 1 {kg/s};
END on_load;
END gas_turbine_test;


(*
	simple model assumes no pressure drop, but heating losses due to
	flue gas temperature
*)
MODEL combustor REFINES air_equipment;
	NOTES
		'block' SELF {Simple combustion chamber model}
	END NOTES;

	inlet.p, outlet.p ARE_THE_SAME;
	Qdot_fuel IS_A energy_rate;
	Qdot IS_A energy_rate;

	Qdot = mdot * (outlet.h - inlet.h);

	eta IS_A fraction;
	Qdot = eta * Qdot_fuel;

	qdot_fuel IS_A specific_energy_rate;
	qdot_fuel * mdot = Qdot_fuel;
END combustor;

MODEL combustor_test REFINES combustor;
	(* nothing here *)
METHODS
METHOD on_load;
	FIX inlet.p;
	FIX inlet.T;
	FIX eta;
	FIX outlet.T;
	FIX mdot;

	inlet.p := 15 {bar};
	inlet.T := 500 {K};

	eta := 0.8;
	outlet.T := 700 {K};
	mdot := 1 {kg/s};
END on_load;
END combustor_test;


(*
	this is really simple (fluid props permitting): just work out the heat
	that must be expelled to get the gas down to a specified temperature
*)
MODEL dissipator REFINES air_equipment;
	NOTES
		'block' SELF {Simple condenser model}
	END NOTES;

	inlet.p, outlet.p ARE_THE_SAME;
	Qdot IS_A energy_rate;

	Qdot = mdot * (outlet.h - inlet.h);
	
END dissipator;

MODEL dissipator_test REFINES dissipator;
	(* nothing here *)
METHODS
METHOD on_load;
	FIX inlet.p, inlet.T;
	FIX outlet.T;
	FIX mdot;

	inlet.p := 1 {bar};
	inlet.T := 500 {K};
	outlet.T := 300 {K};
	mdot := 1 {kg/s};
END on_load;
END dissipator_test;


MODEL brayton;
	NOTES
		'description' SELF {
			This is a model of a simple Brayton cycle with
			irreversible compressor (eta=0.8) and turbine (eta=0.85) operating
			between 300 K and 1300 K, with a compression ratio of 8 and an
			assumed inlet pressure of 1 bar. Based on examples 9-5 and 9-6 from
			Çengel & Boles, 'Thermodynamics: An Engineering Approach', 6th Ed,
			McGraw-Hill, 2008}
	END NOTES;

	CO IS_A compressor;	
	TU IS_A gas_turbine;
	BU IS_A combustor;
	DI IS_A dissipator;

	CO.outlet, BU.inlet ARE_THE_SAME;
	BU.outlet, TU.inlet ARE_THE_SAME;
	TU.outlet, DI.inlet ARE_THE_SAME;
	DI.outlet, CO.inlet ARE_THE_SAME;
	braytonpressureratio IS_A positive_factor;
	braytonpressureratio * CO.inlet.p = TU.outlet.p;

	Wdot_CO ALIASES CO.Wdot;
	Wdot_TU ALIASES TU.Wdot;
	Wdot IS_A energy_rate;
	Wdot = Wdot_CO + Wdot_TU;
	
	Qdot_BU ALIASES BU.Qdot;
	Qdot_DI ALIASES DI.Qdot;

	Qdot IS_A energy_rate;
	Qdot = Qdot_BU + Qdot_DI;

	Edot IS_A energy_rate;
	Edot = Wdot - Qdot;

	eta IS_A fraction;
	eta = Wdot / Qdot_BU;

	r_bw IS_A factor;
	r_bw = -Wdot_CO / Wdot_TU;

	state[1..4] IS_A air_node;
	state[1], CO.inlet ARE_THE_SAME;
	state[2], BU.inlet ARE_THE_SAME;
	state[3], TU.inlet ARE_THE_SAME;
	state[4], DI.inlet ARE_THE_SAME;

	eta_TU ALIASES TU.eta;
	eta_CO ALIASES CO.eta;

METHODS
METHOD on_load;
	FIX CO.eta, TU.eta;
	CO.eta := 0.8;
	TU.eta := 0.85;
	FIX CO.inlet.T, TU.inlet.T;
	CO.inlet.T := 300 {K};
	TU.inlet.T := 1300 {K};
	FIX CO.r;
	CO.r := 8;
	FIX CO.inlet.p;
	CO.inlet.p := 1 {bar};
	FIX CO.inlet.mdot;
	CO.inlet.mdot := 1 {kg/s};
	FIX BU.eta;
	BU.eta := 1;
END on_load;
END brayton;


(*
	Regenerator: air-to-air heat exchanger

	Assumption: fluid on both sides have the same c_p.
*)
MODEL regenerator REFINES air_equipment;
	inlet_hot, outlet_hot IS_A air_node;
	
	inlet.p, outlet.p ARE_THE_SAME;
	inlet_hot.p, outlet_hot.p ARE_THE_SAME;

	inlet_hot.mdot, outlet_hot.mdot ARE_THE_SAME;
	mdot_hot ALIASES inlet_hot.mdot;

	(* for perfect eps=1 case: inlet_hot.T, outlet.T ARE_THE_SAME;*)

	epsilon IS_A fraction;

	Qdot IS_A energy_rate;
	mdot_min IS_A mass_rate;
	mdot_min = inlet.mdot + 0.5*(inlet.mdot - inlet_hot.mdot + abs(inlet.mdot - inlet_hot.mdot));

	Qdot = epsilon * mdot_min * (inlet_hot.h - inlet.h);
	outlet.h = inlet.h + Qdot/inlet.mdot;
	outlet_hot.h = inlet_hot.h - Qdot/inlet_hot.mdot;
END regenerator;

MODEL regenerator_test REFINES regenerator;
METHODS
METHOD on_load;
	FIX inlet.mdot, inlet.p, inlet.T;
	FIX inlet_hot.mdot, inlet_hot.p, inlet_hot.T;
	inlet.mdot := 1 {kg/s};
	inlet.p := 1 {bar};
	inlet.T := 300 {K};
	inlet_hot.mdot := 1.05 {kg/s};
	inlet_hot.p := 15 {bar};
	inlet_hot.T := 500 {K};
	FIX epsilon;
	epsilon := 0.8;
END on_load;
END regenerator_test;


MODEL brayton_regenerative;
	NOTES
		'description' SELF {
			This is a model of a regenerative Brayton cycle with
			irreversible compressor (eta=0.8) and turbine (eta=0.85) operating
			between 300 K and 1300 K, with a compression ratio of 8 and an
			assumed inlet pressure of 1 bar. The regenerator effectiveness is 
			0.8.

			Based on example 9-7 from Çengel & Boles, 'Thermodynamics: An 
			Engineering Approach', 6th Ed, McGraw-Hill, 2008}
	END NOTES;

	CO IS_A compressor;	
	TU IS_A gas_turbine;
	BU IS_A combustor;
	DI IS_A dissipator;
	RE IS_A regenerator;

	CO.outlet, RE.inlet ARE_THE_SAME;
	RE.outlet, BU.inlet ARE_THE_SAME;
	BU.outlet, TU.inlet ARE_THE_SAME;
	TU.outlet, RE.inlet_hot ARE_THE_SAME;
	RE.outlet_hot, DI.inlet ARE_THE_SAME;
	DI.outlet, CO.inlet ARE_THE_SAME;

	Wdot_CO ALIASES CO.Wdot;
	Wdot_TU ALIASES TU.Wdot;
	Wdot IS_A energy_rate;
	Wdot = Wdot_CO + Wdot_TU;
	
	Qdot_BU ALIASES BU.Qdot;
	Qdot_DI ALIASES DI.Qdot;

	Qdot IS_A energy_rate;
	Qdot = Qdot_BU + Qdot_DI;

	Edot IS_A energy_rate;
	Edot = Wdot - Qdot;

	eta IS_A factor;
	eta = Wdot / Qdot_BU;

	r_bw IS_A factor;
	r_bw = -Wdot_CO / Wdot_TU;

	Qdot_RE ALIASES RE.Qdot;

	eta_TU ALIASES TU.eta;
	eta_CO ALIASES CO.eta;
	epsilon_RE ALIASES RE.epsilon;

	braytonpressureratio ALIASES CO.r;
METHODS
METHOD on_load;
	FIX CO.eta, TU.eta;
	CO.eta := 0.88;
	TU.eta := 0.85;
	FIX CO.inlet.T, TU.inlet.T;
	CO.inlet.T := 300 {K};
	TU.inlet.T := 1300 {K};
	FIX CO.r;
	CO.r := 4.5;
	FIX CO.inlet.p;
	CO.inlet.p := 1 {bar};
	FIX CO.inlet.mdot;
	CO.inlet.mdot := 1 {kg/s};
	FIX BU.eta;
	BU.eta := 1;
	FIX RE.epsilon;
	RE.epsilon := 0.8;
END on_load;
END brayton_regenerative;


MODEL brayton_intercool_reheat_regen;
	NOTES
		'description' SELF {
			This is a model of a Brayton cycle with intercooling, reheating, and
			regeneration.

			It has an irreversible compressor (eta=0.8) and turbine (eta=0.85)
			and operates between 300 K and 1300 K, with a compression ratio of 8
			and an assumed inlet pressure of 1 bar. The regenerator
			effectiveness is 0.8.

			In adding the intercooling and reheating stages, we assume an
			intermediate pressure that results in two compression stages of
			equal pressure ratio, and two turbine stages of equal pressure 
			ratio.

			Based on example 9-8 from Çengel & Boles, 'Thermodynamics: An 
			Engineering Approach', 6th Ed, McGraw-Hill, 2008}
	END NOTES;

	CO1, CO2 IS_A compressor;	
	TU1, TU2 IS_A gas_turbine;
	BU IS_A combustor;
	DI IS_A dissipator;
	RE IS_A regenerator;
	IC IS_A dissipator;
	RH IS_A combustor;

	CO1.outlet, IC.inlet ARE_THE_SAME;
	IC.outlet, CO2.inlet ARE_THE_SAME;
	CO2.outlet, RE.inlet ARE_THE_SAME;
	RE.outlet, BU.inlet ARE_THE_SAME;
	BU.outlet, TU1.inlet ARE_THE_SAME;
	TU1.outlet, RH.inlet ARE_THE_SAME;
	RH.outlet, TU2.inlet ARE_THE_SAME;
	TU2.outlet, RE.inlet_hot ARE_THE_SAME;
	RE.outlet_hot, DI.inlet ARE_THE_SAME;
	DI.outlet, CO1.inlet ARE_THE_SAME;

	Wdot_CO1 ALIASES CO1.Wdot;
	Wdot_CO2 ALIASES CO2.Wdot;
	Wdot_TU1 ALIASES TU1.Wdot;
	Wdot_TU2 ALIASES TU2.Wdot;

	Wdot_CO, Wdot_TU, Wdot IS_A energy_rate;
	Wdot_CO = Wdot_CO1 + Wdot_CO2;
	Wdot_TU = Wdot_TU1 + Wdot_TU2;
	Wdot = Wdot_CO + Wdot_TU;
	
	Qdot_BU ALIASES BU.Qdot;
	Qdot_DI ALIASES DI.Qdot;
	Qdot_IC ALIASES IC.Qdot;
	Qdot_RH ALIASES RH.Qdot;

	Qdot IS_A energy_rate;
	Qdot = Qdot_BU + Qdot_DI + Qdot_IC + Qdot_RH;

	Edot IS_A energy_rate;
	Edot = Wdot - Qdot;

	eta IS_A factor;
	eta = Wdot / Qdot_BU;

	CO1.r = CO2.r;
	TU1.r = TU2.r;

	RH.outlet.T = BU.outlet.T;
	IC.outlet.T = DI.outlet.T;

	r IS_A factor;
	r = CO2.outlet.p / CO1.inlet.p;

	r_bw IS_A factor;
	r_bw = -Wdot_CO / Wdot_TU;

	Qdot_RE ALIASES RE.Qdot;

	eta_TU1 ALIASES TU1.eta;
	eta_TU2 ALIASES TU2.eta;
	eta_CO1 ALIASES CO1.eta;
	eta_CO2 ALIASES CO2.eta;
	epsilon_RE ALIASES RE.epsilon;
METHODS
METHOD on_load;
	FIX CO1.eta, CO2.eta, TU1.eta, TU2.eta;
	CO1.eta := 0.8;
	CO2.eta := 0.8;
	TU1.eta := 0.85;
	TU2.eta := 0.85;
	FIX CO1.inlet.T, TU1.inlet.T;
	CO1.inlet.T := 300 {K};
	TU1.inlet.T := 1300 {K};
	FIX r;
	r := 8;
	FIX CO1.inlet.p;
	CO1.inlet.p := 1 {bar};
	FIX CO1.inlet.mdot;
	CO1.inlet.mdot := 1 {kg/s};
	FIX BU.eta, RH.eta;
	BU.eta := 1;
	RH.eta := 1;
	FIX RE.epsilon;
	RE.epsilon := 0.8;
END on_load;
END brayton_intercool_reheat_regen;


1	jpye	2687	(* ASCEND modelling environment
2			Copyright (C) 2007, 2008, 2009, 2010 Carnegie Mellon University
3
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			along with this program. If not, see <http://www.gnu.org/licenses/>.
16			)(
17			This file contains some models of Brayton engines and associated cycles,
18			following the development of Çengel & Boles 'Thermodynamcs: An Engineering
19			Approach, 6th Ed, McGraw-Hill, 2008.
20
21			Author: John Pye
22			*)
23
24			REQUIRE "atoms.a4l";
25			REQUIRE "johnpye/thermo_types.a4c";
26			REQUIRE "johnpye/airprops.a4c";
27			IMPORT "sensitivity/solve";
28
29			(* first some models of air as an ideal gas *)
30
31			MODEL ideal_gas_base;
32			M IS_A molar_weight_constant;
33			c_p IS_A specific_heat_capacity;
34			s IS_A specific_entropy;
35			h IS_A specific_enthalpy;
36			v IS_A specific_volume;
37			T IS_A temperature;
38			p IS_A pressure;
39			R IS_A specific_gas_constant;
40
41			R :== 1{GAS_C} / M;
42			p * v = R * T;
43
44			METHODS
45			METHOD bound_self;
46			s.lower_bound := -5 {kJ/kg/K};
47			END bound_self;
48			END ideal_gas_base;
49
50			(*
51			Ideal air assuming ideal gas and constant cp.
52			*)
53			MODEL simple_ideal_air
54			REFINES ideal_gas_base;
55
56			M :== 28.958600656 {kg/kmol};
57
58			c_p = 1.005 {kJ/kg/K};
59
60			T_ref IS_A temperature_constant;
61			p_ref IS_A pressure_constant;
62
63			s = c_p * ln(T / T_ref) - R * ln(p / p_ref);
64			h = c_p * (T - T_ref);
65
66			T_ref :== 273.15 {K};
67			p_ref :== 1 {bar};
68
69			METHODS
70			METHOD on_load;
71			RUN ClearAll;
72			RUN bound_self;
73			FIX T, p;
74			T := 300 {K};
75			p := 1 {bar};
76			END on_load;
77			END simple_ideal_air;
78
79			(*
80			Ideal air, using a quartic polynomial for c_p as given in Moran & Shapiro
81			4th Ed.
82			*)
83			MODEL ideal_air_poly REFINES ideal_gas_base;
84			M :== 28.958600656 {kg/kmol};
85
86			a[0..4] IS_A real_constant;
87			a[0] :== 3.653;
88			a[1] :== -1.337e-3{1/K};
89			a[2] :== 3.294e-6{1/K^2};
90			a[3] :== -1.913e-9{1/K^3};
91			a[4] :== 0.2763e-12{1/K^4};
92
93			T_ref IS_A temperature_constant;
94			p_ref IS_A pressure_constant;
95			h_ref IS_A real_constant;
96			s_ref IS_A real_constant;
97
98			T_ref :== 300 {K};
99			p_ref :== 1 {bar};
100			h_ref :== -4.40119 {kJ/kg};
101			s_ref :== 0. {kJ/kg/K};
102
103			c_p * M = 1{GAS_C} * SUM[a[i]*T^i \| i IN [0..4]];
104
105			(h - h_ref) * M = 1{GAS_C} * SUM[a[i]/(i+1) * T^(i+1) \| i IN[0..4]];
106
107			s0 IS_A specific_entropy;
108			s0 = R * (a[0]ln(T/T_ref) + SUM[a[i]/i (T - T_ref)^i \| i IN[1..4]]) + 1.294559 {kJ/kg/K} + 0.38191663487 {kJ/kg/K};
109
110			s = s0 - R * ln(p/p_ref);
111
112			METHODS
113			METHOD bound_self;
114			RUN ideal_gas_base::bound_self;
115			s0.lower_bound := -1e20 {kJ/kg/K};
116			END bound_self;
117			METHOD on_load;
118			RUN ClearAll;
119			RUN bound_self;
120			FIX T, p;
121			T := 200 {K};
122			p := 1 {bar};
123			END on_load;
124			END ideal_air_poly;
125
126			IMPORT "johnpye/datareader/datareader";
127			MODEL drconf;
128			filename IS_A symbol_constant;
129			format IS_A symbol_constant;
130			parameters IS_A symbol_constant;
131			END drconf;
132
133			MODEL ideal_air_table REFINES ideal_gas_base;
134			M :== 28.958600656 {kg/kmol};
135
136			dataparams IS_A drconf;
137			dataparams.filename :== 'johnpye/idealair.csv';
138			dataparams.format :== 'CSV';
139			dataparams.parameters :== '1,4,6';
140
141			s0 IS_A specific_entropy;
142			u IS_A specific_energy;
143			p_ref IS_A pressure;
144
145			data: datareader(
146			T : INPUT;
147			u, s0 :OUTPUT;
148			dataparams : DATA
149			);
150
151			h = u + R * T;
152
153			s = s0 - R * ln(p/p_ref);
154			END ideal_air_table;
155
156			(*
157			Thermo properties
158			*)
159			MODEL air_state;
160			air IS_A ideal_air_poly;
161			p ALIASES air.p;
162			T ALIASES air.T;
163			h ALIASES air.h;
164			s ALIASES air.s;
165			v ALIASES air.v;
166			METHODS
167			METHOD default;
168			p := 10{bar};
169			p.nominal := 42 {bar};
170			h := 2000 {kJ/kg};
171
172			T := 400 {K};
173			v.nominal := 10 {L/kg};
174			s := 4 {kJ/kg/K};
175			END default;
176			METHOD solve;
177			EXTERNAL do_solve(SELF);
178			END solve;
179			METHOD on_load;
180			RUN default_all;
181			FIX p, h;
182			END on_load;
183			END air_state;
184
185
186			(* a simple connector that includes calculation of steam properties *)
187			MODEL air_node;
188			state IS_A air_state;
189			p ALIASES state.p;
190			h ALIASES state.h;
191			v ALIASES state.v;
192			T ALIASES state.T;
193			s ALIASES state.s;
194			mdot IS_A mass_rate;
195			METHODS
196			METHOD default;
197			mdot.nominal := 2 {kg/s};
198			END default;
199			METHOD solve;
200			EXTERNAL do_solve(SELF);
201			END solve;
202			METHOD on_load;
203			RUN default; RUN reset; RUN values;
204			FIX p,h;
205			END on_load;
206			END air_node;
207
208			MODEL air_equipment;
209			inlet "in: inlet air stream" IS_A air_node;
210			outlet "out: outlet air stream" IS_A air_node;
211
212			inlet.mdot, outlet.mdot ARE_THE_SAME;
213			mdot ALIASES inlet.mdot;
214			END air_equipment;
215
216
217			MODEL compressor REFINES air_equipment;
218			NOTES
219			'block' SELF {Simple model of a compressor using isentropic efficiency}
220			END NOTES;
221
222			dp IS_A delta_pressure;
223			inlet.p + dp = outlet.p;
224
225			outlet_is IS_A air_state;
226			outlet_is.p, outlet.p ARE_THE_SAME;
227
228			outlet_is.s, inlet.s ARE_THE_SAME;
229			eta IS_A fraction;
230
231			r IS_A factor;
232			r * inlet.p = outlet.p;
233
234			eta_eq:eta * (inlet.h - outlet.h) = (inlet.h - outlet_is.h);
235
236			(* work done on the environment, will be negative *)
237			Wdot IS_A energy_rate;
238			Wdot_eq:Wdot * eta = mdot * (inlet.h - outlet_is.h);
239
240			w IS_A specific_energy;
241			w_eq:w = eta * (outlet.h - inlet.h);
242
243			END compressor;
244
245			MODEL compressor_test REFINES compressor;
246			(* no equations here *)
247			METHODS
248			METHOD on_load;
249			FIX inlet.T;
250			FIX inlet.p;
251
252			inlet.T := 300 {K};
253			inlet.p := 1 {bar};
254
255			FIX r;
256			FIX eta;
257			FIX mdot;
258
259			r := 8;
260			eta := 0.8;
261			mdot := 1 {kg/s};
262			END on_load;
263			END compressor_test;
264
265
266
267			MODEL gas_turbine REFINES air_equipment;
268			NOTES
269			'block' SELF {Simple turbine model}
270			END NOTES;
271
272			dp IS_A delta_pressure;
273			inlet.p + dp = outlet.p;
274
275			outlet_is IS_A air_state;
276			outlet_is.p, outlet.p ARE_THE_SAME;
277			outlet_is.s, inlet.s ARE_THE_SAME;
278
279			eta IS_A fraction;
280			eta_eq:eta * (inlet.h - outlet_is.h) = (inlet.h - outlet.h);
281
282			(* work done on the environment, will be positive *)
283			Wdot IS_A energy_rate;
284			Wedot_eq:Wdot = mdot * (inlet.h - outlet.h);
285
286			w IS_A specific_energy;
287			w_eq:w = inlet.h - outlet.h;
288
289			r IS_A factor;
290			r * outlet.p = inlet.p;
291
292			END gas_turbine;
293
294			MODEL gas_turbine_test REFINES gas_turbine;
295			(* no equations here *)
296			METHODS
297			METHOD on_load;
298			FIX inlet.p;
299			FIX inlet.T;
300			FIX outlet.p;
301			FIX eta;
302			FIX mdot;
303
304			inlet.p := 15 {bar};
305			inlet.T := 1200 {K};
306			outlet.p := 1 {bar};
307			eta := 0.85;
308			mdot := 1 {kg/s};
309			END on_load;
310			END gas_turbine_test;
311
312
313
314
315			(*
316			simple model assumes no pressure drop, but heating losses due to
317			flue gas temperature
318			*)
319			MODEL combustor REFINES air_equipment;
320			NOTES
321			'block' SELF {Simple combustion chamber model}
322			END NOTES;
323
324			inlet.p, outlet.p ARE_THE_SAME;
325			Qdot_fuel IS_A energy_rate;
326			Qdot IS_A energy_rate;
327
328			Qdot = mdot * (outlet.h - inlet.h);
329
330			eta IS_A fraction;
331			Qdot = eta * Qdot_fuel;
332
333			qdot_fuel IS_A specific_energy_rate;
334			qdot_fuel * mdot = Qdot_fuel;
335			END combustor;
336
337			MODEL combustor_test REFINES combustor;
338			(* nothing here *)
339			METHODS
340			METHOD on_load;
341			FIX inlet.p;
342			FIX inlet.T;
343			FIX eta;
344			FIX outlet.T;
345			FIX mdot;
346
347			inlet.p := 15 {bar};
348			inlet.T := 500 {K};
349
350			eta := 0.8;
351			outlet.T := 700 {K};
352			mdot := 1 {kg/s};
353			END on_load;
354			END combustor_test;
355
356
357
358			(*
359			this is really simple (fluid props permitting): just work out the heat
360			that must be expelled to get the gas down to a specified temperature
361			*)
362			MODEL dissipator REFINES air_equipment;
363			NOTES
364			'block' SELF {Simple condenser model}
365			END NOTES;
366
367			inlet.p, outlet.p ARE_THE_SAME;
368			Qdot IS_A energy_rate;
369
370			Qdot = mdot * (outlet.h - inlet.h);
371
372			END dissipator;
373
374			MODEL dissipator_test REFINES dissipator;
375			(* nothing here *)
376			METHODS
377			METHOD on_load;
378			FIX inlet.p, inlet.T;
379			FIX outlet.T;
380			FIX mdot;
381
382			inlet.p := 1 {bar};
383			inlet.T := 500 {K};
384			outlet.T := 300 {K};
385			mdot := 1 {kg/s};
386			END on_load;
387			END dissipator_test;
388
389
390			MODEL brayton;
391			NOTES
392			'description' SELF {
393			This is a model of a simple Brayton cycle with
394			irreversible compressor (eta=0.8) and turbine (eta=0.85) operating
395			between 300 K and 1300 K, with a compression ratio of 8 and an
396			assumed inlet pressure of 1 bar. Based on examples 9-5 and 9-6 from
397			Çengel & Boles, 'Thermodynamics: An Engineering Approach', 6th Ed,
398			McGraw-Hill, 2008}
399			END NOTES;
400
401			CO IS_A compressor;
402			TU IS_A gas_turbine;
403			BU IS_A combustor;
404			DI IS_A dissipator;
405
406			CO.outlet, BU.inlet ARE_THE_SAME;
407			BU.outlet, TU.inlet ARE_THE_SAME;
408			TU.outlet, DI.inlet ARE_THE_SAME;
409			DI.outlet, CO.inlet ARE_THE_SAME;
410			braytonpressureratio IS_A positive_factor;
411			braytonpressureratio * CO.inlet.p = TU.outlet.p;
412
413			Wdot_CO ALIASES CO.Wdot;
414			Wdot_TU ALIASES TU.Wdot;
415			Wdot IS_A energy_rate;
416			Wdot = Wdot_CO + Wdot_TU;
417
418			Qdot_BU ALIASES BU.Qdot;
419			Qdot_DI ALIASES DI.Qdot;
420
421			Qdot IS_A energy_rate;
422			Qdot = Qdot_BU + Qdot_DI;
423
424			Edot IS_A energy_rate;
425			Edot = Wdot - Qdot;
426
427			eta IS_A fraction;
428			eta = Wdot / Qdot_BU;
429
430			r_bw IS_A factor;
431			r_bw = -Wdot_CO / Wdot_TU;
432
433			state[1..4] IS_A air_node;
434			state[1], CO.inlet ARE_THE_SAME;
435			state[2], BU.inlet ARE_THE_SAME;
436			state[3], TU.inlet ARE_THE_SAME;
437			state[4], DI.inlet ARE_THE_SAME;
438
439			eta_TU ALIASES TU.eta;
440			eta_CO ALIASES CO.eta;
441
442			METHODS
443			METHOD on_load;
444			FIX CO.eta, TU.eta;
445			CO.eta := 0.8;
446			TU.eta := 0.85;
447			FIX CO.inlet.T, TU.inlet.T;
448			CO.inlet.T := 300 {K};
449			TU.inlet.T := 1300 {K};
450			FIX CO.r;
451			CO.r := 8;
452			FIX CO.inlet.p;
453			CO.inlet.p := 1 {bar};
454			FIX CO.inlet.mdot;
455			CO.inlet.mdot := 1 {kg/s};
456			FIX BU.eta;
457			BU.eta := 1;
458			END on_load;
459			END brayton;
460
461
462			(*
463			Regenerator: air-to-air heat exchanger
464
465			Assumption: fluid on both sides have the same c_p.
466			*)
467			MODEL regenerator REFINES air_equipment;
468			inlet_hot, outlet_hot IS_A air_node;
469
470			inlet.p, outlet.p ARE_THE_SAME;
471			inlet_hot.p, outlet_hot.p ARE_THE_SAME;
472
473			inlet_hot.mdot, outlet_hot.mdot ARE_THE_SAME;
474			mdot_hot ALIASES inlet_hot.mdot;
475
476			(* for perfect eps=1 case: inlet_hot.T, outlet.T ARE_THE_SAME;*)
477
478			epsilon IS_A fraction;
479
480			Qdot IS_A energy_rate;
481			mdot_min IS_A mass_rate;
482			mdot_min = inlet.mdot + 0.5*(inlet.mdot - inlet_hot.mdot + abs(inlet.mdot - inlet_hot.mdot));
483
484			Qdot = epsilon * mdot_min * (inlet_hot.h - inlet.h);
485			outlet.h = inlet.h + Qdot/inlet.mdot;
486			outlet_hot.h = inlet_hot.h - Qdot/inlet_hot.mdot;
487			END regenerator;
488
489			MODEL regenerator_test REFINES regenerator;
490			METHODS
491			METHOD on_load;
492			FIX inlet.mdot, inlet.p, inlet.T;
493			FIX inlet_hot.mdot, inlet_hot.p, inlet_hot.T;
494			inlet.mdot := 1 {kg/s};
495			inlet.p := 1 {bar};
496			inlet.T := 300 {K};
497			inlet_hot.mdot := 1.05 {kg/s};
498			inlet_hot.p := 15 {bar};
499			inlet_hot.T := 500 {K};
500			FIX epsilon;
501			epsilon := 0.8;
502			END on_load;
503			END regenerator_test;
504
505
506
507			MODEL brayton_regenerative;
508			NOTES
509			'description' SELF {
510			This is a model of a regenerative Brayton cycle with
511			irreversible compressor (eta=0.8) and turbine (eta=0.85) operating
512			between 300 K and 1300 K, with a compression ratio of 8 and an
513			assumed inlet pressure of 1 bar. The regenerator effectiveness is
514			0.8.
515
516			Based on example 9-7 from Çengel & Boles, 'Thermodynamics: An
517			Engineering Approach', 6th Ed, McGraw-Hill, 2008}
518			END NOTES;
519
520			CO IS_A compressor;
521			TU IS_A gas_turbine;
522			BU IS_A combustor;
523			DI IS_A dissipator;
524			RE IS_A regenerator;
525
526			CO.outlet, RE.inlet ARE_THE_SAME;
527			RE.outlet, BU.inlet ARE_THE_SAME;
528			BU.outlet, TU.inlet ARE_THE_SAME;
529			TU.outlet, RE.inlet_hot ARE_THE_SAME;
530			RE.outlet_hot, DI.inlet ARE_THE_SAME;
531			DI.outlet, CO.inlet ARE_THE_SAME;
532
533			Wdot_CO ALIASES CO.Wdot;
534			Wdot_TU ALIASES TU.Wdot;
535			Wdot IS_A energy_rate;
536			Wdot = Wdot_CO + Wdot_TU;
537
538			Qdot_BU ALIASES BU.Qdot;
539			Qdot_DI ALIASES DI.Qdot;
540
541			Qdot IS_A energy_rate;
542			Qdot = Qdot_BU + Qdot_DI;
543
544			Edot IS_A energy_rate;
545			Edot = Wdot - Qdot;
546
547			eta IS_A factor;
548			eta = Wdot / Qdot_BU;
549
550			r_bw IS_A factor;
551			r_bw = -Wdot_CO / Wdot_TU;
552
553			Qdot_RE ALIASES RE.Qdot;
554
555			eta_TU ALIASES TU.eta;
556			eta_CO ALIASES CO.eta;
557			epsilon_RE ALIASES RE.epsilon;
558
559			braytonpressureratio ALIASES CO.r;
560			METHODS
561			METHOD on_load;
562			FIX CO.eta, TU.eta;
563			CO.eta := 0.88;
564			TU.eta := 0.85;
565			FIX CO.inlet.T, TU.inlet.T;
566			CO.inlet.T := 300 {K};
567			TU.inlet.T := 1300 {K};
568			FIX CO.r;
569			CO.r := 4.5;
570			FIX CO.inlet.p;
571			CO.inlet.p := 1 {bar};
572			FIX CO.inlet.mdot;
573			CO.inlet.mdot := 1 {kg/s};
574			FIX BU.eta;
575			BU.eta := 1;
576			FIX RE.epsilon;
577			RE.epsilon := 0.8;
578			END on_load;
579			END brayton_regenerative;
580
581
582
583
584			MODEL brayton_intercool_reheat_regen;
585			NOTES
586			'description' SELF {
587			This is a model of a Brayton cycle with intercooling, reheating, and
588			regeneration.
589
590			It has an irreversible compressor (eta=0.8) and turbine (eta=0.85)
591			and operates between 300 K and 1300 K, with a compression ratio of 8
592			and an assumed inlet pressure of 1 bar. The regenerator
593			effectiveness is 0.8.
594
595			In adding the intercooling and reheating stages, we assume an
596			intermediate pressure that results in two compression stages of
597			equal pressure ratio, and two turbine stages of equal pressure
598			ratio.
599
600			Based on example 9-8 from Çengel & Boles, 'Thermodynamics: An
601			Engineering Approach', 6th Ed, McGraw-Hill, 2008}
602			END NOTES;
603
604			CO1, CO2 IS_A compressor;
605			TU1, TU2 IS_A gas_turbine;
606			BU IS_A combustor;
607			DI IS_A dissipator;
608			RE IS_A regenerator;
609			IC IS_A dissipator;
610			RH IS_A combustor;
611
612			CO1.outlet, IC.inlet ARE_THE_SAME;
613			IC.outlet, CO2.inlet ARE_THE_SAME;
614			CO2.outlet, RE.inlet ARE_THE_SAME;
615			RE.outlet, BU.inlet ARE_THE_SAME;
616			BU.outlet, TU1.inlet ARE_THE_SAME;
617			TU1.outlet, RH.inlet ARE_THE_SAME;
618			RH.outlet, TU2.inlet ARE_THE_SAME;
619			TU2.outlet, RE.inlet_hot ARE_THE_SAME;
620			RE.outlet_hot, DI.inlet ARE_THE_SAME;
621			DI.outlet, CO1.inlet ARE_THE_SAME;
622
623			Wdot_CO1 ALIASES CO1.Wdot;
624			Wdot_CO2 ALIASES CO2.Wdot;
625			Wdot_TU1 ALIASES TU1.Wdot;
626			Wdot_TU2 ALIASES TU2.Wdot;
627
628			Wdot_CO, Wdot_TU, Wdot IS_A energy_rate;
629			Wdot_CO = Wdot_CO1 + Wdot_CO2;
630			Wdot_TU = Wdot_TU1 + Wdot_TU2;
631			Wdot = Wdot_CO + Wdot_TU;
632
633			Qdot_BU ALIASES BU.Qdot;
634			Qdot_DI ALIASES DI.Qdot;
635			Qdot_IC ALIASES IC.Qdot;
636			Qdot_RH ALIASES RH.Qdot;
637
638			Qdot IS_A energy_rate;
639			Qdot = Qdot_BU + Qdot_DI + Qdot_IC + Qdot_RH;
640
641			Edot IS_A energy_rate;
642			Edot = Wdot - Qdot;
643
644			eta IS_A factor;
645			eta = Wdot / Qdot_BU;
646
647			CO1.r = CO2.r;
648			TU1.r = TU2.r;
649
650			RH.outlet.T = BU.outlet.T;
651			IC.outlet.T = DI.outlet.T;
652
653			r IS_A factor;
654			r = CO2.outlet.p / CO1.inlet.p;
655
656			r_bw IS_A factor;
657			r_bw = -Wdot_CO / Wdot_TU;
658
659			Qdot_RE ALIASES RE.Qdot;
660
661			eta_TU1 ALIASES TU1.eta;
662			eta_TU2 ALIASES TU2.eta;
663			eta_CO1 ALIASES CO1.eta;
664			eta_CO2 ALIASES CO2.eta;
665			epsilon_RE ALIASES RE.epsilon;
666			METHODS
667			METHOD on_load;
668			FIX CO1.eta, CO2.eta, TU1.eta, TU2.eta;
669			CO1.eta := 0.8;
670			CO2.eta := 0.8;
671			TU1.eta := 0.85;
672			TU2.eta := 0.85;
673			FIX CO1.inlet.T, TU1.inlet.T;
674			CO1.inlet.T := 300 {K};
675			TU1.inlet.T := 1300 {K};
676			FIX r;
677			r := 8;
678			FIX CO1.inlet.p;
679			CO1.inlet.p := 1 {bar};
680			FIX CO1.inlet.mdot;
681			CO1.inlet.mdot := 1 {kg/s};
682			FIX BU.eta, RH.eta;
683			BU.eta := 1;
684			RH.eta := 1;
685			FIX RE.epsilon;
686			RE.epsilon := 0.8;
687			END on_load;
688			END brayton_intercool_reheat_regen;
689
690
691
692