% Thermodynamic Engine Modelling
% Thomas Holland
% University of Wales TSD

%Constants

Bore = 0.079;                    %Cylinder Bore (m)
Stroke = 0.0814;                 %Piston Stroke (m)
Rod = 0.137;                     %Con Rod Length (m)
Clearance = 39.27;               %Clearance Volume (cc)
LHV = 44e6;                      %Fuel Calorific Value (J/Kg)
a = 5;                    
m = 2;
Ignition = -11;            %Spark Timing (degrees)
IgnDel = 5;                %Ignition Delay (degrees)
gamma = 1.25;              %Ratio of Specific Heats
bdur = 40;                 %Burn Duration (degrees)
R = 287;                    %Specific Gas Constant
rho = 1.167;               %Air Density (kg/m^3)
Combeff = 0.9;             %Combustion Efficency
Voleff = 0.9;              %Volumetric Efficency
AFR = 14.7;                 %Air/Fuel Ratio

%Equations
Vd = (pi/4) * Bore^2 * Stroke;
Air_Mass = rho * Vd * Voleff;
Fuel_Mass = Air_Mass / AFR;
Total_Mass = Air_Mass + Fuel_Mass;
Pcom = 100000;
Vc = Clearance/1000000;
Radius = Stroke/2;
CR = (Vc + Vd) / Vc;
Crank = -180:1:180;
s = Radius*cosd(Crank) + sqrt(Rod^2 - (Radius*sind(Crank)).^2);
Vol = Vc + (pi/4)*Bore^2 * (Radius + Rod - s);

% Calculation Loop for Non-Combustion Cylinder Pressure (n compression = 1.33) & (n expansion = 1.27)
P = zeros(1,length(Crank));
P(1) = Pcom;
for i = 1:length(Crank)-1

    if Crank(i) < 0
        n = 1.311;
    else
        n = 1.365;
    end
    P(i+1) = P(i) * (Vol(i) / Vol(i+1))^n;
end

% Definition of Wiebe Function
xb = zeros(1,length(Crank));
SoB = Ignition + IgnDel;
EoB = SoB + bdur;

for i=1:length(Crank)
    if Crank(i) < SoB
        xb(i) = 0;
    elseif Crank(i) <= EoB
        xb(i) = 1 - exp(-a * ((Crank(i)-SoB)/bdur)^(m+1));
    else
        xb(i) = 1;
    end
end


% Calculation for dQr (Heat Released per Degree)
Qtol = Combeff * Fuel_Mass * LHV;
dQr= zeros(1,length(Crank));
dxb = zeros(1,length(Crank));
Pcomp = P;
for i = 1:length(Crank)
    if Crank(i) >= SoB && Crank(i) <= EoB
        theta_term = (Crank(i) - SoB) /bdur;
        dxb(i) = (a*(m+1)/bdur)*(theta_term^m) * exp(-a*(theta_term^(m+1)));
        dQr(i) = Qtol * dxb(i);
    else
        dxb(i) = 0;
        dqr(i) = 0;
    end
end

% Calculation Loop for Combustion Cylinder Pressure
Pcomb = zeros(1,length(Crank));
Pcomb(1) = P(1);
for i = 1:length(Crank)-1
    if Crank(i) < SoB
        Pcomb(i+1) = P(i+1);
    elseif Crank(i) <= EoB
        dV = Vol(i+1) - Vol(i);
        dQ = dQr(i);
        dP = ((gamma - 1)/Vol(i)) * dQ - (gamma * Pcomb(i)/Vol(i)) * dV;
        Pcomb(i+1) = Pcomb(i) + dP;
    else
        Pcomb(i+1) = Pcomb(i) * (Vol(i)/Vol(i+1))^1.27;
    end    
end



%Calculations for Temperature
Tempnon = (P .* Vol) / (Total_Mass * R);
Tempcomb = (Pcomb .* Vol) / (Total_Mass * R);

%Test in Grader
HeatReleased=Qtol;
PeakP=max(Pcomb)/100000;
PeakT=max(Tempcomb);

%Upon initial comparison between the real data analysis and the simulation
%model above there is one big initial weakness of the simulation. All
%components of the Wiebe function are limited with the simulation
%requiring constants that must be input before running which restricts the
%real world application. This model also assumes a fully ideal scenario
%where there is no loss of energy through either heat transfer or airflow
%turbulence as it assumes a constant laminar flow of air into and out of
%the combustion chamber.

%Improvements upon this model can potentially come from a more active Wiebe
%function implementation where