Category When Is A Flow Compressible?

Elements of Finite-Difference Methods

The method of characteristics described in the previous section legitimately can be considered a part of computational fluid dynamics because it uses discrete algebraic forms of the governing equations (such as Equations (13.17) and (13.18)) which are solved at discrete points in the flow (the characteristic mesh illustrated in Figure 13.5). However, most authors consider that CFD is represented by mainly finite difference and finite volume techniques, such as are discussed in Reference 64, and the method of characteristics is usually not included in the study of CFD. The purpose of this section is to give you the flavor of finite-difference techniques by describing one particular

method that is readily applicable to a number of compressible flow problems. The method discussed here is representative of mainstream CFD, but it is just the tip of the iceberg of CFD. The intensive research in CFD since 1960 has produced a multitude of different algorithms and philosophies, and it is far beyond the scope of this book to go into the details of such work. See Reference 64 for an in-depth presentation of CFD at the introductory level. In addition, you are strongly encouraged to read the current literature in this regard, in particular the AIAA Journal, Computers and Fluids, and the Journal of Computational Physics. Finally, in this chapter we are dealing with numerical solutions of inviscid supersonic flows. See Reference 21 for an expanded discussion of finite difference methods applied to supersonic flows.

How do we use the finite differences obtained here? Imagine that a flow in xy space is covered by the mesh shown in Figure 13.2b. Assume there are N grid points. At each one of these grid points, evaluate the continuity, momentum, and energy equations with their partial derivatives replaced by the finite-difference expressions derived above. For example, replacing the derivatives in Equations (7.40), (7.42a and b), and (7.44) with finite differences, along with Equations (7.1) and (7.6a), we obtain a system (over all N grid points) of 6N simultaneous nonlinear algebraic equations in terms of the 6N unknowns, namely, p, u, v, p, T, and e, at each of the N grid points. In principle, we could solve this system for the unknown flow variables at all the grid points. In practice, this is easier said than done. There are severe problems in solving such a large number of simultaneous nonlinear equations. Moreover, we have to deal with problems associated with numerical instabilities that sometimes cause such attempted solutions to “blow up’’ on the computer. Finally, and most importantly, we must properly account for the boundary conditions. These considerations make

all finite-difference solutions a nontrivial endeavor. As a result, a number of spe­cialized finite-difference techniques have evolved, directed at solving different types of flow problems and attempting to increase computational efficiency and accuracy. It is beyond the scope of this book to describe these difference techniques in detail. However, one technique in particular was widely used during the 1970s and 1980s. This is an approach developed in 1969 by Robert MacCormack at the NASA Ames Research Center. Because of its widespread use and acceptance at the time, as well as its relative simplicity, we will describe MacCormack’s technique in enough detail to give you a reasonable understanding of the method. This description will be carried out in the context of the following example.

Consider the two-dimensional supersonic how through the divergent duct shown in Figure 13.8a. Assume the how is supersonic at the inlet, and that all properties are known at the inlet. That is, the how-field variables at grid points (1, 1), (1, 2), (1, 3),

(1,4) , and (1,5) are known. The duct is formed by a fiat surface at the bottom and a specified contour, ys = f(x), at the top. In addition, assume that the how is inviscid, adiabatic, and steady, and with no body forces. It can be rotational or irrotational—the method of solution is the same. The governing equations are obtained from Equations (7.40), (7.42a and b), (7.44), (7.1), and (7.6a), which yield

[13.35]

Let us express these equations in slightly different form, as follows. Multiplying Equation (13.30) by u, and adding the results to Equation (13.31), we have

Similarly, multiplying Equation (13.30) by v, and adding the result to Equation

(13.32) , we obtain

d(puv) d (pv2 + p)

dx dy

Then, Equations (13.30) and (13.36) to (13.38) become

Equations (13.40) to (13.43) are the continuity, * and у momentum, and energy equations, respectively—but in a slightly different form from those we are used to seeing. The above form of these equations is frequently called the conservation form. Let us now treat F, G, H, and К as our primary dependent variables; these quantities are called flux variables, in contrast to the usual p, p, T, u, v, e, etc., which are called primitive variables. It is important to note that once the values of F, G, H, and К are known at a given grid point, the primitive variables at that point can be found from Equations (13.39a to d) and

p = pRT [13.44]

e = cvT [13.45]

V2 = u2 + v2 [13.46]

That is, Equations (13.39a to d) and (13.44) to (13.46) constitute seven algebraic

equations for the seven primitive variables, p, u, v, p, e, T, and V.

Let us return to the physical problem given in Figure 13.8a. Because the duct diverges, it is difficult to deal with an orthogonal, rectangular mesh; rather, a mesh which conforms to the boundary of the system will be curved, as shown in Figure 13.8a. On the other hand, to use our finite-difference quotients as given in Equation (13.27), (13.28), or (13.29), we desire a rectangular computational mesh. Therefore, we must transform the curved mesh shown in Figure 13.8a, known as the physical plane, to a rectangular mesh shown in Figure 13.8b, known as the computational plane. This transformation can be carried out as follows. Define

|=jc [13.47a]

У

h = —

ys

where ys = f(x) [13.47b]

In the above transformation, t] ranges from 0 at the bottom wall to 1.0 at the top wall. In the computational plane (Figure 13.8b), t] = constant is a straight horizontal line, whereas in the physical plane, t] = constant corresponds to the curved line shown in Figure 13.8. Because we wish to apply our finite differences in the computational plane, we need the governing equations in terms of f and r rather than x and у. To ac­complish this transformation, apply the chain rule of differentiation, using Equations

[13.48]

[13.49]

.43) become

Note in the above equations that the £ derivatives are on the left and the i] derivatives are all grouped on the right.

Let us now concentrate on obtaining a numerical, finite-difference solution of the problem shown in Figure 13.8. We will deal exclusively with the computational plane, Figure 13.8b, where the governing continuity, x and у momentum, and energy equations are given by Equations (13.50) to (13.53), respectively. Grid points (1, 1), (2, 1), (1, 2), (2, 2), etc., in the computational plane are the same as grid points (1, 1), (2, 1), (1, 2), (2, 2), etc., in the physical plane. All the flow variables are known at the inlet, including F, G, H, K. The solution for the flow variables downstream of the inlet can be found by using MacCormack’s method, which is based on Taylor’s series expansions for F, G, Я, and К as follows:

gm-j=gij+&Lm

/дН

