MATLAB code to plot P-V diagram of Otto cycle and calculate it's thermal efficiency

Otto cycle processes:-

Process 1-2 (Isentropic compression)

Process 2-3 (Contant volume heat addition)

Process 3-4 (Isentropic expansion)

Process 4-1 (Contant volume heat rejection)

'main_script' has all the inputs and calculates state variables i.e pressure, volume and temeprature. Then, it calls 'v_at_crank_angle.m' function to get cylinder volume at various given crank angles. Further, after 'v_at_crank_angle.m' returns array of cylinder volume, it calculates in process (1-2) and (3-4) state variables at various crank angles. Finally, it plots the P-V diagram for the Otto cycle, calculates and displays the themal efficiency.

Themal efficiency, `eta = 1 - (1/(r^(gamma-1)))`

% main_script.m


% Calculating thermal efficiency of Otto cycle and
% plotting it's P-V diagram
% Clearing
close all;
clear;
clc;
% Inputs
p_1 = 1; %Pressure at state 1 (atm) (Bottom dead centre (BDC))
t_1 = 500; % Temperature at state 1 (K)
bore_dia = 0.1; %Cylinder bore diameter (m)
l = 0.1; %Cylinder stroke (m)
cr = 12; %Compression ratio
gamma = 1.4;
t_3 = 2300; %Temerature at state 3 that is also the maximum allowable temperature (Kelvin)
conec_rod = 0.35; % Connecting rod length (m)
%Calculation
v_stroke = (double(pi) * bore_dia^2 * l)/4;
v_clearance = v_stroke/(cr-1);
v_1 = v_clearance + v_stroke;
v_2 = v_clearance;
v_3 = v_2;
p_2 = p_1 * (v_1/v_2)^gamma;
t_2 = t_1 * (v_1/v_2)^(gamma-1);
p_3 = p_2 * (t_3/t_2);
v_4 = v_1;
p_4 = p_3 * (v_3/v_4)^gamma;
% Plotting
figure(1);
hold on;
grid on;
plot(v_1 , p_1 , 'r.' , 'MarkerSize' , 10);
plot(v_2 , p_2 , 'r.' , 'MarkerSize' , 10);
plot(v_3 , p_3 , 'r.' , 'MarkerSize' , 10);
plot(v_4 , p_4 , 'r.' , 'MarkerSize' , 10);
xlabel('Volume (m^3)');
ylabel('Pressure (atm)');
title('P-V Diagram of Otto cycle');
% Process 1-2 (Isentropic compression)
v_12 = v_at_crank_angle(v_clearance, l, cr, conec_rod, 180, 360);
p_12 = p_1 * ((v_1./v_12).^gamma);
plot(v_12 , p_12 , 'b.-' , 'MarkerSize' , 3);
% Process 2-3 (Contant volume heat addition)
plot([v_2 v_3] , [p_2 p_3] , 'b.-' , 'MarkerSize' , 3);
% Process 3-4 (Isentropic expansion)
v_34 = v_at_crank_angle(v_clearance, l, cr, conec_rod, 0, 180);
p_34 = p_3 * ((v_3./v_34).^gamma);
plot(v_34 , p_34 , 'b.-' , 'MarkerSize' , 3);
% Process 4-1 (Contant volume heat rejection)
plot([v_4 v_1] , [p_4 p_1] , 'b.-' , 'MarkerSize' , 3);
% Thermal efficiency of Otto cycle
eff = 1 - (1/cr^(gamma-1));
fprintf('Thermal efficiency = %f' , eff);

'v_at_crank_angle.m' is a function that calculates cylinder volume as a function of crank angle.

% v_at_crank_angle.m

function [v] = v_at_crank_angle(v_compression, l, cr, conec_rod, start_angle, end_angle)
%Inputs
% v_compression = Compressed volume = Clearance volume
% cr = compression ratio
% start_angle = starting crank angle
% end_angle = ending crank angle
% At theta = 0 degree v = v_compression
% At theta = 180 degrees v = v_compression + v_stroke
a = l/2; % Crank pin radius (m)
R = conec_rod/a; % Crank ratio
theta = linspace(start_angle, end_angle, 200); %crank angle (degrees)
term1 = 0.5 * (cr - 1);
term2 = R + 1 - cosd(theta);
term3 = (R^2 - (sind(theta)).^2).^0.5;
v = (1 + term1 * (term2 - term3)) .* v_compression;
end

Output:-