Solving Otto Cycle Using MatLab

AIM:

To write a MatLab program to solve an Otto Cycle and make plots for it .

OTTO CYCLE:

The Otto Cycle is named after Nikolaus Otto (1832–1891, German) who is credited as the first creator of a petroleum fuel based internal combustion engine operating under a four stroke cycle. Otto’s engine work during the eve of the industrial revolution was a critical development leading to the dawn of the automotive industry in the late 19th century. 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 Otto cycle is a description of what happens to a mass of gas as it is subjected to changes of pressure, temperature, volume, addition of heat, and removal of heat. The mass of gas that is subjected to those changes is called the system. The system, in this case, is defined to be the fluid (gas) within the cylinder. By describing the changes that take place within the system, it will also describe in inverse, the system's effect on the environment. In the case of the Otto cycle, the effect will be to produce enough network from the system so as to propel an automobile and its occupants in the environment.

Figure.1

Today the Otto Cycle is frequently associated with an ideal cycle analysis of automotive spark ignition gasoline internal combustion. This ideal cycle assumes that the combustion process occurs instantaneously at constant volume. To a first order this Otto Cycle analysis can provide insights into how engine efficiency is affected by changing operating and design parameters. Engines operating under the Otto Cycle are the principal power-plants used in automobiles worldwide at the present time, although in some regions of the world the Diesel Cycle is a very close second in automotive applications. Engines operating under the Otto cycle are relatively inexpensive to manufacture and provide reasonable power to weight ratios with good efficiencies. In recent years engine design improvements coupled with exhaust after-treatment have resulted in vehicles that produce very low levels of pollutants.

In Otto Cycle following processes occurs :-

Figure.2

Figure.3

1. Process 0–1 a mass of air is drawn into piston/cylinder arrangement at constant pressure.
2. Process 1–2 is an adiabatic (isentropic) compression of the charge as the piston moves from bottom dead center(BDC) to top dead center (TDC).
3. Process 2–3 is a constant-volume (isochoric) 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.
4. Process 3–4 is an adiabatic (isentropic) expansion (power stroke).
5. Process 4–1 completes the cycle by a constant-volume (isochoric) process in which heat is rejected from the air while the piston is at bottom dead center.
6. Process 1–0 the mass of air is released to the atmosphere in a constant pressure process.

GOVERNING EQUATION:

Figure.4

Now to calculate State variables i.e. Volume(V), Pressure(R), Temperature(T) at different state points following equations are followed :

1. In an ideal Otto cycle, a fuel/air mixture is compressed adiabatically from a volume, `V_1`, to a volume, `V_2`, during the compression phase. The gas traces out the curve between the points 1 and 2 in the P-V diagram . Since from point 1 to 2 it is isentropic process so

`P*V^gamma =`constan`, therefore

`P_1/P_2 = (V_2/V_1)^gamma` ......... (i)

Where,

`P_1`= Pressure at state point 1

`V_1`= Volume at state point 1

`P_2`= Pressure at state point 2

`V_2`= Volume at state point 2

Now from Ideal Gas Equation,

`(P*V)/T =` constan`, therefore

`(P_1*V_1)/T_1 = (P_2*V_2)/T_2` ..........(ii)

`T_2= T_1((P_2*V_2)/(P_1*V_2))`

By substituting (`P_1/P_2`) from equation (i)

`T_2=T_1(V_1/V_2)^(gamma-1)`

`T_2=T_1*(r)^(gamma-1)` ……….(iii)

Where,

`T_1`= Temperature at state point 1

`T_2` = Temperature at state point 2

r = compression ratio =`V_1/V_2`

2. At the end of the compression phase, Point 2, a spark ignites the mixture, which, ideally, is assumed to burn instantaneously. Both temperature and pressure rise instantaneously, the latter reaching the value indicated in Point 3 of the diagram. There is no change in volume during this heat addition or combustion. Therefore,

`V_3=V_2` ………… (iv)

Where,

`V_3` = Volume at State point 3

Now from Ideal Gas Equation,

`(PV)/T` = Contant, therefore

`(P_2*V_2)/T_2 = (P_3*V_3)/T_3` ..…… (v)

Where,

`T_3` = Temperature at state point 3

`P_3` = Pressure and state point 3

3. The increased high pressure exerts a force on the piston and pushes it towards the BDC. Expansion of working fluid takes place isentropically from Point 3 to Point 4 and work is done by the system on the piston. The volume ratio `(V_4/V_3)` is called the "isentropic expansion ratio". [For the Otto cycle is the same as the compression ratio (V1/V2)]. Mechanically this is the expansion of the hot gaseous mixture in the cylinder known as expansion (power) stroke.

Now here,

`P.V^gamma`= Constant, so

`P_3/P_4 = (V_4/V_3)^gamma` ……… (vi)

`V_4/V_3`= expansion ratio = compression ratio=`V_1/V_2` = r

Where,

`P_4`= Pressure at state point 4

`V_4`= Volume at state point 4

From Ideal Gas Equation,

`(P*V)/T` = Constant,

`(P_3*V_3)/T_3 = (P_4*V_4)/T_4` ……….. (vii)

`T_3= T_4 ((P_3*V_3)/(P_4*V_4))`

By substituting (`P_3/P_4`) from equation (vi)

`T_3 = T_4(V_4/V_3)^(gamma-1)`

`T_3 = T_4*(r)^(gamma-1)` ……….. (viii)

Where,

`T_4`= Temperature at state point 4

4. The piston is momentarily at rest at BDC. The working gas pressure drops instantaneously from point 4 to point 1 during a constant volume process as heat is removed to an idealized external sink that is brought into contact with the cylinder head. In modern internal combustion engines, the heat-sink may be surrounding air (for low powered engines), or a circulating fluid, such as coolant. The gas has returned to state 1. Therefore,

`V_4 = V_1` ……….. (ix)

5. Efficiency can be calculated from the heat supplied to engine and heat rejected by the engine.

Now from processes 2 to 3 combustion takes place heated is added here,

Therefore,

`Q_(IN) = Q_23=m*C_v*(T_3-T_2)` ………….(x)

Where,

`Q_(IN)`= heat supplied in process 2-3

`m` = mass of the air-fuel mixture

`C_v` = specific heat at constant volume.

Now from process 4 to 1 heat is rejected to the surroundings,

Therefore’

`Q_(OUT)=Q_41= m*Cv*(T_4-T_1)` ……….(xi)

Where,

`Q_(OUT)` = heat rejected in process 4-1

Therefore, Net Workdone is

`W_(NET) =Q_(IN) - O_(OUT)` ............(xii)

Where,

`W_(NET)`= Net work done by the gas

Now, thermal efficiency is

`eta`= Wnet/Qin = Qin-Qout/Qin

by substituting equation (xi), (xii) and solving, we will get

`eta=1- (T_4-T_1)/(T_3-T_2)` ………..(xiii)

Where,

`eta` = Thermal Efficiency

from equation (iii) and (viii)

`T_2/T_1= T_3/T_4=r^(gamma-1)`

Threrefore,

`T_2/T_3=T_4/T_1`

By substituting above equation in equation(xiii)

`eta = 1- [(T_1(T_4/T_1 – 1))/(T_2(T_3/T_2 -1))]`

= 1- `T_1/T_2`

`eta = 1- 1/r^(gamma-1)` ..........(xiv)

6. In the above graph distance from origin to point `V_2` is known as clearance volume, this gap between piston and cylinder head and distance from `V_2` to`V_1` is known as swept volume, this where piston is free to move. So total volume of Cylinder is

`V= V_c+V_s` ……….. (xv)

Where,

V= total volume of cylinder

`V_c=` clearance volume =`V_s /(r-1)`

`V_s`= swept volume = `pi/4``*``(`diameter of the bore`)^2``*` length of the storke

6. To find how volume inside combustion chamber is changing ,a formula is used which is derived from geometric relationships e using bore of cylinder, stroke of the piston, compression ratio, connecting rod length of piston and crank angle.

`V/V_c = 1+ 1/2(r-1)[R+1-cosθ-(R^2 – sin^2θ)^(1/2)]` ……… (xvi)

Where,

R=(connecting rod length) `/` (crank pin radius)

Let,

a = crank pin radius = (stroke of piston) `/` 2

Therefore,

R=connecting rod length `/` a

8. Using the above calculations and equations a program is written in MatLab to solve the Otto cycle.

OBJECTIVES:

1. To create a P-V diagram using MatLab Program.
2. To calculate Thermal Efficiency of the engine using MatLab Program.

PROGRAM:

MAIN CODE:

clear all
clc
close all
 
%Inputs
 gamma=1.4;
 t3=2000;
 
%STATE VARIABLES AT STATE POINT 1
 p1=101325;
 t1=500;
 
%Engine geometric parameters
 bore= 0.1;
 stroke= 0.1;
 con_rod = 0.15;
 r = 12;
 
%Calculating the swept volume and the clearance volume
 v_swept=(pi/4)*bore^2*stroke;
 v_clearance= v_swept/(r-1);
 
 v1= v_swept + v_clearance;
 v2= v_clearance;
 
 
% STATE VARIABLES AT STATE POINTS 2
% p2v2^gamma=p3v3^gamma
% p2=p1*(v1/v2)^gamma = p1*r^gamma

 p2=p1*r^gamma;

%p1v1/t1 = p2v2/t2 | t2=p2*v2*t1/(p1*v1) 
 t2=p2*v2*t1/(p1*v1);
 
%STATE VARIABLES AT STATE POINTS 3
 v3=v2;

%p3v3/t3=p2v2/t2 | p3=p2*t3/t2
 p3=p2*t3/t2;
  
%STATE VARIABLES AT STATE POINT 4
 v4=v1;
 
%p4v4^gamma=p3v4^gamma | p4=p3*(v3/v4)^gamma
 p4=p3*(v3/v4)^gamma;
  
 figure(1)
 hold on
 plot(v1,p1,'*' , 'color','r')
 plot(v2,p2,'*' , 'color','r')
 plot(v3,p3,'*' , 'color','r')
 plot(v4,p4,'*' , 'color','r')
  
 figure(2)
 plot([v1 v2 v3 v4],[p1 p2 p3 p4],'color','b')

EXPLANATION:

1. Generally for solving problems on Otto cycle, initial temperature, initial pressure, peak pressure, gamma, compression ratio and geometric parameters are given. So in the above code all these values are entered as input values.
2. In the above code ‘con_rod’ refers to length of connecting rod of piston.
3. First code is started with ‘close all’, ‘clear all’ , ‘clc’ commands and then Swept Volume(v_swept) and Clearance Volume(v_clearance) are entered in the code as per the formulae mentioned in above equation (xv) .
4. Volume at state point 1 and state point 2 are determined using swept volume and clearance volume.
5. Referring to above equation(i) and equation(ii) temperature and pressure at state point 2 are calculated.
6. As peak temperature i.e. is temperature at state point 3 (t3) is given and entered in above code as input , referring to above equation(vi) and equation(v) pressure and volume at state point 3 are determined.
7. Using above equation(vi) and equation(ix) volume and pressure at state point 4 is determined.
8. Using all the volume and pressure values at all state points, P-V graph is plotted using ‘plot()’ command and this command is also used to join the these state points. In the plot command ‘*’ is used to represent the each state point in the graph.
9. Using ‘figure()’command the plots and the cycle diagram are divided into 2 figures.

Figure.5

Figure.6

10. Using above equations all the state points are determined and plotted.
11. After running the code the results are as in the above figures.
12. In Figure.5 we can observe that the plots are plotted accurately according to the values that are shown in command window.
13. In Figure.2 the graph of Otto cycle is not correct. From point 1 to point 2 and from point 3 to point 4 the graph has wrong path curves, so plots should not be joined directly using ‘plot()’ command i.e. there should not be a straight line.
14. Correct path equation has to be entered in order to get correct graph. Thus the code has to be corrected.

CODE OF PATH EQUATION:

clc 
close all 
clear all

%Geometric Inputs
bore = 0.1;
stroke = 0.1;
con_rod = 0.15;
r= 12;
a = stroke/2;
R = con_rod/a;

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

theta = linspace(0,180,180);

term1 = 0.5*(r-1);
term2 = R+1 - cosd(theta);
term3 = (R^2 - sind(theta).^2).^0.5;

V = (1+ term1*(term2 - term3)).*v_c;
plot(theta,V)

EXPLANATION:

1. Here the program is started with the close all, clear all and clc commands and then the geometric parameters are given as inputs.
2. In order to give correct path equation we have used equation(xvi) ,which tells us how volume varies from one state point to another state point in isentropic processes.
3. Clearance volume and Swept volume i.e ‘v_c’ and ‘v_s’ are calculated using the formulae stated in equation(xv).
4. Then crank angles are given here using variable ‘theta’ and it is defined using ‘linspace()’ command. Values of crank are defined from 0 degree to 180 degree that is piston is moving from TDC to BDC.
5. The above volume equation is divided in order for easy understanding and if any errors occurs it will be easy to find them. So whenever there are big formulas it is better to divide them and enter in the code.
6. Here it is divided into 3 parts i.e. ‘term1’ , ‘term2’ and ‘term3’ and then they are used to enter the full equation . Since ‘theta’ here is array so elementary wise operator(.) is used.
7. Then plot command is entered to generate the curve and the code is executed.

Figure.7

8. Above figure is the result for the above code and the curve is absolutely correct i.e graph is perfect.
9. Now we can copy-paste this above code in the previous code to obtain the correct path curves in the P-V graph.
10. But we have to copy this code and paste it twice, one time for isentropic compression and another for isentropic expansion which will make the code lengthy and clumsy, if any error occurs it will difficult to find and resolve them.
11. So instead of copy-pasting the code we can covert this code to a function and use it in the previous code, which will simplify our code and it will be easy to handle the code.

FUNCTION CODE:

function[V] =engine_kinematics(bore, stroke, con_rod, r, start_crank, end_crank)

a = stroke/2;
R = con_rod/a;

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

theta = linspace(start_crank,end_crank,100);

term1 = 0.5*(r-1);
term2 = R+1 - cosd(theta);
term3 = (R^2 - sind(theta).^2).^0.5;

V = (1+ term1*(term2 - term3)).*v_c;

end

EXPLANATION:

1. Here the path equation code is converted into a function code.
2. Here volume(V) is define as the function of bore, stroke, con_rod (i.e. connecting rod length) , r (i.e. compression ratio), start_crank (i.e. crank start angle) and end_crank (i.e. crank ending angle) , using the ‘function’ command and this function is named as ‘engine_kinematics’.
3. For example if we type X= engine_kinematics (0.1,0.1,0.15,12,180,360), MatLab will take the value of bore=0.1, stroke=0.1, con_rod=0.15, r = 12 , crank_start=180 degrees, crank_end=360 degrees , will be taken as input and the below equations will be executed using these inputs values thus giving us the result.
4. Here for ‘theta’ variable, change is made in the inputs of ‘linspace()’ command i.e. the angles for theta are define as to take 100 values between ‘start_crank’ to ‘end_crank’. So there will results for 100 angles between start crank angle and end crank angle.
5. All the inputs in the above code which were given fixed values are now converted to variables which can be assigned any value according to the problem.
6. To close the function ‘End’ is used,which completes our function.
7. Now once we use this function in command window we will get the value of ‘V’ directly.

Figure.8

8. From the above the above figure we can observe that the function that we have created is working perfectly.
9. Now save this function with same name that is given to the function i.e. ‘engine_kinematics’ and save it in the same folder where the main code is saved.
10. Now we can used this function in the main code and obtain the correct curve for the isentropic processes.

FINAL CODE:

clear all
clc
close all

%Inputs
 gamma=1.4;
 t3=2000;

%State Variables at state point 1
 p1=101325;
 t1=500;

%Engine geometric parameters
 bore= 0.1;
 stroke= 0.1;
 con_rod= 0.15;
 r=12;

%Calculating the swept volume and the clearance volume
 v_swept=(pi/4)*bore^2*stroke;
 v_clearance= v_swept/(r-1);
 
 v1= v_swept + v_clearance;
 v2= v_clearance;
 
 
