/* MSW Power Plant
******************************************

Hypotheses:

1. Steady state
2. Negligible heat rejection, negligible kinetic and potential energy variations
3. Negligible pressure drop on pipes
4. Combustion products can be modeled as air being ideal gas

Analyze:

*/

// Incinerator

// Mass balance

rac = 17.24

// mdot1 = 69.8  // kg/s

// alfa = 0.5  // percentual of CH4 in fuel

mdot1 = mdotch + mdotar

mdotar = rac * mdotch

mdotch4 = alfa * mdotch

mdotrsu = (1 - alfa) * mdotch

// MSW lower heating value

LHVrsu = 1300 * 4.189  // kJ/kg

// CH4 lower heating value

LHVch4 = 11940 * 4.189 // kJ/kg

qcomb = mdotrsu * LHVrsu + mdotch4 * LHVch4

// Energy balance

// q0 = 0

q0 + qcomb = mdot1 * hprod

hprod = h_T("Air",T1)

// T1 = T_h("Air", hprod)

// Environmental conditions

T0 = 298  // K

// Incinerator sizing

hext = 1  // external natural convection coefficient

hint = 500 // internal forced convection coefficient + radiation

tp = 0.1  // m - insulation thickness

kiso = 1  // W/mK - insulation thermal conductivity

D = 1  // m - incinerator diameter

// pi = 3.14159

hinc = 2 // m - incinerator height

uinv = 1/hext + tp/kiso + 1/hint

u = 1/uinv  // W/m2K - global heat transfer coefficient

area = pi * D * D / 2 + pi * D * hinc 

q0 = - u * area * (T1 - T0) / 1000 // heat leak to ambient

// Recovery heat exchanger sizing - Boiler

cpr = 1 // kJ/kgK

cag = 4.18  // kJ/kgK

Cprod = mdot1 * cpr

Cmin = Cprod

// mdot1 * (h1 - h2) = effet * Cmin * (T1 - T3)

T1 - T2 = effet * (T1 - T3)

effet = 1 - exp (-ua_gv/Cmin)


// Power generation system

// Conversion of units from English system to SI

	vprodE=2E5 // ft3/min

	 vprod = vprodE * 0.028317 / 60     // ft3/min to m3/s

	 T1E=400  // oF

	T1 = 5 * (T1E - 32)/9 + 273     // K

	p1 = 1.01325       // atm or bar

	T2E = 260  //oF

	T2 = 5 * (T2E - 32)/9 + 273     // K

	// T2 = T3

	p2 = 1.01325     // bar

	mdot3 = 275 * 0.4536 / 60   // kg/s

	p3 = 1.01325 * 40 / 14.7   // bar 

	T3 = 5 * (102 - 32)/9 + 273     // K

	p5E = 1 // psi

	p5 = 1.01325* p5E / 14.7   // bar

	x5 = 0.93

	celet = 0.08  // $/kWh

	rbar = 8.314  // kJ/(kmol.K)

	m = 28.97  // kg/kmol

	p4 = p3

	// Mass balance

	0 = mdot1 - mdot2  // HX products side

	0 = mdot3 - mdot4  // HX water side

	0 = mdot4 - mdot5  // vapor in turbine	

	// Energy balance

	wdot = mdot1 * (h1 - h2) + mdot3 * (h3 - h5)

	// Products mass flow rate

	mdot1 = p1 * 1E2 * vprod / (rbar/m) / T1  
	
	// Enthalpies calculation

	h1 = h_T("Air", T1) 

	h2 = h_T("Air", T2)

	h3 = h_PT("Water/Steam", p3, T3)

	h5 = hsat_Px("Water/Steam", p5, x5) 

	// T4 determination

	0 = mdot1 * (h1 - h2) + mdot3 * (h3 - h4)

	T4 = T_Ph("Water/Steam", p4, h4)

	// Isoentropic efficiency

	s5s = s4

	x5s0 = x_sP("Water/Steam", s5s, p0) 

	h5s0 = hsat_Px("Water/Steam", p0, x5s0)  

	x5s = x_sP("Water/Steam", s5s, p5) 

	h5s = hsat_Px("Water/Steam", p5, x5s)  

	eta_s_t = (h4 - h5) / (h4 - h5s)  // turbine isoentropic efficiency
	
	wextra = eta_s_t * mdot4 * (h5s0 - h5s) // Extra power necessary for vaccum

	// Condenser

	// Energy balance

	DT = 0 // Temperature difference below the working fluid condenser input temperature

	T7 = T5 - DT  // Desired condenser water output temperature

	h7 = h_PT("Water/Steam", p0, T7)  

	0 = mdot5 * (h5 - h3) + mdot6 * (h0ag - h7)

