Plotting a PV diagram & finding efficiency of an Otto cycle

Objective: To write a code that can solve an otto cycle

1. To plot PV diagram &
2. To find thermal efficiency of the engine.

Description:

The four-stroke air standard Otto cycle is made up of the following four internally reversible processes:

The piston moves within the cylinder between TDC and BDC. In order to complete a full combustion cycle, the piston executes four strokes and the crankshaft makes two complete turns. The swept/displaced volume is the volume in which the piston moves, the clearance volume is the volume left in the cylinder when the piston reaches TDC.

For one cylinder, the Swept/Displaced volume : `V_d=S⋅A_c`

The area of the cylinder is: `A_c=π/4⋅B^2`

Clearance volume : `Vc=V_s/(cr−1)`

i.e.. `V c = V _2` & `V s + V c = V _1`

Compression Ratio (cr) is the ratio between the maximum (total) volume of the cylinder (when the piston is at BDC) and the minimum (clearance) volume (when the piston is as TDC).

`cr=V_(max)/V_(min)=V_1/V_2`

Crank angle radius:`a=S/2`

connecting rod length to crank radius ratio: `R=r/a`

As 2-3 & 4-1 processes are constant, we get straight lines which are easy to plot where as the processes 1-2 & 3-4 are isentropic and the curve will not be linear. It is important to plot the PV diagram correctly as the area under a process curve on a p-V diagram is equal to the work performed by air during the process.

Isentropic process follows the equation: `V/V_c = 1 + 1/2*(cr-1)[R + 1 - costheta - (R^2 - sin^2theta)^(1/2)]`

In this project the otto cycle is solved by taking a general problem statement:

In an air standard Otto cycle, the compression ratio is 7 and the compression begins at 35^oC and 0.1MPa. The maximum temperature of the cycle is 1100^oC. Find (a) the temperature and the pressure at various points in the cycle and plot PV diagram, (b) the cycle efficiency.

MATLAB code with explanation:

%Otto cycle
clear all
close all
clc

%Given data in the problem statement
cr = 7;           %compression ratio
t1 = 35+273;      %K
p1 = 100;         %bars
t3 = 1100+273;    %K

%% Need to find t2,t4,p2,p3,p4,v1,v2,v3=v2,v4=v1

Rc = 0.287;       %Characteristic gas constant of air
g  = 1.4;         %Isentropic index gamma(g)
Cv = 0.718;       %Specific Volume
theta = 0:360;    %Crank angle

t2 = t1*cr^(g-1);

p2 = p1/(1/cr)^g; %p1v1^g = p2v2^g
p3 = p2*t3/t2;    %p2/t2 = p3/t3
p4 = p3/(cr)^g;   %p3v3^g = p4v4^g

v1 = Rc*t1/p1;    %p1v1=Rc.t1
v2 = v1/7;
v3 = v2;
v4 = v1;

t4 = p4*v4/Rc;    %p4v4=Rc.t4

%% To find P & V w.r.t crank angle
% for adiabetic process pv^g=constant(k)
k1 = p1*v1^g;
k2 = p3*v3^g;

for i=1:length(theta)
    if i <=180          %compression stroke (0 - 180 degrees)
v_comp(i) = isentropic_process(cr,theta(i),v2);
p_comp(i) = k1/v_comp(i)^g;
    elseif i>180        %expansion stroke (181 - 360 degrees)
v_expa(i) = isentropic_process(cr,theta(i),v2);
p_expa(i) = k2/v_expa(i)^g;
    end
end

%% To find thermal Efficiency
Q_in  = Cv*(t3-t2);
Q_out = Cv*(t4-t1);
Thermal_Efficiency = (Q_in-Q_out)/Q_in

%% To plot PV diagram 
figure(1)
plot ([v2,v3],[p2,p3],'r','linewidth',2)
hold on
plot ([v4,v1],[p4,p1],'r','linewidth',2)
plot (v_comp,p_comp,'r','linewidth',2)
plot (v_expa(181:360),p_expa(181:360),'r','linewidth',2)
xlabel('Volume(V)')
ylabel('Pressure(P)')
title(sprintf('PV diagram with Thermal efficiency %g',Thermal_Efficiency))

Function to be called

function [V] = isentropic_process(cr,theta,v2)
S = 0.1;
r = 0.15;% connecting rod length
a = S/2; % crank pin radius
R = r/a;

term1 = 0.5*(cr-1);
term2 = R + 1 - cosd(theta);
term3 = (R^2 - sind(theta)^2)^0.5;
V = (1 + term1*(term2-term3))*v2;
end

Output:

Thermal_Efficiency = 0.5408