1、 - NASA Technical Paper 1619 Application of a - LOAN COPY: RETURN TO NASA . i AWL TECHNICAL LIBRARV Tp KIRTLAND AFB, N.M. 161 9 ! c.1 I ! Generated Numerically Orthogonal Coordinate System to the Solution- of Inviscid -Axisym.metric . Supersonic Flow. Over Blunt Bodies , “ ., H. Harris Hamilton I1 a
2、nd Randolph: A. Graves, Jr. d MARCH 1980 MSA , ,. Provided by IHSNot for ResaleNo reproduction or networking permitted without license from IHS-,-,-TECH LIBRARY KAFB, NM NASA Technical Paper 1619 Application of a Numerically Generated Orthogonal Coordinate System to the Solution of Inviscid Axisymme
3、tric Supersonic Flow Over Blunt Bodies H. Harris Hamilton I1 and Randolph A. Graves, Jr. Latzgley Research Cerrter Harnptorz, Virgitlia National Aeronautics and Space Administration Scientific and Technical Information Office 1980 Provided by IHSNot for ResaleNo reproduction or networking permitted
4、without license from IHS-,-,-SUMMARY A numerically generated orthogonal coordinate system (with the body surface and shock wave as opposite boundaries) has been applied with a time asymptotic method to obtain steady-flow solutions for axisymmetric inviscid flow over several blunt bodies including sp
5、heres, paraboloids, ellipsoids, hyperboloids, hemisphere cylinders, spherically blunted cones, and a body with a concavity in the stagnation region. Canparisons with experimental data and the results of other computational methods have demonstrated that accurate solutions can be obtained with this a
6、pproach. The numerically generated orthogonal coordinate system used in the present paper should prove useful for applications to complex body shapes, particularly those with concave regions. In addition, the use of the present orthogonal coordinate system simplifies the form of the governing equati
7、ons and simplifies the application of boundary conditions at the body sur- face and shock wave. INTRODUCTION Previous investigations of the direct blunt-body problem have demonstrated the utility of treating the bow shock wave as a discrete discontinuity. (See, for example, refs. 1 to 5.) In each of
8、 these studies, some form of coordinate transformation was used to map the physical space within the shock layer into a rectangular computational domain, with the bow shock wave as one boundary and the body surface as the opposite boundary. These transformations generally pro- duce a nonorthogonal c
9、oordinate system in the physical space which complicates the form of the equations of motion. However, they usually produce a uniform grid system for calculating the flow-field solution and simplify the application of the boundary conditions at the shock wave and the body surface. Recently, Graves (
10、ref. 6) applied a method for numerically generating an orthogonal coordinate system between two arbitrary continuous curves. By allow- ing one of the curves to be the body surface and the other to be the bow shock wave (which can move with time), this method can be used to construct an orthog- onal
11、coordinate for solving the blunt-body problem. An orthogonal coordinate system of this type leads to a much simpler form of the equations of motion and further simplifies the application of boundary conditions at the body sur- face and the shock wave. In addition, this type of coordinate system can
12、be applied to complex body shapes, including bodies with concave regions. The purpose of the present paper is to explore the application of the numerically generated orthogonal coordinate system presented in reference 6 to flaw-field problems. For this purpose, the inviscid, axisymmetric supersonic
13、flow over a blunt-nose body will be investigated. Provided by IHSNot for ResaleNo reproduction or networking permitted without license from IHS-,-,-SYMBOLS speed of sound, ;/Vm ; also length of Z-axis for ellipse -2 a Bb b cP CV d h3 I i L A M P Rb R s t U VS Vn Vt bluntness parameter length of r-ax
14、is for ellipse specific heat at constant pressure specific heat at constant volume diameter of cylinder, i/i metric coefficients in 6- and rl-directions, ;1/i and hz/i, respectively metric coefficient in circumferential direction, i3/i static enthalpy, ;/ern index for Steady-flow solutions are obtai
15、ned by integrating the system of equa- tions (2) to (6) in time, from an assumed set of initial conditions, until steady state is reached. Since the integrated form of the energy equation (eq. (5) is being used, the results are not accurate during the transient phase of the solution. However, equati
16、on (5) becomes more accurate as the 5 Provided by IHSNot for ResaleNo reproduction or networking permitted without license from IHS-,-,-steady state is approached and the correct converged steady-flow solution is obtained Using the energy equation in this form speeds up the computational procedure.
17、Along the stagnation streamlines u = h3 = 0 and ah,/ag = hl; thus the u ah3 V ah3 h2h3 an terms - - and - - appearing in equation (2) are of the indefinite form (O/O). Taking the limit of these terms as 6 -+ 0 by using LH6pitals rule h1h3 35 Thus on the stagnation line (5 = 0), the limiting form of
18、the continuity equa- tion (2) becomes 6 Provided by IHSNot for ResaleNo reproduction or networking permitted without license from IHS-,-,-Equations (2) to (9) are written in nondimensional form using the following definition for the nondimensional variables: Finite-Difference Scheme The solution of
19、the flow field is obtained by integrating the compressible, time-dependent Euler equations described in the previous section. The integra- tion is carried out using the explicit finite-difference scheme of Brailovskaya (ref. 7) . Using the following computational module i- 1 i i+l T“-t“-t j+l I I I
20、I I I . this two-step difference scheme can be represented as follows: 7 Provided by IHSNot for ResaleNo reproduction or networking permitted without license from IHS-,-,-Predictor step: Corrector step: t+At Fi,j where F can represent p, u, or v. The time derivatives appearing in equa- tions (1 l a)
21、 and (1 l b) are obtained from equations (2) to (4) . The spatial deriv- atives in equations (2) to (4) are evaluated using central difference formulas except at the boundaries where three-point, one-sided difference formulas are used. The numerical stability limit for this difference scheme is the
22、well-known CFL condition given by the equation (ref. 8) The local maximum time step at each grid point is calculated by using equa- tion (12), and then a fraction of the local value (usually 50 to 80 percent) is applied at each grid point to advance the solution in time. This leads to an inaccurate
23、representation of the time-dependent nature of the .solution, but it has been found to accelerate convergence by a factor of approximately 1.5 to 2 when compared with using the minimum global time step. Boundary and Initial Conditions The boundary condition of no flow normal to the body surface is s
24、ince the body surface is the inner boundary of an orthogonal coordinate system. The remaining properties at the surface are computed using three-point forward 8 Provided by IHSNot for ResaleNo reproduction or networking permitted without license from IHS-,-,-differences to calculate the spatial deri
25、vatives in both the predictor and cor- rector steps of the finite-difference scheme described in the previous section. At the outflow boundary (5 = 1) and the shock wave (n = 11, three-point backward differences (away from the boundary) are used. At the shock wave, the pressure thus calculated is us
26、ed with the shock-wave relations (described in the next section) to calculate the shock velocity and then to update all other properties at the shock wave. Along the stagnation streamline (5 = 01, the following conditions are set: u=o and the remaining derivatives are computed by taking into account
27、 the asymmetry in the u-component of velocity. To start the solution, a shock-wave shape is assumed and the shock-wave velocities are set equal to zero. Properties at the shock wave are computed from the shock-wave relations. Properties on the body surface are computed from modified Newtonian theory
28、, assuming that the entropy is constant on the surface. Properties at a constant value of 5 across the shock layer are approximated by linearly interpolating between the values at the shock wave and the body surface. Solution for Shock Wave A method similar to that used by Weilmuenster and Hamilton
29、(ref. 9) has been used to track the movement of the shock wave during the transient portion of the solution. The method can be summarized as follows (details are given in the appendix). First the pressure is computed on the downstream side of the shock wave by using the two-step difference scheme de
30、scribed previously. With the pressure at the shock wave ps known, the density is computed from the equation 1 1 2 Im + -(Ps - PJ 9 Provided by IHSNot for ResaleNo reproduction or networking permitted without license from IHS-,-,-and then the shock velocity is computed with the equation where V, and
31、(Vn)oo are defined in the sketches in the appendix. Next, the movement of the shock wave is computed from the following differential equations: dZS - = V, sin Bs dt dr S “ - -vs cos Bs dt where Bs = tanl (drs/dZ) by using the following two-step finite-difference approximations: Predictor step: - (,)
32、t+ = z: + At (vS sin Bs) (rs)t+At = rk - At(vs cos 6,) 10 Provided by IHSNot for ResaleNo reproduction or networking permitted without license from IHS-,-,-Corrector step: The coordinate system generator is used to update the shock wave after the predictor step and the entire coordinate system after
33、 the corrector step. AS the solution approaches convergence, the shock velocity approaches zero and the shock wave remains essentially fixed in space. The shock velocities at convergence are generally less than in magnitude. Convergence and Computational Time The average change in density from one t
34、ime step to the next was used to determine convergence. Based on the results of several test cases, it was found that the solution was converged when the average change in density was less than 2 x 1 Om5. For the cases presented in the present paper, this generally required from 400 to 600 time step
35、s. From timing studies on the Control Data CYBER 175 computer system, it has been found that approximately 3.4 x 1 0-4 sec/point/step were required to compute the flow field and approximately 7.8 x 1 0-4 sec/point/step were required to com- pute the coordinate system. Thus, in the present nonoptimiz
36、ed computer code, approximately 70 percent of the time is spent computing the coordinate system. One obvious method of improving the overall efficiency is to update the coordi- nate system less frequently. For example, it is possible to update the shock wave after each time step and hold the remaind
37、er of the coordinate system fixed for several time steps before recomputing. This would greatly reduce the per- centage of computational time spent on generating the coordinate system. Body Geometry Most of the body geometrics that are considered in the present paper have a projection in the 2,? pla
38、ne that is described by the equation for a general conic section 11 Provided by IHSNot for ResaleNo reproduction or networking permitted without license from IHS-,-,-where Rb eter . The section. I generates is the radius of curvature at z = 0 and Bb is a bluntness param- I bluntness parameter Bb cha
39、racterizes the eccentricity of the conic Its significance is better understood if it is noted that Bb o generates an ellipse (with Bb = 1 for the special case of a circle). For an ellipse, Bb is related to the ratio of major to minor axis (b/a) by the relation Some of the body shapes considered in t
40、he present paper were spherically blunted cones. These body shapes were generated using equation (20) for the nose, with Bb = 1, followed by a straight-line segment tangent to the nose at the appro- pr iate cone angle. One body shape considered had a reverse curvature at the nose matched to a 30 con
41、ical afterbody at the point 9 = 0 and r = 1 and defined by the equation 0 RESULTS AND DISCUSSION In this section, results of computation by the present method are compared with experimental data and with the results of other computational methods to demonstrate the applicability of the numerically g
42、enerated orthogonal coordi- nate system of reference 6 for flow-field calculations. Sphere Pressure distributions over a sphere for y = 1.4 are presented in fig- ure 4 for free-stream Mach numbers of 2, 3, 4, 6.05, and 8.06. Experimental data from reference 10 and computed results from reference 3 a
43、re presented for comparison. The surface pressures calculated by the present method agree well with the other results for each Mach number considered. The shock-layer thickness for these conditions is presented in figure 5. The agreement with the other results is excellent except at M, = 2 where the
44、 experimental data fall slightly below the calculated results. Since shock-layer thickness changes rapidly at low Mach numbers, this behavior could easily be due to a small change in the effective Mach number for the experimental tests which would change the shock-layer thickness but would have litt
45、le effect on the mea- sured pressure distribution. 12 Provided by IHSNot for ResaleNo reproduction or networking permitted without license from IHS-,-,-The calculated shock-wave shape and sonic-line location about a sphere at y = 1.4 are presented in figure 6 for free-stream Mach numbers of 6.8 and
46、4.76. Experimental data from references 11 and 12 are presented for comparison. For both conditions, the agreement with the experimental data is very good. Paraboloid The calculated surface pressure distributions over a paraboloid at y = 1.4 are presented in figure 7 for Mach numbers of 3 and 10. Ca
47、lculated results from references 3 and 13 are included for comparison. The results of the present calculation are in very good agreement with the other calcu- lated results . The shock-layer thicknesses for these conditions are compared in figure 8 where, again, excellent agreement with the results
48、of references 3 and 13 will be noted. Ellipsoid Pressure distributions for an ellipsoid with b/a = 0.5 and y = 1.4 are presented in figure 9 for Mach numbers of 4, 6.05, and 8.06. Experimental data from reference 10 and calculated results from reference 3 are presented for com- parison. The present calculated results agree closely with the other data for each Mach number. Pressure distributions for an ellipsoid with b/a = 1.5 and y = 1.4 are presented in figure 10 for Mach numbers of 3, 4, 6.05, and 8.06. Experimental data from reference 10 and calculated res
