To demonstrate the working of an otto cycle & calculate its thermal efficiency.

Aim:- Write a code that can solve otto cycle and make plots of it, also calculate the thermal efficiency of the engine

The Otto cycle describe how heat engines turn gasoline into motion. It is most used in automobile engines.

1. a) Process 1-2 = Isentropic compression process.
2. b) Process 2-3 = Constant volume heat addition process
3. c) Process 3-4 = Isentropic expansion process
4. d) Process 4-1 Constant volume heat rejection process.

Consider the mass 1 kg

* Heat addition at constant volume is given by: - Qa = mCv(t3-t2) KJ/KG where Cv is specific heat

* Heat rejection at constant volume is given by: - Qr = mCv(t4-t1) KJ/KG where Cv is specific heat

* Work done by engine

Wd = Qa -Qr

Wd = Cv(t3-t2)-Cv(t4-t1)

* Efficiency of otto cycle is given n =Wd/Qa

n = Cv(t3-t2)-Cv(t4-t1)/Cv(t3-t2)

n = 1-((t4-t1)/(t3-t2)) .....eq 1

for the reversible adiabatic process 3-4 and 1-2 we can write

t4/t3=((v3/v4)^(gamma-1)) .....eq 2

t1/t2 = ((v1/v2)^(gamma-1)) .....eq 3

v2=v3 and v4 = v1

t4/t3 = t1/t2

t4-t1/t3-t1

cr = (v1/v2) .....eq 4

Hence,

n = 1-((1/cr)^(gamma-1))

Piston kinematics equation

How its combustion changes as function of crank angle

bore = 0.1

stroke = 0.1

cond_rod = 0.15

v_swept = pi/4*bore^2*stroke

v_clearance = v_swept/(cr-1)

a = stroke/2

R = cond_rod/2

Formula for how volume inside combustion chamber is going to change derived from pure geometric relation.

theta = 0

term = 0.5*(cr-1)

term2 = R + 1 -cosd(theta)

term3 = [R^2 - sind(theta)^2)0.5

V = [1 + term1*(term2-term3)]*v_c

Point to remember

To check the runs code correctly when theta is set to 0 then volume must correspond to clearance volume which is basically TDC volume.

We convert the above formula into function so it colud be reuse for varrying value of theta.

name of the function and the name of the file shouid be same.

function [V] = engine_kinematics(bore, stroke, cond_rod, cr, start_crank, end_crank)

% bore dia in meter
% stroke length in meter
% connecting rod length in meter
% cr = Compression ratio
% start_crank = Starting crank angle
%end_crank = Ending crank angle


a = stroke/2;   % a =ratio of connecting rod to crank pin radius
R = cond_rod/a; 

v_s = pi/4*bore^2*stroke;  % v_s = Swepth volume
v_c = v_s/(cr-1);          % v_c = Clearance volume

theta = linspace(start_crank, end_crank, 100)

% Volume inside chamber calculate by geomatric relation
% V/Vc = 1+1/2((cr-1)[(R+1-cosd(theta)-(R^2-(sind(theta))^2))

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

Using the above code we could plot the pv diagram but yet it will not follow the path which it should be , Because MATLAB code is just joining the point of compression, combustion, expansion, exhaust.

where on other side it should folloew in actual path by taking each stages in consideration, it could be achieved by deternining the constanbt k.

k=pv^gamma is constant along the isentropic process, where pressure times volume is going to be constant

constant_c = p1*v1^gamma

P_compression = constant_c./V_compression

From the above code we can determine P_compression and plot it against V_compression to follow the desired path.

P_expansion = constant_c./V_expansion

From the above code we can determine P_expansion and plot it against V_expansion to follow the desired path.

Final plot following desired paths after adding title 7 label

figure(1);
grid on
grid minor
hold on
plot(v1,p1,'color','r','linew', 2')
plot(V_compression, P_compression, 'Linew',2)
plot(v2,p2,'color','b','linew',2)
line([v3 v3], [p2 p3],'linew',2)
plot(v3,p3,'color','g', 'Linew',2)
plot(V_expansion,P_expansion,'linew',2)
line(v4,p4,'color','y','Linew',2)
plot([v4 v1], [p4 p1],'Linew',2)


%Adding the labels
xlabel('Volume in m^3')
ylabel('Pressure in pa')

%Additing the title
title('PV diagramof an IC engine')
legend('Compression', 'Combustion', 'Expansion', 'Exhaust')

Output plot

Thermal Efficiency

Final code

% Program to perform otto cycle

clear all
close all
clc

% Inputs
gamma = 1.4;
t3 = 2300;

% State variable at point 1
p1 = 101325;   % p1 = pressure at point 1 
t1 = 500;      % t1 = Temp. at point 1

% Engine geometric parameters.
bore = 0.1;          % in meter
stroke = 0.1;        % in meter
cond_rod = 0.15;     % connecting rod length in meter
cr = 12;             % cr = compression ratio       

% Calculate the Swept volume and the clearance volume.

v_s = (pi/4)*bore^2*stroke;       % v_s = Swept volume
v_c = v_s/(cr-1);                 % v_c = Clearance volume

v1 = v_s + v_c                    % v1 = Volume at point 1
v2 = v_c                          % v2 = Volume at point 2

% State variable at point 2
% p2v2^gamma = p1v1^gamma
% p2 = p1*(v1/v2)^gamma = p1*cr^gamma
p2 = p1*cr^gamma                  % p1 & p2 = pressure at point 1 & 2 

% Ideal gas law
% p1v1/t1 = p2v2/t2 t2
t2 = p2*v2*t1/(p1*v1);           % t1 & t2 = pressure at point 1 & 2

constant_c = p1*v1^gamma;
V_compression = engine_kinematics(bore, stroke, cond_rod, cr, 0,180)
P_compression = constant_c./V_compression.^gamma

% State variable at point 3
v3 = v2         % v3 = Volume at point 3

% p3v3/t3 = p2v2/t2  p3 = p2*t3/t1 

p3 = p2*t3/t2          % p3 = pressure at point 3 & t3 = temp. at point 3

constant_c = p3*v3^gamma;
V_expansion = engine_kinematics(bore, stroke, cond_rod, cr, 0,180)
P_expansion = constant_c./V_expansion.^gamma

% State variable at point 4

v4 = v1              % v4 = Volume at point 4

% p3v3^gamma = p4v4^gamma
p4 = p3*(v3/v4)^gamma       % p4 = pressure at point 4 & v4 = Volume at point 4 

figure(1);
grid on
grid minor
hold on
plot(v1,p1,'color','r','linew', 2')
plot(V_compression, P_compression, 'Linew',2)
plot(v2,p2,'color','b','linew',2)
line([v3 v3], [p2 p3],'linew',2)
plot(v3,p3,'color','g', 'Linew',2)
plot(V_expansion,P_expansion,'linew',2)
line(v4,p4,'color','y','Linew',2)
plot([v4 v1], [p4 p1],'Linew',2)


%Adding the labels
xlabel('Volume in m^3')
ylabel('Pressure in pa')

%Additing the title
title('PV diagramof an IC engine')
legend('Compression', 'Combustion', 'Expansion', 'Exhaust')

%Calculate Thermal Efficiency
thermal_efficiency = 1-1./cr^(gamma-1)
fprintf(['The thermal efficiency of the cycle is : ' num2str(thermal_efficiency)])
fprintf('n')

Conclusion

In this report otto cycle was studied and plotted using MATLAB