Code to solve problem on Otto cycle

Code to solve problem on Otto cycle

Based on the challenge for week 2.1 based on air standard cycles and within that, Otto cycle as the cycle of consideration, here's a code to solve problems based on Otto cycle with a given set of inputs.

So, basically here we need to consider an engine with some specific characteristic dimensions(bore,stroke,etc) to work upon. For that matter, I have considered the engine of Bugatti Veyron, which is a W16 engine with a BoreXStroke of 86X86.5 mm and a compression ratio of 9:1. I have considered its connecting rod to be 100mm long. The other set of inputs are:

Intial pressure & temp(correspondong to point 1) - 1 atm & 400K respectively,

Adiabatic index - 1.4,

Maximum temperature (corresponding to point 3) - 2300K

Understanding the problem:

To begin with, we need to specify the engine dimensions and other initial conditions as a set of inputs. Having defined those, we'll have to define swept volume and clearance volume. Swept volume is given by,

v_swept=(pi)/4Xbore^2Xstroke

whereas clearance volume is given by,

v_clearance=v_swept/(cr-1).....cr is compression ratio

Then we'll have to calculate variable values at different state points.

For state point 1, pressure and temperature are known whereas volume corresponding to it is the total volume which is the sum of swept and clearance volume.

For state point 2, volume corresponds to the clearance volume whereas pressure can be calculated from the equation,

p1*v1^(gamma)= p2*v2^(gamma).....for an adiabatic/isentropic process

Temperature can be calculated from the equation,

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

For state point 3, volume is equal to volume at state point 2 as 2-3 is a constant volume process. Temperature at this point is the maximum as combustion takes place around this point, which has been mentioned. Pressure for point 3 can be calculated from,

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

For state point 4, volume is equal to volume at point 1. Pressure can be calculated from,

p3*v3^(gamma)= p4*v4^(gamma).....for an adiabatic/isentropic process

whereas temperature can be calculated from,

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

Now, if we move on to plotting with these conditions, we'll get all the state points connected about by straight lines, which shouldn't be the case. As an adiabatic process traces a curve under the p-v plot, we'll have two curves in our problem: points 1-2 and points 3-4. Now, inorder to trace these curves, we'll have to define an array carrying various points of the curve so as to trace upon them. For this matter we can either define 2 separate variables or define a function and call upon that function twice for the two arrays. For the given problem, I have defined a function namely engine_kinematics. In our code, I have called upon this function twice to define two arrays for compression volume and expansion volume. Also, the pressure values corresponding to all these points can be obtained from,

p*v = constant

and the value of this constant can be obtained from

constant = p1*v1.......for compression pressure

&

constant = p3*v3.......for expansion pressure

Now, for points 1-2 we'll have to name the arrays as compression volume and compression pressure whereas for points 3-4 they'll be expansion volume and expansion pressure, respectively. Having done this, we can move on to plotting all the variables to get the points, lines and curves to form a Otto cycle plot.

Lastly, we have to calculate the thermal efficiency of the cycle as well. It is given by,

eta(thermal) = 1-(t1/t2)

or, eta(thermal) = 1-1/(cr^(gamma-1)).

The Code

% Code to solve problem on otto cycle

% Inputs
bore = 0.08; %m
stroke = 0.086; %m
con_rod = 0.1; %m
cr = 9; %compression ratio
v_swept = (pi)/4*bore^2*stroke; %m^3
v_clearance = v_swept/(cr-1); %m^3
p1=101325; %Pa
t1=400; %K
gamma=1.4; %adiabatic index
t3=2300; %K
v1=v_swept+v_clearance
v2=v_clearance

% State variables at point 2
p2 = p1*(cr)^gamma
t2 = (p2*v2*t1)/(p1*v1)

c1 = p1*v1.^gamma;
v_compression = engine_kinematics(.08,.086,0.1,9,180,0);

p_compression = c1./v_compression.^gamma;

% State variables at point 3
v3=v2
p3 = p2*(t3/t2)

c2 = p3*v3.^gamma;
v_expansion = engine_kinematics(.08,.086,0.1,9,0,180);

p_expansion = c2./v_expansion.^gamma;

% State variables at point 4
v4=v1
p4=p3*(1/cr)^gamma
t4 = (p4*v4*t3)/(p3*v3)

% Plotting
hold on
plot(v1,p1,'*','color','r')
plot(v_compression,p_compression,'linewidth',3,'color','r')
plot(v2,p2,'*','color','r')
line([v2 v3],[p2 p3],'linewidth',3,'color','g')
plot(v3,p3,'*','color','r')
plot(v_expansion,p_expansion,'linewidth',3,'color','r')
plot(v4,p4,'*','color','r')
line([v4 v1],[p4 p1],'linewidth',3,'color','g')

% Thermal efficiency
eta_th = 1-(t1/t2) %or
eta = 1-1/(cr^(gamma-1))


%%
% Functon engine_kinematics

function [V] = engine_kinematics(bore,stroke,con_rod,cr,start_crank,end_crank)
  
  a= stroke/2;
  R = con_rod/a;
  v_swept = (pi)/4*bore^2*stroke;
  v_clearance = v_swept/(cr-1);
  
  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_clearance;
  end

Outputs

Errors faced

1. Incorrectly defining a variable. While defining a variable the variable name must be equated against its value and not the other way round.

2. Using distinct notations while defining and using a variable. term 1 while defining and term1 while using it.