Global stability analysis of axisymmetric boundary layer over a circular cone

This paper presents the linear Global stability analysis of the incompressible axisymmetric boundary layer on a circular cone. The base flow is considered parallel to the axis of cone at the inlet. The angle of attack is zero and hence the base flow is axisymmetric. The favorable pressure gradient develops in the stream-wise direction due to cone angle. The Reynolds number is calculated based on the cone radius (a) at the inlet and free-stream velocity ($U_{\infty}$). The base flow velocity profile is fully non-parallel and non-similar. Linearized Navier-Stokes equations (LNS) are derived for the disturbance flow quantities in the spherical coordinates. The LNS are discretized using Chebyshev spectral collocation method. The discretized LNS along with the homogeneous boundary conditions forms a general eigenvalues problem. Arnoldi's iterative algorithm is used for the numerical solution of the general eigenvalues problem. The Global temporal modes are computed for the range of Reynolds number from 174 to 1046, semi-cone angles $2^o$, $4^o$, $6^o$ and azimuthal wave numbers from 0 to 5. It is found that the Global modes are more stable at higher semi-cone angle $\alpha$, due to the development of favorable pressure gradient. The effect of transverse curvature is reduced at higher semi-cone angles ($\alpha$). The spatial structure of the eigenmodes show that the flow is convectively unstable. The spatial growth rate ($A_x$) increases with the increase in semi-cone angle ($\alpha$) from $2^o$ to $6^o$. Thus, the effect of increase in semi-cone angle ($\alpha$) is to reduces the temporal growth rate ($\omega_i$) and to increases the spatial growth rate ($A_x$) of the Global modes at a given Reynolds number.


I. INTRODUCTION
The laminar-turbulent transition in boundary layers has been a subject of interest to many researchers in past few decades. It is important to understand the onset of transition in boundary layers as the flow pattern and its effects are very different in laminar and turbulent flows. For example, the drag force in a turbulent flows is much higher than that of a laminar flows. At low free-stream levels the boundary layers undergo transition through the classical TS wave mechanism. The amplification of the disturbance waves is the primary step in the transition process and this is studied in linear stability analysis. The results from stability analysis and the prediction of transition onset is very useful in hydrodynamics and aerodynamics applications like submarines, torpedoes, rockets, missiles etc.
The linear stability analysis of shear flows with parallel flow assumption is well understood by the solution of the Orr-Sommerfeld equation [1]. It is known that the stability characteristics in a boundary layer is strongly influenced by various factors such as pressure gradient, surface curvature free-stream turbulence level.
The boundary layer forms over a right circular cone is axisymmetric and it is qualitatively different from a flat-plate boundary layer. The boundary layer on a circular cone develops continuously in spatial directions, hence the parallel flow assumption is not valid. Due to the wall normal velocity component, the boundary layer is non-parallel. The instability analysis of such twodimensional base flow is called Global stability analysis [2]. The literature on the instability analysis of axisymmetric boundary layer is very sparse. The available literature on the axisymmetric boundary layer is limited to local stability analysis. Rao [3] first studied the stability of axisymmetric boundary layer. He found that non-axisymmetric disturbances are less stable than that of two-dimensional disturbances. Tutty [4] investigated that for non-axisymmetric modes critical Reynolds number increases with the azimuthal wave number. The critical Reynolds number found to be 1060 for N = 1 mode and 12439 for N = 0 mode. The axisymmetric mode is found least stable fourth mode. Vinod [5] investigated that higher non-axisymmetric mode N = 2 is linearly stable for a small range of curvature only. The helical mode N=1 is unstable over a significant length of the cylinder, but never unstable for curvature above 1. Transverse curvature has overall stabilizing effect over mean flow and perturbations. Malik [6] studied the effect of transverse curvature on the stability of incompressible boundary layer. They investigated that the body curvature and streamline curvature are having significant damping effects on disturbances. The secondary instability of an incompressible axisymmetric boundary layer is also studied by Vinod [7]. They found that laminar flow is always stable at high transverse curvature to secondary disturbances. The Global stability analysis for hypersonic and supersonic flow over a circular cone is reported in the literature. In supersonic and hypersonic boundary layers there exists higher acoustic instability modes with higher frequencies in addition to first instability mode [8]. It has been confirmed by the experiments that the two dimensional Mack mode dominates in hypersonic boundary layer [9][10][11]. Experimental results show that leading bluntness affects the transition location significantly in circular cone [12]. In case of axisymmetric boundary layer on a circular cylinder, transverse curvature effect increases in the stream-wise direction, which helps in stabilizing the flow. In case of a axisymmetric boundary layer on a circular cone, the body radius of the cone increases at higher rate than a boundary layer thickness (δ) and hence the transverse curvature effect (δ/a) reduces in the stream-wise direction as shown in figure  4. However, the favorable pressure gradient develops in the stream-wise direction due to cone angle (α). Thus, it becomes interesting to study the combined effect of transverse curvature and pressure gradient on the stability characteristics. The two dimensional global modes are also computed for the flat-plate boundary layer by some investigators [13][14][15]. However, this is the first attempt to carry out Global stability analysis of incompressible boundary layer over a circular cone. The main aim of this paper is to study the Global stability characteristics of the axisymmetric boundary layer on a circular cone and the effect of the transverse curvature and pressure gradient on the stability characteristics.

II. PROBLEM FORMULATION
The standard procedure is followed for the derivation of the Linearized Navier-Stokes equations(LNS) for the disturbance flow quantities. The Navier-Stokes (N-S) equations for the base flow and instantaneous flow are written in the spherical coordinate system (r,θ,φ). The equations are non-dimensionalised using free-stream velocity U ∞ and body radius of the cone a at the inlet. The LNS for disturbance flow quantities are obtained by subtraction and subsequent linearization. The base flow is two dimensional and disturbances are three dimensional in nature. This will determine whether the small amplitudes of the disturbances to grow or decay for a given steady laminar flow. The Reynolds number is defined as, where a is the surface radius of cone at inlet and ν is the kinematic viscosity.
The flow quantities are presented as sum of the base flow and the perturb quantities as, U r = U r +u r , U θ = U θ +u θ , U φ = 0+u φ , P = P +p, (2) The disturbances are considered in normal mode form and varying in radial(r) and polar(θ) directions. Thus the disturbances having periodic nature in the azimuthal (φ) direction. where, where,

A. Boundary conditions
No slip and no penetration boundary conditions are considered at the surface of the cone. The magnitude of all the disturbance velocity components are zero at the solid surface of the cone due to viscous effect. u r (r, θ min ) = 0, u θ (r, θ min ) = 0, u φ (r, θ min ) = 0 (9) It is expected to vanish the disturbance flow quantities at free-stream far away from the solid surface of the cone. The Homogeneous Dirichlet conditions are applied to all velocity and pressure disturbances at free-stream.
The boundary conditions are not straight forward in the stream-wise direction for Global stability analysis.
As suggested by Theofilis [16], homogeneous Dirichlet boundary conditions are considered for velocity disturbances at inlet. Here we are interested in the disturbances evolved within the basic flow field only. The boundary conditions at outlet can be applied based on the incoming and outgoing wave information [17]. Such conditions impose wave like nature of the disturbances and so it is more restrictive in nature which is not appropriate from the physical point of view. Even stream-wise wave number α is not known initially in case of the Global stability analysis. Alternatively one can impose such numerical boundary condition which extrapolate the information from the interior of the computational domain. Linear extrapolated conditions are applied by several investigators in such case. The review of the literature on Global stability analysis suggests that the linear extrapolated boundary conditions are least restrictive in nature [16,18]. Thus, we considered linear extrapolated conditions at the outlet for the numerical solution of the general eigenvalues problem.
u rn−2 [r n −r n−1 ]−u rn−1 [r n −r n−2 ]+u rn [r n−1 −r n−2 ] = 0 (12) similarly, one can write extrapolated boundary conditions for polar and azimuthal disturbance components u θ and u φ respectively. The boundary conditions for pressure do not exist physically at the wall. However compatibility conditions derived from the LNS equations are collocated at the solid wall of the cone [16].
The Linearized Navier-Stokes equations are discretized using Chebyshev spectral collocation method. The Chebyshev polynomial generates non uniform grids. By nature it takes more number of points towards the ends. This is favorable arrangement for the boundary value problems.
The gradient of the velocity disturbances are very high near the solid surface of the cone. To improve the spatial resolution in the boundary layer region near the wall, it is necessary to stretch more number of collocation points near the wall. The following algebraic equation is used for grid stretching ( [19]).
In the above grid stretching method half number of the collocation points are concentrated within the θ i angle from the lower surface θ min only. The non uniform nature of the distribution for the collocation points in the stream-wise direction is undesirable. The maximum and minimum distance between the grid points are at the center and at the end respectively in stream-wise direction.
Thus it makes poor resolution at the center of the domain and very small distance between the grids at the end gives rise to Gibbs oscillations. To improve the spatial resolution and to minimize the Gibbs oscillation in the solution, grid mapping is implemented in stream-wise direction using following algebraic equation [20].
The value of α m is chosen carefully to improve the spatial resolution in the stream-wise direction. Very small value of α m keeps the grid distribution almost like a Chebyshev and near to unity almost uniform grid distribution. For detail description of grid mapping readers are suggested to refer [20]. To incorporate the effect of physical dimensions of the domain [L r , L θ ] along with grid stretching and mapping it is required to multiply the Chebyshev differentiation matrices by proper Jacobean matrix. Once all the partial derivatives of the LNS are discretized by spectral collocation method using Chebyshev polynomials, the operator of the differential equations generate matrices A and B. These matrices are square, real and sparse in nature, formulates general eigenvalues problem.
where A and B are square matrix of size 4nm, iω is an eigenvalues and φ is a vector of unknown disturbance amplitude functions u r ,u θ , u φ and p. The above mentioned all the boundary conditions are properly incorporated in the matrix A and B.

B. Solution of general eigenvalues problem
The A and B matrix are very large in size and sparse in nature. The full solution of the general eigenvalue problem is very difficult due to large size of matrix A and B. However, the shear flow becomes unstable due to very few leading eigenmodes only. The eigenmodes with largest imaginary part are important for instability analysis. Hence we are interested in the few eigenvalues and associated eigen vectors only. The QZ algorithm is not the right choice for the solution, because it computes full spectrum. Arnoldi's iterative algorithm is selected to compute the few selected eigenvalues. The Krylov subspace provides possibility of extracting major part of the spectrum using shift and invert strategy. The computations of Krylov subspace along with Arnoldi's algorithm applied to eigenvalues problem becomes easy.
where, λ being the shift parameter and µ is the eigenvalues of the converted problem. Sometimes it is also called spectral transformation, which converts generalized eigenvalues problem to standard eigenvalues problem. The Krylov subspace may be computed by successive resolution of linear system with matrix (A − λB),using LU decomposition. Full spectrum method is employed for this small subspace to get good approximate solution of the original general eigenvalues problem [21]. The major part of the spectrum can be extracted with large size of Krylov subspace. Generally the computed spectrum is always nearby the shift parameter. The good approximation of shift parameter reduces the number of iterations to converge the solution to required accuracy level. However larger size of subspace makes the solution almost independent from the shift parameter λ. We tested the code for several values of shift parameter. The convergence of the solution depends on the value of shift parameter. Good approximation of the shift value needs less number of iterations. However we have taken maximum number of iteration equal to 300, hence convergence of the solution may not be affected by shift parameter. The grid convergence study for the base flow is obtained, using U (x=0.5, r=0.038965 ) and V ( x=0.5, r=0.038965) for U∞=0.1 m/s. The grid refinement ratio (α)in each direction is 1.4142. The relative error ( ) and Grid Convergence Index (GCI) are calculated using two consecutive grid size. The j and j+1 represents course and fine grids respectively. = Given the large subspace size of k = 250,the part of spectrum for our instability analysis could be recovered in the one computation which took about 4 hours on Intel Xenon(R)CPU E5 26500@2.00GHz × 18.

