MATLAB code to solve and generate plot for Otto cycle

OBJECTIVE:

To write a Matlab code to solve and generate plot for Otto cycle. The code should create a PV diagram of the cycle and should also calculate the thermal efficiency of the engine.

SOLUTION:

A Matlab code was written to solve the Otto cycle, to find the pressure, volume and temperature values at the various state points. Further, the code generates a plot of the Otto cycle with respect to the input variables given and also calculates the thermal efficiency of the engine.

Explanation:

The aim of the program is to solve and generate the Otto cycle. In order to perform this task, it is necessary to solve for the pressure, volume and temperature values at the different state points. A theoretical Otto cycle has four state points, with four process paths; isentropic compression, isochoric compression, isentropic expansion and isochoric expansion. In order to solve the cycle, we need to have knowledge of some of the variables. Hence, some variables were defined at the initial stage including p1, t1, t3, γ, bore diameter, stroke length, compression ratio and connecting rod length.

The swept volume and clearance volume was calculated using these variables. As per the cycle, the volume at state point 1, v1 is equal to the sum of swept volume and clearance volume and v2 is equal to clearance volume. Since 2-3 and 4-1 are isochoric processes, v2 = v3 and v4 = v1. All the other state points were calculated using the below relations:

(PV^γ)_a = (PV^γ)_a --- [A]

(PV/T)_a = (PV/T)_b --- [B]

After finding the variables at the four state points, the isentropic compression and expansion paths was simulated. For this, the below relation is used:

where,

r_c - compression ratio

R - ratio of connecting rod length to crank pin radius

Crank pin radius = a = stroke length/2

To simulate the compression and expansion paths, a function named engine kinematics was defined. Inside this function an array of crank angle was introduced and the volume values at each of angle was calculated. In the main program, the pressure values corresponding to volume was found.

The efficiency of the cycle is calculated using the relation W_net/Q_in.

W_net/Q_in = 1 - ( t4 - t1 / t3 - t2 )

Outputs generated:

1. All values found at the state points was printed as a table

2. The otto cycle plot was generated

3. The efficiency of the cycle was calculated and printed

%program to simulate Otto cycle

clear all
close all
clc

%definition of variables

p1 = 101325;
t1 = 500;
gamma = 1.4;
t3 = 2300;

%calculation of intial volume v1

bd = 0.1;
stroke = 0.1;
cr = 12;
con_rod = 0.15;
swept_vol = (pi/4)* bd^2 * stroke;
clearence_vol = swept_vol/(cr-1);

v1 = swept_vol+clearence_vol;

%calculation of second state point values

v2 = clearence_vol;
p2 = p1*(v1/v2)^gamma;
t2 = (p2*v2*t1)/(p1*v1);

%calculation of v_compression

v_compression = engine_kinematics(clearence_vol,stroke,con_rod,cr,180,0);

%calculation of constant

c1 = p1*v1^gamma;

%calculation of p_compression

p_compression = c1./v_compression.^gamma;

%calculation of third state point values

v3 = v2;
p3 = (p2*t3)/t2;

%calculation of fourth state point values

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

v_expansion = engine_kinematics(clearence_vol,stroke,con_rod,cr,0,180);
c2 = p3*v3^gamma;
p_expansion = c2./v_expansion.^gamma;


hold on
plot(v1,p1,'*','color','b');
plot(v2,p2,'*','color','b');
plot(v3,p3,'*','color','b');
plot(v4,p4,'*','color','b');
plot(v_compression,p_compression, 'color', 'r');
plot(v_expansion,p_expansion, 'color', 'r');
plot([v2 v3], [p2 p3],'color','g');
plot([v1 v4], [p1 p4],'color','g');
xlabel('Volume');
ylabel('Pressure');

%calculation of thermal efficiency
eeta1 = 1-((t4-t1)/(t3-t2));

eeta2 = eeta1 *100;

%displaying state point values in the form of a table
Data = [ 1      p1     v1    t1
         2      p2     v2    t2
         3      p3     v3    t3
         4      p4     v4    t4];

ColumnNames = {'state_point', 'pressure (Pa)', 'volume (m^3)', 'temperature (K)'};

T = table(Data(:,1),Data(:,2),Data(:,3),Data(:,4),'VariableNames', ColumnNames);

%Printing final output of the program
disp(T);
fprintf('The thermal efficiency (in percentage) = %d',eeta2);

function [V] = engine_kinematics(clearence_vol,stroke,con_rod,cr,start_crank,end_crank)

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

theta = linspace(start_crank,end_crank,100);

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

V = (1+term1*(term2-term3))*clearence_vol;

end 

OUTPUT: