Air Standard Cycle

Q) Based on the concepts you have learnt during the forward kinematics routine, go ahead and write code that can solve an otto cycle and make plots for it.

Here are the requirements

• Your code should create a PV diagram
• You should output the thermal efficiency of the engine.

ANS:

• Theory

To write a program that plots the PV diagram of an otto cycle, we should first understand the PV diagram. The diagram plots the relation between the pressure and the volume of the air fuel mixture inside an IC engine at different stages during the otto cycle. The first stage is isentropic compression. This occurs during the compression stroke of the engine where the piston travels upwards towards the top dead centre (TDC) to compress the air fuel mixture. Work is done by the piston on the air fuel mixture to compress it. The second stage is called adiabatic heat addition. In this stage, no transfer of work or heat is done. This happens at constant volume. The third stage is called isentropic expansion. This occurs during the power stroke of the otto cycle where the piston travels downwards to the bottom dead centre (BDC) as work is done on it by the air fuel mixture. The last stage is called adiabatic heat rejection. This also happens at constant volume.

With some inputs of temperature and pressure, we can calculate the required information for each stage successively. Initially, we can take the temperature at stage 1, t1, to be 500K and the corresponding pressure, p1 to be 1 atm which is equal to 101325 Pa. These variables have been called 'State Point 1 Variables' in the script. Gamma is another constant that is taken as 1.4. The temperature, t3 is for the third stage and will be used later.

• Inputs

clear all
close all
clc

% Constants

gamma = 1.4;
t3 = 2300

% State Point 1 Variables

p1 = 101325;
t1 = 500;

The first thing we need to find is the engine's geometric parameters. This refers to the engine's bore, stroke, compression ratio, and the length of the connecting rod. These parameters can be assumed as per our wish and I have taken the values shown in the script below. Since, we are dealing with an Otto cycle, we need to choose a compression ratio between 8 and 14.

• Engine Geometric Parameters

% Engine Geometric Parameters

bore = 0.1;
stroke = 0.1;
cr = 12;
con_rod = 0.15;

Next, using this information, we need to now find the parameters of swept volume and clearance volume. The swept volume is the volume that is displaced as the piston travels from the TDC to the BDC in one stroke length. The clearance volume is the volume between the cylinder head and the piston top when the piston is at the TDC position. The formulas to calculate both are shown below:

`V_s = Pi/4xxD^2xxL`

where Vs is the swept volume, D is the bore, and L is the stroke length.

`V_c = V_s / (cr - 1)`

where Vc is the clearance volume and cr is the compression ratio.

I have written these formulas in to the script as shown below:

% Calculation for Swept Volume & Clearance Volume

V_s = (pi/4)*bore^2*stroke;
V_c = V_s / (cr-1);

v1 = V_s + V_c
v2 = V_c

As you can see, I have also calculated 2 new variables, V1 and V2. V1 is the sum of the swept and clearance volumes and tells us the maximum volume of the engine cylinder. V2 is the clearance volume and represents the minimum volume of the cylinder. The ratio of V1 to V2 tells us the compression ratio of the engine cylinder.

After that, we can find the state point variables for state 2 which marks the end of the isentropic compression stage.

% State Point 2 Variables (End of Isentropic Compression)
% p2v2^gamma = p1v1^gamma
% p2 = p1 * (v1/v2)^gamma
% p2 = p1 * (cr)^gamma


p2 = p1 * (cr)^gamma;

In the isentropic compression stage, `P_2V_2^gamma = P_1V_1^gamma`.

In other words, either side of the equation is equal to a constant. Knowing the values of P1, V1, V2 and gamma, we can determine P2, the pressure at the end of isentropic compression.

Because the ratio of V1 to V2 is the compression ratio, I used the variable of 'cr' in the code which was defined earlier in the engine geometric parameters.

Next, we need is the temperature at the end of stage 2, T2. This can be determined by using the Ideal Gas Law:

`PV = nRT`