III. BASE FLOW SOLUTION
The base flow velocity profile is obtained by numerical solution of the steady incompressible axisymmetric Navier-Stokes equations using finite-volume code ANSYS FLUENT. The incoming flow is parallel to the axis of the cone at inlet. The axisymmetric domain is modeled with 1m and 0.5m in axial and radial directions respectively as shown in figure 1.
Appropriate boundary conditions are applied to close the formulation of the above problem. The uniform inlet velocity U ∞ is imposed at the inflow boundary. No-slip and no penetration conditions are applied on the surface of the cone. At far away from the solid surface, slip boundary condition is applied. The outflow boundary conditions are applied at outlet, which consists ∂U ∂x = 0, using cubic spline interpolation. The simulation was started initially with the coarse grid size of 250 and 125 grid points in axial and radial directions respectively. The grid was refined with a factor of 1.4142 in each directions. The discretization error was calculated through Grid Convergence Index (GCI) study on four consecutive refined grids [22]. The monotonic convergence is obtained for all this grids as expected. The error and GCI are computed for two different field values U (x=0.5, r=0.038965) and V (x=0.5, r=0.038965) near the solid boundary where the gradient is higher. The GCI and error were computed for four different grids as shown in table I. The error between Mesh 1 and Mesh 2 is too small. The GCI has also reduced with the refinement of the grids. Thus,the solution has converged monotonically towards the grid independent one. The further refinement in the grids will hardly improve the accuracy  [23]. The governing equations for instability analysis are derived in spherical coordinates. The base velocity field is transformed to spherical coordinates to perform Global stability analysis.

IV. CODE VALIDATION
To validate the Global stability results for the incompressible flow over a circular cone, a blunt cone with the very small semi-cone angle, α = 10 −12 is considered. Thus the geometry of the circular cylinder and a blunt is nearly similar due to very small semi-cone angle α as shown in figure 5. A rectangle epqr shows computational domain for the cylindrical coordinates and efgh for spherical coordinates. L x and L r are the stream-wise domain length for cylindrical and spherical coordinates. The stability equations for the axisymmetric flow over circular cylinder and cone are written in the polar cylindrical and spherical coordinates respectively. The Reynolds number is computed based on the body radius (a) at inflow and free-stream velocity for both the case. The stream-wise and wall normal extent of the domain is also same for both the domain. The homogeneous Dirichlet and linear extrapolated boundary conditions are imposed at inflow and outflow for both the case. No-slip and free-stream boundary conditions are applied as usual in wall normal direction. The Global stability results are already validated for the axisymmetric boundary layer over a circular cylinder [24]. Figure 6 shows the comparison of the eigenspectrum for the Global stability analysis of the axisymmetric boundary layer over a circular cone and cylinder for N = 0 and Re = 1000. The discrete and continuous part of the spectrum is in excellent agreement with each other. Figure 7 shows the comparison of the eigen-functions for the least stable eigenmode at the same stream-wise location. The eigen-functions are also in good agreement with each other.

