Solving an ideal Otto cycle using MATLAB

Aim:

To write code for solving an ideal Otto cycle to obtain PV diagram and to calculate the thermal efficiency of the cycle using MATLAB

Introduction:

Air standard cycles are an approximate representation of the performance of Internal Combustion Engines. Otto cycle is one such air standard cycle that is used for petrol engines. Also known as "Constant Volume Cycle", it is usually represented in the form of a PV diagram (P - Pressure and V - Volume). The curve represents the pressure and volume during the four processes namely isentropic compression (1-2), Constant volume heat addition (2-3), isentropic expansion (3-4) and Constant volume heat rejection (4-1). The volume at point 2 is equal to point 3 and the volume at point 1 is equal to point 4.

Source: www.study.com

The image shows a typical Otto cycle where P and V are Pressure and Volume respectively and Q1 and Q2 indicate the constant volume heat addition and constant volume heat rejection processes.

Formulae:

1. Swept Volume `V_s` = `pi/4 * d^2 * l`

where,

`d` = Bore diameter of the cylinder

`l` = Stroke length of the cylinder

2. Total Volume of the cylinder = Swept Volume + Clearance Volume = `V_s + V_c`

3. Compression ratio `k` = Total Volume/Clearance Volume = `(V_s+V_c)/(V_c) = 1+ (V_s/V_c)`

where,

`V_s` = Swept Volume of the cylinder`

`V_c` = Clearance Volume of the cylinder

4. Clearance Volume `V_c` = `V_s/(k-1)`

5. Thermal Efficiency `eta = 1-(1/(k^(gamma-1)))`

Thermodynamic relations for Otto Cycle:

• At point 2,

`p_1*v_1^gamma = p_2 * v_2^gamma`

`(p_1*v_1)/t_1 = (p_2*v_2)/t_2`

• At point 3,

`v_3 = v_2`

`p_3 / t_3 = p_2/t_2`

• At point 4,

`v_4 = v_1`

`p_4*v_4^gamma = p_3 * v_3^gamma`

where,

`p` , `v` and `t` indicate Pressure, Volume and Temperature at the respective states.

• To estimate the Piston kinematics,

`V / V_c = 1 + 1/2 * (k -1) [R + 1 - cos theta - (R^2 - sin^2 theta)^(1/2)]`

where,

`V` = Volume under expansion or compression

`V_c` = Clearance Volume

`k` = Compression Ratio

`R` = Connecting rod length / `a`

`a` = Crank pin radius = Stroke length / 2 = `l /(2)`

MATLAB Code:

For Otto cycle :

%Air Standard Cycle
clear all 
close all
clc


%Initial Conditions
%Pressure in Pa
%Temperature in K
p1 = 100000;
t1 = 300;

gamma = 1.4;
t3 = 2500;

%Engine parameters
bore = 0.2;
stroke = 0.25;

%Compression ratio 
k = 9;


%Connecting rod length
con_rod_len = 0.2;

%Swept volume and Compressed volume
v_swept = (pi/4)*bore^2*stroke;
v_clearance = v_swept/(k-1);

v1 = v_swept + v_clearance;
v2 = v_clearance;

%At point 2
%p1v1^gamma = p2v2^gamma
p2 = p1*k^gamma;

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

constant = p1*v1.^gamma;
v_comp = Engine_kinematics(bore,stroke,k,con_rod_len,180,0);
p_comp = constant./v_comp.^gamma;

%At point 3
v3 = v2;
p3 = p2*t3/t2;

constant = p3*v3.^gamma;
v_exp = Engine_kinematics(bore,stroke,k,con_rod_len,0,180);
p_exp = constant./v_exp.^gamma;

%At point 4
v4 = v1;
p4 = p3*(v3/v4)^gamma;

%Thermal Efficiency
eff = 1 - (1/k^(gamma-1))

figure(1)
hold on
plot(v1,p1,'*','color','r');
plot(v2,p2,'*','color','r');
plot(v3,p3,'*','color','r');
plot(v4,p4,'*','color','r');
plot(v_comp,p_comp,'color','b');
plot(v_exp,p_exp,'color','b');
plot([v2 v3],[p2 p3],'color','b');
plot([v4 v1],[p4 p1],'color','b');

xlabel('Volume');
ylabel('Pressure');

For engine_kinematics function:

%Engine Kinematics function

function [V_change] = engine_kinematics(bore,stroke,k,con_rod_len,start_angle,end_angle)

%Crank pin radius
a = stroke/2;
R = con_rod_len/a;

%Swept volume and Compressed volume
v_swept = (pi/4)*bore^2*stroke;
v_clearance = v_swept/(k-1);

theta = linspace(start_angle,end_angle,100);

t_1 = 0.5*(k-1);
t_2 = R+1 - cosd(theta);
t_3 = (R^2 - sind(theta).^2).^0.5;

V_change = (1 + t_1*(t_2-t_3)).*v_clearance;

end

Explanation:

The above image is the Workspace after running the program for Otto cycle and the values obtained are displayed.

• The engine parameters such as bore, stroke, connecting rod length, and compression ratio are given as inputs.
• The pressure and temperature at point 1 is assumed as 1 bar (1,00,000 Pa) and 300K respectively. The temperature at point 3 is assumed as 2500K. The value of the ratio of specific heat (gamma) is a constant and its value is 1.4
• The conventional formulae such as for swept volume, ideal gas equation and the relation between pressure and volume are used to solve the unknowns.
• A function called "engine_kinematics" is defined which is used to indicate how the volume changes inside the cylinder. Based on this, the curve varies during compression(process 1-2) and expansion(3-4)

Note: The above values are assumed and depending on the requirement and application, the values can be modified by the user to calculate Thermal efficiency and plot PV diagram, as the code for execution is the same as above.

The variables are explained below:

p - pressure (The numbers indicate the process)

t - temperature (The numbers indicate the process)

k - Compression ratio

con_rod_len - connecting rod length

v_swept - swept volume

v_clearance - clearance volume

v_comp - to calculate and store the values of volume during the compression process

p_comp - to calculate and store the values of pressure during the compression process

v_exp - to calculate and store the values of volume during the expansion process

p_exp - to calculate and store the values of pressure during the expansion process

eff - Thermal efficiency (not in percentage)

Results:

PV Diagram :

The PV diagram obtained for the values used in the code is as above. The red stars indicate the values at each point and the blue curve is the path traced according to the values during each process. As discussed in the introduction part, the volume at points 2 and 3 are the same and volume at points 1 and 4 are the same. The curves during processes 1-2 and 3-4 are obtained due to the engine_kinematics function that defines the variation of volume.

Thermal Efficiency:

The thermal efficiency was calculated as 0.5848 , that is 58.48% (indicated as "eff" in the Workspace and Command Window)

Conclusion:

The MATLAB program to solve an ideal Otto cycle and to calculate Thermal efficiency and plot PV Diagram was created and the explanation was given with the help of screenshots from application window.