AIR STANDARD OTTO CYCLE USING MATLAB

OBJECTIVE:

To plot the Otto cycle and calculate its thermal efficiency using MATLAB

INTRODUCTION:

 Internal combustion engine is a heat engine which converts chemical energy of fuel in to mechanical energy usually in the form of rotating shaft. Chemical energy of fuel is converted into thermal energy first by means of combustion of fuel, which raises pressure and temperature of gases inside the cylinder. This high pressure gases expands by converted into rotating crank shaft by the mechanical linkages of the engine.

IC engine in which the combustion of fuel takes place inside the engine cylinder are non cyclic heat engines in which the working fluid undergoes permanent chemical change due to combustion, and these products of combustion after doing work are thrown out of engine and fresh charge is taken to complete thermodynamic cycle . To simplify the analysis, the thermodynamic cycle is treated as air standard cycle.

Assumptions of Air Standard Cycles: 

1. The working fluid is taken as air which operates in a complete thermodynamic cycle.

2. Air is assumed to behave as an ideal gas, and its specific heats are assumed to be constant.

3. All the processes in the cycle are completely reversible.

4. Combustion and Exhaust processes are treated as Heat Addition and Heat Removal processses which exchange heat with external heat reservoirs.

The Otto cycle:

The Otto cycle is the air standard cycle of Spark Ignition (SI) engine which is used in automobiles. The `P` - `v` and `T` - `s` diagram of the Otto cycle is shown in the following figure

Process 0 - 1: Constant Pressure Fuel Intake: The inlet valve is open, the piston moves to the right, admitting fuel-air mixture into the cylinder at constant pressure. 

Process 1 - 2: Isentropic Compression: Both the valves are closed, the piston, moves from BDC to TDC, compresses the combustible mixture to the minimum volume.

Process 2 - 3: Constant Volume Heat Addition or Combustion:  The mixture is then ignited by means of a spark, combustion takes place, and there is an increase in temperature and pressure.

Process 3 - 4: Isentropic  Expansion: The products of combustion do work on the piston which moves from TDC to BDC, and the pressure and temperature of the gases decrease.

Process 4 - 1: Constant Volume Heat Removal or Blow-down: The exhaust valve opens, and the pressure drops to the initial pressure.
Process 1 - 0: Constant Pressure Exhaust: With the exhaust valve open, the piston moves inwards to expel the combustion products from the cylinder at constant pressure.

The series of processes as described above constitute a mechanical cycle, and not a thermodynamic cycle. To treat it as thermodynamic cycle, intake (process 0 - 1) and exhaust (process 1 - 0) processes are neglected. The cycle is completed in four strokes of the piston.

 INPUT:

Initial Pressure = `P_1` = 1 atm = 101325 kPa

Initial Temperature = `T_1` = 500 K

Maximum Temperature = `T_3` = 3500 K

Adiabatic Index = `gamma` = 1.4

Cylinder bore = `d` = 0.1 m

Cylinder stroke = `l` = 0.1 m

Length of connecting rod = `l_c` = 0.2 m

Compression ratio = `r_c` = 4,6,8

FORMULAE:

1. Swept volume  = `V_s`  = `pi/4``d^2``l`

2. Clearance volume  = `V_c`  = `V_s/(r_c-1`

3. Reversible adiabatic or isentropic process equation = `Pv^gamma = constant`

4. Combined gas law equation = `(Pv)/T = constant`

5. Thermal efficiency of Otto cycle = `eta = 1 - 1/(r_c^(gamma - 1)`

6. Volume trace inside the engine combustion chamber is given by

        `V = (1 + 0.5(r_c - 1[R + 1 - cos theta - sqrt(R^2 - sin^2 theta)))V_c`

     where,

         `theta` - Crank angle

         `R` - Ratio of connecting rod length to crank pin radius `a`

         `a` - stroke length/2

CODE EXPLANATION:

1. Known inputs `P_1`,`T_1`,`T_3`,`gamma`, bore, stroke of cylinder, length of connecting rod and compression ratio are entered.

2. Swept volume and clearance volume are caluclated using above formulae (1) and (2).

3. Initial state 1 volume is calculated as it is sum of swept volume and clearance volume since at state 1, piston is at bottom dead center (BDC) and state 2 volume is taken as clearance volume itself since piston will be at top dead center (TDC) at the end of compression. 

        `v_1 = v_s + v_c`

        `v_2 = v_c`

4. Volumes of state 2 and state 3 are equal as the process 2-3 is constant volume process and similarly volumes of state 1 and state 4 are equal.

        `v_1 = v_4`

        `v_2 = v_3`

5. Process 1-2 is isentropic compression:

Since the process is isentropic, formula (3) is used to calculate pressure `P_2` and formula (4) is used to calculate temperature `T_2`

        `P_2 = P_1r_c^gamma`

        `T_2 = (P_2v_2T_1)/(P_1v_1)`

If (`P_1`, `v_1`) and (`P_2`, `v_2`) are plotted simply using plot command it will give a straight line between state 1 and 2 which actually should be hyperbola. In order to get that plot, volumes at intermediate points between state 1 and 2 are calculated using volume trace equation for a particular crank angle (`theta` is from 180 to 0 degrees) and pressures at these intermediate points are calculated using adiabatic process equation. Using (`P`,`v`) at these intermediate points, process 1-2 is plotted.

6. Process 2-3 is constant volume heat addition

Combined gas law is used between state 2 and 3 to calculate `P_3`

        `P_3 = (P_2T_3)/T_2`

7. Process 3-4 is isentropic expansion 

The process follwed for process 3-4 is similar to process 1-2 since both are isentropic processes and `P - v` plot is plotted using volume trace equation for which crank angle is from 0 to 180 degrees.

        `P_4 = P_3(1/r_c)^gamma`

        `T_4 = (P_3v_3T_4)/(P_4v_4)`

8. Thermal efficiency of Otto cycle is calculated using above formula (5).

9. Finally all states points are plotted to get `P-v` plot

CODE:

% Program to plot otto cycle PV diagram and Calulation of its efficiency
clear all
close all
clc

% Inputs
gamma = 1.4;

% State variables
P1 = 101325;                % in Pa
T1 = 500;                   % in Kelvin
T3 = 3500;                  % in Kelvin

% Cylinder Dimensions
bore = 0.1;                 % in meter
stroke = 0.1;               % in meter
con_rod = 0.2;              % in meter
cr_ratio = 6;               % cr ratio varies from 4 to 8 for SI engine

% Calculating Swept volume and Clearance volume
v_s = (pi/4)*(bore^2)*stroke;
v_c = v_s/(cr_ratio - 1);

% Calculation of volumes at each state
% Both State 1 and State 4 for are when piston is at Bottom Dead Center(BDC)
% At BDC, Total volume = Swept volume + Clearance volume 
v1 = v_s + v_c;
v4 = v1;

% Both State 2 and State 3 for are when piston is at Top Dead Center(TDC)
% At TDC, Total volume = Clearance volume 
v2 = v_c;
v3 = v2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Process 1-2 is isentropic compression.
% P1*v1^gamma = P2*v2^gamma = constant
% P2 = P1*(v1/v2)^gamma = P1*(cr_ratio)^gamma
P2 = P1*(cr_ratio)^gamma;
% Between end states 1 and 2, combined gas law equation is used
% (P1*v1)/T1 = (P2*v2)/T2
T2 = (P2*v2*T1)/(P1*v1);