V. RESULTS AND DISCUSSIONS
In the present analysis Reynolds number is varying from 174 to 1047 and azimuthal wave-numbers from 0 to 5 for different semi-cone angles α. The Reynolds number is computed based on the cone radius(a) at inlet and free-stream velocity (U ∞ ). The stream-wise length is normalized by cone radius(a). The domain size in polar(wall normal) direction is taken as L θ = 12 o . The stream-wise(radial) length is taken 0.75m. The number of collocation points taken in radial and polar directions are n=121 and m=121 respectively. The general eigenvalues problem is solved using Arnoldi's iterative algorithm. The computed eigenvalues are accurate upto three decimal point. The additional lower resolution cases were also run to verify the monotonic convergence of the solution. Arnoldi's iterative algorithm is used to solve the general eigenvalues problem. Heavy sponging is applied at the outflow boundary to prevent spurious reflection. The selected Global eigenmodes are also checked for the spurious mode. A grid convergence study was performed to check the accuracy level of the solution and appropriate grid size. Table II shows the values of two leading eigenvalues computed for Re = 349, N = 1 and semicone angle α = 2 o using three different grid size. The grid resolution was successively improved by a factor of 1.14 in radial and polar direction respectively. The real and imaginary parts of the eigenvalues shows monotonic convergence of the solution with the increased resolution.
In the table II n and m indicates number of collocation points in the radial and polar directions respectively. The relative error is computed between two consecutive grid sizes for real and imaginary parts. The largest associated error among both the eigenmode is considered. The mesh 1 is used for all the computations for stability analysis results reported here.
A. Semi-cone angle α = 2 o Figure 8 shows the eigenspectrum of axisymmetric mode (N = 0) for Re = 349 and semi-cone angle α = 2 o . The eigenmodes marked by square and circle are called stationary and oscillatory mode. The stationary mode has a complex frequency ω = 0 − 0.02667i. The Global mode is temporally stable because ω i < 0 and hence the amplitudes of the disturbances decay in time. Figure 9(a) and (b) presents the two dimensional spatial structure of the eigenmodes for radial (u r ) and polar (u θ ) velocity disturbances. The magnitudes of the disturbance amplitudes are zero at inlet as it is applied inlet boundary condition. The disturbance amplitudes evolve with the time in the flow domain and moves in the stream-wise direction towards downstream. The size and magnitudes of the disturbances grows as they move towards downstream. The magnitudes of the u r disturbances are one order higher then that of u θ . The polar disturbances u θ contaminates the flow field upto large extent than that of u r disturbances, however its magnitudes very small compare to u θ . Figure 10 shows the spatial structure of oscillatory eigenmode. The associated complex frequency for this mode is ω = 0.008154−0.03937i. This oscillatory Global mode is also temporally stable, because ω i < 0. The variation of the amplitude is not monotonic for this eigenmode. The wave-like nature of the disturbances are found for this mode. The oscillatory modes evolved in the flow field, grows in size and magnitude while moving towards the downstream and hence the flow is convec- tively unstable. Figure 11 shows the two dimensional mode structure of the u r and u θ for the ω r = 0.1775. The stream-wise domain length is 0.0000 . The disturbances are observed to evolve near the wall surface and decay exponentially at free-stream. The typical length scale of the wavelet structure decreases with the increase in the frequency ω r . It also demonstrates that the region of contamination reduces with the increases in frequency. The length scale of the wavelet structure is found small for the increased frequency as shown in figure 11. B. Semi-cone angle α = 4 o Figure 12 shows the spectrum for helical mode N=1, Re=523 and semi-cone angle α = 4 o . The most unstable oscillatory mode has an eigenvalue ω = 0.01942 − 0.05512i. The Global mode is temporally stable because largest ω i < 0. Figure 13 shows the two dimensional mode structure for the N=1, Re=523 and semi-cone angle α = 4 o . It has been observed that the disturbances evolved in the flow field at the earlier stage than that of a cone with semi-cone angle α = 2 o . The magnitudes of the disturbance velocity components and the region of contamination in polar direction are higher then that of α = 2 o , which makes the flow convectively more unstable. The largest ω i has reduced with the increased semi-cone angle α for a given Reynolds number, which makes the Global modes more stable. C. Semi-cone angle α = 6 o Figure 14 shows the spectrum for helical mode N=1, Re=698 and semi-cone angle α = 6 o . The most unstable oscillatory mode has an eigenvalue ω = 0.02539 − 0.07707i. The Global mode is temporally stable because largest ω i < 0. Figure 15 shows the two dimensional mode structure for the N=1, Re=698 and semi-cone angle α = 6 o . The structure of discrete part of the spectrum disturbs with the increase in semi-cone angle α. The magnitudes of the u r disturbance amplitudes is one order higher than that of u θ and u Φ . With the increase in semi-cone angle α the disturbances starts to grow at early stage. It has been observed that the disturbances evolved in the flow field at the earlier stage than that of a cone with semi-cone angle α = 2 o . The spatial structure of the disturbances are of similar nature, however damping rate of the Global mode increases with the increase in semi-cone angle.
D. Temporal growth rate Figure 16 shows the temporal growth rate of the least stable eigenmodes for different Reynolds number and semi-cone angles α. The growth rate of eigenmodes increases with the increase in Reynolds number for all azimuthal wave-numbers(N) and semi-cone angles (α). For the Range of Reynolds number and semi-cone angles the least stable eigen modes are negative, hence the Global modes are stable. The Global modes with helical mode, N = 1 are least stable for α = 2 o , α = 4 o and α = 6 o . The damping rate of Global modes with N = 2 is higher than N = 0 for α = 2 0 . At smaller Re, N = 0 is more stable than N = 2, however at higher Re, N = 2 is more stable then N = 0 for semi-cone angle α = 4 0 . The helical modes N = 3, 4 and 5 have larger damping rates then that of N = 0, 1 and 2 for all Re and α. The Global modes are more stable at higher semi-cone angle α. The transverse curvature reduces in the stream-wise direction  with the increase in Reynolds number for same α. The favorable pressure gradient is higher at higher semi-cone angle α. The Global modes are more stable at higher semi-cone angle, proves that favorable pressure gradient has damping effect on the Global modes.

