Otto Cycle Engine Kinematics using Matlab

Aim:

The main aim of this project is to write a Matlab program for Engine parameters of an Otto cycle engine whose variables like Inlet temperature(T1), pressure(p1) and temperature(T3) at the end of expansion are defined and other parameters are computed with respective formulas.

The Engine design parameters like bore, stroke, connecting rod length, compression ratio, start crank angle, end crank angle were provided. The tasks of this program are:-

1.Create a PV (pressure vs velocity) diagram using Matlab

2. Determine the Thermal Efficiency of this cycle.

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. It consists of four processes which occur in it.

  The processes are described as follows: -

1. Process 1–2 is an isentropic compression of the charge as the piston moves from Bottom Dead Centre (BDC) to Top Dead Centre (TDC). In this process, the piston moves from bottom dead centre (BDC) to top dead centre (TDC) position. Air undergoes reversible adiabatic (isentropic) compression. We know that compression is a process in which volume decreases and pressure increases. Hence, in this process, the volume of air decreases from v1 to v2 and pressure increases from p1 to p2. Temperature increases from t1 to t2. 

2. Process 2–3 is a constant-volume heat transfer to the working gas from an external source while the piston is at top dead centre. This process is intended to represent the ignition of the fuel-air mixture and the subsequent rapid burning. Here, piston remains at top dead centre for a moment. Heat is added at constant volume v2 = v3 from an external heat source. Temperature increases from t2 to t3, pressure increases from p2 to p3.  

 3. Process 3–4 is an isentropic expansion (power stroke). In this process, air undergoes isentropic (reversible adiabatic) expansion. The piston is pushed from top dead centre (TDC) to bottom dead centre (BDC) position. Here, pressure decreases from p3 to p4, volume rises from v3 to v4, temperature falls from t3 to t4. 

4. 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 centre. The piston rests at BDC for a moment and heat is rejected at constant volume v4=v1. In this process, the pressure falls from p4 to p1, temperature decreases from t4 to t1.

The following assumptions were made before start of program. These are also called as Input Parameters. which are as follows: -

Pressure at State point 1 (p1) = 101325 Pascals

Temperature at State point 1 (t1) = 500 Kelvin

 Gamma (`gamma`) = 1.4

Temperature at State point 3 (t3) = 2300 Kelvin

Geometric Parameters: -

Bore = 0.1 meters

Stroke = 0.1 meters

Connecting rod = 0.15 meters

Compression Ratio (cr) = 10

swept volume `(v_s) = (pi/4)((b)^2)(stroke)` where b is the bore

clearance volume `(v_c) = v_s/(cr-1)`

Volume at State point 1 (v1) = `v_s + v_c`

State Point 2: -

After the compression process, the volume in the combustion chamber is equal to the clearance volume

v2 = `v_c`

From the isentropic relation `p1v1^gamma = p2v2^gamma`,  the pressure can be computed

`p2 = (p1(v1^gamma))/(v2^gamma)` According to Ideal gas law,  

`(p2v2)/(t2) = (p1v1)/(t1)`  

Lets assume that `(p1v1)/(t1) = rhs`, so

`t2 = (p2v2)/(rhs)`  

State Point 3: -

Since it's constant volume heat addition,

v3 = v2

Similar to the earlier method,  

`p3 = rhs((t3)/(v3))`

State Point 4: -

Similarily,

v4 = v3

`t4 = (p4v4)/(rhs)`  

Inorder to calculate the relation between crack angle and the volume in the combustion chamber a function is defined which calculate the relation between the crank angle and the volume in the combustion chamber. Which is as follows,

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

So basically ,it is necessary to determine the compression volume. With the compression volume, using relations between compression pressure and compression volume the compression pressure can be determined. For compression, crank angle starts from 180 degrees and ends at 0 degrees.

In an adiabatic process PV^gamma = Constant,

From this,  Compression pressure = Constant /v^gamma.

Using the above method, Expansion volume and pressure can be determined, but in this case, the crank angle starts from 0 degrees and ends at 180 degrees. With the calculated values, we plot the following PV diagram

Finally,  the thermal efficiency of this cycle can be determined by the following equation,

`eta = 1 - (1/(cr^(gamma-1)))`

`eta = 1 - (cr^(1-gamma))`

clear all 
close all
clc

%Inputs
gamma = 1.4;
T3 = 2300;

% state variables
 p1 = 101325;
 T1 = 500;
 
 %Engine Geometric Parameters
 
 bore = 0.1;
 stroke = 0.1;
 con_rod = 0.15;
 cr = 12;
 
 %calculating the swept volume and the clearance volume
 
 v_swept = (pi/4)*bore^2*stroke;
 v_cleareance = v_swept/(cr-1);
 
 v1= v_swept + v_cleareance;
 v2= v_cleareance;
 
 %state variables at state 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 = otto_cycle_engine_kinematics(bore,stroke,con_rod,cr,180,0);
 
 P_compression = constant_c ./((V_compression).^gamma) ;
 
 
 
 
 %state variables at point 3
 
 v3 = v2;
 
 %p3v3/t3 = p2v2/t2  |  p3 = p2*t3/t2
 
 p3 = p2*T3/T2;
 
 %state variables at state point 4
 
 v4 = v1;
 
 %p3v3^gamma = p4v4^gamma | p4 = p3*(v3/v4)^gamma
 
 p4 = p3*(v3/v4)^gamma;
 
 constant_c = p3*v3^gamma;
 
 V_expansion = otto_cycle_engine_kinematics(bore,stroke,con_rod,cr,0,180);
 
 P_expansion = constant_c ./((V_expansion).^gamma) ;


% Thermal efficiency
 
 eta = 1 - (cr^(1-gamma));
 
 
 figure(1)
 hold on
 
 plot(v1,p1,'*','color','r')
 plot(V_compression,P_compression)
 plot(v2,p2,'*','color','r')
 plot([v2 v3],[p2 p3])
 plot(v3,p3,'*','color','r')
 plot(V_expansion,P_expansion)
 plot(v4,p4,'*','color','r')
 plot([v4 v1],[p4 p1])
 xlabel('VOLUME')
 ylabel('PRESSURE')
 
 figure(2)
 plot([v1 v2 v3 v4 v1],[p1 p2 p3 p4 p1],'color','b')
 

Function Program represents Volumetric change with respect to crank angle:

function [V]= otto_cycle_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

The thermal efficiency of this Otto cycle is 62.99%.