Since the air fuel mixture is the same and in the same environment, PV/T is a constant because we can disregard the number of moles, n and the molar gas constant, R in the ideal gas equation.

So, we can say mathematically that `P_1V_1/(T_1) = P_2V_2/(T_2)`.

By reorganizing, we can find T2. The code has been entered in the script below:

% Calculation of t2 (using Ideal Gas Law)
% p1v1 / t1 = p2v2 / t2
% t2 = p2v2 * t1 / p1v1

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

% State Point 3 Variables (End of Constant Volume Heat Addn.)

v3 = v2

Next, we can move on the state variables for state point 3 which represents the end of the adiabatic heat addition. Since this stage happens at constant volume between stage 2 and 3, we can say that V3 = V2. So, we can now determine the other variables for this state.

NOTE: The value for T3 was not pre determined, so we can assume a value for it. However, T3 needs to be the highest temperature of all the 4 stages as it has the highest pressure value. Since pressure and temperature are proportional, T3 must be greater than the other temperature values. So T3 is created as a variable and defined as being equal to 2300 in the constants section of the script shown in the script above.

% Calculation of p3 (using Ideal Gas Law)
% p3v3 / t3 = p2v2 / t2 
% p3 = p2v2 * t3 / t2 * v3      |      v3 = v2
% p3 = p2 * t3 / t2

p3 = p2*t3 / t2;

Once again, we will use the Ideal Gas Equation to determine P3. Using the same principles as before, and knowing v2 = v3 (we can cancel out the volumes), P3 is easily determined.

% State Point 4 Variables (End of Isometric Expansion)

v4 = v1;

Finally, we just need to find the information for the fourth and final stage of the Otto cycle which represents the end of the isentropic expansion stage. As we know that the transition between stage 4 back to stage 1 is the adiabatic heat rejection phase, this happens at constant volume. So we can say that V4 = V1.

Then, to find P4, we can use the same formula that we used before for determining P2. This is because we know that during the isentropic expansion stage `P_4V_4^gamma = P_3V_3^gamma`. Knowing V4, V3, and P3, we can find P2. This formula is coded in the script below.

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

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

At last, we just have to find the temperature at the end of stage 4, T4 using the ideal gas equation once more.

% Determining t4 (using Ideal Gas Law)
% p3v3 / t3 = p4v4 / t4
% t4 = p4v4 * t3 / p3v3

t4 = (p4 * v4 * t3) / (p3 * v3);

• Plotting

To plot the PV diagram, we just need to use the plot command, and enter the volume values we calculated as the x coordinates and their respective pressure values as the y coordinates. For instance, the point for the first stage starts at (V1, P1) and vice versa. I have added a marker style of a cross with a marker size of 8 and a black colour so they are clearer. The code for this is shown below.

% Plotting

figure(1)
hold on

plot (v1, p1, 'marker', 'x', 'MarkerSize', 8, 'color', 'black')

plot (v2, p2, 'marker', 'x', 'MarkerSize', 8, 'color', 'black')

plot (v3, p3, 'marker', 'x', 'MarkerSize', 8, 'color', 'black')

plot (v4, p4, 'marker', 'x', 'MarkerSize', 8, 'color', 'black')

This code will only plot 4 points representing the 4 different stages of the otto cycle. It will not give us the full picture. This is the output graph we receive when this code is executed:

• Piston Kinematics

In order to be able to plot the curved sections of the PV diagram during the isentropic compression and expansion stages, we need to understand the engine's piston kinematics. In other words, we need to be able to characterize the volume trace of the piston to see how it changes during the cycle. For this, I started a new script to write a separate program.

As I did before, I started with defining the engine geometric parameters using the same values.

% Engine Geometric Parameters

bore = 0.1;
stroke = 0.1;
cr = 12;
con_rod = 0.15;

The equations for calculating the swept and clearance volumes were also written in the script once more.

% Calculation for Swept Volume & Clearance Volume

V_s = (pi/4)*bore^2*stroke;
V_c = V_s / (cr-1);

• Volume Trace

In order to characterize the volume trace, we need to know how it varies. The volume displaced by the piston changes as the crankshaft rotates in motion. It rotates clockwise from 0 to 360 degrees. We can find the volume at an instant if we know the angle, but this will only tell us a single value. In order to characterize the relation, we need more values for the angle.

So we can generate an array for the angles and assign it to the variable, 'theta' using the linspace command. Ranging from a start crank angle of 0 to and an end crank angle of 180, I have chosen 100 equally spaced values:

theta = linspace (0,180,100);

To relate the volume trace to the crank angle, we need to use a formula which is written below:

`V/V_c = 1 + 1/2 (cr - 1)[R + 1 - costheta-(R^2 - sin^2theta)^(1/2)]`.

In this formula, a new variable is introduced called, 'R' which represents the ratio of the connecting rod length to the crank pin radius. The crank pin radius is another variable and is equal to the stroke length divided by 2.

So we can create 2 new variables, 'R' and 'a'. 'R' can be defined as the connecting rod length, 'con_rod' divided by the crank pin radius, 'a'.

Also we need to define the crank pin radius, 'a' by equating as being the 'stroke' divided by 2. These 2 new variables are included in the script in the inputs section as shown below.

a = stroke/2
R = con_rod / a

This formula is coded in the MATLAB script by splitting it up in to 3 parts with each part being assigned to a different variable named 'term1', 'term2' and 'term3'. This was done in order to make the calculation easier to perform and to reduce the risk of mathematical error.

Once the 3 terms are determined separately, another variable is created called 'V_trace', which will combine all the 3 variables and multiply it with the clearance volume, Vc, in order to calculate the volume trace.

This calculation will output an array of values for the volume trace with each value corresponding to a certain value for 'theta'. If you have chosen the 'theta' array to have 100 values like I have, then we will obtain 100 values for the volume trace. The code for the volume trace calculation is shown below:

% Coding of Volume Trace Eqn

term1 = 0.5*(cr - 1);
term2 = R + 1 - cosd(theta);
term3 = (R^2 - [sind(theta)].^2).^0.5;

V_trace = [1 + term1 * (term2 - term3) ] * V_c;

Error 1

Because we have an array for the 'theta' values, and we are using these values in a calculation to obtain the volume trace, we need to use the elemental operators correctly for the program to execute. For the variable, 'term3' there are 2 exponents in the calculation. So, if we do not use the operator (.), before the exponent 2 and also before the exponent 0.5, the code will not execute.

Initially, I had forgotten to enter the operators at all.

term3 = (R^2 - [sind(theta)]^2)^0.5;

This gave me the error statement:

After including the operators, the spacing was another issue.

term3 = (R^2 - [sind(theta)]. ^2). ^0.5;

The line of code shown above will not execute because there are spaces after the operator. There shouldn't be any gaps. This gave me another error statement:

Going back to the volume trace, when the code is correct and executed successfully, we will obtain an array of values for the volume trace. The screenshot below shows a few of the columns that are generated in the array for the volume trace.

Now that we have sufficient values to characterize how the volume changes with the crank angle, we can plot a graph to see how it actualy varies. This will act as a check to see if our code is correct and giving us the information we need.

I used the plot command shown below to plot the crank angle (x axis) against the volume (y axis).

plot (theta, V_trace)

This gave me the graph shown below as the output:

As we can see, the volume increases as the crank angle changes from being a minimum at 0 to being a maximum at 180 degrees. This means our array of volume trace values is correct. The changing of the crank angle from 0 to 180 degrees occurs during the expansion stroke. So naturally, the volume will be smallest at the start crank angle and highest at the end crank angle.

• Creating a Function

In order to perform the Piston Kinematics calculations more conveniently, we can create a function for determining the volume trace. This function would require us to input the necessary values for the calculations and would execute the code to give us the array of volume trace all with one line of code. This is shown below:

