MATLAB Programme To Solve An Otto Cycle

AIM:

To write a program in MATLAB to plot an otto cycle and calculate its thermal efficiency.

GOVERNING EQUATIONS:

• Compression ratio : Vs+Vc/Vc

Where Vs=swept volume=(pi/4)*d^2*L

Where d=bore

L=stroke

Vc=clearance volume

• Volume change inside the cylinder due to various processes :

• For isoentropic process, P*V^gamma=constant.
• Thermal efficiency of otto cycle=1-(1/cr^(gamma-1))

where cr=compression ratio=10 (assumed)

gamma=specific heat of air=1.4

OBJECTIVE:

To solve the different state variables of otto cycle.

To plot a P-V diagram.

To calculate the thermal efficiency of an otto cycle.

THEORY:

OTTO CYCLE:

The Air Standard-Otto cycle is the idealized cycle for the spark-ignition Internal Combustion engines. In the Process 0-1, the piston moves from Top Dead Centre (TDC) to Bottom Dead Centre (BDC) and the mixture of fuel and air (charge) is sucked inside the combustion chamber.

Process 1-2: Isentropic compression

In this process, the piston moves from BDC to TDC. Air undergoes reversible adiabatic (isentropic) compression. Compression is a process in which volume decreases and pressure increases. Therefore, volume of air decreases from V1 to V2 and pressure increases from p1 to p2 and temperature increases from T1 to T2.

Process 2-3: Constant volume heat addition

It is an isochoric heat addition process. In this process, the piston remains at TDC for a moment. Heat is added at constant volume (V2 = V3) from an external heat source (By creating a spark by means of a spark plug). Temperature further increases from T2 to T3 and pressure also increases from p2 to p3.

Process 3-4: Isentropic Expansion

In this process, air undergoes reversible adiabatic (isentropic) expansion. The piston is pushed from TDC to BDC. Here, pressure decreases from p3 to p4, volume increases from v3 to v4 and temperature falls from T3 to T4.

Process 4-1: Constant volume heat rejection

In this process, the piston rests at BDC for a moment and heat is rejected at constant volume (V4=V1). During this process, the pressure falls from p4 to p1 and temperature further decreases from T4 to T1.

ASSUMPTIONS:

• The initial pressure P1 and temperature T1 is assumed to be 101325 pascal and 500 kelvin.
• The maximum temperature inside the engine cylinder. T3 is assumed to be 2300 kelvin.

EXPLANATION OF MATLAB CODE:

Step 1: Declaring constants, predetermined state variables and engine parameters: Constant gama (ratio of sp. heat at const. pressure to sp. heat at const. volume) is declared as 1.4. p1, t1 are predetermined state variables. Temperature at the end of Isochoric heat addition process is declared as 2300k (t3=2300k). In Engine parameters, variables such as ‘bore’, ‘stroke’, ‘con_rod’ and ‘cr’ are also declared.

Step 2: Calculating the swept volume and clearance volume:

Swept volume and clearance volume are calculated. When the piston is at BDC, the volume inside the combustion chamber (v1) is equal to the summation of swept volume and clearance volume. And when the piston is at TDC, the volume inside the combustion chamber (v2) is equal to clearance volume.

Step 3: Calculating state variables:

1. Process 1-2: p2 is calculated with the help of isentropic law (p1v1^gama = p2v2^gama) and t2 is calculated with the help of gas law (p1v1/t1 = p2v2/t2).
2. Process 2-3: p3 is calculated with the help of Gay-Lussac’s law (p3/p2 = t3/t2). Since this is isochoric process, v3 = v2.
3. Process 3-4: p4 is calculated with the help of isentropic law (p1v1^gama = p2v2^gama) and t4 is calculated with the help of gas law (p3v13/t3 = p4v4/t4). Since Process 4-1 is an isochoric process, therefore, v4 = v1.

Step 4: Computation of Isentropic Processes

A formula for determining volume inside combustion chamber at any crank angle is used.

It is: (v/vc)=1+0.5*(cr−1)[R+1−cos(theta)−(R2−sin(theta)^2)^0.5]

Where,

R = con_rod_length/crank_pin_radius,

crank_pin_radius = stroke/2,

cr is compression ratio,

theta is the crank angle.

This formula is broken into three terms for simplicity:

term1 = 0.5*(cr-1)

term2 = R+1-cosd(theta)

term3 = (R^2 - sind(theta)^2)^0.5

therefore, volume inside combustion chamber at any crank angle can be computed by the formula:

