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

() (Converting Files) (Computing Initial Conditions) (Interpolation) (Computing the Reynolds Number) (Preparing NPARC or XAIR Namelists) (Installation and Execution) (Conclusions) () (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