Thermal Efficiency and P-V Plot for an Otto cycle
A programme to solve an otto cycle and make plots of its state points,its process and also calculate the thermal efficiency of that cycle.
This programme is used to plot the variations of pressure by changing volume. The inputs taken are:
a.) Constants used were gamma=1.4 ,t3=1200, pi=3.141 and geometric parameters of engine were taken as
bore=0.1
stroke=0.1
c_rod=0.15 (connecting rod)
cr=16 (compression ratio)
b.) we calculated the the swept volume and clearance volume by using these formulas:
v_swept=(pi/4)*bore*bore*stroke
v_clearance= v_swept/(cr-1)
c.) As the otto cycle is having 4 process (compression,ignition,expansion,intake), these will lead to 4 quantised states(state 1,2,3,4). Starting with the state 1, we calculated the pressure and volume at all states using the below datas and formula.
state variable 1
p1=101325 (pressure at state 1)
t1=500 (temperature at state 1)
v1=v_clearance + v_swept (volume at state 1)
constant=p1*v1^gamma (determining the constant for the process 1 to 2, using the formula p*v^gamma=constant)
state variable 2
v2= v_clearance (volume at state 2)
p2=p1*(cr^gamma) (pressure at state 2)
t2=(p2*v2*t1)/(p1*v1) (temperature at state 2)
state variable 3
v3=v2 (volume at state 3 is same as that of state 2)
p3=(p2*t3)/t2 (pressure at state 3)
state variable 4
v4=v1 (volume at state 4 is same as that of state 1)
p4=p3*(v3/v4)^gamma (pressure at state 4)
t4=[(p4*v4)*t3]/(p3*v3) (temperature at state 4)
constant=p3*(v3^gamma) (determining the constant for the process 3 to 4, using the formula p*v^gamma=constant)
d.) The above code will give us the value of pressure and volume at all the 4 state points. Joining these states simply wil give us inappropriate plots for expansion and compression processes.Joining these points will give us simply the straight lines, which is not the actual and correct plot. So for solving these problems we took the formula which calculates the variation of volume with changing crank angle(start angle 0 degree at Top Dead Centre and 180 degree at Bottom Dead Centre).
The relation used was:
V/V_clearance = 1+(1/2)*(cr-1)[R+1-cos(theta)-√{R^2−sin(theta)^2}].
For creating the exact plot we define a function names engine_kinematics , which calculates the exact pressure and volume at all the points in the compression and expansion process when called from an another programe. File name and the function name should be same.
we use the function definition to define the function definition as shown below:
function [v]=engine_kinematics(bore,stroke,c_rod,cr,start_crank,end_crank);
this function ends with an "end" keyword.
we again calculated the the swept volume and clearance volume by using these formulas:
v_swept=(pi/4)*bore*bore*stroke
v_clearance= v_swept/(cr-1)
we define an array variable as theta using linspace(start_crank,end_crank,100), where theta is an array of angle starting from start_crank to end_crank angle with 100 divisions(for compression stroke start_crank is 180 and end_crank is 360, whereas for expansion stroke the start_crank is 0 and end_crank is 180).
Then the formula of variation of swept volume with changing crank angle is broken down into three terms for simplification. As given below:
1.) term1=0.5*(cr-1)
2.) term2=R + 1 - cosd(theta)
3.) term3=(R^2 - sind(theta).^2).^0.5
And now the full expression of total volume is assembeled as shown below:
v=(1+ term1*(term2-term3)).*v_clearance
e.) Now the above function is called in the main programme by passing the values of its parameters such as bore,stroke,c_rod,cr,start_crank,end_crank for process 1 to 2 and 3 to 4 to calculate (p_compression ,v_compression) and (p_expansion, v_expansion).
For calculating p_compression and v_compression we call the function as follows:
v_compression=engine_kinematics(bore,stroke,c_rod,cr,180,360);
p_compression=constant./v_compression.^gamma;
And for calculating p_expansion and v_expansion we call the function as follows:
v_expansion=engine_kinematics(bore,stroke,c_rod,cr,0,180);
p_expansion=constant./v_expansion.^gamma;
f.) Now for finding the thermal efficiency of the otto cycle we use the following relations:
etta_main=1-(t1/t2), where etta_main is a variable for the thermal efficiency of the otto cycle, t1&t2 are the therodynamic temperature at state 1 and 2.
But for plotting the variation of thermal effeciency with changing compression ratio,we create an array of compression ratio starting from 0 to 16 with 100 values in between.The code for this calculation is shown below:
cr=linspace(0,16,100)
etta = [1-1./(cr).^(gamma-1)]*100 %this is the relation between thermal efficiency and %compression ratio.
g.) To plot all the 4 state points we use plot command and we also use the plot command to join all the (p_compression ,v_compression) and (p_expansion, v_expansion) points, as shown below:
For plotting the states and actual processes of the otto cycle curve.
figure(1)
hold on
plot(v1,p1,'*','color','r')
plot(v_compression,p_compression)
plot(v2,p2,'*','color','b')
plot(v3,p3,'*','color','g')
plot(v_expansion,p_expansion)
plot(v4,p4,'*','color','y')
plot([v2 v3],[p2 p3],'color','r')
plot([v1 v4],[p1 p4],'color','b')
xlabel('VOLUME [V] in (m^3)') %x-axis is named as VOLUME [V] in (m^3)
ylabel('PRESSURE [P] in (N/m^2)') %y-axis is named as PRESSURE [P] in (N/m^2)
grid on
And for plotting the approx paths of otto cycle.
figure(2)
plot([v1 v2 v3 v4 v1],[p1 p2 p3 p4 p1],'color','r');
grid on
xlabel('VOLUME [V] in (m^3)') %x-axis is named as VOLUME [V] in (m^3)
ylabel('PRESSURE [P] in (N/m^2)') %y-axis is named as PRESSURE [P] in (N/m^2)
x-axis is named as VOLUME and y-axis is named as PRESSURE.
And for plotting the variation of thermal efficiency with compression ratio we use the code below.
figure(3)
hold on
plot(cr,etta,'color','b');
axis([0 20 0 100])
xlabel('COMPRESSION RATIO(cr)')
ylabel('THERMAL EFFICIENCY(etta)')
grid on
CODE FOR PLOTTING THE ACTUAL PV DIAGRAM, APPROX PV DIAGRAM & THERMAL EFFICIENCY VS COMPRESSION RATIO PLOT:
% programme to calculate and plot the otto cycle state points.
close all
clear all
clc
%constants used
gamma=1.4
t3=1200
pi=3.141
% geometric parameters of engine
bore=0.1
stroke=0.1
c_rod=0.15
cr=16
% calculation of of swept volume and clearance volume
v_swept=(pi/4)*bore*bore*stroke
v_clearance= v_swept/(cr-1)
% calculating state variable 1
p1=101325 %setting the initial state variable
t1=500
v1=v_clearance + v_swept
constant=p1*v1^gamma %using p*v^gamma=constant
% calculating state variable 2
v2= v_clearance
p2=p1*(cr^gamma)
t2=(p2*v2*t1)/(p1*v1)
v_compression=engine_kinematics(bore,stroke,c_rod,cr,180,360); %calling the
%engine_kinematics
%function by passing
%values of the
%parameters and
%calculating an array
%of pressure and
%volume
%for process 1 to 2.
p_compression=constant./v_compression.^gamma; %using
%p*v^gamma=constant. For
%compression start angle
%is 180 and end angle is
%360 degree.
% calculating state variable 3
v3=v2
p3=(p2*t3)/t2
% calculating state variable 4
v4=v1
p4=p3*(v3/v4)^gamma
t4=[(p4*v4)*t3]/(p3*v3)
constant=p3*(v3^gamma)
v_expansion=engine_kinematics(bore,stroke,c_rod,cr,0,180); %calling the
%engine_kinematics function
%by passing values of the
%parameters and calculating
%an array of pressure and
%volume for process 3 to 4.
p_expansion=constant./v_expansion.^gamma; %using p*v^gamma=constant.
%For expansion start angle
%is 0 and end angle is 180
%degree.
%calculating thermal efficiency of the otto cycle.
cr=linspace(0,16,100);
etta = [1-1./(cr).^(gamma-1)]*100;
etta_main=1-(t1/t2);
%plotting the states and actual processes of the otto cycle curve.
figure(1)
hold on
plot(v1,p1,'*','color','r')
plot(v_compression,p_compression) %plotting the values of the array of v_compression
%and p_compression.
plot(v2,p2,'*','color','b')
plot(v3,p3,'*','color','g')
plot(v_expansion,p_expansion) %plotting the values of the array of v_expansion
%and p_expansion.
plot(v4,p4,'*','color','y')
plot([v2 v3],[p2 p3],'color','r')
plot([v1 v4],[p1 p4],'color','b')
xlabel('VOLUME [V] in (m^3)') %x-axis is named as VOLUME [V] in (m^3)
ylabel('PRESSURE [P] in (N/m^2)') %y-axis is named as PRESSURE [P] in (N/m^2)
grid on
%plotting the approx paths of otto cycle.
figure(2)
hold on
plot([v1 v2 v3 v4 v1],[p1 p2 p3 p4 p1],'color','r');
grid on
xlabel('VOLUME [V] in (m^3)') %x-axis is named as VOLUME [V] in (m^3)
ylabel('PRESSURE [P] in (N/m^2)') %y-axis is named as PRESSURE [P] in (N/m^2)
%polt of thermal efficiency vs changing compression ratio(cr).
figure(3)
hold on
plot(cr,etta,'color','b');
axis([0 20 0 100]) %setting the co-ordinate axis.
xlabel('COMPRESSION RATIO(cr)') %x-axis is named as COMPRESSION RATIO(cr)
ylabel('THERMAL EFFICIENCY(etta)') %y-axis is named as THERMAL EFFICIENCY(etta)
grid onCODE FOR DEFINING THE engine_kinematics FUNCTION:
% programme to define a function named engine_kinematics and calculate the variation of %swept volume with changing crank angle in an otto cycle.
function [v]=engine_kinematics(bore,stroke,c_rod,cr,start_crank,end_crank);
%constants used
gamma=1.4
t3=1200
pi=3.141
% geometric parameters of engine
a=stroke/2;
R=c_rod/a;
% calculation of of swept volume and clearance volume
v_swept=(pi/4)*bore*bore*stroke
v_clearance= v_swept/(cr-1)
theta=linspace(start_crank,end_crank,100) %theta is an array of angle starting from
%start_crank to end_crank angle.
term1=0.5*(cr-1) %the formula of variation of swept volume with
%changing crank angle is broken down into three terms
%for simplification.
term2=R + 1 - cosd(theta)
term3=(R^2 - sind(theta).^2).^0.5
v=(1+ term1*(term2-term3)).*v_clearance %this expression calculates the variation
%of swept volume with change in the crank
%angle.
endPLOTS OF THE RESULTS:
1.) Plot of the states and actual processes of the otto cycle curve
2.)Plot of the approx paths of otto cycle.
3.)Plot of the variation of thermal efficiency with compression ratio.