Skip navigation links
(Wind-US Documentation Home Page) (Wind-US User's Guide) (GMAN User's Guide) (MADCAP User's Guide) (CFPOST User's Guide) (Wind-US Utilities) (Common File User's Guide) (Wind-US Installation Guide) (Wind-US Developer's Reference) (Guidelines Documents)

(b4wind User's Guide) (Converting Files) (Computing Initial Conditions) (Interpolation) (Computing the Reynolds Number) (Preparing NPARC or XAIR Namelists) (Installation and Execution) (Conclusions) (Computing the Mach Number From the Area Ratio) (Trilinear Interpolation)

Computing the Mach Number From the Area Ratio

Given the equation

At / A = M { [ (γ+1) / 2 ] / [ 1 + ((γ−1) / 2) M2 ] }E/2   (1)

where A is the cross-sectional area, At is the cross-sectional area of the throat, γ is the ratio of specific heats, M is the Mach number, and E = (γ+1)/(γ−1).

Assume everything is known but M. This equation has two solutions. At a point upstream to the throat the flow is subsonic, so there is the subsonic solution giving a M < 1. Downstream of the throat the flow is supersonic giving a solution, M > 1. Karl Kneile devised an efficient method for computing M.

Initial Guess

The equation will be solved using an iterative process which which requires a starting solution. A very good first approximation can be obtained as follows.

The Subsonic Case

The subsonic case will be done first.

[(A / At) M ]2 = { [ 1 + ((γ−1) / 2) M2 ] / [ (γ+1)/2) ] }E

[(A / At) M ]2 = { [ (2 / (γ−1)) + M2 ] / [ (γ+1) / (γ−1) ] }E

[(A / At) M ]2 = { [ (2 / (γ−1)) + M2 ] / E }E

Or,

R X = (P + Q X)E   (2)

where

R = (A / At)2   (3)

X = M2

P = 2 / (γ + 1) = (E − 1) / E

Q = (γ − 1) / (γ + 1) = 1 / E

Note that P + Q = 1 and EQ = 1.

Let

F(X) = (P + Q X)E

Then the derivative is

F'(X) = E Q (P + Q X)E−1 = (P + Q X)E−1

So

F(0) = PE

F(1) = (P + Q)E = 1

F'(1) = (P + Q)E−1 = 1

If F is approximated by a quadratic, then to the accuracy of the quadratic

F(X) = a + b X + c X2

F'(X) = b + 2 c X

Thus

F(0) = a = PE

F(1) = a + b + c = 1

F'(1) = b + 2 c = 1

which has the solution

a = PE = P1/Q

b = 1 − 2a

c = a

The Supersonic Case

Now doing the same thing with the supersonic case let

Y = 1 / X

Equation (2), at the start of the subsonic case, can now be written

R / Y = (P + Q/Y)E   (4)

or,

R YE−1 = (P Y + Q)E

Rewriting and changing the definitions of the symbols,

R X = (P X + Q)E/(E−1)

where

Rnew = Rold1/(E−1) = (A / At)2/(E−1)   (5)

Xnew = 1 / Xold

Following the same process as before (and using the same symbols)

F(X) = (P X + Q)E/(E−1)

F'(X) = [E/(E−1)] P (P X + Q)1/(E−1) = (P X + Q)1/(E−1)

Thus

F(0) = a = QE/(E−1)

F(1) = a + b + c = 1

F'(1) = b + 2 c = 1

which yields

a = QE/(E−1) = Q1/P

b = 1 − 2a

c = a

Equations (2) and (4) become

R X = a + b X + c X2

= a + (1 − 2a) X + a X2

or

X2 + [ (1 − R − 2a) / a ] X + 1 = 0

X2 − 2 [ (R − 1) / (2a) + 1 ] X + 1 = 0

or

X2 − 2 (r + 1) X + 1 = 0

where

r = (R − 1) / (2a)

Note that in both cases (subsonic and supersonic) from Eqs (3) and (5), that since At (the cross-sectional area at the throat) is less than A, R is greater than one; thus r is positive. The quadratic formula yields

X = (1 + r) ± sqrt [ (1 + r)2 − 1 ]

Since in the subsonic case X = M2 and for the supersonic X = 1/M2, in both cases it is desired that X be less than one; thus the minus sign is chosen, giving

X = (1 + r) − sqrt [ (1 + r)2 − 1 ]

From this form it seen that X > 0. Numerically, a better form is

X = 1 / { (1 + r) + sqrt [ r ( r + 2 ) ] }

