Adiabatic Compressible Flow in a Converging Duct

[Twigg]

TL;DR Summary

Conventional wisdom says that when a subsonic gas flow accelerates through a converging duct, it can only achieve Mach 1 at the point of minimal area (throat). The math makes sense but I don't understand the physics that makes this happen on the local scale. I also have incomplete simulation results that I don't understand.

For compressible flow in a duct, mass conservation combined with reversibility (no entropy change) implies $$(1-\frac{u^2}{c^2})\frac{du}{u} = -\frac{dA}{A}$$
where u is the flow velocity of the gas, c is the speed of sound in the gas, and A is the area of the duct. I am assuming a calorically perfect, ideal gas (ideal gas law holds for u = 0 and constant heat capacities, but $γ=\frac{c_p}{c_v}$ need not be 5/3, as it would be for a true ideal gas with no internal structure). Based on these assumptions, the speed of sound is $c=\sqrt{γRT}$ where R is the specific gas constant (universal gas constant divided by molar mass, $R = \frac{\bar{R}}{M}$). Just to be extra clear, this flow is adiabatic, not isothermal, so the speed of sound is not constant.
Now, the argument I see in the textbooks is that when a flow gets accelerated to Mach 1, the left hand side of the above equation goes to zero, therefore $dA=0$. The textbooks interpret this to mean that Mach 1 only ever occurs at the throat. I understand the argument on mathematical grounds. Originally, my physical intuition was that as the gas flows it compresses, and so the temperature rises. As the temperature rises, the speed of sound increases. As a result, the flow can only approach Mach 1 asymptotically, with the speed of sound and flow velocity increasing at almost equal rates. However, this does not agree with my numerical results.

In my simulations, I see the flow velocity approach Mach 1 sharply, followed by pathological behavior of the numerics once Mach 1 is reached. Extrapolating the non-pathological regime suggests that the flow really does reach Mach 1.

I have also done this simulation with a fixed step solver I wrote in C++, but the code is not easy to read and it calls several external tools, so it'd be a real hassle for anyone to either try to deciper it or try to run it on their own system. I'll just say that the qualitative results are identical, namely that the Mach number approaches 1 sharply.

R = 8.314 / (44.01E-3); %specific gas constant for propane
g = 1.13; %heat capacity ratio

r1 = 0.5 / 2 * 2.54e-2; %initial radius in meters (0.5 inches)
r2 = 0.031 / 2 * 2.54e-2; %final radius in meters (0.031 inches)
A1 = pi*r1^2;
A2 = pi*r2^2;

p1 = (30 + 14.7)/14.7 * 101500; %initial pressure in bar (30psi gauge at inlet)
T1 = (40 -32) * (5/9) + 273; %initial temperature in Kelvin (40 Fahrenheit)
u1 = 10; %initial flow velocity in m/s

T0 = (R*T1 + 0.5 * (1-1/g)*u1^2)/R; %stagnation temperature, from Bernoulli equation & ideal gas law
p0 = (T0/T1)^(g/(g-1)) * p1; %stagnation pressure, from adiabatic flow relations
rho0 = p0/R/T0; %stagnation density, from ideal gas law

rho1 = (T1/T0)^(1/(g-1)) * rho0; %initial density, from adiabatic flow relations
c0 = sqrt(g*R*T0); %speed of sound at stagnation point

dudA = @(A,u) - (real(u)/A) / (1 - (real(u)/c0)^2 * (rho0*real(u)*A/u1/A1/rho1)^(g-1));
meh = ode45(dudA,[A1,A2],u1);
% note: the use of "real" is just to force MATLAB to only consider the real part

u = meh.y; A = meh.x; 
rho = (((A1*u1) ./ A) ./ real(u)) * rho1; %local density, from mass conservation
p = rho .* (p0/rho0 - 0.5 * (1- 1/g) * real(u).^2); %local pressure, from Bernoulli equation
T = (rho/rho0).^(g-1) * T0; %local temperature from adiabatic flow relations
c = sqrt(g*R*T); %local speed of sound

%figure(1); clf;
%plot(A,real(u)./c)
%xlabel('Duct Area (m^2)')
%ylabel('Local Mach Number')

% Sanity check: local entropy
% Get the stagnation states for each local thermodynamic state
l_T0 = (R*T + 0.5*(1-1/g)*real(u).^2)/R;
l_p0 = (l_T0./T).^(g/(g-1)) .* p;

for n = 1:length(u)
  s(n) = CoolProp.PropsSI("S","T",double(l_T0(n)),"P",double(l_p0(n)),"Propane"); %lookup mass specific entropy
end

