CC BY
62
Вычислительные технологии
Том 6, № 6, 2001
STABILITY ANALYSIS OF SUPERSONIC ENTROPY LAYERS INCLUDING SHOCK EFFECTS
St. Mählmann Aerodynamisches Institut Aachen Rheinisch-Westfälische Technische Hochschule Aachen, Germany e-mail: stefan@aia.rwth-aachen.de
Локализация перехода ламинарного режима в турбулентный в сверхзвуковых пограничных слоях необходима для точного расчета аэродинамических сил и тепловых режимов при аэродинамическом проектировании гиперзвуковых летательных аппаратов и при управлении ими. Характеристики устойчивости сверхзвуковых пограничных слоев зависят, например, от параметров затупления носовой части, поперечной кривизны, температуры стенок, ударных волн и т. д. Эти параметры могут быть теоретически исследованы с помощью численного моделирования условной устойчивости с заданием на удаленной границе либо нулевых краевых условий, либо условий, описывающих асимптотические возмущения. В этом подходе формирование ударной волны перед передним фронтом затупленного тела игнорируется. Тем не менее для улучшения понимания взаимодействия между волнами неустойчивости, возникающими внутри сверхзвукового пограничного слоя, с волнами, идущими от невязкого энтропийного слоя, наличие этой ударной волны должно приниматься во внимание. В работе представлен метод моделирования устойчивости с физически согласованным включением разрывов. Внешние условия на свободный поток заменяются соответствующими условиями на разрыве. Необходимые уравнения возмущений могут быть выведены из из линеаризованных нестационарных уравнений Рэнкина — Гюгонио с учетом эффекта ударных осцилляций, образующихся вследствие волн возмущения, которые возникают из поля течения с наветренной стороны разрыва.
Introduction
The physical understanding of the stability properties of supersonic boundary layers is based mainly on results obtained with the linear stability theory. In most stability computations of supersonic boundary layers the presence of a shock-wave is ignored and vanishing perturbation conditions are imposed at the far field or the asymptotic form outside the boundary layer. This approach simplifies the upper boundary conditions and allows one to focus on the flow stability due to the viscous boundary layer. However, for supersonic flows over bodies with blunt leading edges, where the inviscid entropy layer between the boundary layer and the shock wave has a significant impact on the global stability behaviour of the flow [1], the free-stream conditions should be replaced by appropriate shock conditions.
© St. Mählmann, 2001.
In stability analysis of supersonic flows the interaction between the shock-wave and instability waves originating inside the boundary-layer and free-stream disturbances has to be taken into account [2]. Furthermore, it is well known that any kind of modes, for instance, acoustic, vorticity or entropy modes, which can generally be found in any real supersonic free stream, produces all three modes when passing the shock wave [3]. One of the first attempts to understand the effect of a shock wave on the boundary-layer stability was Petrov's work [4]. In his study the inviscid eigensolution outside the boundary layer was replaced by the linearized Rankine — Hugoniot jump conditions for the normal momentum equation. By coupling this boundary conditions with asymptotic conditions associated with the viscous eigensolutions in the free stream, he obtained solutions for a two-dimensional flow over a wegde in the hypersonic limit.
We are interested in the investigation of shock effects on the stability of the windward supersonic flow field over blunt flat plates. A set of perturbation equations is derived from the unsteady Rankine — Hugoniot conditions by accounting for the effect of shock velocity due to perturbated waves which originate behind the shock. These perturbation equations are then applied as outer boundary conditions for the quasi-parallel linearized interior stability equations [2].
1. Formulation 1.1. Equations
The compressible flow over a blunt flat plate confined between two boundaries located at the wall (y* = 0) and the shock shape (y* = L*) is analysied for linear stability. In a system of cartesian coordinates the streamwise, wall-normal and spanwise directions are represented by x*, y* and z*, where the superscript * denotes dimensional quantities. A perfect Newtonian gas is assumed. The three dimensional Navier — Stokes equations are
dq* d*
+ qVq*
= -Vp* + V[A*(Vq*)I + ^*(Vq* + Vq*tr), dp*
* * P cp
dT* ~3t*
dt* + q * VT
+ V(p * q *) = 0,
dp*
= V(k * VT*) + dt* + q* Vq* + $ *, = p*R*T *,
(1) (2)
(3)
(4)
where u * is the velocity vector [u,v,w]tr, p * is the density, p* is the pressure, T * is the temperature, R * is the gas constant, cp* is the specific heat at constant pressure, k * the thermal conductivity, ^* the first coefficient of viscosity, and A* is the second coefficent of viscosity. The viscous dissipation function $* is given as
$ *
tri 2
A *(Vu *)2 + y[Vu * + Vu ]
(5)
The flow variables and equations are nondimensionalized by the corresponding values on the boundary-layer edge: velocities by u*, length scales by I* = \Jv*s*/u*, where s* is streamwise coordinate in the body-oriented coordinate system with is origin in the stagnation point, the density by p*, temperature by Te*, pressure by p*eu*2, and the time scale by l*/u*. The Reynolds
*
p
*
number Re1b is defined as
^e *
and the Mach number Mae is
U *
* * 7 *
Re1b = -—*-> (6)
Mae = , e , (7)
where R* = the specific heat at contant volume, and 7 is the ratio of specific
heats. The Prandtl number is defined as Pr = y*c*p/k*. The viscosity coefficient is determined by Sutherland's law,
,i.5. 1 + c
" = yr + c), (8)
where C is a constant. In this paper we set C = 0.375, A = -2/3y, 7 = 1.4, and Pr = 0.715.
1.2. Basic-flow solution
The full Navier — Stokes equations for a two-dimensional, compressible flow are solved using the Advection Upstream Splitting Method (AUSMDV) of Wada and Liou [5] with a modified Limiter of Van Albada [6]. A second order accurray in space is obtained with the MUSCL interpolation of Van Leer [7]. Earlier studies [8] have shown that the solutions for the basic flow are sufficiently accurate to allow reliable instability analysis.
1.3. Linear stability equations
The stability investigation of the laminar basic flow is based on the a normal mode analyis of the linearized perturbationd equations of the three-dimensional Navier — Stokes equations. The governing stability equations are derived be representing as the sum of a mean and a fluctuation quantity, i. e.
u = U7(y) + u'(x, y, z, t), v = V(y) + v' (x, y, z, t), w = w'(x, y, z, t), T = T + T '(x,y,z,t),
p = p + p'(x,y,z,t). (9)
Substituting Eq. (9) into the nondimensional form of the governing Eqs. (1) - (5), and dropping the nonlinear and high-order terms yields to a set of linear differential equations for the perturbations equations. Further details of the linear perturbation equations and other formulations are described in Malik [9]. In the normal mode analysis for linear disturbance scenarios, the fluctuations of the flow quantities are assumed as harmonic waves of the following form:
[u',v' ,w',T',p']tr = [U(y),v(y),W(y),T(y),p(y)]tr (10)
where a and ft are the streamwise and spanwise wavenumbers and u denotes the frequency of the disturbance wave. Substituting Eq. (10) into the linearized perturbation equations leads (to a system of ordinary differential equations:
(AD2 + BD + C)$ = 0, (11)
where D is the differential operator in the wall-normal direction, i.e. D = d/dy,D2 = d2/dy2. In Eq. (11), $ is a vector defined as
$=[u(y),V(y),W(y),T(y),p(y)]tr. (12)
A complete listing of the 5 x 5 matrices A, B and C can be found in [9].
1.4. Boundary conditions for linear stability analysis
The governing stability equations (11) are solved as a boundary value problem. Therefore appropriate boundary condition have to be formulated. At the body surface no-slip conditions are imposed
y* = 0 : $1 = 0, $2 = 0, $3 = 0, $4 = 0. (13)
Here temperature fluctuations are assumed to vanish at the solid wall boundary. This is a reasonable approach for high frequency disturbances where the temperature fluctuactions will not penetrate deep into the solid boundary due to the thermal inertia of the wall. Therefore the wall will appear insulated on the time scale of the mean flow but not on the short time scales of the disturbances. Additionally, Neumann conditions for the pressure eigenfunction are enforced as
dp) dy
dp)
= a 7T
y*=0 dy
= b, (14)
y* =L*
where a and b are evaluated at the two boundaries using the normal momentum equation.
At the bow shock location y = L , shock conditions must hold. To understand the interaction between an oscillating shock wave and instability modes, which originate from the flow field behind the shock, the unsteady motion of the shock wave has to be taken into account. Under the assumption of an inviscid shock nature this unsteady motion can be fully described by the Euler equations. For a given shock position y0 = f (x,z,t) with the time averaged shock position y0 = f 0 and a local shock slope a = df /dx the jump conditions arcoss the shock are
|[Q] + f i*]- F1 + fiGl = °- (15>
The vectors Q, E, F and G are defined as
Q = [p, pu, pv, pw,e]tr,
E = [pu, pu2, puv, puw, (e + p]u)tr,
F = [pv, puv, pv2, pvw, (e + p]v)tr,
G = [pw, puw, pvw, pw2, (e + p)w]tr, (16)
where the total energy e is given by the relation
p , 1 ..2 , „,2
- p(u2 + v2). (17)
Y — 1 2
We note that Eq. (15) represents the unsteady Rankine — Hugoniot relation which governs the motion of a shock wave. For the blunt flat plate we consider the mean flow as "quasi-parallel". Under the assumption of small disturbance amplitudes the time dependent position function of the shock can now be perturbated according to f = f + f'. Although only perturbations behind the shock are considered in our study, Eq. (15) is in general valid for any small disturbancies
across the shock. After introducing harmonic waves into equation (15) we finally get the relations of eigenfunctions for normal mode analysis
i(a[E] + $[G] - w[Qj) + a[E] - [F] = 0,
(18)
with the vectors E and F
E
P2 U + pU2 PU22 + 2p2 U2 U + p
PU2V2 + P2UV2 + P2U2U PU2W2 + P2U2W2 + P2U2W (e2 + P2)U + U2(e + p)
P2V + PV2
PU2V2 + P2UV2 + P2U2V
Pv2 + 2P2UV2 + p
PV2W2 + P2V2W2 + P2V2W
(e2 + P2)V + V2(e + p5)
(19)
where the quantity e can be derived from
Y - 1
+
U22 + V22
P + P2(U2 U + V2V).
(20)
2
2. Numerical methods and validation
For a given set of real-value a and j3, the temporal stability analysis leads to u and $ as eigenvalues and eigenfunctions of the boundary value problem. Various numerical methods for solving the eigenvalue problem were discussed and tested by Malik [9]. The eigenvalue problem can by solved by either global of local methods. In this paper the eigen spectra of the discretized eigenvalue problem is computed with a global method determines yieds all eigenvalues together the associated eigenvectors, taking into bargain the computational more expensive effort compared to the local methods. To approximate the derivatives of the eigenfunctions in the wall-normal direction, we use a Chebychev spectral collocation method. The numerical procedure will be descibed in the this section
2.1. Chebychev spectral collocation method
The N-th order Chebychev polynominals TN are defined on the interval G [-1,1], where denotes the Gauss — Lobatto collocation points
0 = cos( j, J = 0,1, N. (21)
In order to apply the spectral collocation method, an interpolation polynominal is constructed for the dependent variables in terms of their values at the collocations points. A N-order polynominal is
N
) = E Afc (£)$(£fc ), (22)
k=0
where the interpolants \k (£) for the Chebychev scheme is
Ak(Î) = (-!)<«> (i-f ) Ng (23)
where c0 = cN = 2, and ck = 1 for 0 < k < N. From Eq. (22) the first derivative of $(£) can be written as
d$ N
dC
, (24)
j k=0
where Eij are the elements of the derivatives matrix defined as:
Ejj = j , for j = k, (25)
ck — ^k
Ejj =
C
j
E00 = - ENN
2(1 - CJ)'
2N2 + 1
6
With the scaling factor
Sj = f j • j = 0- (26)
the first derivative matrix F in the physical domain may be written as
Fij = Sj Eij, (27)
and the second derivative matrix Gij is
Gij FjmFmk■ (28)
Now the governing system of stability equations can be formulated at the collocation points as
NN
Aj J] Gjk$k + B^ Fjk$k + Cj= 0. (29)
k=0 k=0
To cluster grid points near the wall, algebraic and exponential stretching functions may be employed to transform the physical domain 0 < y < ymax to the computational domain —1 < C < 1. However, algebraic stretchings are more robust for spectral methods. We employ the stretching
y = a 1—1, (30)
where a = yiymax/(ymax — 2yi) and b = 1 + 2a/ymax. Here yi is the location corresponding to C = 0, i.e. half the grid points are clustered in the region 0 < y < yi, and ymax denotes the upper boundary location.
The discretized formulation of the interior stability equations (29) together with the boundary conditions (13), (14), (18) leads to a matrix eigenvalue problem, which is solved numerically with the QZ matrix eigenvalue algorithm.
2.2. Validation of the stability solver
The stability code was validated first by comparing the eigenvalue solutions of a linear stability analysis of a flat-plate compressible boundary-layer with zero pressure gradient with those of Malik [9] and Hu and Zhong [10]. For this test the flow condition are: Mae = 2.5, Re1b = 3000, T0 = 333.3K, and Tw/Tadb = 1. The three methods tested by Malik are a fourth-order compact
finite-difference (4CD) scheme, a single domain spectral collocation (SDSP) method and a multi-domain spectral collocation (MDSP) method. He applied a global method to compute all eigenvalues and purified them with a local iteration precedure. Hu and Zhong used a spectral collocation (SC) global method. Tab. 1 shows that our results agree very well this those of the other authors.
Table 1
Comparison of the eigenvalue solution of complex frequency u for the linear stability of a compressible boundary-layer over a flat plate with zero pressure gradient at Mae = 2.5, Re^ = 3000, To = 333.3K, Tw/Tadb = 1.0, a = 0.06, 0 = 0.1
Methods Grids u
4CD (Malik [9]) 61 (0.0367321, 0.005847)
SDSP (malik [9]) 61 (0.0367339, 0.005840)
MDSP (malik [9]) 61 (0.0367340, 0.005840)
SC (Hu, Zhong [10]) 100 (0.0367337, 0.005845)
4CD (Hu, Zhong [10]) 100 (0.0367338, 0.005840)
SC (Pres. Calculation) 90 (0.0367336, 0.005842)
The second test case concerns the compressible stability equations in the limit of Mach number approaching 0 to simulate the Orr — Sommerfeld results. Calculations were actually done for Mae = 10-6 and Re1b = 580. Some of the computed eigenvalues are compared in Tab. 2 with the results of Mack [11] using an Orr — Sommerfeld (OS) solver and Malik [9]'s MDSP methods for the compressible stability equations, where the real and imaginary parts of the computed phase velocity (c = u/a) are presented. We note, that our results agree very well with those computed by Mack [11] and Malik [9]. Also the second mode in Tab. 2, which is an eigenvalue of the energy equation and therefore cannot be calculated with the Orr — Sommerfeld solver, is well resolved.
Table 2
Comparison of the first 5 modes (phase velocity c = u/a) of the compressible stability equations in the incompressible limit (Mae = 10-6) with the Orr — Sommerfeld modes computed by Mack [11] and Malik at Re1b = 580, a = 0.179, 0 = 0
Pres. Calculation Malik [9] (MDSP) MACK [11] (OS)
(0.3641, 0.0079) (0.2329, —0.1344) (0.2897, —0.2769) (0.4840, —0.1921) (0.5573, —0.3654) (0.3641, 0.0079) (0.2329, —0.1343) (0.2897, —0.2768) (0.4839, —0.1921) (0.5571, —0.3655) (0.3641, 0.0080) (0.2897, —0.2769) (0.4839, —0.1921) (0.5572, —0.3653)
Furthermore, results of a stability investigation of the supersonic flow field over a blunt flat plate with a nose radius of Rn = 2.5 mm at a free-stream Mach number Maw = 2.5, a Reynolds number per unit length Re^/l* = 9-9-106/m and angle of attack a = 70 are compared with those obtained from the spatial theory using the NOLOT (NOnLOcal Transition analysis) stability solver [12] of the DLR1. We remark that the NOLOT code is based on a local eigenvalue iteration procedure. The lower computational effort of this method allows one to cluster more points in critial layers than for the global method. Asymptotic vanishing disturbance amplitudes are assumed at the upper boundary. As the spatial theory yields complex wavenumbers for a given
1Institut fur Stromungsmechanik, Deutsches Zentrum fiir Luft- und Raumfahrt e. V. Germany
real frequency and the temporal theory a complex frequency for given real-value wavenumbers, the frequencies, determined from both theories, can be compared. The reduced frequencies wred = 2nlb/(Relbue) x 105 are shown in Tab. 3 for the first boundary-layer instability (BL) mode and an instability mode of the inviscid entropy-layer (ES). The slight differences in the reduced frequencies can be explained by the limited resolution of the temporal scheme especially outside the boundary-layer. A subsequent local iteration combined with a multi-domain decomposition might improve the accuracy of the global results.
Table 3
Comparison of temporal/spatial solutions of the reduced frequency wred = 2nlb/(Reibue) x 105 for the supersonic flow field over a blunt body of Rn = 2.5 mm at Maœ = 2.5,
Re^/l * = 9.9 ■ 106/m, a = 70
Reib (Re(a), Re(fi)) wred Temp. Theory wred Spat. Theory [12] Type
547.7 (6.7529E-2, 7.4844E-2) 6.2009 6.2006 BL-Mode
593.8 (5.0882E-2, 0) 0.8843 0.8833 ES-Mode
2037.6 (3.3300E-2, 6.3100E-2) 0.6496 0.6483 BL-Mode
2112.5 (1.2092E-2, 0) 0.5440 0.5454 ES-Mode
3. Results
The linear stability of the supersonic flow field over a flat plate with a semicircle leading edge of Rn = 2.5 mm is investigated at Maw = 2.5, Re^/l* = 9.9 • 106/m, a = 70. The shock conditions Eq. (18) are imposed at the upper boundary to study the effect of the shock on the entropy-layer and first boundary-layer instability modes for different positions x1 = s*/Rn in streamwise direction. The wall-normal coordinate in the following figures in non-dimensionalized with characteristic boundary-layer length scale according to n = y*/l* \/(Rezb).
For the entropy-layer mode, we pass the region of high flow gradients around the body slope singularity and start at position x1 = 20.2. Here Fig. 1 shows clearly differences between the eigenfunctions computed using the asymptotic approch and those with imposed shock conditions. At the same streamwise position, Fig. 2 demonstrates, that the amplitude of the pressure eigenfunction has a finite amplitude at the shock wave for the shock conditions in contrast to the asymptotic condition. Further downstream, at position x1 = 101.2 the eigenfunctions for both boundary conditions fall complety into each other, which is documented in Fig. 3.
Furthermore, the shock has no effect on the boundary-layer instability mode at position x1 = 2.95, as shown in Fig. 4. The eigenfunctions decay exponentially and therefore do not notice the presence of the shock. Going downstream, the shock distance increases faster than the critical layer shifts away from the wall, so we expect to find no shock influence on the boundary-layer stability at higher streamwise positions.
4. Concluding remarks
The linearized shock conditions Eq. (18), derived from the unsteady Rankine — Hugoniot relations yield to a more realistic mathematical modelling for stability analysis of supersonic flows over blunt bodies including shock waves compared to the common asymptotic approach. However, far downstream the nose region, the simplified conditions imposed at the upper boundary, are a rational approximation due to the fact that disturbance waves from the flow field windward of the shock vanish before they can interact with the shock wave.
400
!
350 300 250 200 150 100 50
...... 1 i i I u asympt. be.-
\ : fasympt.be.----
v P asympt. be. —
u shock be.........
\ f shock be.
u 1 P shock be. -
gv
!\
\\
-
K
\ A. -----------
i
---------------
| - -35 1 1 1
0 0.2 0.4 0.6 0.8 1
[?/?max][-]
Fig. 1. Wall-normal distribution of the eigenfunction U, T and P of the entropy-layer instability mode at streamwise position x1 = 20.2.
Fig. 2. Scope of the pressure eigenfunction of the entropy-layer instability mode p near the shock wave position.
Fig. 3. Wall-normal distribution of the
eigenfunction U, T and P of the entropy-
layer instability mode at streamwise position x1
-1 - 101.2.
??50 45 40 35 30 25 20 15 10 5
0
1 1 1 1 1 1 1 1 1 »asympt, be.-
T asympt. be,---- Pasympt.be.----- u shock be.........
A T shock be.----" P shock be.----_
i\ -
_
( -1-1-T- l—~rH
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Fig. 4. Wall-normal distribution of the eigenfunction U, T and P of the first boundary-layer instability at streamwise position x1 = 2.95.
References
[1] Lysenko V. I. Influence of the entropy layer on the stability of supersonic shock layer and transition of the laminar boundary layer to turbulence //J. Appl. Mech. Tech. Phys. 1990. Vol. 31, No. 6. P. 868-873.
[2] Chang C. L., Malik M. R. Effects of Shock on the Stability of Hypersonic Boundary Layer. AIAA. Paper 90-1448. 1990.
[3] KovAsznay L.S.G. Turbulence in supersonic flow // J. Aeronaut. Sci. 1953. Vol. 20. P. 657-674.
[4] PETROV G. V. Stability of a thin viscous layer on a wedge in hypersonic flow of a perfect gas: in laminar-turbulent Transition // Proc. of 2nd IUTAM Symp. / Ed. by V. V. Kozlov. 1984.
[5] Wada Y., Liou M.-S. A Flux Splitting Scheme With High-Resolution and Robustness for Discontinuities. NASA TM 106452 (AIAA 94-0083. 1994).
[6] Van Albada G., Van Leer B., Roberts J. A comparative study of computational methods in cosmic gas dynamics // Astron. Astrophys. 1982. Vol. 108. P. 76-84.
[7] Van Leer B. Towards the ultimate conservative difference scheme V. A second-order sequel to Godunov's method // J. of Comput. Phys. 1979. Vol. 32. P. 101-136.
[8] DiETZ G., MEIJERING A. Numerical investigation of the boundary layer instabilities over a blunt flat plate at angle of attack in supersonic flow // Notes on Numer. Methods. 1996. P. 103-110.
[9] Malik M. R. Numerical methods for hypersonic boundary layer stability // J. of Comput. Phys. 1988. Vol. 86, No. 2. P. 376-413.
[10] Hu S. H., ZHONG X. Linear Stability of Hypersonic Couette Flow. AIAA. Paper 97-0432. 1997.
[11] Mack L. M. Boundary-Layer Stability Theory. AGARD. Rep. 709. 1984.
[12] Hein S., BERTOLOTTI F.P., SiMEN M., Hanifi A., Henningson D. Linear Nonlocal Instability Analysis. The Linear NOLOT Code. DLR IB 223-94 A56. 1994.
Received for publication June 22, 2001
