Matlab - Compressible Flow in a Petrol Engine Inlet Tract

1. Introduction
The Matlab Simulink model previously developed that looked at the air flow in the inlet tract, made the assumption that the air flow was in-compressible. This is a convenient simplification as the Bernoulli equation can be applied and as the density is assumed constant the calculations are straightforward. For a more accurate approach to this application the compressible nature of air should be considered. The flow is subsonic but there are three conditions that need to be satisfied as the air flow passes
first: the throttle plate and
secondly: the inlet valve:

2. Approach to Modeling
There are a number of ways of tackling this problem but the method suggested here is to set up a Matlab calculation into which initial conditions are entered and assumed pressure ratio are also put in. The calculations are done, the answers examined, ie the values for the 3 parameters listed above and changes made to the appropriate parameters and the analysis re-run until the values are very close to remaining constant,
ie: g1 almost equal to g2 almost equal to g3
eb1 almost equal to eb2 almost equal to eb3
e1 almost equal to e2 almost equal to e3

3. Matlab Listing
% Compressible flow inlet tract - 15th March 2005
k=1.412; % Ratio of specific heats, cp/cv
r=287; %Universal gas constant
t1=293; % Initial temperature, ambient in deg K
km12=(k-1)/2;
km1k=(k-1)/k;
invkm1k=1/km1k;
maxp=input('Enter power for calculation in watts ');
tractcsa=0.05*0.06; % Inlet tract cross section
openper=input('Enter the throttle opening as a percentatge of full open ');
opencsa=openper*tractcsa/100; % Gives area of inlet tract opened by throttle
p1=100000; % Initial pressure, ambient.
p2p1=input('Enter trial pressure drop ratio ratio accross the throttle plate, p2/p1 ');
p2=p2p1*p1; % Pressure downside of throttle
p3p2=input('Enter trial pressure drop ratio accross the inlet valve, p3/p2 ');
effic=0.3; % Nominal efficiency
numinvalve=1; % Nimber of inlet valves per cylinder
numcyl=4; % Number of cylinders
grossmaxp=maxp/effic; % Gross maximum power
petrolmflow=grossmaxp/46000000; % Petrol mass flow rate
airmflow=petrolmflow*14.7; % Air mass flow rate, kg/s
ro1=p1/(r*t1); % Initial density
v1=airmflow/(ro1*tractcsa); %Initial velocity in the inlet tract
eb1=v1*v1/2+invkm1k*p1/ro1; % Initial bernoulli energy
nm1=v1/(sqrt(1.412*r*t1)); % Initial mach number
nm22=(1+km12*nm1*nm1-exp(km1k*log(p2/p1)))/(km12*exp(km1k*log(p2p1))); % final mack number, squared
nm2=sqrt(nm22); % Mach number going past throttle
t2=t1*(1+km12*nm1*nm1)/(1+km12*nm2*nm2); % Final temperature deg K
v2=v1*(nm2/nm1)*sqrt((1+km12*nm1*nm1)/(1+km12*nm2*nm2)); % Final velocity
ro2=p2/(r*t2); % Final density
eb2=v2*v2/2+invkm1k*p2/ro2; % Final Bernoulli energy
g1=airmflow; % Initial air mass flow
g2=ro2*v2*opencsa; % Calculated final air mass flow rate
v2av=g2/(ro2*tractcsa); % velocity after throttle
nm2av=v2av/(sqrt(1.412*r*t2)); % Mach no. beyond throttle plate
cp=1007; % Specific heat at constant pressure
e1=cp*t1+v1*v1/2; % Initial energy
e2=cp*t2+v2*v2/2; % Final energy Porsche 944 air flow, compressible flow, 90 kW
p3p2=0.96; %Pressure drop ratio accross the inlet valve
p3=p3p2*p2p1*p1; % Absolute value of p3
invalvedia=35; % Inlet valve diameter in mm
openheight=10; % Mean open height of inlet valve
valveopenarea=3.14159*invalvedia*openheight*numinvalve*numcyl/4000000; % Area of valve opening ,assumes valve open 1/4 of time
nm32=(1+km12*nm2av*nm2av-exp(km1k*log(p3/p2)))/(km12*exp(km1k*log(p3p2))); % final mack number, squared
nm3=sqrt(nm32); % Final mach number
t3=t2*(1+km12*nm2av*nm2av)/(1+km12*nm3*nm3); % Final temperature deg K
v3=v2av*(nm3/nm2av)*sqrt((1+km12*nm2av*nm2av)/(1+km12*nm3*nm3)); % Final velocity
ro3=p3/(r*t3); % Final density
eb3=v3*v3/2+invkm1k*p3/ro3; % Final Bernoulli energy
g3=ro3*v3*valveopenarea; % Calculated final air mass flow rate
e2=cp*t2+v2*v2/2; % Final energy Porsche 944 air flow, compressible flow
e3=cp*t3+v3*v3/2; % Final energy
pw=(p1-p3p2*p2p1*p1)*airmflow/(1000*ro1); % Pumping work, kW, 15th March 2005.
a=[g1, g2, g3];
b=[eb1, eb2, eb3, e1, e2, e3];
disp('Air mass flow rates, kg/s ');
disp(a);
disp('Energies, Bernoulli and CpT etc ');
disp(b);
disp('Pumping work, kW ');
disp(pw);

Several of the values are standard and can probably be left unaltered ie: k, r, t1 and p1.

You should also note that the 'pumping work' calculation here is somewhat artificial as it takes no account of the engine speed. This should be incorporated into the calculation.

4. Running the Program, Values to Adjust and Saving the Program
Open Matlab and either type in the above in the command window or start a new .m file (this starts the m file editor - use menu: 'File' 'New') and type in the above. Save it as an .m file.
Initially numinvalve and numcyl, tractcsa, invalvedia and openheight will have to be set.
For a particular maxp, the following will have to be varied until consistent values (see above) are obtained: p2p1, opencsa and p3p2.

Before you log off from the work station, save the .m file in a directory on a floppy disc, or USB storage device. When you next come to use the program on a different workstation, you can either:

Copy the file from the floppy or USB device into the appropriate location - for Matlab 6.5 this is the MATLAB6p5/bin directory and for Matlab 7 this is the Matlab14/work directory.

Alternatively you can set the search path to the folder where you stored the file. This is done from the 'File' menu by choosing 'Set Path' and putting in the correct device and folder name - i.e: A:/folder.

Typing the name of the .m file at the Matlab prompt then runs the program.

5. Application to the Assignment.
Incorporating all of the above into a Simulink model would not be easy. However for the assignment you are in one part looking at full throttle acceleration and in another part full throttle after 1 second of opening the throttle. You could say that you will assume that as far as the pressure drops go you will assume that full throttle conditions apply. This gives a very small pressure drop between P1 and P2 (P2/P1 is 0.999 or more) with virtually all the pressure drop occurring between P2 and P3. You can run the above Matlab model for a range of power outputs with P1/P2 set as 0.999 (or more) and iterate out what the P2/P3 ratios are. You can then determine a relationship between P2, P3 engine speed and power to use in your Simulink model. You should keep in mind that you are trying to simulate a dynamic system which is not in equilibrium. The calculation done above assumes steady state conditions, so obviously this is not a particularly accurate representation of the actual dynamic system.

David J Grieve, 11th August 2005.