% ***********************************************************************
% Uniflow_Engine_Isentropic.m
%
% This model predicts the efficiency of the uniflow engine as described
% in associated thesis. This model assumes isentropic expansion
% using ideal gas laws.
%************************************************************************
 
 
clear all
clc
 
% ---------------------------------Inputs----------------------------------
 
T_cond = 20; % Condenser operating temperature (C)
Pb = 1; % Boiler pressure (bar)
cyl_bore = 50E-3; % Cylinder bore diameter (m)
cyl_swept = 100E-3; % Stroke lenth (m)
ER = 1; % Expansion ratio used in the cylinder
xi = 1; % Vapour fraction of steam supplied by the boiler
gamma = 1.08; % Adiabatic coefficient of expansion
 
 
% --------------------- Downstroke (Power Loss Stroke)-----------------
 
% Calculate the total cylinder volume to be swept (m3)
V_swept = (((pi * (cyl_bore^2)) / 4) * cyl_swept);
 
% Calculate the volume of the cylinder that will be swept with atmospheric
% steam producing boundary work (m3)
V_boundary = V_swept * (1/ER);
 
% Set the volume of steam required equal to that of the boundary volume (m3)
steam_vol_required = V_boundary;
 
% Calculate mass of steam required (kg)
steam_mass_required = XSteam('rhoV_p', Pb) * steam_vol_required;
 
% Calculate boundary work produced from sweeping this cylinder volume
% This does not include work from the excess steam
W_down_b = ((Pb * 1E5) - 101325) * V_boundary; % (J)
 
% Calculate the expansion work produced (J)
if ER > 1
% Set the inital volume for the expansion as that of the boundary volume
V_exp_0 = V_boundary;
 
% Set the final volume after expansion (m3) 
V_exp_1 = V_swept;
 
% Set the starting pressure for the expansion (Pa)
P_exp_0 = Pb * 1E5;
 
% Calculate the final pressure after expansion (Pa)
P_exp_1 = P_exp_0 * ((V_exp_0 / V_exp_1)^gamma);
 
% Save this pressure value for later (bar)
P_cyl_exp = P_exp_1 / (1E5);
 
% Calculate work produced during the expansion (J)
W_exp_e = (1/(gamma - 1)) * ((P_exp_0 * V_exp_0) - (P_exp_1 * V_exp_1));
 
% Calculate losses from atmospheric pressure (J)
W_exp_atm = 101325 * (V_exp_1 - V_exp_0);
 
% Resolve expansion work (J)
W_down_e = W_exp_e - W_exp_atm;      
    
else % i.e. there is no expansion
    P_cyl_exp = Pb; % Pressure at BDC is still equal to boiler pressure (bar)
    W_down_e = 0; % There is no expansion work
end
 
% Find entropy of steam, assumed constant over the expansion (kJ/kg) 
s_in_cyl = (xi * XSteam('sV_P',Pb)) + ((1-xi) * XSteam('sL_P',Pb));
 
% Find the steam outlet vapour fraction after expansion
x_exp = XSteam('x_ps',P_cyl_exp,s_in_cyl);
 
% Calculate steam outlet enthalpy after expansion (kJ/kg)
h_exp = (x_exp * XSteam('hV_P',P_cyl_exp)) + ((1-x_exp) * XSteam('hL_P',P_cyl_exp)); 
 
% Find saturation temperature of outlet steam (deg C)
T_exp_sat = XSteam('Tsat_p', P_cyl_exp);
 
% Resolve boundary, expansion, and back pressure work terms (kJ)
W_down = W_down_b + W_down_e;
 
 
 
% ---------------------- Steam Evacuation --------------------------------
 
% It is assumed that all steam is evacuated from the cylinder in this case.
% The slot has been shown through calculation large enough to allow for
% full evacuation. Any residual steam that would be left is allowed to the
% condenser through the open outlet valve during the upstroke.
% As an approximation the pressure in the cylinder during the upstroke is
% maintained at equal that of the condenser.
 
P_cond = XSteam('psat_T', T_cond);
P_cyl_evac = P_cond;
 
%-------------------Work From Upstroke (Power Stroke)---------------------
 
% In the upstroke the atmospheric pressure external to the piston drives it
% upwards against the negative pressures inside the cylinder. During this
% stroke the residual steam inside the cylinder will be compressed. The
% pressure rise inside the cylinder will reduce the power generated by the
% piston. This is undesirable. In this case an outlet valve maintains the
% pressure inside the cylinder equal to that of the condenser
 
W_up = (101325 - (P_cyl_evac * 1E5)) * V_swept; %(J)
 
 
% ----------------------- Energy Requirements -----------------------------
% The mass of steam entered per stroke is what is required to fill the
% boundary volume.
 
% The required energy input is considered to be that to raise steam
% from the condensate pumped back to the boiler (J)
h_condensate = XSteam('hL_T',T_cond);
h_vapour = XSteam('hV_p',Pb);
Q_in = ((h_vapour - h_condensate) * steam_mass_required) * 1000;
 
 
%--------------------------------Model Outputs-----------------------------
 
% Calculate the total work produced during the whole stroke
W_stroke = W_down + W_up;
   
% Calculate the thermal efficiency of the engine
Thermal_Efficiency = (W_stroke / Q_in) * 100;
disp(['The thermal efficiency of the uniflow engine operating with an ER of ', num2str(ER), ' is ', num2str(Thermal_Efficiency), '%'])