PV Diagram & Thermal Efficiency of Air Standard Otto Cycle

%to create a PV Diagram using required values & state points

clear all
close all
clc

% Input varialbles
gamma = 1.14;
t3 = 2100;

% Engine dimension
bore = 0.1;
stroke = 0.1;
con_rod = 0.15;
cr = 12;

%Calculation of swept & clearance volume
v_swept = (pi/4)*bore*2*stroke;
v_clearance = v_swept/(cr-1);

% Standard varibles at point 1
v1= v_swept + v_clearance;
p1 = 1.01325;
t1 = 300;

% State variables at point 2
v2 = v_clearance;
p2 = p1*(v1/v2)^gamma;
t2 = (p2*v2*t1)/(p1*v1);

%for the values of Pressure & Volume b/w p1,v1 & p2,v2 respectively
constant_c = p1*v1^gamma;
V_compression = engine_kinematics(bore, stroke, con_rod, cr, 180, 0);
P_compression = constant_c./V_compression.^gamma;

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

%for the values of Pressure & Volume b/w p3,v3 & p4,v4 respectively
constant_c = p3*v3^gamma;
V_expansion = engine_kinematics(bore, stroke, con_rod, cr, 0, 180);
P_expansion = constant_c./V_expansion.^gamma;

%State varibles at point 4
v4 = v1;
p4 = p3*(v3/v4)^gamma;
t4 = (p4*v4*t3)/(p3*v3);

%Thermal Efficiency of Otto cycle
%nth = 1-(t4-t1/t3-t2)
%nth = 1-(v2/v1)^gamma-1
%nth = 1-(1/cr^gamma-1)
nth1 = 1-((t4-t1)/(t3-t2))
nth2 = 1-((v2/v1)^(gamma-1))
nth3 = 1-(1/(cr^(gamma-1)))

%plotting the PV diagram
figure(1)
hold on
plot(v1, p1, \'+\', \'color\', \'b\')
plot(V_compression, P_compression, \'color\', \'r\', \'linewidth\', 1.5)
plot(v2, p2, \'+\', \'color\', \'b\')
plot(v3, p3, \'+\', \'color\', \'b\')
plot([v2 v3], [p2 p3], \'color\', \'r\', \'linewidth\', 1.5)
plot(V_expansion, P_expansion, \'color\', \'r\', \'linewidth\', 1.5)
plot(v4, p4, \'+\', \'color\', \'b\')
plot([v4 v1], [p4 p1], \'color\', \'r\', \'linewidth\', 1.5)

%Labelling x and y axes
xlabel(\' Volume \');
ylabel(\' Pressure \');

------------------------------------------------------------------------------------------

%function engine_kinematics

function [V] = engine_kinematics(bore, stroke, con_rod, cr, start_crank, end_crank)

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

v_swept = (pi/4)*bore*2*stroke;
v_clearance = v_swept/(cr-1);

theta = linspace(start_crank,end_crank,100);

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

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

end