% Thermodynamic Engine Modelling
% Dr Owen Williams
% University of Wales TSD

Bore=98;                %Cylinder Bore (mm)
Stroke=66;              %Piston Stroke (mm)
Rod=120;                %Con Rod Length (mm)
Clearance=58.8;         %Clearance Volume (cc)
LHV=44e6;               %Fuel Calorific Value (J/Kg)
a=9;                    
m=2.1;
Ignition=-5;           %Spark Timing (degrees)
IgnDel=0;               %Ignition Delay (degrees)
gamma=1.25              %Ratio of Specific Heats
bdur=35                 %Burn Duration (degrees)
rho=1.167               %Air Density (kg/m^3)
Combeff=0.9             %Combustion Efficency
Voleff=0.9              %Volumetric Efficency
AF=14.7                 %Air/Fuel Ratio

Air=rho*(Stroke*((pi)*(Bore^2))/4)*Combeff/1000000      %Air Mass (kg)
Fuel=Air/AF                                             %Fuel Mass (kg)
Pcomp=1;                                                %Ambient Air Pressure (Bar)   
Pcom=100000;                                            %Ambient Air Pressure (N/m^2)
Vc=Clearance*1000                                       %Clearance Volume (mm^3)
Rad=Stroke/2                                            %Crank Radius (mm)
Comp=((Stroke*((pi)*(Bore.^2))/4)+Vc)/Vc                %Compression Ratio
Crank=-180:1:180;                                       %Crank Angle Range
s=Rad*cosd(Crank)+((Stroke.^2)-(Rad^2)*sind(Crank).*sind(Crank)).^(1/2);  %s  
Vol=(Vc+((((pi)*(Bore.^2))/4))*(Stroke+Rad-s))/1000;    %Cylinder Volume (cc)

% Calculation Loop for Non-Combustion Cylinder Pressure
for i=1:1:360;
    if (i<180)
    n=1.33;
    else 
    n=1.27;
    end
    Pcomp(i+1)=Pcomp(i)*(Vol(i)/Vol(i+1))^n;
end

% Definition of Wiebe Function
Wiebe=1-exp(-a*((Crank-(Ignition+IgnDel))/bdur).^(m+1));
for i=1:1:360;
    if (i<(180+(Ignition+IgnDel)))
    stop=0;
    else if (i>(180+((Ignition+IgnDel)+bdur)))
    stop=0;
    else 
    stop=1;
        end
    MFB(i)=(Wiebe(i+1)-Wiebe(i))*stop;                  %Calculation of Mass Fraction Burned
    end
end

% Plot Wiebe Function
figure(1)
clf
plot(Crank,Wiebe,'k')
axis([-20 40 0 1])
xlabel('Angle (degrees)');
ylabel('Wiebe');
legend('Wiebe Function');

% Calculation for dQr (Heat Released per Degree)
dQr=MFB*LHV*Fuel*Voleff; 

% Calculation Loop for Combustion Cylinder Pressure
for i=1:1:360;
    if (i<180) %180 degrees will be TDC
    n=1.33;
    else 
    n=1.27;
    end
    Pcom(i+1)=(Pcom(i)*(Vol(i)/Vol(i+1))^n)+(((dQr(i)*(gamma-1))/((Vol(i+1)/1000))));
    Combustion=Pcom/100000;                              %Convert from N/m^2 to Bar
end

% Plot figures for Combustion and Without Combustion
figure(2)
plot(Crank,Combustion,'k',Crank,Pcomp,'b--')
xlabel('Crank Angle (degrees)');
ylabel('Pressure (Bar)');
legend('Combustion','No Combustion')