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.