Obtaining the PV diagram and thermal efficiency of an Otto cycle spark ignition piston engine using MATLAB

A MATLAB code has been written to obatain the PV diagram and thermal efficiency of an Otto cycle spark ignition piston engine.

The Otto cycle is a description of what happens to a mass of gas as it is subjected to changes of pressure, temperature, volume, addition of heat, and removal of heat.

The Inputs that will be given are:

1.) t1 and p1 which are the temperature and pressure values at the beginin of the compression stoke before the isentropic compression takes place

2.) t3, v3 and p3 which are the temperature, volume and pressure values at ignition before isentropic expansion or power stroke takes place.

3.) The specific heat ratio, gamma.

4.) The geometric parameters of the engine like: the bore diameter, the stroke length, the connecting rod length and the compression ratio (cr).

The outputs obatined are:

1.) The PV diagram of the cycle

2.) The thermal efficiency of the engine

3.) The values temperature, pressure and volume at points 2 and 4.

The code uses a function called engine_kinamatics that used to find the volume at compression and expansion. The code and the function are given below:

clear all
close all
clc

%Inputs
gamma = 1.4; %Specific Heat ratio Cp/Cv (i.e: Specific Heat at Constant Pressure/ Sprcific Heat at Constant Volume)
t3 = 2300; %K

%State Variables
p1 = 101325; %Pa
t1 = 500; %K

%Engine geometeric parameters in meters
bore = 0.1;
stroke =  0.1;
con_rod = 0.15; %Connecting Rod Length
cr = 10; %Compression Ratio

%Calculating the swept volume  and the clearence volume which are giving by
%the formaulas shown below:
v_swept = (pi/4)*bore.^2*stroke;
v_clearence = v_swept/(cr-1);

v1 = v_swept + v_clearence;
v2 = v_clearence;

%State variables at state point 2 
%p2v2.^gamma = p1v1.^gamma (As it is a reversible adiabatic process)
%p2 = p1*(v1/v2)^gamma = p1*(cr)^gamma
p2 = p1*cr.^gamma;

%p1v1/t1 = p2v2/t2 | t2 = p2*v2*t1/(p1*v1)|From ideal gas law PV=nRT=const
t2 = p2*v2*t1/(p1*v1); 

constant_c = p1*v1.^gamma;

V_compression = engine_kinematics(bore,stroke,con_rod,cr,180,0); %During compression piston moves from 180 to 0 degrees (Refer to the funcion explanation)

P_compression = constant_c./V_compression.^gamma;

% State Variable at state point 3
v3 = v2; %Contant volume combustion/heat addition process (V2 and V3 are on the same line in the PV Diagram)

%p3v3/t3 = p2v2/t2 | p3 = p2*t3/t2 | From ideal gas law
p3 = p2*t3/t2;

constant_c = p3*v3.^gamma;
V_expansion = engine_kinematics(bore,stroke,con_rod,cr,0,180); % During exapnsion piston moves from 0 to 180 degrees.

P_expansion = constant_c./V_expansion.^gamma;

%State variables at state point 4 | Contant volume heat rejection process
v4 = v1;

%p3v3.^gamma = p4v4.^gamma | p4 = p3(v3/v4)^gamma
p4 = p3*(v3/v4)^gamma;

figure(1)
hold on
plot([v1 v4], [p1 p4],'color','b','linewidth',2)
plot([v2 v3], [p2 p3],'color','b','linewidth',2)
plot(V_compression, P_compression,'color','b','linewidth',2)
plot(V_expansion, P_expansion,'color','b','linewidth',2)
title('Otto cycle PV Diagram')
xlabel('Volume')
ylabel('Pressure')
grid on 

Now for the function engine_kinematics, that is used to find the volume inside the bore at compression and expansion.

The formula given above detemines how the volume inside the combustion chamber changes. It has been derived from purely geometric relationships.

theta,θ = Crank angle
rc = compression ratio,cr
R = connecting rod length / crank pin radius
Vc = clearence volume

This formula in the function code is split up into 3 terms, term 1, term 2 and term 3.

function [V] = engine_kinematics(bore,stroke,con_rod,cr,start_crank_angle,end_crank_angle)

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

v_s = (pi/4)*bore.^2*stroke; %Swept Volume
v_c = v_s/(cr-1); %Clearence Volume

theta = linspace(start_crank_angle,end_crank_angle,100);

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

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

This function is called in the main code to find the volume at compression and expansion.

The thermal efficiency of an Otto cycle is given by:

where CR is the compression ratio and K is gamma, the specific heat ratio.

Errors:

I had intially tried to plot the entire PV diagram using in a single plot command as shown below:

The resulting graph obatined was as shown below:

The above graph has unnecessary straight lines joining points 3 to 4 and 1 to 2

This mistake was rectified using seperate plot commands as shown in the code provided.

Results: