# 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 aﬀect 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 coeﬃcient on a ﬂat plate is estimated to be an order of magnitude higher than the corresponding laminar value.27 Despite signiﬁcant research eﬀorts 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 ﬂow 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 ﬂat plate or crossﬂow 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-ﬂight 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 diﬀerent path to transition, namely, by-pass transition.19,20

The ﬁrst application of the theoretical transient growth framework to compressible boundary layers is due to Haniﬁ 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 ﬂow, and that the optimal initial perturbations for the compressible case are found to be streamwise vortices generated by the lift-up eﬀect. 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 diﬀerent 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 eﬀects. 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 ﬂows 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 modiﬁcations 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 eﬀorts, transient growth is not yet conclusively linked to the measured onset of transition over blunt-body conﬁgurations.

In the last two decades, aided by the rapid increase in computational power, classical linear stability theory has been extended to base ﬂows 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 ﬂows at an increased computational cost. Almost all of the transient growth studies available nowadays that deal with compressible basic ﬂows 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 ﬂow, using two-dimensional linear stability theory (2D-LST). The results reveal a signiﬁcant 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 ﬂows, a very interesting application of 2D-LST concerns the wake induced behind isolated roughness elements in supersonic and hypersonic ﬂow. Such three-dimensional elements with heights comparable to the local boundary-layer thickness generally induce counter-rotating streamwise vortices in the wake ﬂow, 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 conﬁgurations, 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 ﬂow.

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 ﬁrst 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 ﬂuid. The primitive ﬂow 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 ﬂuid, with e being the speciﬁc internal energy. The viscous stress tensor τi j is deﬁned 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 ﬂuid and δi j is the Kronecker delta. The conductive heat ﬂux qi is modeled

by means of Fourier’s law of heat conduction, qi = −k (∂T/∂xi), for a ﬂuid 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 speciﬁc gas constant, cv is the speciﬁc heat at constant volume and γ is the ratio of speciﬁc 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 ﬂow variables q = [u, v, w, ρ, T ]T is split into a steady reference state q¯ , also known as base ﬂow, and a small unsteady perturbation ﬁeld q˜ :

q = q¯ + q˜ ,

(7)

with 1. The base ﬂow 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 deﬁnition, 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 deﬁned 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 diﬀerential equations

can be written in the following compact form

Aqˆ = ωBqˆ ,

(9)

where A and B are complex and nonsymmetric diﬀerential 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 deﬁned

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 ﬂow ﬁeld 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 ﬁrst can be achieved by the deﬁnition 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 Haniﬁ 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, deﬁned 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 deﬁned.

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 ﬁrst 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 coeﬃcients that determine the individual contribution of each mode to the transient response. By employing the inner product deﬁned 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 deﬁnite matrix D, that is, D = FHF, which is deﬁned 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 diﬀerential 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 deﬁned 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 deﬁned by the Chebyshev-Gauss-Lobatto collocation points is diﬀerent than the mesh employed to obtain the base ﬂow solution. As a consequence, before solving the eigenvalue problem, the base ﬂow 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 ﬂuctuation is determined by means of a compatibility condition. In the wall-normal far-ﬁeld boundary, the perturbations are forced to decay by imposing a Dirichlet boundary condition as well, with the exception of ρˆ, which once again satisﬁes a compatibility condition. Regarding the spanwise boundaries, periodic boundary conditions are speciﬁed. This implies that the value of each disturbance quantity and its ﬁrst 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 Haniﬁ 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 Haniﬁ 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 Haniﬁ 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 ﬂow 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 proﬁles 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 ﬁrst case evaluates the self-similar proﬁle at a Reynolds number of Re = 300, with a streamwise wavenumber α = 0 and a spanwise wavenumber β = 0.1. The second conﬁguration 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 ﬁrst test case, the boundary layer is stable to linear modal perturbations. The energy reaches a maximum transient growth at a ﬁnite time and then progressively decays to zero. On the contrary, the second conﬁguration features a linearly unstable mode, and as a result the energy grows exponentially as time tends to inﬁnity. Nevertheless, for short times the disturbance energy still undergoes a signiﬁcant 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 ﬂows 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 ﬂow wall-normal velocity is set to zero in all the 2D-LST computations.

5.1 Case 1: Re = 3000, α = 0.06

The ﬁrst 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 speciﬁc 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 diﬀerent 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 diﬀerent spanwise domain lengths at the conditions of case 1.

for diﬀerent 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 diﬀerent. Now, all modes with β = n0.05 are included, and this results in a diﬀerent 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 ﬁrst-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 diﬀerence 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 ﬂows, 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 Haniﬁ 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 Haniﬁ 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 diﬀerent Reynolds numbers and spanwise wavenumbers have been considered, resulting ﬁrst in a case where the boundary layer is linearly unstable to modal disturbances but undergoes a signiﬁcant amount of transient growth for small times, and second in a linearly stable base ﬂow 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 ﬂows 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 ﬂuid ﬂow (Part 1). Acta Mech., 1(3):215–234, 1965. [4] N. De Tullio, P. Paredes, N. D. Sandham, and V. Theoﬁlis. Laminar-turbulent transition induced by a discrete

roughness element in a supersonic boundary layer. J. Fluid Mech., 735:613–646, 2013. [5] A. Haniﬁ, P. J. Schmid, and D. S. Henningson. Transient growth in compressible boundary layer ﬂow. 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 eﬀects on transition on spherical capsules. J. Spacecraft Rockets, 56(2):388–404, 2019. [7] H. R. Quintanilha Jr., V. Theoﬁlis, and A. Haniﬁ. Global transient-growth analysis of hypersonic ﬂow 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 ﬂat 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 aﬀect 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 coeﬃcient on a ﬂat plate is estimated to be an order of magnitude higher than the corresponding laminar value.27 Despite signiﬁcant research eﬀorts 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 ﬂow 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 ﬂat plate or crossﬂow 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-ﬂight 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 diﬀerent path to transition, namely, by-pass transition.19,20

