% *********************************************************************** % Two_Stage_Engine_Simulation_V1.m % % This model predicts the increase in efficiency achieved by using heat % recovery to run a second sub-atmospheric cylinder in the condensing % engine. % % Model uses XSteam database. Holmgren, M., 2007 Steam:Thermodynamic Properties of Water and Steam % [Online] Available at: https://uk.mathworks.com/matlabcentral/fileexchange/9817-x-steam % thermodynamic-properties-of-water-and-steam [Accessed: 21 July 2021]. %************************************************************************ clear all clc % --------------------------- Inputs ------------------------------------- % Inputs for the first stage have a notation '1', and inputs for the second % stage a notation '2'. m1 = 1; % (kg) Mass basis of initial stream to perform the analysis on x1i = 0.9; % Initial steam dryness fraction x2i = 1; % Initial steam dryness fraction T1b = 100; % (bar) First stage boiler pressure T1he = 50; % (C) Define first stage heat exchanger operating temperature T1cond = 20; % (Celsius) Temperature of the first stage condenser T2cond = 20; % (Celsius) Temperature of the second stage condenser he_r = 0.7; % Fraction of available energy that is recovered in HE ER1 = 5; % Expansion Ratio used in the first stage Ni = 100; % Number of iterations to be used in work calculation gamma = 1.08; % Value of adiabatic coefficient of expansion % The heat exchanger steam outlet dryness values to be called if a relevant % heat exchanger inlet temperature is used. % The vectors give the dryness values for a range of ERs from 1 to 10. % The temperature of the heat exchanger limits the range of first stage ERs % that are able to have energy recovered in the heat exchanger (NaN). % These values are read manually from a T-s chart with the condensation % process plotted. New values can be added for any temperature case % required. If this is done, an associated 'elseif' must be added to the % available energy calculation code below to utilise this data. % If 'Use_graphical' is set to 1 then the graphical values input below will % be called. If this variable is set to 0 then the density method will be % used instead. Set this to 0 when the density method is required. This could % be for general testing or if there are no graphical values input % for the chosen heat exchanger operating temperature. Use_Graphical = 1; HEDryness50 = [0.26,0.32,0.40,0.45,0.65,0.76,NaN,NaN,NaN,NaN]; HEDryness60 = [0.40,0.50,0.64,0.70,NaN,NaN,NaN,NaN,NaN,NaN]; HEDryness70 = [0.52,0.70,0.85,NaN,NaN,NaN,NaN,NaN,NaN,NaN]; % -------------- First Stage Cylinder Calculations ------------------------ % Set the boiler pressure based on heat source temperature (bar) P1b = XSteam('psat_T', T1b); % (m3/kg) Find first stage initial vapour and liquid specific volumes v_V1i = XSteam('vV_P',P1b); v_L1i = XSteam('vL_P',P1b); % Calculate the volume of space the inital vapour will occupy in the first cylinder V1i = ((x1i * v_V1i) + ((1- x1i)*v_L1i)) * m1; % (m3) % Calculate boundary work from initially flooding the cylinder W1_b = (P1b * 10^2) * V1i; % (kJ) % Calculate the volume at end of expansion and the volume step to iterate % with V1e = V1i * ER1; % Volume at end of expansion (m3) V1s = (V1e - V1i) / Ni; % Volume step to use in for loop (m3) % Initialise enthalpy value of initial steam (kJ/kg) h1c = (x1i * XSteam('hV_P',P1b)) + ((1-x1i) * XSteam('hL_P',P1b)); % Initialise entropy value of initial steam (kJ/kg) s1c = (x1i * XSteam('sV_P',P1b)) + ((1-x1i) * XSteam('sL_P',P1b)); % Initialise first stage expansion work value (kJ) W1_e = 0; % Calculate expansion work in the first cylinder if ER1 == 1 % i.e. if there is no steam expansion P1e = P1b; % Oulet pressure is equal to boiler pressure (bar) % W1_e remains at zero (kJ) else % i.e. if there is steam expansion % Iterate between initial volume and final volume using calculated volume step for V = V1i : V1s : V1e dW1_e = W1_e / m1; % Scale work done to that of 1kg h = h1c - dW1_e; % (kJ) New enthalpy value based on work energy removed from steam s = s1c; % (kJ) Assume isentropic process P_bar = XSteam('p_hs',h,s); % Calculate new pressure for this iteration P_Pa = P_bar * 1E5; % Convert pressure value to Pascal for calculation % Calculate expansion work (kJ) W1_e = W1_e + ((P_Pa * V1s)/1000); % Save the final pressure value for the iteration (bar) P1e = P_bar; end end % Fing final vapour fraction after expansion at final pressure and entropy x1e = XSteam('x_ps',P1e,s1c); % Calculate outlet enthalpy (kJ/kg) h1e = (x1e * XSteam('hV_P',P1e)) + ((1-x1e) * XSteam('hL_P',P1e)); % Find the operating pressure of the first stage condenser (bar) P1cond = XSteam('psat_T', T1cond); % Calculate back pressure work during stroke (kJ) W1_l = (P1cond * 10^2) * V1e; % Calculate work done in first stage as sum of boundary and expansion work % values minus losses due to back pressure (kJ) W1 = W1_b + W1_e - W1_l; % ---------------- First Stage Heat Exchanger Calculations ---------------- % Set critical value of pressure, below which energy cannot be recovered % (bar) P1crit = XSteam('psat_T',T1he); % Calculate the value of dryness of the outlet vapour from the heat % exchanger and the corresponding energy that can be recovered and sent to % the second stage boiler (E1he). Three possible calculations exist: % (1) If the vapour outlet pressure after expansion is less than the % critical pressure (i.e. the saturation pressure at heat exchanger % operating temerature) then no condensation will occur and E1he is zero. % (2) In the case that the the variable 'Use_Graphical' has been input % as zero then the density ratio method is used to determine dryness values. % (3) If 'Use_Graphical' has been input as 1 and the operating temperature % of the heat exchanger matches one of the temperatures which has pre % defined dryness values calculated by hand graphically and entered as % a vector input then these are called for the relevant ER1 and used to % calculate E1he. if round(P1e,3) <= round(P1crit,3) E1he = 0; % If exit pressure after expansion is less than the critical then % recovered energy from the heat exchanger is zero elseif Use_Graphical == 0 % Find the density of vapour leaving the cylinder after expansion % (kg/m3) rho_1e = XSteam('rhoV_p', P1e); % Find the density of the vapour leaving the first heat exchanger % (kg/m3) rho_1he = XSteam('rhoV_T', T1he); % Calculate the vapour fraction ratio for the heat exchanger outlet x1he = rho_1he / rho_1e; % Check that dryness fraction does not increase through the heat % exchanger if round(x1he,3) > round(x1e,3) x1he = x1e; E1he = 0; else % Calculate the enthalpy value of the wet vapour leaving the heat % exchanger (kJ/kg) h1he = (x1he * XSteam('hV_T',T1he)) + ((1-x1he) * XSteam('hL_T',T1he)); % Calculate available energy for recovery in the heat exchanger (kJ) E1he_a = (h1e - h1he) * m1; % Calculate energy actually recovered based on assumed efficiency (kJ) E1he = he_r * E1he_a; end elseif T1he == 50 % Call the relevant vapour fraction ratio for the heat exchanger outlet x1he = HEDryness50(ER1); % Calculate recovered energy as done above (kJ) h1he = (x1he * XSteam('hV_T',T1he)) + ((1-x1he) * XSteam('hL_T',T1he)); E1he_a = (h1e - h1he) * m1; E1he = he_r * E1he_a; elseif T1he == 60 % Repeat method commented on for 50C case above x1he = HEDryness60(ER1); h1he = (x1he * XSteam('hV_T',T1he)) + ((1-x1he) * XSteam('hL_T',T1he)); E1he_a = (h1e - h1he) * m1; E1he = he_r * E1he_a; elseif T1he == 70 % Repeat method commented on for 50C case above x1he = HEDryness70(ER1); h1he = (x1he * XSteam('hV_T',T1he)) + ((1-x1he) * XSteam('hL_T',T1he)); E1he_a = (h1e - h1he) * m1; E1he = he_r * E1he_a; end % Calculate the enthalpy of the wet vapour sent to the first stage % condenser (kJ/kg) h1cond = h1e - (E1he / m1); % --------------- First Stage Condenser Calculations ---------------------- % Find the enthalpy of the final condensate at condenser operating % temperature (kJ/kg) h1L = XSteam('hL_T', T1cond); % Find the energy that must be removed by the condenser to fully condense % the remaining vapour (kJ) E1cond = (h1cond - h1L) * m1; % Calculate the volume of condensate sent to the air pump (m3) V1cond = XSteam('vL_T', T1cond) * m1; % ----------------- Second Stage Boiler Calculations ---------------------- % Set the pressure of the second stage boiler (bar) P2b = XSteam('psat_T', T1he); % Set the pressure of the second stage condenser (bar) P2cond = XSteam('psat_T', T2cond); % Find the enthalpy of vapour leaving the boiler (kJ/kg) hvap_2 = (x2i * XSteam('hV_p',P2b)) + ((1 - x2i) * XSteam('hL_p', P2b)); % Find the enthalpy of the final condensate at condenser operating % temperature (kJ/kg) h2L = XSteam('hL_T', T2cond); % Energy input required per kg is difference between that of the condensate % returned to the boiler and the final vapour (kJ/kg) Q2in = hvap_2 - h2L; % Calculate the mass of steam that will be produced by the recovered energy % (kg) m2 = E1he / Q2in; % ------------------- Second Stage Cylinder Calculations ------------------ if m2 == 0 % i.e. no heat has been transferred from the 1st to 2nd stage W2 = 0; % then no work will be produced in the 2nd stage else % If heat is transferred then calculate what work is produced in the 2nd stage % Calculate the maximum ER usable in the second cylinder such that cylinder % final pressure does not fall below that of the condenser. for ER2i = 1:10 % Calculate the final pressure at given expansion ratio using ideal gas % law (bar) P2ei = P2b * ((1 / ER2i)^gamma); if P2ei >= P2cond ER2 = ER2i; elseif P2ei < P2cond % Once this criteria is satisfied further expasnion is not possible ER2 = ER2i - 1; % Set the chosen 2nd stage ER as the previous value break; % Breaks the for loop now the maximum ER is known end end % Find second stage initial vapour specific volume (m3/kg) v_V2i = XSteam('vV_P',P2b); % Find second stage initial liquid specific volume (m3/kg) v_L2i = XSteam('vL_P',P2b); % Calculate the volume of space the inital vapour will occupy in the second % cylinder (m3) V2i = ((x2i * v_V2i) + ((1 - x2i) * v_L2i)) * m2; % Calculate boundary work from initially flooding the cylinder (kJ) W2_b = (P2b * 10^2) * V2i; % (kJ) % Calculate the volume at the end of expansion and volume step to iterate V2e = V2i * ER2; % Volume at end of expansion (m3) V2s = (V2e - V2i) / Ni; % Volume step to use in for loop (m3) % Initialise enthalpy value of initial steam (kJ/kg.K) h2c = (x2i * XSteam('hV_P',P2b)) + ((1-x2i) * XSteam('hL_P',P2b)); % Initialise entropy value of initial steam (kJ/kg) s2c = (x2i * XSteam('sV_P',P2b)) + ((1-x2i) * XSteam('sL_P',P2b)); % Initialise second stage expansion work value (kJ) W2_e = 0; if ER2 == 1 % i.e. there is no steam expansion P2e = P2b; else % i.e. there is steam expansion % Iterate between initial volume and final volume using calculated volume step for V = V2i : V2s : V2e % This uses the same emthod as the first stage calculation dW2_2 = W2_e / m2; h = h2c - dW2_2; s = s2c; P_bar = XSteam('p_hs',h,s); P_Pa = P_bar * 1E5; W2_e = W2_e + ((P_Pa * V2s)/1000); % (kJ) P2e = P_bar; end end % Find final vapour fraction after expansion at final pressure and entropy x2e = XSteam('x_ps',P2e,s2c); % Calculate outlet enthalpy (kJ/kg) h2e = (x2e * XSteam('hV_P',P2e)) + ((1-x2e) * XSteam('hL_P',P2e)); % Calculate the back pressure work losses (kJ) W2_l = (P2cond * 10^2) * V2e; % (kJ) % Calculate work done in second stage as sum of boundary and expansion work (kJ) W2 = W2_b + W2_e - W2_l; % Include check which avoids sending negative W2 values to the output % calculations if W2 < 0 W2 = 0; end % ------------------ Second Stage Condenser Calculations ------------------ % Find the energy that must be removed by the condenser to fully condense % the remaining vapour (kJ) E2cond = (h2e - h2L) * m2; % Calculate the volume of condensate sent to the air pump (m3) V2cond = XSteam('vL_T', T2cond) * m2; end % ------------------ Calculate Thermal Energy Required -------------------- % There is only one thermal energy input required, in the first stage % boiler % Found as difference between enthalpies of the vapour leaving the boiler and % the condensate water that returns to the boiler Q1in = (h1c - h1L) * m1; % (kJ) % -------- Calculate Percentage that Energy Recovered is of -------------- % -------------- Total Released During Condensation ----------------------- % Calculate what the total energy released during condensation is (kJ) E1_total = ((h1e - h1L) * m1); % % Calculate percentage that actual recovered energy is of this (%) recover_of_total = ((E1he / E1_total)*100); % ----------------------- Output Calculations ----------------------------- OneStageTE = (W1 / Q1in) * 100; TwoStageTE = ((W1 + W2) / Q1in) * 100; TEIncrease = ((TwoStageTE - OneStageTE) / OneStageTE) * 100; SecondStageTE = (W2 / E1he) * 100; disp(['The one stage thermal efficiency is ', num2str(OneStageTE), '%']) disp(['The two stage thermal efficiency is ', num2str(TwoStageTE), '%']) disp(['The increase in thermal efficiency is ', num2str(TEIncrease), '%']) disp(['The second stage thermal efficiency is ', num2str(SecondStageTE), '%']) disp(['The energy recovered is ', num2str(E1he), 'kJ']) disp(['The energy recovered is ', num2str(recover_of_total), '% of the total released during condensation'])