MATLAB program to solve an Otto cycle and to plot PV diagram of it

Otto Cycle:

Otto cycle is a thermodynamic work producing cycle. Spark ignition or petrol engine works on this cycle. This is a closed cycle which has 4 diffrent processes as follows:

Process 1-2 : Isentropic compression

Process 2-3 : Constant volume heat addition

Process 3-4 : Isentropic expansion

Process 4-5 : Constant volume heat rejection

To plot the PV diagram and calculate thermal efficiency of otto cycle I had written code in MATLAB.

So, PV diagram can be plotted if we know the values of pressure and volume at all states and during the all processes.

In this code we have to know some variables i.e.

1) Pressure and Temperature at intake

2) Maximum temperature of cycle , value of 'γ' (gamma)

3) Engine geometric parameters i.e. bore, stroke length, compression ratio, connecting rod length

MATLAB Code:

% program to plot otto cycle
clear all
close all
clc

% Given inputs
gamma = 1.4;      
t3 = 1300;         % Maximum temperature of cycle is 1300K

% Initial given state variables for state 1
p1 = 101325;       % 1 atm pressure at state 1
t1 = 200;          % temperature of intake air is 300K

% Engine geometric parameters

bore = 0.1;                   % bore diameter of cylinder is 100mm
stroke_length = 0.1;          % considering bore=stroke, hence stroke = 100mm
connecting_rod = 0.15;        % connecting rod length is 150mm
cr = 13;                      % compression ratio is 13

% Calculating swept volume
v_swept = (pi/4)*bore*bore*stroke_length;

% Calculating clearance volume
v_clearance = v_swept/(cr-1);

v1 = v_clearance + v_swept;     % volume at state 1 is sum of clearance and swept volume
v2 = v_clearance;               % volume at state 2 is equal to clearance volume

% Calculating state variables at state 2
p2 = p1*((cr)^gamma);           % Using relation p1*v1^gamma = p2*v2^gamma or p2 = p1*(cr)^gamma
t2 = (p2*v2*t1)/(p1*v1);        % using relation (p1v1/t1) = (p2v2/t2)

constant_c = p2*(v2^gamma);     %pressure*volume^gamma = constant

% calling volume_trace function to calculate different values of volume during compression 
v_compression = volume_trace(bore,stroke_length,connecting_rod,cr,180,0); 
p_compression = constant_c./(v_compression.^gamma); % Calculating different values of pressure during compression

% Calculating state variables at state 3
v3 = v2;                  % Constant volume heat addition
p3 = (p2*t3)/t2;

constant_c1 = p3*(v3^gamma);     %pressure*volume^gamma = constant

% calling volume_trace function to calculate different values of volume during expansion 
v_expansion = volume_trace(bore,stroke_length,connecting_rod,cr,0,180); 
p_expansion = constant_c1./(v_expansion.^gamma); % Calculating different values of pressure during expansion


% Calculating state variables at state 4
v4 = v1;                  % Constant volume heat rejection
p4 =p3*((v3/v4)^gamma);   % using relation p3*v3^gamma = p4*v4^gamma

% Plotting state points at all states
figure(1)
hold on
plot(v1,p1,'*','color','b');
plot(v2,p2,'*','color','b');
plot(v3,p3,'*','color','b');
plot(v4,p4,'*','color','b');

% Plotting curves to join all atate
plot(v_compression,p_compression);    % Joining state 1 to state 2
plot([v2,v3],[p2,p3]);                % Joining state 2 to state 3
plot(v_expansion,p_expansion);        % Joining state 3 to state 4
plot([v1,v4],[p1,p4]);                % Joining state 4 to state 1

xlabel('Volume (m^3)');
ylabel('Pressure (Pa)');
title('PV diagram of Otto cycle');

% Calculation of thermal efficiency of the engine

thermal_efficiency = (1 - (1/(cr)^(gamma - 1)))*100

volume_trace function to calculate various values of volume during expansion and compression

function [V] = volume_trace(bore,stroke,connecting_rod,cr,start_crank_angle,end_crank_angle)


crank_radius = stroke/2;         
R = connecting_rod/crank_radius;

vs = (pi/4)*bore*bore*stroke;     % Calculating swept volume
vc = vs/(cr-1);                   % Calculating clearance volume

theta = linspace(start_crank_angle,end_crank_angle,180);  % Creating 180 different values of theta between start and end crank angle

term_1 = 0.5*(cr-1);
term_2 = R + 1 - cosd(theta);
term_3 = (R^2 - sind(theta).^2).^0.5;

V = (1 + term_1*(term_2 - term_3)).*vc;  % volume calculations


end

OUTPUT:

PV diagram of Otto cycle

Thermal efficiency calculated by program :

By observing the plot we can see that heat addition process is far more longer than heat rejection process. Ideally heat rejection and heat addition should happen at constant volume but in practical case to have a constant volume process is not possible.

Ideal efficiency of otto cycle is 100% but due to many losses it is not possible to achieve this in practical case.