Xw-j=Xij+(wL&t

In Equations (13.54a to d), F, G, Я, and К at point (i, j) are considered known, and these equations are used to find F, G, Я, and К at point (/ + 1, j) assuming that we can calculate the values of (9F/9£)ave, (9G/9£)ave, etc. The main thrust of MacCormack’s method is the calculation of these average derivatives. Examining

Equations (13.54a to d), we find that this finite-difference method is clearly a “down­stream marching” method; given the flow at point (i, j) we use Equations (13.54a to d) to find the flow at point (i + l, j). Then the process is repeated to find the flow at point (i + 2, j), etc. This downstream marching is similar to that performed with the method of characteristics.

The average derivatives in Equations (13.54a to d) are found by means of a straightforward “predictor-corrector” approach, outlined below. In carrying out this approach, we assume that the flow properties are known at grid point (i, j), as well as at all points directly above and below (г, j), namely, at (г, j + 1), (г, j +2), (г, j — 1), (ij -2), etc.

Two-Dimensional Poiseuille Flow

Consider the parallel flow between two horizontal plates separated by the distance D as sketched in Figure 16.12. In this case, unlike Couette flow where one of the plates is in motion, we consider both plates to be stationary. Recall that the driving force which established Couette flow was the shear stress between the moving plate and the fluid. In the present case sketched in Figure 16.12, what is the driving force; that is, what makes the fluid move? In Chapter 1, we emphasized that the only way that nature can exert a force on a fluid is by means of shear stress and pressure distributions. In the present problem, since the walls are not moving, there is no shear stress to drive the flow. Hence, the only other possibility is the pressure distribution. Indeed, to establish the flow shown in Figure 16.12, there must be a pressure gradient acting

on the gas. Moreover, in Figure 16.12 the flow extends to infinity in both directions along the x axis. As in the case of Couette flow, this implies that the velocity и is independent of x; that is, и = и (у). Since the streamlines are parallel, v = 0. This flow is called two-dimensional Poiseuille flow, named after the French physician J. L. M. Poiseuille, who studied similar flows in pipes.

Let us examine the Navier-Stokes equations in light of the problem just outlined. For simplicity, we will consider only steady, incompressible flow. First, return to the continuity equation for an incompressible flow, given by Equation (3.39). In cartesian coordinates, this is

du dv dx 9y

Since the flow is parallel, v = 0; hence, dv/dy = 0. From Equation (16.99), then du/dx = 0; this confirms that и is constant with x; that is, и is a function of у only. From the у-momentum equation, Equation (15.19b), we have

dp

— =0 [16.100]

9y

Hence, p varies only in the x direction; p = p(x). From the x-momentum equation, Equation (15.19a), we have

On the left-hand side of Equation (16.102), p is a function of x only. On the right-hand side of Equation (16.102), и is a function of у only. Hence, the left-hand and right- hand sides of Equation (16.102) must be equal to the same constant. This confirms an important aspect of this flow, namely, the pressure gradient is constant along the flow direction. Once again, we emphasize that it is this pressure gradient that drives the flow. The pressure gradient must be provided by an outside mechanism, that is, some source of high pressure toward the left and low pressure toward the right.

The velocity profile across the flow is obtained by solving Equation (16.102). For convenience, and to emphasize that p = p(x), we write the pressure derivative as an ordinary derivative, dp/dx. Integrating Equation (16.102) twice across the flow, we have

^ y2 + ay + b [16.103]

where a and b are constants of integration. Evaluating Equation (16.103) at у = 0, where и = 0, we have

b = 0

Evaluating Equation (16.103) at у = D, where и = 0, we have

Hence, Equation (16.103) becomes

At у = 0, Equation (16.106) yields

/Зи D /dp

Эу/ш 2/л dx

Hence, the wall shear stress is

Note the interesting fact from Equation (16.108) that zw does not depend on the viscosity coefficient (i, but rather only on the separation distance of the walls and on the pressure gradient. Clearly, this flow is a force balance between the pressure gradient acting toward the right on the gas and the shear stress at the walls acting toward the left on the gas.

This flow is sometimes called fully developed flow, for the following reason. Consider an actual flow in the laboratory wherein a uniform flow enters a channel, such as shown in the photograph in Figure 16.13. Here, velocity profiles in water flow are made visible by the hydrogen bubble method, where the bubbles are generated by electrolysis on a fine wire used as a cathode at the entrance of the channel. Near the entrance, the flow is uniform over a large portion of the distance across the channel; the viscous effects are limited to a thin boundary layer at the walls. However, as the flow progresses downstream, the viscous effects are felt over a larger portion of the flow. Finally, after the flow has covered a sufficient distance through the channel, the velocity profile is totally dominated by viscosity; a parabolic velocity profile is achieved, and the real flow becomes essentially the Poiseuille flow studied in this section. When this type of real flow is reached in the channel (at the right of Figure 16.13), it is called “fully developed flow.”

16.6 Summary

The parallel flows discussed in this chapter illustrate features common to many more complex viscous flows, with the added advantage of lending themselves to a relatively straightforward solution. The purpose of this discussion has been to introduce many of the basic concepts of viscous flows in a fashion unencumbered by fluid dynamic complexities. In particular, we have studied Couette and two-dimensional Poiseuille flows and found the following.

Supersonic Wind Tunnels

Return to the road map given in Figure 10.3. The material for the left and right branches is covered in Sections 10.3 and 10.4, respectively. In turn, a mating of these two branches gives birth to the fundamental aspects of supersonic wind tunnels, to be discussed in this section.

Imagine that you want to create a Mach 2.5 uniform flow in a laboratory for the purpose of testing a model of a supersonic vehicle, say, a cone. How do you do it? Clearly, we need a convergent-divergent nozzle with an area ratio Ae/A* = 2.637 (see Appendix A). Moreover, we need to establish a pressure ratio, pa/pe — 17.09, across the nozzle in order to obtain a shock-free expansion to Me = 2.5 at the exit. Your first thought might be to exhaust the nozzle directly into the laboratory, as sketched in Figure 10.16. Here, the Mach 2.5 flow passes into the surroundings as a “free jet.’’ The test model is placed in the flow downstream of the nozzle exit. In order to make certain that the free jet does not have shock or expansion waves, the nozzle exit pressure pe must equal the back pressure pB, as originally sketched in Figure 10.14c. Since the back pressure is simply that of the atmosphere surrounding the free jet, pa — pe = 1 atm. Consequently, to establish the proper isentropic expansion through the nozzle, you need a high-pressure reservoir with po — 17.09 atm at the inlet to the nozzle. In this manner, you would be able to accomplish your objective, namely, to produce a uniform stream of air at Mach 2.5 in order to test a supersonic model, as sketched in Figure 10.16.

Figure 10.l 6 Nozzle exhausting directly to the atmosphere.

In the above example, you may have a problem obtaining the high-pressure air supply at 17.09 atm. You need an air compressor or a bank of high-pressure air bottles—both of which can be expensive. It requires work, hence money, to create reservoirs of high-pressure air—the higher the pressure, the more the cost. So, can you accomplish your objective in a more efficient way, at less cost? The answer is yes, as follows. Instead of the free jet as sketched in Figure 10.16, imagine that you have a long constant-area section downstream of the nozzle exit, with a normal shock wave standing at the end of the constant-area section; this is shown in Figure 10.17. The pressure downstream of the normal shock wave is p2 = Рв = 1 atm. At M = 2.5, the static pressure ratio across the normal shock is рг/ре = 7.125. Hence, the pressure upstream of the normal shock is 0.14 atm. Since the flow is uniform in the constant-area section, this pressure is also equal to the nozzle exit pressure; that is, pe = 0.14 atm. Thus, in order to obtain the proper isentropic flow through the nozzle, which requires a pressure ratio of po/Pe — 17.09, we need a reservoir with a pressure of only 2.4 atm. This is considerably more efficient than the 17.09 atm required in Figure 10.16. Hence, we have created a uniform Mach 2.5 flow (in the constant-area duct) at a considerable reduction in cost compared with the scheme in Figure 10.16.

In Figure 10.17, the normal shock wave is acting as a diffuser, slowing the air originally at Mach 2.5 to the subsonic value of Mach 0.513 immediately behind the shock. Hence, by the addition of this “diffuser,” we can more efficiently produce our uniform Mach 2.5 flow. This illustrates one of the functions of a diffuser. However, the “normal shock diffuser” sketched in Figure 10.17 has several problems;

1. A normal shock is the strongest possible shock, hence creating the largest total pressure loss. If we could replace the normal shock in Figure 10.17 with a weaker shock, the total pressure loss would be less, and the required reservoir pressure po would be less than 2.4 atm.

2. It is extremely difficult to hold a normal shock wave stationary at the duct exit; in real life, flow unsteadiness and instabilities would cause the shock to move somewhere else and to fluctuate constantly in position. Thus, we could never be certain about the quality of the flow in the constant-area duct.

Nozzle exit Normal shock

Figure 1 0*1 7 Nozzle exhausting into a constant-area duct, where a normal shock stands at the exit of the duct.

3. As soon as a test model is introduced into the constant-area section, the oblique

waves from the model would propagate downstream, causing the flow to become

two – or three-dimensional. The normal shock sketched in Figure 10.17 could not

exist in such a flow.

Hence, let us replace the normal shock in Figure 10.17 with the oblique shock diffuser shown in Figure 10.15/?. The resulting duct would appear as sketched in Figure 10.18. Examine this figure closely. We have a convergent-divergent nozzle feeding a uniform supersonic flow into the constant-area duct, which is called the test section. This flow is subsequently slowed to a low subsonic speed by means of a diffuser. This arrangement—namely, a convergent-divergent nozzle, a test section, and a convergent-divergent diffuser—is a supersonic wind tunnel. A test model, the cone in Figure 10.18, is placed in the test section, where aerodynamic measurements such as lift, drag, and pressure distribution are made. The wave system from the model propagates downstream and interacts with the multireflected shocks in the diffuser. The pressure ratio required to run the supersonic tunnel is ро/рв – This can be obtained by making po large via a high-pressure reservoir at the inlet to the nozzle or by making pB small via a vacuum source at the exit of the diffuser, or a combination of both.

The main source of total pressure loss in a supersonic wind tunnel is the diffuser. How does the oblique shock diffuser in Figure 10.18 compare with the hypothetical normal shock diffuser in Figure 10.17? Is the total pressure loss across all the reflected oblique shocks in Figure 10.18 greater or less than across the single normal shock wave in Figure 10.17? This is an important question, since the smaller the total pressure loss in the diffuser, the smaller is the pressure ratio ро/рв required to run the supersonic tunnel. There is no pat answer to this question. However, it is usually true that progressively reducing the velocity of a supersonic flow through a series of oblique shocks to a low supersonic value, and then further reducing the flow to subsonic speeds across a weak normal shock, results in a smaller total pressure loss than simply reducing the flow to subsonic speeds across a single, strong normal shock wave at the initially high supersonic Mach number. This trend is illustrated by Example 9.4. Therefore, the oblique shock diffuser shown in Figures 10.15/? and 10.18 is usually more efficient than the simple normal shock diffuser shown in Figure 10.17. This is not always true, however, because in an actual real-life oblique shock diffuser, the

shock waves interact with the boundary layers on the walls, causing local thickening and even possible separation of the boundary layers. This creates an additional total pressure loss. Moreover, the simple aspect of skin friction exerted on the surface generates a total pressure loss. Hence, actual oblique shock diffusers may have efficiencies greater or less than a hypothetical normal shock diffuser. Nevertheless, virtually all supersonic wind tunnels use oblique shock diffusers qualitatively similar to that shown in Figure 10.18.

Notice that the supersonic wind tunnel shown in Figure 10.18 has two throats: the nozzle throat with area At i is called the first throat, and the diffuser throat with area A, 2 is called the second throat. The mass flow through the nozzle can be expressed as m — pu A evaluated at the first throat. This station is denoted as station 1 in Figure 10.18, and hence the mass flow through the nozzle is mi — pxuAt = p*a*Att j. In turn, the mass flow through the diffuser can be expressed as m = pu A evaluated at station 2, namely, m2 = p2u2A,2- For steady flow through the wind tunnel, m = rhj. Hence,

раА, л = р2и2А, л

Since the thermodynamic state of the gas is irreversibly changed in going through the shock waves created by the test model and generated in the diffuser, clearly p2 and possibly u2 are different from p* and a*, respectively. Hence, from Equation (10.33), the second throat must have a different area from the first throat; that is, Лі2 ф A, |.

Question: How does A, 2 differ from A, ? Let us assume that sonic flow occurs at both stations 1 and 2 in Figure 10.18. Thus, Equation (10.33) can be written as

A,,2 = /Offlf Af, 1 P2a2

Recall from Section 8.4 that a* is constant for an adiabatic flow. Also, recall that the flow across shock waves is adiabatic (but not isentropic). Hence, the flow throughout the wind tunnel sketched in Figure 10.18 is adiabatic, and therefore a* = a*2. In turn, Equation (10.34) becomes

At,2 _

At, і P2

Recall from Section 8.4 that T* is also constant throughout the adiabatic flow of a calorically perfect gas. Hence, from the equation of state,

РІ = p*/RT* = p p2 p2/RTf P2

Substituting Equation (10.36) into (10.35), we have

At,2 _ p

Аг, 1 P2

From Equation (8.45), we have

Substituting the above into Equation (10.37), we obtain

Examining Figure 10.18, the total pressure always decreases across shock waves; therefore, po,2 < Po. i – In turn, from Equation (10.38), A,2 > А, л. Thus, the second throat must always be larger than the first throat. Only in the case of an ideal isentropic diffuser, where p0 — constant, would A,2 = A, and we have already discussed the impossibility of such an ideal diffuser.

Equation (10.38) is a useful relation to size the second throat relative to the first throat if we know the total pressure ratio across the tunnel. In the absence of such information, for the preliminary design of supersonic wind tunnels, the total pressure ratio across a normal shock is assumed.

For a given wind tunnel, if At 2 is less than the value given by Equation (10.38), the diffuser will “choke”; that is, the diffuser cannot pass the mass flow coming from the isentropic, supersonic expansion through the nozzle. In this case, nature adjusts the flow through the wind tunnel by creating shock waves in the nozzle, which in turn reduce the Mach number in the test section, producing weaker shocks in the diffuser with an attendant overall reduction in the total pressure loss; that is, nature adjusts the total pressure loss such that po. i/Au = Рол/Рв satisfies Equation (10.38). Sometimes this adjustment is so severe that a normal shock stands inside the nozzle, and the flow through the test section and diffuser is totally subsonic. Obviously, this choked situation is not desirable because we no longer have uniform flow at the desired Mach number in the test section. In such a case, the supersonic wind tunnel is said to be unstarted. The only way to rectify this situation is to make Al2/A, a large enough so that the diffuser can pass the mass flow from the isentropic expansion in the nozzle, that is, so that Equation (10.38) is satisfied along with a shock-free isentropic nozzle expansion.

As a general concluding comment, the basic concepts and relations discussed in this chapter are not limited to nozzles, diffusers, and supersonic wind tunnels. Rather, we have been discussing quasi-one-dimensional flow, which can be applied in many applications involving flow in a duct. For example, inlets on jet engines, which diffuse the flow to lower speeds before entering the engine compressor, obey the same principles. Also, a rocket engine is basically a supersonic nozzle designed to optimize the thrust from the expanded jet. The applications of the ideas presented in this chapter are numerous, and you should make certain that you understand these ideas before progressing further.

In Section 1.2, we subdivided aerodynamics into external and internal flows. You are reminded that the material in this chapter deals exclusively with internal flows.

Example 10.4 I For the preliminary design of a Mach 2 supersonic wind tunnel, calculate the ratio of the diffuser throat area to the nozzle throat area.

Solution

Assuming a normal shock wave at the entrance of the diffuser (for starting), from Appendix B, po. i/Poj = 0.7209 for M = 2.0. Hence, from Equation (10.38),

d^ = m = _L_ =| 1.387

At, і Рол 0.7209 —’——-

Hypersonics and Computational Fluid Dynamics

The design of hypersonic vehicles today is greatly dependent on the use of computa­tional fluid dynamics, much more so than the design of vehicles for any other flight regime. The primary reason for this is the lack of experimental ground test facilities that can simultaneously simulate the Mach numbers, Reynolds numbers, and high – temperature levels associated with hypersonic flight. For such simulation, CFD is the primary tool. Reflecting once again on the philosophy illustrated in Figure 2.44, in the realm of hypersonic flow the three partners are not quite equal. Pure experimental work in hypersonics usually involves tests at either the desired Mach number, the desired Reynolds number, or the desired temperature level, but not all at the same time nor in the same test facilities. As a result, experimental data for the design of hypersonic vehicles is a patchwork of different data taken in different facilities under different conditions. Moreover, the data are usually incomplete, especially for the high-temperature effects, which are difficult to simulate in a wind tunnel. The de­signer must then do his or her best to piece together the information for the specified design conditions. The next partner shown in Figure 2.44, pure theory, is greatly hampered by the nonlinear nature of hypersonic flow, hence making mathematical solutions intractable. In addition, the proper inclusion of high-temperature chemi­cally reacting flows in any pure theory is extremely difficult. For these reasons, the third partner shown in Figure 2.44, computational fluid dynamics, takes on a domi-

nant role. The numerical calculation of both inviscid and viscous hypersonic flows, including all the high-temperature effects discussed in Section 14.2, has been a major thrust of CFD research and design application since the 1960s. Indeed, hypersonics has paced the development of CFD since its beginning.

As an example of CFD applied to a hypersonic flight vehicle appropriate to this chapter, consider the space shuttle shown in Figure 14.15. A numerical solution of the three-dimensional inviscid flow field around the shuttle was carried out by Maus et al. in Reference 78. They made two sets of calculations, one for a perfect gas with у — 1.4, and one assuming chemically reacting air in local chemical equilibrium. The freestream Mach number was 23 in both cases. The CFD technique used for these calculations involved a time-dependent solution of the flow in the blunt nose region, patterned after our discussion in Section 13.5, and starting beyond the sonic line a downstream marching approach patterned after our discussion in Section 13.4. The calculated surface pressure distributions along the windward centerline of the space shuttle for both the perfect gas case (the circles) and the chemically reacting case (the triangles) are shown in Figure 14.16. The expansion around the nose, the pressure plateau over the relatively flat bottom surface, and the further expansion over the slightly inclined back portion of the body, are all quite evident. Also note that there is little difference in the pressure distributions between the two cases; this is an example of the more general result that pressure is usually the flow variable least affected by chemically reacting effects.

It is interesting to note, however, that a flight characteristic as mundane as the vehicle pitching moment coefficient is affected by chemically reacting flow effects. Close examination of Figure 14.16 shows that, for the chemically reacting flow, the

Figure 14.15 Space shuttle geometry.

pressures are slightly higher on the forward part of the shuttle, and slightly lower on the rearward part. This results in a more positive pitching moment. Since the moment is the integral of the pressure through a moment arm, a slight change in pressure can substantially affect the moment. This is indeed the case here, as shown in Figure 14.17, which is a plot of the resulting calculated pitching moment as a function of angle of attack for the space shuttle. Clearly, the pitching moment is substantially greater for the chemically reacting case. The work by Maus et al. was the first to point out this effect on pitching moment, and it serves to reinforce the importance of high-temperature flows on hypersonic aerodynamics. It also serves to reinforce the importance of CFD in the analysis of hypersonic flows. The predicted pitching mo­ment used for the space shuttle design came from “cold-flow” wind tunnel tests which did not simulate the high-temperature effects, that is, the designers used data for a perfect gas with у = 1.4 obtained in the wind tunnel. This is represented by the lower curve in Figure 14.7. The early flight experience with the shuttle indicated a much higher pitching moment at hypersonic speeds than predicted, which required that the body flap deflection for trim to be more than twice that predicted—an alarming situa­tion. The reason for this is now known; the actual flight environment encountered by the shuttle at high Mach numbers was that of a high-temperature chemically reacting flow—the situation reflected in the upper curve in Figure 14.17. The difference in the pitching moment between the two curves in Figure 14.17 is enough to account for the unexpected extra body flap deflection required to trim the shuttle. Although these CFD results were obtained well after the design of the shuttle, they serve to underscore the importance of CFD to present and future hypersonic vehicle designs.

Figure 1 4.1 7 Predicted pitching moment coefficient for the space shuttle; comparison between a calorically perfect gas and equilibrium air calculations. {Source: Maus et al., Reference 78.)

14.8 Summary

Only a few of the basic elements of hypersonic flow are presented here, with special emphasis on newtonian flow results. Useful information on hypersonic flows can be extracted from such results. We have derived the basic newtonian sine-squared law;

Cp = 2 sin2 9 [14.4]

and used this result to treat the case of a hypersonic flat plate in Section 14.4. We also obtained the limiting form of the oblique shock relations as Mx —»■ oo, that is, the hypersonic shock relations. From these relations, we were able to examine the significance of newtonian theory more thoroughly, namely, Equation (14.4) becomes an exact relation for a hypersonic flow in the combined limit of Мж —»• oo and у —у 1. Finally, these hypersonic shock relations illustrate the existence of the Mach number independence principle.

Problems

1. Repeat Problem 9.13 using (ia) Newtonian theory (b) Modified newtonian theory

Compare these results with those obtained from exact shock-expansion theory (Problem 9.13). From this comparison, what comments can you make about the

accuracy of newtonian and modified newtonian theories at low supersonic Mach numbers?

2. Consider a flat plate at a = 20° in a Mach 20 freestream. Using straight newto­nian theory, calculate the lift – and wave-drag coefficients. Compare these results with exact shock-expansion theory.

3. Consider a hypersonic vehicle with a spherical nose flying at Mach 20 at a standard altitude of 150,000 ft, where the ambient temperature and pressure are 500°R and 3.06 lb/ft2, respectively. At the point on the surface of the nose located 20° away from the stagnation point, estimate the: (a) pressure, (b) temperature,

(c) Mach number, and (d) velocity of the flow.

Stagnation Point Aerodynamic Heating

Contrary to what you might think, even though the flow velocity is zero at a stagnation point, the boundary layer at the stagnation point can be defined and has a finite thickness. The flow conditions at the edge of the stagnation point boundary layer are given by the inviscid solution for a stagnation point; in particular, at the boundary-layer

edge, the velocity is zero and the temperature is the total temperature, that is, ue = 0 and Te = To. This is shown in Figure 18.10. Moreover, along the vertical line in the ^-direction shown in Figure 18.10, и = 0 at every point inside the boundary layer. However, the ratio (и/и,,) — (0/0) is an indeterminant form that has a finite value at each point in the boundary layer. As in the case of the flat plate solutions discussed in Sections 18.2 and 18.3, we define a function /(r?) such that (u/ue) = and

/’ has a definite profile through the boundary layer. Indeed, we can define the edge of the boundary layer as the point where (u/ue) = /'(fl) = 0.99. Finally, we note that the shear stress at the wall at the stagnation point (point A in Figure 18.10) is zero. This not only comes out of the solution of the boundary layer equations, but it is obvious by inspection. Along the wall above point A the shear stress acts upward, and below point A it acts downward. Hence right at point A the shear stress must go through zero.

If the above discussion sounds rather theoretical, the temperature profile through the stagnation point boundary layer is easier to visualize. The temperature at the outer edge is the total temperature 7b. The temperature at the wall at r] — 0 is Tw. Hence, there is a temperature profile that exists in the normal direction through the stagnation point boundary layer. The heat transfer at the stagnation point is given by

the temperature gradient at point A, namely,

[18.58]

The practical purpose of a stagnation point boundary-layer solution is to calculate the heat transfer, qw.

The boundary-layer equations, Equation (17.28)—(17.31) applied at the stagna­tion point region are transformed using a version of the transformation described in Section 18.3, namely,

where h is the static enthalpy (since и = 0, the static and total enthalpies are the same). This leads to the stagnation point boundary layer equations given below. For a detailed derivation of these equations, see for example Chapter 6 of Reference 55.

[18.64]

where C = (pfi/pefie). Equations (18.63) and (18.64) are the governing equations for a compressible, stagnation-point boundary layer. Examining these equations, we see no § -dependency. Hence, the stagnation point boundary layer is a self-similar case.

Numerical solutions to Equations (18.63) and (18.64) can be obtained by the “shooting technique” as described earlier in the flat plate case. There is nothing to be gained in going through the details at this stage of our discussion. Instead, we simply state the result of solving Equations (18.63) and (18.64), correlated in the following expression obtained from Reference 82:

Cylinder: qw = 0.57 Рґ0’6(редг)І/2./^(/і„ – hw) [18.65]

V dx

If we had considered an axisymmetric body, the original transformation given by Equations (18.59) and (18.60) would have been slightly modified as follows:

1 = peueper2 dx Jo

where r is the vertical coordinate measured from the centerline, as shown in Fig­ure 18.10. Equations (18.66) and (18.67) lead to equations for the axisymmetric stagnation point almost identical to Equation (18.63) and (18.64), namely,

[18.68]

[18.69] where C = (p/і/pe/ie). In turn, the resulting heat transfer expression is (Refer­ence 82):

[18.70]

Compare Equation (18.65) for the two-dimensional cylinder with Equation

(18.70) for the axisymmetric sphere. The equations are the same except for the leading coefficient, which is higher for the sphere. Everything else being the same, this demonstrates that stagnation point heating to a sphere is larger than to a two­dimensional cylinder. Why? The answer lies in a basic difference between two – and three-dimensional flows. In a two-dimensional flow, the gas has only two directions to move when it encounters a body—up or down. In contrast, in an axisymmetric flow, the gas has three directions to move—up, down, and sideways—and hence the flow is somewhat “relieved,” that is, in comparing two – and three-dimensional flows over bodies with the same longitudinal section (such as a cylinder and a sphere), there is a well-known three-dimensional relieving effect for the three-dimensional flow. As a consequence of this relieving effect, the boundary-layer thickness 8 at the stagnation point is smaller for the sphere than for the cylinder. In turn, the temperature gradient at the wall, (dT/dy)w, which is of the order of (Te/8), is larger for the sphere. Since qw = k(dT/dy)w, then qw is larger for the sphere. This confirms the comparison between Equations (18.65) and (18.70).

The above results for aerodynamic heating to a stagnation point have a stunning impact on hypersonic vehicle design, namely, they impose the requirement for the vehicle to have a blunt, rather than a sharp, nose. To see this, consider the velocity gradient, due/dx, which appears in Equations (18.65) and (18.70). From Euler’s equation applied at the edge of the boundary layer

Assuming a newtonian pressure distribution over the surface, we have from Equation

(14.4)

Cp = 2 sin2 в

where 9 is defined as the angle between a tangent to the surface and the freestream direction. If we define ф as the angle between the normal to the surface and the freestream, then Equation (14.4) can be written as

Cp =2 cos2 ф

From the definition of Cp, Equation (18.73) becomes

Differentiating Equation (18.74), we obtain

dPe Л Л ■ л^

‘ = —4<7oc cos ф sin 0 —

dx ax

Combining Equations (18.72) and (18.75), we have

Equation (18.76) is a general result which applies at all points along the body. Now consider the stagnation-point region, as sketched in Figure 18.11. In this region, let Ax be a small increment of surface distance above the stagnation point, corresponding to the small change in ф, Аф. The inviscid velocity variation in the stagnation region can be shown to be

Also, in the stagnation region ф is small, hence, from Figure 18.11,

cos ф ~ 1

[18.78]

Ax

[18.79]

sin ф ~ ф % Аф ~ —–

R

dф 1

[18.80]

dx R

where R is the local radius of curvature of the body at the stagnation point. Finally, at the stagnation point, Equation (18.73) becomes

Substituting Equations (18.77)—(18.81) into (18.76), we have

/duA2 _ 2(pe – poo) /Ax / J dx ) pe Ay R ) R)