The ﬁrst application of the theoretical transient growth framework to compressible boundary layers is due to Haniﬁ 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 ﬂow, and that the optimal initial perturbations for the compressible case are found to be streamwise vortices generated by the lift-up eﬀect. 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 diﬀerent 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 eﬀects. 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 ﬂows 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 modiﬁcations 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 eﬀorts, transient growth is not yet conclusively linked to the measured onset of transition over blunt-body conﬁgurations.

In the last two decades, aided by the rapid increase in computational power, classical linear stability theory has been extended to base ﬂows 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 ﬂows at an increased computational cost. Almost all of the transient growth studies available nowadays that deal with compressible basic ﬂows 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 ﬂow, using two-dimensional linear stability theory (2D-LST). The results reveal a signiﬁcant 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 ﬂows, a very interesting application of 2D-LST concerns the wake induced behind isolated roughness elements in supersonic and hypersonic ﬂow. Such three-dimensional elements with heights comparable to the local boundary-layer thickness generally induce counter-rotating streamwise vortices in the wake ﬂow, 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 conﬁgurations, 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 ﬂow.

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 ﬁrst 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 ﬂuid. The primitive ﬂow 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 ﬂuid, with e being the speciﬁc internal energy. The viscous stress tensor τi j is deﬁned 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 ﬂuid and δi j is the Kronecker delta. The conductive heat ﬂux qi is modeled

by means of Fourier’s law of heat conduction, qi = −k (∂T/∂xi), for a ﬂuid 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 speciﬁc gas constant, cv is the speciﬁc heat at constant volume and γ is the ratio of speciﬁc 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 ﬂow variables q = [u, v, w, ρ, T ]T is split into a steady reference state q¯ , also known as base ﬂow, and a small unsteady perturbation ﬁeld q˜ :

q = q¯ + q˜ ,

(7)

with 1. The base ﬂow 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 deﬁnition, 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 deﬁned 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 diﬀerential equations

can be written in the following compact form

Aqˆ = ωBqˆ ,

(9)

where A and B are complex and nonsymmetric diﬀerential 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 deﬁned

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 ﬂow ﬁeld 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 ﬁrst can be achieved by the deﬁnition 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 Haniﬁ 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, deﬁned 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 deﬁned.

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 ﬁrst 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 coeﬃcients that determine the individual contribution of each mode to the transient response. By employing the inner product deﬁned 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 deﬁnite matrix D, that is, D = FHF, which is deﬁned 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 diﬀerential 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 deﬁned 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 deﬁned by the Chebyshev-Gauss-Lobatto collocation points is diﬀerent than the mesh employed to obtain the base ﬂow solution. As a consequence, before solving the eigenvalue problem, the base ﬂow 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 ﬂuctuation is determined by means of a compatibility condition. In the wall-normal far-ﬁeld boundary, the perturbations are forced to decay by imposing a Dirichlet boundary condition as well, with the exception of ρˆ, which once again satisﬁes a compatibility condition. Regarding the spanwise boundaries, periodic boundary conditions are speciﬁed. This implies that the value of each disturbance quantity and its ﬁrst 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 Haniﬁ 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 Haniﬁ 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 Haniﬁ 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 ﬂow 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 proﬁles 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 ﬁrst case evaluates the self-similar proﬁle at a Reynolds number of Re = 300, with a streamwise wavenumber α = 0 and a spanwise wavenumber β = 0.1. The second conﬁguration 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 ﬁrst test case, the boundary layer is stable to linear modal perturbations. The energy reaches a maximum transient growth at a ﬁnite time and then progressively decays to zero. On the contrary, the second conﬁguration features a linearly unstable mode, and as a result the energy grows exponentially as time tends to inﬁnity. Nevertheless, for short times the disturbance energy still undergoes a signiﬁcant 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 ﬂows 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 ﬂow wall-normal velocity is set to zero in all the 2D-LST computations.

5.1 Case 1: Re = 3000, α = 0.06

The ﬁrst 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 speciﬁc 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 diﬀerent 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 diﬀerent spanwise domain lengths at the conditions of case 1.

for diﬀerent 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 diﬀerent. Now, all modes with β = n0.05 are included, and this results in a diﬀerent 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 ﬁrst-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 diﬀerence 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 ﬂows, 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 Haniﬁ 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 Haniﬁ 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 diﬀerent Reynolds numbers and spanwise wavenumbers have been considered, resulting ﬁrst in a case where the boundary layer is linearly unstable to modal disturbances but undergoes a signiﬁcant amount of transient growth for small times, and second in a linearly stable base ﬂow 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 ﬂows 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 ﬂuid ﬂow (Part 1). Acta Mech., 1(3):215–234, 1965. [4] N. De Tullio, P. Paredes, N. D. Sandham, and V. Theoﬁlis. Laminar-turbulent transition induced by a discrete

roughness element in a supersonic boundary layer. J. Fluid Mech., 735:613–646, 2013. [5] A. Haniﬁ, P. J. Schmid, and D. S. Henningson. Transient growth in compressible boundary layer ﬂow. 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 eﬀects on transition on spherical capsules. J. Spacecraft Rockets, 56(2):388–404, 2019. [7] H. R. Quintanilha Jr., V. Theoﬁlis, and A. Haniﬁ. Global transient-growth analysis of hypersonic ﬂow 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 ﬂat 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