figure(1); clf;
subplot(2,1,1)
plot(A,real(u)./c)
xlabel('Duct Area (m^2)')
ylabel('Local Mach Number')
subplot(2,1,2)
plot(A,s/1000)
xlabel('Duct Area (m^2)')
ylabel('Specific entropy (kJ/(kg*K))')

I put the entropy in there so you know where to stop trusting the numerics.

 

[Twigg]

I've started having ideas about this. I'm starting to wonder if there's an absolute upper limit on the flow velocity into the duct at the inlet.

What I've started doing is using the compressible Bernoulli equation, the mass conservation equation, and the adiabatic flow equation (pressure vs. density), and seeing if they have any solution for my initial conditions and geometry. Turns out, they don't. There is no solution until the inlet velocity is less than approximately 1 m/s (given 30 psi gauge pressure at inlet, $\gamma$ = 1.13 for propane, 40F inlet temperature, 0.5 inch inlet diameter, and 0.031 inch outlet diameter). So now, I'm thinking that what's going on is that 10 m/s into a duct with these parameters does not give you a steady state flow, and the inlet velocity will drop in time until it reaches a speed (~1 m/s) where you get Mach 1 at the outlet. Does that sound reasonable?

 

[boneh3ad]

Consider what would happen if you had a point that reaches Mach 1 somewhere upstream of the throat. The walls are still converging, so does the flow downstream of that accelerate or decelerate? If the Mach number increased past that sonic point, then the walls are converging and it is contradictory since a supersonic flow should slow down under that condition. If the Mach number decreased past that point, then the flow would be subsonic but the walls are converging so it should accelerate. That's another contradiction. The only way you get past that is by requiring Mach 1 to occur only at the throat.

Regarding your numerics, I haven't had a chance to comb through them, but the bottom line is that not every random combination of inlet velocity and pressure ratios has a solution.

 

[Twigg]

This is starting to make a lot more sense. Thanks, boneh3ad!

 

[Klystron]

This description of transonic wind tunnel design with attention to the nozzles supports boneh3ad's explanation and may provide additional food for thought.

NASA Ames Unitary Plan Wind Tunnel System combines subsonic, transonic and hypersonic regimes in a single complex using different nozzles and baffle designs. One can film transitions through the Schlieren windows in the 9'x7' test section.

General Schlieren image. 9'x7' windows (notice the model mount through left window.)

 

[Twigg]

This is definitely some useful food for thought Klystron! Some of the aerodynamics jargon is new to me though. Could you explain to me what a baffle is or point me to an introductory resource? I've only heard of baffles in heat exchangers. Is it the same principle?

Also, which part do you mean by the model mount?

Lastly, in the Schlieren image on the left, based on all the sharp changes in index, isn't that a bunch of shock waves? How should I compare that situation to adiabatic internal flow? Thanks for bearing with me

 

[Klystron]

This is definitely some useful food for thought Klystron! Some of the aerodynamics jargon is new to me though. Could you explain to me what a baffle is or point me to an introductory resource? I've only heard of baffles in heat exchangers. Is it the same principle?

Thank you for sharing my enthusiasm. Not sure of the principle. I would rather leave those explanations to our aerodynamicists. Baffles in this context refer to guide structures built into the interior wind tunnel walls that perform a variety of functions to optimize smooth laminar air flow through the test section.

Also, which part do you mean by the model mount?

Unlike their enormous cousins -- the full-scale wind tunnels that can enclose an entire aircraft fuselage -- the Unitary Plan tunnels test scale models, often fabricated at a nearby facility. The model mount partly visible through the circular Schlieren window of the 9x7' hypersonic tunnel is a metal structure installed normal to the tunnel floor with a narrower angled strut pointing in the horizontal direction of the air flow. Engineers attach scale models to the horizontal member; both model and mount loaded with a variety of sensors, transducers and gauges.

My task as software engineer included programming computers to sample, compare, collate and collect data during operations. The model mounts were themselves marvels of engineering capable of measuring even slight movements as the models react to the intense air flow.

Lastly, in the Schlieren image on the left, based on all the sharp changes in index, isn't that a bunch of shock waves? How should I compare that situation to adiabatic internal flow? Thanks for bearing with me

Again I defer to specialists for accurate explanations. I understand that actual transonic and hypersonic effects are not visible to the unaided eye, particularly in the rarefied atmosphere at high altitudes and induced in pressurized closed wind tunnels; hence the use of Schlieren (literally "streamer") techniques to visualize air movement.

As analogy, one can compare a pressurized wind tunnel to the cloud chambers used to detect subatomic particles. We cannot see most particles* but can photograph their tracks in the cloud chambers.

*According to some physics textbooks the dark adapted human eye can detect photons.