1. Introduction
The studies of three-dimensional sweep boundary layers date back to years ago. Most of them are based on local models: the attachment-line models and three-dimensional cross-flow models. The sweep flow at the leading edge is often modelled by the sweep Hiemenz configuration, and the most unstable mode is symmetric along the chordwise direction perpendicular to the attachment line (Lin & Malik Reference Lin and Malik1996; Theofilis Reference Theofilis1998; Obrist & Schmid Reference Obrist and Schmid2003). Further downstream, because of the non-alignment in three-dimensional inviscid streamlines and pressure gradients, an inflection point appears in the velocity profile of the three-dimensional boundary layer which leads to cross-flow instability (Reed & Saric Reference Reed and Saric1989; Saric, Reed & White Reference Saric, Reed and White2003). Two types of cross-flow instability have been identified: travelling and stationary modes. The stationary mode plays an important role in the roughness-induced transition while travelling modes are related to external and unsteady perturbations. The dominance of either type of mode depends on the specific configurations and the disturbance environment.
The global stability analyses (GSA) around a swept parabolic body in hypersonic flows were performed by Mack, Schmid & Sesterhenn (Reference Mack, Schmid and Sesterhenn2008) and Mack & Schmid (Reference Mack and Schmid2011). They were able to uncover a global spectrum containing boundary-layer modes, acoustic modes and wave-packet modes. For the first time, these authors not only addressed the attachment-line instabilities, but also showed the connection of attachment-line modes with cross-flow instability through GSA, though Bertolotti (Reference Bertolotti1999) had furnished strong evidence for such connection in local stability analysis. Recently, Meneghello, Schmid & Huerre (Reference Meneghello, Schmid and Huerre2015) performed global stability, receptivity and sensitivity analyses for incompressible flows around a Joukowski airfoil. They found that the global eigenfunction was the most responsive to forces applied in the vicinity of the attachment line. However, the analyses of such flows in the hypersonic region are quite limited.
The receptivity process describes the procedures of penetration of external perturbations into the boundary layer and excitation of modes inside the boundary layer, which plays an important role in the boundary-layer transition. Based on theoretical methods, such as finite-Reynolds-number methods (Choudhari Reference Choudhari1994) and triple-deck theory (Ruban Reference Ruban1984), the bi-orthogonal eigenfunction system was found to be an effective tool for the local receptivity analyses (Hill Reference Hill1995; Fedorov & Khokhlov Reference Fedorov and Khokhlov2002; Tumin Reference Tumin2007). Recently, a comprehensive review of this method is given by Tumin (Reference Tumin2020). However, there exists no application of this approach to complex hypersonic flows.
There are two objectives in the present work. Firstly, it is to understand the stability features of hypersonic attachment-line instability over a large range of sweep angles. Secondly, it is to highlight the receptivity behaviours in the vicinity of the leading edge of a swept blunt body in hypersonic flows with the multi-dimensional bi-orthogonal eigenfunction system. This paper is organized as follows. In § 2, the governing equations are introduced and the bi-orthogonal eigenfunction system for the global stability system is established for solving the receptivity problem. The results for global stability/receptivity analysis are presented in § 3 and the conclusion is given in § 4.
2. Theoretical framework
The stability analyses of hypersonic steady flows is performed around the front part of a blunt body, as sketched in figure 1. In general, the compressible Navier–Stokes equations for variable $\boldsymbol {Q} = (\rho ,u,v,w,T)$, denoting the density, velocity components and temperature, can be written as
where $\boldsymbol{\varGamma} $ is the coefficient matrix between conservative and primitive variables. The linearized Navier–Stokes equations (LNSE) around a stationary state $\boldsymbol {Q}$, with $\boldsymbol {\mathscr {R}}(\boldsymbol {Q}) = 0$, can be represented by the combination of a linearized operator $\boldsymbol {\mathscr {L}} = \partial {\boldsymbol {\mathscr {R}}}/\partial \boldsymbol {Q}$ and a perturbation $\boldsymbol {p}=(\hat {\rho },\hat {u},\hat {v},\hat {w},\hat {T})$ field. This process forms a homogeneous system as
The solution of the above problem (2.2) is obtained with the method of Laplace transform.
The Laplace transform $\hat {\boldsymbol {p}}(\sigma )$ of function $\boldsymbol {p}(t)$ in (2.2) is calculated from the following integral
with $\sigma$ being a complex parameter and $\sigma = -\textrm {i} \omega$. From $\hat {\boldsymbol {p}}(\sigma )$, the original function $\boldsymbol {p}(t)$ can be restored with the Bromwich integral:
The inversion integral (2.4) suggests that the solution of LNSE (2.2) can be represented as the superposition of normal modes, which can be written as
where c.c. is the complex conjugate. Substitution of (2.5) into (2.2) leads to a standard global stability problem (Theofilis Reference Theofilis2011):
The detailed expressions of these operators were given in the authors’ previous work (Xi, Ren & Fu Reference Xi, Ren and Fu2021).
On the wall surface and the bow shock, all perturbations except for the density are set to zero ($(\hat {u},~\hat {v},~\hat {w},~\hat {T}) = 0$). Along the chordwise direction, at the exits, high-order extrapolations are applied for all perturbations. Because the operators are non-self-adjoint and the eigenfunctions are not orthogonal, a complete solution of LNSE needs the help of adjoint eigenfunctions. An inner-product for two arbitrary vectors $\boldsymbol {a}$ and $\boldsymbol {b}$ is defined as $(\boldsymbol {a},\boldsymbol {b}) = \iint _{V}\boldsymbol {a}\boldsymbol {\cdot }\boldsymbol {b}^{\textrm {T}}\,\textrm {d}x\,\textrm {d}y$, where the superscript T stands for the transpose and $V$ represents the integration region over $x\text{--}y$ plane. Based on the inner-product, the adjoint problem of (2.6) can be expressed as
where the superscript ${\dagger}$ stands for the adjoint variables and operators. It should be noted that boundary terms can be eliminated with proper boundary conditions in the integration by parts. In addition, based on the definition of inner-product, the adjoint and the original problems have the same eigenvalue spectrum.
A bi-orthogonal eigenfunction system $\{\hat {\boldsymbol {p}}_{a},\hat {\boldsymbol {p}}^{{\dagger} }_{b}\}$ is now formed with the help of eigenfunctions from the original problem $\hat {\boldsymbol {p}}_a$ and the adjoint problem $\hat {\boldsymbol {p}}^{{\dagger} }_b$. The direct and adjoint equations lead to the following expression:
where $\omega _a$ stands for the eigenvalue of the original problem and $\omega_b^{{\dagger}}$ stands for the eigenvalue of the adjoint problem. Based on the definition of the adjoint problem, then
With the above equations (2.8) and (2.9), the following relationship can be achieved:
If $\omega ^{{\dagger} }_b = \omega _a$ and one can further define that
where $C_0^{a}$ represents a normalization constant for a specific mode $a$. If $\omega ^{{\dagger} }_b \ne \omega _a$, the eigenfunctions are orthogonal to each other, so
Relations (2.11) and (2.12) form a bi-orthogonality condition in a two/three-dimensional domain. An inverse Laplace transform (2.4) is adopted for the return to physical space. Also, generally, the solution of linear Navier–Stokes equations can be expanded in a combination of continuous and discrete modes as
where $m$ and $n$ represent the indices of discrete modes and continuous branches, respectively; $k$ stands for the integration parameter along continuous branches, $\boldsymbol {C}_m^{d}$ and $\boldsymbol {C}_{n}^{c}$ represent the amplitude of those discrete modes and branches, respectively. With relationships (2.6)–(2.13), for a specific discrete mode $m$, then
A generic external force $\boldsymbol {f}$, representing the source term in the flow field, together with the external disturbances $\boldsymbol {g}$ standing for any surface/free-stream perturbations, can be added to system (2.2). The amplitudes of these additional perturbations are assumed to be small hence the whole system can still be described by the linear approximation. Thus, the inhomogeneous system can be expressed as
Similar to the previous process, a normal-modes assumption is adopted with respect to the force terms and boundary conditions:
The relative system in phase space can also be written as
with the inhomogeneous boundary conditions $\hat {\boldsymbol {p}}=\hat {\boldsymbol {g}}$ on boundary lines.
Considering a dot product between an adjoint eigenvector $\hat {\boldsymbol {p}}^{{\dagger} }$ and (2.17) for any discrete mode $m$, with integration over the whole domain, we have
Based on (2.7), the following adjoint homogeneous equation can be obtained:
Substituting (2.19) into (2.18), the following relationship is then derived for a specific discrete mode $m$:
Here $B.C.$ represents concomitant boundary terms due to the inhomogeneous boundary conditions determined by using Green formulas (detailed expression is given in Appendix A). With (2.4), (2.14) and (2.20), the following integral equation is derived for the physical perturbation $\boldsymbol {p}_m$:
By closing the Bromwich integral (2.21) in the complex $\omega$-plane, the integral is obtained as the residue value at the point $\omega = \omega _m$:
The amplitude $\boldsymbol {C}_m^{d}$ for a specific discrete mode $m$ can be recovered as
For the fact that $C_0^m$ is a constant for a specific discrete mode, the adjoint field represents a scaled Green's function for the receptivity problem. In fact, formula (2.23) is consistent with those obtained from incompressible flows (Giannetti & Luchini Reference Giannetti and Luchini2007). It can thus be considered as a general form for both incompressible and compressible flow systems.
3. Global stability and receptivity of the leading-edge region
The test model is a two-dimensional parabolic body swept with a sweep angle $\varLambda$. The model geometry and computational domain are shown in figure 1. The free-stream Reynolds number $Re_{\infty }$, the sweep Reynolds number $Re_s$, the free-stream Mach number $M_{\infty }$, the sweep Mach number $M_s$ and the recovery temperature $T_r$ are defined as
Here, $\xi _w$ is a constant for specific free-stream conditions ($M_{\infty }$ and $\varLambda$) and determined based on the study of Reshotko & Beckwith (Reference Reshotko and Beckwith1958); $R=0.1\ \textrm {m}$ represents the radius of the leading edge; ${\boldsymbol {V}}_{\infty }$ stands for the free-stream velocity vectors with $U_{\infty },V_{\infty }$ and $W_{\infty }$ along $x,y$ and $z$ direction, respectively. Here $T_{\infty }$ of $50.93\ \textrm {K}$ and $T_{0}$ stand for the free-stream and stagnation temperature, respectively. The parameters $c_{\infty }$ and $c_s$ are the speeds of sound before and after the leading shock, $\nu _r$ represents the kinematic viscosity at $T_r$. The viscosity length scale $\delta ^*$ is defined as $\delta ^* = \sqrt {{\nu _r R}/{2 U_2}}$, where $U_2$ represents chordwise velocity behind the shock. The Prandtl number $Pr$ of $0.71$ and the specific heat ratio $\gamma$ of $1.4$ are set following the ideal gas assumption of air. Table 1 lists the parameters of the chosen cases. Here, $M_s$ is used to define the supersonic and hypersonic configurations as mentioned below, and $T_w = T_r$ is the surface temperature, which is used for all cases. The spanwise wavenumber $\beta = 2{\rm \pi} \delta ^*/ \lambda ^*_z$ is adopted corresponding to the physical wavelength $\lambda ^*_z$ of 11.7 mm.
Firstly, a high-order shock fitting method (Zhong Reference Zhong1998) is employed to obtain the basic flow fields over the largest domain $R3$ as shown in figure 1. In the shock fitting method, the shock is modelled as a boundary of the computational domain; smooth fields can be achieved for the adjoint calculation without considering shock discontinuities (Giles & Pierce Reference Giles and Pierce2001). For all cases, 1201 grid points are used along the wall surface and 201 grid points in the wall-normal direction (at least 51 points are clustered inside the boundary layer). Compared with the previous DNS study (Mack et al. Reference Mack, Schmid and Sesterhenn2008) with $255\times 128$ grids, the basic flows are adequately resolved in the present grid resolution. A matrix-based high-order global stability analysis is then performed to solve the direct and adjoint problems. A Krylov–Shur method (Stewart Reference Stewart2002a,Reference Stewartb), based on PETSc (http://www.mcs.anl.gov/petsc) and SLEPc (http://slepc.upv.es) with shift-and-invert spectral transformation is adopted to recover a eigenvalue window (20–50) of interest. Sparse linear algebra packages, MUMPS (http://mumps.enseeiht.fr) and SuperLU (https://github.com/xiaoyeli/superlu_dist) are applied to undertake the inverse of the matrix during the spectral transformations. A complete review of multi-dimensional and global stability analysis techniques is given in Theofilis (Reference Theofilis2011). More details of the validations/verifications for the present solvers can be found in Xi et al. (Reference Xi, Ren and Fu2021). A standard sixth-order finite difference scheme is used to discretize both directions. Equally spaced grid points are adopted along the wall surface and the following equation is used to cluster grids at the wall:
Here, $\eta$ is the uniform grid in computational domain $[-1,1]$, $H_s$ stands for the distance between the shock and the wall, $h$ represents the wall-normal grid in the physical domain and $\eta _h$, set to 0.025, is a location parameter for clustering control. To achieve a grid-independent solution, the number of discretization points is checked for convergence. A large number of grid points up to $5401\times 401$ are employed to identify the structure of direct and adjoint eigenfunctions with a minimum of 20 points per wavelength in each direction for different sizes of computational domains.
The spectrum of the attachment-line modes for CASE P2 is shown in figure 2(a). The branches of attachment-line modes are marked by black filled dots. Similar to the incompressible cases (Lin & Malik Reference Lin and Malik1996; Meneghello et al. Reference Meneghello, Schmid and Huerre2015), symmetric ($S1,S2,\ldots$) and antisymmetric ($A1,\ldots$) modes alternate from the most unstable to the most stable. The leading symmetric mode $S1$ has the highest growth rate. All the attachment-line modes are travelling in the spanwise direction at a nearly constant phase speed around $0.5$. For the high-sweep Mach number case (CASE P8), however, as shown in figure 2(b), only one symmetric discrete mode is obtained that is marked with $S$. Based on the recent study (Xi et al. Reference Xi, Ren and Fu2021), for large sweep Mach numbers, the attachment-line mode is inviscid in nature, while for lower sweep Mach numbers, the attachment-line instability exhibits the behaviours of viscous Tollmien–Schlichting waves.
The leading attachment-line modes are calculated through three different domains (see figure 1). For all domains, the same leading eigenvalues and consistent eigenfunctions are obtained as shown in figure 3. It indicates that the characteristics of attachment-line modes are not much affected by the computation domain size. From figure 3(a), the behaviour of the surface-density perturbations, which are normalized with the values at the attachment line, are of great difference in the downstream region $s/R>1$ for these two cases. Based on eigenfunctions in figure 3(b,c), it is seen that the downstream perturbations in CASE P2 are of the cross-flow type while those in CASE P8 are of the second Mack-mode type, as highlighted more clearly in figure 4(e) (CASE P2) and figure 4(f) (CASE P8).
Figure 4(a,b) display the structures of direct and adjoint eigenvectors obtained from the largest R3 domain for the leading eigenvalues of CASES P2 and P8, respectively. The structures are visualized with the iso-surfaces of the real part of the direct spanwise velocity $\hat {w}$ and the adjoint spanwise velocity $\hat {w}^{{\dagger} }$. For both cases, red/blue colours stand for positive/negative values $(\pm 10^{-6}$) for direct eigenfunctions and blue/white colours represent positive/negative values ($\pm 10^{-3}$) for adjoint eigenfunctions. Figure 4(a,b) display the typical features for the direct global eigenfunctions exhibitting attachment-line instability near the leading edge, while either cross-flow instability (CASE P2) or second Mack mode (CASE P8) dominate downstream flows. For smaller $M_s$ (CASE P2), from the attachment-line plane, marked $a1$ in figure 4(a) and shown in figure 4(c), to the downstream plane, marked $a2$ in figure 4(a) and shown in figure 4(e), perturbations evolve away from wall surfaces forming cross-flow vortices aligned with external streamlines, as also reported by Mack et al. (Reference Mack, Schmid and Sesterhenn2008) and Meneghello et al. (Reference Meneghello, Schmid and Huerre2015). For larger $M_s$ (CASE P8), despite the similar mode structure at the attachment-line plane, marked $b1$ in figure 4(b) and shown in figure 4(d), the perturbations at the downstream plane, marked $b2$ in figure 4(b) and shown in figure 4(f), exhibit the behaviour of the second Mack mode instability with perturbations mainly located below the sonic line. In contrast to the feature that the perturbations of direct mode cover a large area, the adjoint mode appears only in the vicinity of the attachment line as shown in figure 4(a,b,g,h). Figure 4(g,h) show the contours of adjoint spanwise velocity $\hat {w}^{{\dagger} }$ for CASE P2 and P8 at the attachment-line planes, respectively. The basic structures of adjoint fields for both cases are similar and the bulk of adjoint fields moves away from the wall surface for large sweep Mach numbers. Moreover, based on the theoretical receptivity analysis on external force $\boldsymbol {f}$, this region is the most responsive to this discrete global mode. In other words, even though the direct global mode covers a large region further away from the attachment line, it is still possible to control/excite the mode by introducing forces at a relatively small region close to the attachment line. This finding is consistent with those from incompressible flows by Meneghello et al. (Reference Meneghello, Schmid and Huerre2015), indicating that the leading-edge control might be effective to the hypersonic flows over swept blunt bodies.
The characteristics of the leading attachment-line mode on sweep angles and the relative sweep Mach numbers are shown in figure 5. As the sweep angle varies from $20^{\circ }$ to $70^{\circ }$, the growth rate of the leading attachment-line mode shows a rising and then declining trend, reaching the locally highest point around $40^{\circ }$. When the sweep angle further increases, the growth rate of the leading attachment-line mode exhibits the features of monotonous increase and becomes positive again at large sweep Mach numbers. Moreover, this finding is consistent with two facts: the first is that the plane stagnation flow is known to be linearly stable for three-dimensional perturbations, as can be seen as the limit of zero sweep angle for the present cases; the other is that the attachment-line instability is found to become unstable for large sweep angles (Gaillard, Benard & Alziary de Roquefort Reference Gaillard, Benard and Alziary de Roquefort1999).
Based on (2.23), the amplitudes of the modes excited by various boundary perturbations are mainly determined by the concomitant boundary term ($B.C.$). From the expression in Appendix A, it is clearly shown that the receptivity is evaluated from the gradients of adjoint variables and physical perturbations on the boundary. Moreover, this procedure permits the extraction of the receptivity amplitude pertaining to any discrete modes for any type of boundary perturbations at any location. Taking the surface perturbations in CASES P2 and P8, for instance, two types of surface perturbation (detailed expressions can be found in Appendix B) are calculated, as shown in figure 6. It demonstrates that the strong receptive regions of the leading attachment-line modes to surface perturbations are in the vicinity of attachment line and the surface perturbations at the exact attachment line are the most effective. Also, the surface vibrations are found to be a little more effective than suction/blowing, because of the higher excited amplitudes. Moreover, as the delta-function form perturbations are used and the basis for representation of LNSE cover the whole domain, the distributed surface perturbations receptivity analyses can be easily performed by integrations along the finite length. Further, as the present framework is from the basis of LNSE, the receptivity problem for the excitation of discrete modes of free-stream disturbances with surface inhomogeneities can be also taken into account.
As mentioned by Tumin (Reference Tumin2020), the eigenfunction expansion methods for the two- and three-dimensional problems are the natural extensions of the widely used local cases. The analysis of discrete modes is similar for both local and global cases and can be easily performed. However, problems become more complicated with the continuous spectra since the mathematical properties of ordinary and partial differential operators are very different. Although uncertainties associated with the continuous spectra for global cases still exist, the bi-orthogonal eigenfunction system is still a powerful tool for solving global receptivity problems of discrete modes for complex high-speed boundary layers.
4. Conclusions
Global instabilities and receptivities around the leading edge of a swept blunt body are studied in the hypersonic-flow region. The characteristics of the leading attachment-line mode under the variation of sweep angles from $20^{\circ }$ to $70^{\circ }$ are obtained for the first time. As the sweep angle increases, the growth rates of the leading attachment-line modes exhibit a rising then declining trend and then rising again in which the modes also show the transformation from the features of the cross-flow instability to the second Mack-mode instability further downstream. Moreover, a general bi-orthogonal eigenfunction system for the hypersonic global stability system is established to address receptivity problems to any external forces and boundary perturbations. The receptivity analyses indicate that the leading attachment-line mode is the most responsive to external disturbances in the vicinity of the leading edge, in the hypersonic region. Furthermore, it illustrates that the present framework can be easily applied to both local and non-local linear stability/receptivity analyses for real geometry.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2021.217.
Acknowledgements
We appreciate Dr M. Choudhari and Professor A. Hanifi for helpful discussion on the receptivity analysis.
Funding
This work received support from the National Key Project GJXM92579, the National Sci. & Tech. Major Project (2017-II-0004-0016), NSFC Grant 12072173 and the National Key Research and Development Plan through project No. 2019YFA0405201.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Concomitant boundary terms
The concomitant boundary terms are determined by using Green formulas. Details for derivation of these terms can be found in supplementary material available at https://doi.org/10.1017/jfm.2021.217. Taking the detailed explicit forms of the operators $\boldsymbol {\mathscr {L}}$, the $B.C.$ term can be expressed as
where the superscript ${\dagger}$ and subscript $bc$ represent the adjoint variables and variables at the boundary line $\varGamma$. The direction of the boundary line $\varGamma$ is defined to ensure the computation domain is always on the left.
Appendix B. Boundary conditions for surface perturbations
The boundary conditions for perturbations can be expressed as
where $H(x,y)$ is the shape function of perturbations and is chosen as the Dirac delta function
in the present study. Here, $S_s$ stands for the surface distance away from the attachment line, and $S_l$ represents a specific location along the surface, where the perturbations are embedded. When a surface is subject to small vibrations, the boundary conditions with the frequency $\omega _r$ can be specified as
where $\boldsymbol {n} = \left (\boldsymbol {n}_x,\boldsymbol {n}_y\right )$ stands for the surface normal direction. In the circumstance of wall-normal blowing/suction, the boundary conditions are given by