Air standard cycle

AIM- To write a Matlab code to find out thermal efficiency of an otto cycle and plot a P-V graph .

THEORY-

fig 1 P-V plot of an otto cycle

PROCESS DESCRIPTION-

1.PROCESS 1-2 - Adiabatic or Isentropic Compression

(which follows the equation `PV^(gamma)=C`)

2.PROCESS 2-3 - Constant volume heat addition

(which follows the equation `(PV)/T=C`)

3.PROCESS 3-4 - Adiabatic or Isentropic Expansion

(which follows the equation `PV^(gamma)=C`)

4.PROCESS 4-1 - Constant volume heat rejection

(which follows the equation `(PV)/T=C`)

Compression ratio = volume before compression / volume after compression

= `=V_1/V_2`

`= (V_s+V_c)/V_C`

`V_s`= Swept volume

`V_c`= Clearance volume

ENGINE PARAMETERS-

bore (bore diameter)=0.1m
stroke=0.1m
con_rod(connecting rod length)=0.15
comp_ratio (compression ratio)=12

INPUT PARAMETRES-

gamma(`gamma`)=1.4

t1(temperature at state point 1)= 500k

t3(temperature at state point 3) = 200k

p1(preassure at state point 1)=101325 pascals

Software used- MatlabR2018a

INPUT PROGRAM-

clear all 
close all
clc
% inputs
gamma=1.4
 t3=2300

% state variables 
p1=101325
t1=500
% parameters of engine 
bore=0.1;
stroke=0.1;
con_rod=0.15;
comp_ratio=12;

% swept volume and clearence volume caliculation
v_swept=(pi/4)*bore^2*stroke;
v_clearance=(v_swept)/(comp_ratio-1)

% state points at 1
v1=v_swept+v_clearance


% state variables at point 2
% p1v1^gamma=p2v2^gamma
% p2= (p1v1/v2)^gamma = p1*comp_ratio^gamma
p2=p1*comp_ratio^gamma;
v2=v_clearance

% p1v1/t1=p2v2/t2 \ t2=(p2*v2*t1)/(p1*v1)
 t2=(p2*v2*t1)/(p1*v1);
 constant_c=p1*v1^gamma;
 V_compression=engine_kinematics(bore,stroke,con_rod,comp_ratio,180,0)
 P_compression=constant_c./V_compression.^gamma;

 % state variables at point 3
 v3=v2;
 % p3v3/t3=p2v2/t2 \ p3=p2t3/t2
 p3=p2*t3/t2;
 
 constant_c=p3*v3^gamma;
 V_expansion=engine_kinematics(bore,stroke,con_rod,comp_ratio,0,180)
 P_expansion=constant_c./V_expansion.^gamma;

  % state variables at point 4
  v4=v1
  % p3v3^gamma=p4v4^gamma \ p4=p3(v3/v4)^gamma
  p4=p3*(v3/v4)^gamma;
  efficiency=1-1/(comp_ratio^(gamma-1))
  
  figure(1)
  hold on
  plot(v1,p1,'*','color','b')
  plot(V_compression,P_compression ,'color','g')
  plot(v2,p2,'*','color','b')
  plot([v2,v3],[p2,p3] ,'color' ,'r')
  plot(v3,p3,'*' ,'color' ,'b')
  plot(V_expansion,P_expansion,'color','g')
  plot(v4,p4,'*' ,'color' ,'b')
  plot([v4,v1],[p4,p1],'color','r')
  grid on
 xlabel('volume')
 ylabel('pressure')
 
  figure(2)
  plot([v1 v2 v3 v4 v1],[p1 p2 p3 p4 p1] ,'color','b')
 xlabel('volume')
 ylabel('pressure')
  grid on

MATLAB CODE EXPLANATION-

//

clear all 
close all
clc
% inputs
gamma=1.4
 t3=2300

% state variables 
p1=101325
t1=500
% parameters of engine 
bore=0.1;
stroke=0.1;
con_rod=0.15;
comp_ratio=12;

Initially parameters such as at state point 1,t3 and engine parameters are assigned as mentioned above.

//

% swept volume and clearence volume caliculation
v_swept=(pi/4)*bore^2*stroke;
v_clearance=(v_swept)/(comp_ratio-1)

Further now swept volume and clearence volume are caliculated using the above formula and assigned to arrays v_swept and v_clearance.

//

% state points at 1
v1=v_swept+v_clearance

% state variables at point 2
% p1v1^gamma=p2v2^gamma
% p2= (p1v1/v2)^gamma = p1*comp_ratio^gamma
p2=p1*comp_ratio^gamma;
v2=v_clearance
% p1v1/t1=p2v2/t2 \ t2=(p2*v2*t1)/(p1*v1)
 t2=(p2*v2*t1)/(p1*v1)

