Creating a PV diagram for an Otto cycle

Here the objective is to plot PV diagram for an otto cycle given some assumptions as inputs in calculations. Some of the initial and final state variables are assumed for the ease of calculation which are given in the code. From the bore, stroke, and compression ratio, the engine cyclinder's swept volume and clearance volume is calculated. Using the ideal gas laws and adiabatic process laws, the state variables for all states 1,2,3,4 are calculated as shown in the code. This gives the 4 points in the PV diagram.

From the piston kinematic geometrical equation (`V/V_c=1+1/2(cr-1)[R+1-costheta-(R^2-sin^2theta)^(1/2)]`), the adiabatic volume variation with respect to the crank angle can be derived as shown in the piston kinematics function code. Using this volume and adiabatic process law, the adiabatic pressure variation can be calculated as shown in the main code.

The graph with all the state variables, the adiabatic processes and isochoric processes are plotted in a single plot. The thermal efficiency of this engine is calculated using the standard formula for efficiency with compression ratio and gamma(adiabatic gas constant).

Below, the main code is given.

%program to represent PV diagram of an otto cycle 
clear all
close all
clc

%inputs
bore = 0.1;        %bore diameter of the engine in m
stroke = 0.1;      %stroke length of the engine in m
con_rod = 0.15;    %connecting rod length of the engine in m
cr = 12;           %compression ratio of the engine

%assumptions
p1 = 101325;       %pressure at state 1 in Pa
t1 = 500;          %temperature at state 1(before compression) in degree K
gamma = 1.4;       %adiabatic gas constant for dry air at normal temperature
t3 = 2300;         %temperature at state 3(before expansion) in degree K

%calculations of state variables
V_s = (pi/4)*(bore^2)*stroke;        %swept volume of the engine in m^3
V_c = V_s/(cr-1);                    %clearance volume of the engine in m^3

v1 = V_s+V_c;      %total volume of the engine aka volume at state 1(before compression) in m^3
v2 = V_c;          %volume at state 2(after compression) also equal to top dead centre(clearance volume) in m^3

%from p1*(v1^gamma)=p2*(v2^gamma) and cr=v1/v2
p2 = p1*(cr^gamma);   %pressure at state 2(after compression) in Pa

%from ideal gas law: p*v/t= constant
t2 = p2*v2*t1/(p1*v1);  %temperature at state 2(after compression) in deg K

v3 = v2;           %from isochoric combustion volume at state 3 in m^3

%from ideal gas law: p*v/t= constant
p3 = p2*(t3/t2);   %pressure at state 2(after compression) in Pa

v4 = v1;           %from isochoric exhaustion volume at state 4 in m^3

%from p3*(v3^gamma)=p4*(v4^gamma)
p4 = p3*((v3/v4)^gamma);   %pressure at state 2(after compression) in Pa

%calculating volume and pressure change during adiabatic compression and expansion
V_compression = piston_kinematics(bore,stroke,con_rod,cr,180,0);   %in m^3
V_expansion = piston_kinematics(bore,stroke,con_rod,cr,0,180);     %in m^3

%from p1*v1^gamma=constant for adiabatic process
const_c1 = p1*(v1^gamma);
P_compression = const_c1./(V_compression.^gamma);                  %in Pa

const_c2 = p3*(v3^gamma);
P_expansion = const_c2./(V_expansion.^gamma);                      %in Pa

%plotting PV diagram and outputting thermal efficiency
figure(1)
hold on
plot(v1,p1,'*','color','r');
plot(v2,p2,'*','color','r');
plot(v3,p3,'*','color','r');
plot(v4,p4,'*','color','r');
xlabel('volume, V in m^3');
ylabel('pressure, P in Pa');
plot(V_compression,P_compression,'color','b');   %adiabatic compression curve
plot(V_expansion,P_expansion,'color','b');       %adiabatic expansion curve
plot([v2 v3],[p2 p3],'color','b');               %isochoric combustion curve
plot([v4 v1],[p4 p1],'color','b');               %isochoric exhaust curve

%from definition of thermal efficiency=1-((v2/v1)^(gamma-1)) and cr = (v1/v2) 
efficiency = 1-(1/(cr^(gamma-1)))

Below the function code for piston kinematics is given.

function [V] = piston_kinematics(bore,stroke,con_rod,cr,start_crank,end_crank)

 %defining terms of piston kinematic equation derived through piston geometry
 a = stroke/2;      %a factor in piston kinematic equation
 R = con_rod/a;     %a factor in piston kinematic equation
 V_s = (pi/4)*(bore^2)*stroke;        %swept volume of the engine in m^3
 V_c = V_s/(cr-1);                    %clearance volume of the engine in m^3

 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;

 %output volume in m^3
 V = V_c.*(1+term1*(term2-term3));

end

The output with the plot and efficiency value of 62.99% is shown below.

The errors faced during coding is as given below.

1. While calculating the adiabatic pressure variation, I forgot to put the element-wise operator for the power of gamma.

2. The factor R was not defined in the piston kinematics function code, so was not executing it. then I defined both 'R' and 'a' inside the function code.

3. Similarly I had a problem with the variable 'V_c' which then I defined in both main code and the function code since both codes required 'V_c' during their execution.

4. I had one more error with 'V_c' with respect to the multiplication operator. I had forgotten to put the operator itself.