% Since process 1-2 is an isentropic process, we cannot join directly end
% states 1 and 2 in the plot which gives straight line which is wrong.
% Hence we use formula which determines how volume inside the engine combustion
% chamber changes wrt to crank angle, cylinder dimensions & compression ratio
constant_c = P1*v1^gamma;
v_compression = engine_kinematics(bore,stroke,con_rod,cr_ratio,180,0);
% Calculation of pressure at each point of v_compression value using
% isentropic process equation (P1*v1^gamma = P2*v2^gamma = constant)
P_compression = constant_c./(v_compression.^gamma);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Process 2-3 is Constant Volume Heat Addition.
% (P2*v2)/T2 = (P3*v3)/T3 and v2 = v3
P3 = (P2*T3)/T2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Process 3-4 is isentropic expansion.
% P3*v3^gamma = P4*v4^gamma
% P4 = P3*(v3/v4)^gamma = P3*(1/cr_ratio)^gamma
P4 = P3*(1/cr_ratio)^gamma;
% Between end states 1 and 2, combined gas law equation is used
% (P1*v1)/T1 = (P2*v2)/T2
T4 = (P4*v4*T3)/(P3*v3);

% Since process 3-4 is an isentropic process, we cannot join directly end
% states 3 and 4 in the plot which gives straight line which is wrong.
% Hence we use formula which determines how volume inside the engine combustion
% chamber changes wrt to crank angle, cylinder dimensions & compression ratio
constant_c = P3*v3^gamma;
v_expansion = engine_kinematics(bore,stroke,con_rod,cr_ratio,0,180);
% Calculation of pressure at each point of v_expansion value using
% isentropic process equation (P3*v3^gamma = P4*v4^gamma = constant)
P_expansion = constant_c./(v_expansion.^gamma);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Process 4-1 is Constant Volume Heat Rejection.
% (P1*v1)/T1 = (P4*v4)/T4 and v1 = v4
% Variables at State 1 i.e P1, v1, T1 are already known.

%%%%%%%%%%%%%%%%%%%%%%%%%%% THERMAL EFFICIENCY %%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Thermal efficiency of Otto cycle depends on only compression ratio
Thermal_Efficiency = 1-(1/(cr_ratio^(gamma-1)))

% Ploting
figure(1)
hold on
% plot(v1,P1,'*','color','r')
% plot(v2,P2,'*','color','r')
% plot(v3,P3,'*','color','r')
% plot(v4,P4,'*','color','r')

plot(v_compression, P_compression, 'linewidth', 2, 'color', 'r')
plot([v2 v3], [P2 P3], 'linewidth', 2, 'color', 'b')
plot(v_expansion, P_expansion, 'linewidth', 2, 'color', 'g')
plot([v4 v1], [P4 P1], 'linewidth', 2, 'color', 'k')

axis([0 (v1+5e-05) -5e05 (P3+5e05)])
title('P-v plot for compression ratio of 6')
xlabel('Volume')
ylabel('Pressure')
grid on
legend('Isentropic Compression','Constant Volume Heat Addition' , 'Isentropic Expansion', 'Constant Volume Heat Removal')

function [V] = engine_kinematics(bore,stroke,con_rod,cr_ratio,start_crank,end_crank)
    a = stroke/2;
    R = con_rod/a;

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

    theta = linspace(start_crank,end_crank,200);
    % Relation between crank angle and volume of the combustion chamber
    % V = v_c * (1 + 0.5*((cr_ratio - 1)*(R + 1 - cos(theta) - (R^2 - sin(theta)^2)^0.5)
    term1 = 0.5*(cr_ratio - 1);
    term2 = R + 1 - cosd(theta);
    term3 = (R^2 - sind(theta).^2).^0.5;

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

OUTPUT:

Thermal efficiency of Otto cycle depends on compression ratio only

For given above inputs, the value of thermal efficiency for Otto cycle is 

Compression ratio                  Thermal efficiency

         4                                         0.4257

         6                                         0.5166

         8                                         0.5647

It can be observed that as the value of compression ratio inceases, the thermal efficiency of engine also increases.

P-v plot of Otto cycle:

ERRORS:

Made mistake in entering the formula for volume trace during piston movement. Because of thet the final output is shown during plotting without proper completion of the compression and expansion processes.