Otto cycle

%Air Standard Cycle(Otto cycle)

clear all
close all
clc

%Inputs

gamma = 1.4;
t3 = 2300;

%State variables

p1 = 101325;
t1 = 500;

%Engine geometric parameters
bore = 0.1;
stroke = 0.1;
con_rod = 0.15;
cr = 12;

%claculating thermal efficiency

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

%Calculating the swept 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 at state point 2
%p1v1^gamma = p2v2^gamma
%p2 = p1*(v1/v2)^gamma = p1*cr^gamma

p2 = p1*cr^gamma;

%p1v1/t1 = p2v2/t2 ~ t2 = p2*v2*t1/(p1*v1)

t2 = p2*v2*t1/(p1*v1);
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 state point 3

v3 = v2;

%p2v2/t2 = p3v3/t3 ~ p3 = (p2*v2*t3)/(v3*t2)

p3 = (p2*v2*t3)/(v3*t2);
constant_c = p3*v3^gamma;
V_expansion = engine_kinematics(bore, stroke, con_rod, cr, 0, 180);
P_expansion = constant_c./V_expansion.^gamma;

%State variables at state point 4

v4 = v1;

%p3v3^gamma = p4v4^gamma ~ p4 = p3*(v3/v4)^gamma

p4 = p3*(v3/v4)^gamma;

t4 = (p4*v4*t3)/(p3*v3);

%plotting

figure(1)
hold on
plot(v1,p1,\'*\',\'color\',\'r\')
plot(V_compression,P_compression,\'linewidth\',2)
plot(v2,p2,\'*\',\'color\',\'r\')
plot(v3,p3,\'*\',\'color\',\'r\')
plot(V_expansion,P_expansion,\'linewidth\',2)
plot(v4,p4,\'*\',\'color\',\'r\')
plot([v1,v4],[p1,p4],\'color\',\'g\',\'linewidth\',2)
plot([v2,v3],[p2,p3],\'color\',\'g\',\'linewidth\',2)

%labeling

xlabel(\'volume\')
ylabel(\'pressure\')

.....................................................................................................................

%engine_kinematics function.....

function [V] = engine_kinematics(bore, stroke, con_rod, cr, start_angle, end_angle)

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

theta = linspace(start_angle,end_angle,500);

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

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

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

end