Challenge Week 2 - Otto Cycle PV Diagram and thermal efficiency of engine

% function file

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

% calculates values of V and stores them in arrays

% engine geometric parameters

a = stroke/2; % crank pin radius
R = conn_rod/a;

% swept volume
Vs = (pi/4)*bore^2*stroke;

% clearance volume
Vc = Vs/(cr - 1);

% change of volume in a combustion chamber
theta = linspace(start_crank,end_crank,100); % for rotation during the complete cycle of 4 strokes

V = (1 + (cr-1)*(R + 1 - cosd(theta) - (R^2 - sind(theta).^2).^(1/2))./2).*Vc;

end

% main file

clear all
close all
clc

%iinputs
gamma = 1.4;

% state variables
P1 = 101325;
T1 = 500;
T3 = 2300; % assumed

% engine geometric parameters
bore = 0.1;
stroke = 0.1;
conn_rod = 0.15;
cr = 12;

% swept volume
Vs = (pi/4)*bore^2*stroke;

% clearance volume
Vc = Vs/(cr - 1);

% total volume of cylinder at state 1
V1 = Vs + Vc;

% total volume of cylinder at state 2
V2 = Vc;

% isentropic compression
P2 = P1*(V1/V2)^gamma;

% ideal gas law
T2 = P2*V2*T1/(P1*V1);

V_compression = piston_kinematics(bore,stroke,conn_rod,cr,180,0);

Constant1 = P1*V1^gamma; % isentropic compression process
P_compression = Constant1./V_compression.^gamma;

% due to constant volume heat addition,
V3 = V2;

% ideal gas law
P3 = P2*T3/T2;

V_expansion = piston_kinematics(bore,stroke,conn_rod,cr,0,180);

Constant2 = P3*V3^gamma; % isentropic expansion process
P_expansion = Constant2./V_expansion.^gamma;

% volume at state 4
V4 = V1;

% isentropic expansion
P4 = P3*(V3/V4)^gamma;

hold on
plot(V1,P1,\'*\',\'color\',\'r\')
plot(V_compression,P_compression,\'linewidth\',1,\'color\',\'b\')
plot(V2,P2,\'*\',\'color\',\'r\')
plot([V2 V3],[P2 P3],\'linewidth\',1,\'color\',\'b\')
plot(V3,P3,\'*\',\'color\',\'r\')
plot(V_expansion,P_expansion,\'linewidth\',1,\'color\',\'b\')
plot(V4,P4,\'*\',\'color\',\'r\')
plot([V1 V4],[P1 P4],\'linewidth\',1,\'color\',\'b\')
grid on

xlabel(\'Volume\')
ylabel(\'Pressure\')

thermal_efficiency = 1 - (1/(V4/V3)^(gamma-1))

% I thought we have to create a single file for added difficulty since we made exactly the same thing by watching the video :D