p-v diagram of an otto cycle and its thermal efficiency

        CREATING THE OTTO CYCLE P-V PLOT AND DETERMINING THR THERMAL EFFICIENCY


%1) Your code should create a PV diagram
%Otto cycle plot

clear all
close all
clc

%Inputs

gamma = 1.4;
t3 = 2300

%State variables

p1 = 101325;
t1 = 500;

%Engine geomatric parameter

bore = 0.1;
stroke = 0.1;
con_rod = 0.15;
cr = 12;


%Calculating the swept volume and the clearence volume

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

v1 = v_swept + v_clearence
v2 = v_clearence

%State variable at point 2
%p2v2^gamma=p1v1^gamma

p2 = p1*cr^gamma

%p1v1/t1 = p2v2/t2 | t2 = p2v2t1/p1v1

t2 = p2*v2*t1/(p1*v1)
v_compression = engine_kinematics(bore,stroke,con_rod,cr,180,0)
c_const1 = p1*v1^gamma;
%p_compression*v_vompression^gamma=c_const

p_compression = c_const1./v_compression.^gamma;


%State variable at point 2

v3 = v2

%p3v3/t3 = p2v2/t2 | p3 = p2v2t3/t2v3 

p3 = p2*v2*t3/(t2*v3)

%State variable at point 4

v4 = v1

%p3v3^gamma = p4v4^gamma | p4 = p3v3/v4

p4 = p3*(v3/v4)^gamma
v_expansion = engine_kinematics(bore,stroke,con_rod,cr,0,180)
c_const2 = p3*v3^gamma;
%p_compression*v_vompression^gamma=c_const

p_expansion = c_const2./v_expansion.^gamma;

%2) You should output the thermal efficiency of the engine.
%Thermal efficiency

thermal_efficiency = 1-(1/cr^(gamma-1))



figure(1)
hold on 
plot(v1,p1,\'*\',\'color\',\'r\') %Isentropic compression
plot(v_compression,p_compression)
plot(v2,p2,\'*\',\'color\',\'r\')%Constant volume heat addition
plot([v2 v3],[p2 p3],\'color\',\'r\')
plot(v3,p3,\'*\',\'color\',\'r\') %Isentropic expansion
plot(v_expansion,p_expansion)
plot(v4,p4,\'*\',\'color\',\'r\') %Constant volume heat rejection
plot([v4 v1],[p4 p1],\'color\',\'c\')
xlabel(\'volume\')
ylabel(\'pressure\')
title(\'P-V diagram of otto cycle\')


                            CREATING THE FUNCTION


%Pston kinematics

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


%Inputs



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

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

theta = linspace(start_crank,end_crank,360);

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

%Final formula

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

end