function [V_trace] = Piston_Kinematics(bore, stroke, con_rod, cr, start_crank, end_crank)

NOTE: This line of code should be at the top of the script before anything else and an end command should be given at the very end of the script.

I chose the function's name as 'Piston_Kinematics' as this name should match the name of the script. The function calculates the volume trace, so we need to make sure the same variable name for the volume trace is used in the function code. Before, I used 'V_trace' as the name of the variable when performing the calculations. So the same name is used.

The inputs for the function are the bore, stroke, connecting rod length, compression ratio, the start crank angle and the end crank angle.

So, we can remove the lines of code for the bore, stroke, compression ratio and the connection rod length from the inputs since they are no longer fixed values. They can take any value we want to give us the required results for the volume trace.

• Calling the Function

Calling the function means to execute the function to do what we want it to do. To do this, we just have to type the name of the function in the command window. After the name, we have to insert open parentheses and enter the required inputs and end with a closed parentheses.

NOTE 1: If there are 5 inputs required (as stated in the function line code), we need to enter 5 inputs. No more, no less.

NOTE 2: Because we have defined the start and end crank angles as 0 and 180 degrees in the function, the same values don't need to be entered in the 'theta' array. Before, the 'theta' array was directly defined with start and end crank angles.

theta = linspace (start_crank, end_crank, 100);

Now, the 'theta' array has been coded as a template. Whatever start and end crank angle we input in the function, will follow through with the 'theta' array to give us 100 outputs just as before because I have still chosen 100 equally spaced values.

The function is called as shown below:

I have entered the same values for the bore, stroke, connecting rod length and the compression ratio to show that we still obtain the same outputs for the volume trace as when before the function was implemented. So we now know that the function works.

Notice how the values for 'a' and 'R' are also calculated correctly from the inputs to the function.

• Integration with Otto Cycle

Now, we can copy the function code in to the script for the Otto cycle. Because the volume trace values we received were for the angle range of 0 to 180 (expansion stage), we can assign the array of values to a new variable, 'V_expansion'. The line below can be copied and pasted in the Otto cycle script after the T4 calculation.

% Calculation of P_expansion

V_expansion = Piston_Kinematics(0.1, 0.1, 0.15, 12, 0, 180);

From the above code, it is seen that we want to calculate 'P_expansion' which are the corresponding array of values for the pressure during the expansion stage. We can calculate this array using the 'V_expansion' array. But we also need another piece of data.

NOTE: We know that at the end of the isentropic compression and expansion stages, the product of the pressure and the volume raised to gamma, is equal to a constant. In other words:

`PV^gamma = c`.

Using the P3 and V3 values, we can find the constant, c2.

NOTE: This line of code below should come before the 'V_expansion' function code.

% Determining constant, c_2

c_2 = p3*v3^gamma;

Now we can use it and the 'V_expansion' array to find the 'P_expansion' array. By reordering the equation so that:

`P = c / V^gamma`.

P_expansion = c_2 ./ V_expansion.^gamma;

NOTE: The above line of code must come straight after the line of function code.

Error 2

As with Error 1, the same error with element operator occured again. 'V_expansion' is an array, so when using it in a calculation, we need to use the (.) operator before the exponent and also before the division symbol.

• P_expansion array

This screenshot shows the generation of the 'P_expansion' array from the calculation. As we have 100 values in the 'V_expansion' array, we will have 100 corresponding values for the pressure.

• P_compression

Similarly, we repeat the same process to find the pressure array during the compression stage.

Firstly, we change the piston kinematics function inputs so that the start and end crank angles are 180 and 360 degrees respectively (compression stage). Second, we find the respective constant during this stage.

NOTE: This line of code appears after the calculation of P2.

% Determining constant, c_1

c_1 = p2*v2^gamma;

After this, we can copy the function code below it but this time assigning the output values of the function under a new variable, 'P_compression'.

