Plotting Otto cycle using MatLab

OTTO CYCLE

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.

(Source: https://en.wikipedia.org/wiki/Otto_cycle)

Generally, an otto cycle is a typical graph used to visualise the change in pressure, temperature, volume during the combustion process of fuel(Petrol) in a SI engine.

An ideal otto cycle consists of mainly four thermodynamic processes :

• 1 → 2: Isentropic Compression process.
• 2 → 3: Constant Volume heat addition process.
• 3 → 4: Isentropic Expansion process.
• 4 → 1: Constant volume heat rejection process.

Otto Cycle

(Source: https://en.wikipedia.orghttps://en.wikipedia.org/wiki/Otto_cycle#/media/File:P-V_Otto_cycle.svg)

Otto cycle helps to understand how the heat energy present in petrol is converted into motion.

Each process follows different thermal processes, so the way of approach varies for each process.

1 → 2: Isentropic Compression process:

In this process, the air is compressed to a particular pressure isentropically, where the change in entropy is zero i.e., there will be no transfer of heat and matter.

For most of the problems, initial conditions provided and either max.pressure or temperature is provided or while experimenting we can measure the pressure and temperature of the air going into the cylinder and maximum temperature can be determined by the material properties of metal used as engine cylinder lining and lubricant used.

In this stage,

Initially, the piston will be at BDC(Bottom Dead Centre). So it is the sum of swept volume (`V_s`) by the piston (BDC to TDC or vice versa) and clearance volume (`V_c`).

The Clearance volume(`V_c`) is provided to accommodate the compressed volume of air and fuel injected into the chamber.

`V_s = pi/4 *d^2 *stroke` .

d = Bore diameter.

stroke = Stroke length of the piston.

Compression ratio

`C_r = (V_1 / V_c)`.

where volume at state-1 is,

`V_1= V_c + V_s`

`C_r = (V_c+V_s)/V_c`

`=> C_r = 1+(V_s/V_c)`

`therefore` `V_c = V_s/(Cr-1)`

Let `P_1, V_1, T_1` be pressure, volume and temperature at state-1and similarly `P_2, V_2, T_2` are at state-2, `P_3, V_3, T_3` for state-3 and `P_4, V_4, T_4` for state-4.

As it is an Isentropic process, the process of compression of air follows the below gas equation.

`p*v^gamma=Const`

i.e., `p_1* v_1^gamma = p_2*v_2^gamma`.

where `gamma` = 1.4 for air medium

Here `P_1, V_1, T_1`, (`V_2 = V_c`) values are known.

To calculate pressure `P_2`

`P_1* V_1^gamma = P_2*V_2^gamma`.

`=>P_2 = P_1 * (V_1/V_2)^gamma`

`therefore P_2 = P_1 * (C_r)^gamma`→ (.1.)

To calculate the temperature `T_2`

By applying the ideal gas equation at state-1,

`P_1* V_1 = M* R* T_1`

Here M = mass of intake air (assumed as constant.)

R = Gas constant.

Upon rearranging the above equation,

`P_1 = (M* R *T_1) / V_1`

Similarly at state-2,

`P_2 = (M* R* T_2) / V_2`

Substituting P1 and P2 in equation 1.

`(M*R*T_2) / V_2 = ( (M*R*T_1) / V_1 ) * (C_r)^gamma`

`(M*R) * T_2 = (M*R) * (T_1*(V_2/V_1)) * C_r^gamma`

`implies` `T_2 = (T_1/C_r)*C_r^gamma` `[because (V_1/V_2) = C_r ]`

`therefore T_2 = T_1 * C_r ^ ( gamma-1)`

Therefore, `P_2` and `T_2` can be found out using the above equations.

2 → 3: Constant Volume heat addition process:

In this process, the piston is at TDC until state-3, where fuel is injected into the clearance volume(Vc).

The pressure is increased by maintaining volume constant.

During this process, gas follows the isochoric process i.e., constant volume process.

So the gas follows below gas equation,

`(P*V)/T = const`

Now the above equation can be written as

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

`P_2/T_2 = P_3/T_3` `[because V_2 = V_3]`

As we know `T_3`,` P_2`, `T_2 `

By rearranging,

`therefore P_3 = P_2* (T_3/T_2)`

By the end of this, we have all state variables at both states 2 and 3.

By the end of this process, the Spark plug in the engine will ignite the pressurized air and fuel mixture. Which leads to a sudden increase in temperature and pressure which is `T_3` and `P_3` respectively.

3 → 4: Isentropic Expansion process:

At the start of the process, the Piston is pushed towards BDC due to high-pressure accumulation in less volume by the isentropic expansion process.

This process is similar to `1->3`, but the exact reverse will happen i.e., instead of compression, expansion takes place.

As we know state variables at state-3 and volume at state-4 `(V_4 = V_1)`.

Now to calculate remaining state variables at state-4 is by

We can use, the relations derived during the process `1 ->2`

Here, we know the values of `T_3, P_3, V_3, V_4`.

To calculate the pressure at state-4 from the equation-1,

`P_4 = P_3 * (C_r)^gamma`

Similarly temperature,

`T_4 = T_3 * C_r^(gamma-1)`

Plotting the Otto cycle:

Now we have all state variables at all the states in the process.

To plot the Pressure Vs Volume graph, Plot (x,y) function can be used.

Which gives a graph as shown in the figure below,

The above graph shows, the state points at each state are connected by a straight line.

But processes `1->2 & 3->4` are isentropic process. where the curve equation is based on

`P*V^gamma = const`

`=> P = (const) /V^gamma`

`=> P prop (1/V^gamma)`

which is not a linear curve.

So to obtain proper curve, we need to calculate the pressure and the temperature at each point of both isentropic processes.

Now to calculate those points, we need to know the relationship between crank angle and volume or pressure.

Because as we discussed earlier piston moves during compression and expansion

Which is the below equation:

`(V__comp )/V_c = [1 + 0.5 * (C_r-1) * (R + 1 - cos(theta)- sqrt(R^2 - sind(theta)^2))]`

Here R= (connecting rod length)/(Crank pin radius).
Crankpin radius = stroke length/2.

By using this equation we will have an array of volume data for each specific theta subtended by the crank.

Then by using `[P_1*V_1^gamma = P__comp*(V__comp)^gamma]`, P_comp can be calculate

Now, v_comp and p_comp are arrays calculated from the above formula.

Similarly, we calculate the p_exp and v_exp.

As this equation might make the code uncoordinated, Using function helps.

Now using these values graph can be plotted using plot(x,y). as shown in fig. below,

Calculation of efficiency:

The efficiency of a system is the desired output by the desired input.

For the Otto cycle efficiency is

`eta` = Work done/Heat added.

Heat added `(Q_1)` = `M*C_v *(T_3 - T_2)` [`because 2->3` Constant volume heat addition process.]

Heat rejected `(Q_2)`= `M*C_v*(T_4 - T_1)` [`because 4->1` Constant volume heat addition process.]

Here `C_v` is volume constant.

`C_v = 0.718 (KJ) / (Kg K)`

Workdone (W) = Heat added - Heat rejected.

`W = Q_1 - Q_2`

`eta = (Q_1 - Q_2) / Q_1`

`= {((M*C_v *(T_3 - T_2)) - (M*C_v*(T_4 - T_1))) / (M*C_v *(T_3 - T_2))}`

`=( (T_3- T_2) - (T_4 - T_1) )/( (T_3 - T_2))`

`= 1 - ((T_4 - T_1)/(T_3 - T_2))`

`= 1 - {[T_1* ( (T_4/T_1)-1) ]/[T_2 * ((T_3/T_2)-1)]}` [`because` taking out `T_1` from numerator and `T_2 form denominator]

`= 1- {(T_1/T_2) * [((T_4/T_1)-1)/((T_3/T_2)-1)]}`

`=> 1 - {(1/C_r^(gamma-1)) * [1]}` `[because` ` (T_2/T_1)=C_r^(gamma-1) & (T_3/T_4) =C_r^(gamma-1)`

` => (T_2/T_1)=(T_3/T_4) `

` => (T_4/T_1) = (T_3/T_2) `

`]`

`therefore eta = 1 - (1/C_r^(gamma-1) )`

CODE:

Main program to Plot Otto cycle

close all
clear all
clc

%Constants
    gamma = 1.4; %cp/cv = 1.4 for air.

%compression ratio
    %cr = V1/V2 
    cr = 12;

%Maximum temperature of the cycle in Kelvin
    T3 = 2500;

% Dimensions of cylinder
    Br = 0.1;  %Bore diameter in m
    Stroke_L =  0.1;   %Stroke Length in m
    con_rod = 0.15; %length of connecting rod in m

%calculation of clerance volume and swept volume
    vs = pi/4*Br^2 * Stroke_L ;  %swept volume by pistion
    vc = vs/(cr-1);  %clearance volume 

%State variables at state 1 (AFTER SUCTION)
    p1 = 101325;  %pressure in KN/m^2
    T1 = 290 ;  %temperature in Kelvin
    V1 = vs + vc; % suction volume in m^3 (VMAX)

%State variables at state 2 (AFTER COMPRESSION)
    %p1*V1^g = p2*V2^gamma | p2 = (V1/V2)^gamma * p1 | cr = V1/V2
    p2 = p1 * cr^gamma;
    V2 = vc ;
    T2 = T1 * cr^(gamma-1);

%State variables at state 3 (AFTER COMBUSTION)
    % p3/p2 = T3/T2 | p3 = (T3/T2)*p2
    p3 = (T3/T2)* p2;
    V3 = V2;

%State variables at state 4 (AFTER EXPANSION)
    V4 = V1;
    % p3/p4 = (V4/V3)^gamma | p4 = p3 * (V2/V1)^gamma
    p4 = p3 * cr^(-gamma);
    % T3/T4 = (V4/V3)^(g-1) | T4 = T3 * (V2/V1)^(gamma-1)
    T4 = T3 * cr^(1-gamma);

% calculating variation of volume and pressure with respect to crank angle.
   
  %Variation of volume and pressure form 1 to 2 (COMPRESSION PROCESS-BDC to TDC)
    constant1 = p1 * V1^gamma;
    V_comp = Engine_kinematics(Br,Stroke_L,con_rod,cr,180,360);
    P_comp = p1 .* (V1./V_comp).^gamma;
   
  % For process 2 to 3 (HEAT ADDITION)
  %Volume is constant and pressure varies linearly.

  %Variation of volume and pressure form 3 to 4 (EXPANSION PROCESS-TDC to BDC)
    constant2 = p3 * V3^gamma;
    V_exp = Engine_kinematics(Br,Stroke_L,con_rod,cr,0,180);
    P_exp = p3 .* (V3./V_exp).^gamma;

   %For process 3 to 4(HEAT REJECTION)  
   %Volume is constant and pressure varies linearly.

%EFFCIENCY OF THE OTTO CYCLE
    
    n = 1 - (1/cr^(gamma-1)); %efficieny of Otto cycle (Derivation in report)
    
   
%plotting the OTTO CYCLE
    figure(1)
    hold on

    plot([V2 V3],[p2,p3],'linewidth',1,'color','r') %Const. volume heat addition
    plot([V4 V1],[p4,p1],'linewidth',1,'color','r') %Const. volume heat rejection

    plot(V_comp,P_comp,'linewidth',1,'color','b') %Iscentopic compression
    plot(V_exp,P_exp,'linewidth',1,'color','b')   %Isentropic expansion

    xlabel('Volume')
    ylabel('pressure')


    plot(V1,p1,'*','color','b')
    plot(V2,p2,'*','color','b')
    plot(V3,p3,'*','color','b')
    plot(V4,p4,'*','color','b')
    
    figure(2)
    plot([V1 V2 V3 V4],[p1 p2 p3 p4],'linewidth',2,'color','r') 

    xlabel('Volume')
    ylabel('pressure')

 
 

Program to calculate volume array at different position of crank.

function [V] = Engine_kinematics(Br,Stroke_L,con_rod,cr,start_crank,end_crank)




theta = linspace(start_crank,end_crank,180); %Crank angle

crank_pin = Stroke_L/2;  %Crank pin radius
R = con_rod/crank_pin;
 
%calculation of clerance volume and swept volume
vs = pi/4*Br^2 * Stroke_L ;  %swept volume by pistion
vc = vs/(cr-1);  %clearance volume

%calculation of volume in terms of crank angle
term1 = 0.5 * (cr-1);
term2 = R + 1 - cosd(theta);
term3 = (R^2 - sind(theta).^2).^0.5;

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

end