Otto Cycle Plot in MATLAB

Title - Otto Cycle Plot in MATLAB

Objectives - 1. To plot PV diagram of otto cycle in MATLAB

2. To determine thermal efficiency of cycle in MATLAB

3. To learn application of MATLAB for PV plots

Function Code -

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

a = stroke/2;
R = con_rod/a;

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

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

end

Main Code -

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 = 6;

% calculating the swept volume and the 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
% p2v2^gamma = p1v1^gamma
% p2 = p1*(v1/v2)^gamma = p1*cr^gamma
p2 = p1*cr^gamma;

% 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 state point 3
v3 = v2;

% p3*v3/t3 = p2*v2/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 state point 4
v4 = v1;
p4 = p3*(v3/v4)^gamma;

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([v1 v4],[p1 p4])
plot([v2 v3],[p2 p3])

% Efficiency of otto engine

thermal_efficiency = (1 - (1/(cr^(gamma-1))))*100

Output -

Command Window -

Graph -

Code Explanation -

Function for volume at different stages of isentropic expansion and contraction is made. The formula for this is v/clearance volume = 1 + 0.5*(cr - 1)(R + 1 - Cos(θ) - (R^2 - (Sinθ)^2)^0.5), which is derived mathematically. So, the function evaluated the volume by intaking different values required in the function. This function will be used in the main code.

In the main code, different stages of otto cycle is computed. Given, is gamma of air = 1.4, compression ratio of stroke, stroke bore, length of connecting rods, diameter of piston. Also, given are p1, t1 and t3. swept volume = (pi/4)*bore^2*stroke and clearance volume = swept volume (compression ratio - 1)

First cycle is isentropic contraction, hence pv^gamma = constant. from this we will calculate, p2 which is p1*cr^gamma, since v1/v2 = cr. From ideal gas equation, we can get t2 since we know t1, p1, v1, p2 and v2. Now, Array of volume in this stage will be calculated by using volume function.

Second cycle, is isochoric heat addition. Thus v3 = v2. We know t2, t3 and p2, thus we can calculate p3 from ideal gas equation for isochoric process

Third cycle is isentropic expansion. v4 = v1 (since last cycle is isochoric). We know v3, t3, p3, v4, now we can calculate p4 from isentropic equation p4 = p3*(v3/v4)^gamma.

Now, we have all state variables and interval values for isentropic functions. We can now plot all state variables. For isentropic process, we will join the array of values returned from the function. For isochoric processes, we will join with straight lines.

Errors encountered - Following error was encountered

" Undefined function or variable 'con_rod'.

Error in Otto_Cycle_error (line 35)
V_compression = engine_kinematics(bore,stroke,con_rod,cr,180,0); "

Error Screenshot - Following is the error screenshot

Steps to overcome error -

The error ocurred due to mismatch of variable name defined at the top and in the function argument. At the top of the code it was "connecting_rod", here it is con_rod. After the correction, code ran smoothly.

Thanks!!