% Calculation of P_compression

V_compression = Piston_Kinematics(0.1, 0.1, 0.15, 12, 180, 360);

P_compression = c_1./ V_compression.^gamma;

• P_compression array

We will have 100 values for the pressure during the compression stage also.

Error 3

The last error I found was in the 'P_expansion' and the 'P_compression' sections. I had written the code for the calculation of 'P_expansion' and 'P_compression' before their respective function codes. This will not work since we need the functions to execute first to obtain the volume traces 'V_expansion' and 'V_compression'. Then these values will be used to find the pressure during compression and expansion.

This is the error output if the function code for 'V_compression' and 'V_expansion' appear after the pressure calculation.

• Plotting PV Diagram

Now that we know how the volume changes during the compression and expansion stages, we can plot their values on to the PV diagram.

% Plotting

figure(1)
hold on

plot (v1, p1, 'marker', 'x', 'MarkerSize', 8, 'color', 'black')
plot (V_compression, P_compression, 'linewidth', 2, 'color', 'red')
plot (v2, p2, 'marker', 'x', 'MarkerSize', 8, 'color', 'black')
line ([v2 v3], [p2 p3], 'linewidth', 2, 'color', 'black')
plot (v3, p3, 'marker', 'x', 'MarkerSize', 8, 'color', 'black')
plot (V_expansion, P_expansion, 'linewidth', 2, 'color', 'blue')
plot (v4, p4, 'marker', 'x', 'MarkerSize', 8, 'color', 'black')
line ([v4 v1], [p4 p1], 'linewidth', 2, 'color', 'black')

title ('PV Diagram of an Otto Cycle')
xlabel ('Volume (m^3)')
ylabel ('Pressure (Pa)')

From the above code, after the lines of code for the plotting of the first and third point, I entered a plot command to plot the V_compression and P_compression and the V_expansion and P_expansion respectively. These 2 lines of code will give us the curves between stages 1 and 2 and stages 3 and 4.

The straight line during adiabatic heat addition and heat rejection, is plotted using the 2 line commands. The first line command comes after the plotting of the second point and the second line command comes after the plotting of the fourth point.

The line commands simply connect the x and y coordinates of the first and fourth coordinates as well as that of the second and third. Line commands were used for these 2 stages since they happen at constant volume. So a vertical line is required to illustrate these stages.

Finally, a title and labels were added for presentation.

This is the PV diagram I obtained:

2) Thermal Efficiency

To find the thermal efficiency of this Otto cycle, we need to find the specific internal energy, u (kJ/kg) which is dependent on the temperature of air in Kelvin. From all the temperature values we have for all the 4 stages, we can find their corresponding values for the specific internal energy from a data book. In my case, I found the values online. Whenver the temperatures were in between ranges, I took average values.

For T1 = 500K, u1 = 359.8kJ/kg

For T2 = 1351K, u2 = 1068kJ/kg

For T3 = 2300K, u3 = 1970kJ/kg

For T4 = 851K, u4 = 633.5kJ/kg

These values for specific internal energy are entered as variables and defined by the values found.

% Specific Internal Energy at each respective temperature

u1 = 359.8
u2 = 1068
u3 = 1970
u4 = 633.5

To find the thermal efficiency, the formula is:

`eta = 1 - (u_4 - u_1)/(u_3 - u_2)`

So using the variables, I entered the formula to find the thermal efficiency to which I assigned the variable, 'eta'. The thermal efficiency was calculated as a percentage also. The coding is shown below and the output value along with it.

% Calculation of Thermal Efficiency, eta

eta = 1 - [(u4 - u1) / (u3 - u2)]
eta_percentage = eta*100

As we can see, the thermal efficiency came to be calculated as 0.6966 or 69.7%.

Conclusion

The air standard cycle was solved using the Forward Kinematics routine on MATLAB and its PV diagram was plotted. The thermal efficiency of the engine in the cycle was also determined.