From this form (r is positive) it is seen that X < 1.

Summary

Summarizing,

P = 2 / (γ + 1)

Q = 1 − P

E = 1 / Q

R = (A / At)2  (Subsonic)   (6a)

R = (A / At)2Q/P  (Supersonic)   (6b)

a = P1/Q  (Subsonic)

a = Q1/P  (Supersonic)

r = (R − 1) / (2a)

X = 1 / { (1 + r) + sqrt [ r ( r + 2 ) ] }   (7)

M = sqrt (X)  (Subsonic)   (8a)

M = 1 / sqrt (X)  (Supersonic)   (8a)

This gives a good first approximation for an iterative solution to

(P + Q X)ER X = 0  (Subsonic)

(P X + Q)E/(E−1)R X = 0  (Supersonic)

These equations can be rewritten

(P + Q X)1/QR X = 0  (Subsonic)

(P X + Q)1/PR X = 0  (Supersonic)

Notice that the above two equations are of identical form except that P and Q are interchanged. Therefore it is sufficient to the first of these two.

Newton-Raphson

Looking first at Newton-Raphson, from Eq. (2) let

f(X) = (P + Q X)1/QR X = 0   (9)

Taking the derivative

f'(X) = (P + Q X)1/Q − 1R   (10)

Also (though not needed for Newton-Raphson),

f''(X) = P (P + Q X)1/Q − 2   (11)

Newton-Raphson yields

Xnew = Xf(X) / f'(X)

= X − [ (P + Q X)1/QR X ] / [ (P + Q X)1/Q − 1R ]

= [ X (P + Q X)1/Q − 1R X − (P + Q X)1/Q + R X ] / [ (P + Q X)1/Q − 1R ]

= [ X (P + Q X)1/Q − 1 − (P + Q X)1/Q ] / [ (P + Q X)1/Q − 1R ]

= [ (XPQ X) (P + Q X)1/Q − 1 ] / [ (P + Q X)1/Q − 1R ]

= (XPQ X) / [ 1 − R (P + Q X)(1 − 1/Q) ]

Recall that

1 − Q = P

and

1 − 1 / Q = (Q − 1) / Q = −P / Q

so

Xnew = P (X − 1) / [ 1 − R (P + Q X)P/Q ]

which is the result for Newton-Raphson.

A Higher Order Method

Better, let

D = XnewX

Then

f(Xnew) = f(X) + f'(X) D + f''(X) D2 / 2 = 0

or

f'' D2 + 2f' D + 2f = 0

Solving for D

D = [−f' + s sqrt (f'2 − 2f f'') ] / f''

where s is ±1, yet to be determined. Multiplying numerator and denominator by [ −f's sqrt (f'2 − 2f f'') ] and simplfying,

D = −2f / [ f' + s sqrt (f'2 − 2f f'') ]   (12)

Making some observations, from Eqs. (6), since At (the throat cross-sectional area) is less than any other cross-sectional area, R > 1. From Eq. (9), f(0) = P(1/Q) > 0 and F(1) = 1 − R < 0. Therefore there is a root between 0 and 1. From Eq. (10), since X < 1, f' < (P + Q)1/Q − 1R = 1 − R < 0. So f' is always negative and by Eq. (11) it is seen that f'' is always positive.

In the neighborhood of a root, it is necessary that D approach zero as f approaches zero in order for the process to converge. In Eq. (12), choosing s to be +1 and letting f approach zero, since sqrt(f'2) = abs(f'), which is positive, and since f' is negative, it is seen that there is a zero over zero. Therefore the limit must be obtained by L'Hospital's rule. Taking the derivative of the numerator and the denominator with respect to f results in

sqrt(f'2 − 2f f'') / f''

which in the limit as f approaches zero is abs(f') / f'' which is not zero. Therefore s must be chosen as −1. Choosing s as −1

D = −2f / [ f' − sqrt (f'2 − 2f f'') ]

Taking the limit as f approaches zero gives

D = −2f / [ f' − abs (f') ]

and since f' is negative this gives

D = −f / f'

which is Newton-Raphson and is seen to be zero in the limit as required.

Therefore the root is found by iterating

Xnew = X − 2f / [ f' − sqrt (f'2 − 2f f'') ]

and the Mach number is given by Eqs. (8). Obtaining the first approximation by Eq. (7), then iterating the above equation, and then comparing the results with a gas table, it was found that only two iterations were necessary to obtain the accuracy of the table.


Last updated 5 May 2003