P-V Diagram for Otto cycle using MATLAB.

Aim:

Writing a MATLAB code to plot Otto cycle.

Objectives:

• To plot P-V diagram for Otto cycle.
• To calculate thermal efficiency of cycle.

Theory:

Assumptions:

• Cycle is reversible.
• Working fluid is air.
• Air behaves as a perfect gas.
• Specific heat of working fluid remains constant.

Stroke volume is the displacement volume of the piston.

Clearance volume is a volume between the cylinder head and the piston top when the piston is at top dead centre (TDC).

These are given blow,

`Vs = pi/4xx(D)^2xxLs` where, D = bore; Ls = stroke.

`Vc = (Vs)/(cr-1)`

`V1 = Vs+Vc` ; `V2 = Vc`

1-2 is isentropic compression process for which,

`P2 = P1*(cr)^gamma`

using ideal gas equation,

`T2 = (P2*V2*T1)/(P1*V1)`

at point3,

`V3 = V2` ; `P3 = (P2*T3)/(T2)`

similarly, at point4

`V4 = V1` ; `P4 = (P3)/(cr)^gamma`

Thermal efficiency of cycle is,

`eta = (1 - (T1)/(T2))`

Body of Content:

The programming is done in MATLAB. Specific heat ratio is taken as 1.4, Temperature at poin3 is also taken as input. In addition to calculating state variables at four points as given in PV diagram, we need to plot variation of pressure and volume and for compression and expansion process volume is a function of crank angle. So, we need to calculate volume for different values of crank angle, to care of this we first defined the function which takes bore, stroke, connecting rod length, compression ratio, start and end crank angle as inputs. Then we defined various term to calculate variation in volume according to crank angle.

function [V] = ch4_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; %this V and V in the function is same.

end

In main program the inputs are defined as given. From which we calculated stroke volume and clearance volume first, it forms the basis to calculate state variables. We calculated the state variables at point 2,3&4 using formulas given above. To get volume and pressure array for different values of crank angle for process 1-2 and 3-4 we incorporated the function. Then finally we plotted the state variables and process to form a closed cycle 1-2-3-4.

clear all 
close all 
clc

gamma = 1.4;
t3 = 2300;  %temperature(in K) at point3

%State Variables
p1 = 101325;    %Pressure(Pa) at point1
t1 = 500;       %temperature(in K) at point1

%Engine geometric parameters
bore = 0.08;    %Diameter(m) of piston
stroke = 0.1;   %lenght(m) of travel of piston
con_rod = 0.06; %connecting rod lenght(m)
cr = 8;        %compression ratio

%swept and clearance volume
v_s = (pi/4)*bore^2*stroke; %swept volume
v_c = v_s/(cr-1);           %clearance volume

v1 = v_s + v_c              %Volume(m^3)at point1.
v2 = v_c

%State variable at point2
p2 = p1*(cr)^gamma;
t2 = p2*v2*t1/(p1*v1)

constant_c = p1*v1^gamma;
V_compression = ch4_engine_kinematics(bore, stroke, con_rod, cr, 180, 0); %to get volume array on compression curve using array of theta.
P_compression = constant_c./(V_compression).^gamma;

%State variable at point3
v3 = v2;
p3 = p2*t3/t2

constant_c = p3*v3^gamma;
V_expansion = ch4_engine_kinematics(bore, stroke, con_rod, cr, 0, 180); 
P_expansion = constant_c./(V_expansion).^gamma;                     %to get pressure array on expansion curve using array of theta.


%State variable at point4
v4 = v1;
p4 = p3/(cr)^gamma

%Thermal Efficiency
effi = (1 - t1/t2)

%plotting
figure(1)   
hold on     %to execute every plot command below one after another without suppressing any.
plot(v1,p1, '*','color','r')    %plotting state1.
plot(V_compression, P_compression,'color', 'k') %plotting process 1-2
plot(v2,p2, '*','color','r')
plot(v3,p3, '*','color','r')
plot(V_expansion, P_expansion,'color','k')
plot(v4,p4, '*','color','r')
plot([v2 v3],[p2 p3],'color','k')   %plotting process 2-3
plot([v4 v1],[p4 p1], 'color', 'k')
xlabel('Volume', 'color', 'b')
ylabel('Pressure', 'color', 'b')

Final Output:

PV plot for Otto cycle

Output in the command window:

Errors:

I forgot to include element wise operator in pressure array.