Examine Equations (18.65) and (18.70) in light of Equation (18.82). We see that

[18.83]

This states that stagnation-point heating varies inversely with the square root of the nose radius; hence, to reduce the heating, increase the nose radius. This is the reason why the nose and leading edge regions of hypersonic vehicles are blunt; otherwise, the severe aerothermal conditions in the stagnation region would quickly melt a sharp leading edge.

Return to Section 1.1 and review our qualitative discussion contrasting the aero­dynamic heating for slender and blunt reentry vehicles. There we argued on a qual­itative basis that to minimize aerodynamic heating a blunt nose must be used. We have now quantitatively proven this fact with the derivation of Equation (18.83).

The fact that qw is inversely proportional to ~/R is experimentally verified in Figure 18.12, obtained from Reference 83. Here, various sets of experimental data for Ся at the stagnation point are plotted versus Reynolds number based on nose diameter; the abscissa is essentially proportional R. This is a log-log plot, and the data exhibits a slope of —0.5, hence verifying that qw ос 1 //R.

Figure 18.13 Stagnation point Stanton number versus Re based on nose radius. (Source: Koppenwallner, Reference 83.)

Historical Note: High-Speed Airfoils—Early Research and Development

Twentieth-century aerodynamics does not have the exclusive rights to the observation of the large drag rise on bodies flying at near the speed of sound; rather, in the eigh­teenth century the Englishman Benjamin Robins, inventor of the ballistic pendulum, reported that “the velocity at which the body shifts its resistance (from a V2 to a Vі relation) is nearly the same with which sound is propagated through air.” His statement was based on a large number of experiments during which projectiles were fired into his ballistic pendulum. However, these results had little relevance to the early aerodynamicists of this century, who were struggling to push aircraft speeds to 150 mi/h during and just after World War I. To these people, flight near the speed of sound was just fantasy.

With one exception! World War I airplanes such as the Spad and Nieuport had propeller blades where the tips were moving at near the speed of sound. By 1919, British researchers had already observed the loss in thrust and large increase in blade drag for a propeller with tip speeds up to 1180 ft/s—slightly above the speed of sound. To examine this effect further, F. W. Caldwell and E. N. Fales, both engineers at the U. S. Army’s Engineering Division at McCook Field near Dayton, Ohio (the forerunner of the massive Air Force research and development facilities at Wright – Patterson Air Force Base today), conducted a series of high-speed airfoil tests. They designed and built the first high-speed wind tunnel—a facility with a 14-in-diameter test section capable of velocities up to 675 ft/s. In 1918, they conducted the first wind-tunnel tests involving the high-speed flow over a stationary airfoil. Their results showed large decreases in lift coefficient and major increases in drag coeffiient for the thicker airfoils at angle of attack. These were the first measured “compressibility effects” on an airfoil in history. Caldwell and Fales noted that such changes occurred at a certain air velocity, which they denoted as the “critical speed”—a term that was to evolve into the critical Mach number at a later date. It is interesting to note that Orville Wright was a consultant to the Army at this time (Wilbur had died prematurely in 1912 of typhoid fever) and observed some of the Caldwell and Fales tests. However, a fundamental understanding and explanation of this critical-speed phenomenon was completely lacking. Nobody at that time had even the remotest idea of what was really happening in this high-speed flow over the airfoil.

Members of the National Advisory Committee for Aeronautics were well aware of the Caldwell-Fales results. Rather than let the matter die, in 1922 the NACA contracted with the National Bureau of Standards (NBS) for a study of high-speed flows over airfoils, with an eye toward improved propeller sections. The work at NBS included the building of a high-speed wind tunnel with a 12-in-diameter test section, capable of producing a Mach number of 0.95. The aerodynamic testing was performed by Lyman J. Briggs (soon to become director of NBS) and Hugh Dryden (soon to become one of the leading aerodynamicists of the twentieth century). In addition to the usual force data, Briggs and Dryden also measured pressure distributions over the airfoil surface. These pressure distributions allowed more insight into the nature of the flow and definitely indicated flow separation on the top surface of the airfoil. We now know that such flow separation is induced by a shock wave, but these early researchers did not at that time know about the presence of such shocks.

During the same period, the only meaningful theoretical work on high-speed airfoil properties was carried out by Ludwig Prandtl in Germany and Hermann Glauert in England—work which led to the Prandtl-Glauert compressibility correction, given by Equation (11.51). As early as 1922, Prandtl is quoted as stating that the lift coefficient increased according to (1 — M^)-1^2; he mentioned this conclusion in his lectures at Gottingen, but without written proof. This result was mentioned again 6 years later by Jacob Ackeret, a colleague of Prandtl, in the famous German series Handbuch der Physik, again without proof. Subsequently, in 1928 the concept was formally established by Hermann Glauert, a British aerodynamicist working for the Royal Aircraft Establishment. (See Chapter 9 of Reference 21 for a biographical sketch of Glauert.) Using only six pages in the Proceedings of the Royal Society,

vol. 118, p. 113, Glauert presented a derivation based on linearized small-perturbation theory (similar to that described in Section 11.4) which confirmed the (1 – A/^)-1”2 variation. In this paper, entitled “The Effect of Compressibility on the Lift of an Airfoil,” Glauert derived the famous Prandtl-Glauert compressibility correction, given here as Equations (11.51) to (11.53). This result was to stand alone, unaltered, for the next 10 years.

Hence, in 1930 the state of the art of high-speed subsonic airfoil research was characterized by experimental proof of the existence of the drag-divergence phe­nomenon, some idea that it was caused by flow separation, but no fundamental under­standing of the basic flow field. In turn, there was virtually no theoretical background outside of the Prandtl-Glauert rule. Also, keep in mind that all the above work was paced by the need to understand propeller performance, because in that day the only component of airplanes to encounter compressibility effects was the propeller tips.

All this changed in the 1930s. In 1928, the NACA had constructed its first rudi­mentary high-speed subsonic wind tunnel at the Langley Aeronautical Laboratory, utilizing a 1 – ft-diameter test section. With Eastman Jacobs as tunnel director and John Stack as the chief researcher, a series of tests was run on various standard airfoil shapes. Frustrated by their continual lack of understanding about the flow field, they turned to optical techniques, following in the footsteps of Ernst Mach (see Section 9.9). In 1933, they assembled a crude schlieren optical system consisting of 3-in-diameter reading- glass-quality lenses and a short-duration-spark light source. In their first test using the schlieren system, dealing with flow over a cylinder, the results were spectacular. Shock waves were seen, along with the resulting flow separation. Visitors flocked to the wind tunnel to observe the results, including Theodore Theodorsen, one of the ranking theoretical aerodynamicists of that period. An indicator of the psychology at that time is given by Theodorsen’s comment that since the freestream flow was sub­sonic, what appeared as shock waves must be an “optical illusion.” However, Eastman Jacobs and John Stack knew differently. They proceeded with a major series of airfoil testing, using standard NACA sections. Their schlieren pictures revealed the secrets of flow over the airfoils above the critical Mach number. (See Figure 131b and its at­tendant discussion of such supercritical flow.) In 1935, Jacobs traveled to Italy, where he presented results of the NACA high-speed airfoil research at the fifth Volta Confer­ence (see Section 7.1). This is the first time in history that photographs of the transonic flow field over standard-shaped airfoils were presented in a large public forum.

During the course of such work in the 1930s, the incentive for high-speed aero­dynamic research shifted from propeller applications to concern about the airframe of the airplane itself. By the mid-1930s, the possibility of the 550 mi/h airplane was more than a dream—reciprocating engines were becoming powerful enough to consider such a speed regime for propeller-driven aircraft. In turn, the entire airplane itself (wings, cowling, tail, etc.) would encounter compressibility effects. This led to the design of a large 8-ft high-speed tunnel at Langley, capable of test-section velocities above 500 mi/h. This tunnel, along with the earlier 1 – ft tunnel, established the NACA’s dominance in high-speed subsonic research in the late 1930s.

In the decade following 1930, the picture had changed completely. By 1940, the high-speed flow over airfoils was relatively well understood. During this period,

Stack and Jacobs had not only highlighted the experimental aspects of such high­speed flow, but they also derived the expression for C;)cr as a function of Mcr given by Equation (11.60), and had shown how to estimate the critical Mach number for a given airfoil as discussed in Section 11.6. Figure 11.24 shows some representative schlieren photographs taken by the NACA of the flow over standard NACA airfoils. Although

NACA 64A006 NACA 64A009 NACA 64A012

Figure 1 1.24 Schlieren pictures and pressure distributions for transonic flows over several NACA airfoils. These pictures were taken by the NACA in 1949. (Source: John V. Becker, "The High-Speed Frontier," NASA SP-445, 1980.)

these photographs were taken in 1949, they are similar to the results obtained by Stack and Jacobs in the 1930s. Superimposed on these photographs are the measured pressure distributions over the top (solid curve) and bottom (dashed curve) surfaces of the airfoil. Study these pictures carefully. Moving from bottom to top, you can see the influence of increasing freestream Mach number, and going from left to right, you can observe the effect of increasing airfoil thickness. Note how the shock wave moves downstream as Мж is increased, finally reaching the trailing edge at = 1.0. For this case, the top row of pictures shows almost completely supersonic flow over the airfoil. Note also the large regions of separated flow downstream of the shock waves for the Mach numbers of 0.79, 0.87, and 0.94—this separated flow is the primary reason for the large increase in drag near Mach 1. By 1940, it was well understood that the almost discontinuous pressure increase across the shock wave creates a strong adverse pressure gradient on the airfoil surface, and this adverse pressure gradient is responsible for separating the flow.

The high-speed airfoil research program continues today within NASA. It led to the supercritical airfoils in the 1960s (see Sections 11.9 and 11.12). It has produced a massive effort in modem times to use computational techniques for theoretically solving the transonic flow over airfoils. Such efforts are beginning to be successful, and in many respects, today we have the capability to design transonic airfoils on the computer. However, such abilities today have roots which reach all the way back to Caldwell and Fales in 1918.

For a more detailed account of the history of high-speed airfoil research, you are encouraged to read the entertaining story portrayed by John V. Becker in The High-Speed Frontier, NASA SP-445, 1980.

Negligible Viscous Dissipation

To some extent, we have already discussed this case in regard to the local heat flux at any point within the flow. If ue is very small, hence г is very small, then the amount of viscous dissipation is negligibly small, and Equation (16.21) becomes

/he – hw Pr V D

Clearly, for this case, the heat flux is constant across the flow. Moreover, the enthalpy profile given by Equation (16.16) becomes

h = hw + {he — hw)~ [16.27]

Since h = cpT, the temperature profile is identical to the enthalpy profile:

T = Tw + (Te — Tw)~ [16.28]

Note that the temperature varies linearly across the flow, as sketched in Figure 16.4. The case shown here is for the upper wall at a higher temperature than the lower wall. The heat transfer at the lower wall is obtained from Equation (16.24) with a negligible

The heat transfer at the upper wall is similarly obtained as

Equations (16.29) and (16.30) are identical; this is no surprise, since we have already shown that the heat flux is constant across the flow, as shown by Equation (16.26), and therefore the heat transfer at both walls should be the same. Equations (16.29) and (16.30) can also be written in terms of temperature as

[16.31]

Examining Equations (16.29) to (16.31), we can make some conclusions which can be generalized to most viscous flow problems, as follows:

1. Everything else being equal, the larger the temperature difference across the viscous layer, the greater the heat transfer at the wall. The temperature difference (Te — Tw) or the enthalpy difference (he — hw) takes on the role of a “driving potential” for heat transfer. For the special case treated here, the heat transfer at the wall is directly proportional to this driving potential.

2. Everything else being equal, the thicker the viscous layer (the larger D is) the smaller the heat transfer is at the wall. For the special case treated here, qw is inversely proportional to D.

3. Heat flows from a region of high temperature to low temperature. For negligible viscous dissipation, if the temperature at the top of the viscous layer is higher than that at the bottom, heat flows from the top to the bottom. In the case sketched in Figure 16.4, heat is transferred from the upper plate into the fluid, and then is transferred from the fluid to the lower plate.

Flow over a Complete Airplane

We end this section by noting a history-making calculation. In Reference 57, a solution of the complete Navier-Stokes equations for the flow field over an entire airplane was reported—the first such calculation ever made. In this work, Shang and Scherr carried out a time-dependent solution using MacCormack’s method—just as we have discussed it in Section 16.4. See Reference 57 for the details. Also, a lengthy

Figure 20.3 Curvilinear, boundary-fitted finite-

difference grid for the solution of the flow over an airfoil. (Source: Kothari and Anderson, Reference 56.)

description of this work can be found in Chapter 8 of Reference 55. Shang and Scherr applied their calculation to the hypersonic viscous flow over the X-24C hypersonic test vehicle. To illustrate the results, the surface streamline pattern is shown in Figure 20.5, as calculated in Reference 57. In reality, since the flow velocity is zero at the surface in a viscous flow (the no-slip condition), the lines shown in Figure 20.5 are the surface shear stress directions.

Oblique Shock Relations

Consider the oblique shock wave sketched in Figure 9.6. The angle between the shock wave and the upstream flow direction is defined as the wave angle, denoted by fl. The upstream flow (region 1) is horizontal, with a velocity V) and Mach number Mi – The downstream flow (region 2) is inclined upward through the deflection angle в and has velocity V2 and Mach number М2. The upstream velocity Vt is split into components tangential and normal to the shock wave, w and ji, respectively, with the associated tangential and normal Mach numbers Mtt and Mn,, respectively. Similarly, the downstream velocity is split into tangential and normal components W2 and /І2, respectively, with the associated Mach numbers M,2 and M„,2.

Consider the control volume shown by the dashed lines in the upper part of Figure 9.6. Sides a and d are parallel to the shock wave. Segments b and c follow the upper streamline, and segments e and / follow the lower streamline. Let us apply the integral form of the conservation equations to this control volume, keeping in mind that we are dealing with a steady, inviscid, adiabatic flow with no body forces. For these assumptions, the continuity equation, Equation (2.48), becomes

pV • dS = 0

s

This surface integral evaluated over faces a and d yields — pU A + p2u2A2, where A1 = A2 — area of faces a and d. The faces b, c, e, and / are parallel to the velocity,

Oblique Shock Relations

and hence contribute nothing to the surface integral (i. e., V • dS = 0 for these faces). Thus, the continuity equation for an oblique shock wave is

Подпись: or Подпись: PU — P2U2 Подпись: [9.2]

PUA + P2U2A2 = 0

Keep in mind that u and U2 in Equation (9.2) are normal to the shock wave.

The integral form of the momentum equation, Equation (2.64), is a vector equa­tion. Hence, it can be resolved into two components, tangential and normal to the shock wave. First, consider the tangential component, keeping in mind the type of flow we are considering:

Подпись: (pV • dS)wПодпись: sOblique Shock RelationsПодпись: [9.3](P d. S’)tangCntia|

In Equation (9.3), w is the component of velocity tangential to the wave. Since dS is perpendicular to the control surface, then (p dShangentiai over faces a and d is zero. Also, since the vectors p dS on faces b and / are equal and opposite, the pressure integral in Equation (9.3) involves two tangential forces that cancel each other over faces b and /. The same is true for faces c and e. Hence, Equation (9.3) becomes

Подпись: [9.4]-(puxAx)w + (P2U2A2)w2 = 0

Dividing Equation (9.4) by Equation (9.2), we have

Подпись: u>i = W2[9.5]

Equation (9.5) is an important result; it states that the tangential component of the flow velocity is constant across an oblique shock.

The normal component of the integral momentum equation is, from Equation (2.64),

Подпись: [9.6]’jlj’ (pV • dS)M = (pdS)normal

s s

Here, the pressure integral evaluated over faces a and d yields the net sum —pA + P2A2. Once again, the equal and opposite pressure forces on b and / cancel, as do those on c and e. Hence, Equation (9.6) becomes, for the control volume shown in Figure 9.6,

-(ріИіАі)И] + (P2U2A2)U2 = -(-P1A1 + P2A2)

Подпись: Pl + Pu] = P2 + P2U Подпись: [9.7]

Since А і = A2, this becomes

Again, note that the only velocities appearing in Equation (9.7) are the components normal to the shock.

Finally, consider the integral form of the energy equation, Equation (2.95). For our present case, this can be written as

Подпись: pV-dSOblique Shock Relations[9.8]

Again noting that the flow is tangent to faces b, c, /, and e, and hence V • dS = 0 on these faces, Equation (9.8) becomes, for the control volume in Figure 9.6,

~P (єі + у ^ щАх + P2 (e2 + y^ u2A2 = —{—pUA + P2U2A2) [9.9]

Подпись: or Oblique Shock Relations Oblique Shock Relations Подпись: = 0 Подпись: [9.10]

Collecting terms in Equation (9.9), we have

Dividing Equation (9.10) by (9.2), we have

V2 V2

hi + – j – = h2+ – j – [9.11]

Since h + V2/2 = h0, we have again the familiar result that the total enthalpy is constant across the shock wave. Moreover, for a calorically perfect gas, ho — cpTo‘, hence, the total temperature is constant across the shock wave. Carrying Equation

(9.11) a bit further, note from Figure 9.6 that V2 = u2 + w2. Also, from Equation

(9.5) , we know that иц = w2. Hence,

V2 — V22 = (и2 + w2) — (и2 + w) = u — u

Thus, Equation (9.11) becomes

Oblique Shock Relations[9.12]

Let us now gather our results. Look carefully at Equations (9.2), (9.7), and (9.12). They are the continuity, normal momentum, and energy equations, respectively, for an oblique shock wave. Note that they involve the normal components only of velocity и і and u2, the tangential component w does not appear in these equations. Hence, we deduce that changes across an oblique shock wave are governed only by the component of velocity normal to the wave.

Again, look hard at Equations (9.2), (9.7), and (9.12). They are precisely the governing equations for a normal shock wave, as given by Equations (8.2), (8.6), and

(8.10) . Hence, precisely the same algebra as applied to the normal shock equations in Section 8.6, when applied to Equations (9.2), (9.7), and (9.12), will lead to identical expressions for changes across an oblique shock in terms of the normal component of the upstream Mach number Mn . Note that

Подпись: [9.13]Mn. = M sin P

Подпись: Pi Pi El pi Oblique Shock Relations Подпись: [9.14] [9.15] [9.16]

Hence, for an oblique shock wave, with Mn. given by Equation (9.13), we have, from Equations (8.59), (8.61), and (8.65),

The temperature ratio T2/T follows from the equation of state:

Подпись:II – ElEL

T pi p2

Note that Mn 2 is the normal Mach number behind the shock wave. The downstream Mach number itself, M2, can be found from Mn2 and the geometry of Figure 9.6 as

Подпись:Mn, 1

sin(/l – 9)

Examine Equations (9.14) to (9.17). They state that oblique shock-wave proper­ties in a calorically perfect gas depend only on the normal component of the upstream Mach number Mn.. However, note from Equation (9.13) that Mn. depends on both M and f. Recall from Section 8.6 that changes across a normal shock wave depend on one parameter only—the upstream Mach number M ,. In contrast, we now see that changes across an oblique shock wave depend on two parameters—say, M and

ft. However, this distinction is slightly moot because in reality a normal shock wave is a special case of oblique shocks where ft = ir/2.

Equation (9.18) introduces the deflection angle в into our oblique shock analysis; we need в to be able to calculate М2. However, 9 is not an independent, third parameter; rather, 9 is a function of M and ft, as derived below. From the geometry of Figure 9.6,

Mi

tan fi = —

[9.19]

Wi

and

М2

tan(£ – 0) = —

[9.20]

W2

Dividing Equation (9.20) by (9.19), recalling that wx = w2, and invoking the conti­nuity equation, Equation (9.2), we obtain

tan(/l – 9) u2 pi

tan ft и і p2

Подпись: [9.22]Combining Equation (9.21) with Equations (9.13) and (9.15), we obtain tan(/i – 9) 2+ (y — l)M2sin2/J

tan ft (y + l)Mf sin2 ft

Oblique Shock Relations Подпись: [9.23]

which gives в as an implicit function of M and ft. After some trigonometric substi­tutions and rearrangement, Equation (9.22) can be cast explicitly for 9 as

Equation (9.23) is an important equation. It is called the в-fi-M relation, and it specifies 9 as a unique function of M and ft. This relation is vital to the analysis of oblique shock waves, and results from it are plotted in Figure 9.7 for у — 1.4. Examine this figure closely. It is a plot of wave angle versus deflection angle, with the Mach number as a parameter. The results given in Figure 9.7 are plotted in some detail—this is a chart which you will need to use for solving oblique shock problems.

Figure 9.7 illustrates a wealth of physical phenomena associated with oblique shock waves. For example:

1. For any given upstream Mach number M1, there is a maximum deflection angle 0max. If the physical geometry is such that 9 > 0max, then no solution exists for a straight oblique shock wave. Instead, nature establishes a curved shock wave, detached from the corner or the nose of a body. This is illustrated in Figure

9.8. Here, the left side of the figure illustrates flow over a wedge and a concave corner where the deflection angle is less than 0max for the given upstream Mach number. Therefore, we see a straight oblique shock wave attached to the nose of the wedge and to the comer. The right side of Figure 9.8 gives the case where the deflection angle is greater than 0max; hence, there is no allowable straight oblique shock solution from the theory developed earlier in this section. Instead, we have

Oblique Shock Relations

Figure 9.7 Oblique shock properties: у = 1.4. The в-fi-M diagram. (Source: NACA Report 1 135, Ames Research Staff, "Equations, Tables and Charts for Compressible Flow," 1953.)

Shock-wove angle 0, degree*

Figure 9.7 (continued) (Source: NACA Report 1 135, Ames Research Staff, "Equations, Tables and Charts for Compressible Flow," 1953.)

a curved shock wave detached from the nose of the wedge or from the comer. Return to Figure 9.7, and note that the value of 0max increases with increasing M. Hence, at higher Mach numbers, the straight oblique shock solution can exist at higher deflection angles. However, there is a limit; as M approaches infinity, #max approaches 45.5° (for у — 1.4).

2. For any given 9 less than 6>max, there are two straight oblique shock solutions for a given upstream Mach number. For example, if M = 2.0 and в = 15°, then from Figure 9.7, f5 can equal either 45.3 or 79.8°. The smaller value of f) is called the weak shock solution, and the larger value of fi is the strong shock solution. These two cases are illustrated in Figure 9.9. The classifications “weak” and “strong” derive from the fact that for a given M, the larger the wave angle, the larger the normal component of upstream Mach number МпЛ, and from Equation (9.16) the larger the pressure ratio pifp. Thus, in Figure 9.9, the higher-angle shock wave will compress the gas more than the lower-angle shock wave, hence the terms “strong” and “weak” solutions. In nature, the weak shock solution usually prevails. Whenever you see straight, attached oblique shock waves, such as sketched at the left of Figure 9.8, they are almost always the weak shock solution. It is safe to make this assumption, unless you have specific information to the contrary. Note in Figure 9.7 that the locus of points connecting all the values of 0max (the curve that sweeps approximately horizontally across the middle of Figure 9.7) divides the weak and strong shock solutions. Above

/

/

Figure 9.9 The weak and strong shock cases.

this curve, the strong shock solution prevails (as further indicated by the 9-fi-M curves being dashed); below this curve, the weak shock solution prevails (where the 9-fi-M curves are shown as solid lines). Note that slightly below this curve is another curve which also sweeps approximately horizontally across Figure 9.7. This curve is the dividing line above which М2 < 1 and below which М2 > 1. For the strong shock solution, the downstream Mach number is always subsonic M2 < 1. For the weak shock solution very near 0max, the downstream Mach number is also subsonic, but barely so. For the vast majority of cases involving the weak shock solution, the downstream Mach number is supersonic М2 > 1. Since the weak shock solution is almost always the case encountered in nature, we can readily state that the Mach number downstream of a straight, attached oblique shock is almost always supersonic.

3. If 9 = 0, then fi equals either 90° or /1. The case of /3 = 90° corresponds to a normal shock wave (i. e., the normal shocks discussed in Chapter 8 belong to the family of strong shock solutions). The case of fi — /z corresponds to the Mach wave illustrated in Figure 9.3b. In both cases, the flow streamlines experience no deflection across the wave.

4. (In all of the following discussions, we consider the weak shock solution exclu­sively, unless otherwise noted.) Consider an experiment where we have super­sonic flow over a wedge of given semiangle 9, as sketched in Figure 9.10. Now assume that we increase the freestream Mach number M. As M increases, we observe that /1 decreases. For example, consider 9 — 20° and M = 2.0, as shown on the left of Figure 9.10. From Figure 9.7, we find that fi = 53.3°. Now assume M is increased to 5, keeping 9 constant at 20°, as sketched on the right of Figure 9.10. Here, we find that fi = 29.9°. Interestingly enough, although this shock is at a lower wave angle, it is a stronger shock than the one on the left. This is because Mn is larger for the case on the right. Although fi is smaller, which decreases Mn , the upstream Mach number M is larger, which increases M„ 1 by an amount which more than compensates for the decreased /3. For example, note the values of Mnд and Рг! P given in Figure 9.10. Clearly, the Mach 5 case on the right yields the stronger shock wave. Hence, in general for attached shocks with a fixed deflection angle, as the upstream Mach number M increases, the wave angle fi decreases, and the shock wave becomes stronger.

Going in the other direction, as Mi decreases, the wave angle increases, and the shock becomes weaker. Finally, if M is decreased enough, the shock wave will become detached. For the case of в = 20° shown in Figure 9.10, the shock will be detached for Mi < 1.84.

5. Consider another experiment. Here, let us keep M fixed and increase the de­flection angle. For example, consider the supersonic flow over a wedge shown in Figure 9.11. Assume that we have M = 2.0 and в — 10°, as sketched at the left of Figure 9.11. The wave angle will be 39.2° (from Figure 9.7). Now assume that the wedge is hinged so that we can increase its deflection angle, keeping Mi constant. In such a case, the wave angle will increase, as shown on the right of Figure 9.11. Also, Mn will increase, and hence the shock will become stronger. Therefore, in general for attached shocks with a fixed upstream Mach number, as the deflection angle increases, the wave angle fi increases, and the shock becomes stronger. However, once в exceeds 0max, the shock wave will

become detached. For the case of M = 2.0 in Figure 9.11, this will occur when в > 23°.

The physical properties of oblique shocks just discussed are very important. Before proceeding further, make certain to go over this discussion several times until you feel perfectly comfortable with these physical variations.

Consider a supersonic flow with M = 2, p = 1 atm, and T = 288 K. This flow is deflected at a compression comer through 20°. Calculate M, p, T, p0, and T0 behind the resulting oblique shock wave.

Solution

From Figure 9.7, for Mi = 2 and в = 20°, ft = 53.4°. Hence, M„д = M sin ft = 2 sin 53.4° = 1.606. From Appendix B, for M„л = 1.60 (rounded to the nearest table entry),

Note: For oblique shocks, the entry for ро. г/P in Appendix В cannot be used to obtain po 2′, this entry in Appendix В is for normal shocks only and is obtained directly from Equation (8.80). In turn, Equation (8.80) is derived using (8.77), where M2 is the actual flow Mach number, not the normal component. Only in the case of a normal shock is this also the Mach number normal to the wave. Hence, Equation (8.80) holds only for normal shocks; it cannot be used for oblique shocks with Mi replaced by M„,i. For example, an incorrect calculation would be to use pop/Pi = 3.805 for МпЛ = 1.60. This gives />0 2 = 3.805 atm, a totally incorrect result compared with the correct value of 7.00 atm obtained above.

Example 9.2 | Consider an oblique shock wave with a wave angle of 30°. The upstream flow Mach number is 2.4. Calculate the deflection angle of the flow, the pressure and temperature ratios across the shock wave, and the Mach number behind the wave.

1. This is a fairly weak shock wave—only a 51 percent increase in pressure across the wave. Indeed, examining Figure 9.7, we find that this case is close to that of a Mach wave, where д = sin~'(l/M) = sin*1!^) = 24.6°. The shock­wave angle of 30° is not much larger than д; the deflection angle of 6.5° is also small—consistent with the relative weakness of the shock wave.

2. Only two properties need to be specified in order to define uniquely a given oblique shock wave. In this example, M and fi were those two properties. In Example 9.1, the specified M and 9 were the two properties. Once any two properties about the oblique shock are specified, the shock is uniquely defined. This is analogous to the case of a normal shock wave studied in Chapter 8. There, we proved that all the changes across a normal shock wave were uniquely defined by specifying only one property, such as M. However, implicit in all of Chapter 8 was an additional property, namely, the wave angle of a normal shock wave is 90°. Of course, a normal shock is simply one example of the whole spectrum of oblique shocks, namely, a shock with /1 = 90°. An examination of Figure 9.7 shows that the normal shock belongs to the family of strong shock solutions, as discussed earlier.

Note: Once again, the oblique shock is uniquely defined by two properties, in this case fi and Pi/Pi-

Consider a Mach 3 flow. It is desired to slow this flow to a subsonic speed. Consider two separate ways of achieving this: (1) the Mach 3 flow is slowed by passing directly through a normal shock wave; (2) the Mach 3 flow first passes through an oblique shock with a 40° wave angle, and then subsequently through a normal shock. These two cases are sketched in Figure 9.12. Calculate the ratio of the final total pressure values for the two cases, that is, the total pressure behind the normal shock for case 2 divided by the total pressure behind the normal shock for case 1. Comment on the significance of the result.

Solution

For case 1, at M = 3, we have, from Appendix B,

() = 0.3283

V Po, ) case 1

For case 2, we have МпЛ = M, sin ft = 3 sin 40° = 1.93. From Appendix B,

= 0.7535 and Mn 2 = 0.588

Po, ‘

From Figure 9.7, for M, = 3 and ft = 40°, we have the deflection angle в = 22°. Hence,

From Appendix B, for a normal shock with an upstream Mach number of 1.9, we have Р03/Po2 = 0.7674. Thus, for case 2,

Figure 9.12 Illustration for Example 9.4.

The result of Example 9.4 shows that the final total pressure is 76 percent higher for the case of the multiple shock system (case 2) in comparison to the single normal shock (case 1). In principle, the total pressure is an indicator of how much useful work can be done by the gas; this is described later in Section 10.4. Everything else being equal, the higher the total pressure, the more useful is the flow. Indeed, losses of total pressure are an index of the efficiency of a fluid flow—the lower the total pressure loss, the more efficient is the flow process. In this example, case 2 is more efficient in slowing the flow to subsonic speeds than case 1 because the loss in total pressure across the multiple shock system of case 2 is actually less than that for case 1 with a single, strong, normal shock wave. The physical reason for this is straightforward. The loss in total pressure across a normal shock wave becomes particularly severe as the upstream Mach number increases; a glance at the P0.2/P0.1 column in Appendix В attests to this. If the Mach number of a flow can be reduced before passing through a normal shock, the loss in total pressure is much less because the normal shock is weaker. This is the function of the oblique shock in case 2, namely, to reduce the Mach number of the flow before passing through the normal shock. Although there is a total pressure loss across the oblique shock also, it is much less than across a normal shock at the same upstream Mach number. The net effect of the oblique shock reducing the flow Mach number before passing through the normal shock more than makes up for the total pressure loss across the oblique shock, with the beneficial result that the multiple shock system in case 2 produces a smaller loss in total pressure than a single normal shock at the same freestream Mach number.

A practical application of these results is in the design of supersonic inlets for jet engines. A normal shock inlet is sketched in Figure 9.13a. Here, a normal shock forms ahead of the inlet, with an attendant large loss in total pressure. In contrast, an oblique shock inlet is sketched in Figure 9. 3b. Here, a central cone creates an oblique shock wave, and the flow subsequently passes through a relatively weak normal shock at the lip of the inlet. For the same flight conditions (Mach number and altitude), the total pressure loss for the oblique shock inlet is less than for a normal shock inlet. Hence, everything else being equal, the resulting engine thrust will be higher for the oblique shock inlet. This, of course, is why most modem supersonic aircraft have oblique shock inlets.

Predictor Step

First, predict the value of Fi+j by using a Taylor series where 3E/3f is evaluated at point (i, j). Denote this predicted value by Fi+j:

Fi+ij = F, j + ( ^ ) A*

In Equation (13.55), (3E/3f )jj is obtained from the continuity equation, Equation

(13.50) , using/orwa/г/ differences for the p derivatives; that is,

In Equation (13.56), all quantities on the right-hand side are known and allow the calculation of (3E/3f ), j which is, in turn, inserted into Equation (13.55). A similar procedure is used to find predicted values of G, H, and K, namely, G,+ij, Hi+ij, and Ki+ij, using forward differences in Equations (13.51) to (13.53). In turn, predicted values of the primitive variables, Pi+j, Рі+ij, etc., can be obtained from Equations (13.39a to d) and (13.44) to (13.46).

13.4.3 Corrector Step

The predicted values obtained above are used to obtain predicted values of the deriva­tive (3E/3f );+ij, using rearward differences in Equation (13.50):

/ЭЕ) _ ( Ц^уА fi’+i. i ~ Fi+ij-1 1 (pv)i+u ~ (pv)i+ij-1

/ ~ y, dx)i+lj Art ys Ap

[13.57]

In turn, the results from Equations (13.56) and (13.57) allow the calculation of the average derivative

Finally, this average derivative is used in Equation (13.54a) to obtain the corrected value of Fi+j. The same process is followed to find the corrected values of G, t L/, Ні+ij, and A’^i, using rearward differences in Equations (13.51) to (13.53) and calculating the average derivatives (9G/9§)ave, etc., in the same manner as Equation

(13.58) .

The above finite-difference procedure allows the step-by-step calculation of the flow field, marching downstream from some initial data line. In the flow given in Figure 13.8, the initial data line is the inlet, where properties are considered known. Although all the calculations are carried out in the transformed, computational plane, the flow-field results obtained at points (2, 1), (2, 2), etc., in the computational plane are the same values at points (2, 1), (2, 2), etc., in the physical plane.

There are other aspects of the finite-difference solution which have not been described above. For example, what values of At] and A§ in Equations (13.54a to d), (13.55), (13.56), and (13.57) are allowed in order to maintain numerical stability? How is the flow-tangency condition at the walls imposed on the finite-difference calculations? These are important matters, but we do not take the additional space to discuss them here. See chapter 11 of Reference 21 for details on these questions. Our purpose here has been to give you only a feeling for the nature of the finite-difference method.