Ploting the PV diagram of an Otto Cycle and also Calculating the thermal efficiency

%Program to trace the otto cycle of a petrol engine and also calculate the
%thermal efficiency

clear all
close all
clc

%Inputs given 
gamma=1.4;
t3=2300;

%State variables (1)
p1=101325;
t1=500;

%Engine Geometric Parameters
bore=0.1;                             %Bore diameter
stroke=0.1;                           %Stroke length
con_rod=0.15;                         %Connecting Rod
cr=12;                                %Compression Ratio

%Calculating the swept volume and the clearance volume
v_swept=(pi/4)*bore^2*stroke;
v_clearance=v_swept/(cr-1);

%As summing of the swept volume and the clearance volume we can get the V1 and the cleance volume is equal to V2 
v1=v_swept + v_clearance;
v2=v_clearance;

%State variables at state point 2
%P2V2^gamma =P1V1^gamma
%P2=P1*(V1/V2)^gamma=P1*CR^gamma

 p2= p1*cr.^gamma;

%P1V1/T1= P2v2/T1 | T2=P2*V2*T1/(P1*V1)
 t2=p2*v2*t1/(p1*v1);
 
%Calculating the value of constant using PV^gamma=C equation 
 constant_c= p1*(v1^gamma);
 
%Calculating the Compression pressure and volume using the engine
%kinematics function
 v_compression= engine_kinematics(bore,stroke,con_rod,cr,180,0);
 p_compression=(constant_c./(v_compression.^gamma));


%State Variables at State Point(3)
 v3=v2;
 
%P3V3/T3=P2V2/T2 | P3=P2*T3/T2
  p3=p2*t3/t2;
  
%Calculating the value of constant using PV^gamma=C equation 
  constant_c= p3*(v3^gamma);
  
%Calculating the Expansion pressure and volume using the engine
%kinematics function
  v_expansion= engine_kinematics(bore,stroke,con_rod,cr,0,180);
  p_expansion= constant_c./v_expansion.^gamma;

 %State Variables at State Point(4)
 v4=v1;

 
 %P3V3^gamma=P4V4^gamma | P4=P3(V3/V4)^gamma
 p4=p3.*(v3/v4).^gamma;
 t4=p4*t1/p1;
 
 %To trace the isochoric curves taking a spacing of 1000 variables and
 %using the linspace command
 sp = 1000;         %Spacing variable
 %For state point 2 and 3
 y= linspace(p2, p3, sp);       
 %For state point 1 and 4
 y11 = linspace(p1, p4, sp);
 
 %Plotting the PV diagram
 figure
 hold on
 plot(v2,y,'.','color','g')                      %Tracing the isochoric curve (1)
 plot(v1,y11,'.','color','b')                    %Isochoric curve (2)
 plot(v1,p1,'*','color','r')                     %State point 1
 up = plot(v_compression,p_compression,'k');     %Tracing the compression curve
 up(1).LineWidth = 2;                            %Setting the width of the curve
 plot(v2,p2,'*','color','r')                     %State point 2
 plot(v3,p3,'*','color','r')                     %State point 3
 down = plot(v_expansion,p_expansion,'r');       %Tracing the expansion curve
 down(1).LineWidth = 2;                          %Setting the width of the curve
 plot(v4,p4,'*','color','r')                     %State point 4
 
%Labelling and Naming the curve
 xlabel('volume')
 ylabel('Pressure')
 legend('P-V Diagram')

%Calculating the thermal efficiency of the otto cycle using the thermal function 
 eff= thermal(v1,v2,gamma);
 disp('\n The thermal efficiency is:-');
 disp(eff)
 
 
%Function to calculate the thermal efficiency 

function [t]= thermal (vol1,vol2,gammaN)

%Taking the t variable for storing the thermal efficiency
t= 1-(vol2/vol1)^(gammaN-1); 

end

%Function to calculate the compression or expansion volume of the otto cycle.

function [v]=engine_kinematics(bore,stroke,con_rod,cr,start_crank,end_crank)

%Crank pin dimensions
a=stroke/2;                                
R=con_rod/a;

%Volume dimensions
v_s=(pi/4)*(bore.^2)*stroke;
v_c=v_s/(cr-1);

%Calculating the theta values for a range between start_crank and End_crank
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;

%Compression or Expansion volume=
v=(1+ term1*(term2 - term3)).*v_c;


end

above is the code

PROCEDURE:

Here in this Matlab code, we have traced the PV diagram of the Otto cycle and also find out the thermal efficiency of the Otto cycle. First of all, we will take the engine geometric parameters bore diameter, stroke length, connecting rod dimensions and compression ratio as inputs. Then we will calculate the swept and clearance volume inside the engine cylinder using the formulas. The state variables of the state point 1 ie pressure, volume, and temperature are given and the state variables of the state point 2 can be calculated by the ideal gas law.

PV=nRT

The expansion and the compression curves of the Otto cycle can be traced by engine kinematics functions ie an user-defined function. In the engine kinematics function, we are solving an equation through which the compression volume can be calculated. The formula is given below:-

`V/(Vc)=(1+1/2(CR-1)[R+1-cos(theta)-(R^2-(sin(theta)^2)^(1/2)])`

The engine kinematics function will be taking all the engine geometric parameters and also the starting and ending crank angle. In that function, we will be calculating the crank pin dimensions and then using a linspace command we will be taking a range of 100 values of theta between the starting and ending crank angles and then implementing the formula by dividing it into 3 terms. As we know the ideal gas equation for compression and expansion:-

`PV^gamma=C`

The value of constant ‘c’ using the formula above. Then by calling the engine kinematics function we can calculate the compression volume and thereby pressure at compression time. Similarly, we can calculate the pressure and volume during the expansion time. And also the state variables at 3 and 4 can find out by using the ideal gas equation. Then we will plot the state variables in the plot and join them by using the suitable commands. In order to trace the expansion and compression curves, we can use the values of the pressure and expansion at that points in the plot command. In order to join the isochoric curves, we will use the linspace command to get a range of 1000 values between the respective pressure points. And we will plot in by fixing the x-axis parameter to the particular value of the volume and the y-axis parameter as the array of 1000 value i.e. used calculated earlier. And the width of the expansion and compression curve can be fixed by linewidth command. We will label the x-axis and y-axis of the PV diagram as ’Volume’ and ‘Pressure’. For calculating the thermal efficiency we have taken a separate function ‘thermal’ and we will pass the volume at state points 1 and 2 and the gamma value. The thermal efficiency is calculated and output is given.

ERROR:-

Here, the error occurs due to matrix dimensions not satisfying. it can be overcome by putting './ ' after the constant_c variable in the p_expansion and p_compression line because v_compression is of 1*100 double array so that p_compression must satisfy.

it is a logical error at deciding the value of the spacing variable('sp') because a large value to the variable will consume more time in running the code and also more crowded dots in the isochoric curve. Hence looking at these problems, sp=1000 is taken so that it can overcome all those problems.

OUTPUT:-