Otto cycle

close all
clear all
clc

%inputs
gamma =1.4;
t3=2300;
p1=101325;
t1=500;

%engine Geometry
bore=.1;
con_rod=.15;
stroke=.1;
cr=12;

%calculate swept volume and clearance volume

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

%state variables at point 1 

v1= v_swept+v_clearance;
v2=v_clearance;


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

%state variables at point 3
constant_c= p1*v1^gamma;
v_compress= piston_kinematics(.1,.1,.15,12,180,360)
p_compress= constant_c./(v_compress).^gamma;

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

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

constant_c= p3*v3^gamma;
v_expansion= piston_kinematics(.1,.1,.15,12,360,540)
p_expansion= constant_c./(v_expansion).^gamma;

figure(1)
hold on 
plot(v1,p1,'*','color','r');
plot(v_compress,p_compress ,'color','r');
plot(v_expansion,p_expansion ,'color','r');
plot(v2,p2,'*','color','r');
plot(v3,p3,'*','color','r');
plot(v4,p4,'*','color','r');

plot([v2 v3 ],[p2 p3]);
plot([v4 v1 ],[p4 p1]);

thermal_efficiency= 100*(1-(t1/t2))
End of main Program



Now the function code is as follows:

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


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



v_swept= (bore^2)*stroke*(pi/4);
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