Otto cycle program in MATLAB

% Program to get PV diagram for Otto cycle 
% Roop Sannigrahi
clear all
close all
clc

%inputs
gamma = 1.4;
t3 = 2300;
 
%State variables
p1 = 101325;
t1 = 500;

%Engine geometric parameters
bore = 0.1;
stroke = 0.1;
connecting_rod = 0.15;
cr = 12;

%Calculating the swept volume and clearance volume
v_swept = (pi/4)*(bore^2)*stroke;
v_clearance = v_swept/(cr-1);

v1 = v_swept + v_clearance;
v2 = v_clearance;

%state variables at state point 2
%p2*v2^gamma = p1*v1^gamma (along the isnetropic process pv^gamma is
%constant)
%p2 = p1*(v1/v2)^gamma = p1*(cr)^gamma

%calculation
p2 = p1*(cr^gamma);

%Using ideal gas eqn 
%(p1*v1)/t1 = (p2*v2)/t2  | t2 = ((p2*v2)*t1)/(p1*v1)
t2 = ((p2*v2)*t1)/(p1*v1);

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


%state variables at point 3
v3 = v2;

%Using ideal gas eqn 
% (p3*v3)/t3 = (p2*v2)/t2  | p3 = p2*(t3/t2)
p3 = p2*(t3/t2);

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

%state variables at point 4 
v4 = v1;

% p3*(v3^gamma)=p4*(v4^gamma)   | p4 = p3*(v3/v4)^gamma
%(p3*v3)/t3 = (p4*v4)/t4  | t4 = ((p4*v4)*t3)/(p3*v3)
p4 = p3*((v3/v4)^gamma); 
t4 = ((p4*v4)*t3)/(p3*v3);

% Thermal efficiency for otto cylce is given by 
thermal_efficiency = (1 - ((t4-t1)/(t3-t2)))*100;
T = ['Thermal Efficiency = ',num2str(thermal_efficiency),'%'];
disp(T)
figure(1)
hold on
title('PV diagram of Otto cycle')
xlabel('Volume')
ylabel('Pressure')
plot(v1,p1,'*','color','r')
plot(V_compression,P_compression)
plot(v2,p2,'*','color','r')
plot([v2 v3],[p2 p3])
plot(v3,p3,'*','color','r')
plot(V_expansion,P_expansion)
plot(v4,p4,'*','color','r')
plot([v4 v1],[p4 p1])



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


a = stroke/2;
R = connecting_rod/a;

v_swept = (pi/4)*(bore^2)*stroke;
v_clearance = v_swept/(cr-1);

theta = linspace(start_crank,end_crank,100);

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)).*v_clearance);

end

Otto cycle program

In the challenge it was asked to write a program that would give the output as the PV diagram of the otto cycle and the overall thermal efficiency of the cycle.

So first we need some inputs to start the program , here I have considered the constant inputs as

Gamma(γ) = 1.4 and temperature T₃= 2300K

Then we had to define some state variable (these variables represent a particular condition)

So P₁=101325 Pa i.e atmospheric pressure and T₁=500K

Then we moved on to the geometric parameters of the engine this will be very important in calculating the volume function. We took basic parameters such as bore =0.1 , stroke = 0.1 , connecting rod length = 0.15 and compression ratio (cr) = 12 . Using these parameters we can obtain clearance volume , swept volume , V₁ and V₂ i.e volumes at state point 1 and 2.

Now to calculate state variables at state point 2 we will use the following relation P₂*(V₂)^γ= P₁*(V₁)^γ. Using this relation we get the value of P₂. For calculating T₂ we will use the following relation

Now at this point we need compression volume and pressure. So for compression volume (V_compression) we will use the following formula:

V_c is clearance volume

r_c is compression ratio

θ is the crank angle

R is the ratio of connecting rod length/crank pin radius

The above relation gives us the changes in the volume in combustion chamber.

Here I have used the function command and created a function known as engine_kinematics(). What this does is, it gives an ease of passing different values for the input parameters of the function.

In the function engine_kinematics I used the above formula and broke it down in different terms to simplify it. For getting θ I have used linspace command to get the angles from start of the crank to the end of crank rotation. so the θ that we take is from 180° to 0° as it is compressing so the volume is reducing. so we obtain V_compression

At state point 3 we use this relation to get the state variable P₃ . V_{2 =}V₃ as it is an isochoric process i.e no change in volume.

Also along the isentropic process there is a constant P₁*(V₁)^γ = constant this is a constant in isentropic process we can compression pressure by P_compression = constant/(V_compression)^γ

We can use the engine_kinematics function to find V_{expansion a} as we did in V_compression but here the θ will be from 0° to 180° as there is expansion happening . Similarly for P_expansion we will use the relation P₃*(V₃)^γ =constant

P_expansion = constant/(V_expansion)^γ

To find the state variables at state 4 we use the relation P₃*(V₃)^γ= P₄*(V₄)^γ we get the value of P₄ using we find the value of T₄.

For thermal efficiency we can use the formula :

All the values of T_1, T_2, T_3, T₄ were obtained from above relations , this helps us finding the thermal efficiency of the cycle. So the thermal efficiency that I obtained was 62.9893% or 0.629893. for displaying this thermal efficiency I stored in a variable and used num2string() to convert it into a string to display the result.

For the plot of the PV diagram I used the different values of Pressure( P₁, P₂, P_compression , P₃ , P_expansion ,P₄) Volume (V₁, V₂, V_compression , V₃ , V_expansion , V₄) that was obtained using various relations

OUTPUT:

PV diagram:

Thermal efficiency of the cycle: