Minimum radiative transfer geometries

Computers & Fluids, Vol. 1, pp. 379-398. Pergamon Press, 1973. Printed in Great Britain

MINIMUM

RADIATIVE

TRANSFER

GEOMETRIES*

H. J. GruELING1"and J. R. BARON:~ Department of Aeronautics and Astronautics, Massachusetts Institute of Technology, Cambridge, Massachusetts, U.S.A. (Received 13 December 1972) Abstract---Optimal control methods are applied to minimize the radiative energy transfer to a two-dimensional body with imposed isoperimetric constraints on base height and surface arclength. Primary emphasis is on the effect of changes in body geometry on the radiative energy transfer to the surface. The radiation field is specified by the differential approximation and simultaneous evaluation of the coupled flow and radiation fields is necessary. The fields are based on a one-strip integral method with a radiation closure condition applied at the base to circumvent the considerable difficulty associated with the attainment of radiative equilibrium in the large. A first-order gradient algorithm is applied to the system describing the state of the shock layer gas using body surface curvature as the control variable. A major problem proves to be numerical instabilities associated with the approximate, stiff system of ordinary differential equations. Some evidence indicates that the one-strip formulation leads to unrealistic behavior in the influence functions (adjoint variables) for the optimization procedure. Optimal solutions, obtained for a sequence of nose angles, demonstrate that nose angle has a crucial effect on the total radiative transfer to the surface for the assumed geometric constraints and flight conditions (Moo ~ 30). Physically, the reduced transfer is related to the shock layer volume control achieved with optimal shapes; however, average temperature also tends to decrease for the optimal sequence. An overall minimum radiative transfer geometry requires a cusped nose; nevertheless, the total energy transfer rates are less by an order of magnitude relative to a reference wedge level for nose angles as large as 25 per cent of the wedge angle. NOMENCLATURE Bo

Boltzmann number, (Puh/Q)~ dependent variable derivative, equation (20) G ~f/~Ko, equation (27) h enthalpy IA, Io, I~ angularly averaged intensity; surface and shock values J performance index, equation (12) body surface curvature, -dO/dx Ko L reference length M Mach number pressure P qx, qy radiative flux components radiative energy transfer rate atotal reference radiative energy transfer rate, Table 1 atotal(r©f) rB, R a body thickness, Fig. 1 ; body base height T temperature U, /3 velocity components V~t, Vac,cf) shock layer volume; reference volume, Table 1 X, X t arc-length coordinate; overall body arc-length dependent variable, equation (5) Y, Y~

f,f,

* Sponsored by U.S.A.F. Office of Scientific Research, Contract F44620-71-C-0009. t Present address: Los Alamos Scientific Laboratory of the University of California, Los Alamos, New Mexico. 379

380 ~Z

EB, ~w

0, 0~ ;~, AA, An P TL~ TX~ T6

H . J . GIBELING

and J.

R. BARON

volumetric absorption coefficient shock wave angle, Fig. 1 ratio of specific heats shock layer thickness, Fig. l base region emissivity,equation (l l); wall emissivity body surface angle, Fig. 1; nose angle influence functions, equation (24, 25) density reference optical depth, cqL; longitudinal depth; transverse depth terminal constraints, equation (22, 23)

Subscripts o0 W

0 F ¢ S

avg ()

free stream wall gas adjacent to wall shock wave designates flow field variable designates 'constraint' variable reference value for normalization averaged across shock layer consistent with integral method quantity normalized by wedge solution reference value 1. INTRODUCTION

Superorbital entry into planetary atmospheres can result in significant transfer of energy by thermal radiation, which may then dictate the gross features of vehicle geometry. Minimization of radiative energy transfer is then of some importance, and methods adapted from modern optimal control theory are considered here to investigate appropriate configurations to accomplish that goal. Minimum drag bodies based on variational techniques have received considerable attention in the past[l]. In the majority of such cases, simplifications afforded by Newtonian flow or linearization assumptions enabled the drag to be expressed as an integral over body length of a known algebraic function of body radius and slope. Unfortunately, a realistic expression for radiative heat transfer is not generally available, and instead necessarily involves simultaneous determination of the coupled flow and radiative fields in the shock layer. The problem of determining a minimum heat transfer geometry (radiation and/or conduction) has received relatively little attention. The geometry for minimum convective heat transfer was determined by Aihara[2] using the well known result of Lees[3]. Allen et al.[4] used a semi-empirical estimate for radiative heat transfer rates to show that cones (halfangle near 30 °) were more desirable than extremely blunt shapes when speeds become sufficiently high. The analysis neglected radiative cooling, geometrical complications, and selfabsorption so that the quantitative conclusions are somewhat questionable in detail. More recently, Tauber and Wakefield[5] showed that slightly blunted cones (half-angle near 60 °) minimize the total mass loss for entry into Jupiter's atmosphere. Neither Allen[4] nor Tauber and Wakefield[5] employed optimization procedures, but considered only conical shapes over a range of half-ar/gles. Furey[6] considered minimum energy transfer hypersonic nose and leading edge contours using a gradient method formulation and the Pontryagin maximum principle. For radiative heating, consideration was limited to moderate hypersonic Mach numbers and an optically thin gas; effectively, radiative energy losses were neglected in determining all shock layer properties except temperature. The latter distribution was then based upon an approximate energy equation including radiation cooling, and the flux to the body surface was specified by means of an optically thin tangent-slab approximation. The optimization procedure

Minimumradiativetransfer geometries

381

actually required that average values of temperature and velocity be taken at each body station, and the implicit assumption was then made that such averages were independent of local slope. Multidimensional radiative transfer is quite complex due to the integro-differential nature of energy conservation. However, it has been demonstrated that differential approximations give at least qualitatively reasonable one-dimensional results. Moreover, differential approximations have been employed with some success in multi-dimensional situations[7,8] despite uncertain accuracy. The differential approximation appears to be both useful and realistic, and is employed here as a first basis for radiative flux distributions. 2. FORMULATION For simplicity we assume the gas to be inviscid, nonconducting, and in local thermodynamic equilibrium. In the absence of transport properties optimum shapes are expected to be pointed, and only attached Rankine-Hugoniot shock waves are considered. A steady flow is assumed, as is a gray volumetric absorption coefficient to minimize numerical computation. The basic conservation equations and differential approximation (radiation moment) relations are readily available[9,10]. A practical optimization scheme virtually necessitates a reduction of the governing partial differential system to an approximate set of the ordinary kind. The speeds necessary for significant thermal radiation imply geometrically thin shock layers, so that the method of integral relations[l 1,12] is realistic for that purpose. The governing partial differential equations[13] are based in the body-oriented (x, y) Eulerian coordinate system illustrated in Fig. 1, with a one-strip integral method as a first approximation in view of its demonstrated

Lv~~"~qy Sh°ck

gOX Fig. 1. Geometryand coordinatesystem. success[12]. Reduction of the partial differential equations to the desired form involves the assumption of specific y-dependences for certain groups of flow and radiation field variables. Specifically the one-strip approximation assumes linear profiles for p u v , p u [ h + ½(u2 + v2)], u, qx and qy. A quadratic form is permissible for pressure, since the pressure gradient normal to the surface follows from y-momentum conservation. Similarly a quadratic profile may be constructed for the angularly averaged intensity IA. Using the differential approximation, the nonadiabatic contribution is - V • q = ct(IA - 4 a T 4)

C A F Vol. 1 N o . 4---E

(l)

382

H. J. GIBELINGand J. R. BARON

in which the Planck mean absorption coefficient for temperatures less than 20,000 ° K is reasonably represented by ~,-~pTS.[14] The temperature nonlinearity in equation (1) proves quite troublesome during the optimization procedure. Some caution is required when evaluating the average radiative flux divergence across the shock layer, since the high temperature behind the shock wave (T6) improperly dominates the energy balance. An approximate quadratic form constructed for T ( x , y ) consistent with the one-strip integral method approximation proves necessary and adequate[13]. For x # 0, the resulting conservation equations are of the form IYll

Ad dx

y~ =K(y)

@5 = (1 + Y5 Ko)tan(y4 - 0) dx d B~x

(2)

Y3 Y4

Y7 Ys

= L(y)

(3) (4)

Y9

where yr = (Uo, P0, Po, fl, 6, qx~, qy~, qxo, qro)

(5)

and q:,~ = q~6 - q~o, qr~, = qy, - qyo" The coefficients A i j , K i , Bij and L i are nonlinear functions of the dependent variables, and are summarized in the appendix. 3. R A D I A T I O N F I E L D C L O S U R E C O N D I T I O N

Consistent boundary conditions for radiation moment methods have been discussed by Finkleman[15] and Traugott[16]. For one-dimensional problems a Mark-type boundary condition from neutron transport theory[17] is consistent with full-range angular moments and has been employed here. Appropriate radiation boundary conditions at the surface and shock wave are, respectively, 2 - ew

Io(x) + - ~w

x/3qro(x) = 4aT~

I6(x) - ~ / 3 [qy,(x)cos(fl -7 0) - q~ sin(fl - 0)] = 4aT~

(6) (7)

with the upstream gas assumed to radiate as a black body at T~. Both conditions should remain valid as x ---,0, since the surface interaction is a local phenomenon. In addition, application of the method of integral relations at the apex and the requirement that the derivatives of the dependent variables be regular yields a relationship between surface and shock wave variables, and in a two-dimensional situation each field variable is identical at the apex for either shock or surface. It is clear from equations (6), (7) that one radiation variable must be guessed at x = 0 in order to establish the other two. A two-point boundary value problem results for the coupled flow and radiation fields, which is to be expected in view of the attempt of the differential approximation to simulate the global nature of the radiation field. In principle, one should

Minimum radiative transfer geometries

383

allow for eventual attainment of radiative equilibrium far downstream; however, this cannot be accomplished with complete consistency unless upstream absorption is included[18]. Moreover, for practical flight conditions and finite lengths, the transverse shock layer optical depth will rarely be large and evaluation must be continued beyond the base of the body. It would be extremely difficult, if not impossible, to obtain a solution in such a manner. However, at the base an expansion must be anticipated, and downstream temperatures decrease rather rapidly. This suggests that a precise closure boundary for radiation is of lesser importance. Since an approximate closure condition is not unique, an easily understood physical condition applicable at the base of the body was sought. Hence the ' closing boundary' is taken to be the normal to the surface at the base, and the radiation condition is applied at this boundary to account for radiant energy entering the layer from downstream. With one radiation variable guessed at x = 0, only a single condition remains to be specified at the base and follows from the Mark-type closure condition averaged across the shock layer in a manner consistent with the one-strip integral method. The result is ¼(IB - x/3q~a) = Q~ where IB = c~( i)

1 qxB --"

~(xs)

f~(xl)

I'*(xI ' y) dy

[~(~,) -o

qx(x I , y) dy

(8) (9) (10)

and Q~, the average one-sided radiant energy flux entering the shock layer normal to the closing boundary, depends on the temperature distribution downstream of the closing boundary as follows 1

f6(xs) a T g ( x l , Y) dy.

(11)

T h e ' base region emissivity' en accounts approximately for the radiation from the gas downstream of the closing boundary. 4. F O R M U L A T I O N OF T H E O P T I M I Z A T I O N P R O B L E M

The objective is to minimize the total rate of radiative energy transfer to a body subject to appropriate constraints. Define a scalar 'performance index' as J = -

q,o(x) d x = - Q,o,,~

(12)

which is the negative radiative heat transfer rate per unit depth for the two-dimensional body under consideration. The control variable is taken to be surface curvature, K o = - dO/dx, in view of its appearance in the differential state equations (2--4). It then becomes necessary to consider the slope, 0, as a state variable. In order to specify a meaningful optimization problem, two geometric or 'isoperimetric' constraints must be imposed. Here, the arc-length (xs), or surface area per unit depth, along the surface and the base height are specified. If the reference length, L, is taken to be the specified arc-length, the nondimensional geometrical constraints are x I = 1;

C

sin O(x) d x = R e , specified.

(13)

384

H . j . GIBELING a n d J. R. BARON

The integral base height constraint is most easily treated by introducing a new state variable, r s, such that dr B

sin O(x)

dx

(14)

with the boundary conditions r~(0)=0, rn(xf)=RB

(15)

which may now be treated as a terminal constraint in the context of optimal control theory [13,19]. In optimization theory the body nose angle, 0~ --- O(x = 0), may be either fixed or free; in contrast, a solution of the system equations requires specification of 01 just as Ko(x) must be specified. The essential point is that optimum geometries may be found for fixed nose angle in addition to seeking an overall minimum corresponding to no such constraint. The dependent variable vector must be enlarged to include both O(x) and rB(x). Symbolically, the new state vector may be written as yj

(16)

_-iy, ol ( r~° ) \Y11/

,17

Y= where YF is defined by equation (5), and

Yc

dyc (-Ko) ~xx = f c = sin0 "

(18)

The system of conservation equations (2-4) is easily reduced to the 'normal form,' dyF dx

-

fF(YF, 0, K0)

(19)

so that the complete system of state equations becomes

dy_ (fF) dx f[y(x), Ko] = fc "

(20)

The apex shock wave angle, fl(x = 0), follows from frozen shock wave relations; the initial conditions for equation (20) then follow from the frozen shock wave relations and equations (6), (7), and (15).

yj(x=O) =yj(00 yj(x = O) = 0 ys(x=O)=qx. y9(x=O)=y9(qx., 01) ylo(x=O)=01.

for j - 1, 2, 3, 4 for j = 5, 6, 7, 11 (21)

Minimum radiative transfer geometries

385

Reference states for normalization were taken to be those immediately downstream of the shock wave at the apex. An appropriate reference value for the radiation variables (q~, qy and IA) is unknown a priori, so that a general reference flux Qs has been chosen (note that the black body flux normalization, Qs = aTe, may be inappropriate in problems of this kind). With this normalization procedure it is necessary to consider the dimensional initial conditions equation (21) in order to carry out the optimization scheme. Terminal constraints for optimization purposes consist of the base height specification (15) and the radiation field closure condition (8). It is convenient to define a terminal constraint vector

I//A

(22)

where ~A(x s) = r,,(x s) -- Rx, 1

~B(xl) = q~, + "-'75 (4Q~ - In).

(23a) (23b)

The functional forms of In, qxB, and Q~ follow from equations (9-11) upon substitution of the appropriate interpolation polynomials for l , ( x , y), q,,(x, y) and T4(x, y). Of course, ~,a(xf) and ~O~(xs) may not vanish in general, since at an arbitrary step in the optimization procedure the terminal conditions (23) may not be satisfied precisely. The performance index J, equation (12), is to be minimized subject to the system equations (20) with the associated initial (21) and terminal (23) conditions. It proves possible to satisfy equation (23b) prior to the starting the optimization procedure by simply solving the twopoint boundary value problem for qx.; the geometric constraint equation (23a) can be satisfied by the proper choice of a ' nominal' curvature distribution Ko(x) for a given nose angle 01 . It should be realized that only the optimal' solution' is required to satisfy all of the terminal conditions for either fixed or free nose angle. Nevertheless, it is preferable to start with a ' nominal solution' which satisfies such terminal conditions. Here the term ' solution' refers to any consistent combination of Ko(x), 01 and q,., and the resulting uniquely determined state variable distribution y(x), 0 < x <_ x f . A first-order gradient algorithm has been employed to determine the desired optimal solutions; the specific procedure is outlined in [13], while a somewhat different algorithm is presented in ref. 19. Briefly, several sets of linear influence functions or adjoint variables are constructed in order to predict the effect of small changes in the dependent (state) variables on the performance index J and the terminal constraints q(xf); i.e., ~.(x) is an influence vector for changes in J, equation (12), due to changes in y(x), 6J 2i(x)~ayi(x).

(24)

Similarly, Aa(x) and Aa(x) are influence vectors for changes in g'A and g'B,

&,(x) ~ aO,(x,) 6yi(x) ' ct = A or B.

(25)

386

n. J. GIBELINGand J. R. BARON

Such influence functions permit construction of an 'impulse response function' 8H/dKo representing the variation in J due to a unit impulse in K o at the position x for fixed initial conditions y(x = 0) while satisfying the system equations (20), but not necessarily the terminal constraints (23). Here H is the Hamiltonian,

H(x) = -qyo(x) + ~rf

=

-Y9

+ ~rf

(26)

and ~K----~

~ o --=--~'rG'

(27)

Further, A~G and Anr G are impulse response functions representing the variations in Oa(X:) and On(X:), respectively, due to a unit impulse in Ko(x ). The first-order expressions for the changes in J, 0a(X:) and OB(x:) due to weak variations in the dependent variables y(x) and the control variable Ko(x ) are simplified considerably by the introduction of ~., AA and AB [13]. It can be shown that ~ satisfies a linear nonhomogeneous adjoint differential equation, d2i dx

-

-

Fji2j + ~i9

(28)

where Fij = ~fffc~yj may be termed the influence matrix for the system (20) and c5i9 is a Kronecker delta. AA and An satisfy the homogeneous adjoint equation dA~, _ dx

FjiAaj , :~ = A or B.

(29)

Boundary conditions at xf for ~,, AA and AB (transversality relations) follow from the variational formulation. The special form of the state equations, terminal constraints and boundary conditions permits certain simplifications in the gradient formulation[13]. Actually, AA(X) may also be determined quite simply from geometrical considerations, and F o is singular at x = 0 due to the vanishing shock layer thickness. On physical grounds the influence functions themselves must remain finite, so that in practice a polynomial extrapolation was used to obtain values for ~ and AB at the apex. Preferably the optimization procedure should maximize the predicted change IdJI in the performance index at any given stage. Since a first-order gradient approach is being employed, a quadratic penalty[19] is imposed on the size of the optimization step which may be taken, thereby yielding the optimal corrections 8Ko(x), 6qx, and 601 [13]. 5. DISCUSSION OF NUMERICAL PROCEDURES The optimization algorithm was not completed without encountering numerical difficulties. Of major importance were numerical instabilities associated with the system equations, and there is some evidence that the one-strip integral method formulation can lead to unrealistic behavior in the influence functions if care is not taken. These topics and the related question of optimization step size control merit discussion. It has been observed[15,20] that radiation moment equations imply unstable roots that may cause divergence when integrating toward radiative equilibrium. In one-dimensional problems, this occurs when the optical depth measured from the boundaries becomes large. Finkleman[15] points out that the difference formulas of Curtiss and Hirshfelder[21] force radiative equilibrium in the far field, which may not be suitable for self-consistent problems.

Minimum radiative transfer geometries

387

As mentioned earlier, the attainment of radiative equilibrium is unlikely for the geometries considered here. More important, however, is that the omission of upstream absorption and application of a one-strip integral method may in fact preclude the attainment of equilibrium, The stability of a system of ordinary differential equations is discussed in some detail by Lomax[22]. Consider the system dy = Cy + b dx

(30)

where C is a constant matrix. Only the eigenvalues of C are important in determining the stability and accuracy of any numerical method, and the system is ' inherently unstable' if the real part of any eigenvalue is positive. For the nonlinear, nonconstant coefficient system (20), the eigenvalues of the influence matrix Fij = ~fi/ay~ must be considered. The numerical stability of an inherently unstable system is virtually never discussed except in the context of relative stability (i.e., exponential growth of the solution due to the prescribed initial or boundary conditions). Unfortunately, the present system of differential equations [specifically the subset of radiation moments, equation (4)] is inherently unstable. The existence of exponentially growing solutions causes solutions based on conventional difference formulas to diverge, despite initial conditions that imply no such behavior. The influence matrix of the system possesses positive eigenvalues proportional to [6(x)]-1 ; near the apex the instability is quite severe since 6 ~ 0. Clearly a physically reasonable solution does not exhibit exponential growth and the difficulty appears to be associated with application of the integral method to the differential approximation in the presence of a pointed shock layer. Cheng and Vincenti[8], e.g., applied a series truncation method tO the differential approximation for a stagnation region flow and mentioned no stability problems. An integral method also aggravates the instabilities due to inherent differencing of nearly equal quantities near the apex. A Taylor series solution in the nose region serves to alleviate this otherwise troublesome loss in accuracy. A simple but illustrative example of the numerical problem and its solution was presented by Curtiss and Hirshfelder[21 ] for a single first-order differential equation of the form

dy dx

--

=

y(x) - Z(x) c(x, y)

(31)

They indicate t h e ' stiffness' corresponds to [c(x, y)/Ax I ,~ 1, where Ax is the desired integration step size. Only for a special solution y = Y(x), close to Z(x), does the solution not grow exponentially regardless of the initial conditions. For a system of equations, the situation described by Lomax[22] holds true. Curtiss and Hirshfelder[21] show that the implicit forward difference scheme (first-order in Ax) Y~--Yn-t__ Y ~ - Z n

Ax

c(x,, y,)

(32)

converges to the 'correct' solution independent of the initial conditions, but involves iteration when c depends on y. Otherwise, c ~ c(y), an interesting asymptotic solution can be obtained[21]. If there is inherent instability, i.e., c(x) > 0, it can be shown that a minimum step size exists below which the solution will be unstable. The implication is that for sufficiently large gradients, the numerical method must fail due to unacceptably large truncation error associated with the minimum step size.

388

H.J. GIBELINGand J. R. BARON

The present system may be classified as ' stiff' because of the wide variation in magnitude of the influence matrix eigenvalues[21,23]. However, the implicit forward finite difference formulas[21] do control the spurious numerical errors due to unwanted exponentially growing solutions, and can often provide stable, accurate solutions. Gear[23] has further developed the higher-order, multistep versions of this stiff approach, and has suggested a special form that allows both variable integration step size and variable order. He also showed that use of a Newton iteration process to solve the basic corrector formula leads to convergence for larger step sizes than with standard corrector iteration procedures. The stiff nature of the radiation moment equations for pointed shock layers along with the problem of numerical instability suggested use of the Gear approach. The optimization iteration procedure is initiated by assuming a nominal curvature distribution Ko(X) such that both terminal constraints (23) are satisfied. Using a Taylor series starting solution[13] the state equations are integrated from x = 0 to xf, and the entire state variable solution stored for that range. The total heat transfer rate follows from equation (14). After the first optimization step the actual change in J, dJ a, is determined, and adjustment of the step size is made by comparison of dJa with the change, dJ~, specified prior to taking the step. The adjoint differential equations (28) and (29) are then integrated backward from x = x I to 0, and the impulse response functions GT~., GrA~ and GTAA stored for 0 _< x_< x I . The optimal corrections 6Ko(x), 13qx,~and fi01 are specified by the first-order variational formulation[13]. Briefly, fiqxn and ~01 are expressed as linear combinations of the influence functions at x = 0; 13Ko(x) is expressed as a linear combination of the impulse response functions. The coefficients in these expressions are Lagrange multipliers determined such that both a specified change in the terminal constraints (23) and the specified step size parameter, dJs, are satisfied to first order. The quantities K0(x), q~ and 01 can then be modified by the optimal corrections, and the sequence repeated. On physical grounds the individual adjoint variables 21(x) and ABe(x) would be expected to exhibit smooth variations. Figure 2 (Case 3, Table 1), however, indicates the influence function, AB4, for variations in shock angle can have a large amplitude oscillatory behavior that decays as integration proceeds toward the nose, and is not due to numerical instability. In fact such behavior may be traced to the assumed temperature profile. A quadratic profile for T(x, y) shows appreciable improvement (Fig. 2) over linear profiles for ~(x, y)[ = pT 5] and T4(x, y) used initially. Physically, the source of the difficulty is a poor representation of the radiative cooling term ~T~[ = pT 9] in the integral method, since only surface and shock wave values are adjustable. The result is an unreasonably strong dependence of (~T4)avg on shock wave angle, and is compounded by the highly nonlinear relationship between (~T4)avg and fl(x). Further evidence is the relatively small oscillations found in ).~(x), which is related to J being an integrated quantity and therefore insensitive to local changes in shock angle. Simultaneously, the resultant oscillations introduced in Ko(x) remain unacceptably large, and eventually induce divergence of the optimization iteration. The cost of implementing a multi-strip integral method would be excessive and probably not eliminate the entire problem. A practical approach is to approximate the impulse response functions GTAB and Grk by low order least squares polynomials such that the mean behavior of the oscillating quantities is well represented[13]. The values of the influence functions themselves are required only at the apex, and appear to be sufficiently accurate for present purposes. It is essential to realize that the calculated geometries based upon the smoothed gradient information (i.e., GrAB and Grk) always yield a valid flow and radiation field (with iteration for the appropriate q.~, if necessary), independent of the manner in which the geometry

Minimum radiative transfer geometries

7°f

]

l

I

6.0 5.0

-I

389

--

-/

Basis: a(y) ond T4(y) linear

--~

4.°

.:

Bos,. : .,.,

\ \ A

.oo°,o,,0

i

A

00~ . . . . I0

0 Arc- length x

Fig. 2. Transverse profile effects on shock angle influence function, case 3. may have been obtained. Without a smoothing procedure for GrAB the linear prediction for 6qx, proved to be quite accurate, but led to divergence during the optimization iteration ; with the smoothing procedure convergence results but with some associated error in the predicted value of 6qx,. In order to satisfy the closure condition (23b) using the predicted 6q,,. the necessary optimization step size is unduly small, and hence the correct value of q~, is best determined by iteration when the closure condition error ffn(xs) is greater than a specified value[13]. This procedure is necessary only if' solutions' during optimization iteration are simultaneously to yield valid flow and radiation fields as well. The specified change in the performance index dJs and the actual change dJ a after a given optimization step remain in good agreement even with the smoothing procedure in effect, as long as the linear relations for 6Ko(x ), 6qx, and 60, remain valid. Of course, in approaching the optimum solution, the step size parameter dJ~ must decrease to prevent divergence, since a linear step along the local gradient in state space is employed. However, only the total radiation heat transfer rate to the surface being of physical interest, the optimization procedure was terminated when the actual change in J was less than 0.5 per cent. 6. RESULTS AND CONCLUSIONS The specific example of radiating wedge flow furnishes a realistic reference flow field. All results must necessarily be viewed as somewhat qualitative in any case, since real gas effects will be important for the range of conditions considered and perfect gas shock wave relations are employed in the present work to facilitate the computation of the partial derivatives with respect to shock angle necessary for both the integral method and the optimization procedure.

390

H . J, G I B E L I N G

and

J. R. B A R O N

In partial compensation for perfect gas (? = 1.4) simplicity, flight conditions have been chosen to achieve appropriate temperature levels. The free stream properties (p~, T~) for the cases presented correspond to an altitude of 60 kin. In particular, a wedge of half-angle 0~ = 0.61 radians at Mach number Moo = 28.5 has an implied shock layer temperature of 19950 ° K at the nose. To achieve this temperature in air with equilibrium chemistry, an actual free stream velocity of about 34 km/sec is necessary[24]. The effective emissivity eB [equation (11)] can be shown to be quite small in practical situations, and the value eB = 0 was employed in all cases considered with an approximate error in total heat transfer rate of at most a few per cent. Both the shock layer transverse optical depth %(x) and an average longitudinal optical depth r~(x) are relevant for the shock layer: 6

%(x) = "cL/o ~(x, y) dy

L,(x)

'= r4x') = zL fo ~ - ~ dx'

and are presented in Fig. 3. The shock layer proves to be optically thin in both transverse and longitudinal directions, despite an appreciable reference optical depth (rL = 0'62). 0.07[

1

V~-

I

J

o,oe ~ [

0.05

i

I

i

'1

0.04~

r 0.03 ~-

i ,i -c8 x I 0

0.02 -~

Ji I

i

o.ofl

oE_ 0

L

02

l

0.4. Arc-lenqth

l

0,6

. . . . .

0.8

i

1.0

x

Fig. 3. Longitudinal and transverse optical depth distributions, case 1 ( w e d g e f l o w , 01 = 0.61 r a d . , Moo ~ 28"5, Ts = 1 9 9 5 0 ° K ) .

Figure 4 shows wedge flow values of Uo, u~, Po, P~, To, T~, Po and 5. Since the wall temperature is relatively low (1000 ° K) the radiation field is emission-dominated, i.e., Ia ~ 4aT4~ Q~ as expected, and the variations are consistent with a strongly radiating flow. In particular the temperatures To, T6 and density Po show the effects of significant radiative cooling.

Minimum radiative transfer geometries

I.a~,

I

'

]

'os

I

I. 0

'

I

'

L

--3.0

Uo

~ 0,8

391

,../

.'"

-

T.,PSJ

- 2.5

0,6 0,4 0.2

_ _ / /

0

--

0,2

0.4

0.6

Arc-

leflqth

0.8

1.5

1.0

x

Fig. 4. Gas property distributions along surface and shock, case 1 (wedge flow, 01 = 0'61 rad., M® = 28.5, T, = 19950 ° K).

The components of the radiative flux at both the surface and shock, qxo, qx~, qyo and qy~, are shown in Fig. 5. Also shown is the a posteriori calculation of t.he surface normal flux component using a local tangent slab approximation and the calculated temperature profile consistent with the integral method and differential approximation. Since IA ~ 4GT4/Q~ the 0.81"0 ~

0,6

qY8x I0~

•~0--}

Differential al~ocoxlmation

A postltiOti" IO(:OItangent

0.4

slab

0.2

~x° x I0 z - -

0

~xax I0 z _

-02 -0.4 -0.6 -0.8 -1.0 -t.2

~yox aO 3 I

I

02

0.4

I 0.6 Arc - length

I

o,a

I

l.o

•

Fig. 5. Surface and shock wave radiative flux component distributions, case 1 (wedge flow, 01 = 0.61 rad., Moo = 28"5, T~ = 19950 ° K).

392

rl. J. GIBELINGand J. R. BARON

temperature distribution is not strongly dependent on I o and I~; thus the tangent slab comparison indicates the consistency of the differential approximation for qro(X). The remaining discrepancy is deemed acceptable in view of both the inexactness of the tangent slab result and the differential approximation used in an optically thin situation. The surface longitudinal flux component qxo is somewhat large compared to qx~, which may be attributed to use of the integral method. Due to averaging the flux divergence across the layer, dqxo/dX is sensitive to errors in ~[Ia - 4aT4/Qs] as a result of uncertainty in ~T 4 mentioned previously. Even though Iq.... ] is roughly five times larger than either Iqrol or lqr6], there is little longitudinal energy transport[13]. Errors in dqxo/dx do not have a large influence on other flow and radiation variable derivatives. Comparison of qx,,,, and I, vg (--- x/3qy,) reveals that the radiation field is quite anisotropic over most of the shock layer as was to be anticipated. The benefits of optimization are illustrated in Fig. 6 for two Mach numbers. The total 20

'

-'1

-r

]

r

I

Constant K n

(Moo= 32)~'~//p

1,6=. --

~

'

/(3

/ /

--

/

Optimol ,~.....~//

(M~ = 3z)

a'

1.2 o°

-¢

/

-

/

Constant

d'" ~ / 7 /,*

/

°'7 / ol0

~3

/,o [ 0,2

~

! I 0.4

Nose Angle

I 0.6

01

,

I 0.8

(rod)

Fig. 6. T o t a l surface radiative energy transfer rate vs. nose angle f o r optimal and constant

curvature geometries. radiative energy transfer rate to a surface is displayed for a range of nose angles corresponding to specific optimal curvature distributions as in Fig. 7. For reference there are also shown constant curvature body results for the identical geometric constraints (Cases 1-5). A tabulation of key values for cases 1-14 appears in Table 1. The reduction in energy transfer is quite appreciable. The shape distributions corresponding to the actual geometries (Fig. 7) indicate relatively little change in curvature over the last half of the body, and reflect the importance of the large (negative) initial curvature as the nose angle decreases[13]. In nearly all examples, a substantial portion of the decrease in I Qtot,l[ is due to a reduction in shock layer volume. In effect, the strong compression buildup due to the initial surface curvature causes the shock wave strength to increase rapidly (Fig. 8), but over a length sufficient to permit appreciable radiative losses from the

T a b l e 1. S u m m a r y

Case No.

01 (radians)

Q,ota, Qtotal(ref)

1" 2* 3* 4* 5* 6 7 8 9 10 11 12 13 14

0"610 0"560 0"500 0"434 0"350 0"610 0"573 0"500 0"400 0.269 0"149 0"100 0"050 O'OlO

1"0 0-905 0-837 0"820 0"849 0"922 0"815 0"641 0"477 0.277 0.145 0"109 0"083 0"068

Vs~ Vs/(reO 1"0 ----0"991 0"905 0"810 -0"444 0"230 0.151 -0-0875

TL 0-620 0"287 0-101 0-275 x 0"377 x 0"620 0"351 0"101 0"130 x 0"328 x 0"189 x 0"918 x 0-214 x 0-838 ×

10 -1 10 - z

10 -1 10 - 3 10 - s 10 - 7 10 - s 10 - l °

rL B--"O 23"1 6"35 1"14 1"35 x 5"42 x 23"1 8"91 1"14 3-98 x 1"09 x 3"60x 2-89 x 1"43 × 1-43 x

10 -1 10 - 3

10 - 2 10 - 4 10 - 6 10 - 7 10 - s 10 - 8

of results, cases 1-14

Tx(x I ) X 10 0'66 ----0'64 0"61 0"56 -0"42 0-37 0"38 -0"41

To(xl) X 10 z 0"29

T,

--

1-0 0-856 0"875

--0"40 0"43 0"34 -0'079 0-043 0"043 -0'027

0"948 1"09 1"0 0"914 0"905 0"875 0.887 0.845 0"830 0"825 0"881

-

-

(m -1 )

T0(max)

0.355 0-164 0.580 x 0-157 x 0.216 x 0"355 0.201 0.580 x 0.742 x 0.188 x 0-108 x 0.526 x 0.123 x 0"480 x

M = 28"5, ~ = 1"4, R n = 1"0 m , L = 1"7456 m , T~, = 1 0 0 0 ° K . , ew = 1 '0, Qtotaltref) = - 0 . 1 4 1 4 x l 0 s W , Vs.,ef) = 0 . 1 2 7 m 3 * C o n s t a n t c u r v a t u r e c a s e n u m b e r s ( 1 , 2 , 3 , 4 , 5 ) c o r r e s p o n d t o Ko = ( 0 , - - 0 - 1 0 0 , - - 0 " 2 2 3 , - - 0 " 3 6 0 , - - 0 " 5 3 7 )

10 -1 10 - 1 10 - 2

10-1 10 - 2 10 - 2 10 - s 10 - 7 10 - 8 10 - 1 °

(°K.) 19950 17110 13900 10720 7220 19950 17810 13900 9230 4450 1620 913 467 286

as

( W / m 2) 0.898 x 0.485 × 0.212 x 0.749 x 0.154 x 0.898 x 0.571 x 0-212 x 0.412 x 0.223 x 0-435 x 0-350 x 0.253 x 0.180 x

101° 101° 101° 109 109 101° 10 l ° 101° 109 10 s 108 10 s 10 s 109

394

H.J. GIBELING and J. R. BARON 06~--

3

]

7---i

I ---~

I~---1---T Case 8

9

1 --FII 12-

SJJ../

, ~°~ i

/yjj

03

~ i

~

0.2

~,~"

OI

/~ 0

0

~

/.. / 1/

0

// O.I

~J 0.2

0.3

0.4

0.5

0.6

0.7

0.8

Axiol Length xox

Fig. 7. Optimum two-dimensional geometries.

~0~

I

I

[

°~~.~

I

0.41

i

O. 02~--/

--

I I

°'1~ II

o? 0

~ J

'

0.2

I

_. . . .

0.4

L___~ 0.8

06

Arc - lenglh

i

i.O

x

Fig. 8. Shock wave temperature distributions.

layer. The resulting average temperatures across the layer are decreased markedly, with corresponding higher densities reflected in thinner layers (Fig. 9). The large initial curvature also aids in satisfying the geometric base height constraint as the nose angle decreases. A comparison of the surface normal flux distribution with the a posteriori local tangent slab result (Fig. 10) shows the slab model to be extremely sensitive to the local temperature distribution (Fig. 8), in contrast to the relatively insensitive differential approximation. Near the nose the slab model underpredicts Jqyo ] in view of the one-sided dominance of the hot downstream regions. Similarly, near the maximum shock wave temperature station, the

395

Minimum radiative transfer geometries 0.08

I

I

I

0.07 - -

i I 0.06 --

i 0.05 - -

to

0.04-

0.03

Cose

--i

--

0.02

O.01 14

o ~-------~i--~ 0

0.2

I

]

I

0.4

0.6

0.8

1.0

Arc - l e n g t h x

Fig. 9. Shock layer thickness distributions.

-0.;

- O.c

X

o

-0.6

+

S':S'r,°S: tangent slab

-

-0.8

-I.O

-I.2

0

I

I

I

0.2

0.4

0.6

I 0.8

1.0

Fig. 10. Surface transverse radiative flux component distributions. flux is over predicted due to the smaller radiative flux that actually emanates from the cooler adjoining regions. Nevertheless, a comparison of the heat transfer rate Qtotal using either basis (Fig. 11) indicates close agreement for the total amount of energy radiated from shock layer to surface. It is interesting to compare these results with Furey[6] who on a tangent

396

H.J. GIBELINGand J. R. BARON

I4

l

I

I

[

I

1 -~

1,2~-

}

!

0.8 ~~" ~

o ~'4~

0

A posteriori locol ~-On ~ b ~"~.

01

~

0.2

i

0.3

~_

0.4

,., ///fie ~

~

05

, -~ |

,

0.6

j

0.7

Nose Angle 01 (rod)

Fig. 11. Total surface radiative energy transfer rate versus nose angle (M~ = 28.5).

slab basis concluded o p t i m u m axisymmetric bodies were conical but with cusped noses of extremely large curvature. In summary, large reductions in radiative energy transfer to a surface m a y be realizable with appropriate changes in nose angle a n d surface shape. Sharp nose regions followed by d o w n s t r e a m expansion a n d recompression as required by geometric constraints are the qualitative features associated with such configurations. Physically, the reduced energy transfer experienced by the surface is attributed largely to control o f shock layer volume.

REFERENCES 1. A. Miele, (ed.), Theory of Optimum Aerodynamic Shapes, Academic Press, New York (1965). 2. Y. Aihara, Optimum body geometries of minimum heat transfer at hypersonic speeds, AIAA J. 6, 21872189 (1968). 3. L. Lees, Laminar heat transfer over brunt-nosed bodies at hypersonic flight speeds, Jet Propul. 26, 259269 (1956). 4. H. Allen, A. Seiff, and W. Winovich, Aerodynamic heating of conical entry vehicles at speeds in excess of earth parabolic speed, NASA TR R-185 (1963). 5. M. E. Tauber and R. M. Wakefield, Heating environment and protection during Jupiter entry, AIAA Paper No. 70-1324 (1970). 6. R. J. Furey, Minimum energy hypersonic nose and leading edge shapes, AIAA Paper No. 70-825 (1970). 7. P. Cheng, Dynamics of a radiating gas with application to flow over a wavy wall, AIAA J. 4, 238-245 (1966). 8. P. Cheng and W. G. Vincenti, Inviscid radiating flow over a blunt body, J. Fluid Mech. 27, 625-646 (1967). 9. W. G. Vincenti and C. H. Kruger, Jr., Introduction to Physical Gas Dynamics, Wiley, New York (1965). 10. S. S. Penner and D. B. Olfe, Radiation and Reentry, Academic Press, New York (1968). 11. A. A. Dorodnitsyn, A contribution to the solution of mixed problems of transonic aerodynamics, Advances in Aeronautical Sciences, Vol. 2, Pergamon Press, Oxford (1959). 12. J. C. South, Jr., Application of the method of integral relations to supersonic nonequilibriumflow past wedges and cones, NASA TR R-205 0964). 13. H. J. Gibeling, Optimum radiative transfer geometries in hypersonic flow, Ph.D. Thesis, M.I.T. (June 1972).

Minimum radiative transfer geometries

397

14. S. C. Traugott, Shock structure in a radiating, heat conducting and viscous gas, Martin Company Report No. RR-57 (1964). 15. D. Finkleman, Self-consistent solution of nonlinear problems in unsteady radiation gasdynamics, M.I.T. Aerophysics Laboratory, TR 147 (1967). 16. S. C. Traugott, An improved differential approximation for radiative transfer with spherical symmetry, AIAA J. 7, 1825-1832 (1969). 17. B. Davison, Neutron Transport Theory, Oxford University Press, London (1958). 18. D. B. Olfe, A modification of the differential approximation for radiative transfer, AIAA J. 5, 638-643 (1967). 19. A. E. Bryson, Jr. and Y. C. Ho, Applied Optimal Control, Blaisdell, Waltham, Mass. (1969). 20. M.P. Sherman, Moment methods in radiative transfer problems, J. Quant. Spectr. Radiative Trans. 7, 89-109 (1967). 21. C. F. Curtiss and J. O. Hirschfelder, Integration of stiff equations, Proc. Nat. Acad. ofSci., Vol. 38 (1952). 22. H. Lomax, An operational unification of finite difference methods for the numerical integration of ordinary differential equations, NASA TR R-262 (1967). 23. C. W. Gear, The automatic integration of stiff ordinary differential equations, Information Processing 68, Proc. IFIP Congr. 68, Vol. 1, Morrell, A. J. H. (ed.), North-Holland Publishing Co., Amsterdam (1969). 24. P. V. Marrone, Normal shock waves in air: equilibrium composition and flow parameters for velocities from 26,000 to 50,000 ft/sec., Cornell Aeronautical Laboratory Report No. AG-1729-A-2 (1962).

7. A P P E N D I X The coefficients A lj, Kl, B~ and L~ are presented in very compact form; further details may be found in

[131. A 1 3 = A 2 2 = A 2 4 = A 3 1 ~ A 3 2 = A 3 3 ~ 0; A 2 3 = I ; A a a = po

A12 = A,*3= Uo; A,4 = ~fl(pu)~; A2, : yM~(OU)o (ho + 3

A.

=

A42

+

Bl2 = B14 = B31 = B32 = B33 = B41 ~ B43 = 0 B13 Blx =B21 =B2a=B34 = --f : 1 B22 = B24 = --cot(fl -- 0); B42 = S2(x); B44 = Sl(x) KI=-(~

1 d~

+ Ko)(pv), + ~-~x[(PU),--(pu)o]

K2=0

Id8

Ko)

1 d3

1

(pu)o [ho

CAF Vol. I No. 4--F

.

398

H. J. GIBELING a n d s. a. BARON

d3x _ 2q~,a) _ 2(Ko qya +'rtPR) LI = g 1~ [q xa~x [, 3rL(l +

Ko ~) ~,-- dfl Ko] ~x --

L2 = [qxn cot(fl -- 0) + qy,] IV/3 cos (fl -- O)

[Sd~o + ~ O d ~ ) + 2 R 2 ~ + 3R,--(1 + Ko 8)~,[qxn +qy, tan(fl--O,]) La = T

X/~ (2--ew t

3rL

~

x-TS-~ / + " T ~o o

it

d~

~__2__(2- ,. 1

L" = 3 tRe -~x + ~/~r L ~ e.

/ L3 -- (l + Ko S)~dq~, + q,, tan(fl-- O)]} -- Dlq, o -- D2 qv,

where qxa

qxn

=

- -

qxo; qya = qr6

- -

q~o

dSk Dk= dx RI='~

(l+yKo)o:qxdy=Slqxo+ S2qxa+goS(S2qxo+Saqx~)

1

R2 = ~ fo otqydy = $1 qm + $2 qyA ea=g

o~(l + y K o ) E d y =

k=OEk[Sk+t +(KoS)Sk+2]

E----Ia-- \ Q~]4h 4=k=o ~ Ek(x)

8

h=ho+,--3ho+4hl--h,,; + 2(ho-2hl + h,,(~) 2= [,~=~o't(~)t~]TM

[

]

t = p u h + - ' - ~ M~(uZ + v 2) ; # = p u ; s=puv.

Recall that linear profiles are assumed for t, #, s and u.

1 y

s,=~jo=~

lye*-1

)

)

-

s dy: =_Yolk" = tj+~ + ~ ] ,

k->,.