Matlab - 'Isothermal' Flow in a Petrol Engine Inlet Tract |
For high velocity flows of gases in nozzles losses are often small and the flow can be approximated as a reversible adiabatic process, commonly called an isentropic (constant entropy) process. In this process there are assumed to be no heat transfer gains or losses and no friction losses.
However flow across an the edge of a (throttle) plate will induce turbulence and the above model is not accurate. A better model is the 'isothermal orifice' where the assumption is made that the behaviour of the fluid can be separated as follows:
2. Approach to Modeling
There are a number of ways of tackling this problem but the method suggested here is
to develop a Matlab calculation (such as the one shown below) in which most parameters need only setting up once.
To run a calculation the (net) power output, the percentage that the throttle is
open, the assumed pressure ratios across the throttle plate (p2p1) and across the
inlet valve (p3p2) are entered when prompted.
The calculations are done and the three values of calculated air mass flow rate examined,
and changes made to the p2p1 and p3p2 pressure ratios and the analysis re-run, until the three
air mass flow rates calculated are the same.
3. Matlab Listing
% Isentropic flow inlet tract - 10th August 2005
r=287; % Universal gas constant
t1=293; % Initial temperature, ambient in deg K
Cd=0.6; % Assumed coefficient of discharge around throttle plate
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 percentage 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 across the throttle plate, p2/p1 ');
p2=p2p1*p1; % Pressure downside of throttle
p3p2=input('Enter trial pressure drop ratio across 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
ro2=p2/(r*t2); % Final density
g1=airmflow; % Initial air mass flow to provide gross power.
m2dot=sqrt(2*p2p1*(1-p2p1))*Cd*opencsa*p1/(sqrt(v1*r)); % Calculated air mass flow rate after throttle
p3=p3p2*p2p1*p1; % Absolute value of p3
invalvedia=30; % Inlet valve diameter in mm
openheight=8; % 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
m3dot=sqrt(2*p3p2*(1-p3p2))*Cd*valveopenarea*p2/(sqrt(v1*r)); % Calculated air mass flow rate after throttle
pw=(p1-p3p2*p2p1*p1)*airmflow/(ro1*1000); % Pumping work, kW, 10th August 2005.
a=[g1, m2dot, m3dot];
disp('Air mass flow rates, kg/s ');
disp(a);
disp('Pumping work, kW ');
disp(pw);
Several of the values are standard and can probably be left unaltered ie: 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 power and percentage of throttle opening, p2p1 and p3p2 will have to be varied until
consistent values for air mass flow rate are obtained.
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:
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.
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.