Otto cycle simulator using MATLAB

AIM: To visualize and solve an Otto cycle using MATLAB.

GOVERNING EQUATION: The equations to solve an Otto cycle are as follows:

`PV^gamma` = c

`(PV)/T = c`

`V_c = frac{pi}{4} .D^2 .L`

`V_s = frac{V_c}{r_c - 1}`

`r_{c} = frac{V_2}{V_1} = frac{V_3}{V_4}`

`V_theta = [1 + frac{1}{2}(r_c - 1) sqrt(1 + R - costheta - sqrt(R^2 - sintheta^2))].V_c`

`eta_{o} = 1 - frac{1}{(r_c)^(gamma- 1)) = 1 - frac{T_4 - T_1]{T_3 - T_2}`

`eta_{c} = 1 - frac{T_min}{T_max} = 1 - frac{T_1}{T_3}`

where;

`P` = pressure (Pa)

`V` = volume (m³)

`T` = temperature (K)

`c` = constant

`gamma` = isentropic/adiabatic index

`V_c` = clearance volume (m³)

`V_s` = swept volume (m³)

`theta` = crank angle (degree)

`eta_o` = otto cycle efficiency

`eta_c` = carnot cycle efficiency

OBJECTIVE:

1. To plot the P-V diagram of an Otto cycle.

2. To calculate the thermal efficiency

3. To study piston kinematics by using of an Otto cycle. as a function in the isentropic processes.

THEORY:

An Otto cycle is an idealized thermodynamic cycle that describes the functioning of a typical spark ignition piston engine. It is the thermodynamic cycle most commonly found in automobile engines. The processes involved during this cycle are as follows; Process 1–2: Adiabatic (isentropic) compression of the charge as the piston moves from the bottom dead center (BDC) to top dead center (TDC) Process 2–3: Constant-volume heat transfer to the working gas from an external source while the piston is at top dead center. This process is intended to represent the ignition of the fuel-air mixture and the subsequent rapid burning Process 3–4: Adiabatic (isentropic) expansion (power stroke) Process 4–1: Completes the cycle by a constant-volume process in which heat is rejected from the air while the piston is at bottom dead center

Piston Kinematics: Piston kinematics is the study of the motion of the piston relative to the rotation of the crankshaft. A reciprocating engine derives its power from the vertical motion of a piston within individual cylinders of the engine. The up-down motion of the piston, also known as reciprocating motion, is converted into rotary motion by use of a crankshaft.

SOLUTION:

Programming Language: MATLAB

Procedure/Steps:

1. Assign suitable variables to store input parameters for the links.

2. Calculate swept volume and clearance volume with the help of given equations.

3. Using gas laws and other suitable process relations, determine the pressure, temperature and volume at all state points.

4. An array is incorporated to find V_comp, P_comp and V_exp, P_exp.

5. At this point, we have calculated all state points. So Otto and Carnot efficiencies are calculated.

6. Finally, all the state points are plotted.

Code: The MATLAB program to plot the P-V diagram is shown below.

% Program to plot P-V Diagram and calculate efficiency of an Otto cycle.

% Inputs:
% (NOTE: P = pressure (Pa), V = volume (m^3), T = temperature (K))

gamma = 1.4;   % isentropic index
cr = 10;       % compression ratio

% Engine Geometric Parameters
bore = 0.1;      % bore radius (m)
stroke = 0.1;    % stroke lenght (m)
con_rod = 0.15;  % connecting rod length (m)

% State 1
P1 = 101325;
T1 = 300;

% State 3
T3 = 1600;


% Solution:

% State 1
V_s = (pi/4)*(bore^2)*stroke;  % swept volume
V_c = V_s/(cr-1);              % clearance volume
V1 = V_s + V_c;

% State 2
V2 = V_c;
% P1*V1^gamma = P2*V2^gamma | P2 = P1*(V1/V2)^gamma
P2 = P1*(cr^gamma);
% P1*V1/T1 = P2*V2/T2
T2 = (P2*V2*T1)/(P1*V1);

% Parametrs during Isentropic Compression
c1 = P1*V1^gamma;
V_comp = otto_cycle_kinematics(bore,stroke,con_rod,cr,0,180);
P_comp = c1./(V_comp.^gamma);


% State 3
V3 = V2;
% P2*V2/T2 = P3*V3/T3
P3 = (P2*T3)/T2;

% Parametrs during Isentropic Expansion
c3 = P3*V3^gamma;
V_exp = otto_cycle_kinematics(bore,stroke,con_rod,cr,180,0);
P_exp = c3./(V_exp.^gamma);

% State 4
V4 = V1;
% P3*V3^gamma = P4*V4^gamma | P4 = P3*(V3/V4)^gamma
P4 = P3*(V3/V4)^gamma;
% P3*V3/T3 = P4*V4/T4
T4 = (P4*V4*T3)/(P3*V3);

% Efficiency
otto_eff = (1 - ((T4-T1)/(T3-T2)))*100
carnot_eff = (1 - (T1/T3))*100
txt = ['Efficiency = ', num2str(otto_eff) , '%'];


% Plot:

figure(1)
hold on
plot(V_comp,P_comp,'color','m')
plot([V2 V3],[P2 P3],'color','g')
plot(V_exp,P_exp,'color','b')
plot([V1 V4],[P1 P4],'color','k')
plot(V1,P1,'x','color','r')
plot(V2,P2,'x','color','r')
plot(V3,P3,'x','color','r')
plot(V4,P4,'x','color','r')
text(V1,P1,'1','HorizontalAlignment','left','VerticalAlignment','top')
text(V2,P2,'2','HorizontalAlignment','right','VerticalAlignment','top')
text(V3,P3,'3','HorizontalAlignment','right','VerticalAlignment','bottom')
text(V4,P4,'4','HorizontalAlignment','left','VerticalAlignment','bottom')
text(0.5*(V1+V2),0.5*(P1+P3),txt)
xlabel('Volume (m^3)')
ylabel('Pressure (Pa)')
title('P-V diagram of Otto Cycle')
legend('Isentropic Compression','Constant Volume Heat Addition','Isentropic Expansion','Constant Volume Heat Rejection')
grid on

Output: The efficiencies are calculated and the problem is solved as shown below.

Plot: All four state points are plotted and connected according to their respective process nature as shown below.

Function: The MATLAB program as a volume function [V] is shown below. We call the function otto_cycle_kinematics to find the volume trace as a function of crank angle.

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

% Inputs:

a = stroke/2;   % crank pin radius
R = con_rod/a;  % crank ratio


% Solution:

V_s = (pi/4)*(bore^2)*stroke;  % swept volume
V_c = V_s/(cr-1);              % clearance volume

th = linspace(start_crank, end_crank, 100);  % crank angle

x1 = (cr-1)/2;
x2 = 1 + R - cosd(th);
x3 = (R^2 - sind(th).^2).^0.5;

V = (1 + x1*(x2-x3)).*V_c;   % volume in terms of crank angle


% Plot:

plot(V)
xlabel('Crank angle (deg)')
ylabel('Volume (m^3)')
grid on

end

Plot: The nature of volume with respect to crank angle rotation is shown below.

Error: Problem in obtaining an accurate curve due to merging the state points in a single command.

Solution: A function was written which included crank angle and volume during the compression and expansion process. Also, separate plot commands were given for each state point so that they are displayed according to their process nature.

CONCLUSION: From the above results, we can infer that the Otto cycle efficiency is less as compared to Carnot cycle efficiency.

`eta_o < eta_c`

Hence, our calculations are correct and the Otto cycle is successfully plotted.

REFERENCES:

1. Otto Cycle

2. Piston Kinematics