V = (1 + term1*(term2-term3))*Vc

This formula can be used inside a function

[v] = engine_kinematics(bore,stroke,con_rod,cr,start_crank,end_crank) where bore, stroke, con_rod, cr, start_crank and end_crank are the independent variables and volume (v) is the dependent variable. This command can be used as template to call the function (engine_kinematics) and compute the volume inside combustion chamber at any given crank angle.

This command is very helpful in computing the pressure and volume during isentropic compression and expansion.

Step 5: Calculating Thermal efficiency:

Thermal efficiency is calculated by the formula =(1-(1/(cr)^(gamma-1)))*100

where cr= compression ratio.

Step 6: PV Diagram

Plotting of Isentropic and Isochoric paths:

Commands such as "plot(v_compression,p_compression,'linewidth','3')" and "plot([v2 v3],[p2 p3],'color','m')" to plot the graph for isentropic and iscochoric process respectively.

Plotting of state points:

Coomands such as "plot(v1,p1,'*','color','r')" are used to plot all the state points.

MAIN CODE:

%To calculate the state variables in an otto cycle and to calculate its thermal efficiency
clear all
close all
clc
%Input variables

%gamma for air=1.4
gamma=1.4;
%maximum temperature inside engine cylinder unit-kelvin
T3=2300;

%State variables at state point 1

%Atmospheric pressure unit-pascal
P1=101325;
%Temperature unit-kelvin (assumed)
T1=500;

%Engine geometric parameters
bore=0.1; %bore length unit=metre
stroke=0.1; %stroke length unit=metre
con_rod=0.15; %connecting rod length unit=metre
cr=10; %compression ratio

%Calculating swept and clearance volume
V_swept=(pi/4)*bore^2*stroke;
V_clearance=V_swept/(cr-1);

V1=V_swept+V_clearance;
V2=V_clearance;

%State variable at state point 2

%P2*V2^gamma=P1*V1^gamma | P2=P1*(V1/V2)^gamma | P2=P1*cr^gamma
P2=P1*(cr^gamma);

%P1*V1/T1=P2*V2/T2 | T2=(P2*V2*T1)/(P1*V1)
T2=(P2*V2*T1)/(P1*V1);

C_const=P1*(V1^gamma);
V_compression=engine_kinematics(bore,stroke,con_rod,cr,180,0);
P_compression=C_const./(V_compression.^gamma);


%State variable at state point 3

V3=V2;

%P3*V3/T3=P2*V2/T2 | P3=(P2*V2*T3)/(T2*V3)
 P3=(P2*V2*T3)/(T2*V3);
 
 C_const=P3*(V3^gamma);
 V_expansion=engine_kinematics(bore,stroke,con_rod,cr,0,180);
 P_expansion=C_const./(V_expansion.^gamma);
 
%State variable at state point 4

V4=V1;

%P3*V3^gamma=P4*V4^gamma | P4=P3*(V3/V4)^gamma
P4=P3*(V3/V4)^gamma;

%Calculating thermal efficiency
therm_eff=(1-(1/(cr)^(gamma-1)))*100;
fprintf('The thermal efficiency of the otto cycle is : %f percentage',therm_eff);

%plotting state variables
hold on
plot(V1,P1,'*','color','r')
text(V1,P1,'Point 1')
plot(V_compression,P_compression,'linewidth',3)
plot(V2,P2,'*','color','r')
text(V2,P2,'Point 2')
plot(V3,P3,'*','color','r')
line([V2 V3],[P2 P3],'linewidth',3)
plot(V_expansion,P_expansion,'linewidth',3)
text(V3,P3,'point 3')
plot(V4,P4,'*','color','r')
text(V4,P4,'point 4')
line([V4 V1],[P4 P1],'linewidth',3)

CODE FOR ENGINE KINEMATICS FUNCTION:

%Function for engine kinematics
function [V]=engine_kinematics(bore,stroke,con_rod,cr,start_crank,end_crank)
  %Crank pin radius
  a=stroke/2;
  %Crank radius
  r=con_rod/a;
  
  %calculating swept volume and clearance volume
  V_swept=(pi/4)*bore^2*stroke;
  V_clearance=V_swept/(cr-1);
  
  %how volume changes inside a engine cylinder
  
  %making linearly spaced vector of crank angle
  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_clearance;
  
end

OUTPUT:

CONCLUSION:

Thus we plotted the P-V diagram from the given assumed information and the thermal efficiency was found to be 60.189283 %