MATLAB PROGRAMMING FOR CALCULATING EFFICIENCY, ALL STATE POINTS AND GENERATION OF P - V DIAGRAM OF AN OTTO CYCLE

AIM : To calculate Efficiency and to plot a p - v diagram by calculating all the state point variables

Figure : 

1-2 suction

2-3 Isentropic compression

3-4 Isochoric or constant volume heat addition

4-5 Isentropic expansion

5-6 Isochoric or constant volume heat rejection

MATLAB Code :

For simplification in analysing the code it is divided in to 2 halfs. 

CODE 1 :

%calculation of volume variation wrt to crank angle
function volume=engine_kinematics(bore,stroke,con_rod,cr,start_crank,end_crank,n)
%crank pin radius
a=(stroke/2);
%Ratio of connecting rod to crank pin radius
R=(con_rod/a);
%swept volume
v_s=(pi/4)*bore^2*stroke;
%clearance volume
v_c=v_s/(cr-1);
%crank angle
theta=linspace(start_crank,end_crank,n);
term1=0.5*(cr-1);
term2=R+1-cosd(theta);
term3=(R^2-sind(theta).^2).^0.5;
%volume variation wrt to crank angle
volume=(1+term1*(term2-term3)).*v_c;
end

 CODE 2 :

%matlab program to calculate all state variables and to plot  p-v diagram of otto cycle
clear all
close all
clc
%Engine bore diameter in (m) 
bore=0.1;
%stroke length in (m)
stroke=0.1;
%compression ratio
cr=8;
%Length of connecting rod in (m)
con_rod=0.12;
%Ratio of specific heat capacities for air
gamma=1.4;
%No of crank angle divisions
n=180;
%Maximum temperature in cylinder chamber(K)
t3=1600;
%swept volume in (m^3)
v_s=(pi/4)*bore^2*stroke;
%clearance volume in (m^3)
v_c=v_s/(cr-1);
%state 1 properties
%Volume at the start of compression in (m^3) 
v1=v_c+v_s;
%Temperature at the start of compression in (k)
t1=300;
%pressure at the start of compression in (KPa)
p1=101.325;
%state 2 properties
%Volume at the end of compression in (m^3)
v2=v_c;
%pressure at the end of compression in (Kpa)calculated by p*v1^gamma=p2*v2^gamma
p2=p1*(cr)^gamma;
%Temperature at the end of compression is calculated by t1*v1^(gamma-1)=t2*v2^(gamma-1)
t2=t1*(cr)^(gamma-1);
%volume curve generation by calling the Engine kinematics function
V_compression=engine_kinematics(bore,stroke,con_rod,cr,180,0,n);
constant=p1*v1^gamma;
p_compression=constant./(V_compression.^gamma);
%state 3 properties
%volume before expansion
v3=v2;
%maximum pressure after heat addition
p3=(p2/t2)*t3;
%volume curve generation by calling the Engine kinematics function
V_expansion=engine_kinematics(bore,stroke,con_rod,cr,0,180,180);
constant_1=p3*v3^gamma;
p_expansion=(constant_1)./(V_expansion.^gamma);
%state 4 properties
%volume after expansion
v4=v1;
%pressure after expansion
p4=p3*(1/cr)^gamma;
hold on
axis([0 0.001 0 4500])
%plotting pv diagram
plot(v1,p1,'o','color','g','linewidth',3)
plot(V_compression,p_compression,'linewidth',2,'color','b')
plot(v2,p2,'o','color','g','linewidth',3)
plot(v3,p3,'o','color','g','linewidth',3)
plot(V_expansion,p_expansion,'linewidth',2,'color','b')
plot(v4,p4,'o','color','g','linewidth',3)
line([v2 v3],[p2 p3],'linewidth',2,'color','b')
line([v4 v1],[p4 p1],'linewidth',2,'color','b')
xlabel('volume(m^3)')
ylabel('pressure(kpa)')
Efficiency=(1-(1/cr^(gamma-1)))*100;
disp(['Efficiency of otto cycle in % =',num2str(Efficiency)])

