Air standard Cycle

Aim: To write MATLAB program that can solve an otto cycle and make plots for it.

Objective:

1. Creating a PV diagram.
2. Calculating the thermal efficiency of the engine.

Description:

To make a plot of a P-V diagram we need to understand it first. The P-V diagram of an Otto cycle represents the pressure and volume of air-fuel mixture inside an internal combustion engine.

The Otto cycle consists of four processes:-

Process 1-2: Isentropic Compression: In this process the piston moves from BDC to TDC and Air undergoes isentropic compression.

Process 2-3: Constant volume heat addition: In this process, heat is added at constant volume and piston rests at TDC.

Process 3-4: Isentropic Expansion: In this process piston moves from TDC to BDC and air undergoes isentropic expansion.

Process 4-1: Constant volume heat rejection: In this process, the heat is rejected at constant volume and piston rests at BDC.

To find the volume and temperature change behavior at each process we need to first understand the engine kinematics.

Engine geometric parameters:

1. Bore: Diameter of each cylinder.
2. Stroke: Movement of engine cycle from BDC to TDC and vice versa
3. Connecting Rod Length: Distance between crankshaft pin and piston pin.
4. Compression Ratio (cr): Ratio of volume at BDC and volume at TDC.

1. Function code for Engine Kinematics:

function [V] = engine_kinematics(bore, stroke, con_rod, cr, start_crank, end_crank)

a = stroke/2;
R = con_rod/a;

v_s = pi/4*bore^2*stroke;
v_c = v_s/(cr-1);

%V/V_c = 1+0.5(rc-1) [R+1-cosd(theta)-(R^2-sind(theta).^2]^0.5

theta = linspace(start_crank,end_crank,100);

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

1.1 Output: (Validation of function)

2. Code to visualize the P-V diagram of an Otto cycle:

% P-V Diagram of an Otto cycle

clear all
close all
clc

%Inputs
gamma = 1.4;
t3 = 2300;

%Engine geometric parameters
bore = 0.1;
stroke = 0.1;
con_rod = 0.15
cr = 12;

%Calculating the swept volume and  the clearance volume

V_swept = pi/4 * bore^2 * stroke;

V_clearance = V_swept/(cr-1);

v1 = V_swept + V_clearance
v2 = V_clearance

%State variables at state point 1
p1 = 10325
t1 = 500

%State varibales at state point 2
%p2v2^gamma = p1v1^gamma
%p2 = p1 * (v1/v2)^gamma = p1*cr^gamma

p2 = p1*cr^gamma;

%Calculating temperature at state point using Ideal gas law
%p1v1/t1 = p2v2/t2 | t2 = p2*v2*t1/(p1*v1)

t2 = p2*v2*t1/(p1*v1)

constant_c = p1*v1^gamma;

V_compression = engine_kinematics(bore, stroke, con_rod, cr, 180, 0)

P_compression = constant_c./V_compression.^gamma;

%State variable at point 3

v3 = v2

%p3v3/t3 = p2v2/t2 | p3 = p2*t3/t2

p3 = p2*t3/t2

constant_c = p3*v3^gamma;

V_expansion = engine_kinematics(bore, stroke, con_rod, cr, 0, 180)

P_expansion = constant_c./V_expansion.^gamma

%State variables at point 4

v4 = v1

%p3v3^gamma = p4v4^gamma | p4 = p3(v3/v4)^gamma

p4 = p3*(v3/v4).^gamma

%Plotting

figure(1)
hold on

plot(v1, p1, '*', 'color', 'r')
plot(V_compression, P_compression)
plot(v2, p2, '*', 'color', 'r')
line([v2 v3], [p2 p3], 'color', 'black')
plot(v3, p3, '*', 'color', 'r')
plot(V_expansion, P_expansion)
plot(v4, p4, '*', 'color', 'r')
line([v4 v1], [p4 p1], 'color', 'black')
title('P-V Diagram of an Otto cycle') 
xlabel ('Volume (m^3)')
ylabel ('Pressure (Pa)')

2.1 Output:

3. Code for calculating the thermal efficiency of the engine:

% Calculating the thermal Efficiency of an Otto engine

%Input
cr = 12
gamma = 1.4

%Thermal efficiency of an Otto cycle

Eta = 1 - [1/cr^(gamma-1)]

Eta_percentage = Eta*100

3.1 Output:

4. Explanation of Code

*Understanding the code to generate the P-V diagram of an Otto cycle:

% P-V Diagram of an Otto cycle

clear all
close all
clc

%Inputs
gamma = 1.4;
t3 = 2300;

%Engine geometric parameters
bore = 0.1;
stroke = 0.1;
con_rod = 0.15
cr = 12;

Here,

• Clear all, close all, clc command used to delete the variables created in workspace and command window before writing a new program.
• Inputs gamma and t3 already defined (assumed).
• Engine geometric parameters are already defined.

%Calculating the swept volume and  the clearance volume

V_swept = pi/4 * bore^2 * stroke;

V_clearance = V_swept/(cr-1);

v1 = V_swept + V_clearance
v2 = V_clearance

• Swept volume: It is the volume of the cylinder as the piston moves from TDC to BDC in one stroke.

Mathematically,

V_swept = (Pi/4)*bore^2*stroke

• Clearance Volume: It is defined by the volume between the cylinder head and the piston top when the piston is at TDC.

Mathematically,

V_clearance = V_swept/(cr-1)

Here two variables, v1 and v2. v1 is the sum of swept volume and clearance volume which represents the maximum volume of the cylinder. And v2 is the clearance volume which represents the minimum volume of the engine cylinder.

%State variables at state point 1
p1 = 10325
t1 = 500

Here, the state variables are predefined. Assuming the p1 = 10325 Pa and t1 = 500k.

%State varibales at state point 2
%p2v2^gamma = p1v1^gamma
%p2 = p1 * (v1/v2)^gamma = p1*cr^gamma

p2 = p1*cr^gamma;

%Calculating temperature at state point using Ideal gas law
%p1v1/t1 = p2v2/t2 | t2 = p2*v2*t1/(p1*v1)

t2 = p2*v2*t1/(p1*v1)

constant_c = p1*v1^gamma;

V_compression = engine_kinematics(bore, stroke, con_rod, cr, 180, 0)

P_compression = constant_c./V_compression.^gamma;

State variables at state point 2: It is the end of the isentropic compression stage.

In this stage, p2v2^gamma = p1v1^gamma

With the help of the above equation, we can determine the pressure at stage 2 (p2)

p2 = p1*(v1/v2)^gamma

The ratio of v1 and v2 is compression ratio (cr)

So, p2 = p1*(cr)^gamma

We need to calculate the temperature at stage 2 (t2).

Using Ideal gas law,

Pv = nRT

Since the air-fuel mixture is the same in stage 1 and 2 we can say p1v1 = p2v2.

N and R are constant.

Mathematically, p1v1/t1 = p2v2/t2

So, t2 = p2*v2*t1/(p1*v1)

Next,

V_compression is the function I created earlier, which I will explain later. It is the volume trace in the compression stage. In the compression stage, the starting crank angle is 180 and end crank angle is 0 because during compression the volume decreases.

In isentropic compression, p1*v1^gamma is constant which is represented by constanct_c.

In other words, P_compression*V_compression^gamma = constant_c

So, we have already know the V_compression and constrant_c so we can compute the P_compression by the above expression.

P_compression = constant_c/V_compression^gamma

%State variable at point 3

v3 = v2

%p3v3/t3 = p2v2/t2 | p3 = p2*t3/t2

p3 = p2*t3/t2

constant_c = p3*v3^gamma;

V_expansion = engine_kinematics(bore, stroke, con_rod, cr, 0, 180)

P_expansion = constant_c./V_expansion.^gamma

State variables at state point 3: It represents the endpoint of the constant volume compression process or starting point of the expansion process. And it happens at process 2 and 3

So, v3 = v2

or, p3v3/t3 = p2v2/t2

So, p3 = p2*v2*t3/(t2*v3)

Since v3 = v2

p3 = p2*t3/t2

Next,

V_expansion is the function I created earlier, which I will explain later. V_expansion is the volume trace in the expansion stage. In the expansion stage, the starting crank angle is 0 and end crank angle is 180 because the volume increases in this stage.

In isentropic expansion, p3*v3^gamma is constant which is represented by constanct_c.

In other words, P_expansion*V_cexpansion^gamma = constant_c

So, we have already know the V_expansion and constrant_c so we can compute the P_expansion by the above expression.

P_expansion = constant_c/V_expansion^gamma

%State variables at point 4

v4 = v1

%p3v3^gamma = p4v4^gamma | p4 = p3(v3/v4)^gamma

p4 = p3*(v3/v4).^gamma

State variables at state point 4: Isentropic expansion stage the volume of stage 4 is equal to stage 1.

Hence, v4 = v1

During isentropic expansion stage p4 = v4

p4v4^gamma = p3v3^gamma

Or, p3v3^gamma = p4v4^gamma

By using above expression we can calculate the p4.

P4 = p3*(v3/v4)^gamma

Now that we have determined the value of p1, p2, p3, p4 we can generate a plot out of it.

%Plotting

figure(1)
hold on

plot(v1, p1, '*', 'color', 'r')
plot(V_compression, P_compression)
plot(v2, p2, '*', 'color', 'r')
line([v2 v3], [p2 p3], 'color', 'black')
plot(v3, p3, '*', 'color', 'r')
plot(V_expansion, P_expansion)
plot(v4, p4, '*', 'color', 'r')
line([v4 v1], [p4 p1], 'color', 'black')
title('P-V Diagram of an Otto cycle') 
xlabel ('Volume (m^3)')
ylabel ('Pressure (Pa)')

• Using figure command we can display the graphical output.
• Using plot(v1, p1, '*', 'color', 'r') command to plot the volume and pressure at stage 1 with color red. And the same for stage 2, stage 3 and stage 4.
• Using plot(V_compression, P_compression) and plot(V_expansion, P_expansion) commands to show the volume and pressure at compression and expansion stage 2 and 3.
• Using line([v2 v3], [p2 p3], 'color', 'black') and line([v4 v1], [p4 p1], 'color', 'black') commands to join the volume and pressure at stage 2 and 3 and stage 4 and 1 respectively.
• Using title('P-V Diagram of an Otto cycle') command to give title on the graph.

• Using xlabel ('Volume') and ylabel ('Pressure') commands to give name to the x and y axis.

*Now Understanding the code to calculate the efficiency of an Otto engine:

% Calculating the thermal Efficiency of an Otto engine

%Input
cr = 12
gamma = 1.4

%Thermal efficiency of an Otto cycle

Eta = 1 - [1/cr^(gamma-1)]

Eta_percentage = Eta*100

The thermal efficiency of an Otto engine is,

η = 1 - [1/cr^(gamma-1)]

We already know the clearance ratio and gamma, so we can calculate the efficiency of an Otto cycle using the above expression. We can determine the percentage by multiplying it to 100.

*Now Understanding the Function code:

function [V] = engine_kinematics(bore, stroke, con_rod, cr, start_crank, end_crank)

a = stroke/2;
R = con_rod/a;

v_s = pi/4*bore^2*stroke;
v_c = v_s/(cr-1);

In order to connect the state variables, we need to understand piston kinematics and how the volume changes at each cycle. For this, I have created a separate function[V] called engine_kinematics. A function is a code that can be re-used by giving the correct parameter in a different code.

Here, the function [V] is volume trace named engine_kinematics and includes the following parameters: bore, stroke, connecting rod length, clearance ratio, starting crank angle and end crank angle.

a is the crankpin radius which is stroke divided by 2.

R is the ratio of connecting rod length to the crankpin radius.

We have already calculated the swept volume (v_s) and clearance volume (v_c) before.

%V/V_c = 1+0.5(rc-1) [R+1-cosd(theta)-(R^2-sind(theta).^2]^0.5

theta = linspace(start_crank,end_crank,100);

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

The theta is starting crank angle from stage 1 to the end of the crank angle in one complete cycle (stage 4) which is 360. To generate the array for the angle I have assigned it to linspace command and chosen to give 100 values at equal spaces.

To relate the volume trace [V] to the crank angle we need to use the given formula:

V/V_c = 1+0.5(rc-1) [R+1-cosθ-(R^2-sinθ.^2]^0.5

We are using this formula that determines how the volume inside the combustion chamber changes. And we are breaking it into term1, term2 and term3 to avoid mathematical errors.

So, the volume trace is going to be V = [1+term1*(term2-term3)] and multiplied by v_c

So, this exact function we have used in our code to visualize the P-V diagram of an Otto cycle.

Conclusion:

The Air standard cycle was solved by the forward kinematics routine and the P-V diagram generated. The thermal efficiency was also calculated.