Air standard Cycle (MATLAB)

AIM- A MATLAB code for Air Standard Cycle(otto cycle):-

The otto cycle approximates the processes in petrol or diesel engines. At constant volume, Heat is transferred to the working fluid. It consists of four internally reversible processes

1. Isentropic compression.
2. Constant volume heat addition.
3. Isentropic expansion.
4. Constant volume of heat rejection.

Main Code:-

• Initially provide the all required values (bore, stroke, connecting rod, and compression ratio etc.). 
• The values of swept volume and clearance volume is calculated with the following formula-

            $v_s=\frac{\pi }{4}\cdot b\text{or}{e}^2\cdot stroke$   and

            $v_c=\frac{v_s}{cr-1}$

        For the isentropic compression, we know 

            $p2\cdot v{2}^{\gamma }=p1\cdot v{1}^{\gamma }$

            $p2=p1\cdot {\left(\frac{v1}{v2}\right)}^{\gamma }$

              so  $p2=p1\cdot c{r}^{\gamma }$

         For an ideal gas  $\frac{p1\cdot v1}{t1}=\frac{p2\cdot v2}{t2}$ 

             so $t2=\frac{p2\cdot v2\cdot t1}{p1\cdot v1}$

• As the values of theta and other variables vary in a function of           

            $p\cdot v^{\gamma }=cons\tan t$

• due to the adiabatic process, a function 'engine_keinematics' is created to generate a non-linear p-v curve and also get the value of V_compression, and to calculate P_compression we use the above equation.       
• At state point-3 the volume is constant $(v3=v2)$

            $\frac{p3v3}{t3}=\frac{p2v2}{t3}$ 

       so  $p3=\frac{p2\cdot t3}{t2}$

• During the expansion process, the curve will follow the expression($v4=v1$)

           $p3\cdot v{3}^{\gamma }=p4\cdot v{4}^{\gamma }$

           $p4=p3\cdot {\left(\frac{v3}{v4}\right)}^{\gamma }$

• Similarly, by using the function 'engine_keinematics' we can get the value of V_expansion  and P_expansion.
• Finally use the plot() command and plot all the curves. 

% creating a PV diagram otto cycle


  clear all 
  close all
  clc
  
% given values
  bore= 1;
  stroke= 1;
  connecting_rod= 1.5; 
  gamma=1.5;
% compression ratio 
  cr = 5;
% swept volume
  v_swept= pi/4*bore^2*stroke;
% clearance volume
  v_clearance= v_swept/(cr-1);
% at state point 1
  p1 = 1000;
  t1 = 500;
  v1 = v_clearance + v_swept;
  v2 = v_clearance;
  
% Compression
  
  p2 = p1*cr^gamma;
  t2 = p2*v2*t1/p1*v1;
   V_compression = engine_keinematics(bore, stroke, connecting_rod, cr, 180 , 0);
   constant_c = p1* v1 ^gamma;
   P_compression = constant_c./V_compression.^gamma;
   
   
% at state point 2-3 
  v3=v2;
  t3 = 2000;
  p3= p2*t3/t2;
  
% Expansion  
   constant_c = p3* v3 ^gamma;
   V_expansion = engine_keinematics(bore, stroke, connecting_rod, cr, 0 ,180);
   P_expansion = constant_c./V_expansion.^gamma;
   
  v4 = v1;
  p4 = p3*(v3/v4)^gamma;
  t4 = p4*v4*t3/p3*v3;
  
  
  hold on
   plot(v1,p1, '*' , 'color', 'g' )
  plot( V_compression ,P_compression)
   plot(v2,p2, '*' , 'color', 'g' )
   plot(v3,p3, '*' , 'color', 'g' )
  plot( V_expansion   ,P_expansion ) 
   plot(v4,p4, '*' , 'color', 'g' )
   plot([v2 v3],[p2 p3], 'color', 'b')
   plot([v4 v1],[p4 p1], 'color', 'b')
   xlabel('VOLUME')
   ylabel('PRESSURE')
   
   
% Thermal Efficency
 Th_eff = (1-(1/cr^(gamma-1)))* 100;

Code for function calling:

% calling a function
function [V] = engine_keinematics(bore, stroke, connecting_rod, cr, start_cranc, end_cranc)

   a = stroke/2;
   r = connecting_rod/a;
   v_swept= pi/4*bore^2*stroke;
   v_clearance= v_swept/(cr-1);
   
   theta= linspace(start_cranc,end_cranc,180);
   
   term1 = 0.5*(cr-1);
   term2 = r + 1 - cosd(theta);
   term3 = (r^2 - sind(theta)).^0.5;
  
   V = (1 + term1*(term2 - term3)).* v_clearance ;
    
 end

Results-

Thermal efficiency = 55.27864

Error:-

Wrong closing of brackets in the following expression

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