E. Spatial amplification rate
The global temporal modes exhibits growth/decay in the stream-wise direction while moving towards the downstream. This growth/decay at different stream-wise station can be quantified by spatial amplification rate (A x ). The spatial amplification rate (A x ) shows the growth of all the disturbances together in the streamwise direction. Figure 17 shows the A x for different azimuthal wave-numbers(N) and semi-cone angles(α) for Re = 349. The spatial growth rate(A x ) increases with the increase in semi-cone angle(α) for all the azimuthal wave-numbers for given Reynolds number.We computed it for Re = 174 to Re = 1047, however the result are presented for Re = 349 only. As semi-cone angle(α) increases the growth rate of the disturbances increases in the stream-wise direction and makes the flow convectively unstable. Hence, at higher semi-cone angle flow becomes convectively more unstable. For N = 0, 1 and 2 the disturbances exhibits higher spatial growth rate at small Reynolds number, however for N = 3, 4 and 5 spatial growth rate is high for higher Reynolds number.

F. Conclusions
The linear global stability analysis of boundary layer forms on a circular i cone is performed. The combined effect of transverse curvature and pressure gradient has been studied. The Global temporal modes are computed for the axisymmetric boundary layer on a circular cone for the range of Reynolds number from 174 to 1046 with azimuthal wave-number, N from 0 to 5 and semi-cone angle α=2 o ,4 o & 6 o . The largest imaginary part (ω i ) of the computed Global modes are negative for all the Reynolds numbers(Re), azimuthal wave-numbers(N) and semi-cone angles(α). Hence, the Global modes are temporally stable. The wave-like behaviour of the eigenmodes are found in the stream-wise direction. The wavelength of the wavelet structure reduces with the increase in frequency(ω r ). The 2D spatial structure of the Global modes show that the size and magnitude of the disturbance amplitudes increases while moving towards downstream, which proves that the flow is convectively unstable. The damping rate of the disturbances increases with the increase in semi-cone angle(α) from 2 o to 6 o . Thus, Global modes are more stable at the higher semi-cone angles(α). At the same time with the increase in (α) from 2 o to 6 o , the spatial growth rate (A x ) also increases at a given Reynolds number, thus the flow becomes convectively more unstable at the higher semi-cone angles.
The azimuthal wave-numbers N = 3, 4 and 5 have less temporal growth(ω i ) and spatial growth(A x ) compare to N = 0, 1 and 2. The azimuthal wave-number N = 1 is found to be least stable one for all the Re and α. The increase in the semi-cone angle(α) develops favorable pressure gradient and reduces transverse curvature effect. Thus, the favorable pressure gradient stabilizes the Global modes and effect of transverse curvature reduces. Thus, the role of favorable pressure gradient is more effective than that of transverse curvature on flow stability.