Determination of PV diagram and thermal efficiency of an engine using MATLAB

The Otto cycle is a description of what happens to a mass of gas as it is subjected to changes of pressure, temperature, volume, addition of heat, and removal of heat. The mass of gas that is subjected to those changes is called the system. The system, in this case, is defined to be the fluid (gas) within the cylinder. By describing the changes that take place within the system,In the case of the Otto cycle, the effect will be to produce enough net work from the system so as to propel an automobile and its occupants in the environment.

The below fig(1) shows that how the piston movements and time taken to complete the process by the piston with in the cylinder

fig(1)

The pv diagram show that how the pressure and velocity will change in the otto cycle

Process:

• process A-B=reversible adiabatic compression
• process B-C=constant volume heat addition
• process C-D=reversible adiabatic expansion
• process D-A=constant volume heat rejection

Mathematical equations:

• thermal efficiency

`eta=1-(1/r^gamma-1)`

• swept volume

`v_s=(pi/4)*d*l`

• clerance volume

`v_c=v_s/(r_s-1)`

• initial volume

`v_1=v_s+v_c`

• isontropic compression

`p_1*v_1^gamma=p_2*v_2^gamma`

• isontropic expansion

`p_3*v_3^gamma=p_4*v_4^gamma`

• ideal gas equation

`(p_1v_1)/t_1=(p_2v_2)/t_2=(p_3v_3)/t_3=(p_4v_4)/t_4=cons`

• compression ratio

`c_r=v_1/v_2`

• mathamatical formula to find the volume that the piston can move

`v=(1+1/2(r_c-1)(R+1-costheta-(R^2-sin^2theta)^(1/2))*vc`

where R=ratio of connecting rod length to crank pin radius

Piston geometri:

Programming part:

To srtart the program initial inputs are

initialpressure=p1

initial temperature=t1

diameter of the piston=bore

length of the cylinder=stroke

length of the connecting rod=con_rod

compression ratio=cr=12

gamma=1.4

clear all
close all
clc

%inputs to run the program 
p1=101325;
t1=500;
bore=.1;
stroke=.1;
con_rod=.15;
cr=12;
gamma=1.4;
t3=3600

%calulate the swept volume & clerance volume
vs=(pi/4)*bore^2*stroke;
vc=vs/(cr-1); 
v1=vs+vc
v2=vc
v3=v2

% to identify p2
p2=p1*cr^gamma

%to calculate (t2)
t2=(p2*v2*t1)/(p1*v1)

constant=p1*v1^gamma;
v_compression=piston(bore,stroke,cr,con_rod,180,0)
p_compression=constant./v_compression.^gamma;


% TO CALUCLATE THE p3
p3=(p2*t3)/t2

constant=p3*v3^gamma;
v_expansion=piston(bore,stroke,cr,con_rod,0,180);
p_expansion=constant./v_expansion.^gamma;
v4=v1

%to calc p4
p4=p3*(v3/v4)^gamma

figure(1)
hold on
plot(v1,p1,'*','color','r')
plot(v_compression,p_compression)
plot(v2,p2,'*','color','b')
plot(v3,p3,'*','color','g')
line([v2 v3],[p2 p3])
plot(v_expansion,p_expansion)
plot(v4,p4,'*')
line([v1 v4],[p1 p4])
xlabel('volume ')
ylabel('pressure ')


Efficincy=1-(1/(v1/v2)^(gamma-1))

function syntax:

function [V]=piston(bore,stroke,cr,con_rod,start_crank,end_crank)
v_s=pi/4*bore^2*stroke;
vc=v_s/(cr-1); 


a=stroke/2;
r=con_rod/a;
theta=linspace(start_crank,end_crank,100);

term1=(cr-1)*.5;
term2=r+1-cosd(theta);
term3=((r^2-sind(theta).^2).^0.5);
V=(1+term1*(term2-term3)).*vc;
 

end

in the first program on line 28 and 36 "piston(bore,stroke,cr,con_rod,start_crank,end_crank)" is used to call the program using name of the funcation piston from 2nd program to 1st program

results:

The pv diagram is as shown below

the thermal efficiency:

p1 =101325 p2 =3.2853e+06 p3 =8754480 p4 =2.7001e+05

t1 =500 t2 =1.3510e+03 t3 =3600

v1 =8.5680e-04 v2 =7.1400e-05 v3 =7.1400e-05 v4 =8.5680e-04

Efficincy =0.6299