- Email: [email protected]
Spectral solution of inviscid supersonic flows over wedges and axisymmetric cones
Computers Fluids Vol. 21, No. 2, pp. 247-266, 1992 Printed in Great Britain. All rights reserved
SPECTRAL SOLUTION OVER WEDGES
0045-7930/92 $5.00+ 0.00 Copyright © 1992 Pergamon Press plc
OF INVISCID SUPERSONIC FLOWS AND AXISYMMETRIC CONES DAVID A. K O P R I V A
Department of Mathematics and Supercomputer Computations Research Institute,The Florida State University, B-186, Tallahassee, FL 32306-4052, U.S.A. (Received 1 March 1991; in revised form 16 July 1991; received for publication 24 October 1991) Abstract--We use a shock-fitted multidomain spectral collocation method to solve both steady and unsteady inviscid supersonic flows over bodies. New aspects of the method include two subdomain interface types and a zonal solution procedure to get efficient convergence to steady-state in supersonic regions. Steady-state examples include flow over a sharp cone, a hyperbolic cone and a hyperbolic wedge. For the cone, the exact flow solution is used to show that the method is spectrally accurate. As an example of an unsteady flow, we present a calculation of the interaction of a free-stream hot ring with the flow over a sharp cone.
1. I N T R O D U C T I O N The use of spectral collocation methods to solve multidimensional compressible flow problems is less than a decade old [1]. The earliest papers [e.g. 2, 3] show only the feasibility of the method and date from the early 1980s. Spectral methods were introduced with the hope that the high-order ("spectral") accuracy that proved to be important in the solution of many incompressible flows could be used to advantage to understand the physics of some compressible flow problems. Because they are somewhat new, techniques associated with the more established methods are only now being applied to spectral methods. When first applied to compressible flows, spectral methods had several inherent limitations. One problem was how to compute shock waves. Without special care, spectral methods reduce globally to low-order polynomial accuracy when shocks are present. To address this problem, both shock-capturing and shock-fitting methods have been studied [4-7]. A second limitation was the inflexibility of the methods. In multiple space dimensions, tensor product interpolations are used to define the global spectral approximation. This means that the domain of interest must be simple enough to be mapped onto a square. Given that mapping, the grid points must then be placed according to a specific, typically nonuniform, distribution. The result is that when steep gradients are present, large oscillations appear in the solution unless many grid points are used. To get past the second limitation, multidomain techniques (among others) were introduced in the mid 1980s [8-10]. In a multidomain technique, like the zonal techniques associated with finite difference methods, the physical domain of interest is divided into multiple subdomains. Local mappings of the subdomains allow the solution of flows on complicated domains. Also, the careful placement of subdomains, and the choice of the number of grid points within them, allows for local refinement of the grid. An additional limitation of spectral methods comes from the stiffness in time of the spatially discretized equations. This problem has yet to be resolved satisfactorily. The derivative matrices are full, so explicit time marching procedures are usually used. The time steps associated with the commonly used Chebyshev spectral methods vary as I/N 2 in one space dimension, where N is the number of grid points. This can be compared with the O(1/N) variation typically associated with finite difference methods. A multidomain discretization helps this somewhat [8, 10]. For instance, a 10-domain decomposition with 10 points in each can take a time step roughly 10 times larger than can be taken on a single domain with 100 points. Also, the work required to compute the derivatives in that multiple domain decomposition will be 1/10 that of the single domain. Except for the overheads associated with the interface procedures, the result is an overall factor of 100 advantage for the 10-domain decomposition. 247
248
DAVIDA. KOPRIVA
In this paper, we extend and apply the spectral multidomain technique of Ref. [10] to solve some inviscid supersonic flows over bodies. The extensions are needed because the flows are more complicated than those computed previously. Unlike the flows computed in Ref. [10], the flows considered here will be rotational and not homentropic. Also, a bow shock is present ahead of the bodies. This bow shock is fitted so that the flow is smooth in the computational region. The multidomain decomposition is used to resolve narrow layers in the solutions and to reduce the stiffness of the complete problem. To extend the method, we introduce two new interface conditions. The first condition arises at the intersection of the fitted shock boundary and an interface boundary. This allows the shock to be divided into multiple segments. We also introduce "shear interfaces" along which the grids slide as the calculation progresses. With shear interfaces, subdomains do not need to be fixed relative to each other. Finally, to make the calculation of steady-state flows more efficient, particularly those presented in this paper, we group subdomains into logical "zones". In supersonic regions of a flow, the solution is marched zone by zone from upstream to downstream. To apply the method, we consider steady and unsteady supersonic flows over two simple geometries. We first compute the flow over a sharp nosed cone (the "Taylor-Maccoll" [11] problem). An exact solution to this problem can be found, and we observe spectral accuracy in the computed solution. The second problem is that of flow over a hyperbolic wedge or cone. For long bodies, the significant computational difficulty is the presence of a very narrow layer near the body in which the entropy and vorticity vary strongly. The finite difference approximation of entropy layers was considered in Ref. [12]. We show that a multidomain discretization allows the spectral method to resolve the high gradients in these thin layers. Finally, we consider the interaction of a free-stream hot ring with the flow around a sharp cone as an example of a time-dependent flow. This paper is organized as follows. In Section 2 we present the equations and the mappings. We also derive the shock acceleration formula for time-dependent free-stream variations as an extension to the work in Ref. [7]. In Section 3, we review the multidomain method. We then introduce the new shock-interface condition, the shear interface condition and the zonal approach used to compute the steady-state solutions. Section 4 contains examples. Finally, some conclusions are drawn in Section 5. 2. T H E E Q U A T I O N S , T H E M A P P I N G S AND S H O C K F I T T I N G In this paper, we solve flows around bodies placed in a supersonic free stream. The two geometries considered are shown in Figs l(a, b). The first situation, Fig. 1(a), is used for the solution of the sharp cone problems. We avoid the numerical difficulties at the sharp nose and specify supersonic inflow conditions at a plane located downstream of the nose and perpendicular to the body. Figure l(b) represents the axisymmetric blunt body problem for which single domain spectral results were presented in Ref. [7]. In this problem, a portion of the flow between the body and the shock is subsonic. Also, entropy and vorticity are generated near the nose where the shock is strongly curved.
Shockw
~
M>I (
(b) Fig. 1. Diagrams of the computational regions. (a) Sharp cone geometry. (b) Blunt body geometry.
Inviscid supersonic flows over bodies
249
The inviscid flow between the bow shock and the body can be represented by the Euler gas-dynamics equations:
P,+UPx+VPy+ 7(Ux+Vy+~)
=0,
u, + UUx + vuy + "rPx = O, vt + UVx + vvy + "rPy = O, St+uSx+vSy=O.
(1)
Here, P = In(p) is the logarithm of the pressure. The velocity is q = (u, v) and r = a2/y is the temperature. The entropy, S, is scaled to Cv, and a is the sound speed. We take 7 = 1.4. The pressure and density are scaled to the free-stream values. Velocities are scaled to ~ . For the blunt nose problems, the lengths are scaled to the nose radius. The parameter ~ = 0 for the planar case, and ~ = 1 for the axisymmetric case. A transformation to computational coordinates is done in two steps: (x,y, t)~--~(r,s, t) (X, Y, T). The first mapping is to a coordinate system where s approximates the arclength along the body (with s = 0 at the nose), and r is the normal distance from the body. In (r, s) space, the region between the shock and the body is divided into multiple subregions (subdomains). Then these subdomains are mapped individually to the unit square, (X, Y)E [0, 1] x [0, 1]. The mapping (x, y ) ~ (r, s) is made by the relations x = Xb(S) + r~" ~, y = yb(S) + r~ .y; where (Xb, yb) represents the body as a function of s. As an example of this mapping, we consider the hyperboloid (or hyperbolic wedge) of Section 4. There Yb = tan(~b)x/(x + x0) 2 -- e2, where • = 1/tanZ(~b) and x0 = • + 1. This defines a hyperbola with unit radius of curvature at the nose and center at x = 0. Asymptotically, as x --* oo, the body is straight and with angle q~ to the horizontal• A reasonable parametrization of the distance along the hyperbola is in terms of y, so we define yb(S) = S, Xb(S) = .~/(s /m ) ~ + a z -- x0; where m = tan(q~) is the asymptotic slope of the body. After some algebra, we have the mapping r
x = xb(s) +
,/ l+\ay/ (ox y' Oxb
r ay
where CqXb
S
Oy
m2[xb(s)+ x0]
In (r, s) coordinates, we create multiple nonoverlapping subdomains, Gk, which collectively cover the region between the shock and the body (Fig. 2). Each subdomain is bounded by four curves. In r, the bounds are rmin(s, t, k) and rr~(S, t, k). Time-dependent subdomain boundary positions are included to allow the shock boundary to move with time. In a later paper, we will also consider the use of time periodic body motion to generate acoustic waves at the body surface. For convenience only, we require the subdomains to be bounded in s by straight lines normal to the body surface. CAF 2|/2--H
DAVID A. KOPRIVA
250
Point
Fig. 2. D i a g r a m of the d o m a i n decomposition. Each s u b d o m a i n , Gk, is b o u n d e d by four curves: r = rmi,(s, t, k), r = rmax(S, t, k), s = stain(k), s = smx(k ). T h e new interface points introduced in this paper are (c) and (d).
Given the boundary curves of a subdomain, a linear mapping from (r, s, t) ~ (X, Y, T) is used: r -- rmi.(s, t, k) rmax (S , t, k) - rmin(S ,
(2a) t, k ) '
S - - Stain ( k )
r~=
(2b)
Smax ( k ) - - Stain ( k ) '
T=t.
(2c)
Under the complete transformation from a subdomain to a square, (x, y, t) ~ (X, Y, T) [0, 1] × [0, 1], the equations become
Pr+ UPx+ VPy+ 7[Xxux+ X, vx+ Yxur+ Yyv~+~]=O, a 2
u r + Uux+ Vu~+-- (XxPx+ YxPr) = 0,
7
a 2
Vr + Uvx + Vvr + - - (XyPx + YyPr) = O, Sr + USx + VSr = O;
(3)
where
u = x, + UXx + VXy, V=uYx+vYy. The X, term appears because of the moving shock. Finally, the metric coefficients can be written in terms o f the intermediate metric terms as
Xx= XsSx + X, rx, Xy = Xssy + X~G,
Inviscid supersonic flows over bodies
251
and Yx = Yssx+ Y,r~, Yy = Yssy + Y, ry.
To ensure that the spectral approximation is not used to approximate derivatives across the shock, the shock is fitted as a boundary. The use of shock fitting with spectral methods is based on the work of Moretti [13]. It was shown to be viable in Refs [2, 5]. Kopriva et al. [7] have shown that the simplified shock-fitting procedure used in the early papers is not stable when used with spectral methods to compute steady flows. Instead, they derived a fully characteristic approximation. With the new shock-fitting formula, spectral accuracy of the stagnation pressure was shown for the blunt body problem. In this paper, we include the possibility of a time-varying free stream. So, we rederive the shock acceleration formula presented in Ref. [7] to account for free-stream variations. To fit the shock as a boundary, we set the shock position, denoted rsh(S, t), to be a boundary of the flow. Thus, for a subdomain that is bounded by the shock, rmax(S, t, k) = rsh(S, t). As the flow develops, the shock moves with velocity w = w~, where ~ is a unit vector in the r direction and w = Orsh/Ot is the shock speed. The shock speed is computed by integrating the shock acceleration in time. The shock position is updated by integrating the shock speed in time. To obtain the shock acceleration formula, we differentiate in time the Rankine-Hugoniot relations for the jump of the normal velocity and pressure across the shock. In what follows, let subscript 1 represent the low-pressure (known) quantities upstream of the shock; subscript 2 will represent the downstream or unknown flow quantities. For a time-varying free stream, the Rankine-Hugoniot relations are: ,n~----~--- TI)
(4a)
and ),-1 2a~ 1 ($2 = ~ ($1+ ~ 7+ 1 (st--'
(4b)
where ($i = q i - N - w(~-N) is the normal gas velocity relative to the shock and N = N ~ + Nz.9 is the shock normal. Since equations (4a, b) hold for all time, T, we differentiate them with respect to T to give P: = Pl + G31 + H
(5a)
(~2= F61 + E.
(5b)
and
Here, we use the dots to represent partial time derivatives, (') = O/OT. The coefficients E, F, G and H are related to quantities on the low-pressure (prescribed) side: E ~
-
27~ -
(~ + l)a,'
F=7- 1
27T1
~+l
(~ + 1)($~' 2($1
G=
and I H= -
7-- 1 2
a~
+ - -1
~-~ -- --~"-
71 "~1
252
DAVIDA. KOPRIVA
The derivative of thc relative velocity in terms of the shock accelaration, ~., is 42" N = w(~. N)(1 - F) + F(q,. ~ + ti," N) + w~. l~l(l - F) - q2" i'q + E.
(6)
The pressure variation is related to the shock acceleration by P2 = Pl + G[ql "~ + ¢h .N - w~.(~] - w~.NG + H.
(7)
To relate equations (6) and (7), we use the compatibility condition for the waves hitting the shock from the high-pressure side:
P2 - L~:. N + C = 0,
(8)
a2
where C = (qz - aN)" V P + Rp - a
(N,R. + NzR~),
V P = (XxPx+ YxPy).¢c + ( X y P x + YyPy).f,
Rp=,(Xxux+
Xyvx+ Yxuy+ Y y v r + ~ ) +
X, Px.
R. = Uux + Vuy and
R~ = Uvx + Vv y. Finally, equations (6) and (7) are substituted into equation (8) and the result is solved for the shock acceleration: -
aPl + (q~. Iq + dh'N ) (aG - ~F) - w ~. N(aG - 7F + 7) + ~ q2" N + a ( C + H) + yE ~. N(aG - 7F + ~')
(9)
The upstream quantities, P~, ~h and ~ used in equation (9) come from the time variation of the upstream flow. We define the flow upstream of the shock to be the uniform free stream plus a perturbation, e.g. P~ = P~ + eP'(x, y, t'), where the perturbation is defined in its own rest frame. Since the flow is supersonic, the perturbation is simply advected along with the flow until it interacts with the shock. In the rest frame of the shock, then, the time derivatives become
0,
PI
= E/-~TT.. + (q, + w ) ' V q '
LotF~q' L c3t'
]
1
(ql + w).VP' ,
and [-c3z'
e~ = , 1 _ ~ +
,
] (ql + w ) . W '
.
3. THE M U L T I D O M A I N S P E C T R A L M E T H O D 3.1. Discretization and interface conditions To discretize the system of equations (3), we start with the multidomain Chebyshev spectral method of Kopriva [10]. The method is general in the sense that there can be any number of subdomains defined in any order. A description of the data structures is given in Ref. [9]. Of practical importance is that local refinement of the grid is possible. Subdomains can be of any size. The number of grid points within a subdomain is independent of the number of points in its neighbors. Finally, subdomains can be overlapped, if necessary, but we do not use this feature in this paper.
Inviseid supersonic flows over bodies
253
Within each subdomain, the spatial derivatives in the equations are approximated by a standard Chebyshev spectral collocation method [1]. The result of the approximation is a system of ordinary differential equations for the unknowns at each grid point. The system is advanced in time by one of two time integration procedures. For the time-accurate calculation, we use the second-order midpoint rule. For unsteady problems, the low order of the time differencing procedure is not a serious problem. The stiffness of the discrete equations in time already requires the use of small time steps. For the steady-flow problems we are not interested in time accuracy. So, we use the first-order accurate four-stage Runge-Kutta method [1] commonly used in finite volume integrations. This allows a larger time step and more efficient time integration than does the midpoint rule, yet it has the same storage requirements. Numerical experiments show it is stable for a Courant number of at least 3, based on the smallest grid spacing. The approximation of the shock normals and their derivatives needed for equation (9) is done spectraUy. The two outward shock normal components can be written as /6~rsh
~/
/c~r~h
\ /
and
where
A
/(ar,h sx-rx )2
v
The shock slope, Or,h = -_r,,, -~ y,, ds dY is computed by the one-dimensional Chebyshev derivative of the shock position array. Thus, it is computed to the same accuracy as are the derivatives of the dependent flow quantities. Similarly, the time derivative of the shock normals requires the computation of dw/as, which is also computed spectrally. The derivation of interface conditions between subdomains for various types of interface points is described in detail in Ref. [1(3]. We only summarize those results here. The equations are solved on each subdomain to get the interior point solutions and a (provisional) value of the flow quantities on the boundaries. Each interface point has at least two provisional quantities, one from each subdomain to which it is a neighbor. These provisional quantities are then projected onto inflow and outflow components. To project the solutions, four generalized Riemann variables are computed from the infinite number represented by P + (~/a)q. ~ where g is a unit vector. The four choices are for ~ = + N and ~ = + X, where N and Z represent the normal and tangent vectors to the interface. To obtain the correct domain of dependence, the Riemann variables are then selected according to the subdomain from which their associated bicharacteristic comes. If the grid points do not match across the interface, the incoming quantities are interpolated using the natural Chebyshev interpolation that defines the method. (See Ref. [9] for details.) The resulting Riemann variables are then recombined to get corrected values of the flow quantities. As an example, we consider the computation of the flow quantities at the intersection point between the three domains GI, G2 and G3 shown in Fig. 3. In Ref. [10], this was called a "T" intersection. Superimposed on the grid, we show a projection of the bicharacteristic cone that affects the interface point. We represent the base of the bicharacteristic cone by the circle. Since the length of the fluid velocity vector, q, is shown as smaller than the sound speed represented by the radius of the circle, Fig. 3 shows a subsonic case. Four bicharacteristics have been chosen corresponding to ~ = ___~, ___)3.According to this choice, R + --/3 + (v/a)a and S + = 13 + (v/a)g, computed from the provisional values/3, a and g at the point in subdomain G3, and the sound speed evaluated at the previous time level. We compute S - from provisional values in GI. R - is computed from values interpolated from G2 to the interface point,
254
DAvID A. KoPRIVA
g-4--/ l q \ IIGI
F3!
! !!
Fig. 3. Determination of the domain of dependence in a subsonic flow at the intersection of three subdomains.
The interface condition requires that the true Riemann variables at the interface equal those computed from the provisional quantities [10]. Thus, we have four equations
p + 7-u=R+, a p -Tu=R-, a p + aY-v=S +, and
P - Ya- v = S - , from which to calculate the three quantities, P, u and v. Our choice of the four bicharacteristics makes the calculation of the velocities simple: a
u = Z (R + - R - )
and a
v =T~ (s÷ - s-).
The pressure is computed by combining all four Riemann variables and subtracting the value computed in the upstream subdomain, G3 (see Ref. [!0] for details): e =~(R + +n-
+s + +s-
- 2P3).
Finally, the entropy is advected along streamlines so its domain of dependence is given by the vector q shown in Fig. 3. Thus, the corrected value of the entropy is just the value computed in G3, i.e. S = ~ . To compute the flows considered here, we introduce two new types of interface points. These points are represented in Fig. 2. The first we call "shock-interface points" and they occur where the shock is divided at an interface. Many examples of such points are shown in Section 4. We remark that shock-fitted multidomain solutions have been presented in the past [14, 15]. See also Ref. [1, Chap. 13]. The subdomains in those papers were always strips along the shock direction. With the shock-fitting and interface procedures used in those papers, shock-interface points were not stable. At shock-interface points we use the shock normal and characteristic information from the upstream subdomain. The difficulty at shock-interface points is that there are two values of the shock normal, one from each subdomain, from which the correct one must be chosen. The shock acceleration formula (9) uses the compatibility equation (8) for the normal bicharacteristic. Since that bicharacteristic originates upstream of the shock-interface point, we use the normal and shock acceleration computed on the upstream side of the interface to determine the flow.
Inviseid supersonic flowsover bodies
255
We call the second type of interface a "shear interface" (Fig. 2). It represents a straightforward extension to the interface treatment in Ref. [10]. At a shear interface, neighboring subdomains are allowed to slide relative to each other as a function of time. A shear interface is convenient to use to solve the blunt body problem, for example. The narrow entropy layer that needs to be resolved near the body appears far downstream of, not near, the nose. The use of the shear interface allows for a change in the grid topology from one that does not need to resolve a narrow boundary layer to one that does. To define the shear interface, we use the hyperbolic nature of the problem to consider each subdomain as independent of the others. Then, the interface conditions are nothing more than boundary conditions on the incoming characteristic quantities. If the subdomains do not move relative to each other, the grid connections and interpolations for these characteristic quantities are defined once, at the start of the computation [9]. For the shear interface, the connections and interpolations are updated at each time step. The shear interface is the most expensive of the interface types and is used sparingly. The main expense is in the calculation of the Lagrange interpolating coefficients as described in Ref. [9]. To make this calculation as fast as possible, we no longer compute these coefficients directly using the definition of the Chebyshev polynomial, T, (x) = cos[n cos- l(x)]. They are now defined recursively by the standard three-term recursion for the Chebyshev polynomials [1]. 3.2. Acceleration to steady state To decrease the computer time needed to calculate the steady-state solutions, we introduce the notion of zones that are groups of subdomains to be computed together. This notion is convenient for the problems computed here since the solutions develop from the nose to the downstream and most of the flow is supersonic. In the supersonic portions of the flow, interface conditions are only needed on the upstream side of a subdomain. If the solution in a subdomain upstream of a given subdomain is converged to an acceptable level, there is no need to continue to compute on the upstream subdomain. In other words, subdomains within upstream "zones" can be turned "off" and not be computed further. Finally, subdomains in downstream zones do not have to be computed until the zones upstream of them are converged. Zones are defined in increasing numerical order from upstream to downstream. They define a purely logical grouping device that tells which subdomains are to be computed simultaneously. The interface conditions are not affected in any way. Note that all subdomains that contain subsonic points must be grouped in a single zone and solved together to convergence. For the calculations that we present here, the subsonic region is only a small fraction of the total flow. There are several advantages to be gained by solving by zones. First, only a small fraction of the total number of grid points needs to be solved at a given time. Next, the time step is regulated only by the subdomains that are turned on. If all subdomains are computed simultaneously, the time step must be determined by the minimum over all the subdomains. The minimum is usually the time step at the stagnation point. It can be an order of magnitude smaller than the local time step in the supersonic parts of the flow. Finally, the concept of zones helps to avoid the need for a highly accurate initial shock shape. For long bodies, a small error in the initial shock shape near the nose can translate into a large error far downstream. This has typically lead to a kink instability in the shock as the computation proceeds. To avoid the problem, the solution in a supersonic zone is "bootstrapped" from the previous zone by linear extrapolation. Far downstream, the shock is nearly straight so this is a reasonable initial approximation. 4. EXAMPLES We now present three examples of the use of the multidomain method described in the previous section. The first is the Taylor-Maccoll problem [11] for which an exact solution can be computed. The second is the steady flow over a blunt hyperbolic wedge and cone. The unique feature of the blunt body flow solutions is the narrow entropy and vorticity layer that develops downstream of the nose. Finally, we include a time dependent calculation that models the interaction of a free stream hot ring with the flow around the sharp nosed cone.
256
DAVIDA. KOPRIVA
Fig. 4. Geometry and grid topology for the solution of the Taylor-Maccoll [11] problem.
4.1. The Taylor-Maccoll [11]problem The first problem we consider is the steady flow over a sharp nosed cone at zero angle of attack. We assume that the semi-vertex angle is small enough so that the shock is attached to the cone. Then, the shock is straight and the flow solution is conical. Also, the entropy is constant and the vorticity is zero in the region between the shock and the cone. To find the exact solution to this problem, we solve a system of ordinary differential equations [16]. A Richardson error estimate is used to ensure that enough grid points are used to obtain machine accuracy. To study the error, the problem was solved on the three subdomain grid topology shown in Fig. 4 for several values of the free-stream Mach number and semi-vertex angle. The grid was chosen to show the accuracy of the method when shock-interface points are used. The upper boundary corresponds to the fitted shock wave and the lower boundary is the body. Along the body, the usual charateristic boundary conditions were used [cf. 7]. The inflow boundary at the lower left is supersonic, so all flow variables were prescribed from the exact solution computed at the inflow points. The outflow is also supersonic so no boundary conditions were applied there. As the initial condition, we specified the exact solution. The solution was advanced in time until the shock speed converged. Spectral accuracy was observed for the Taylor-MaccoU [11] problem. As an example, we present the results for the M = 6 flow over a 10° cone. Figure 5 shows the maximum errors of the shock slope, the pressure, the horizontal and the vertical velocity as a function of the number of grid points in the direction normal to the body. The solution is analytic, so we expect exponential decay in the error [1]. All quantities do show exponential decay in the error.
4.2. The blunt body problem We consider next the solution of flows over hyperbolic wedges and cones. The (single domain) spectral solution of similar blunt body problems was described in Ref. [7]. Here we are most interested in the flow over bodies that are long compared with the nose radius. The computational difficulty is to resolve the narrow entropy and vorticity layers that form along the body [12]. To
o
_4,0 ~
,~: ~ ' ~
"~
-~-Slope I ~Pressurel
~
uJ
E
-s.o
E
I I I
~
-6.0 o 0 -J
-7.0 -8.0 ~
3
4
~
5
6 N
....
7
8
J
9
Fig. 5. Error decay of the shock slope, pressure and velocity components for the Taylor-Maccoll [11] problem.
Inviscid supersonic flows over bodies
257
resolve these layers, we place narrow subdomains along the body. We remark that this is similar in spirit to numerical approaches suggested in Ref. [12]. To start the calculations, a reasonable initial solution was used. A parabolic shock shape was assumed near the nose. The Rankine-Hugoniot relations were used to specify the flow variables along the shock• Along the body, the initial condition assumed that the Math number varied as sin(s) up to the sonic point and constant beyond that. The sonic point was placed along the body at a point near the nose. The interior flow variables were then defined by linear interpolation. In all the cases presented in this section, the solutions were converged to maximum shock speeds between 10 -5 and 10 -4 depending on the free-stream Mach number. Solutions were computed for a wide range of parameters. Asymptotic body angles ranged from 10° to 30° for both wedge and cone geometries. Mach numbers varied between 4 and 6. Finally, body lengths ranged from about 2 to 225 nose radii. When possible, solutions were compared with the single domain calculations using the single domain spectral code described in Ref. [7]. As an example, we first present a solution in the subsonic portion of the M = 4 flow over a 10° hyperbolic wedge on a grid with three subdomains. Figure 6 shows the vortieity contours and grid for the initial and converged solutions. The subdomains were chosen so that a shear interface and shock-interface point would occur within the subsonic portion of the flow. Notice that the grid lines have changed relative to each other as the solution converged.
rid
I n t t• "
Conv
Fig. 6. Vorticity contours and grid for the solution near the nose of the M = 4 flow over a hyperbolic wedge. The sonic line is indicated with the vorticity contours.
258
DAVID A. KOPRIVA
A severe test for any numerical method, including this one, is to see how smooth the vorticity contours are. Vorticity is a quantity computed after the fact from the derivatives of the velocities. Thus, the errors in the vorticity will be larger than the errors in the velocity components alone. To compute the vorticity, we use the usual Chebyshev differentiation of the velocities. The results of this calculation are the contours shown in Fig. 6. The contours were not smoothed by the plotting routine. Thus, we see that the multidomain method gives a solution that is not only continuous, but smooth through the interfaces and throughout the flow. As a final measure of the quality of the solution, we note that the vorticity along the wall, which should be 0, is at most 4 × 10 -3. The shape of the vorticity layer further down the body is determined by the distribution of the shock-generated vorticity in the nose region shown in Fig. 6. Since the shock becomes nearly straight in the neighborhood of the symmetry line, the vorticity approaches zero. In that same region, the entropy is nearly constant. Then Kelvin's theorem applies and the vorticity lines will follow the streamlines. Away from the symmetry line, the shock becomes strongly curved, the entropy varies strongly (cf. Fig. 7) and the vorticity generation reaches its maximum. Finally, away from the nose, the shock straightens out and the vorticity again approaches zero. Note, that, in the region of the curved shock and strongly varying entropy, Kelvin's theorem does not apply and
Pressure
Mach
Entropy
Fig. 7(a)
Inviscid supersonic flows over bodies
Fig.
259
7(b)
Fig. 7. The grid contours of the pressure, M a t h number, entropy and vorticity for the steady M = 6 flow over a 25 ° hyperbolic wedge.
the vorticity contours do not follow the streamlines. Thus, we can predict that profile between the body and the shock along a normal to the body of the vorticity layer far downstream of the nose. The vorticity will rapidly drop from zero at the body to a minimum and then more slowly rise to zero at the shock. Figure 7 shows the solution contours and grid for an M = 6 flow over a 50 nose radius long 25 ° hyperbolic wedge. Eight subdomains and a total of 800 grid points were used. The grid included both a shear interface and several shock-interface points in the supersonic portion of the flow. A narrow subdomain layer was used along the body outside the nose region to ensure uniform resolution of the entropy and vorticity layer. Figure 8 shows the variation of the solutions along the normal direction at the outflow boundary of the calculation. The subdomain interface appears at the distance r = 2 from the body. The pressure varies by < 2% across the flow so it is clearly over-resolved by this grid. More structure appears in the velocity [=(u s +v2)~/2], the entropy and the vorticity. For instance, the shock slope is zero at the symmetry axis. So, the derivative of the velocity and the entropy should be zero at the body. This feature is captured by the method. Also, the vorticity is smooth and has the structure suggested by the variation of the vorticity generation near the nose. It is zero at the wall and ends at zero at the shock. In this calculation the vorticity error along the wall is only 6 x 10 -3.
DAVID A. KOPRIVA
260
14
14
. . . .
I
. . . .
I
. . . .
I
. . . .
I
12
10
10
8
8
6
6
4
4
2
2
0 13.15
0 0.5
14
.
13.2
.
.
.
I
.
.
13.25 13.3 Pressure
1 3 . 3 5 13.4
.
I
.
I
.
.
.
.
.
.
.
4
.
(b)
12
~ - - ~ i'
(a)
12
,
. . . .
0.7
P
12
t
*
I
. . . .
0.9
1.1 Entropy
F
. . . .
~
. . . .
I
. . . .
I
1.3
. . . .
~
. . . .
.5
r
. . . .
v '
''
(d)
10
10 '
8
8
6 4
4
2
2
0
'
4
w
'
4.5
5 Velocity
5.5
0
6
'
'
-
-0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 Vorticity
0
0.1
Fig. 8. The variation of the solutions along the outflowplane from Fig. 7. A more difficult calculation is presented in Figs 9 and 10 which show the M = 4 flow over a 10° hyperbolic cone of length 100. This time, 18 subdomains and over 1600 grid points were used. As before, the calculation includes a shear interface and multiple shock-interface points. Compared to the distance between the shock and the body, the entropy and vorticity layers are much narrower than in the previous example. Figure 11 shows the variation of the solution at the outflow plane of the calculation. Three subdomain divisions were used in the direction normal to the body to resolve the structure of the entropy layer. As before, the pressure variation shows little of interest near the body. The velocity, entropy and vorticity, on the other hand, show strong gradients there. Note that, even with this grid, we have not resolved all the structure of the layer. Missing from the velocity and entropy curves are the turnover to zero slope at the body. The vorticity graph has missed the zero vorticity along the wall. These features were missed in part due to the lack of refinement of the grid near the nose. What we see in this calculation, however, is the ability of the method to allow the computation and resolution of very steep, though continuous, features. 4.3. The interaction o f a hot ring with the flow over the sharp nosed cone
The interface conditions, being characteristic, are also applicable to time-dependent problems, particularly in the important acoustic regime. An example of a linear problem was considered in Ref. [10]. For nonhomentropic flows, however, the sound speed used to compute the Riemann variables is linearized in time about the previous time level. As described in Ref. [10], this linearization introduces errors of the same order as the time stepping scheme used. The interface conditions are not applicable, again because they are of characteristic type, to the passage of strong nonlinear waves with unresolved steep gradients such as shocks. In these situations, numerical experiments show that the waves are stopped at the interfaces.
Inviscid supersonic flows over bodies
261
Grid
Fig. 9. Solution contours and grid for the M = 4 flow over a 10° hyperbolic cone.
Our final example uses the method to compute a time-dependent rotational and nonhomentropic flow. We present a calculation of the interaction of a hot ring with the M = 4 flow over a 30 ° cone. Figure 12 shows the geometry and grid. The section of the cone ranged from s = 1.0 to s = 2.0. Eight subdomains and a total of 1736 grid points were used. To begin, we used the cone geometry of Section 4.1 and established a steady state. Then the temperature in the free stream ahead of the shock was modified to r = 1 + E exp( -~t[(x - x0) 2 + (y - y0)2]}. The temperature ring was allowed to advect with the free stream until it encountered the shock wave. This problem lacks an exact solution. However, the numerical solution of a similar s h o c k temperature spot interaction problem was obtained in Ref. [17] with a shock-fitted finite difference scheme. It was solved again in Ref. [15] by a multidomain spectral method that could not split the shock by a subdomain interface. Neither paper included the complications of a wall boundary. Both computations showed the generation of a sound wave, an entropy spot and a pair of counterrotating vortices. The shock-fitted spectral method described here is suitable for the study of the interaction of the shock with low-amplitude free-stream perturbations. High-amplitude perturbations create large variations in the shock shape and the possibility of additional shocks that cannot be captured by the spectral method. Therefore, calculations were performed with amplitude, e, between 0.05 and 0.2. The temperature perturbation was centered on (x0, Y0) --- (0.9, 1.05). The coefficient ~t = 50 was chosen so that at t = 0, the flow just ahead of the shock was consistent with the initially straight shock. Also, with these parameters, the spot was localized enough not to affect the shock at the inflow shock point on the upper left of Fig. 12.
262
DAVID A. KOPRIVA
Pressure
Mach
Fig. 10. Blowup of the nose region of Fig. 9.
The results show the same qualitative features as the previous calculations in Refs [15, 17]. The time history of the result of the interaction for E = 0.05 is shown in Figs 13-15. First, in Fig. 13, we see the generation of a sound wave that strikes the body and moves upstream along it. Figure 14 shows the generation of an entropy spot at the shock. This spot, which is flattened in the direction normal to the shock, breaks away and moves upstream along the streamlines. For our purposes here, we are most interested to see that the entropy is advected smoothly through the interfaces and with no visible dissipation or dispersion. Dissipation would show up as a decay in the maximum entropy and an increase in the size of the core of the spot. Over the 600 time steps taken between t = 0.156 to 0.250, the maximum entropy varies by at most 0.14%. Dispersion effects would be seen as waves trailing (or possibly leading) the entropy spot and are not visible in the plots. The time history of the vorticity is shown in Fig. 15. We see the generation of the counter-rotating vortices represented by positive and negative spots of vorticity. In this nearly linear regime, they are of nearly equal strength. The worst difference in the maximum and minimum vorticity is 1.3% at time 0.188. The cores of the vortex contours do grow with time, so we see that the vortices show the effect of more dissipation than does the entropy. However, they do not show any significant trailing waves that are characteristic of large phase errors, and they propagate smoothly through the interfaces with no visible reflection.
Inviscid supersonic flows over bodies 0
. . . .
I
. . . .
L
. . . .
I
. . . .
I
0
. . . .
_(a)
16
....
, ....
, ....
,
,
. , ....
, • . ,. , ....
(c)
16
12
263
12 t..
8
8
4
0 1.7
1.75
1.8
1.85
1.9
.95
0
0.1
0.2
0.3
0
20
_(b)
16
0.4
0.5
0.6
0.7
0.8
Entropy
Pressure
. ,
. ,
• ,
.
(d)
16
12
. ,
12 t..
8
8
4
4
0
, 3.8
'
'o~ 4
'
~- ' 4.2
^ ~ 4.4
Velocity
0
4.6
4.8
~
-8
.
.
~
-6
.
.
.
I
.
-
-4
.
t
-2
Vorficity
Fig. 1I. The variation of the solutions along the outflow plane from Fig. 9.
Fig. 12. The grid used for the calculation of the hot-ring calculation at the initial time.
~
,
~
.
0
.
•
264
DAVID A. KOPRIVA
TIME = 9.387E-02
= •
TIME =
.125
=.
TIME =
.156
=.
Fig. 13. Pressure contours for the hot-ring calculation.
TIME = .188
T
m
-
•
Fig. 14. Entropy contours for the hot-ring calculation.
Inviscid supersonic flows over bodies
265
"1"1
i Fig. 15. Vorticity contours for the hot-ring calculation. The solid lines represent positive vorticity. The dashed lines indicate negative vorticity.
5. C O N C L U S I O N S We have used a multidomain spectral method to solve some inviscid supersonic flows over some simple bodies. Unlike the solutions previously presented in Ref. [10], we have computed flows that have entropy and vorticity variations. Two interface procedures were introduced. One allowed the fitted shock to be intersected by subdomain interfaces. The other allowed interfaces to slide along each other. To accelerate the convergence to steady state and to use locally determined (by subdomain) time steps, subdomains were grouped into zones and marched in time by zone. The solutions obtained by the method were accurate and smooth. First, we solved the Taylor-Maccoll [1 l] problem. It has an exact, analytic solution. As expected, the errors decayed exponentially fast. The more complex problems such as the flow over the hyperbolic bodies had entropy and vorticity generated by the curved shock. By careful choice of the subdomain decomposition, it was possible to resolve steep gradients in these quantities. Though the interface conditions specify only the continuity of the characteristic variables, we see that even the vorticity is smooth across subdomain interfaces. The time-dependent calculation showed the generation of a counter-rotating vortex pair and an entropy ring. These features also propagated smoothly through the subdomain interfaces. Acknowledgements--The author would like to thank Drs C. L. Streett and T. A. Zang for helpful discussions. This research was supported in part by the National Aeronautics and the Space Administration under Grant No. NAGI-862 and by the U.S. Department of Energy through Contract No. DE-FC05-85ER250000.
REFERENCES 1. C. Canuto, M. Y. Hussaini, A. Quarteroni and T. A. Zang, Spectral Methods in Fluid Mechanics. Spdnger-Verlag,New York 0957). 2. M. D. Salas, T. A. Zang and M. Y. Hussaini, Shock fitted Euler solutions to shock-vortex interactions. In Proceedings of the 8th International Conference on Numerical Methods in Fluid Dynamics (Edited by E. Krause), pp. 461-467. Springer-Vedag, New York (1982). 3. D. Gottlieb, L. Lustman and C. $treett, Spectral methods for two-dimensional shocks. In Spectral Methodsfor Partial Differential Equations (Edited by Voigt, D. Gottlieb and M. Y. Hussaini), pp. 79-95. SIAM, Philadelphia, PA 0984). CAF 21/2--I
266
DAVID A. KOPRIVA
4. M. Y. Hussaini, D. A. Kopriva, M. D. Salas and T. A. Zang, Spectral methods for the Euler equations: Part l--Fourier methods and shock capturing. AIAA Jl 23, 64 (1985). 5. M. Y. Hussaini, D. A. Kopriva, M. D. Salas and T. A. Zang, Spectral methods for the Euler equations: Part II--~hebyshev methods and shock-fitting. AIAA Jl 23, 234 (1985). 6. W. Cai, D, Gottlieb and C. W. Shu, Essentially non-oscillatory spectral Fourier method for shock wave calculations. Maths Comput. 52, 389 (1989). 7. D. A. Kopriva, T. A. Zang and M. Y. Hussaini, Spectral methods for the Euler equations: the blunt body problem revisited. AIAA Jl 29, 1458 (1991). 8. D. A. Kopriva, A spectral multidomain method for the solution of hyperbolic systems. Appl. Numer Math. 2, 221 (1986). 9. D. A. Kopriva, Solution of hyperbolic equations on complicated domains with patched and overset Chebyshev grids. SIAM Jl Scient. Statist. Comput. 10, 120 (1989). 10. D. A. Kopriva, Multidomain spectral solution of the Euler gas-dynamics equations. J. Comput. Phys. 96, 428 (1991). 11. G. I. Taylor and J. W. Maccoll, The air pressure on a cone moving at high speeds. Proc. R. Soc. A139, 278 (1933). 12. G. Moretti and M. Padolfi, Entropy layers. Computers Fluids 1, 19 (1973). 13. G. Moretti, Thoughts and afterthoughts about shock computations. PIBAL Report 72-37, Polytechnic Institute of Brooklyn, New York (1972). 14. D. A. Kopriva, A multidomain spectral collocation approximation of the sound generated by a shock-vortex interaction. In Computational Acoustics: Algorithms and Applications (Edited by D. Lee et al.), pp. 377-388. Elsevier, Amsterdam (1988). 15. D. A. Kopriva and M, Y. Hussaini, Multidomain spectral solution of shock-turbulence interactions. In Domain Decomposition Methods (Edited by T. F. Chan, R. Glowinski, J. Periaux and O. B. Widlund), pp. 340-350. SIAM, Philadelphia, PA (1989). 16. M. J. Zucrow and J. D. Hoffmann, Gas Dynamics, Vol. II. Wiley, New York (1977). 17. T. A. Zang, M. Y. Hussaini and D. M. Bushnell, Numerical computations of turbulence amplification in shock wave interactions. Presented at the AIAA 20th Aerospace Sciences Mtg, Orlando, FL, AIAA Paper 82-0293 (1982).
SPECTRAL SOLUTION OVER WEDGES
0045-7930/92 $5.00+ 0.00 Copyright © 1992 Pergamon Press plc
OF INVISCID SUPERSONIC FLOWS AND AXISYMMETRIC CONES DAVID A. K O P R I V A
Department of Mathematics and Supercomputer Computations Research Institute,The Florida State University, B-186, Tallahassee, FL 32306-4052, U.S.A. (Received 1 March 1991; in revised form 16 July 1991; received for publication 24 October 1991) Abstract--We use a shock-fitted multidomain spectral collocation method to solve both steady and unsteady inviscid supersonic flows over bodies. New aspects of the method include two subdomain interface types and a zonal solution procedure to get efficient convergence to steady-state in supersonic regions. Steady-state examples include flow over a sharp cone, a hyperbolic cone and a hyperbolic wedge. For the cone, the exact flow solution is used to show that the method is spectrally accurate. As an example of an unsteady flow, we present a calculation of the interaction of a free-stream hot ring with the flow over a sharp cone.
1. I N T R O D U C T I O N The use of spectral collocation methods to solve multidimensional compressible flow problems is less than a decade old [1]. The earliest papers [e.g. 2, 3] show only the feasibility of the method and date from the early 1980s. Spectral methods were introduced with the hope that the high-order ("spectral") accuracy that proved to be important in the solution of many incompressible flows could be used to advantage to understand the physics of some compressible flow problems. Because they are somewhat new, techniques associated with the more established methods are only now being applied to spectral methods. When first applied to compressible flows, spectral methods had several inherent limitations. One problem was how to compute shock waves. Without special care, spectral methods reduce globally to low-order polynomial accuracy when shocks are present. To address this problem, both shock-capturing and shock-fitting methods have been studied [4-7]. A second limitation was the inflexibility of the methods. In multiple space dimensions, tensor product interpolations are used to define the global spectral approximation. This means that the domain of interest must be simple enough to be mapped onto a square. Given that mapping, the grid points must then be placed according to a specific, typically nonuniform, distribution. The result is that when steep gradients are present, large oscillations appear in the solution unless many grid points are used. To get past the second limitation, multidomain techniques (among others) were introduced in the mid 1980s [8-10]. In a multidomain technique, like the zonal techniques associated with finite difference methods, the physical domain of interest is divided into multiple subdomains. Local mappings of the subdomains allow the solution of flows on complicated domains. Also, the careful placement of subdomains, and the choice of the number of grid points within them, allows for local refinement of the grid. An additional limitation of spectral methods comes from the stiffness in time of the spatially discretized equations. This problem has yet to be resolved satisfactorily. The derivative matrices are full, so explicit time marching procedures are usually used. The time steps associated with the commonly used Chebyshev spectral methods vary as I/N 2 in one space dimension, where N is the number of grid points. This can be compared with the O(1/N) variation typically associated with finite difference methods. A multidomain discretization helps this somewhat [8, 10]. For instance, a 10-domain decomposition with 10 points in each can take a time step roughly 10 times larger than can be taken on a single domain with 100 points. Also, the work required to compute the derivatives in that multiple domain decomposition will be 1/10 that of the single domain. Except for the overheads associated with the interface procedures, the result is an overall factor of 100 advantage for the 10-domain decomposition. 247
248
DAVIDA. KOPRIVA
In this paper, we extend and apply the spectral multidomain technique of Ref. [10] to solve some inviscid supersonic flows over bodies. The extensions are needed because the flows are more complicated than those computed previously. Unlike the flows computed in Ref. [10], the flows considered here will be rotational and not homentropic. Also, a bow shock is present ahead of the bodies. This bow shock is fitted so that the flow is smooth in the computational region. The multidomain decomposition is used to resolve narrow layers in the solutions and to reduce the stiffness of the complete problem. To extend the method, we introduce two new interface conditions. The first condition arises at the intersection of the fitted shock boundary and an interface boundary. This allows the shock to be divided into multiple segments. We also introduce "shear interfaces" along which the grids slide as the calculation progresses. With shear interfaces, subdomains do not need to be fixed relative to each other. Finally, to make the calculation of steady-state flows more efficient, particularly those presented in this paper, we group subdomains into logical "zones". In supersonic regions of a flow, the solution is marched zone by zone from upstream to downstream. To apply the method, we consider steady and unsteady supersonic flows over two simple geometries. We first compute the flow over a sharp nosed cone (the "Taylor-Maccoll" [11] problem). An exact solution to this problem can be found, and we observe spectral accuracy in the computed solution. The second problem is that of flow over a hyperbolic wedge or cone. For long bodies, the significant computational difficulty is the presence of a very narrow layer near the body in which the entropy and vorticity vary strongly. The finite difference approximation of entropy layers was considered in Ref. [12]. We show that a multidomain discretization allows the spectral method to resolve the high gradients in these thin layers. Finally, we consider the interaction of a free-stream hot ring with the flow around a sharp cone as an example of a time-dependent flow. This paper is organized as follows. In Section 2 we present the equations and the mappings. We also derive the shock acceleration formula for time-dependent free-stream variations as an extension to the work in Ref. [7]. In Section 3, we review the multidomain method. We then introduce the new shock-interface condition, the shear interface condition and the zonal approach used to compute the steady-state solutions. Section 4 contains examples. Finally, some conclusions are drawn in Section 5. 2. T H E E Q U A T I O N S , T H E M A P P I N G S AND S H O C K F I T T I N G In this paper, we solve flows around bodies placed in a supersonic free stream. The two geometries considered are shown in Figs l(a, b). The first situation, Fig. 1(a), is used for the solution of the sharp cone problems. We avoid the numerical difficulties at the sharp nose and specify supersonic inflow conditions at a plane located downstream of the nose and perpendicular to the body. Figure l(b) represents the axisymmetric blunt body problem for which single domain spectral results were presented in Ref. [7]. In this problem, a portion of the flow between the body and the shock is subsonic. Also, entropy and vorticity are generated near the nose where the shock is strongly curved.
Shockw
~
M>I (
(b) Fig. 1. Diagrams of the computational regions. (a) Sharp cone geometry. (b) Blunt body geometry.
Inviscid supersonic flows over bodies
249
The inviscid flow between the bow shock and the body can be represented by the Euler gas-dynamics equations:
P,+UPx+VPy+ 7(Ux+Vy+~)
=0,
u, + UUx + vuy + "rPx = O, vt + UVx + vvy + "rPy = O, St+uSx+vSy=O.
(1)
Here, P = In(p) is the logarithm of the pressure. The velocity is q = (u, v) and r = a2/y is the temperature. The entropy, S, is scaled to Cv, and a is the sound speed. We take 7 = 1.4. The pressure and density are scaled to the free-stream values. Velocities are scaled to ~ . For the blunt nose problems, the lengths are scaled to the nose radius. The parameter ~ = 0 for the planar case, and ~ = 1 for the axisymmetric case. A transformation to computational coordinates is done in two steps: (x,y, t)~--~(r,s, t) (X, Y, T). The first mapping is to a coordinate system where s approximates the arclength along the body (with s = 0 at the nose), and r is the normal distance from the body. In (r, s) space, the region between the shock and the body is divided into multiple subregions (subdomains). Then these subdomains are mapped individually to the unit square, (X, Y)E [0, 1] x [0, 1]. The mapping (x, y ) ~ (r, s) is made by the relations x = Xb(S) + r~" ~, y = yb(S) + r~ .y; where (Xb, yb) represents the body as a function of s. As an example of this mapping, we consider the hyperboloid (or hyperbolic wedge) of Section 4. There Yb = tan(~b)x/(x + x0) 2 -- e2, where • = 1/tanZ(~b) and x0 = • + 1. This defines a hyperbola with unit radius of curvature at the nose and center at x = 0. Asymptotically, as x --* oo, the body is straight and with angle q~ to the horizontal• A reasonable parametrization of the distance along the hyperbola is in terms of y, so we define yb(S) = S, Xb(S) = .~/(s /m ) ~ + a z -- x0; where m = tan(q~) is the asymptotic slope of the body. After some algebra, we have the mapping r
x = xb(s) +
,/ l+\ay/ (ox y' Oxb
r ay
where CqXb
S
Oy
m2[xb(s)+ x0]
In (r, s) coordinates, we create multiple nonoverlapping subdomains, Gk, which collectively cover the region between the shock and the body (Fig. 2). Each subdomain is bounded by four curves. In r, the bounds are rmin(s, t, k) and rr~(S, t, k). Time-dependent subdomain boundary positions are included to allow the shock boundary to move with time. In a later paper, we will also consider the use of time periodic body motion to generate acoustic waves at the body surface. For convenience only, we require the subdomains to be bounded in s by straight lines normal to the body surface. CAF 2|/2--H
DAVID A. KOPRIVA
250
Point
Fig. 2. D i a g r a m of the d o m a i n decomposition. Each s u b d o m a i n , Gk, is b o u n d e d by four curves: r = rmi,(s, t, k), r = rmax(S, t, k), s = stain(k), s = smx(k ). T h e new interface points introduced in this paper are (c) and (d).
Given the boundary curves of a subdomain, a linear mapping from (r, s, t) ~ (X, Y, T) is used: r -- rmi.(s, t, k) rmax (S , t, k) - rmin(S ,
(2a) t, k ) '
S - - Stain ( k )
r~=
(2b)
Smax ( k ) - - Stain ( k ) '
T=t.
(2c)
Under the complete transformation from a subdomain to a square, (x, y, t) ~ (X, Y, T) [0, 1] × [0, 1], the equations become
Pr+ UPx+ VPy+ 7[Xxux+ X, vx+ Yxur+ Yyv~+~]=O, a 2
u r + Uux+ Vu~+-- (XxPx+ YxPr) = 0,
7
a 2
Vr + Uvx + Vvr + - - (XyPx + YyPr) = O, Sr + USx + VSr = O;
(3)
where
u = x, + UXx + VXy, V=uYx+vYy. The X, term appears because of the moving shock. Finally, the metric coefficients can be written in terms o f the intermediate metric terms as
Xx= XsSx + X, rx, Xy = Xssy + X~G,
Inviscid supersonic flows over bodies
251
and Yx = Yssx+ Y,r~, Yy = Yssy + Y, ry.
To ensure that the spectral approximation is not used to approximate derivatives across the shock, the shock is fitted as a boundary. The use of shock fitting with spectral methods is based on the work of Moretti [13]. It was shown to be viable in Refs [2, 5]. Kopriva et al. [7] have shown that the simplified shock-fitting procedure used in the early papers is not stable when used with spectral methods to compute steady flows. Instead, they derived a fully characteristic approximation. With the new shock-fitting formula, spectral accuracy of the stagnation pressure was shown for the blunt body problem. In this paper, we include the possibility of a time-varying free stream. So, we rederive the shock acceleration formula presented in Ref. [7] to account for free-stream variations. To fit the shock as a boundary, we set the shock position, denoted rsh(S, t), to be a boundary of the flow. Thus, for a subdomain that is bounded by the shock, rmax(S, t, k) = rsh(S, t). As the flow develops, the shock moves with velocity w = w~, where ~ is a unit vector in the r direction and w = Orsh/Ot is the shock speed. The shock speed is computed by integrating the shock acceleration in time. The shock position is updated by integrating the shock speed in time. To obtain the shock acceleration formula, we differentiate in time the Rankine-Hugoniot relations for the jump of the normal velocity and pressure across the shock. In what follows, let subscript 1 represent the low-pressure (known) quantities upstream of the shock; subscript 2 will represent the downstream or unknown flow quantities. For a time-varying free stream, the Rankine-Hugoniot relations are: ,n~----~--- TI)
(4a)
and ),-1 2a~ 1 ($2 = ~ ($1+ ~ 7+ 1 (st--'
(4b)
where ($i = q i - N - w(~-N) is the normal gas velocity relative to the shock and N = N ~ + Nz.9 is the shock normal. Since equations (4a, b) hold for all time, T, we differentiate them with respect to T to give P: = Pl + G31 + H
(5a)
(~2= F61 + E.
(5b)
and
Here, we use the dots to represent partial time derivatives, (') = O/OT. The coefficients E, F, G and H are related to quantities on the low-pressure (prescribed) side: E ~
-
27~ -
(~ + l)a,'
F=7- 1
27T1
~+l
(~ + 1)($~' 2($1
G=
and I H= -
7-- 1 2
a~
+ - -1
~-~ -- --~"-
71 "~1
252
DAVIDA. KOPRIVA
The derivative of thc relative velocity in terms of the shock accelaration, ~., is 42" N = w(~. N)(1 - F) + F(q,. ~ + ti," N) + w~. l~l(l - F) - q2" i'q + E.
(6)
The pressure variation is related to the shock acceleration by P2 = Pl + G[ql "~ + ¢h .N - w~.(~] - w~.NG + H.
(7)
To relate equations (6) and (7), we use the compatibility condition for the waves hitting the shock from the high-pressure side:
P2 - L~:. N + C = 0,
(8)
a2
where C = (qz - aN)" V P + Rp - a
(N,R. + NzR~),
V P = (XxPx+ YxPy).¢c + ( X y P x + YyPy).f,
Rp=,(Xxux+
Xyvx+ Yxuy+ Y y v r + ~ ) +
X, Px.
R. = Uux + Vuy and
R~ = Uvx + Vv y. Finally, equations (6) and (7) are substituted into equation (8) and the result is solved for the shock acceleration: -
aPl + (q~. Iq + dh'N ) (aG - ~F) - w ~. N(aG - 7F + 7) + ~ q2" N + a ( C + H) + yE ~. N(aG - 7F + ~')
(9)
The upstream quantities, P~, ~h and ~ used in equation (9) come from the time variation of the upstream flow. We define the flow upstream of the shock to be the uniform free stream plus a perturbation, e.g. P~ = P~ + eP'(x, y, t'), where the perturbation is defined in its own rest frame. Since the flow is supersonic, the perturbation is simply advected along with the flow until it interacts with the shock. In the rest frame of the shock, then, the time derivatives become
0,
PI
= E/-~TT.. + (q, + w ) ' V q '
LotF~q' L c3t'
]
1
(ql + w).VP' ,
and [-c3z'
e~ = , 1 _ ~ +
,
] (ql + w ) . W '
.
3. THE M U L T I D O M A I N S P E C T R A L M E T H O D 3.1. Discretization and interface conditions To discretize the system of equations (3), we start with the multidomain Chebyshev spectral method of Kopriva [10]. The method is general in the sense that there can be any number of subdomains defined in any order. A description of the data structures is given in Ref. [9]. Of practical importance is that local refinement of the grid is possible. Subdomains can be of any size. The number of grid points within a subdomain is independent of the number of points in its neighbors. Finally, subdomains can be overlapped, if necessary, but we do not use this feature in this paper.
Inviseid supersonic flows over bodies
253
Within each subdomain, the spatial derivatives in the equations are approximated by a standard Chebyshev spectral collocation method [1]. The result of the approximation is a system of ordinary differential equations for the unknowns at each grid point. The system is advanced in time by one of two time integration procedures. For the time-accurate calculation, we use the second-order midpoint rule. For unsteady problems, the low order of the time differencing procedure is not a serious problem. The stiffness of the discrete equations in time already requires the use of small time steps. For the steady-flow problems we are not interested in time accuracy. So, we use the first-order accurate four-stage Runge-Kutta method [1] commonly used in finite volume integrations. This allows a larger time step and more efficient time integration than does the midpoint rule, yet it has the same storage requirements. Numerical experiments show it is stable for a Courant number of at least 3, based on the smallest grid spacing. The approximation of the shock normals and their derivatives needed for equation (9) is done spectraUy. The two outward shock normal components can be written as /6~rsh
~/
/c~r~h
\ /
and
where
A
/(ar,h sx-rx )2
v
The shock slope, Or,h = -_r,,, -~ y,, ds dY is computed by the one-dimensional Chebyshev derivative of the shock position array. Thus, it is computed to the same accuracy as are the derivatives of the dependent flow quantities. Similarly, the time derivative of the shock normals requires the computation of dw/as, which is also computed spectrally. The derivation of interface conditions between subdomains for various types of interface points is described in detail in Ref. [1(3]. We only summarize those results here. The equations are solved on each subdomain to get the interior point solutions and a (provisional) value of the flow quantities on the boundaries. Each interface point has at least two provisional quantities, one from each subdomain to which it is a neighbor. These provisional quantities are then projected onto inflow and outflow components. To project the solutions, four generalized Riemann variables are computed from the infinite number represented by P + (~/a)q. ~ where g is a unit vector. The four choices are for ~ = + N and ~ = + X, where N and Z represent the normal and tangent vectors to the interface. To obtain the correct domain of dependence, the Riemann variables are then selected according to the subdomain from which their associated bicharacteristic comes. If the grid points do not match across the interface, the incoming quantities are interpolated using the natural Chebyshev interpolation that defines the method. (See Ref. [9] for details.) The resulting Riemann variables are then recombined to get corrected values of the flow quantities. As an example, we consider the computation of the flow quantities at the intersection point between the three domains GI, G2 and G3 shown in Fig. 3. In Ref. [10], this was called a "T" intersection. Superimposed on the grid, we show a projection of the bicharacteristic cone that affects the interface point. We represent the base of the bicharacteristic cone by the circle. Since the length of the fluid velocity vector, q, is shown as smaller than the sound speed represented by the radius of the circle, Fig. 3 shows a subsonic case. Four bicharacteristics have been chosen corresponding to ~ = ___~, ___)3.According to this choice, R + --/3 + (v/a)a and S + = 13 + (v/a)g, computed from the provisional values/3, a and g at the point in subdomain G3, and the sound speed evaluated at the previous time level. We compute S - from provisional values in GI. R - is computed from values interpolated from G2 to the interface point,
254
DAvID A. KoPRIVA
g-4--/ l q \ IIGI
F3!
! !!
Fig. 3. Determination of the domain of dependence in a subsonic flow at the intersection of three subdomains.
The interface condition requires that the true Riemann variables at the interface equal those computed from the provisional quantities [10]. Thus, we have four equations
p + 7-u=R+, a p -Tu=R-, a p + aY-v=S +, and
P - Ya- v = S - , from which to calculate the three quantities, P, u and v. Our choice of the four bicharacteristics makes the calculation of the velocities simple: a
u = Z (R + - R - )
and a
v =T~ (s÷ - s-).
The pressure is computed by combining all four Riemann variables and subtracting the value computed in the upstream subdomain, G3 (see Ref. [!0] for details): e =~(R + +n-
+s + +s-
- 2P3).
Finally, the entropy is advected along streamlines so its domain of dependence is given by the vector q shown in Fig. 3. Thus, the corrected value of the entropy is just the value computed in G3, i.e. S = ~ . To compute the flows considered here, we introduce two new types of interface points. These points are represented in Fig. 2. The first we call "shock-interface points" and they occur where the shock is divided at an interface. Many examples of such points are shown in Section 4. We remark that shock-fitted multidomain solutions have been presented in the past [14, 15]. See also Ref. [1, Chap. 13]. The subdomains in those papers were always strips along the shock direction. With the shock-fitting and interface procedures used in those papers, shock-interface points were not stable. At shock-interface points we use the shock normal and characteristic information from the upstream subdomain. The difficulty at shock-interface points is that there are two values of the shock normal, one from each subdomain, from which the correct one must be chosen. The shock acceleration formula (9) uses the compatibility equation (8) for the normal bicharacteristic. Since that bicharacteristic originates upstream of the shock-interface point, we use the normal and shock acceleration computed on the upstream side of the interface to determine the flow.
Inviseid supersonic flowsover bodies
255
We call the second type of interface a "shear interface" (Fig. 2). It represents a straightforward extension to the interface treatment in Ref. [10]. At a shear interface, neighboring subdomains are allowed to slide relative to each other as a function of time. A shear interface is convenient to use to solve the blunt body problem, for example. The narrow entropy layer that needs to be resolved near the body appears far downstream of, not near, the nose. The use of the shear interface allows for a change in the grid topology from one that does not need to resolve a narrow boundary layer to one that does. To define the shear interface, we use the hyperbolic nature of the problem to consider each subdomain as independent of the others. Then, the interface conditions are nothing more than boundary conditions on the incoming characteristic quantities. If the subdomains do not move relative to each other, the grid connections and interpolations for these characteristic quantities are defined once, at the start of the computation [9]. For the shear interface, the connections and interpolations are updated at each time step. The shear interface is the most expensive of the interface types and is used sparingly. The main expense is in the calculation of the Lagrange interpolating coefficients as described in Ref. [9]. To make this calculation as fast as possible, we no longer compute these coefficients directly using the definition of the Chebyshev polynomial, T, (x) = cos[n cos- l(x)]. They are now defined recursively by the standard three-term recursion for the Chebyshev polynomials [1]. 3.2. Acceleration to steady state To decrease the computer time needed to calculate the steady-state solutions, we introduce the notion of zones that are groups of subdomains to be computed together. This notion is convenient for the problems computed here since the solutions develop from the nose to the downstream and most of the flow is supersonic. In the supersonic portions of the flow, interface conditions are only needed on the upstream side of a subdomain. If the solution in a subdomain upstream of a given subdomain is converged to an acceptable level, there is no need to continue to compute on the upstream subdomain. In other words, subdomains within upstream "zones" can be turned "off" and not be computed further. Finally, subdomains in downstream zones do not have to be computed until the zones upstream of them are converged. Zones are defined in increasing numerical order from upstream to downstream. They define a purely logical grouping device that tells which subdomains are to be computed simultaneously. The interface conditions are not affected in any way. Note that all subdomains that contain subsonic points must be grouped in a single zone and solved together to convergence. For the calculations that we present here, the subsonic region is only a small fraction of the total flow. There are several advantages to be gained by solving by zones. First, only a small fraction of the total number of grid points needs to be solved at a given time. Next, the time step is regulated only by the subdomains that are turned on. If all subdomains are computed simultaneously, the time step must be determined by the minimum over all the subdomains. The minimum is usually the time step at the stagnation point. It can be an order of magnitude smaller than the local time step in the supersonic parts of the flow. Finally, the concept of zones helps to avoid the need for a highly accurate initial shock shape. For long bodies, a small error in the initial shock shape near the nose can translate into a large error far downstream. This has typically lead to a kink instability in the shock as the computation proceeds. To avoid the problem, the solution in a supersonic zone is "bootstrapped" from the previous zone by linear extrapolation. Far downstream, the shock is nearly straight so this is a reasonable initial approximation. 4. EXAMPLES We now present three examples of the use of the multidomain method described in the previous section. The first is the Taylor-Maccoll problem [11] for which an exact solution can be computed. The second is the steady flow over a blunt hyperbolic wedge and cone. The unique feature of the blunt body flow solutions is the narrow entropy and vorticity layer that develops downstream of the nose. Finally, we include a time dependent calculation that models the interaction of a free stream hot ring with the flow around the sharp nosed cone.
256
DAVIDA. KOPRIVA
Fig. 4. Geometry and grid topology for the solution of the Taylor-Maccoll [11] problem.
4.1. The Taylor-Maccoll [11]problem The first problem we consider is the steady flow over a sharp nosed cone at zero angle of attack. We assume that the semi-vertex angle is small enough so that the shock is attached to the cone. Then, the shock is straight and the flow solution is conical. Also, the entropy is constant and the vorticity is zero in the region between the shock and the cone. To find the exact solution to this problem, we solve a system of ordinary differential equations [16]. A Richardson error estimate is used to ensure that enough grid points are used to obtain machine accuracy. To study the error, the problem was solved on the three subdomain grid topology shown in Fig. 4 for several values of the free-stream Mach number and semi-vertex angle. The grid was chosen to show the accuracy of the method when shock-interface points are used. The upper boundary corresponds to the fitted shock wave and the lower boundary is the body. Along the body, the usual charateristic boundary conditions were used [cf. 7]. The inflow boundary at the lower left is supersonic, so all flow variables were prescribed from the exact solution computed at the inflow points. The outflow is also supersonic so no boundary conditions were applied there. As the initial condition, we specified the exact solution. The solution was advanced in time until the shock speed converged. Spectral accuracy was observed for the Taylor-MaccoU [11] problem. As an example, we present the results for the M = 6 flow over a 10° cone. Figure 5 shows the maximum errors of the shock slope, the pressure, the horizontal and the vertical velocity as a function of the number of grid points in the direction normal to the body. The solution is analytic, so we expect exponential decay in the error [1]. All quantities do show exponential decay in the error.
4.2. The blunt body problem We consider next the solution of flows over hyperbolic wedges and cones. The (single domain) spectral solution of similar blunt body problems was described in Ref. [7]. Here we are most interested in the flow over bodies that are long compared with the nose radius. The computational difficulty is to resolve the narrow entropy and vorticity layers that form along the body [12]. To
o
_4,0 ~
,~: ~ ' ~
"~
-~-Slope I ~Pressurel
~
uJ
E
-s.o
E
I I I
~
-6.0 o 0 -J
-7.0 -8.0 ~
3
4
~
5
6 N
....
7
8
J
9
Fig. 5. Error decay of the shock slope, pressure and velocity components for the Taylor-Maccoll [11] problem.
Inviscid supersonic flows over bodies
257
resolve these layers, we place narrow subdomains along the body. We remark that this is similar in spirit to numerical approaches suggested in Ref. [12]. To start the calculations, a reasonable initial solution was used. A parabolic shock shape was assumed near the nose. The Rankine-Hugoniot relations were used to specify the flow variables along the shock• Along the body, the initial condition assumed that the Math number varied as sin(s) up to the sonic point and constant beyond that. The sonic point was placed along the body at a point near the nose. The interior flow variables were then defined by linear interpolation. In all the cases presented in this section, the solutions were converged to maximum shock speeds between 10 -5 and 10 -4 depending on the free-stream Mach number. Solutions were computed for a wide range of parameters. Asymptotic body angles ranged from 10° to 30° for both wedge and cone geometries. Mach numbers varied between 4 and 6. Finally, body lengths ranged from about 2 to 225 nose radii. When possible, solutions were compared with the single domain calculations using the single domain spectral code described in Ref. [7]. As an example, we first present a solution in the subsonic portion of the M = 4 flow over a 10° hyperbolic wedge on a grid with three subdomains. Figure 6 shows the vortieity contours and grid for the initial and converged solutions. The subdomains were chosen so that a shear interface and shock-interface point would occur within the subsonic portion of the flow. Notice that the grid lines have changed relative to each other as the solution converged.
rid
I n t t• "
Conv
Fig. 6. Vorticity contours and grid for the solution near the nose of the M = 4 flow over a hyperbolic wedge. The sonic line is indicated with the vorticity contours.
258
DAVID A. KOPRIVA
A severe test for any numerical method, including this one, is to see how smooth the vorticity contours are. Vorticity is a quantity computed after the fact from the derivatives of the velocities. Thus, the errors in the vorticity will be larger than the errors in the velocity components alone. To compute the vorticity, we use the usual Chebyshev differentiation of the velocities. The results of this calculation are the contours shown in Fig. 6. The contours were not smoothed by the plotting routine. Thus, we see that the multidomain method gives a solution that is not only continuous, but smooth through the interfaces and throughout the flow. As a final measure of the quality of the solution, we note that the vorticity along the wall, which should be 0, is at most 4 × 10 -3. The shape of the vorticity layer further down the body is determined by the distribution of the shock-generated vorticity in the nose region shown in Fig. 6. Since the shock becomes nearly straight in the neighborhood of the symmetry line, the vorticity approaches zero. In that same region, the entropy is nearly constant. Then Kelvin's theorem applies and the vorticity lines will follow the streamlines. Away from the symmetry line, the shock becomes strongly curved, the entropy varies strongly (cf. Fig. 7) and the vorticity generation reaches its maximum. Finally, away from the nose, the shock straightens out and the vorticity again approaches zero. Note, that, in the region of the curved shock and strongly varying entropy, Kelvin's theorem does not apply and
Pressure
Mach
Entropy
Fig. 7(a)
Inviscid supersonic flows over bodies
Fig.
259
7(b)
Fig. 7. The grid contours of the pressure, M a t h number, entropy and vorticity for the steady M = 6 flow over a 25 ° hyperbolic wedge.
the vorticity contours do not follow the streamlines. Thus, we can predict that profile between the body and the shock along a normal to the body of the vorticity layer far downstream of the nose. The vorticity will rapidly drop from zero at the body to a minimum and then more slowly rise to zero at the shock. Figure 7 shows the solution contours and grid for an M = 6 flow over a 50 nose radius long 25 ° hyperbolic wedge. Eight subdomains and a total of 800 grid points were used. The grid included both a shear interface and several shock-interface points in the supersonic portion of the flow. A narrow subdomain layer was used along the body outside the nose region to ensure uniform resolution of the entropy and vorticity layer. Figure 8 shows the variation of the solutions along the normal direction at the outflow boundary of the calculation. The subdomain interface appears at the distance r = 2 from the body. The pressure varies by < 2% across the flow so it is clearly over-resolved by this grid. More structure appears in the velocity [=(u s +v2)~/2], the entropy and the vorticity. For instance, the shock slope is zero at the symmetry axis. So, the derivative of the velocity and the entropy should be zero at the body. This feature is captured by the method. Also, the vorticity is smooth and has the structure suggested by the variation of the vorticity generation near the nose. It is zero at the wall and ends at zero at the shock. In this calculation the vorticity error along the wall is only 6 x 10 -3.
DAVID A. KOPRIVA
260
14
14
. . . .
I
. . . .
I
. . . .
I
. . . .
I
12
10
10
8
8
6
6
4
4
2
2
0 13.15
0 0.5
14
.
13.2
.
.
.
I
.
.
13.25 13.3 Pressure
1 3 . 3 5 13.4
.
I
.
I
.
.
.
.
.
.
.
4
.
(b)
12
~ - - ~ i'
(a)
12
,
. . . .
0.7
P
12
t
*
I
. . . .
0.9
1.1 Entropy
F
. . . .
~
. . . .
I
. . . .
I
1.3
. . . .
~
. . . .
.5
r
. . . .
v '
''
(d)
10
10 '
8
8
6 4
4
2
2
0
'
4
w
'
4.5
5 Velocity
5.5
0
6
'
'
-
-0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 Vorticity
0
0.1
Fig. 8. The variation of the solutions along the outflowplane from Fig. 7. A more difficult calculation is presented in Figs 9 and 10 which show the M = 4 flow over a 10° hyperbolic cone of length 100. This time, 18 subdomains and over 1600 grid points were used. As before, the calculation includes a shear interface and multiple shock-interface points. Compared to the distance between the shock and the body, the entropy and vorticity layers are much narrower than in the previous example. Figure 11 shows the variation of the solution at the outflow plane of the calculation. Three subdomain divisions were used in the direction normal to the body to resolve the structure of the entropy layer. As before, the pressure variation shows little of interest near the body. The velocity, entropy and vorticity, on the other hand, show strong gradients there. Note that, even with this grid, we have not resolved all the structure of the layer. Missing from the velocity and entropy curves are the turnover to zero slope at the body. The vorticity graph has missed the zero vorticity along the wall. These features were missed in part due to the lack of refinement of the grid near the nose. What we see in this calculation, however, is the ability of the method to allow the computation and resolution of very steep, though continuous, features. 4.3. The interaction o f a hot ring with the flow over the sharp nosed cone
The interface conditions, being characteristic, are also applicable to time-dependent problems, particularly in the important acoustic regime. An example of a linear problem was considered in Ref. [10]. For nonhomentropic flows, however, the sound speed used to compute the Riemann variables is linearized in time about the previous time level. As described in Ref. [10], this linearization introduces errors of the same order as the time stepping scheme used. The interface conditions are not applicable, again because they are of characteristic type, to the passage of strong nonlinear waves with unresolved steep gradients such as shocks. In these situations, numerical experiments show that the waves are stopped at the interfaces.
Inviscid supersonic flows over bodies
261
Grid
Fig. 9. Solution contours and grid for the M = 4 flow over a 10° hyperbolic cone.
Our final example uses the method to compute a time-dependent rotational and nonhomentropic flow. We present a calculation of the interaction of a hot ring with the M = 4 flow over a 30 ° cone. Figure 12 shows the geometry and grid. The section of the cone ranged from s = 1.0 to s = 2.0. Eight subdomains and a total of 1736 grid points were used. To begin, we used the cone geometry of Section 4.1 and established a steady state. Then the temperature in the free stream ahead of the shock was modified to r = 1 + E exp( -~t[(x - x0) 2 + (y - y0)2]}. The temperature ring was allowed to advect with the free stream until it encountered the shock wave. This problem lacks an exact solution. However, the numerical solution of a similar s h o c k temperature spot interaction problem was obtained in Ref. [17] with a shock-fitted finite difference scheme. It was solved again in Ref. [15] by a multidomain spectral method that could not split the shock by a subdomain interface. Neither paper included the complications of a wall boundary. Both computations showed the generation of a sound wave, an entropy spot and a pair of counterrotating vortices. The shock-fitted spectral method described here is suitable for the study of the interaction of the shock with low-amplitude free-stream perturbations. High-amplitude perturbations create large variations in the shock shape and the possibility of additional shocks that cannot be captured by the spectral method. Therefore, calculations were performed with amplitude, e, between 0.05 and 0.2. The temperature perturbation was centered on (x0, Y0) --- (0.9, 1.05). The coefficient ~t = 50 was chosen so that at t = 0, the flow just ahead of the shock was consistent with the initially straight shock. Also, with these parameters, the spot was localized enough not to affect the shock at the inflow shock point on the upper left of Fig. 12.
262
DAVID A. KOPRIVA
Pressure
Mach
Fig. 10. Blowup of the nose region of Fig. 9.
The results show the same qualitative features as the previous calculations in Refs [15, 17]. The time history of the result of the interaction for E = 0.05 is shown in Figs 13-15. First, in Fig. 13, we see the generation of a sound wave that strikes the body and moves upstream along it. Figure 14 shows the generation of an entropy spot at the shock. This spot, which is flattened in the direction normal to the shock, breaks away and moves upstream along the streamlines. For our purposes here, we are most interested to see that the entropy is advected smoothly through the interfaces and with no visible dissipation or dispersion. Dissipation would show up as a decay in the maximum entropy and an increase in the size of the core of the spot. Over the 600 time steps taken between t = 0.156 to 0.250, the maximum entropy varies by at most 0.14%. Dispersion effects would be seen as waves trailing (or possibly leading) the entropy spot and are not visible in the plots. The time history of the vorticity is shown in Fig. 15. We see the generation of the counter-rotating vortices represented by positive and negative spots of vorticity. In this nearly linear regime, they are of nearly equal strength. The worst difference in the maximum and minimum vorticity is 1.3% at time 0.188. The cores of the vortex contours do grow with time, so we see that the vortices show the effect of more dissipation than does the entropy. However, they do not show any significant trailing waves that are characteristic of large phase errors, and they propagate smoothly through the interfaces with no visible reflection.
Inviscid supersonic flows over bodies 0
. . . .
I
. . . .
L
. . . .
I
. . . .
I
0
. . . .
_(a)
16
....
, ....
, ....
,
,
. , ....
, • . ,. , ....
(c)
16
12
263
12 t..
8
8
4
0 1.7
1.75
1.8
1.85
1.9
.95
0
0.1
0.2
0.3
0
20
_(b)
16
0.4
0.5
0.6
0.7
0.8
Entropy
Pressure
. ,
. ,
• ,
.
(d)
16
12
. ,
12 t..
8
8
4
4
0
, 3.8
'
'o~ 4
'
~- ' 4.2
^ ~ 4.4
Velocity
0
4.6
4.8
~
-8
.
.
~
-6
.
.
.
I
.
-
-4
.
t
-2
Vorficity
Fig. 1I. The variation of the solutions along the outflow plane from Fig. 9.
Fig. 12. The grid used for the calculation of the hot-ring calculation at the initial time.
~
,
~
.
0
.
•
264
DAVID A. KOPRIVA
TIME = 9.387E-02
= •
TIME =
.125
=.
TIME =
.156
=.
Fig. 13. Pressure contours for the hot-ring calculation.
TIME = .188
T
m
-
•
Fig. 14. Entropy contours for the hot-ring calculation.
Inviscid supersonic flows over bodies
265
"1"1
i Fig. 15. Vorticity contours for the hot-ring calculation. The solid lines represent positive vorticity. The dashed lines indicate negative vorticity.
5. C O N C L U S I O N S We have used a multidomain spectral method to solve some inviscid supersonic flows over some simple bodies. Unlike the solutions previously presented in Ref. [10], we have computed flows that have entropy and vorticity variations. Two interface procedures were introduced. One allowed the fitted shock to be intersected by subdomain interfaces. The other allowed interfaces to slide along each other. To accelerate the convergence to steady state and to use locally determined (by subdomain) time steps, subdomains were grouped into zones and marched in time by zone. The solutions obtained by the method were accurate and smooth. First, we solved the Taylor-Maccoll [1 l] problem. It has an exact, analytic solution. As expected, the errors decayed exponentially fast. The more complex problems such as the flow over the hyperbolic bodies had entropy and vorticity generated by the curved shock. By careful choice of the subdomain decomposition, it was possible to resolve steep gradients in these quantities. Though the interface conditions specify only the continuity of the characteristic variables, we see that even the vorticity is smooth across subdomain interfaces. The time-dependent calculation showed the generation of a counter-rotating vortex pair and an entropy ring. These features also propagated smoothly through the subdomain interfaces. Acknowledgements--The author would like to thank Drs C. L. Streett and T. A. Zang for helpful discussions. This research was supported in part by the National Aeronautics and the Space Administration under Grant No. NAGI-862 and by the U.S. Department of Energy through Contract No. DE-FC05-85ER250000.
REFERENCES 1. C. Canuto, M. Y. Hussaini, A. Quarteroni and T. A. Zang, Spectral Methods in Fluid Mechanics. Spdnger-Verlag,New York 0957). 2. M. D. Salas, T. A. Zang and M. Y. Hussaini, Shock fitted Euler solutions to shock-vortex interactions. In Proceedings of the 8th International Conference on Numerical Methods in Fluid Dynamics (Edited by E. Krause), pp. 461-467. Springer-Vedag, New York (1982). 3. D. Gottlieb, L. Lustman and C. $treett, Spectral methods for two-dimensional shocks. In Spectral Methodsfor Partial Differential Equations (Edited by Voigt, D. Gottlieb and M. Y. Hussaini), pp. 79-95. SIAM, Philadelphia, PA 0984). CAF 21/2--I
266
DAVID A. KOPRIVA
4. M. Y. Hussaini, D. A. Kopriva, M. D. Salas and T. A. Zang, Spectral methods for the Euler equations: Part l--Fourier methods and shock capturing. AIAA Jl 23, 64 (1985). 5. M. Y. Hussaini, D. A. Kopriva, M. D. Salas and T. A. Zang, Spectral methods for the Euler equations: Part II--~hebyshev methods and shock-fitting. AIAA Jl 23, 234 (1985). 6. W. Cai, D, Gottlieb and C. W. Shu, Essentially non-oscillatory spectral Fourier method for shock wave calculations. Maths Comput. 52, 389 (1989). 7. D. A. Kopriva, T. A. Zang and M. Y. Hussaini, Spectral methods for the Euler equations: the blunt body problem revisited. AIAA Jl 29, 1458 (1991). 8. D. A. Kopriva, A spectral multidomain method for the solution of hyperbolic systems. Appl. Numer Math. 2, 221 (1986). 9. D. A. Kopriva, Solution of hyperbolic equations on complicated domains with patched and overset Chebyshev grids. SIAM Jl Scient. Statist. Comput. 10, 120 (1989). 10. D. A. Kopriva, Multidomain spectral solution of the Euler gas-dynamics equations. J. Comput. Phys. 96, 428 (1991). 11. G. I. Taylor and J. W. Maccoll, The air pressure on a cone moving at high speeds. Proc. R. Soc. A139, 278 (1933). 12. G. Moretti and M. Padolfi, Entropy layers. Computers Fluids 1, 19 (1973). 13. G. Moretti, Thoughts and afterthoughts about shock computations. PIBAL Report 72-37, Polytechnic Institute of Brooklyn, New York (1972). 14. D. A. Kopriva, A multidomain spectral collocation approximation of the sound generated by a shock-vortex interaction. In Computational Acoustics: Algorithms and Applications (Edited by D. Lee et al.), pp. 377-388. Elsevier, Amsterdam (1988). 15. D. A. Kopriva and M, Y. Hussaini, Multidomain spectral solution of shock-turbulence interactions. In Domain Decomposition Methods (Edited by T. F. Chan, R. Glowinski, J. Periaux and O. B. Widlund), pp. 340-350. SIAM, Philadelphia, PA (1989). 16. M. J. Zucrow and J. D. Hoffmann, Gas Dynamics, Vol. II. Wiley, New York (1977). 17. T. A. Zang, M. Y. Hussaini and D. M. Bushnell, Numerical computations of turbulence amplification in shock wave interactions. Presented at the AIAA 20th Aerospace Sciences Mtg, Orlando, FL, AIAA Paper 82-0293 (1982).