working of Air standard Cycle (otto cycle)
In this project I will be calculating thermalefficiency of the engine and I will draw a PV diagram.
First of all, lets gather the inputs.
Given Inputs: gamma for air=G=1.4
t3=2300 k
%state variables at point 1; p1= 101325 and t1=500k
as we know p1 and t1, we can calculate v1.
But first have a look at engine dimension, bore=0.1; stroke=0.1; connecting rod = cn_rod= 0.15; compression ratio= cr= 12;
Calculating swept and clearance volume, swept_volume = v_s=(pi/4)*bore^2*stroke; clearance_volume = v_c= v_s/(cr-1); As cylinder moves from TDC to BDC, the air get sucked into the cylinder, so the volume at point 1 is: v1=v_s + v_c;
Now, the piston moves from BDC to TDC, so volume at point 2 will be equal to clearance volume: v2=v_c;
State variables at point 2: p2v2^G=p1v1^G p2=p1*(v1/v2)^G=p1*cr^G p2=p1*cr^G;
By using Ideal Gas law: p1v1/t1=p2v2/t2 t2=(p2*v2*t1)/(p1*v1) t2=(p2*v2*t1)/(p1*v1)
State variables at point 3: constant volume heat addition Therefore, v3=v2 p2v2/t2=p3v3/t3 p3=p2*(v2/v3)*(t3/t2) p3=p2*(v2/v3)*(t3/t2)
State variable at point 4: constant volume heat rejection Therefore, v4=v1
To find pressure at point 4: p4v4^G=p3v3^G p4=p3*(v3/v4)^G=p3*(v3/v4)^G
Now we can plot the pv diagram,
hold on plot([v1 v2],[p1 p2]) plot([v2 v3],[p2 p3]) plot([v3 v4],[p3 p4]) plot([v4 v1],[p4 p1])
check the code:
close all
clear all
clc
%inputs
G=1.4;
t3=2300;
%state variables at point 1
p1= 101325;
t1= 500;
%Engine dimensions
bore=0.1;
stroke=0.1;
cn_rod= 0.15;
cr=12;
%Calculating swept and clearance volume
v_s=(pi/4)*bore^2*stroke;
v_c= v_s/(cr-1);
v1=v_s + v_c;
v2=v_c;
%state variables at point 2
%p2v2^g=p1v1^G | p2=p1*(v1/v2)^G=p1*cr^G
p2=p1*cr^G;
%Ideal Gas law: p1v1/t1=p2v2/t2 | t2=(p2*v2*t1)/(p1*v1)
t2=(p2*v2*t1)/(p1*v1);
% state variables at point 3
% constant volume heat addition
v3=v2;
%p2v2/t2=p3v3/t3 | p3=p2*(v2/v3)*(t3/t2);
p3=p2*(v2/v3)*(t3/t2);
%state variable at point 4
v4=v1;
p4=p3*(v3/v4)^G;
hold on
plot([v1 v2],[p1 p2])
plot([v2 v3],[p2 p3])
plot([v3 v4],[p3 p4])
plot([v4 v1],[p4 p1])
but we will get this pv diagram,
We can clearly see that, From point 1 to 2 the line we are getting is straight. But we know that the process is Isentropic expansion process. So to get ride of this issue we have solution.
Find v_compression and p_compression separately, We have values of bore, stroke,compression ratio, connecting rod lenght, so we can calculate crank pin radius and ratio of connecting rod length to crank pin radius i.e. R. So, crank pin radius=a=stroke/2 R=cn_rod/a
Now, we have this formula to calculate volume by using different values of theta. In this case we have to find v_compression, so values of theta are gonne be 0 to 180.
theta=linspace(0,180,180)
Equation:
v/vc=1+(1/2)(cr-1)[R+1-cosd(theta)-(R^2-sind(theta)^2)^0.5]
Here, v= v_vompression
As the equation is very complex, so for our convinience I will break this into three terms.
term1=0.5*(cr-1) term2=R+1-cosd(theta) term3=(R^2-sind(theta).^2).^0.5 v_compression/v_c= (1+ term1*(term2-term3)) v_compression= (1+ term1*(term2-term3)).*v_c
close all
clear all
clc
%Engine dimensions
bore=0.1;
stroke=0.1;
cn_rod= 0.15;
cr=12;
%crank pin radius a=stroke/2
a=stroke/2;
R= cn_rod/a;
v_s=(pi/4)*bore^2*stroke;
v_c= v_s/(cr-1);
% v/vc=1+(1/2)(cr-1)[R+1-cosd(theta)-(R^2-sind(theta)^2)^0.5]
theta=linspace(0,180,180);
term1=0.5*(cr-1);
term2=R+1-cosd(theta);
term3=(R^2-sind(theta).^2).^0.5;
v_compression= (1+ term1*(term2-term3)).*v_c;
plot(theta,v_compression)
We can see that, as this is compression stroke, volume increases as theta increases. Using v_compression, we can find p_compression, From the formula, p1v1^G=p2v2^G, here p1v1^G is constant, So, constant for compression=constant_c=p1v1^G
Therefore, p_compression=constant_c./v_compression.^G
Similarly, we can find v_expansiona and p_expansion using theta 180 to 360. Expansion stroke.
But, this is very bulky and difficult to understand if we put this commands in our code. So to simplify we can use function command.
Please check the function command,
function [v]=engine_kinematics(bore,stroke,cr,cn_rod,start_crank,end_crank)
%crank pin radius a=stroke/2
a=stroke/2;
R= cn_rod/a;
v_s=(pi/4)*bore^2*stroke;
v_c= v_s/(cr-1);
% v/vc=1+(1/2)(cr-1)[R+1-cosd(theta)-(R^2-sind(theta)^2)^0.5]
theta=linspace(start_crank,end_crank,180);
term1=0.5*(cr-1);
term2=R+1-cosd(theta);
term3=(R^2-sind(theta).^2).^0.5;
v= (1+ term1*(term2-term3)).*v_c;
end
From the above function command, we can easily calculate v_compression and v_expansion.
Volume is the function of bore,stroke,cr,con_rod and thetas.
So by simply calling the function using different start_crank and end_crank angle we can find volume of any point.
Use equation of v_compression and v_expansion in their respective parts in the final code.
Now the thermal efficiency is calcluated by using:
Thermal_efficiency=work done/ heat supplied
Work done= heat addition -heat dissipated = (t3-t2) - (t4-t1)
heat supplied= t3-t2
Therefore, Thermal_efficiency= 1 - (t4-t1)/(t3-t2)
Pleas echck the code.
close all
clear all
clc
%inputs
G=1.4;
t3=2300;
%state variables at point 1
p1= 101325;
t1= 500;
%Engine dimensions
bore=0.1;
stroke=0.1;
cn_rod= 0.15;
cr=12;
%Calculating swept and clearance volume
v_s=(pi/4)*bore^2*stroke;
v_c= v_s/(cr-1);
v1=v_s + v_c;
v2=v_c;
%state variables at point 2
%p2v2^g=p1v1^G | p2=p1*(v1/v2)^G=p1*cr^G
p2=p1*cr^G;
%Ideal Gas law: p1v1/t1=p2v2/t2 | t2=(p2*v2*t1)/(p1*v1)
t2=(p2*v2*t1)/(p1*v1);
v_compression=engine_kinematics(bore,stroke,cr,cn_rod,0,180)
constant_c=p1*v1^G
p_compression=constant_c./v_compression.^G
% state variables at point 3
% constant volume heat addition
v3=v2;
%p2v2/t2=p3v3/t3 | p3=p2*(v2/v3)*(t3/t2);
p3=p2*(v2/v3)*(t3/t2);
%state variable at point 4
v4=v1;
v_expansion=engine_kinematics(bore,stroke,cr,cn_rod,180,360)
p_expansion=p3.*(v3./v_expansion).^G
%p4v4^G=p3v3^G | p4=p3*(v3/v4)^G
p4=p3*(v3/v4)^G;
%Ideal Gas law: p3v3/t3=p24v4/t4 | t4=(p4*v4*t3)/(p3*v3)
t4=(p4*v4*t3)/(p3*v3);
hold on
plot(v1,p1)
plot(v_compression,p_compression)
plot([v2 v3],[p2 p3])
plot(v3,p3)
plot(v_expansion,p_expansion)
plot([v4 v1],[p4 p1])
xlabel('volume')
ylabel('pressure')
%Thermal_efficiency=work done/ heat supplied
Thermal_efficiency= 1 - (t4-t1)/(t3-t2)
Thank you |