Plotting Otto cycle using MATLAB

AIM:
Writing a program in Matlab to calculate the efficiency of otto cycle and to plot PV diagrams for the same.

THEORY:
An Otto cycle is an idealized thermodynamic cycle that describes the functioning of a typical spark ignition piston engine. It is the thermodynamic cycle most found in automobile engines. In Otto cycle, enough work is produced by combusting the fuel to propel an automobile and its occupants in the environment.

• Process 0–1 a mass of air is drawn into piston/cylinder arrangement at constant pressure.
• Process 1–2 is an adiabatic (isentropic) compression of the charge as the piston moves from bottom dead centre (BDC) to top dead centre (TDC).
• Process 2–3 is a constant-volume heat transfer to the working gas from an external source while the piston is at TDC. This process is intended to represent the ignition of the fuel-air mixture and the subsequent rapid burning.
• Process 3–4 is an adiabatic (isentropic) expansion also known as power stroke.
• Process 4–1 completes the cycle by a constant-volume process in which heat is rejected.
• Process 1–0 the mass of air is released to the atmosphere in a constant pressure process.

FORMULA:

• p1v1^gamma=p2v2^gamma then p2=p1(v1/v2)^gamma
• v1/v2= Compression ratio
• p1v1/t1=p2v2/t2
• p2v2/t2=p3v3/t3
• p3v3^gamma=p4v4^gamma
• p3v3/t3=p4v4/t4
• Efficiency=1-(1/cr^(gamma-1))

All the above-mentioned formulae are used for writing a code in Matlab.

MATLAB PROGRAMME:

%Inputs to the problem
p1=101325;
t1=500;
t3=2300;
gamma=1.4;

%Engine parameters
b=0.1; %bore
s=0.1; %stroke
con_rod=0.15; %connecting rod length
cr=12; %compression ratio

vs=(pi/4)*b^2*s; %Swept volume
vc=vs/(cr-1); %Clearance volume

v1=vc+vs;
v2=vc;
p2=p1*(cr)^gamma; %p1v1^gamma=p2v2^gamma then p2=p1(v1/v2)^gamma and v1/v2=cr
t2=(p2*v2*t1)/(p1*v1); %p1v1/t1=p2v2/t2
constant_c=p1*v1^gamma;
[V_compression]=engine_kinematics(b,s,con_rod,cr,180,0);
P_compression=constant_c./V_compression.^gamma;

v3=v2; %constant volume process
p3=(p2*v2*t3)/(t2*v3); %p2v2/t2=p3v3/t3
constant_c1=p3*v3^gamma;
[V_expansion]=engine_kinematics(b,s,con_rod,cr,0,180);
P_expansion=constant_c1./V_expansion.^gamma;

v4=v1; %constant volume process
p4=p3*(v3/v4)^gamma; %p3v3^gamma=p4v4^gamma
t4=(p4*v4*t3)/(p3*v3); %p3v3/t3=p4v4/t4

%plotting
hold on
plot(v1,p1,'*','color','r');
plot(V_compression,P_compression,'color','b','linewidth',1.5);
plot(v2,p2,'*','color','r');
plot([v2,v3],[p2,p3],'color','b','linewidth',1.5);
plot(v3,p3,'*','color','r');
plot(V_expansion,P_expansion,'color','b','linewidth',1.5);
plot(v4,p4,'*','color','r');
plot([v4,v1],[p4,p1],'color','b','linewidth',1.5);
xlabel('Volume');
ylabel('Pressure');
title('PV diagram of Otto Cycle');

%efficiency calculation
eff=1-(1/cr^(gamma-1));
efficiency=eff*100+"%";

Volume trace function:

function [V]=engine_kinematics(b,s,con_rod,cr,start_crank,end_crank)

a=s/2;
R=con_rod/a;

vs=(pi/4)*b^2*s; %Swept volume
vc=vs/(cr-1); %Clearance volume

theta=linspace(start_crank,end_crank,180);
t1=0.5*(cr-1);
t2=R+1-cosd(theta);
t3=(R^2-sind(theta).^2).^0.5;

V=(1+t1*(t2-t3)).*vc;
end

RESULTS:

DISCUSSION OF RESULTS/CONCLUSION:
It is seen that the volume changes according to the engine kinematics relations and not linearly. The efficiency of the otto cycle was calculated and it is shown below.