(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)E − R X = 0 (Subsonic)
(P X + Q)E/(E−1) − R X = 0 (Supersonic)
These equations can be rewritten
(P + Q X)1/Q − R X = 0 (Subsonic)
(P X + Q)1/P − R 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/Q − R X = 0 (9)
Taking the derivative
f'(X) = (P + Q X)1/Q − 1 − R (10)
Also (though not needed for Newton-Raphson),
f''(X) = P (P + Q X)1/Q − 2 (11)
Newton-Raphson yields
Xnew = X − f(X) / f'(X)
= X − [ (P + Q X)1/Q − R X ] /
[ (P + Q X)1/Q − 1 − R ]
= [ X (P + Q X)1/Q − 1 − R X − (P + Q X)1/Q + R X ] /
[ (P + Q X)1/Q − 1 − R ]
= [ X (P + Q X)1/Q − 1 − (P + Q X)1/Q ] /
[ (P + Q X)1/Q − 1 − R ]
= [ (X − P − Q X) (P + Q X)1/Q − 1 ] /
[ (P + Q X)1/Q − 1 − R ]
= (X − P − Q 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 = Xnew − X
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 − 1 − R = 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