PV diagram visualizer in MATLAB using engine kinematics

Objective: To plot the PV diagram of the otto cycle and to determine the thermal efficiency of the Otto cycle.

Background:

Processes in Otto Cycle:

Process 1-2: Isentropic compression

In this process, the piston moves from the bottom dead center (BDC) to top dead center (TDC) position. Air undergoes reversible adiabatic (isentropic) compression. We know that compression is a process in which volume decreases and pressure increases. Hence, in this process, the volume of air decreases from V₁ to V₂ and pressure increases from p₁ to p₂. Temperature increases from T₁ to T₂. As this an isentropic process, entropy remains constant

Process 2-3: Constant Volume Heat Addition:

Process 2-3 is an isochoric (constant volume) heat addition process. Here, piston remains at the top dead center for a moment. Heat is added at constant volume (V₂ = V₃) from an external heat source. Temperature increases from T₂ to T₃, pressure increases from p₂ to p₃ and entropy increases.

Process 3-4: Isentropic expansion

In this process, air undergoes isentropic (reversible adiabatic) expansion. The piston is pushed from the top dead center (TDC) to bottom dead center (BDC) position. Here, pressure decreases from p₃ to p₄, volume rises from v₃ to v₄, the temperature falls from T₃ to T₄ and entropy remains constant

Process 4-1: Constant Volume Heat Rejection

The piston rests at BDC for a moment and heat is rejected at constant volume (V₄=V₁). In this process, the pressure falls from p₄ to p₁, temperature decreases from T₄ to T₁ and entropy falls.

CODE:

clear all
close all
clc 

% inputs
gamma = 1.4
T3 = 2300 % kelvin
Cv = 45500 % calorific value KJ/Kg
mf = 1 % mass of fuel in kg

% state variables at state point 1
P1 = 101325 % Pa
T1 = 500    % kelvin

% engine geometric parameters
bore = 0.1; % in m
stroke = 0.1; % stroke length in m
con_rod = 0.15;% connecting rod in m
cr = 12; % compression ratio

%computing the values of swept volume and clearance volume
V_swept = (pi/4)*bore^2*stroke;
V_clearance = V_swept/(cr-1);

V1 = V_swept + V_clearance
V2 = V_clearance

% state variables at state point 2
% wkt P1* V1^gamma = P2* V2^gamma, P2 = (P1*V1^gamma)/V2^gamma =
% P2 = P1*(cr)^gamma
P2 = P1*(cr)^gamma;
%((P1*V1)/T1)= ((P2*V2)/T2)
T2 = (P2/P1)*T1*(1/cr);

constant_c = P1*V1^gamma
V_compression = engine_kinematics(bore,stroke,con_rod,cr,180,360)
P_compression = constant_c./ V_compression.^gamma;

% state variables at state point 3
V3 = V2;

%((P2*V2)/T2)= ((P3*V3)/T3)
P3 = (T3/T2)*P2*(V2/V3)

constant_c = P3*V3^gamma
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
% wkt P3* V3^gamma = P4* V4^gamma, P4 = (P3*V3^gamma)/V4^gamma =
% P4 = P3*(1/cr)^gamma
P4 = P3*(1/cr)^gamma;


figure(1)
hold on
plot(V1,P1, 'color', 'r','linewidth', 4)
plot(V_compression,P_compression, 'color', 'r','linewidth', 4)
plot([V2 V3], [P2 P3], 'color', 'b', 'linewidth', 4)
plot(V_expansion,P_expansion, 'color', 'r', 'linewidth', 4)
plot([V4 V1], [P4 P1], 'color', 'b', 'linewidth', 4)
xlabel('Volume in m^3');
ylabel('Pressure in Pa');
grid on
fill([V1 V_compression V2 V3 V_expansion V4 V1], [P1 P_compression P2 P3 P_expansion P4 P1],'g')% The colored area represents the work done

Otto_cycle_efficiency = (1 - (1/(cr)^(gamma-1)))*100
WD = Otto_cycle_efficiency*mf*Cv %Work done in KJ
WD_in_MJ = WD/1000

EXPLANATION:

The inputs are mentioned as gamma = 1.4(ratio of CP/CV) for air, Temperature at state point 3 (T3=2300K) and calorific value (Cv) = 45.5 MJ/Kg and mass of fuel (mf) to determine the Work done by the cycle.

The additional values of bore diameter(bore)= 0.1m, stroke length(stroke) = 0.1m, connecting rod length(con_rod) of 0.15m, compression ratio(cr) = 12 and the initial state point conditions at P1= 101.32 Kpa atm and T1 =500k are established.

V1 and V2 are determined by using the formula:

V_swept = (pi/4)*bore^2*stroke
V_clearance = V_swept/(cr-1)

V1 = V_swept + V_clearance
V2 = V_clearance

Varaiables at state point 2, P2 and T2 are determined by the equations:

P1* V1^gamma = P2* V2^gamma,
P2 = P1*(cr)^gamma
((P1*V1)/T1)= ((P2*V2)/T2)
T2 = (P2/P1)*T1*(1/cr)

Followed by the compression volume determined by calling a function engine_kinematics as mentioned in the code below:

This code defines the geometrical aspects to determine the compression volume using the formula

V = (1+(1/2)(cr-1)*((R+1-cos(theta)-(R^2-sin(theta)^2)^1/2) * V_clearance

function [V] = engine_kinematics(bore,stroke,con_rod,cr,start_angle,end_angle)

a = stroke/2; %crank radius
R = con_rod/a; %ratio of connecting rod and crank radius

V_s = (pi/4)*bore^2*stroke;
V_c = V_s/(cr-1);

theta = linspace(start_angle,end_angle,180)

t1 = 0.5*(cr-1)
t2 = R+1-cosd(theta);
t3 = (R^2 - sind(theta).^2).^0.5

V = (1+t1*(t2-t3)).* V_c;

end

The formula is split into 3 separate terms for ease of calculation and eliminating errors. The function output is V and as for the input it is bore, stroke, connecting rod, compression ratio, start angle, and end angle of the crank defined for compression between the angle 180 to 360 deg.

WKT P1*V1^gamma= constant
V_compression is determined from the function defined. followed by
P_compression = constant/ V_compression^gamma

state variables at state point 3 are determined by
V3 = V2

and

((P2*V2)/T2)= ((P3*V3)/T3)
P3 = (T3/T2)*P2*(V2/V3)

As for the expansion stroke the function engine kinematics is called in again to compute the variation of the curve

P3*V3^gamma= constant
P_expansion = constant/ V_expansion^gamma

state variables at state point 4. Determined from the following equations.
V4 = V1
P4 = P3*(1/cr)^gamma

RESULTS :

Results are plotted using the plot command at each state point and for the corresponding expansion and compression strokes.

The plot command is defined with appropriate linewidth and color to accentuate the variation.

The fill command is provided to highlight the area enclosed in the cycle to notify the work done by the otto cycle.

Further the equations to determine the thermal efficiency

Thermal efficiency = 1- (1/cr^(gamma-1)) and the work done is also computed for the Cv and mass of fuel.

The results can be viewed and validated below for a better understanding of the output received.