Air standard Cycle

%Air standard Cycle [Otto cycle] 
%Create a PV diagram of otto cycle
%Thermal efficiency of the engine for various compression ratio vales

close all
clear all
clc

%inputs
gamma=1.4
t3=2000

%physical parameters
bore_dia=0.1
stroke_leng=0.1
con_rod=0.1
cr=12

%calculate swept and clearance volume
v_swept=(pi/4)*bore_dia^2*stroke_leng
v_clearance=v_swept/(cr-1)

v1=v_swept+v_clearance
v2=v_clearance

%for point1
p1=101325
t1=500

constant=p1*v1^gamma
v_compression=engine_kinematics(bore_dia,stroke_leng,con_rod,cr,180,0)
p_compression=constant./v_compression.^gamma

%for point2
p2=p1*(cr)^gamma
t2=p2*v2*t1/(p1*v1)


%for point3
v3=v2
p3=p2*t3/t2

Constant=p3*v3^gamma
v_expansion=engine_kinematics(bore_dia,stroke_leng,con_rod,cr,0,180)
p_expansion=Constant./v_expansion.^gamma

%for point4
v4=v1
p4=p3*(1/cr )^gamma
t4=p4*v4*t3/(p3*v3)

%create plot of PV Diagram
figure(1)
hold on
plot(v1,p1,'*')
plot(v_compression,p_compression)
plot(v2,p2,'*')
plot(v3,p3,'*')
plot(v_expansion,p_expansion)
plot(v4,p4,'*')
plot([v2,v3],[p2,p3])
plot([v4,v1],[p4,p1])
xlabel('volume')
ylabel('pressure')
%create plot of compression ratio vs thermal efficiency
figure(2)
%inputs
comp_ratio=1:1:20;
thermal_eff=(1-(1./(comp_ratio).^(gamma-1)))*100
plot(comp_ratio,thermal_eff)
xlabel('compression ratio')
ylabel('thermal efficiency')

Subprogram used

function [v]=engine_kinematics(bore_dia,stroke_leng,con_rod,cr,start_crank,end_crank)
a=stroke_leng/2;
R=con_rod/a;
%swept and clearance volume
v_swept=(pi/4)*bore_dia^2*stroke_leng;
v_clearance=v_swept/(cr-1);

theta =linspace(start_crank,end_crank,180);
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

In the below screenshot of output

figure 1 represents pv diagram

figure 2 represents plot of compression ratio vs thermal efficiency