To obtain the P-V diagram of Otto Cycle using Matlab/Octave programming and calculate the thermal efficiency of the cycle.

AIM

To write a code in Matlab/Octave to obtain the P-V diagram of Otto Cycle and thereby calculate the thermal efficiency of the cycle.

OTTO CYCLE

An Internal Combustion Engine or an IC Engine is generally referred to as the heart of the automobile. It is the region where combustion of fuel occurs, which produces energy for motion of the vehicle. The combustion of fuel takes place with the help of an oxidizer, usually air, in a chamber, which is an important part of the fluid flow circuit. In an IC Engine, chemical energy is converted into mechanical energy. The combustion in an IC engine is intermittent, which means the combustion happens inside a closed environment.

There are mainly four strokes on which an engine generates power to run an automobile which is shown in diagram

1. Suction Stroke – It is the start of the Engine cycle. During this stroke, the inlet valve opens, and the piston which is at the top dead center (TDC) moves towards bottom dead center (BDC) as the air-fuel mixture enters the cylinder.
2. Compression Stroke – The piston now moves from the bottom dead center to the top dead center which increases the pressure of the cylinder due to compression of the air-fuel mixture.
3. Expansion/Power Stroke – This is the power generating stroke of the engine cycle. The combustion of fuel takes place due to spark generation, and the piston moves from the top dead center to bottom dead center due to the enormous power generated during sudden combustion of the fuel.
4. Exhaust Stroke – After power stroke, the piston again moves towards top dead center, followed by the opening of the exhaust valve, through which the gases produced during combustion escapes the cylinder thereby making space for fresh air-fuel mixture to enter and the cycle continues.

A typical petrol engine works on the Otto cycle, which is an ideal thermodynamic cycle which elaborates the working of a general spark ignition piston system. Most automobile engines work on the Otto cycle. The ideal representation of the Otto cycle is shown in diagram. The Otto cycle is a basic description of the changes occurring in a mass of gas when subjected to parameters such as temperature, pressure, volume, heat addition and heat removal and the system’s effect on the environment. The ideal Otto cycle consists of the following phases.

Process 0-1 – It is a constant pressure process where air is drawn into the piston.

Process 1-2 – Isentropic compression phase as piston moves from BDC to TDC.

Process 2-3 – Constant Volume Heat Transfer/Addition phase when the piston is at TDC.

Process 3-4 – Isentropic Expansion or Power Stroke

Process 4-1 – Constant Volume Heat Rejection Phase when the piston is at BDC.

Process 1-0 – Air released into atmosphere at Constant Pressure.

GOVERNING EQUATIONS

The equations are with respect to the otto cycle diagram given above and its various points.

`1)Vs= pi/4*B^2*S`

where,

Vs = Swept volume

B = Bore

S = Stroke length

`2) Vc= (Vs)/(Rc-1)`

where,

Vc = Clearance volume [ Volume present when the engine cylinder is at Top Dead Center (TDC)]

Rc = Compression ratio [ Vs/Vc ]

3) `PV^gamma=C`

where,

P = Pressure of gas

V = Volume of gas

C = Constant

AT POINTS 1 AND 2 OF THE CYCLE

`V1=Vs+Vc`

`V2=Vc`

From the equation `P1V1^gamma = P2V2^gamma`, `P2=P1*((V1)/(V2))^gamma`

