PV Diagram and thermal efficiency of the Otto Cycle using MATLAB

PV Diagram and thermal efficiency of the Otto Cycle using MATLAB

Objective:

To write code in MATLAB to plot the PV Diagram for Otto cycle and calculating the thermal efficiency of the same.

Theory:

Otto cycle is the thermodynamic cycle which forms the basis for the working for spark-ignition engines.

The process involved in the cycle is described in the below figures.

IC Engine Nomenclature:

Schematic representation of engine geometry parameters

Inputs:

1. Temperature(T1) and pressure (P2) at the beginning of the compression stroke. This usually corresponds to the ambient temperature and pressure.
2. Temperature(T3) at the end of the heat addition/ignition process. This temperature is usually limited by the metallurgical limits of the material used.
3. Specific heat ratio(gamma)
4. Engine geometric parameters like bore diameter, stroke length, connecting rod length and compression ratio.

Code:

%Program for plotting PV diagram and finding thermal efficiency of the otto cycle
clear all
close all 
clc

%Inputs
gamma = 1.4
T3 = 2300

%State Variables
P1 = 101325
T1 = 298

%Engine Geometric Parameters
bore = 0.1;
stroke = 0.1;
con_rod = 0.15
cr = 12;

%calculating 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 = P1*cr^gamma
T2 = P2*T1/(P1*V1)

constant_c = P1*V1^gamma;
V_compression = engine_kinematics(bore,stroke,con_rod,cr,180,0);
P_compression = constant_c./V_compression.^gamma;

%State variables at state point 3
V3 = V2
P3 = P2*T3/T2

constant_c = P3*V3^gamma;
V_expansion = engine_kinematics(bore,stroke,con_rod,cr,0,180);
P_expansion = constant_c./V_expansion.^gamma;

%State variables at state point 4
V4 = V1
P4 = P3*(V3/V4)^gamma
T4 = P4*T1/P1

figure(1)
hold on
plot(V1,P1,'.','color','r');
plot(V_compression,P_compression,'color','r');
plot(V2,P2,'.','color','r');
plot([V2 v3],[P2 P3],'.','color','r');
plot(V3,P3,'.','color','r');
plot(V_expansion,P_expansion,'color','r');
plot(V4,P4,'.','color','r');
plot([V1 v4],[P1 P4],'.','color','r');
xlabel('volume');
label('Pressure');
grid on

%Thermal efficiency
eta = 1 - (1/cr^(gamma-1));
fprintf('Thermal efficiency of the cycle is: %f',eta*100)

Explanation:

The initial inputs gamma, P1, T1, T3 are initialized. Units for the quantities are

Pressure - Pascal, Temperature - Kelvin, Volume - cubic meters, Other geometric parameters - meters.

The swept/stroke volume and clearance are found using the engine geometric parameters

Volumes V1 and V2 are assigned.

The pressure at state point 2 can be calculated using the isentropic relation

The temperature at state point 2 can be calculated by

Finding the state variables at intermediate points between 1 & 2

Since the volume change from 1 to 2 is not a linear one, we need to calculate the pressure and volume at intermediate points for plotting the graph. This can be done using the following equation

This equation is based on the engine parameters and gives the cylinder volume during the compression as well as the expansion process.

R = Connecting rod length/ crank radius

θ = crank angle

rc = compression ratio

The Matlab code calls for the function engine_kinematics for finding out the intermediate cylinder volumes which take input arguments in the order bore, stroke, connecting rod length, compression ratio, starting crank angle and ending crank angle. The starting and ending crank angles for the compression process are 180 deg and 0 deg respectively. The function returns the intermediate volume values as an array. Pressure values at the intermediate points are found out using the isentropic relation.

Code for the function engine_kinematics.m

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

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

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

The variables at point 3 V3 & P3 are found out using the same engine_kinematics function the difference being the starting and ending crank angles are 0 deg and 180 deg respectively. The state variables at point 4 are also found out.

The graph is plotted using the plot command.

The thermal efficiency of the cycle is found out using the formula

Results:

Output

gamma = 

     1.4000


T3 = 

     2300


P1 = 

    101325


T1 = 

   298


V1 = 

    8.5680e-004


V2 = 

   7.1400e-0.005
            

P2 = 

   3.2853e+006


T2 = 
  
   805.1722


V3 = 

   7.1400e-005


P3 = 

   8.5680e-004


P4 = 

   2.8944e+005


T4 = 

   851.2465


The thermal efficiency of the cycle is: 62.989283>>

The PV diagram for the otto cycle is plotted.

The thermal efficiency of the cycle = 62.9%