Analysis of otto cycle using MATLAB
Aim
To analyse Otto cycle and plot its P-V diagram
Governing equations
- Ideal gas equation: `(P*V)/T` = constant. This equation is true for ideal gases and since Otto cycle is an air standard cycle, this equation holds true.
Here, P = pressure
T = temperature
V = volume
- For adiabatic process, `P*V^gamma` = constant
`gamma` = `C_P/C_v`
Value of `gamma` for air = 1.4
- Swept volume, `V_(swept)` = `pi/4*d^2*l`
Here, d = bore or diameter of piston cylinder
l = stroke or distance traversed by piston
- Clearance volume = `V_(cl earance)` = `V_(swept)/(r-1)`
- `V/V_c = 1 + 1/2(r-1)[R+1-costheta-(R^2-sin^2theta)^0.5]`
- Efficiency of Otto cycle: `eta` = work done/heat added
= `(Q_(i n) - Q_(out))/Q_(i n) = 1 - Q_(out)/Q_(i n)`
= `1-(m*C_v*(T_4-T_1))/(m*C_v*(T_3-T_2)`
= 1 – `(T_4-T_1)/(T_3-T_2)`
The above equation can be written in the terms of compression ratio, r as
`eta = 1-1/(r^(gamma-1))`
Objective
The objective of this program is to analyse the Otto cycle, plot its P-V curve and determine its thermal efficiency. An Otto cyclee is an air standard thermodynamic cycle that describes working of a typical spark ignition piston engine in which a mass of gas (usually air) is subjected to changes in pressure, temperature, volume, heat addition and rejection. There are usually 4 processes involved in a 2 stroke Otto cycle.
The 4 processes involved in an Otto cycle
1) Process 1-2: It is the adiabatic or isentropic (constant entropy) compression of charge in which the piston moves from the Bottom Dead centre (BDC) to Top Dead Centre (TDC). For this process, `P*V^gamma`= constant. Therefore, `P_1*V_1^gamma = P_2*V_2^gamma`
2) Process 2-3: It is isochoric (constant volume) heat addition while piston remains at TDC. It causes fuel ignition and rapid burning take place. Therefore, `V_2 = V_3`.
3) Process 3-4: It is the adiabatic (isentropic) expansion of charge in which piston moves from TDC to BDC. Again, `P*V^gamma`=constant. This implies `P_3*V_3^gamma = P_4*V_4^gamma`.
4) Process 4-1: It is the isochoric heat rejection while the piston remains at BDC. Thus, `V_4 = V_1`.
Given the state variables at point 1 and the temperature after the heat addition process, the state variables at all the states are determined using the suitable governing equations.
The nature of volume change during the adiabatic processes is determined using the governing equation and then the plot between pressure and volume is generated.
At last, thermal efficiency of the cycle is calculated using the formula.
Program
%Program for analysis of air compression cycle
%inputs
gamma = 1.4; %gamma = C_p/C_v
t3 = 2500;
%geometric parameters
d = 0.2; %bore
l = 0.4; %stroke
c_r = 0.25; %length of connecting rod
r = 6; %compression ratio
%state variables
p1 = 101325; %1atm = 101325 pascal
t1 = 500;
%calculating swept volume and clearance volume
v_swept = (pi/4)*d^2*l;
v_clearance = v_swept/(r-1);
v1 = v_swept + v_clearance;
v2 = v_clearance;
%state 2 variables
%p1*v1^gamma = p2*v2^gamma
p2 = p1*r^gamma; %v1/v2 = r
t2 = p2*v2*t1/(p1*v1); %ideal gas equation
constant_c = p1*v1^gamma;
v_compression = otto_cycle(d,l,c_r,r,180,0);
p_compression = constant_c./v_compression.^gamma;
%state 3 variables
v3 = v2;
p3 = p2*t3/t2; %ideal gas equation
constant_c = p3*v3^gamma;
v_expansion = otto_cycle(d,l,c_r,r,0,180);
p_expansion = constant_c./v_expansion.^gamma;
%state 4 variables
v4 = v1;
%p3*v3^gamma = p4*v4^gamma
p4 = p3/r^gamma; %r = v3/v4
t4 = p4*v4*t3/(p3*v3); %ideal gas equation
%plotting
hold on
plot(v1,p1,'*','color','r') %plotting the point
text(v1,p1,'1','FontSize',14,'color','r') %to label the points on the plot
plot(v_compression,p_compression,'linewidth',2,'color','b') %plotting the curve
plot(v2,p2,'*','color','r')
text(v2,p2,'2','FontSize',14,'color','r')
plot([v2,v3],[p2,p3],'linewidth',2,'color','b')
plot(v3,p3,'*','color','r')
text(v3,p3,'3','FontSize',14,'color','r')
plot(v_expansion,p_expansion,'linewidth',2,'color','b')
plot(v4,p4,'*','color','r')
text(v4,p4,'4','FontSize',14,'color','r')
plot([v4,v1],[p4,p1],'linewidth',2,'color','b')
xlabel('volume')
ylabel('pressure')
%efficiency of otto cycle
% eta = work done/heat added
%eta = (heat added - heat rejected)/heat added
eta = 1 - (t4 - t1)/(t3 - t2)
Initially, the required inputs are provided which include the geometric parameters (stroke, bore, connecting rod length), compression ratio and required state variables.
Then state variables for different state points are calculated, For plotting, hold on command is used to keep various plots and points corresponding to various states are plotted.
In order to determine the nature of curve for the adiabatic (isentropic) processes, a function is defined which is then called in the main program to get the accurate plot. Given below is the function defined to get actual change in volume during adiabatic compression and expansion processes.
function[v] = otto_cycle(d,l,c_r,r,start_crank,end_crank)
a = l/2; % crank pin radius
R = c_r/a;
v_s = (pi/4)*d^2*l;
v_c = v_s/(r-1);
theta = linspace(start_crank,end_crank,180);
% v/v_c = 1 + (1/2)(r-1)[R + 1 - cos(theta) - (R^2 - Sin^2(theta))^0.5]
term1 = 0.5*(r-1);
term2 = R + 1 - cosd(theta);
term3 = (R^2 - sind(theta).^2).^0.5;
v = (1 + term1*(term2 - term3)).*v_c;
endA function nammed otto_cycle is created with volume (v) as output a various inputs like bore, stroke, connecting rod length, compression ratio, start crank angle and end crank angle. Additional parameters like crank pin radius (a) and R (ratio of connecting rod lngth and crank pin radius) are defined.
An equation is formed to calculate ctual volume change using a geometric relation. Whenever this function would be called in the main program, it would calculate the volume (`V_(compression) and V_(expansion)`) in this case.
After `V_(compresion), P_(compression), V_(expansion) and P_(expansion)` are calculated, P-V diagram is plotted to get the actual P-V curve for Otto cycle.
Output
Conclusion
The P-V diagram for the Otto cycle is made and the thermal efficiency for the given state conditions comes out to be 0.5116.