Clarification of Answer by
mathtalk-ga
on
25 Jan 2003 20:55 PST
/*/ obliqshk.c (mathtalk-ga : 2003-01-25)
These routines are related to formation of attached shock waves
formed by steady supersonic flow past a wedge-shaped obstacle.
For the isentropic flow of an ideal gas, we have the equation:
cot(A) = tan(S) [ (g + 1)M^2/(M^2 sin^2(S) - 1) - 1 ]
where:
A = angle of symmetric wedge surface wrt downstream
S = angle of symmetric shock boundry wrt downstream
M = Mach number ahead of shock wave (assumed > 1)
g = gamma, ratio of specific heats in ideal gas law
This equation gives A more or less explicitly as a function of
M, g, and S. Our general goal is to provide a means to solve
for S given M, g, and A.
The following should be considered:
1) Angle S is naturally constrained to between:
arcsin(1/M) < S < pi/2
2) On this range the right hand tends to +oo at both endpoints,
with one interior minimum. Angle A is thus constrained:
0 < A <= A_max
with A_max being the inverse cotangent of the RHS minimum
depending only on M and g.
3) For particular A and given M,g there may be no solution S
(if A > A_max), one solution S (if A = A_max), or two (if
A < A_max). When two solutions are available, the smaller
one is the stable solution and ordinarily the physically
realized one. The smaller angle corresponds to what is
called a "weak" shock wave, the larger angle to what is
called a "strong" shock wave. The strong shock wave has
a large increase in pressure and a correspondingly large
drop in Mach speed to subsonic flow. The weak shock wave
has a smaller increase in pressure and in most cases the
flow beyond the shock remains supersonic.
Taking these observations into account, the design of these
routines aims to provide for the smaller angle solution and
for determination of A_max, given M and g, so that existence
of a solution can be determined beforehand.
/*/
#include <math.h>
#include <assert.h>
#include <stdio.h>
/*/
The expressions involving angle S can be reduced to a
rational form involving u = cot(S), leading actually to
a cubic equation to solve for u in terms of M, g, and
tan(A). However those routines below which directly
take u as an argument are for "internal use" and are
thus marked with the static keyword.
tan(A) = u (R - u^2)/(T - R + (T+1)u^2)
where:
R = M^2 - 1
T = M^2 (g + 1)
Since this is C code, we prefer lowercase for variables.
Hence Mach number becomes m instead of M, and the angles
A and S are represented in code with a and s, resp.
One of the cheesier aspects of the code below is the
free overwriting of calling arugments, passed by value,
on the way to quick computation of the return result.
/*/
static double amgu(double m, double g, double u)
{
g = (g + 1.0) * (m *= m); // overwrites g with T
m -= 1.0; // finish overwriting m with R
return atan(u*(m - u*u)/(g - m + (g+1.0)*u*u));
}
double wedgeangle(double m, double g, double s)
// compute angle A from M, g, and S
{
return amgu(m,g,1.0/tan(s)); // pass in u = cot(S)
}
/*/
We note above that tan(A) can be expressed in terms
of a rational function of u:
f(u) = u (R - u^2)/(T + (T+1)u^2)
which is reciprocal to the expression for cot(A) once
u is identified with cot(S). Now f(u) is bounded on
[0,sqrt(R)] and odd as a function of u.
So the function f(u) attains a maximum on this range
where f'(u) vanishes (note that f(u) is zero at both
endpoints and otherwise positive in the given interval).
But now f'(u) is an even function of u, and hence can
be expressed in terms of w = u^2. In fact f'(u) = 0
if and only if:
(T+1)w^2 + (R(T+1) + 3(T-R))w - R(T-R) = 0
a quadratic in w with exactly one positive root.
/*/
static double maxu(double m, double g)
// find u that maximizes expression f(u)
{
double w;
g = (g + 1.0) * (m *= m); //overwrite g with R
m -= 1.0; // finish overwriting m with T
w = m*(g + 1.0) + 3.0*(g - m);
w += sqrt( w*w + 4.0*m*(g - m)*(g + 1.0));
return sqrt(2.0*m*(g - m)/w);
}
double wedgemax(double m, double g)
// compute A_max from M & g
{
double u;
u = maxu(m,g);
return amgu(m,g,u);
}
/*/
Finally we turn to solving for u = cot(S) given the
values of M, g, and A. Letting v = cot(A) we see
that the rational function f(u) = v gives a cubic:
u^3 + v(T+1)u^2 - Ru + v(T-R) = 0
which, by Descarte's rule of signs, must have two
positive real roots or none. We assume v > 0 and:
0 < M^2 - 1 = R < M^2 <= T = (g+1)M^2
The interesting point here is that in order to get
the smaller angle S, we want to find the root u of
this equation between maxu(m,g) and sqrt(R). This
is because cot(S) is a decreasing function, and so
the larger root in terms of u gives the smaller
solution for angle S.
A number of root finding strategies or tactics are
available to us. At one extreme we might apply a
variant of Cardano's formula; at the other one may
try exhaustive search.
I recently read some ideas on explicit solutions
to cubic equations that I'd like to try, but for
now Newton's method seems most attractive. Note
that the cubic polynomial is concave up when:
0 <= u <= sqrt(R)
and takes a positive value at u = sqrt(R) since:
u^3 - Ru = 0
there. So u = sqrt(R) makes a suitable starting
value, provided a positive root lies below it,
i.e provided the cubic polynomial evaluated at
maxu(m,g) is negative.
To minimize dependencies I decided not to call
maxu(m,g) to directly test this feasibility. My
alternative relies on the polynomial's convexity
on [0,sqrt(R)]. Using u = sqrt(R) as an initial
guess, we'd approach the larger of any positive
root montonically from above. Consequently the
deriviatives would remain positive in this case.
However if no positive real roots exist, we will
eventually slip around the local minimum in such
a manner that either u or the slope becomes non-
positive. Asserts are placed to guard against
these possibilities.
/*/
static double umga(double m, double g, double a)
// find root for u using tailored Newton's method
{
double u,v,p,q;
g = 1.0 + g*(m *= m); // overwrite g with T-R
v = tan(a); // overwrite a with v = tan(A)
// p = u^3 + v(T+1)u^2 - Ru + v(T-R)
// q = dp/du = 3u^2 + 2v(T+1)u - R
a = v*(g + m); // overwrite a with v(T+1)
m -= 1.0; // finish overwriting m with R
g *= v; // and overwrite g with v(T-R)
u = sqrt(m);
do {
p = ((u + a)*u - m)*u + g;
q = (3.0*u + a + a)*u - m;
assert(q > 0.0);
u -= p/q;
assert(u > 0.0);
} while (p > 1.0e-6); // Newton iteration
return u;
}
double shockangle(double m, double g, double a)
// compute angle S from M, g, and A
{
return atan(1.0/umga(m,g,a)); // u = cot(S)
}
/*/
test data (angles in degrees):
M A S gamma
2.5 15.0 36.94490 0.2
/*/
int main(int argc, char* argv[])
{
double a,s;
const double pi = 3.14159265358979323846;
a = wedgemax(2.5,0.2)*180.0/pi;
printf("Max. wedge if Mach 2.5, gamma 0.2 is %f degrees.\n", a);
if (a > 15.0)
{
s = shockangle(2.5,0.2,pi/12.0); // pi/12 is 15 degrees
printf("Shock angle for 15 degree wedge is %f degrees.\n",
s*180.0/pi);
// s = (36.9449/180.0)*pi;
a = wedgeangle(2.5,0.2,s)*180.0/pi;
printf("Corresponding wedge angle is %f degrees.\n", a);
}
else printf("Maximum wedge angle is below 15 degrees.\n");
return 0;
}
,
0 Comments
)