Plotting P-V diagram of Otto cycle using Matlab

Project Description: An otto cycle is a thermodynamic cycle that signifies the functional working of a typical spark ignition engine. A spark ignition engine intakes air fuel mixture, compresses it, ignites the compress charge producing heat from the chemicals burnt and then uses this energy to push down the piston, hence creating mechanical energy. The piston then due to inertia again moves upwards and pushes out the burnt charge and is ready for repeating the cycle. To plot a P-V diagram of this otto cycle and calculate its thermal efficiency by using Matlab, is the soul aim of this project.

Processes: An otto cycle consists of four processes viz.,

1. isentropic compression(1-2): The compression process in S.I. engines takes place isentropically. i.e., with constant entropy. Thishappens when the piston moves from BDC to TDC, being full with the intakecharge and compresses it till TDC.
2. constant volume head addition(2-3): The heat addition processs starts with the spark plug. As the spark is generated the mixed air fuel charge gets ignited and generates heat energy from chemical. This head addition happens at constant volume.
3. isentropic expansion(3-4): After the burn, the heat generated forces the piston to move downwards from TDC to BDC. This is called as power stroke because here the heat energy is being actually converted into mechanical power. As the piston moves down the volume goes on increasing and expansion process takes place. This expansion occurs at constant entropy.
4. constant volume heat rejection(4-1): after the expansion process the pressure in the system is reduced and hence unused heat is rejected from the system at constant volume.

PLOTING THE CURVE:
Firstly, to plot P-V diagram we will need the (volume,pressure) coordinates at the four cycle points. For determining the process curves (2-3) and (4-1) of constant volume, we can obtain straight lines easily. But, for isentropic curves (1-2) and (3-4) we need to calculate the adiabatic constant at finite values of the pressure and volume. Let us divide the adiabatic curves in hundred equal parts. This can be easily done by dividing the crank angle in 100 for a range of 180-360 for compression curve and 0-180 for expansion curve. But to calculate these 100 values of the volume and pressure we will use function command in Matlab.

Code to plot the P-V diagram:

clear all
close all
clc

% code to plot P-V diagram

%inputs
bore=0.1;
stroke=0.1;
con_rod=0.15;
c_r=12;
gamma=1.4;
%inlet pressure and temperature of the cycle 
p1=101325;
t1=500;
%peak temperature of the cycle
t3=2300;
 

%swept volume and clearance volume
v_swept=pi/4*bore^2*stroke;
v_clearance=v_swept/(c_r-1);
v1=v_swept+v_clearance;
v2=v_clearance;
v3=v2;
v4=v1;
%calculating pressure at 2
% p1v1^gamma=p2v2^gamma
p2=p1*c_r^gamma;

%calculating intermediate volumes and pressures for compression curve
constant_c = p1*v1^gamma;
v_compression = engine_kinematics(0.1,0.1,0.15,12,180,360);
p_compression=constant_c./v_compression.^gamma;

%calculating temperature at t2
%p1v2/t1=p2v2/t2
t2=p2*v2*t1/(p1*v1);

%calculating pressure at 3
p3=p2*t3/t2;

%calculating intermediate volumes and pressures for expansion curve
constant_c = p3*v2^gamma;
v_expansion=engine_kinematics(0.1,0.1,0.15,12,0,180);
p_expansion=constant_c./v_expansion.^gamma;

%calculating pressure at 4
%p4*v46gamma=p3*v36gamma
p4=p3*(v3/v4)^gamma;


figure(1)
hold on
 
plot(v1,p1,'*','color','r');
plot(v_compression,p_compression,'color','r');
plot(v2,p2,'*','color','r');
plot([v2 v3],[p2 p3],'color','r')
plot(v3,p3,'*','color','r');
plot(v_expansion,p_expansion,'color','r');
plot(v4,p4,'*','color','r');
plot([v4 v1],[p4 p1],'color','r');
 
%efficiency of the cycle
n=1-1/c_r^(gamma-1);

For this code we need the minimum following input values :
bore=inner diameter of cylinder;
stroke=length of the piston travel;
con rod length;
crank pin radius=stroke/2;
R=connecting rod length/crank pin radius;
CR=compression ratio;
gamma=1.4;
inlet temperature and peak temperature=t1 and t3 resp.;
inlet pressure=p1;

Function to find varying volume at varying crank angles:

%function to calculate volume at different crank angles
function [V]= engine_kinematics(bore,stroke,con_rod,c_r,start_angle,end_angle)

a=stroke/2;
R=con_rod/a;
v_clearance=7.1400e-05;
theta=linspace(start_angle,end_angle,100);

term1=0.5*(c_r-1);
term2=R+1-cosd(theta);
term3=(R^2-sind(theta).^2).^0.5;

V=(1+((term1)*(term2-term3))).*v_clearance

end

From this function we can obtain different volume values from v1 to v2 during compression that are decreasing and increasing volume values from v3 to v4.
The formula used here is V=(1+0.5(CR-1)[R+1-cos(theta)-(R^2-sin(theta)^2)^0.5])/Vc
where,
V= cylinder volume
Vc=clearance volume
CR=comprssion ratio
theta=crank angle
R=connecting rod length/crank pin radius
crank pin radius = stroke/2

Errors:
The most common error faced during programming was not including ''.'' when multipication of matrix was done.
Another error was wrong syntax during plotting the lines of constant volume processes.

This above error is due to incorrect syntax given to plot command while plotting constant volume lines.

Result:
The P-V plot for otto cycle was successfully plotted in Matlab using function command spcifically to get the piston kinematics and then plot adiabatic curves.
This code is specifically best at calculating varying volumes at varying crank angles.
The efficiency of this cycle is 62.99=63%.

The above mentioned plot is tha P-V diagram for otto cycle.