Air Standard Cycle Challenge

The following MATLAB code is for plotting PV diagram for an otto cycle and finding its thermal efficency for a given set of input parameters. The following code utilises a function defined in a separate MATLAB script to find the volume during compression stroke and volume during expansion stroke which in turn is returned to the main script after passing certain inputs in correct order.

clear all
close all
clc

%Input
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;

%Calculating 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
p2 = p1*(cr^gamma);
t2 = (t1/cr)*(p2/p1);

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;
p3 = p2*(t3/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;
p4 = p3/(cr^gamma);
v4 = v1;
t4 = t3*cr*(p4/p3); 

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

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

% Code for \'Engine_kinematics\' function
% This code has to be pasted in a separate MATLAB script and must be named as \"Engine_kinematics.m\".This script has to be saved and run atleast once before running the main code.

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


% crank pin radius
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,100);

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

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

end