From the equation, `(P1V1)/(T1)=(P2V2)/(T2)``, `T2=(P2*V2*T1)/(P1*V1)`

AT POINTS 3 AND 4 OF THE CYCLE

`V2=V3`

`V4=V1`

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

From the equation `P3V3^gamma = P4V4^gamma`, `P4=P3*((V3)/(V4))^gamma`

`(P3V3)/(T3)=(P4V4)/(T4)`,`T4=(P4*V4*T3)/(P3*V3)`

PISTON KINEMATICS EQUATION

`V/(Vc)=1+1/2(Rc-1)[R+1-costheta-(R^2 - sin^2theta)^(1/2)]`

where,

V = Volume variation with respect to crank angle

Vc = Clearance volume

R = `(Crank Radius)/(S/2)`

OBJECTIVES OF THE PROJECT

The major objective of this project is to implement functions in the Matlab/octave code with respect to piston kinematics equation and call the function when needed to plot the Otto cycle and calculate the thermal efficiency. The piston kinematics equation gives us the variation of volume with respect to the crank angle within the engine during compression and expansion. This calculation helps us get an isentropic curve for both compression and expansion in the Otto cycle.

EXPLANATION OF THE CODE

The following code is written in Octave

There are two codes in this project, one for the Otto cycle and the other one for obtaining the isentropic curves. Both codes shall be explained below.

1) Piston Kinematics Code

function[S]=piston1(bore,stroke,con_rod_length,r_c,start_crank_angle,end_crank_angle) – The function command is used in Octave to specify required values during the time of running the program, so that variable values can be tested in the running of the program instead of a standard set of values. Here, the various engine parameters like bore, stroke, length of connecting rod, compression ratio and the starting and ending crank angles are the parameters that can be specified in the command window. ‘piston1’ is the name of the saved file from which the function will be called. ‘S’ is a variable which stores the value of the volume obtained as a factor of changing crank angles.

theta = linspace(start_crank_angle,end_crank_angle,50); - The linspace command is used to divide the starting and ending values into equal number of parts as specified by the programmer. Here, the starting and ending crank angles are varied during the compression and expansion strokes, and 50 values are considered between them in this case.

crank_pin_radius = stroke/2;

R=r_c/crank_pin_radius;

v_swept = (pi/4)*bore*bore*stroke;

v_clear = v_swept/(r_c - 1); - The above statements are used to calculate the swept and clearance volumes to calculate the volume variation during compression and expansion strokes.

term_1 = 0.5*(r_c -1);

term_2 = R+1-cosd(theta);

term_3 = ((R^2)-(sind(theta).*sind(theta))).^0.5; - As mentioned above, the governing equation for the variation of volume is a very huge equation. The three terms bear three portions of the equation which are then merged to form the final equation

S = (1+(term_1*(term_2 - term_3)))*v_clear; - The merger of the three terms to form the final equation.

Endfunction – This command is used to close the function command

2) Otto Cycle code

p1=101325;

t1=600;

bore=0.15;

stroke=0.17;

con_rod_length=0.2;

r_c=12;

gamma=1.4;

t3=2500; - The above statements are used to input the constant values that are required to obtain the thermal efficiency of the Otto cycle. The value of gamma is constant at 1.4 for a polytropic process. The values of bore, stroke and connecting rod lengths are assumed. Similarly, the values of initial pressure and temperatures are also assumed.

vswept=(pi/4)*bore*bore*stroke

vclear=vswept/(r_c-1) – The above statements are used to calculate the swept and clearance volumes based on the input engine parameters.

v1=vswept+vclear

v2=vclear – These statements are used to calculate the values of v1 and v2 which are theoretically obtained ones.

p2=p1*(v1/v2).^gamma – The pressure at point 2 is calculated using standard theoretical formula.

C=p1*(v1^gamma); - This command is used to find the constant value ‘C’ in the equation “ PV^gamma=C” which in turn is used to graph the compression isentropic curve.

V_compression=piston1(bore,stroke,con_rod_length,r_c,180,0);

[P_compression]=C./(V_compression.^gamma); - These statements are used to call the values from the piston kinematics equation code. The value of volumes at different crank angles during compression stroke(V_compression) is first calculated and then the corresponding pressures(P_compression) are calculated to obtain an isentropic curve between points 1 and 2 of the Otto cycle for crank angles from 180 to 0 degrees.

t2=(p2*v2*t1)/(p1*v1)

v3=v2

p3=(p2*v2*t3)/(t2*v3) – These commands are used to find temperature, volume and pressure at point 3 respectively.

C1=p3*(v3^gamma); - This command is used to find the constant value ‘C’ in the equation “ PV^gamma=C” which in turn is used to graph the expansion isentropic curve.

V_expansion=piston1(bore,stroke,con_rod_length,r_c,0,180);

[P_expansion]=C1./(V_expansion.^gamma); - These statements are used to call the values from the piston kinematics equation code. The value of volumes at different crank angles during expansion stroke(V_expansion) is first calculated and then the corresponding pressures(P_expansion) are calculated to obtain an isentropic curve between points 3 and 4 of the Otto cycle for crank angles from 0 to 180 degrees.

v4=v1

p4=p3*(v3/v4).^gamma

t4=(p4*v4*t3)/(p3*v3) - These commands are used to find temperature, volume and pressure at point 4 respectively.

figure(1)

plot(v1,p1,'r')

hold on

plot(V_compression,P_compression,'r')

hold on

plot(v2,p2,'r')

hold on

plot(v3,p3,'r')

hold on

plot([v2,v3],[p2,p3],'r')

plot(V_expansion,P_expansion,'r')

hold on

plot(v4,p4,'r')

plot([v1,v4],[p1,p4],'r') – These statements are used to plot the Otto cycle. The hold on statement is used so that all points and plots appear on the same graph. We can observe that the isentropic compression and expansion curves are also formed as in the ideal Otto cycle.

Thermal_efficiency = (1 - (1/(r_c.^(gamma-1)))) – This statement is used to calculate the thermal efficiency of the Otto cycle as per theoretical formula.

ERRORS

1)

The above error appears because the value of theta has a range of 50 values while the values of the result will store only one value. This problem can be rectified by using the element wise operator ‘.’ before the mathematical symbols as shown in the code.

2)

This error does not shoot up in the command window but the Otto cycle will not appear correctly. This occurs due to not using ‘hold on’ command while plotting. The hold on command ensures that all plots and points appear in the same graph. This problem can be rectified using the hold on command while plotting the points as shown in the code. The rectified code and graph are shown below.

RESULT

The plot of the Otto cycle is shown below and the value of thermal efficiency is obtained from the command window. The thermal efficiency of the cycle from the entered data is 0.62909 or 62.9 percent.

CODE

The below code is for Piston Kinematics Equation where the function is later called in the main program. 'piston1' is the name of the saved file which is called in the main program.

function [S]=piston1(bore,stroke,con_rod_length,r_c,start_crank_angle,end_crank_angle)
  theta = linspace(start_crank_angle,end_crank_angle,50);
  crank_pin_radius = stroke/2;
  R=r_c/crank_pin_radius;
  v_swept = (pi/4)*bore*bore*stroke;
  v_clear = v_swept/(r_c - 1);
  
  
  term_1 = 0.5*(r_c -1);
  term_2 = R+1-cosd(theta);
  term_3 = ((R^2)-(sind(theta).*sind(theta))).^0.5;
  
  
  S = (1+(term_1*(term_2 - term_3)))*v_clear;
  
  
endfunction

The below code is the main program for obtaining the Otto cycle and calculating the thermal efficiency.

% Otto cycle plot formation and calculation of thermal efficiency

% Standard Inputs

p1=101325;
t1=600;
bore=0.15;
stroke=0.17;
con_rod_length=0.2;
r_c=12;
gamma=1.4;
t3=2500;

% Calculation of swept and clearance volumes based on engine parameters

vswept=(pi/4)*bore*bore*stroke
vclear=vswept/(r_c-1)

v1=vswept+vclear
v2=vclear

p2=p1*(v1/v2).^gamma

C=p1*(v1^gamma); % Calculation of constant C " P1V1^gamma = C " for compression curve

% Calling the function from Piston kinematics code for curve between points 1 and 2

V_compression=piston1(bore,stroke,con_rod_length,r_c,180,0);
[P_compression]=C./(V_compression.^gamma); % P


t2=(p2*v2*t1)/(p1*v1)
v3=v2
p3=(p2*v2*t3)/(t2*v3)


C1=p3*(v3^gamma); % Calculation of constant C1 " P3V3^gamma = C1 " for expansion curve

%  Calling the function from Piston kinematics code for curve between points 3 and 4

V_expansion=piston1(bore,stroke,con_rod_length,r_c,0,180);
[P_expansion]=C1./(V_expansion.^gamma);


v4=v1
p4=p3*(v3/v4).^gamma
t4=(p4*v4*t3)/(p3*v3)

% Plotting the Otto cycle 

figure(1)
plot(v1,p1,'r')
hold on
plot(V_compression,P_compression,'r')
hold on
plot(v2,p2,'r')
hold on
plot(v3,p3,'r')
hold on
plot([v2,v3],[p2,p3],'r')
plot(V_expansion,P_expansion,'r')
hold on
plot(v4,p4,'r')
plot([v1,v4],[p1,p4],'r')

% Calculation of thermal efficiency

Thermal_efficiency = (1 - (1/(r_c.^(gamma-1))))