Transient growth analysis of a compressible boundary layer
Transcript Of Transient growth analysis of a compressible boundary layer
DOI: 10.13009/EUCASS2019-591
8TH EUROPEAN CONFERENCE FOR AERONAUTICS AND AEROSPACE SCIENCES (EUCASS)
DOI: ADD DOINUMBER HERE
Transient growth analysis of a compressible boundary layer using two-dimensional linear stability theory
Iván Padilla Montero† and Fabio Pinna Aeronautics and Aerospace Department, von Karman Institute for Fluid Dynamics
Chaussée de Waterloo 72, 1640 Rhode-Saint-Genèse, Belgium [email protected] · [email protected] †Corresponding author
Abstract
Transient growth results are presented for a compressible boundary layer at Mach 2.5 by performing a singular value decomposition of the two-dimensional linear stability operator. It is shown that the nonmodal growth predicted by two-dimensional linear stability theory corresponds to the envelope of all the transient growth curves associated to each of the one-dimensional linear stability (LST) computations at the spanwise wavenumbers that can be resolved by the two-dimensional discretization. Emphasis is made on the choice of the domain size in the spanwise direction and how this can affect the resulting energy gain.
1. Introduction
Boundary-layer transition from a laminar to a turbulent regime is a critical driver in the design of high-speed vehicles. Turbulent supersonic and hypersonic boundary layers are usually characterized by big heat transfer rates that lead to strong aerothermodynamic loads. For example, for a Reynolds number of 107 and a moderate supersonic Mach number, the turbulent heat-transfer coefficient on a flat plate is estimated to be an order of magnitude higher than the corresponding laminar value.27 Despite significant research efforts in the last decades, the physical mechanisms governing high-speed boundary layer transition are not well understood yet, and the existing transition prediction tools for practical applications still rely heavily on empirical correlations and extensive wind tunnel testing.18
In order to improve the transition prediction capabilities for future high-speed applications, new methodologies involving a higher degree of flow physics for each particular case should be developed. To achieve this goal, theoretical research in laminar-turbulent transition has traditionally focused on the study of small-amplitude disturbances undergoing exponential growth based on the mechanisms predicted by linear modal stability theory. Examples include the development of Tollmien-Schlichting waves on a flat plate or crossflow instabilities in a three-dimensional boundary layer. However, this scenario of modal growth does not always provide a satisfactory explanation for the occurrence of transition encountered in experimental and in-flight observations. Probably, the most notable example of this situation is the so-called blunt-body paradox, which makes reference to the occurrence of transition observed in blunt bodies at conditions which are stable to modal perturbations. In view of the need of an alternative explanation, nonmodal growth, also known as transient growth, has emerged as a possible mechanism leading to a different path to transition, namely, by-pass transition.19,20
The first application of the theoretical transient growth framework to compressible boundary layers is due to Hanifi et al.,5 who employed a singular value decomposition of the linearized evolution operator based on classical temporal linear stability theory (LST). They showed that a large amount of transient growth can take place over a wide range of parameter values for compressible flow, and that the optimal initial perturbations for the compressible case are found to be streamwise vortices generated by the lift-up effect. Tumin & Reshotko26 later extended the same analysis to spatially growing disturbances. Building upon this work, Reshotko & Tumin21 developed a correlation for roughness-induced transition based on spatial transient growth results which is able to reproduce transition trends observed experimentally, giving weight to the idea that transient growth could be a plausible mechanism for early transition due to distributed surface roughness. A different theoretical approach was independently developed by Andersson et al.2 and Luchini9 to quantify the spatial transient energy growth in boundary layers including non-parallel effects. This technique is based on the parabolized stability equations (PSE) and consists in an optimization marching procedure in which the disturbance energy over a given spatial extent is maximized. Recently, by employing this methodology, Paredes and
Copyright © 2019 by all authors. Published by the EUCASS association with permission.
DOI: 10.13009/EUCASS2019-591
TRANSIENT GROWTH IN A COMPRESSIBLE BOUNDARY LAYER USING 2D-LST
co-authors15 have extended the current body of results for optimal nonmodal disturbance growth in boundary layer flows to the hypersonic Mach number regime. Additionally, Paredes et al.14 have revisited the blunt-body paradox using an improved transient growth framework, leading to some modifications of the Reshotko-Tumin correlation but at the same time highlighting the need to further investigate the optimal-growth criterion underlying this correlation. Similarly, Hein et al.6 have investigated the disturbance growth in the wake of a roughness patch mounted on the forebody of a reentry capsule at Mach 5.9, for which the onset of transition is observed experimentally. Nevertheless, no evidence of modal or nonmodal growth could be found in that study. Despite these efforts, transient growth is not yet conclusively linked to the measured onset of transition over blunt-body configurations.
In the last two decades, aided by the rapid increase in computational power, classical linear stability theory has been extended to base flows featuring two and even three inhomogeneous spatial directions, giving birth to twoand three-dimensional linear stability theories, also known as BiGlobal and TriGlobal theory,24 respectively. These theories allow the study of more complex flows at an increased computational cost. Almost all of the transient growth studies available nowadays that deal with compressible basic flows are based either on classical linear stability theory or parabolized stability equations, in both cases assuming a strong inhomogenous character only in one spatial direction. Very recently, Quintanilha et al.7 have performed temporal transient growth computations on the HIFiRE-5 elliptic cone model in hypersonic flow, using two-dimensional linear stability theory (2D-LST). The results reveal a significant amount of transient growth taking place at short times, after which modal unstable perturbations take the lead and the growth becomes exponential. In the context of high-speed flows, a very interesting application of 2D-LST concerns the wake induced behind isolated roughness elements in supersonic and hypersonic flow. Such three-dimensional elements with heights comparable to the local boundary-layer thickness generally induce counter-rotating streamwise vortices in the wake flow, giving rise to low-velocity streaks that are surrounded by regions of high-detached shear. Several fundamental studies exist nowadays concerning the linear modal instability of these configurations, see for instance De Tullio et al.,4 Theiss et al.,23 Padilla Montero & Pinna.13 These analyses have revealed that the roughness wake supports the growth of sinuous and varicose instability modes that develop in the high-shear regions introduced by the counter-rotating vortex pair, and that these modes can undergo a substantial growth during the linear stages of the transition process. However, to the best of the authors’ knowledge, no studies currently exist which analyze the potential for nonmodal growth in the wake behind isolated roughness elements in high-speed flow.
The aim of this work is to highlight the particularities on the transient growth analysis of compressible boundary layers by means of two-dimensional linear stability theory and validate its implementation within the von Karman Institute stability code known as VESTA toolkit.17 This, in turn, sets the base for future transient growth computations on the wake induced by a discrete roughness element inside a hypersonic boundary layer. On first place, the mathematical formulation of the linear stability problem and the associated transient growth theoretical background are summarized. Next, some details of the current numerical implementation are discussed and a validation of the transient growth solver is presented. Finally, results are shown for 2D-LST transient growth computations in a self-similar boundary layer at Mach 2.5, emphasizing the comparison with the nonmodal growth predicted by classical linear stability theory.
2. Mathematical formulation
2.1 Governing equations
The equations governing the stability solutions obtained in this work are the Navier-Stokes equations for a compress-
ible Newtonian fluid. The primitive flow variables are denoted by density ρ, pressure p, temperature T and velocity
components ui, with u1 = u, u2 = v and u3 = w. In a Cartesian reference frame, the system of equations can be written
in conservation form as
∂ρ + ∂(ρui) = 0,
(1)
∂t ∂xi
∂(ρui) + ∂(ρuiu j) + ∂p − ∂τi j = 0, (2)
∂t
∂xj
∂xi ∂x j
∂(ρE) + ∂(ρEui + pui) + ∂qi − ∂(uiτi j) = 0, (3)
∂t
∂xi
∂xi ∂x j
where t is the time coordinate, xi is the ith spatial coordinate (x1 = x, x2 = y and x3 = z) and E = e + uiui/2 is the total energy of the fluid, with e being the specific internal energy. The viscous stress tensor τi j is defined as
∂ui ∂u j 2 ∂uk τi j = µ ∂x j + ∂xi − 3 ∂xk δi j , (4)
2
DOI: 10.13009/EUCASS2019-591
TRANSIENT GROWTH IN A COMPRESSIBLE BOUNDARY LAYER USING 2D-LST
where µ is the dynamic viscosity of the fluid and δi j is the Kronecker delta. The conductive heat flux qi is modeled
by means of Fourier’s law of heat conduction, qi = −k (∂T/∂xi), for a fluid with thermal conductivity denoted by k. A
calorically perfect gas is assumed, so the system of equations is closed using the perfect gas equation of state and two
thermodynamic relationships:
R
p = ρRT,
e = cvT and
cv
=
γ
−
, 1
(5)
where R is the specific gas constant, cv is the specific heat at constant volume and γ is the ratio of specific heats. Air
is considered with the values γ = 1.4 and R = 287 J/(kg·K). Sutherland’s law is used to model the variation of the
dynamic viscosity with temperature
T 3/2 Tre f + S µ
µ = µre f Tre f
T +Sµ ,
(6)
where S µ = 110.4 K is the Sutherland temperature constant and µre f = 1.716 × 10−5 kg/(m·s) and Tre f = 273.15 K are respectively the reference dynamic viscosity and static temperature. Regarding the thermal conductivity, its value is determined from the assumption of a constant Prandtl number Pr, that is, k = γcvµ/Pr. A value of Pr = 0.7 is employed for all the cases shown in this study.
2.2 Linear stability theory
According to classical linear stability theory, the vector of primitive flow variables q = [u, v, w, ρ, T ]T is split into a steady reference state q¯ , also known as base flow, and a small unsteady perturbation field q˜ :
q = q¯ + q˜ ,
(7)
with 1. The base flow is assumed to be locally parallel in the streamwise (x) direction, so that q¯ = q¯ (y, z) at a given
x coordinate. The perturbations are assumed to have an inhomogeneous character in the spanwise and wall-normal
directions, so that the amplitude of the perturbations is considered to be a function of both y and z coordinates. These
assumptions lead to a two-dimensional linear stability theory (2D-LST), for which the ansatz of the modal perturbations
can be written as
q˜ (x, y, z, t) = qˆ (y, z) exp[i(αx − ωt)] + c.c.,
(8)
where qˆ is the two-dimensional amplitude function, α is the wavenumber along the streamwise direction, ω is the
angular frequency and c.c. denotes the complex conjugate. For temporal transient growth analysis, the temporal sta-
bility framework is considered. The quantity α is real and represents the streamwise wavenumber of the perturbations.
On the contrary, ω is complex, with the real part ωr = {ω} being the angular frequency of q˜ and the imaginary
part ωi = {ω} its temporal growth rate. With this definition, a negative value of ωi means a temporal decay of the
amplitude function whereas a positive value implies temporal growth.
The equations governing the linear stability problem are obtained by substituting the splitting given by equa-
tion (7) together with the ansatz defined by equation (8) into the Navier-Stokes system, equations (1) to (3), and
neglecting the non-linear terms, which are of order O( 2). The resulting system of linear partial differential equations
can be written in the following compact form
Aqˆ = ωBqˆ ,
(9)
where A and B are complex and nonsymmetric differential matrix operators. After discretization of equation (9), a twodimensional algebraic generalized eigenvalue problem (GEVP) is obtained, which can be solved to obtain the angular frequency (eigenvalue) and the amplitude function (eigenvector) of each eigenmode supported by the linear stability problem.
The stability equations are made nondimensional using the boundary layer edge properties, denoted by the subscript e. The following dimensionless quantities are defined
t† = tue , x† = xi , ρ† = ρ , u† = ui , p† = p , E† = E , T † = T , µ† = µ and k† = k , (10)
le i le
ρe i ue
ρeu2e
u2e
Te
µe
ke
where le is the local Blasius length scale at the streamwise location x where the stability analysis is performed, given by le = µe x/ρeue. This choice of reference quantities leads to a set of equations where the similarity parameters are γ and the boundary layer edge Reynolds, Mach and Prandtl numbers, respectively denoted by Re, Me and Pr. In the remaining sections of this manuscript, the nondimensional equations are considered. The superscript † is dropped
unless otherwise noted.
3
DOI: 10.13009/EUCASS2019-591
TRANSIENT GROWTH IN A COMPRESSIBLE BOUNDARY LAYER USING 2D-LST
2.3 Transient growth analysis
In order to quantify the potential for transient energy growth in a given flow field satisfying the assumptions discussed
in the previous section, the orthogonality of the eigenfunctions and the size of the perturbations obtained by linear
stability theory have to be evaluated. The first can be achieved by the definition of an inner product and the second by an associated energy norm. In this work, the disturbance energy norm introduced by Chu,3 Mack10 and later rederived by Hanifi et al.5 is adopted, given by
E(t) = 1 ρ¯ u˜2 + v˜2 + w˜ 2 + T¯ ρ˜2 +
ρ¯
T˜ 2 dxdydz
(11)
2V
γρ¯ Me2 γ (γ − 1) T¯ Me2
in a volume V. For transient growth analysis, it is convenient to express the disturbance energy in terms of the timedependent amplitude functions, defined as qˇ (y, z, t) = qˆ (y, z) exp(−iωt), so that the ansatz (equation (8)) becomes
q˜ (x, y, z, t) = qˇ (y, z, t) exp(iαx) + c.c.
(12)
Additionally, since the perturbations are assumed to be periodic along the streamwise direction, the disturbance energy can be integrated over one wavelength in x without loss of generality, see Weder et al.28 Substituting equation (12) into equation (11) and integrating along x between 0 and 2π/α, the following expression is obtained
E(t) = 2π ρ¯ (uˇ∗uˇ + vˇ∗vˇ + wˇ ∗wˇ ) + T¯ ρˇ∗ρˇ +
ρ¯
Tˇ ∗Tˇ dydz,
(13)
αΩ
γρ¯ Me2
γ (γ − 1) T¯ Me2
where ∗ denotes the complex conjugate and Ω is the space in which the two-dimensional eigenfunctions are defined.
equation (13) can be rewritten in a compact form as
2απ E(t) = ||qˇ ||2E = (qˇ , qˇ )E , (14)
where the subscript E denotes the energy norm. The associated inner product for a given pair of amplitude functions is then expressed as
(qˇ k, qˇ l)E = qˇ kHMqˇ l dydz,
(15)
Ω
where the superscript H refers to the complex conjugate transpose and the matrix M is given by
T¯
ρ¯
M = diag ρ¯, ρ¯, ρ¯, γρ¯ M2 , γ (γ − 1) T¯ M2 .
(16)
e
e
After a choice on the inner product and the associated energy norm has been made, the temporal transient growth analysis requires an eigenvector expansion of the time dependent amplitude function. Denoting the first K eigenvalues and eigenvectors obtained from the solution of the eigenvalue problem given by equation (9) as ωk and qˆ k, respectively, the vector qˇ can be expanded as
K
qˇ (y, z, t) = κk(t)qˆ k(y, z)
(17)
k=1
with
κk(t) = κk(0) exp(−iωkt),
(18)
or in matrix form
qˇ (y, z, t) = Qκ(t) with κ(t) = Λκ(0) and Λ = diag exp(−iω1t), ..., exp(−iωKt) ,
(19)
where Q is a matrix whose columns contain the eigenvectors qˆ k and the vector κ(0) = [κ1(0), ..., κK(0)]T contains the initial expansion coefficients that determine the individual contribution of each mode to the transient response. By employing the inner product defined by equation (15) together with the expansion in equation (17), the following
relationship can be found between the energy norm and the L2-norm
(qˇ k, qˇ l)E = (Qκk, Qκl)E = κk∗QHMQκl dydz = κk∗Dκl = κk∗FHFκl = (Fκk, Fκl)2 ,
(20)
Ω
4
DOI: 10.13009/EUCASS2019-591
TRANSIENT GROWTH IN A COMPRESSIBLE BOUNDARY LAYER USING 2D-LST
where the subscript 2 denotes the L2-norm. The matrix F is the Cholesky decomposition of the Hermitian and positive definite matrix D, that is, D = FHF, which is defined element-wise as
Dkl = qˆ kHMqˆ l dydz.
(21)
Ω
The relation given by equation (20) allows to rewrite the disturbance energy norm in terms of the L2-norm as follows
||qˇ (t)||2E = ||Fκ(t)||22.
(22)
Finally, the maximum transient energy growth G(t), optimized for all the possible initial conditions, is given by (see Schmid & Henningson22)
G(t) = max ||qˇ (t)||2E = max ||Fκ(t)||22 = max ||FΛF−1Fκ(0)||22 = ||FΛF−1||2.
(23)
||qˇ (0)||2E ||Fκ(0)||22 ||Fκ(0)||22 2
The L2-norm of matrix S = FΛF−1 can be computed by means of singular value decomposition (SVD). The value of G(t) is then given by the square of the largest singular value of matrix S. The optimal initial conditions giving
rise to the maximum transient growth can be computed via the right singular eigenvector r associated to the largest singular value of matrix S, that is, κ(0) = F−1r.
3. Numerical methodology
The numerical solution of the linear stability problem and the posterior transient growth analysis are performed using the von Karman Institute Extensible Stability & Transition Analysis (VESTA) toolkit, originally developed by Pinna.16,17 The discretization of the partial differential eigenvalue problem is performed employing the Chebyshev collocation method.25 This technique is based on a Lagrange polynomial interpolation constructed in a structured grid with a non-uniform point distribution given by the Chebyshev-Gauss-Lobatto collocation points. These points are defined on a transformed domain with spanwise and wall-normal coordinates respectively denoted by ξ and η, with ξ, η ∈ [−1, 1]. This transformed space does not coincide with the physical domain of interest. Therefore, geometrical mappings have to be used. In the computations presented in this work, the mapping originally introduced by Malik11 is applied in the wall-normal direction, which allows placing half of the collocation points below a given coordinate in each direction, whereas a uniform mapping is used in the spanwise direction, which simply translates the transformed domain to the physical one. In this way, the grid used for solving the stability problem can be clustered towards the boundary layer. On the other hand, despite this process transforms the stability grid into the physical domain, the mapped grid defined by the Chebyshev-Gauss-Lobatto collocation points is different than the mesh employed to obtain the base flow solution. As a consequence, before solving the eigenvalue problem, the base flow data is interpolated on the collocation grid by means of a cubic spline interpolation in each spatial direction.
The boundary conditions imposed on the perturbation quantities are as follows. At the wall, the velocity and the temperature perturbations are set to zero via a homogeneous Dirichlet condition, whereas the density fluctuation is determined by means of a compatibility condition. In the wall-normal far-field boundary, the perturbations are forced to decay by imposing a Dirichlet boundary condition as well, with the exception of ρˆ, which once again satisfies a compatibility condition. Regarding the spanwise boundaries, periodic boundary conditions are specified. This implies that the value of each disturbance quantity and its first derivative normal to the boundary are set to be equal at both spanwise boundaries.
For the transient growth analysis to be physically meaningful, a large number of eigenmodes must be included in the construction of the matrix S, see for instance Hanifi et al.5 In this work, the QZ algorithm12 (also known as generalized Schur decomposition) is used to compute the full spectrum of eigenvalues and eigenvectors associated to the discrete generalized eigenvalue problem, employing subroutines belonging to the LAPACK1 library. The integral in equation (21) is evaluated using the Chebyshev integration weights derived in Hanifi et al.5 This ensures that the order of accuracy of the discretization is preserved after integration. The calculation of the largest singular value of matrix S is also performed using the LAPACK library.
4. Validation
In order to validate the implementation of the transient growth analysis methodology within the von Karman Institute (VKI) stability code (VESTA), two test cases investigated by Hanifi et al.5 are reproduced here based on classical
5
DOI: 10.13009/EUCASS2019-591
TRANSIENT GROWTH IN A COMPRESSIBLE BOUNDARY LAYER USING 2D-LST
Figure 1: Temporal evolution of transient growth for a self-similar boundary layer at the conditions described in section 4. Comparison of the results obtained between VESTA and the reference5 data for one-dimensional linear stability theory.
linear stability (LST). Both cases employ a base flow obtained by solving the compressible self-similar boundary layer equations, see for instance Anderson.8 These equations are solved as well within VESTA by means of a fourth-order Runge-Kutta method. The self-similar profiles considered for validation have a boundary layer edge Mach number of Me = 2.5, a total temperature of T0 = 333 K and assume an adiabatic wall. The first case evaluates the self-similar profile at a Reynolds number of Re = 300, with a streamwise wavenumber α = 0 and a spanwise wavenumber β = 0.1. The second configuration assumes a Reynolds number of Re = 3000 and wavenumbers α = 0.06 and β = 0.1.
Figure 1 shows the energy growth as a function of time for each of the two cases as computed by the current solver. A very good agreement with the reference results is obtained. At the conditions of the first test case, the boundary layer is stable to linear modal perturbations. The energy reaches a maximum transient growth at a finite time and then progressively decays to zero. On the contrary, the second configuration features a linearly unstable mode, and as a result the energy grows exponentially as time tends to infinity. Nevertheless, for short times the disturbance energy still undergoes a significant transient growth, showing a faster increase than the exponential growth associated to the unstable mode.
5. Transient growth results based on two-dimensional linear stability theory
This section presents results obtained for transient growth computations using eigenmodes solved by means of twodimensional linear stability theory and their comparison against the one-dimensional counterpart. The base flows considered in this analysis have conditions that are similar to the validation test cases shown in the previous section, namely, a self-similar boundary layer at Me = 2.5, with a total temperature T0 = 333 K and an adiabatic wall. For comparison purposes between LST and 2D-LST results, the base flow wall-normal velocity is set to zero in all the 2D-LST computations.
5.1 Case 1: Re = 3000, α = 0.06
The first case considered corresponds to the same Reynolds number and streamwise wavenumber as for the second validation test case, that is, Re = 3000 and α = 0.06. Before proceeding further, it is important to note that for the 2DLST computations a choice has to be made regarding the spanwise size of the stability domain, denoted by Lz. Together with the use of periodic boundary conditions, this size determines the specific wave-numbers that are resolved by the computation. For a given value of Lz, only modes with spanwise wanumbers that are smaller multiples of, or equal to, 2π/Lz can be computed. Additionally, the number of grid points in the spanwise direction, denoted by Nz, limits the number of spanwise multiple modes that can be retrieved for each LST mode. For a given discretization, Nz − 2 multiple modes are computed. For instance, if a value of Lz = 2π/0.1 is chosen, this means that, in addition to all the modes obtained by the classical LST operator for β = 0.1, multiple modes with spanwise wavenumbers β = n0.1, with n = 2, ..., Nz − 2, are also found in the 2D-LST spectrum.
Figure 2 displays the transient growth results obtained for the current case by employing two different spanwise domain lengths, namely, Lz = 2π/0.1 and Lz = 2π/0.05 and a grid resolution of Nz × Ny = 61 × 61 points. LST results
6
DOI: 10.13009/EUCASS2019-591
TRANSIENT GROWTH IN A COMPRESSIBLE BOUNDARY LAYER USING 2D-LST
Figure 2: Comparison of transient growth results between LST and 2D-LST for two different spanwise domain lengths at the conditions of case 1.
for different spanwise wavenumbers are also shown for comparison, which have been computed using 151 collocation points in the wall-normal direction. As it can be observed, the transient energy growth predicted by means of 2D-LST is equal to the envelope of the transient growth curves associated to each of the LST computations at the spanwise wavenumbers that can be resolved by the two-dimensional discretization. In other words, when the domain length of Lz = 2π/0.1 is used, the 2D-LST transient growth is the envelope of the LST transient growth curves for β = n0.1, with n = 1, ..., 59 in this case. The four multiple spanwise wavenumbers for which the LST operator undergoes the highest transient growth are also included in the left plot of Figure 2 to illustrate this fact. As it can be seen, a linearly unstable wave is only present at β = 0.1 for this case, so the 2D-LST energy growth converges to the exponential growth given by the growth rate of this particular mode in the limit of large times. However, for smaller times, the other wavenumbers exhibit a larger nonmodal growth, and as a result the 2D-LST curve also undergoes a larger growth accordingly. It is important to note that when a domain length twice as big is considered, the spanwise wavenumbers that can be resolved are different. Now, all modes with β = n0.05 are included, and this results in a different envelope mainly due to the transient growth behaviour of the LST operator at β = 0.15, for which the least stable mode has a growth rate that is very close to 0 and the energy decays very slowly.
To better illustrate the relationship between the one-dimensional and two-dimensional linear stability analyses, Figure 3 presents a comparison of the LST and the 2D-LST spectra obtained for the current case. Figure 3a) shows an overview of the relevant part of the temporal spectrum as computed with LST for β = 0.1 and Ny = 81 points and with 2D-LST using Nz = 31, Ny = 81 points and Lz = 2π/0.1. The spectrum contains a vertical continuous branch at ωr = α = 0.06 which consists of entropy and vorticity waves and two horizontal continuous branches which contain acoustic waves. An additional group of modes is resolved in the positive imaginary part on top of the vertical continuous branch. These are numerical spurious modes that appear as part of the numerical discretization. The discrete physical eigenmodes are found in the diagonal-like structures that emerge from the vertical continuous branch. A single discrete unstable is found at these conditions, which corresponds to a first-mode disturbance. All the modes present in this plot are included in the transient growth calculation for both the LST and the 2D-LST theories, with the exception of the spurious numerical modes. Figure 3b) displays a closer view of the same spectra plotted in a) where the difference between the one and two-dimensional computations can be clearly appreciated. As it can be noticed, for each mode appearing in the LST spectrum, there are Nz − 2 multiple modes resolved in the 2D-LST spectrum. Figure 3c) is included as a proof of convergence when the number of discretization points in the spanwise direction is doubled. Finally, Figure 3d) illustrates the comparison between the 2D-LST spectra obtained for both spanwise domain sizes considered. When a domain twice as big is employed, additional discrete modes can be found in the spectrum, respectively corresponding to the spanwise wavenumbers β = 0.05, 0.15, 0.25 etc., which are not resolved by the other domain length.
When performing 2D-LST transient growth analyses of more complex boundary layer flows, care should then be taken to choose an appropriate domain size. Ideally, a previous analysis based on LST should be performed to identify the spanwise wavenumbers that produce the higher contributions to the transient growth envelope.
7
DOI: 10.13009/EUCASS2019-591
TRANSIENT GROWTH IN A COMPRESSIBLE BOUNDARY LAYER USING 2D-LST
Figure 3: Comparison of the temporal stability spectrum between LST and 2D-LST for the conditions analyzed in case 1.
5.2 Case 2: Re = 300, α = 0 A second case is also examined which corresponds to conditions at which the optimal transient growth is found for a Reynolds number of Re = 300. From LST-based computations, see Hanifi et al.,5 the optimal growth is found for α = 0 and β = 0.235. A comparison of the transient growth results obtained for this case is depicted in Figure 4. In this case, the LST transient growth evolution for β = 0.235 dominates for all times. As a result, the 2D-LST curve when Lz = 2π/0.235 is identical to the LST one. The advantage of this case is that it can be used to compare the optimal disturbances computed with both theories as well. Figure 5 shows the optimal perturbations evaluated at t = tmax, where tmax denotes the time at which maximum G is achieved, which for this case corresponds to t = 980. Reference values for the streamwise velocity and temperature perturbations from Hanifi et al.5 are also added for validation, for which a very good agreement is found with the current implementation. The initial optimal disturbance for this case corresponds to streamwise vortices that develop into streamwise streaks as time evolves. The amplitude functions of the optimal perturbations shown here correspond to the signature of these streak structures. It is interesting to note that the relative magnitude of the wall-normal and spanwise velocity perturbations is much smaller compared to the ones associated to the streamwise velocity and the temperature. Similarly, Figure 6 presents the optimal perturbations computed by means of the transient growth computation based on 2D-LST modes. As it can be observed, the shape of the optimal disturbances as well as their relative magnitude coincides with the results obtained with LST. Additionally, they feature a spanwise wavelength equal to the size of the domain for the spanwise direction, which means that β = 0.235 also for the optimal perturbations computed with the two-dimensional eigenmodes.
The results shown in this section constitute a validation of the current transient growth implementation based on two-dimensional linear stability theory. 8
DOI: 10.13009/EUCASS2019-591
TRANSIENT GROWTH IN A COMPRESSIBLE BOUNDARY LAYER USING 2D-LST
Figure 4: Comparison of transient growth results between LST and 2D-LST for the conditions of case 2.
Figure 5: Optimal disturbances at t = tmax computed with LST for the conditions studied in case 2.
6. Conclusions
Transient growth computations based on two-dimensional linear stability theory have been performed for a self-similar boundary layer at Mach 2.5 by means of singular value decomposition of the linearized temporal stability operator. Two different Reynolds numbers and spanwise wavenumbers have been considered, resulting first in a case where the boundary layer is linearly unstable to modal disturbances but undergoes a significant amount of transient growth for small times, and second in a linearly stable base flow at the conditions for which the optimal transient growth is encountered. It has been found that the transient energy growth predicted by two-dimensional linear stability theory is equal to the envelope of the nonmodal classical linear stability results evaluated at all the possible spanwise wavenumbers that can be resolved by a given spanwise discretization of the two-dimensional domain. Care should be taken to size the two-dimensional stability domain appropriately in the spanwise direction in order to resolve all the modes with a wave-like behavior in the spanwise direction that have the largest contribution to the energy growth.
The optimal disturbances have been computed for the second studied case at the time instant in which the maximum nonmodal growth is present. A comparison of the optimal perturbations computed between LST and 2D-LST theories shows that the two-dimensional disturbances have the same shape and magnitude as the one-dimensional ones and feature a spanwise wavenumber equal to the optimal one as well. These results constitute a validation of the current transient growth implementation for two-dimensional eigenmodes, and sets the base for the application of this methodology to more complex base flows with two-inhomogenous spatial directions, such as the wake induced behind a discrete roughness element inside a high-speed boundary layer.
9
DOI: 10.13009/EUCASS2019-591
TRANSIENT GROWTH IN A COMPRESSIBLE BOUNDARY LAYER USING 2D-LST
Figure 6: Optimal disturbances at t = tmax computed with 2D-LST for the conditions studied in case 2.
Acknowledgments
This work is part of a project that has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 675008.
References
[1] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen. LAPACK Users’ Guide. Society for Industrial and Applied Mathematics, Philadelphia, PA, third edition, 1999.
[2] P. Andersson, M. Berggren, and D. S. Henningson. Optimal disturbances and bypass transition in boundary layers. Phys. Fluids, 11(1):134–150, 1999.
[3] B. Chu. On the energy transfer to small disturbances in fluid flow (Part 1). Acta Mech., 1(3):215–234, 1965. [4] N. De Tullio, P. Paredes, N. D. Sandham, and V. Theofilis. Laminar-turbulent transition induced by a discrete
roughness element in a supersonic boundary layer. J. Fluid Mech., 735:613–646, 2013. [5] A. Hanifi, P. J. Schmid, and D. S. Henningson. Transient growth in compressible boundary layer flow. Phys.
Fluids, 8(3):826–837, 1996. [6] S. Hein, A. Theiss, A. Di Giovanni, C. Stemmer, T. Schilden, W. Schröder, P. Paredes, M. M. Choudhari, F. Li,
and E. Reshotko. Numerical investigation of roughness effects on transition on spherical capsules. J. Spacecraft Rockets, 56(2):388–404, 2019. [7] H. R. Quintanilha Jr., V. Theofilis, and A. Hanifi. Global transient-growth analysis of hypersonic flow on the HIFiRE-5 elliptic cone model. In 2019 AIAA Aerospace Sciences Meeting, AIAA SciTech Forum, San Diego, CA, USA, 7-11 January 2019, AIAA 2019-2148. [8] J. D. Anderson Jr. Hypersonic and high-temperature gas dynamics. American Institute of Aeronautics and Astronautics, second edition, 2006. [9] P. Luchini. Reynolds-number-independent instability of the boundary layer over a flat surface: optimal perturbations. J. Fluid Mech., 404:289–309, 2000. [10] L. M. Mack. Boundary layer stability theory. Technical report, Jet Propulsion Laboratory Report 900-277, 1969. [11] M. R. Malik. Numerical methods for hypersonic boundary layer stability. J. Comput. Phys., 86(2):376–413, 1990. [12] C. B. Moler and G. W. Stewart. An algorithm for generalized matrix eigenvalue problems. SIAM J. Numer. Anal., 10(2):241–256, 1973.
10
8TH EUROPEAN CONFERENCE FOR AERONAUTICS AND AEROSPACE SCIENCES (EUCASS)
DOI: ADD DOINUMBER HERE
Transient growth analysis of a compressible boundary layer using two-dimensional linear stability theory
Iván Padilla Montero† and Fabio Pinna Aeronautics and Aerospace Department, von Karman Institute for Fluid Dynamics
Chaussée de Waterloo 72, 1640 Rhode-Saint-Genèse, Belgium [email protected] · [email protected] †Corresponding author
Abstract
Transient growth results are presented for a compressible boundary layer at Mach 2.5 by performing a singular value decomposition of the two-dimensional linear stability operator. It is shown that the nonmodal growth predicted by two-dimensional linear stability theory corresponds to the envelope of all the transient growth curves associated to each of the one-dimensional linear stability (LST) computations at the spanwise wavenumbers that can be resolved by the two-dimensional discretization. Emphasis is made on the choice of the domain size in the spanwise direction and how this can affect the resulting energy gain.
1. Introduction
Boundary-layer transition from a laminar to a turbulent regime is a critical driver in the design of high-speed vehicles. Turbulent supersonic and hypersonic boundary layers are usually characterized by big heat transfer rates that lead to strong aerothermodynamic loads. For example, for a Reynolds number of 107 and a moderate supersonic Mach number, the turbulent heat-transfer coefficient on a flat plate is estimated to be an order of magnitude higher than the corresponding laminar value.27 Despite significant research efforts in the last decades, the physical mechanisms governing high-speed boundary layer transition are not well understood yet, and the existing transition prediction tools for practical applications still rely heavily on empirical correlations and extensive wind tunnel testing.18
In order to improve the transition prediction capabilities for future high-speed applications, new methodologies involving a higher degree of flow physics for each particular case should be developed. To achieve this goal, theoretical research in laminar-turbulent transition has traditionally focused on the study of small-amplitude disturbances undergoing exponential growth based on the mechanisms predicted by linear modal stability theory. Examples include the development of Tollmien-Schlichting waves on a flat plate or crossflow instabilities in a three-dimensional boundary layer. However, this scenario of modal growth does not always provide a satisfactory explanation for the occurrence of transition encountered in experimental and in-flight observations. Probably, the most notable example of this situation is the so-called blunt-body paradox, which makes reference to the occurrence of transition observed in blunt bodies at conditions which are stable to modal perturbations. In view of the need of an alternative explanation, nonmodal growth, also known as transient growth, has emerged as a possible mechanism leading to a different path to transition, namely, by-pass transition.19,20
The first application of the theoretical transient growth framework to compressible boundary layers is due to Hanifi et al.,5 who employed a singular value decomposition of the linearized evolution operator based on classical temporal linear stability theory (LST). They showed that a large amount of transient growth can take place over a wide range of parameter values for compressible flow, and that the optimal initial perturbations for the compressible case are found to be streamwise vortices generated by the lift-up effect. Tumin & Reshotko26 later extended the same analysis to spatially growing disturbances. Building upon this work, Reshotko & Tumin21 developed a correlation for roughness-induced transition based on spatial transient growth results which is able to reproduce transition trends observed experimentally, giving weight to the idea that transient growth could be a plausible mechanism for early transition due to distributed surface roughness. A different theoretical approach was independently developed by Andersson et al.2 and Luchini9 to quantify the spatial transient energy growth in boundary layers including non-parallel effects. This technique is based on the parabolized stability equations (PSE) and consists in an optimization marching procedure in which the disturbance energy over a given spatial extent is maximized. Recently, by employing this methodology, Paredes and
Copyright © 2019 by all authors. Published by the EUCASS association with permission.
DOI: 10.13009/EUCASS2019-591
TRANSIENT GROWTH IN A COMPRESSIBLE BOUNDARY LAYER USING 2D-LST
co-authors15 have extended the current body of results for optimal nonmodal disturbance growth in boundary layer flows to the hypersonic Mach number regime. Additionally, Paredes et al.14 have revisited the blunt-body paradox using an improved transient growth framework, leading to some modifications of the Reshotko-Tumin correlation but at the same time highlighting the need to further investigate the optimal-growth criterion underlying this correlation. Similarly, Hein et al.6 have investigated the disturbance growth in the wake of a roughness patch mounted on the forebody of a reentry capsule at Mach 5.9, for which the onset of transition is observed experimentally. Nevertheless, no evidence of modal or nonmodal growth could be found in that study. Despite these efforts, transient growth is not yet conclusively linked to the measured onset of transition over blunt-body configurations.
In the last two decades, aided by the rapid increase in computational power, classical linear stability theory has been extended to base flows featuring two and even three inhomogeneous spatial directions, giving birth to twoand three-dimensional linear stability theories, also known as BiGlobal and TriGlobal theory,24 respectively. These theories allow the study of more complex flows at an increased computational cost. Almost all of the transient growth studies available nowadays that deal with compressible basic flows are based either on classical linear stability theory or parabolized stability equations, in both cases assuming a strong inhomogenous character only in one spatial direction. Very recently, Quintanilha et al.7 have performed temporal transient growth computations on the HIFiRE-5 elliptic cone model in hypersonic flow, using two-dimensional linear stability theory (2D-LST). The results reveal a significant amount of transient growth taking place at short times, after which modal unstable perturbations take the lead and the growth becomes exponential. In the context of high-speed flows, a very interesting application of 2D-LST concerns the wake induced behind isolated roughness elements in supersonic and hypersonic flow. Such three-dimensional elements with heights comparable to the local boundary-layer thickness generally induce counter-rotating streamwise vortices in the wake flow, giving rise to low-velocity streaks that are surrounded by regions of high-detached shear. Several fundamental studies exist nowadays concerning the linear modal instability of these configurations, see for instance De Tullio et al.,4 Theiss et al.,23 Padilla Montero & Pinna.13 These analyses have revealed that the roughness wake supports the growth of sinuous and varicose instability modes that develop in the high-shear regions introduced by the counter-rotating vortex pair, and that these modes can undergo a substantial growth during the linear stages of the transition process. However, to the best of the authors’ knowledge, no studies currently exist which analyze the potential for nonmodal growth in the wake behind isolated roughness elements in high-speed flow.
The aim of this work is to highlight the particularities on the transient growth analysis of compressible boundary layers by means of two-dimensional linear stability theory and validate its implementation within the von Karman Institute stability code known as VESTA toolkit.17 This, in turn, sets the base for future transient growth computations on the wake induced by a discrete roughness element inside a hypersonic boundary layer. On first place, the mathematical formulation of the linear stability problem and the associated transient growth theoretical background are summarized. Next, some details of the current numerical implementation are discussed and a validation of the transient growth solver is presented. Finally, results are shown for 2D-LST transient growth computations in a self-similar boundary layer at Mach 2.5, emphasizing the comparison with the nonmodal growth predicted by classical linear stability theory.
2. Mathematical formulation
2.1 Governing equations
The equations governing the stability solutions obtained in this work are the Navier-Stokes equations for a compress-
ible Newtonian fluid. The primitive flow variables are denoted by density ρ, pressure p, temperature T and velocity
components ui, with u1 = u, u2 = v and u3 = w. In a Cartesian reference frame, the system of equations can be written
in conservation form as
∂ρ + ∂(ρui) = 0,
(1)
∂t ∂xi
∂(ρui) + ∂(ρuiu j) + ∂p − ∂τi j = 0, (2)
∂t
∂xj
∂xi ∂x j
∂(ρE) + ∂(ρEui + pui) + ∂qi − ∂(uiτi j) = 0, (3)
∂t
∂xi
∂xi ∂x j
where t is the time coordinate, xi is the ith spatial coordinate (x1 = x, x2 = y and x3 = z) and E = e + uiui/2 is the total energy of the fluid, with e being the specific internal energy. The viscous stress tensor τi j is defined as
∂ui ∂u j 2 ∂uk τi j = µ ∂x j + ∂xi − 3 ∂xk δi j , (4)
2
DOI: 10.13009/EUCASS2019-591
TRANSIENT GROWTH IN A COMPRESSIBLE BOUNDARY LAYER USING 2D-LST
where µ is the dynamic viscosity of the fluid and δi j is the Kronecker delta. The conductive heat flux qi is modeled
by means of Fourier’s law of heat conduction, qi = −k (∂T/∂xi), for a fluid with thermal conductivity denoted by k. A
calorically perfect gas is assumed, so the system of equations is closed using the perfect gas equation of state and two
thermodynamic relationships:
R
p = ρRT,
e = cvT and
cv
=
γ
−
, 1
(5)
where R is the specific gas constant, cv is the specific heat at constant volume and γ is the ratio of specific heats. Air
is considered with the values γ = 1.4 and R = 287 J/(kg·K). Sutherland’s law is used to model the variation of the
dynamic viscosity with temperature
T 3/2 Tre f + S µ
µ = µre f Tre f
T +Sµ ,
(6)
where S µ = 110.4 K is the Sutherland temperature constant and µre f = 1.716 × 10−5 kg/(m·s) and Tre f = 273.15 K are respectively the reference dynamic viscosity and static temperature. Regarding the thermal conductivity, its value is determined from the assumption of a constant Prandtl number Pr, that is, k = γcvµ/Pr. A value of Pr = 0.7 is employed for all the cases shown in this study.
2.2 Linear stability theory
According to classical linear stability theory, the vector of primitive flow variables q = [u, v, w, ρ, T ]T is split into a steady reference state q¯ , also known as base flow, and a small unsteady perturbation field q˜ :
q = q¯ + q˜ ,
(7)
with 1. The base flow is assumed to be locally parallel in the streamwise (x) direction, so that q¯ = q¯ (y, z) at a given
x coordinate. The perturbations are assumed to have an inhomogeneous character in the spanwise and wall-normal
directions, so that the amplitude of the perturbations is considered to be a function of both y and z coordinates. These
assumptions lead to a two-dimensional linear stability theory (2D-LST), for which the ansatz of the modal perturbations
can be written as
q˜ (x, y, z, t) = qˆ (y, z) exp[i(αx − ωt)] + c.c.,
(8)
where qˆ is the two-dimensional amplitude function, α is the wavenumber along the streamwise direction, ω is the
angular frequency and c.c. denotes the complex conjugate. For temporal transient growth analysis, the temporal sta-
bility framework is considered. The quantity α is real and represents the streamwise wavenumber of the perturbations.
On the contrary, ω is complex, with the real part ωr = {ω} being the angular frequency of q˜ and the imaginary
part ωi = {ω} its temporal growth rate. With this definition, a negative value of ωi means a temporal decay of the
amplitude function whereas a positive value implies temporal growth.
The equations governing the linear stability problem are obtained by substituting the splitting given by equa-
tion (7) together with the ansatz defined by equation (8) into the Navier-Stokes system, equations (1) to (3), and
neglecting the non-linear terms, which are of order O( 2). The resulting system of linear partial differential equations
can be written in the following compact form
Aqˆ = ωBqˆ ,
(9)
where A and B are complex and nonsymmetric differential matrix operators. After discretization of equation (9), a twodimensional algebraic generalized eigenvalue problem (GEVP) is obtained, which can be solved to obtain the angular frequency (eigenvalue) and the amplitude function (eigenvector) of each eigenmode supported by the linear stability problem.
The stability equations are made nondimensional using the boundary layer edge properties, denoted by the subscript e. The following dimensionless quantities are defined
t† = tue , x† = xi , ρ† = ρ , u† = ui , p† = p , E† = E , T † = T , µ† = µ and k† = k , (10)
le i le
ρe i ue
ρeu2e
u2e
Te
µe
ke
where le is the local Blasius length scale at the streamwise location x where the stability analysis is performed, given by le = µe x/ρeue. This choice of reference quantities leads to a set of equations where the similarity parameters are γ and the boundary layer edge Reynolds, Mach and Prandtl numbers, respectively denoted by Re, Me and Pr. In the remaining sections of this manuscript, the nondimensional equations are considered. The superscript † is dropped
unless otherwise noted.
3
DOI: 10.13009/EUCASS2019-591
TRANSIENT GROWTH IN A COMPRESSIBLE BOUNDARY LAYER USING 2D-LST
2.3 Transient growth analysis
In order to quantify the potential for transient energy growth in a given flow field satisfying the assumptions discussed
in the previous section, the orthogonality of the eigenfunctions and the size of the perturbations obtained by linear
stability theory have to be evaluated. The first can be achieved by the definition of an inner product and the second by an associated energy norm. In this work, the disturbance energy norm introduced by Chu,3 Mack10 and later rederived by Hanifi et al.5 is adopted, given by
E(t) = 1 ρ¯ u˜2 + v˜2 + w˜ 2 + T¯ ρ˜2 +
ρ¯
T˜ 2 dxdydz
(11)
2V
γρ¯ Me2 γ (γ − 1) T¯ Me2
in a volume V. For transient growth analysis, it is convenient to express the disturbance energy in terms of the timedependent amplitude functions, defined as qˇ (y, z, t) = qˆ (y, z) exp(−iωt), so that the ansatz (equation (8)) becomes
q˜ (x, y, z, t) = qˇ (y, z, t) exp(iαx) + c.c.
(12)
Additionally, since the perturbations are assumed to be periodic along the streamwise direction, the disturbance energy can be integrated over one wavelength in x without loss of generality, see Weder et al.28 Substituting equation (12) into equation (11) and integrating along x between 0 and 2π/α, the following expression is obtained
E(t) = 2π ρ¯ (uˇ∗uˇ + vˇ∗vˇ + wˇ ∗wˇ ) + T¯ ρˇ∗ρˇ +
ρ¯
Tˇ ∗Tˇ dydz,
(13)
αΩ
γρ¯ Me2
γ (γ − 1) T¯ Me2
where ∗ denotes the complex conjugate and Ω is the space in which the two-dimensional eigenfunctions are defined.
equation (13) can be rewritten in a compact form as
2απ E(t) = ||qˇ ||2E = (qˇ , qˇ )E , (14)
where the subscript E denotes the energy norm. The associated inner product for a given pair of amplitude functions is then expressed as
(qˇ k, qˇ l)E = qˇ kHMqˇ l dydz,
(15)
Ω
where the superscript H refers to the complex conjugate transpose and the matrix M is given by
T¯
ρ¯
M = diag ρ¯, ρ¯, ρ¯, γρ¯ M2 , γ (γ − 1) T¯ M2 .
(16)
e
e
After a choice on the inner product and the associated energy norm has been made, the temporal transient growth analysis requires an eigenvector expansion of the time dependent amplitude function. Denoting the first K eigenvalues and eigenvectors obtained from the solution of the eigenvalue problem given by equation (9) as ωk and qˆ k, respectively, the vector qˇ can be expanded as
K
qˇ (y, z, t) = κk(t)qˆ k(y, z)
(17)
k=1
with
κk(t) = κk(0) exp(−iωkt),
(18)
or in matrix form
qˇ (y, z, t) = Qκ(t) with κ(t) = Λκ(0) and Λ = diag exp(−iω1t), ..., exp(−iωKt) ,
(19)
where Q is a matrix whose columns contain the eigenvectors qˆ k and the vector κ(0) = [κ1(0), ..., κK(0)]T contains the initial expansion coefficients that determine the individual contribution of each mode to the transient response. By employing the inner product defined by equation (15) together with the expansion in equation (17), the following
relationship can be found between the energy norm and the L2-norm
(qˇ k, qˇ l)E = (Qκk, Qκl)E = κk∗QHMQκl dydz = κk∗Dκl = κk∗FHFκl = (Fκk, Fκl)2 ,
(20)
Ω
4
DOI: 10.13009/EUCASS2019-591
TRANSIENT GROWTH IN A COMPRESSIBLE BOUNDARY LAYER USING 2D-LST
where the subscript 2 denotes the L2-norm. The matrix F is the Cholesky decomposition of the Hermitian and positive definite matrix D, that is, D = FHF, which is defined element-wise as
Dkl = qˆ kHMqˆ l dydz.
(21)
Ω
The relation given by equation (20) allows to rewrite the disturbance energy norm in terms of the L2-norm as follows
||qˇ (t)||2E = ||Fκ(t)||22.
(22)
Finally, the maximum transient energy growth G(t), optimized for all the possible initial conditions, is given by (see Schmid & Henningson22)
G(t) = max ||qˇ (t)||2E = max ||Fκ(t)||22 = max ||FΛF−1Fκ(0)||22 = ||FΛF−1||2.
(23)
||qˇ (0)||2E ||Fκ(0)||22 ||Fκ(0)||22 2
The L2-norm of matrix S = FΛF−1 can be computed by means of singular value decomposition (SVD). The value of G(t) is then given by the square of the largest singular value of matrix S. The optimal initial conditions giving
rise to the maximum transient growth can be computed via the right singular eigenvector r associated to the largest singular value of matrix S, that is, κ(0) = F−1r.
3. Numerical methodology
The numerical solution of the linear stability problem and the posterior transient growth analysis are performed using the von Karman Institute Extensible Stability & Transition Analysis (VESTA) toolkit, originally developed by Pinna.16,17 The discretization of the partial differential eigenvalue problem is performed employing the Chebyshev collocation method.25 This technique is based on a Lagrange polynomial interpolation constructed in a structured grid with a non-uniform point distribution given by the Chebyshev-Gauss-Lobatto collocation points. These points are defined on a transformed domain with spanwise and wall-normal coordinates respectively denoted by ξ and η, with ξ, η ∈ [−1, 1]. This transformed space does not coincide with the physical domain of interest. Therefore, geometrical mappings have to be used. In the computations presented in this work, the mapping originally introduced by Malik11 is applied in the wall-normal direction, which allows placing half of the collocation points below a given coordinate in each direction, whereas a uniform mapping is used in the spanwise direction, which simply translates the transformed domain to the physical one. In this way, the grid used for solving the stability problem can be clustered towards the boundary layer. On the other hand, despite this process transforms the stability grid into the physical domain, the mapped grid defined by the Chebyshev-Gauss-Lobatto collocation points is different than the mesh employed to obtain the base flow solution. As a consequence, before solving the eigenvalue problem, the base flow data is interpolated on the collocation grid by means of a cubic spline interpolation in each spatial direction.
The boundary conditions imposed on the perturbation quantities are as follows. At the wall, the velocity and the temperature perturbations are set to zero via a homogeneous Dirichlet condition, whereas the density fluctuation is determined by means of a compatibility condition. In the wall-normal far-field boundary, the perturbations are forced to decay by imposing a Dirichlet boundary condition as well, with the exception of ρˆ, which once again satisfies a compatibility condition. Regarding the spanwise boundaries, periodic boundary conditions are specified. This implies that the value of each disturbance quantity and its first derivative normal to the boundary are set to be equal at both spanwise boundaries.
For the transient growth analysis to be physically meaningful, a large number of eigenmodes must be included in the construction of the matrix S, see for instance Hanifi et al.5 In this work, the QZ algorithm12 (also known as generalized Schur decomposition) is used to compute the full spectrum of eigenvalues and eigenvectors associated to the discrete generalized eigenvalue problem, employing subroutines belonging to the LAPACK1 library. The integral in equation (21) is evaluated using the Chebyshev integration weights derived in Hanifi et al.5 This ensures that the order of accuracy of the discretization is preserved after integration. The calculation of the largest singular value of matrix S is also performed using the LAPACK library.
4. Validation
In order to validate the implementation of the transient growth analysis methodology within the von Karman Institute (VKI) stability code (VESTA), two test cases investigated by Hanifi et al.5 are reproduced here based on classical
5
DOI: 10.13009/EUCASS2019-591
TRANSIENT GROWTH IN A COMPRESSIBLE BOUNDARY LAYER USING 2D-LST
Figure 1: Temporal evolution of transient growth for a self-similar boundary layer at the conditions described in section 4. Comparison of the results obtained between VESTA and the reference5 data for one-dimensional linear stability theory.
linear stability (LST). Both cases employ a base flow obtained by solving the compressible self-similar boundary layer equations, see for instance Anderson.8 These equations are solved as well within VESTA by means of a fourth-order Runge-Kutta method. The self-similar profiles considered for validation have a boundary layer edge Mach number of Me = 2.5, a total temperature of T0 = 333 K and assume an adiabatic wall. The first case evaluates the self-similar profile at a Reynolds number of Re = 300, with a streamwise wavenumber α = 0 and a spanwise wavenumber β = 0.1. The second configuration assumes a Reynolds number of Re = 3000 and wavenumbers α = 0.06 and β = 0.1.
Figure 1 shows the energy growth as a function of time for each of the two cases as computed by the current solver. A very good agreement with the reference results is obtained. At the conditions of the first test case, the boundary layer is stable to linear modal perturbations. The energy reaches a maximum transient growth at a finite time and then progressively decays to zero. On the contrary, the second configuration features a linearly unstable mode, and as a result the energy grows exponentially as time tends to infinity. Nevertheless, for short times the disturbance energy still undergoes a significant transient growth, showing a faster increase than the exponential growth associated to the unstable mode.
5. Transient growth results based on two-dimensional linear stability theory
This section presents results obtained for transient growth computations using eigenmodes solved by means of twodimensional linear stability theory and their comparison against the one-dimensional counterpart. The base flows considered in this analysis have conditions that are similar to the validation test cases shown in the previous section, namely, a self-similar boundary layer at Me = 2.5, with a total temperature T0 = 333 K and an adiabatic wall. For comparison purposes between LST and 2D-LST results, the base flow wall-normal velocity is set to zero in all the 2D-LST computations.
5.1 Case 1: Re = 3000, α = 0.06
The first case considered corresponds to the same Reynolds number and streamwise wavenumber as for the second validation test case, that is, Re = 3000 and α = 0.06. Before proceeding further, it is important to note that for the 2DLST computations a choice has to be made regarding the spanwise size of the stability domain, denoted by Lz. Together with the use of periodic boundary conditions, this size determines the specific wave-numbers that are resolved by the computation. For a given value of Lz, only modes with spanwise wanumbers that are smaller multiples of, or equal to, 2π/Lz can be computed. Additionally, the number of grid points in the spanwise direction, denoted by Nz, limits the number of spanwise multiple modes that can be retrieved for each LST mode. For a given discretization, Nz − 2 multiple modes are computed. For instance, if a value of Lz = 2π/0.1 is chosen, this means that, in addition to all the modes obtained by the classical LST operator for β = 0.1, multiple modes with spanwise wavenumbers β = n0.1, with n = 2, ..., Nz − 2, are also found in the 2D-LST spectrum.
Figure 2 displays the transient growth results obtained for the current case by employing two different spanwise domain lengths, namely, Lz = 2π/0.1 and Lz = 2π/0.05 and a grid resolution of Nz × Ny = 61 × 61 points. LST results
6
DOI: 10.13009/EUCASS2019-591
TRANSIENT GROWTH IN A COMPRESSIBLE BOUNDARY LAYER USING 2D-LST
Figure 2: Comparison of transient growth results between LST and 2D-LST for two different spanwise domain lengths at the conditions of case 1.
for different spanwise wavenumbers are also shown for comparison, which have been computed using 151 collocation points in the wall-normal direction. As it can be observed, the transient energy growth predicted by means of 2D-LST is equal to the envelope of the transient growth curves associated to each of the LST computations at the spanwise wavenumbers that can be resolved by the two-dimensional discretization. In other words, when the domain length of Lz = 2π/0.1 is used, the 2D-LST transient growth is the envelope of the LST transient growth curves for β = n0.1, with n = 1, ..., 59 in this case. The four multiple spanwise wavenumbers for which the LST operator undergoes the highest transient growth are also included in the left plot of Figure 2 to illustrate this fact. As it can be seen, a linearly unstable wave is only present at β = 0.1 for this case, so the 2D-LST energy growth converges to the exponential growth given by the growth rate of this particular mode in the limit of large times. However, for smaller times, the other wavenumbers exhibit a larger nonmodal growth, and as a result the 2D-LST curve also undergoes a larger growth accordingly. It is important to note that when a domain length twice as big is considered, the spanwise wavenumbers that can be resolved are different. Now, all modes with β = n0.05 are included, and this results in a different envelope mainly due to the transient growth behaviour of the LST operator at β = 0.15, for which the least stable mode has a growth rate that is very close to 0 and the energy decays very slowly.
To better illustrate the relationship between the one-dimensional and two-dimensional linear stability analyses, Figure 3 presents a comparison of the LST and the 2D-LST spectra obtained for the current case. Figure 3a) shows an overview of the relevant part of the temporal spectrum as computed with LST for β = 0.1 and Ny = 81 points and with 2D-LST using Nz = 31, Ny = 81 points and Lz = 2π/0.1. The spectrum contains a vertical continuous branch at ωr = α = 0.06 which consists of entropy and vorticity waves and two horizontal continuous branches which contain acoustic waves. An additional group of modes is resolved in the positive imaginary part on top of the vertical continuous branch. These are numerical spurious modes that appear as part of the numerical discretization. The discrete physical eigenmodes are found in the diagonal-like structures that emerge from the vertical continuous branch. A single discrete unstable is found at these conditions, which corresponds to a first-mode disturbance. All the modes present in this plot are included in the transient growth calculation for both the LST and the 2D-LST theories, with the exception of the spurious numerical modes. Figure 3b) displays a closer view of the same spectra plotted in a) where the difference between the one and two-dimensional computations can be clearly appreciated. As it can be noticed, for each mode appearing in the LST spectrum, there are Nz − 2 multiple modes resolved in the 2D-LST spectrum. Figure 3c) is included as a proof of convergence when the number of discretization points in the spanwise direction is doubled. Finally, Figure 3d) illustrates the comparison between the 2D-LST spectra obtained for both spanwise domain sizes considered. When a domain twice as big is employed, additional discrete modes can be found in the spectrum, respectively corresponding to the spanwise wavenumbers β = 0.05, 0.15, 0.25 etc., which are not resolved by the other domain length.
When performing 2D-LST transient growth analyses of more complex boundary layer flows, care should then be taken to choose an appropriate domain size. Ideally, a previous analysis based on LST should be performed to identify the spanwise wavenumbers that produce the higher contributions to the transient growth envelope.
7
DOI: 10.13009/EUCASS2019-591
TRANSIENT GROWTH IN A COMPRESSIBLE BOUNDARY LAYER USING 2D-LST
Figure 3: Comparison of the temporal stability spectrum between LST and 2D-LST for the conditions analyzed in case 1.
5.2 Case 2: Re = 300, α = 0 A second case is also examined which corresponds to conditions at which the optimal transient growth is found for a Reynolds number of Re = 300. From LST-based computations, see Hanifi et al.,5 the optimal growth is found for α = 0 and β = 0.235. A comparison of the transient growth results obtained for this case is depicted in Figure 4. In this case, the LST transient growth evolution for β = 0.235 dominates for all times. As a result, the 2D-LST curve when Lz = 2π/0.235 is identical to the LST one. The advantage of this case is that it can be used to compare the optimal disturbances computed with both theories as well. Figure 5 shows the optimal perturbations evaluated at t = tmax, where tmax denotes the time at which maximum G is achieved, which for this case corresponds to t = 980. Reference values for the streamwise velocity and temperature perturbations from Hanifi et al.5 are also added for validation, for which a very good agreement is found with the current implementation. The initial optimal disturbance for this case corresponds to streamwise vortices that develop into streamwise streaks as time evolves. The amplitude functions of the optimal perturbations shown here correspond to the signature of these streak structures. It is interesting to note that the relative magnitude of the wall-normal and spanwise velocity perturbations is much smaller compared to the ones associated to the streamwise velocity and the temperature. Similarly, Figure 6 presents the optimal perturbations computed by means of the transient growth computation based on 2D-LST modes. As it can be observed, the shape of the optimal disturbances as well as their relative magnitude coincides with the results obtained with LST. Additionally, they feature a spanwise wavelength equal to the size of the domain for the spanwise direction, which means that β = 0.235 also for the optimal perturbations computed with the two-dimensional eigenmodes.
The results shown in this section constitute a validation of the current transient growth implementation based on two-dimensional linear stability theory. 8
DOI: 10.13009/EUCASS2019-591
TRANSIENT GROWTH IN A COMPRESSIBLE BOUNDARY LAYER USING 2D-LST
Figure 4: Comparison of transient growth results between LST and 2D-LST for the conditions of case 2.
Figure 5: Optimal disturbances at t = tmax computed with LST for the conditions studied in case 2.
6. Conclusions
Transient growth computations based on two-dimensional linear stability theory have been performed for a self-similar boundary layer at Mach 2.5 by means of singular value decomposition of the linearized temporal stability operator. Two different Reynolds numbers and spanwise wavenumbers have been considered, resulting first in a case where the boundary layer is linearly unstable to modal disturbances but undergoes a significant amount of transient growth for small times, and second in a linearly stable base flow at the conditions for which the optimal transient growth is encountered. It has been found that the transient energy growth predicted by two-dimensional linear stability theory is equal to the envelope of the nonmodal classical linear stability results evaluated at all the possible spanwise wavenumbers that can be resolved by a given spanwise discretization of the two-dimensional domain. Care should be taken to size the two-dimensional stability domain appropriately in the spanwise direction in order to resolve all the modes with a wave-like behavior in the spanwise direction that have the largest contribution to the energy growth.
The optimal disturbances have been computed for the second studied case at the time instant in which the maximum nonmodal growth is present. A comparison of the optimal perturbations computed between LST and 2D-LST theories shows that the two-dimensional disturbances have the same shape and magnitude as the one-dimensional ones and feature a spanwise wavenumber equal to the optimal one as well. These results constitute a validation of the current transient growth implementation for two-dimensional eigenmodes, and sets the base for the application of this methodology to more complex base flows with two-inhomogenous spatial directions, such as the wake induced behind a discrete roughness element inside a high-speed boundary layer.
9
DOI: 10.13009/EUCASS2019-591
TRANSIENT GROWTH IN A COMPRESSIBLE BOUNDARY LAYER USING 2D-LST
Figure 6: Optimal disturbances at t = tmax computed with 2D-LST for the conditions studied in case 2.
Acknowledgments
This work is part of a project that has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 675008.
References
[1] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen. LAPACK Users’ Guide. Society for Industrial and Applied Mathematics, Philadelphia, PA, third edition, 1999.
[2] P. Andersson, M. Berggren, and D. S. Henningson. Optimal disturbances and bypass transition in boundary layers. Phys. Fluids, 11(1):134–150, 1999.
[3] B. Chu. On the energy transfer to small disturbances in fluid flow (Part 1). Acta Mech., 1(3):215–234, 1965. [4] N. De Tullio, P. Paredes, N. D. Sandham, and V. Theofilis. Laminar-turbulent transition induced by a discrete
roughness element in a supersonic boundary layer. J. Fluid Mech., 735:613–646, 2013. [5] A. Hanifi, P. J. Schmid, and D. S. Henningson. Transient growth in compressible boundary layer flow. Phys.
Fluids, 8(3):826–837, 1996. [6] S. Hein, A. Theiss, A. Di Giovanni, C. Stemmer, T. Schilden, W. Schröder, P. Paredes, M. M. Choudhari, F. Li,
and E. Reshotko. Numerical investigation of roughness effects on transition on spherical capsules. J. Spacecraft Rockets, 56(2):388–404, 2019. [7] H. R. Quintanilha Jr., V. Theofilis, and A. Hanifi. Global transient-growth analysis of hypersonic flow on the HIFiRE-5 elliptic cone model. In 2019 AIAA Aerospace Sciences Meeting, AIAA SciTech Forum, San Diego, CA, USA, 7-11 January 2019, AIAA 2019-2148. [8] J. D. Anderson Jr. Hypersonic and high-temperature gas dynamics. American Institute of Aeronautics and Astronautics, second edition, 2006. [9] P. Luchini. Reynolds-number-independent instability of the boundary layer over a flat surface: optimal perturbations. J. Fluid Mech., 404:289–309, 2000. [10] L. M. Mack. Boundary layer stability theory. Technical report, Jet Propulsion Laboratory Report 900-277, 1969. [11] M. R. Malik. Numerical methods for hypersonic boundary layer stability. J. Comput. Phys., 86(2):376–413, 1990. [12] C. B. Moler and G. W. Stewart. An algorithm for generalized matrix eigenvalue problems. SIAM J. Numer. Anal., 10(2):241–256, 1973.
10