Code Explanation :

In code 1 a function is taken and relation for volume in terms of crank angle is used and the output volume path is calculated

In code 2 the initial volume is calculated by bore,stroke length and all the input state variables are taken and later all the properties of pressure ,volume and temperature are calculated at each and every state point. The function which is written in code 1 is called in code 2 for volume curve generation.Later required plots and efficiencies are calculated.

Single complex code:

%matlab program to calculate all state variables and to plot  p-v diagram of otto cycle
clear all
close all
clc
%Engine bore diameter in (m) 
bore=0.1;
%stroke length in (m)
stroke=0.1;
%compression ratio
cr=8;
%Length of connecting rod in (m)
con_rod=0.12;
%Ratio of specific heat capacities for air
gamma=1.4;
%No of crank angle divisions
n=180;
%Maximum temperature in cylinder chamber(K)
t3=1600;
%swept volume in (m^3)
v_s=(pi/4)*bore^2*stroke;
%clearance volume in (m^3)
v_c=v_s/(cr-1);
%state 1 properties
%Volume at the start of compression in (m^3) 
v1=v_c+v_s;
%Temperature at the start of compression in (k)
t1=300;
%pressure at the start of compression in (KPa)
p1=101.325;
%state 2 properties
%Volume at the end of compression in (m^3)
v2=v_c;
%pressure at the end of compression in (Kpa)calculated by p*v1^gamma=p2*v2^gamma
p2=p1*(cr)^gamma;
%Temperature at the end of compression is calculated by t1*v1^(gamma-1)=t2*v2^(gamma-1)
t2=t1*(cr)^(gamma-1);
%volume curve generation by calling the Engine kinematics function
V_compression=engine_kinematics(bore,stroke,con_rod,cr,180,0,n);
constant=p1*v1^gamma;
p_compression=constant./(V_compression.^gamma);
%state 3 properties
%volume before expansion
v3=v2;
%maximum pressure after heat addition
p3=(p2/t2)*t3;
%volume curve generation by calling the Engine kinematics function
V_expansion=engine_kinematics(bore,stroke,con_rod,cr,0,180,180);
constant_1=p3*v3^gamma;
p_expansion=(constant_1)./(V_expansion.^gamma);
%state 4 properties
%volume after expansion
v4=v1;
%pressure after expansion
p4=p3*(1/cr)^gamma;
hold on
axis([0 0.001 0 4500])
%plotting pv diagram
plot(v1,p1,'o','color','g','linewidth',3)
plot(V_compression,p_compression,'linewidth',2,'color','b')
plot(v2,p2,'o','color','g','linewidth',3)
plot(v3,p3,'o','color','g','linewidth',3)
plot(V_expansion,p_expansion,'linewidth',2,'color','b')
plot(v4,p4,'o','color','g','linewidth',3)
line([v2 v3],[p2 p3],'linewidth',2,'color','b')
line([v4 v1],[p4 p1],'linewidth',2,'color','b')
xlabel('volume(m^3)')
ylabel('pressure(kpa)')
Efficiency=(1-(1/cr^(gamma-1)))*100;
disp(['Efficiency of otto cycle in % =',num2str(Efficiency)])
%calculation of volume variation wrt to crank angle
function volume=engine_kinematics(bore,stroke,con_rod,cr,start_crank,end_crank,n)
%crank pin radius
a=(stroke/2);
%Ratio of connecting rod to crank pin radius
R=(con_rod/a);
%swept volume
v_s=(pi/4)*bore^2*stroke;
%clearance volume
v_c=v_s/(cr-1);
%crank angle
theta=linspace(start_crank,end_crank,n);
term1=0.5*(cr-1);
term2=R+1-cosd(theta);
term3=(R^2-sind(theta).^2).^0.5;
%volume variation wrt to crank angle
volume=(1+term1*(term2-term3)).*v_c;
end

Output :

Error : No specific errors are generated while executing the code

Note :May be we should look after the axes limits for better view of p-v diagram and V_compression and V_expansion are matrices so . operations should be used in calculations.