MATLAB code to solve efficiency of Otto cycle and generate Pressure-Volume plot

Title : MATLAB code to solve efficiency of Otto cycle and generate Pressure-Volume plot.

Aim : To calculate the efficiency of Otto cyle.

Introduction : The spark ignition piston engines works on Otto cycle. Such engines are use in automobile. The cycle basically shows the behaviour of the gas when it undergoes compression, heat addition, expansion and heat rejection.

Otto cycle is also called as constant volume cycle. Otto cycle consists of two constant volume processes and two isentropic processs.The code to generate the efficiency and pressure- volume plot is shown below

Code:

%program to calculate the otto cycle efficiency
clc
clear all
close all
%inputs
%engine specifications
d = 0.1
l = 0.15
stroke = 0.1
cr = 9.9
gamma = 1.4
%at stage 1
p1 = 101325
t1 = 500
t3 = 2500
a = stroke/2
R = l/a

%volumes

stroke_volume = (pi/4)*d^2*stroke
clearance_volume = stroke_volume/(cr-1)
v1 = stroke_volume + clearance_volume
v2 = clearance_volume

%At stage 2
% p1v1^gamma = p2v2^gamma
%p2 = p1*(cr)^gamma

p2 = p1*(cr)^gamma

%t2 = (p2*v2*t1)/(p1*v1)

t2 = (p2*v2*t1)/(p1*v1)

% p1*v1^gamma = constant_c
%therefore, p_compression*v_compression^gamma = constant_c
%(It is required in order to trace pressure volume relation during compression process)

constant_c = p1*v1^gamma  
v_compression = engine_kine(l,d,cr,stroke,180,360);
p_compression = constant_c./(v_compression.^gamma)

%At stage 3
v3 = v2
p3 = p2*t3/t2
%p1*v1^gamma = constant_c
%therefore, p_expansion*v_expansion^gamma = constant_c
%(It is required in order to trace pressure volume relation during expansion process)

constant_c = p3*v3^gamma 
v_expansion = engine_kine(l,d,cr,stroke,0,180);
p_expansion = constant_c./(v_expansion.^gamma);

%At stage 4
v4 = v1
p4 = p3*(v2/v1)^gamma
t4 = (p4*v1*t3)/(p3*v2)
figure(1)
hold on
plot(v1,p1,'*','color','r')
plot(v_compression,p_compression)
plot(v2,p2,'*','color','r')
plot(v3,p3,'*','color','r')
plot(v_expansion,p_expansion)
plot(v4,p4,'*','color','r')
plot([v2 v3],[p2 p3])
plot([v4 v1],[p4 p1])
axis([0 0.001 0 7*10^6])
xlabel('Volume V')
ylabel('Pressure P')
%estimation of efficiency
efficiency = [1-((t4-t1)/(t3-t2))]*100

Code 2 : To generate function.

Explanation of codes : The program starts with the input parameters where,

d = cylinder bore diameter

l = connecting rod length

cr = compression ratio

a = crank radius = stroke/2

R = l/a

Process 1-2 : Isentropic compression process. We know the atmospheric pressure (p1) and temperature (t1). To calculate the pressure at point 2, use the law p1v1^gamma = p2v2^gamma. Similarly to calculate temperature (t2) use (p1v1/t1) = (p2v2/t2) as shown in code where,

v1 = clearance volume + stroke volume

v2 = clearance volume

Process 2-3 : Constant volume process hence, v2 = v3. To calculate p3 use (p2v2/t2) = (p3v3/t3) as v2 = v3 the equation becomes p3/t3 = p2/t2. Here, t3 is given and p2 gets calculated.

Process 3-4 : Isentropic expansion process. Here the values of p4 and t4 are calculated same as that of process 1-2.

Process 4-1 : The values of p1 and t1 are given

Now, the above code will just give the points. But while joing this points, the nature of lines representing compression and expansion will not be straight line. Hence a seperate function should be generated nwhich shows the volume change at evey instant.Thus the code 2 is generated to determine the volume change.

Errors faced and corrective actions :

1) As the variable constant_c is matrix it will require dot product. So the syntax must include './' instead of '/'.

2) Improper plot. This occured due to incorrect formula of isentropic process.

Corrected action : The syntax and formula is corrected.

Result :

1) Efficiency is estimated to be 60.0289

2) The p-v plot,