Otto Cycle and Piston Kinematics using Matlab

An Otto cycle of an Air Standard cycle consists of two adiabatic processes and two isochoric processes.

1 - 2: Adiabatic Compression
2 - 3: Constant Volume heat addition
3 - 4: Adiabatic Expansion
4 - 1: Constant Volume heat rejection

Following are state variables:

P = Pressure, V = Volume, T = Temperature, cr = Compression ratio, γ = gamma(ratio of heat capacity), η = Thermal efficiency of the engine

Geometric parameters such as bore, stroke and connecting rod length have been assumed in the following program.

close all
clear all
clc

%Input variables
gamma = 1.4;
t3 = 2300;

%To Plot Otto cycle
%Engine Geometric parameters
bore = .12;
stroke = .08;
cr = 12;
con_rod = .15;

%State variables
p1 = 101325;
t1 = 300;

%Calculating Swept Volume and Clearance Volume
v_swept = (pi/4)*bore^2*stroke;
v_clearance = v_swept/(cr-1);

v1 = v_swept + v_clearance;
v2 = v_clearance;


%State variables for state point 1 p1,v1,t1 are known

%State variables for state point 2
p2 = p1*(cr)^gamma;
%p1*v1/t1 = p2*v2/t2 | t2 = (p2*v2*t1)/(p1*v1)
t2 = (p2*v2*t1)/(p1*v1);

constant_c = p1*v1^gamma
v_compression = engine_kinematics(bore,stroke,cr,con_rod,180,0);
p_compression = constant_c./v_compression.^gamma;


%State variables for state point 3
v3 = v2;
%p2*v2/t2 = p3*v3/t3 | p3 = p2*(t3/t2) 
p3 = p2*(t3/t2);

constant_c = p3*v3^gamma
v_expansion = engine_kinematics(bore,stroke,cr,con_rod,0,180);
p_expansion = constant_c./v_expansion.^gamma;


%State variables for state point 4
v4 = v1;
%p3*v3^gamma = p4*v4^gamma | p4 = p3*(v3/v4)^gamma
p4 = p3*(v3/v4)^gamma;

%Plotting
figure(1)
hold on
plot(v1,p1,\'*\',\'color\',\'r\')
plot(v_compression,p_compression)
plot(v2,p2,\'*\',\'color\',\'r\')
plot([v2,v3],[p2,p3])
plot(v3,p3,\'*\',\'color\',\'r\')
plot(v_expansion,p_expansion)
plot(v4,p4,\'*\',\'color\',\'r\')
plot([v4,v1],[p4,p1])
axis([0 10^-3 0 10^7])
xlabel(\'Volume\')
ylabel(\'Pressure\')
title(\'PV Diagram of Otto Cycle\')



%Calculating thermal efficiency for otto cycle using expression eta = 1 - cr^(1-gamma)
eta = 1 - cr^(1-gamma)

Engine Kinematics function called in the above program is follows:

function [V] = engine_kinematics(bore,stroke,cr,con_rod,crank_start,crank_end)


a = stroke/2;
R = con_rod/a;

v_s = (pi/4)*bore^2*stroke;
v_c = v_s/(cr-1);


theta = linspace(crank_start,crank_end,180);

term1 = .5*(cr-1);
term2 = (R + 1 - cosd(theta));
term3 = (R^2 - sind(theta).^2).^.5;

V = (1+ term1*(term2-term3)).*v_c;

plot(theta,V,\'linewidth\',3)

end