Gas Dynamics and exponents

[cspeed]

In an attempt to learn how to code in LV, I'm attempting to program the Mach tables into, and I'm having a little trouble.  One of the problems that I have run into is raising a number (other than 10, 2 and e) to a power (uaually a decimal number).

I could do this with a formula node, but I have to use this same equation in a numeric solver that I've programmed up, so I don't want to slow things down.  What I've done instead is use basic algebra to turn the expression into exp and ln terms.  Specifically, I've turned x^n into exp(n*ln(x)) which works just fine.  I'm wondering if I've really gained anything though, since transendental functions tend to be approximated numerically which can be kind of expensive.  Does anyone out there know which is more expensive (a formula node vs, exp&ln functions)?

Now let's talk about the numeric solver.  First, I tried using the Newton-Rahpson method, but it didn't converge quickly when it did converge which was rare (in fact, none of the open methods that I tried worked).  When I say that it didn't converge quickly I mean that after 1E6 iterations, I still wasn't within 4 significant figures of the accepted value(s).

Since open methods weren't working, I decided to try bracketing,  The False Position method was the first one I tried.  This is the one where the program basis the next guess on the intersection between a line connecting the guesses and the function evaluated at the guesses.  Between the two intervals, the program looks for a sign change and then it assumes that there is a zero within that interval.

The False Position method didn't work either.  It actually gave me worse results than the Newton-Raphson.  The next one I tried was the simplest bisection method around.  This one basis the next two guesses on a guess halfway between the first two guess, and looks for a sign change within the two intervals.  This one actually converged (with double precision) after about 50 iterations.  This was much better that 1E6, but not very good if I have to do it 5000000 times.  Since this subVI will eventually be used to find where a normal shock is, this is a very likely occurance.

I was wondering is anyone had programmed up a numeric solver that didn't require a string as an input.  I've got two variables, one which is given (Cp0/Cv0) and of course A/A*.  So the built-in solvers in LV won't do.  I'm also using the base package, which won't parse anything.

Anyone have any wisdom on this subject?

Thanks,

--CS

[altenbach]

---

@cspeed wrote:

In an attempt to learn how to code in LV, I'm attempting to program the Mach tables into, and I'm having a little trouble.  One of the problems that I have run into is raising a number (other than 10, 2 and e) to a power (uaually a decimal number).

---

LabVIEW has x^y in the "numeric...logarithmic palette". It is called "power of x".

Newton Raphson should converge pretty quickly unless the function is somehow pathological. What are you trying to solve?

---

@cspeed wrote:

I was wondering is anyone had programmed up a numeric solver that didn't require a string as an input.  I've got two variables, one which is given (Cp0/Cv0) and of course A/A*. 

---

Could you provide more details?

[cspeed]

"LabVIEW has x^y in the "numeric...logarithmic palette". It is called "power of x"."

I haven't had a chance to confirm this yet, but I'll feel pretty silly if it's true 🙂

"Newton Raphson should converge pretty quickly unless the function is somehow pathological. What are you trying to solve?"

It's a little difficult to read/write in text, but here goes:

A/A*=(1/M)[(2/(g+1))(1+((g-1)/2)M^2]^(g+1)/(2(g-1))

Where g is gamma, the ratio between the constant pressure specific heat and the constant volume specific heat.  Gamma depends on the fluid.  For air, gamma is around 1.4, for hydrogen, gamma is 4/3.  M is the Mach Number (which is always positive).

The equation itself is extremely useful for modeling compressible flows because all it requires is A* (the area that would cause the Mach Number to be one).  This is nice because once A* is found, one can find the Mach Number (and all other fluid properties) purely by knowing the geometry of the system.  It's also useful for finding where a normal shock wave is, a process that can only be done iteratively.

Thanks

--CS

[altenbach]

CS,

OK, your formula has some mismatched parenthesis and other unclear points. I have implemented it as follows:

A/A*=(1/M) [(2/(g+1)) (1+((g-1)/2)) M^2]^[(g+1)/(2(g-1))]

Is this correct? It could also be:

A/A*=(1/M) [(2/(g+1)) (1+((g-1)/2) M^2)]^[(g+1)/(2(g-1))]

Could you clarify which one is right? Is it also correct to assume the last two blue ones or is the exponent only (g+1)? Is there a website describing this formula?

A very simple Newton-Raphson implementation with numerical derivatives finds the correct M for any A/A* within less than 50 iterations, typically less than 10. (I tried between A/A*=0.0001 to 100000). with a tolerance of 1e-8).

(It works equally well with the second formula with the limitation that it produces two possible solutions, depending on input)

Attached is a simple Newton-Raphson implementation to solve M as a function of A/A* (with a given g) using formula 1 above. (LabVIEW 7.1. Note: There could be bugs please verify formula. There should also be a limit on the total iteration if a solution cannot be found).

[cspeed]

---

@altenbach wrote:

CS,

OK, your formula has some mismatched parenthesis and other unclear points. I have implemented it as follows:

A/A*=(1/M) [(2/(g+1)) (1+((g-1)/2)) M^2]^[(g+1)/(2(g-1))]

Is this correct? It could also be:

A/A*=(1/M) [(2/(g+1)) (1+((g-1)/2) M^2)]^[(g+1)/(2(g-1))]

Could you clarify which one is right? Is it also correct to assume the last two blue ones or is the exponent only (g+1)? Is there a website describing this formula?

---

The second one is the correct one, sorry about the abiguity.  This is touched on briefly at a http://me.lsu.edu/~meniki/me4621/manual/Lab%20Manual_v4.pdf on page 23.  Sorry about the .pdf, I had a hard time finding the equations.

---

A very simple Newton-Raphson implementation with numerical derivatives finds the correct M for any A/A* within less than 50 iterations, typically less than 10. (I tried between A/A*=0.0001 to 100000). with a tolerance of 1e-8).

(It works equally well with the second formula with the limitation that it produces two possible solutions, depending on input)

---

This sounds about right, one is the subsonic solution (M<1) the other is the supersonic solution (M>1,of course).  The subsonic solution ideally happens on the way into a converginig-diverging nozzle, while the supersonic solution happens on the way out.  Either way, there sould be two answers, both positive.  One number is between 0 and 1, the other is bigger than 1.  The only time this isn't true is when A/A*=1 (since A* represents the cross-sectional area needed to accelerate the flow to Mach 1).

---

Attached is a simple Newton-Raphson implementation to solve M as a function of A/A* (with a given g) using formula 1 above. (LabVIEW 7.1. Note: There could be bugs please verify formula. There should also be a limit on the total iteration if a solution cannot be found).

---

Thank you very much.  I will use it to figure out where I have gone wrong.  I often make a sign error when I am using the Newton- Raphson method.  I always try to add instead of subtract.

Thanks again.

--CS