PV Diagram of Otto cycle and Efficiency of Otto cycle.

AIM: To write MATLAB code that can solve an Otto Cycle, make plots and find efficiency of that.

Otto Cycle:

Otto cycle is description of thermodynamic cycle which is based on the function of spark ignition piston engine. Name of this cycle is dedicated to the Nicolaus Otto, German Engineer who build the first working four stroke engine. Most of the automobile engines follows this cycle for process. There are 4 process in the cycle.

Process 1-2: Isentropic compression

Process 2-3: Constant volume heat addition

Process 3-4: Isentropic expansion

Process 4-1: Constant volume heat rejection

To generate above mentioned P-V diagram, we need state points, that can be achieved at different states. This state points are basically P(Pressure), V(Volume) and T(Temperature). There are total 4 stages for that first we need to find the value of above-mentioned state points.

Step-by-step explanation of MATLAB code to obtain all state points:

1. Firstly, we need to get basic inputs and value of γ (gamma) from command window. These inputs are P₁, T₁, and T₃. Where P₁: Initial Pressure, T₁: Initial Temperature, T₃: Temperature after heat addition and γ = 1.4 constant.
2. Secondly, we need some engine parameters, which can be get from command window, too. These parameters are Bore (Engine bore diameter), Stroke (Stroke length), Connecting rod (Length of connecting rod) and r (Engine Compression Ratio).
3. We need to obtain value of V_c (Clearance Volume) and V₅ (Swept Volume) using engine parameters.
4. With the addition of V_c and V_s, we can get V₁ (Initial Volume), we already have other inputs (P₁ and T₁) for State 1.
5. For values of state points of State 2,
   1. V₂ = V_c (Clearance Volume)
   2. P₂ = by using `P_1 V_1^γ=P_2 V_2^γ`
   3. T₂ = by using Ideal gas equation`(P_1 V_1)/T_1 =(P_2 V_2)/T_2`
6. For values of state points of State 3,
   1. V₃ = V2 (Constant Volume)
   2. P₃ = by using Ideal gas equation
   3. T₃= already taken from Command Window`(P_2 V_2)/T_2 =(P_3 V_3)/T_3`
7. For values of state points of State 4,
   1. V₄ = V1 (Get back to initial volume)
   2. P₄ = by using `P_3 V_3^γ=P_4 V_4^γ`
8. To find the efficiency, we use following equation:

Efficiency = `1-(1/r^((γ-1)) )` where r = Engine Compression Ratio

γ = Gamma (is 1.4 in our case)

9.At the end, we have all state points values at States and Efficiency.

With the Help of state points, now we can plot all values on graph. Procedure for that plotting is as following:

1. After plotting all four point on graph, we can’t just connect them with plot command, otherwise we can-not get output as required. So, we have to create the function to calculate changes occurs in V during the compression and expansion process.
2. We can divide connection of lines by each process:
   1. 1 – 2: Isentropic Compression Process
   2. 2 – 3: Constant Volume Heat Addition
   3. 3 – 4: Isentropic Expansion Process
   4. 4 – 1: Constant Volume Heat Rejection
3. By using Engine Kinematics Function V,

`V/V_c = (1+1/2 (r-1)[R+1-cos⁡〖θ- (R^2-〖sin⁡θ〗^2 )^(1⁄2) 〗 ])`

Where, V indicates the rate of change in Volume from one state to other state

This function can use to determine How value of V changes in compression and expansion Process.

4. Now, For the compression process we calculate the V_compression using this function. For isentropic process `PV^γ=Constant`. So, we calculate first the value of constant with the value of P₁ and V₁.V_compression indicates how the volume changes during compression.
5. Then we determine the value of P_compression according to V_compression.
6. For the expansion process, we calculate the V_expansion and P_expansion same as V_compression and P_compression.
7. Now, we simply connect the state point 2 and 3 with plot command. Also connect the state point 4 and 1.
8. So, we obtain complete P-V diagram.

The errors faced while programming & steps taken to overcome those errors:

1. In the line 32, I forgot the use of the element wise operator for the constant_c. So, it is not calculating the p_compression.

Updates: I put the elementwise operator(.) .

2. In the line 23 I put the wrong formula for V₁ so the graph for the compression is wrong.

Updates: I correct the formula of V₁.

Screenshots showing errors thrown in command window

1.

2.

Program for the plot the Otto Cycle

%Write a program  to generate the PV diagram & efficiency of otto cycle
clear all 
close all
clc

%input
p1=input('input the initial pressure in Pa=')
T1=input('input the initial temperature in K=')
T3=input('input the temperature after heat addition in K=')
gamma=1.4

%engine parameters
bore=input('input the engine bore diameter in meter=')
stroke=input('input the engine stroke length in meter=')
con_rod=input('input the engine connecting rod length in meter=')
cr = input('input the engine compression ratio=')

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

%state variable at state 1
v1 = v_s + v_c
%state variable at state 2
v2 = v_c
p2=p1*cr^gamma          %according to p1*v1^gamma=p2*v2^gamma
T2 = p2*v2*T1/(p1*v1)         %p1*v1/T1=p2*v2/T2 (ideal gas eqation)

constant_c = p1*v1^gamma;
v_compression = engine_kinematics(bore,stroke,con_rod,cr,180,0);
p_compression = constant_c./v_compression.^gamma;

%state variable at state 3
v3=v2                   %constant volume
p3=p2*T3/T2           %p2*v2/T2=p3*v3/T3 (ideal gas eqation)

constant_c = p3*v3^gamma
v_expansion = engine_kinematics(bore,stroke,con_rod,cr,0,180);
p_expansion = constant_c./v_expansion.^gamma;

%state variable at state 4
v4=v1
p4=p3*(v3/v4)^gamma      %according to p3*v3^gamma=p4*v4^gamma

%efficiency of otto cycle
efficiency=1-(1/cr^(gamma-1))

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

Program of subfunction

function [V]= engine_kinematics(bore,stroke,con_rod,cr,start_angle,end_angle)

a=stroke/2
R=con_rod/a

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

theta=linspace(start_angle,end_angle,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

Input of the program in command window

input the initial pressure in Pa=95000

p1 =

       95000

input the initial temperature in K=288

T1 =

   288

input the temperature after heat addition in K=1473

T3 =

        1473


gamma =

    1.4000

input the engine bore diameter in meter=0.1

bore =

    0.1000

input the engine stroke length in meter=0.15

stroke =

    0.1500

input the engine connecting rod length in meter=0.2

con_rod =

    0.2000

input the engine compression ratio=8

cr =

     8

Output of the above program

1.Efficiency of otto cycle

2.PV diagram