Matlab- Plot the Otto Cycle PV graph and find out thermal efficiency

Please find the attached code, Explanation,result and errors faced during programming.

Code

%To Create an Otto Air Standard Cycle Plot and to find its Thermal efficiency.
clear all
close all 
clc

% Inputs & Known values

p1=101325
t1=500
t3=2300
stroke=0.1;
bore=0.1;
con_rod=0.15;
cr=12;
gamma=1.4;

swept_volume=(pi/4)*(bore^2)*stroke;
clearance_volume=swept_volume/(cr-1);
v1=swept_volume+clearance_volume
v2=clearance_volume

% for isentropic compression p1*v1^gamma=p2*v2^gamma. 
% also cr=v1/v2

p2=p1*cr^gamma

%(p1v1)/t1=(p2v2)/t2. Therefore t2=(p2v2t1)/(p1v1)

t2=(p2*v2*t1)/(p1*v1)
c=p1*v1^gamma;

%Aaray of values to plot during compresssion
v_compression = engine_kinematics(bore,stroke,con_rod,cr,0,180);
p_compression= (c)./v_compression.^gamma;

%process 2-3 is isochoric process. Therefore v2=v3

v3=v2

% To find p3,p3/t3=p2/t2

p3=p2*t3/t2

% variables at state 4. V4-V1 is an isochoric variable,therefore

v4=v1

%Process 3-4 is isentropic in nature
%Therefore

p4=p3*(v3/v4)^gamma
c=p3*v3^gamma;

%Array of values to plot during expansion
v_expansion= engine_kinematics(bore,stroke,con_rod,cr,180,360);
p_expansion=(c)./v_expansion.^gamma;

%Plotting
figure (2)
hold on 
plot (v1,p1,'*+*','color','r')
plot(v_compression,p_compression)
plot([v2 v3],[p2 p3])

plot (v2,p2,'*+*','color','r')

plot (v3,p3,'*+*','color','r')
plot(v_expansion,p_expansion)

plot (v4,p4,'*+*','color','r')
plot([v4 v1],[p4 p1])

% for finding out the value of t4. process 4-1 is an isochoric process. therefore
%p1/t1=p4/t4,t4=p4t1/p1

t4=(p4*t1)/p1

% Calculating the thermal efficiency of the engine.
% the equation we are going to us is Thermal efficiency=1-((u4-u1)/(u3-u2))

% Here we are considereing the gas as an idela gas. Hence we can obtain the specific internal energy 
% from thermodynamic tables corresponding to different Temperatures. Therefore
% internal enegry of air corresponding to t1 500k=u1 =359.9 KJ/kg.Similarly 
% u2=1064 Kj/Kg, u3=1930 Kj/Kg ,u4= 632.15 Kj/Kg
therm_efficiency=(1-((632.15-359.9)/(1930-1064)))*100

Code Explanation

We take in required values. These include the bore,stroke,connecting rod length,t1,p1 etc.

Using these values the swept volume and clearance volume is calculated.

V1 and v2 are also found using this values.

p2 is found from the isentropic process equation p2v2^gamma=p1v1^gammat2 is found by the ideal gas equation

The values required to plot the points between states 1 and 2 are found out by calling the functiion engine_kinematics and storing in array v_compression.p_compression is found out by dividing v_compression elementwise with constant c. c=p1v1^gamma

p3 is also found out by applying ideal gas equation to the process 2-3. This is found out easily as it is an isochoric process.

v4 is equal to v1

The values required to join states 3 and 4 are found out by calling the function engine_kinematics and storing the resulting values in an array v_expansion.

p_expansion is found out by dividing elementwise with constant c.

After all the required values for plotting are found out the graph is plotted.

we are using the plot function for this. We are marking the state points using '*'.

Following this we are finding out the thermal efficiency. This is done by the equation thermal efficiency=1-((u4-u1)/(u3-u2)). Here u represents the specific internal energy at that temperature. This is found from the thermodynamic tables.

This is the code of the function used to find out the array of values used during plotting

%This function is used to help plot the PV diagram of an otto cycle,After already having obtained the state points.

 function[Volume]= engine_kinematics(bore,stroke,con_rod,CR,start_crank,end_crank) 
 clc

 %The equation we are going to be using is as follows:
 %(V/Vc)=1+0.5(CR-1)(R+1-cos(theta)-((R^2)-(Sin(theta)^2))^0.5)

 
 Vs=(pi/4)*(bore^2)*stroke;
 Vc=Vs/(CR-1);
 a=stroke/2;
 R=con_rod/a;
 theta=linspace(start_crank,end_crank,180);
 term1=0.5*(CR-1);
 term2=(R+1-cosd(theta));
 term3=((R^2)-sind(theta).^2).^0.5; 
 Volume=(1+term1*(term2-term3))*Vc;
 
 

  end
 

Code Explanation

The values required for computing the swept volume and clerance volume are received from the main program.t

he crank pin radius and the value of R is also found.

The value of ending and starting angle are also received from the main program.

The equation required for finding out the volume (V/Vc)=1+0.5(CR-1)(R+1-cos(theta)-((R^2)-(Sin(theta)^2))^0.5)

Using this the value of volume is saved as an array.

Results

Errors Faced.

Errors faced during programmingHad some issues with forgeting element wise operator in the term3 in the function. Worked perfectly after checking.

Not sure whether the method i have chosen for finding out efficiency is coorect. Assumed that i was supposed to find ouut efficiency from chart hence this method. 68 perc approx seems too high an efficiency. Please do comment.