Now volume at state point 1 is calculated using the formula v1=v_swept+v_clearance.

Since from state point 1 to state point 2 the process is Isentropic compression it follows the equation `p_1v_1^gamma=p_2v_2^gamma`,p2 ,t2 and v2 are caliculated using the above mentioned formulas.

//

constant_c=p3*v3^gamma;
 V_expansion=engine_kinematics(bore,stroke,con_rod,comp_ratio,0,180)
 P_expansion=constant_c./V_expansion.^gamma;

We know that from state point 1 to state point 2 isentropic expansion occurs , if we join these two state points a linear curve occurs but it is wrong, inorder to obtain a realistic curve we use engine kinematics equation and obtain volume during expansion stroke.

The evaluation of engine kinematics is done in a function program named engine_kinematics and the obtained results are stored in an array V_expansion.

//

Function program explanation-

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


comp_ratio=12;
a=stroke/2;
R=con_rod/a;
v_swept=(pi/4)*bore^2*stroke;
v_clearance=v_swept/(comp_ratio-1);

theta=linspace(start_crank,end_crank ,100);

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

V=(1+term1*(term2-term3))*v_clearance
end 

In order to obtain volume we have an equation from engine kinematics i.e `V/V_c= 1+1/2(r_c-1)[R+1-costheta-(R^2-sin^2theta)^(1/2)]`

since the equation is pretty long we have divided it into 3 terms as term1,term2,term3 and arranged them to obtain volume[V].

For the above mentioned function program input values of bore,stroke,con_rod values are taken from the main program and the values of start_crank and end_crank values are taken from main program of respective compression and expansion processes.

Between these starting and end crank 100 values are taken.

v_clearence is also calculated using the formula mentioned in the program.

Now all the values of volume are stored in array V.

//

% state variables at point 3
 v3=v2;
 % p3v3/t3=p2v2/t2 \ p3=p2t3/t2
 p3=p2*t3/t2;
 
 constant_c=p3*v3^gamma;
 V_expansion=engine_kinematics(bore,stroke,con_rod,comp_ratio,0,180)
 P_expansion=constant_c./V_expansion.^gamma;

Now at at state point 3 the parameters are obtained using the above mentioned formulas in the program.

Since the process from 3 to 4 is expansion the values are caliculated for volume and pressure as V_expansion and P_expansion,the crank angles are given as 0 to 180 since the crank moves from TDC to BDC.

//

% state variables at point 4
  v4=v1
  % p3v3^gamma=p4v4^gamma \ p4=p3(v3/v4)^gamma
  p4=p3*(v3/v4)^gamma;
  efficiency=1-1/(comp_ratio^(gamma-1))

The values of pressure and volume are obtained at state poin 4.

Then the efficiency is caliculated using the above mentioned formula in the program.

//

 figure(1)
  hold on
  plot(v1,p1,'*','color','b')
  plot(V_compression,P_compression ,'color','g')
  plot(v2,p2,'*','color','b')
  plot([v2,v3],[p2,p3] ,'color' ,'r')
  plot(v3,p3,'*' ,'color' ,'b')
  plot(V_expansion,P_expansion,'color','g')
  plot(v4,p4,'*' ,'color' ,'b')
  plot([v4,v1],[p4,p1],'color','r')
  grid on
 xlabel('volume')
 ylabel('pressure')

Now a plot labeld figure(1) is created , then obtained values of v1,p1, v2,p2, v3,p3, v4,p4 are plotted as points coloured in blue usin the plot command.

Later the process 1 to 2 are joined using plot command for the values of V_compression,P_compression.

Since for the process 2 to 3 it varies linearly the points are joined using plot command.

Now the same process is repeated for point 3 to 4 an4 to 1.

X-axis label and Y-axis label are mentioned using the command xlabel and ylabel, using grin on command the grid of the plot is obtained.

//

  figure(2)
  plot([v1 v2 v3 v4 v1],[p1 p2 p3 p4 p1] ,'color','b')
 xlabel('volume')
 ylabel('pressure')
  grid on

Now a plot named finure(2) is created , this plot is just to show how the plot is obtained if we just join all the points obtained at the end of each and every process using plot command without using engine kinematics concept, the obtained plot would be represented in a wrong manner hence we have used the concept of engine kinematics to determine a correct looking plot.

OUTPUT- The following are the outputs obtained after successful execution of the Matlab code.

fig 2 obtained PV plot of an Otto cycle

fig 3 wrong PV plot

fig 4 showing efficiency obtained

CONCLUSION- Hence efficiency of 62% and the PV plots of Otto cycle are determined successfully.