MATLAB PROGRAM: TO CREATE AN OTTO CYCLE PV DIAGRAM USING PISTON KINEMATICS AND ITS EFFECIENCY

To write a code to generate the pv diagram and to calculate the effeciency of otto cycle using matlab:


Otto cycle:
1.Otto cycle is a thermodynamic cycle which is most commonly used in automobile engines.
2.It describes of what happends to the mass of gas when it is subjected to change of pressure,volume,temperature ,addition of heat and removal of heat.
3.The mass of gas that is subjected to those changes are called system.
4.In case of otto cycle the effect will be to produce enough net work from the system so as to propel an automobile along the environment.

Various process in otto cycle:

The various processes in Otto cycle are :

process 1-2: Isentropic compression process
It is an adiabatic compression process in which the piston moves from BDC to TDC
process 2-3: constant volume heat addition process
combustion stroke in which the heat gets transfered form external source to fluid.
process 3-4:Isentropic expansion process
It is an adiabatic expansion process or power stroke.
process 4-1:constant volume heat reduction
The heat is rejected to the idealized external sink (air).

Given:
Geometric parameters of engine with assigned values:
Bore
Compression ratio
Stroke
Connecting rod length &
initial state variable like t1,t3,p1

% Piston kinematic equation tells us how piston moves from TDC to BDC and how combustion chamber volume changes as a function of crank angle.

PROGRAM CODES WITH EXPLANATIONS:

clear all
close all
clc

%program to stimulate pv diagram and thermal effeciency of otto cycle;

%given inputs
gamma=1.4;
t3=2300

%required variables at state 1:
 p1=101325
 t1=500

%engine geometric parameters:
 bore=0.1;
 stroke=0.1;
 con_rod=0.15;
 cr=12;
 
%inorder to get v1 and v2 we need swept_volume & clearance_volume
%to calculate swept volume and clearance volume

 v_swept=(pi/4)*stroke*bore^2;
 v_clearance= v_swept/(cr-1);
 
 v1=v_swept + v_clearance
 
%v2 is nothing but volume at TDC is our clearance_volume

 v2= v_clearance
 
% Required variables at state 2:
% With isentropic compression process which represents
% p1v1^gamma=p2v2^gamma inorder to find p2?

% p2= p1*(v1/v2)^gamma or p2=p1*cr^gamma
  p2=p1*cr^gamma
  

  
  constant_c=p1*v1^gamma;
% During compression stroke start_crank is 180& end_crank is 0;
 V_compression=eng(bore,stroke,con_rod,cr,180,0);
 p_compression=constant_c./V_compression.^gamma

% To compute t2 we use ideal gas law, since pv = constant
% p1v1/t1=p2v2/t2
% t2=p2*v2*t2/p1*v1

  t2=p2*v2*t1/(p1*v1);
  
% variables at state point 3:
% state pont 3 represents the end of const volume combustion process
  v3=v2;
% p3v3/t3=p2v2/t2
% p3=p2t3/t2 & we have assumed the maximum temp to be 2300k
  p3=p2*(t3/t2);
  
   constant_c=p3*v3^gamma;
% During expansion stroke start_crank is 0 & end_crank is 180;
 V_expansion=eng(bore,stroke,con_rod,cr,0,180);
 p_expansion=constant_c./V_expansion.^gamma
  
 
% state variable at point 4:
  v4=v1
% p3v3^gamma = p4v4^gamma
% p4=p3(v3/v4)^gamma

  p4=p3*(v3/v4)^gamma
  
%p4*t4/t4=p1*v1/t1
%t4=(p4/p1)*t1
  t4=(p4*t1/p1);
  
%               (2)To find the effeciency of an otto cycle
%                  Effeciency=1-((t4-t1)/(t3-t2))
      Effeciency=1-((t4-t1)/(t3-t2))
      
  figure (1)
  hold on
  plot(v1,p1,'*','color','g')
 plot(V_compression,p_compression,'linewidth',2,'color','r')
  plot(v2,p2,'*','color','g')
  plot(v3,p3,'*','color','g')
  plot(V_expansion,p_expansion,'linewidth',2,'color','b')
  plot(v4,p4,'*','color','g')
  xlabel('volume')
  ylabel('pressure')
  
line ([v2 v3],[p2 p3])
line([v4 v1],[p4 p1])

Once we execute the corresponding program it will generate a graph connecting our state points in a straight line .

Inorder to represent it with the curved line we need to call a function which contains the formula how the volume inside the combustion chamber changes can be derived from the geometric equation:



Calling a function in the main program:
PROGRAM CODE:

 function[V]=eng(bore, stroke, con_rod,cr,start_crank,end_crank)
 
 a=stroke/2;
% Crank pin radius R
 R=con_rod/a;
 
 v_s = (pi/4)*stroke*bore^2;
 v_c = v_s/(cr-1);
 
 % Formula to find how the volume inside the combustion chamber changes
 % derived from geometric relationship;
 % Piston kinematic equation tells us how piston moves from TDC to BDC &
 % how combustion chamber volume changes as a function 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_c
 
 end

output graph plot:

output values:

Effeciency :

We got the effeciency about: 63%