// Condenser sizing

Cmin_c = mdot6 * cag

T5 = Tsat_P("Water/Steam", p5)  

T7 - T0 = effet_c * (T5 - T0)

effet_c = 1 - exp (-ua_c/Cmin_c)

FF = ua_c + mdot6

	// Annual value 

	Vanual = wdot * 8000 * celet  // $/ano

	// Exergetic analysis

	// T0 = 298  // K

	p0 = 1.01325  // bar

	h0ar = h_T("Air", T0)  // kJ/kg

	s0ar = s_TP("Air", T0, p0)  // kJ/(kg.K)

	h0ag = h_PT("Water/Steam", p0, T0)

	s0ag = s_PT("Water/Steam", p0, T0)

	// Specific flow exergies calculation

	s1 = s_TP("Air", T1, p1)

	ex1 = (h1 - h0ar) - T0 * (s1 - s0ar)

	s2 = s_TP("Air", T2, p2)

	ex2 = (h2 - h0ar) - T0 * (s2 - s0ar)

	s3 = s_PT("Water/Steam", p3, T3)

	ex3 = (h3 - h0ag) - T0 * (s3 - s0ag)

	s4 = s_PT("Water/Steam", p4, T4)

	ex4 = (h4 - h0ag) - T0 * (s4 - s0ag)

	s5 = ssat_Px("Water/Steam", p5, x5)

	ex5 = (h5 - h0ag) - T0 * (s5 - s0ag)

	s7 = s_PT("Water/Steam", p0, T7)

	ex7 = (h7 - h0ag) - T0 * (s7 - s0ag)

	// Calculate system destroyed exergy and 2nd law efficiency

	wdot = mdot1 * (ex1 - ex2) + mdot3 * (ex3 - ex5) - EDsist_orig  // ED in original system

	wdot = mdot1 * (ex1 - ex2) + mdot6 * (- ex7) - EDsist // ED in complete system

	etaIIsist = (wdot - wbomba) / (mdot1*ex1) // later deduct necessary vaccum extra power

	// Power consumed by pump

	v3 = v_PT("Water/Steam", p3, T3)

	wbomba = mdot4 * v3 * (p3 - p5) * 1E2 // kW

	// Is the temperature rise at the pump output significant? Answer: energy balance - see it is negligible

	wbomba = mdot3 * (h3-h8)

	T8 = T_Ph("Water/Steam", p5, h8) 

	// Thermoeconomic analysis

	// Steam generator - boiler

	Zgv = 0.01  // $/s = CI and OM cost rate

	// Products cost

	c1 = 51E-7  // $/kJ

	// Feed water cost

	rwtwb = wbomba/wdot

	c3 = cw 

	c2 = c1

	// Thermoeconomics accounting equation

	c2 * mdot2 * ex2 + c4 * mdot4 * ex4 = c1 * mdot1 * ex1 + c3 * mdot3 * ex3 + Zgv

	C4 = c4 * mdot4 * ex4 // $/s

	C4_1ano = C4 * 365 * 24 * 3600 // $/year

	// turbine

	Zt = 0.08 * wdot * 1E-4  // $/s

	c5 = c4

	c5 * mdot5 * ex5 + cw * wdot = C4 + Zt

	cwconv = cw * 3600  // $/kWh