Otto cycle using matlab

Aim:

To write a matlab code for solving otto cycle and to generate plot of p-v diagram.

Objective:

To write a code for solving otto cycle using matlab, also to calculate the thermal efficiency of engine and to plot PV diagram of otto cycle.

Introduction:

The Otto cycle is the ideal cycle for spark-ignition reciprocating engines. This cycle is also known as constant volume cycle as the heat is received and rejected at constant volume. The cycle consists of two adiabatic processes and two constant volume processes as shown in P-v and T-s diagrams.

Process 1-2:

In this air is compressed isentropically from V₁ to V₂ Temperature increases from T₁ to T₂. Since this is an adiabatic process heat rejected is zero

Process 2-3:

In this air is heated at constant volume and temperature increases from T₂ to T₃^.

Process 3-4:

In this air is expanded isentropically from V₃ to V₄ and temperature decreases from T₃ to T₄. Since this is an adiabatic process, the heat supplied is zero.

Process 4-1:

In this air is cooled at constant volume and temperature decreases from T4 to T1.

Efficiency of the cycle is given by,

`eta=1+frac{t1}{t2}`

Geometric parameters engine are:

Bore =0.2m

Stroke =0.1m

Connecting rod length = 0.15m

Compression ratio =7

Equations used are:

Swept volume =(pi/4)*bore^2*stroke;

Clearance volume =swept volume/(cr-1);

Volume trace is given by,

`frac{V1}{V2}=1+frac{1}{2}(cr-1)[R+1-costheta-R^2-(sin^2theta)^frac{1}{2}]`

Explanation of code:

Here in this matlab program, Otto cycle has been solved and thermal efficiency is found, also a PV diagram is plotted. In the beginning inputs and geometric parameters of engine are provided (i.e. values of gamma, t3, bore, stroke, connecting rod length, compression ratio). Then swept volume and clearance volume are calculated using formulas. The state variable at state1 (i.e. pressure, volume, temperature) are given and state variables at state2 is calculated by using ideal gas law, PV=nRT,

P2 is calculated using formula, p2 = p1*cr^gamma,

T2 is calculated using formula, t2 = p2*v2*t1/(p1*v1),

To calculate constant C ideal gas equation .

Constant c is calculated using formula, constant_c = p1*v1^gamma.

Expansion and compression curves are made by using engine kinematics function i.e. a user defined function. In this function a equation is solved using a formula, formula is given by,

`frac{V1}{V2}=1+frac{1}{2}(cr-1)[R+1-costheta-R^2-(sin^2theta)^frac{1}{2}]`

In this engine parameters are taken as inputs. This function will be used to calculate crank pin dimensions and then a range of 100 values of theta from starting to ending crank angles is crated using linspace command and then implementing formula by dividing the formula into 3 terms.

Then by calling the function engine kinematics compression volume and compression pressure are calculated. Similarly volume and pressure at expansion time are calculated. Then state variables at state 3 and state 4 are calculated. Then efficiency of engine is calculated. After that plot for PV diagram will be plotted.

code for otto cycle:

clear all
close all 
clc

%inputs
gamma = 1.4;
t3 = 1273;

%state variables
p1 =101325;
t1 = 300;
%engine geometric parameter
bore = 0.2
stroke = 0.1
con_rod = 0.15
cr = 7

%calculating the swept volume & the clearence 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;

%p1v1/t1=p2v2/t2|t2=p2*v2*t1/(p1*v1)
t2 = p2*v2*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;
%p3v3/t3=p2v2/t2|p3=p2*t3/t2
p3 = p2*t3/t2;
constant_c = p3*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;
%p3v3^gamma=p4v4^gamma|p4=p3*v3/v4^gamma
p4 = p3*v3^gamma/v4^gamma;

%efficiency of cycle
E = 1-(t1/t2)

figure(1)
hold on 
plot(v1,p1,'*','color','r')
plot(v_compression,p_compression,'color','m')
plot(v2,p2,'*','color','r')
plot([v2 v3],[p2 p3],'color','m')
plot(v3,p3,'*','color','r')
plot(v_expansion,p_expansion,'color','m')
plot(v4,p4,'*','color','r')
plot([v4 v1],[p4 p1],'color','m')
xlabel('volume')
ylabel('pressure')
grid on

code for engine kinematics function:

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

plot of PV diagram is as below :

Errors:

Here the errors was due to logical values, in the line 48 formula was wrongly written as p4=p3(v3/v4)^gamma which was wrong. Then it is corrected to right form which is p4=p3*v3/v4^gamma.

Conclusion:

Thermal efficiency is determined as 54.08% by the calculation.

A plot of PV diagram is created.