%State variables at state point 2
% p2v2^gamma=p3v3^gamma
% p2=p1*(v1/v2)^gamma = p1*r^gamma
 p2=p1*r^gamma;
 
%p1v1/t1 = p2v2/t2 | t2=p2*v2* 
 t2=p2*v2*t1/(p1*v1);
 constant_c = p1*v1^gamma;
 V_compression = engine_kinematics(0.1,0.1,0.15,12,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*v3^gamma;
 V_expansion = engine_kinematics(0.1,0.1,0.15,12,0,180);
   
 P_expansion = constant_c./V_expansion.^gamma;
  
%State Variables at state point 4
 v4=v1;
 
%p4v4^gamma=p3v4^gamma | p4=p3*(v3/v4)^gamma
 p4=p3*(v3/v4)^gamma;
  
 figure(1)
 hold on
 plot(v1,p1,'*' , 'color','r')
 plot(V_compression,P_compression,'color','r')
 plot(v2,p2,'*' , 'color','r')
 plot(v3,p3,'*' , 'color','r')
 plot(V_expansion,P_expansion,'color','g')
 plot(v4,p4,'*' , 'color','r')
 plot([v2 v3],[p2 p3], 'color','b')
 plot([v4 v1],[p4 p1], 'color','y')
 xlabel('Volume')
 ylabel('Pressure')
 
 %Calculating Thermal Efficiency
 eta= 1- (1/r^(gamma-1))
 efficeincy_in_percent= eta*100

EXPLANATION :

1. In the above code we have used the function that we have created previously for defining volume path.
2. Here we have assigned a variable (V_compression) and defined it with function that we have created by giving it necessary inputs.
3. This represents that the change in volume in isentropic compression process e. from point 1 to point 2.
4. Now we also have to define how pressure varies during this process, so for this we have used equation (pv^gamma = constant). This variable is assigned as ‘P_compression’
5. Now here we don’t know the value of ‘constant’ in the above equation, so for calculating this we used two end points of the process 1-2 [i.e. (v1,p1) and (v2,p2)] and used the same above equation (i.e p1v1^gamma=constant) thus value of ‘constant’ is calculated.
6. Now we have value of ‘V_compression’ and ‘constant’ we can find the value of ‘P_compression’ using equation ‘ P_compression*V_compression^gamma=constant’. Elementary vise operator(.) is used since ‘V_compression’ is an array.
7. Similar method is followed for the process 3-4 which is isentropic expansion ,but the difference is that the variables are given different names i.e. ‘V_expansion’ and ‘P_expansion’ for defining change in Volume and Pressure respectively during this process. For finding ‘constant’ value, here end points [(v3,p3) and (v4,p4)] are used.
8. And also the starting crank angles and ending crank angles for both the processes are different. So always check the angles entered are correct. For compression process the crank angles are from 180 degrees to 0 degree and for expansion process crank angles are 0 degree to 180 degree.
9. Now using the values of V_compression and P_compression, the path is defined using ‘plot()’ command. Similarly using values of V_expansion and P_expansion a curve is plotted.
10. Now a ‘hold on’ command is also used before plot command so the plot commands do not get overwritten by last plot command.
11. Figure (2) is removed from code as it a wrong representation of the graph. Only state point 2 and state point 3, state point 4 and state point1 are joined using plot command as they are straight lines.
12. Each curve is given different color using ‘color’ command for better representation of the graph, and it differentiates different processes that are occurring in Otto cycle.
13. Now to find efficiency of the above Otto cycle ,formula for thermal efficiency is entered according to equation(xiv) and then it is multiplied by hundred to obtain the value in percentage.

Figure.9

14. The above figure shows us the correct P-V graph and we can observe how volume and pressure changes in the isentropic processes.

Figure.10

15. We have got the value of thermal efficiency of the Otto cycle in the command window after we executed the code.

CONCLUSIONS:

1. Thus thermal efficiency and P-V graph of an Otto cycle are obtained using a MatLab program.
2. Similar approach can be used to solve any other Air-Standard cycles.
3. Different input values can be given according to our problem to solve the Otto cycle.