ABSTRACT

Introduction
Boundary layer (BL) transition can occur via different mechanisms. If external disturbances, such as incoming tur- bulence, are small, modal mechanisms, in the form of grow- ing instability waves are usually responsible for transition. In incompressible flows, the only unstable mode is the Tollmien- Schlichting (TS) wave. In supersonic flows, at sufficiently high Mach number, there is also a family of unstable modes, under- stood as acoustic waves trapped in the BL, which were first identified by Mack (1963).  The receptivity to acoustic dis- turbances has been studied both for TS (Goldstein, 1985) and Mack waves (Fedorov & Khokhlov, 1992). If larger disturbances are present, non-modal mecha- nisms are often responsible for triggering transition.  They are characterized by a transient amplification associated with a linear combination of stable modes, a phenomena known as transient-growth, and which has been characterized for both incompressible (Andersson et al. , 1999) and compress- ible (Hanifi et al. , 1996) flows. Which of the two, modal or non-modal mechanisms, are responsible for transition depends on the amplitude of the in- cident disturbances and the receptivity of each mechanism. There are several mechanisms by which free-stream distur- bances enter the boundary layer. For low turbulence levels, the dominant one is likely localized in the leading edge (Brandt et al. , 2004).  The receptivity magnitude will thus provide the initial amplitude of waves involved in modal and non- modal mechanisms.  It is also important for scenarios where transition is observed at the leading edge.  Due to the strong favourable pressure gradient, which stabilizes TS waves, and the relatively low local Mach number, due to slower veloci- ties and higher temperatures near the stagnation point, transi- tion at the leading edge is believed to be caused by non-modal mechanisms (Paredes et al. , 2017).  Most studies considering this possibility have been mainly concerned with transition in- duced by wall roughness. Aiming to complement previous studies that focused on the receptivity of TS, Mack, and cross-flow instability waves, we here focus on the receptivity of streaks in blunt bodies in non-swept configurations, using a paraboloid body as a test case.  A large range of the parameter space was studied, in- cluding varying Reynolds and Mach numbers.  A description of the method used, results obtained, and their discussion are presented next. Numerical method The linearized Navier-Stokes equations can be compactly written as,
(1) where q is a vector containing the flow density, velocities, and temperature, and f contain forcing terms. Assuming cylindri- cal symmetry and using the ansatz q(x,t) = q(z,r)e-icx+mq (2) where r , φ , and z are the cylindrical coordinates, (1) is written as 1 where A m is the compressible linearized Navier-Stokes opera- tor for the m -th azimuthal component.  An input-output relation is obtained re-writing (3) as
where R m ( ω ) = ( iω − A ) − 1 is the resolvent operator. A forcing term is said to be optimal if it maximizes the Rayleigh quotient
where ˆ q and ˆ f satisfy (4). The value σ is usually referred as the gain associated with the forcing/response pair. Suboptimal forcing components maximize the same ratio, with the added constraint of being orthogonal to the leading forcing mode. Gains and forcing/response mode pairs are typically ob- tained from a singular value decomposition (SVD) of the weighted resolvent operator R m ′ (ω) = W ¹ / 2 R m ′ (ω)W − 1/2 , where W is a matrix representing an energy norm and, for non-uniform meshes, integration weights. However, for two- dimensional systems, the matrix inversion in (3) and the SVD of R ′ m become prohibitively costly. Instead, the leading singu- lar values and modes are obtained using an iterative algorithm, such as the power iteration or the Arnoldi algorithms, based on solutions of the linear problem (3). These solutions are typi- cally obtained using an LU decomposition of R m to reduce the overall cost. We use a code developed by Schmidt & Rist (2011) to study the modal stability of streamwise corner-flows, which was later extended to non-modal stability analysis (Schmidt & Rist, 2014), and resolvent analysis in Schmidt et al . (2018).  An important difference with respect to the original code is, that here, the geometry cannot be described by a cartesian mesh, for which the use of finite-differences schemes is most natu- ral. Instead, here an original rectangular domain, described by the coordinates ξ i , is mapped onto the desired geometry, re- sulting in curvilinear coordinates x j = x j (ξ j ). The mapping is illustrated in figure 1. In order to model the system in the physical domain ( x j ), differentiation matrices need to be constructed based on the finite differences used in the original domain ( ξ i ). First order derivatives with respect to x j are obtained from the derivatives in the original domain as
where J ij = ∂x j /∂ξ j represents the mapping Jacobian, g is a scalar field, and the Einstein summation convention is used. The Jacobian J is computed numerically and then inverted to obtain J-1. Second derivatives can be obtained by re-applying (6). However, when using this approach the effective stencil for second derivatives doubles, which is an undesirable feature when used to numerically construct the operator A m , as it be- comes less sparse, increasing the computational cost. Instead, second derivatives are computed as is computed using the matrix identity
The term ∂ x j J contains information about the second derivatives of the mapping, and reads
This term is also computed numerically. With this formulation, the effects of the mapping is com- puted explicitly, and the differentiation matrix retains its spar- sity. Throughout this work, a 4-th order differentiation scheme is used to numerically compute the derivatives, using a cen- tered scheme in the interior the domain and biased stencils on the boundaries.
Geometry and baseflow
The body surface is given by
The free-stream flow moves towards +z with unit velocity. The problem is made non-dimensional with respect to the freestream density, viscosity and velocity and the diameter of curvature of the leading edge. The inviscid and incompressible flow around the body is obtained analytically using parabolic coordinates, (σ,τ,φ), given by
The body surface corresponds to σ = cylindrical coordinates are recovered as z = σ 2 2τσ . The flow potential function, ψ , is given by the solution of the Laplace equation, netration boundary condition at the body surface, for σ = σ 0 , and the requirement that far from the body to flow corresponds to a uniform right-moving flow, i.e. ψ = z = (σ 2 − τ 2 ) at σ → ∞. The solution is given by
From (15), the influence of the paraboloid body on the flow scales decays roughly with the logarithm of the wall distance.  Accurately capturing this effect numerically would require a large computational domain if uniform free-stream boundary conditions were used as inflow of a viscous incom- pressible simulation domain. To reduce the required computa- tional domain needed, the inviscid solution was thus used as a boundary condition for the viscous problem. Viscous incompressible baseflows were obtained by time- marching the incompressible Navier-stokes equations in time until the norm of the velocity time-derivative becomes smaller than 10-8. No-slip boundary conditions were applied at the body surface and outflow conditions on the right-most edge of the domain. Integration was carried out using the the axisym- metric formulation of the spectral code Nek5000 (Fischer & Lottes, 2005).  The inviscid and the Re = 100 baseflows are shown in figure 2. The baseflows for higher values of Re look similar, but with thinner boundary layers.  The computational domain and mesh discretization used to obtain the baseflows for each Reynolds number were chosen to guarantee a domain considerably larger than the boundary layer around the body and to have sufficient grid points inside it. The boundary layer displacement thickness (δ*) and boundary layer edge velocities ( U ε ) as a function of the dis- tance to the leading edge and the Reynolds number is pre- sented in figure 3. The BL thickness for Re > 5000 is seen to scale with Re−1/2 , as typically observed for flat plates. The solution at the boundary layer edge closely follows the invis- cid solution at the surface and is nearly independent of the Reynolds number
Reynolds number effects
To study the effect of the Reynolds number on the lead- ing forcing/response modes close to the stagnation point, the physical domain was initially limited downstream to z = 1. As will be soon reported, the forcing mode extends further upstream for increasing Reynolds numbers, and thus the up- stream boundary has been moved further from the body for larger values of Re . An extended domain, truncated at z = 2, was used to verify the impact of the domain truncation. Figure 4a shows the gains of the optimal forcing/response pair as a function of Reynolds number and the azimuthal mode. Figure 4b shows that maximum gains and their azimuthal wave numbers increase with the Reynolds number. Figure 4c shows that the same gain and azimuthal wavenumber trends are ob- served when the domain is increased, although at higher nomi- nal values. To explore these trends, we first illustrate the struc- ture of forcing and response modes. Figure 5 shows the leading forcing and response modes for the flow with Re = 60000 and azimuthal number m = 26. The forcing mode is mainly located upstream of the body and mainly generates vortical disturbances aligned in the radial di- rection, while the response is concentrated on the thin bound- ary layer, where the velocity perturbations are aligned with the base flow. The dynamics of the system can be described as follows: radially-aligned vortical disturbances are generated upstream of the body and are convected into the BL close to the stag- nation point.  Once in the boundary layer, these vortical dis- turbances are aligned with the base flow, as, near the stagna- tion point, the velocity profile of the boundary layer is mainly radial, and thus triggers the lift-up mechanism which gen- erates wall-aligned streaks.  The trend of the forcing mode having negligible azimuthal components was observed for all Reynolds numbers at the peak azimuthal mode. The sub-optimal modes are shown in figure 6. Their struc- ture is analogous, with sub-optimal modes showing radial os- cillations in the forcing components, and amplifications 3 to 5 times lower. The basic mechanism for all Reynolds numbers remains approximately the same, with minor variations.  Figure 7 shows the forcing energy for the leading mode for different Reynolds numbers.  If, on one hand, it could be expected that an increase in the Reynolds number would lead to an up- stream extension of the forcing support, as the lower dissipa- tion levels may favour a spreading of the forcing terms, on the other, higher azimuthal wave number reduces the distance between vortical disturbances in the azimuthal direction, in- creasing shear and thus dissipation. These two trends seem to partially compensate for each other, and the forcing domain remains approximately the same.  One notable difference is that for higher Reynolds numbers the forcing mode slightly approaches the r = 0 axis. This can be interpreted as perturba- tions impinging closer to the stagnation to penetrate the thinner BLs.
Compressibility effects
In order to investigate the influence of compressibility on the lift-up mechanism, baseflows for Mach numbers 0.3, 0.6 and 0.9 were obtained using the SU2 code (Economon et al. , 2016). Figure 8 shows the leading gains for these conditions using Chu's and the kinetic energy norms.  In the latter, only forcing associated with the momentum equations and flow ve- locities are used in the norm. Gains obtained using Chu's energy norm are higher than those obtained using only the kinetic energy norm.  For the Mach 0.3 and Mach 0.6 cases, the gains obtained with the ki- netic energy norm are approximately the same as those ob- tained with Chu's norm for a lower Mach number. In partic- ular, this suggests that the main effect of compressibility for Mach numbers up to 0.3 is to couple the perturbations of den- sity and temperature to those of velocities, with negligible in- fluence of the baseflow. The preferred azimuthal number is not affected by compressibility for Mach numbers up to 0.6. For the Mach 0.9, a larger discrepancy is observed, and the optimal disturbances are found at a lower azimuthal num- ber. The difference cannot be attributed to local variations in viscosity, as at the stagnation point it increases only by ≈ 15%, and thus is likely due to more complex phenomena, possibly related to density variation, which is 40% greater at the stag- nation point than on the freestream.
Conclusions
Non-normal mechanisms for perturbation amplification near an axisymmetric leading edge of the blunt body were studied. Using a model paraboloid body, identification of lead- ing forcing and response modes reveals that the lift-up mech- anism is optimally excited upstream of the body, creating vor- tical disturbances that are convected into the boundary layer at the stagnation point. This mechanism is different from those studied by Andersson et al. (2001), where streamwise vortical disturbances are generated directly inside the boundary layer, and Brandt et al. (2004), where it is excited by non-linear in- teractions of the free-stream turbulence, and that of Paredes et al. (2017), where optimal perturbations inside the boundary layer are computed, which could be excited by wall roughness elements. The scaling of the optimal gain and associated azimuthal number found in this study contrasts with those found on de- veloped flat-plate boundary layers. On the latter, optimal gains are seen to scale with Re ∞ Re 2 . , and the spanwise wavelength is proportional to the BL thickness.  On the paraboloid body, the gains are found to scale with Re 1/2 ∞ Re 8 -, thus showing a considerably slower growth with Re . The azimuthal wavenum- ber, which for a given position at the body is inversely propor- tional to the spanwise wavelength, scales with Re 1/4 ∞ δ × 1/2 . The lower gains are explained by the domain truncation: in developed boundary layers the streamwise streaks are con- tinuously forced by vortical disturbances.  Both these struc- tures are damped by viscosity at a rate of Re − 1 , and thus have a total length ∞ Re . For larger Re vortical disturbances excite the streaks via the lift-up mechanisms for a longer time, resulting in the classical Re 2 gain scaling. The finite domain used here reduces the efficacy of the lift-up mechanism, and thus reduc- ing the overall gain. This suggests that leading edge transition, for which the relevant dynamics necessarily occurs in a short distance, will exhibit this slower scaling with respect to Re and/or free-stream turbulent intensity. The relatively low mag- nitude of the gains can be explained by two different mecha- nisms. The first one is the flow curvature, which pushes fast streaks away from the body and slow streaks closer to it, coun- teracting the generated vortical disturbances. This mechanism stabilizes the flow, which should further reduce the gains. The second is the radial evolution of the boundary layer. The con- servation of the azimuthal wavenumber, as opposed to that of the span-wise wavenumber in the flat-plate case, means that a given mode has different spanwise lengths at each body lo- cation. Thus, the length, measured in boundary layer thick- nesses, that maximizes the lift-up effect can only be attained at specific locations. Although this is also the case in spatially evolving boundary layers, the evolution of the BL thickness is usually small, contrasting with the evolution of the paraboloid body radius, which varies very fast near the leading edge. Finally, compressibility effects for subsonic flows were investigated. It is observed that compressibility increases per- turbations amplification and reduces the optimal azimuthal wavenumber when the flow approaches sonic speeds.
Acknowledgements
This work is part of the project TRANSITION sup- ported by Région Nouvelle Aquitaine under grant 2018- 1R10220. Additionally, this work was granted access to the HPC re- sources of IDRIS under the allocations 2020-A0092A10868 and 2021-A0112A10868 made by GENCI.