Piecewise and Evaluated If Not Equivalent (at least not within
- To: mathgroup at smc.vnet.net
- Subject: [mg104034] Piecewise and Evaluated If Not Equivalent (at least not within
- From: blamm64 <blamm64 at charter.net>
- Date: Fri, 16 Oct 2009 07:20:13 -0400 (EDT)
Hi All,
I've had Mathematica for a couple of years, but until about a month
and a half ago used it for nothing but deriving special "shape
functions" for a finite element analysis program (in Fortran 2003,
now) for solid mechanics.
Anyway, over the past month and a half I've been devoting as much time
as I can to steep myself in Mathematica syntax/semantics in an effort
to solve a pneumatic/hydraulic system problem, and as in most all
other projects I attempt I build the problem and its solution(s)
incrementally.
Briefly, the system is a pressurized bottle connected to a movable
piston through a torturous gas passage. The piston is resisted by
hydraulic fluid on the face opposite to the gas side. The bottle is
initially highly pressurized (3000 psig) and the analysis I'm
attempting is the bottle is suddenly "uncorked", gas (N2) flows
through the passage into the chamber housing the piston, driving the
piston against the hydraulic fluid as well as an external load on the
piston.
The first increment in the derivation was to simply model the
isentropic compressible flow out of the bottle through a small hole
out to atmosphere. I won't get into too much detail here, but I
generate a two level list (.csv file), with a separate application,
containing thermodynamic properties specific internal energy, density,
and sonic velocity as a function of pressure and (system)
temperature. The list is generated using a constant value of entropy
and ... well, a fairly detailed description is in the attached
notebooks, as well as the table the notebook Imports is attached as
well. These values are for real, not ideal, gas behavior, which is
critically important where the system temperature (the temperature to
which the entire system is soaked prior to initiation) is quite low.
But this sample analysis is run at 70 degrees F.
The fact the flow out of the bottle is initially velocity choked
throws quite a twist in the analysis, but absolutely must be
captured. In this first analysis, in fact, the flow is velocity
choked for almost the entire process.
So I wrote conservation of mass and conservation of power equations,
and "folded" conservation of mass into the conservation of power
equation, yielding one ordinary first order differential equation with
bottle pressure as the time dependent variable, onto which I applied
NDSolve.
Now, the potential for the flow to be velocity choked necessitates
evaluations for the velocity which can be captured either as piecewise
functions, or evaluated If constructs. This is because where the flow
is velocity choked the pressure at the outlet is not the external
atmospheric pressure, but is rather a certain percentage of the
current bottle pressure (consult most any thermodynamics and/or fluid
mechanics text). As the bottle pressure ("stagnation" pressure)
bleeds down, the sonic velocity of the gas exiting the bottle also
changes, further changing the sonic velocity, the exit pressure, the
exit density, the exit specific internal energy, etc.
Physically, the bottle pressure cannot decrease below the external
(atmospheric) pressure. That is important if you read on.
Finally, we arrive at my dilemma: When I use evaluated If constructs
( Evaluate //@ If[ ... ) to branch to the proper exit pressure, exit
density, exit specific internal energy, as well as to branch to either
the sonic velocity interpolation (on pressure) or conservation of mass
to determine the exit velocity, NDSolve solves the problem without
incident (except for a minor one on converging to a stopping
criterion, and I have another story for that I won't get into here).
However, when I use mathematically equivalent piecewise constructs
replacing all the evaluated If constructs, NDSolve goes out of bounds
in the interpolations for specific internal energy, density, and sonic
velocity. NDSolve in this case below the physically possible lowest
bottle pressure, and, again, with the evaluated If constructs NDSolve
never goes "out of bounds".
That is very disturbing to me and I would like to know more about why
that is occurring. Apparently, even though the If constructs are
mathematically equivalent to the piecewise constructs, they are not
Mathematica 'lly equivalent, at least within NDSolve. That being the
case, a dark cloud descended upon me in Mathematica's voracity as a
whole.
If anyone could take a look I would be extremely appreciative to know
what is going on, at least inside NDSolve, where piecewise and
equivalent If constructs are not equivalent in Mathematica's NDSolve.
I see now I don't have an option to attach the files, so I will look
for a way to do that. If I cannot and you are gracious enough to look
into this and enlighten me, I receive the group via email, so I could
send them to you.
Thanks Very Much for your Patience!
Brian L.