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 (). 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 , , and azimuthal wave numbers from 0 to 5. It is found that the Global modes are more stable at higher semi-cone angle , due to the development of favorable pressure gradient. The effect of transverse curvature is reduced at higher semi-cone angles (). The spatial structure of the eigenmodes show that the flow is convectively unstable. The spatial growth rate () increases with the increase in semi-cone angle () from to . Thus, the effect of increase in semi-cone angle () is to reduces the temporal growth rate () and to increases the spatial growth rate () of the Global modes at a given Reynolds number.
pacs:Valid PACS appear here
Also working at] Dept. of Mechanical Engineering, Faculty of Engineering, Marwadi Education Foundations Groups of Institutions, Rajkot, INDIA.
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 Drazin and Raid (2004). 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 two-dimensional base flow is called Global stability analysis Theofilis (2003). 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 Rao (1974) first studied the stability of axisymmetric boundary layer. He found that non-axisymmetric disturbances are less stable than that of two-dimensional disturbances. Tutty Tutty and Price (2002) investigated that for non-axisymmetric modes critical Reynolds number increases with the azimuthal wave number. The critical Reynolds number found to be 1060 for mode and 12439 for mode. The axisymmetric mode is found least stable fourth mode. Vinod Vinod (2006) investigated that higher non-axisymmetric mode 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 Malik and Poll (1985) 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 Vinod and Govindarajan (2012). 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 Maslov (1984). It has been confirmed by the experiments that the two dimensional Mack mode dominates in hypersonic boundary layer Stetson and Kimmel (1993a, b); Maslov et al. (2002). Experimental results show that leading bluntness affects the transition location significantly in circular cone Horvath et al. (2002). 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 () 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 Akervik et al. (2008); Alizard and Robinet (2007); Ehrenstein and Gallaire (2005). 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
and body radius of the cone 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
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 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,
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.
r = radial coordinate
= polar coordinate
= azimuthal coordinate
= circular frequency
= azimuthal wave-number
ii.1 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.
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 Theofilis (2003), 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 Fasel et al. (1990). 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 Swaminathan et al. (2011); Theofilis (2003). Thus, we considered linear extrapolated conditions at the outlet for the numerical solution of the general eigenvalues problem.
similarly, one can write extrapolated boundary conditions for polar and azimuthal disturbance components and 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 Theofilis (2003).
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 (Malik (1990)).
In the above grid stretching method half number of the collocation points are concentrated within the angle from the lower surface 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 Costa et al. (2007).
The value of is chosen carefully to improve the spatial resolution in the stream-wise direction. Very small value of 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 Costa et al. (2007). To incorporate the effect of physical dimensions of the domain [, ] 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 and are square matrix of size , is an eigenvalues and is a vector of unknown disturbance amplitude functions ,, and . The above mentioned all the boundary conditions are properly incorporated in the matrix A and B.
ii.2 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 ,using LU decomposition. Full spectrum method is employed for this small subspace to get good approximate solution of the original general eigenvalues problem Theofilis (2011). 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. Given the large subspace size of ,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@.
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 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 , and . The steady Navier-Stokes equations are solved using SIMPLE algorithm with under-relaxation to get the stable solution. The spatial discretization of the N-S equations are done using a Quick scheme, which is a weighted average of second order upwind and second order central scheme. Thus obtained base velocity profile is interpolated to a spectral grids 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 Roache (1994). 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 1. 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 of the solution while increases the time for the computations. The distribution of grid is geometric in both the directions. The computed convergence order for U and V are 1.98 and 2.01 are in good agreement with the second order discretization scheme used in the finite volume code ANSYS FLUENT. Sufficient large domain height is selected in radial direction, i.e 0.5m. Thus the Mesh 1 is used in all the results presented in to calculate velocity profile for basic state. Figure 2 and 3 show the contour plot and velocity profile on the transformed coordinates and . The velocity profile is qualitatively similar to that of Garratt(2006) Garret and Peak (2007). 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, 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. and 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 Bhoraniya and Vinod (2016).
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 and . 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 (). The stream-wise length is normalized by cone radius(a). The domain size in polar(wall normal) direction is taken as . 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 2 shows the values of two leading eigenvalues computed for , and semi-cone angle 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 2 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.
v.1 Semi-cone angle
Figure 8 shows the eigenspectrum of axisymmetric mode () for and semi-cone angle . The eigenmodes marked by square and circle are called stationary and oscillatory mode. The stationary mode has a complex frequency . The Global mode is temporally stable because 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 () and polar () 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 disturbances are one order higher then that of . The polar disturbances contaminates the flow field upto large extent than that of disturbances, however its magnitudes very small compare to .
Figure 10 shows the spatial structure of oscillatory eigenmode. The associated complex frequency for this mode is . This oscillatory Global mode is also temporally stable, because . 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 convectively unstable.
Figure 11 shows the two dimensional mode structure of the and for the . 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 . 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.
v.2 Semi-cone angle
Figure 12 shows the spectrum for helical mode N=1, Re=523 and semi-cone angle . The most unstable oscillatory mode has an eigenvalue . The Global mode is temporally stable because largest . Figure 13 shows the two dimensional mode structure for the N=1, Re=523 and semi-cone angle . 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 . The magnitudes of the disturbance velocity components and the region of contamination in polar direction are higher then that of , which makes the flow convectively more unstable. The largest has reduced with the increased semi-cone angle for a given Reynolds number, which makes the Global modes more stable.
v.3 Semi-cone angle
Figure 14 shows the spectrum for helical mode N=1, Re=698 and semi-cone angle . The most unstable oscillatory mode has an eigenvalue . The Global mode is temporally stable because largest . Figure 15 shows the two dimensional mode structure for the N=1, Re=698 and semi-cone angle . The structure of discrete part of the spectrum disturbs with the increase in semi-cone angle . The magnitudes of the disturbance amplitudes is one order higher than that of and . 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 . 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.
v.4 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, are least stable for , and . The damping rate of Global modes with is higher than for . At smaller Re, is more stable than , however at higher Re, is more stable then for semi-cone angle . The helical modes , 4 and 5 have larger damping rates then that of , 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.
v.5 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 (). The spatial amplification rate () shows the growth of all the disturbances together in the stream-wise direction. Figure 17 shows the for different azimuthal wave-numbers(N) and semi-cone angles() for . The spatial growth rate() increases with the increase in semi-cone angle() for all the azimuthal wave-numbers for given Reynolds number.We computed it for to , however the result are presented for 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 , 1 and 2 the disturbances exhibits higher spatial growth rate at small Reynolds number, however for , 4 and 5 spatial growth rate is high for higher Reynolds number.
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 =, & . The largest imaginary part () 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 wave-length of the wavelet structure reduces with the increase in frequency(). 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 to . Thus, Global modes are more stable at the higher semi-cone angles(). At the same time with the increase in () from to , the spatial growth rate () also increases at a given Reynolds number, thus the flow becomes convectively more unstable at the higher semi-cone angles. The azimuthal wave-numbers , 4 and 5 have less temporal growth() and spatial growth() compare to , 1 and 2. The azimuthal wave-number 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.
- preprint: APS/123-QED
- footnotetext: Where, , ,
- footnotetext: The dimensions of computational domain in spherical coordinates are and for Global stability computations.
- P. Drazin and W. Raid, Cambridge University Press (2004).
- V. Theofilis, Progress in Aerospace Sciences 39, 3229 (2003).
- G. Rao, Journals of applied Mathematics and Physics 25, 251 (1974).
- O. Tutty and W. Price, Physics of Fluid 14, 628 (2002).
- N. Vinod, Ph.D. Thesis, Engineering Mechanics Unit, JNCASR (2006).
- M. R. Malik and D. I. A. Poll, AIAA Journal 23, 1362 (1985).
- N. Vinod and R. Govindarajan, J.Fluid Eng. 134 (2012).
- L. M. Maslov, AGARD Report 709 (1984).
- K. F. Stetson and R. L. Kimmel, AIAA J 92, 0737 (1993a).
- K. F. Stetson and R. Kimmel, AIAA J 93, 0896 (1993b).
- A. A. Maslov, S. Mironov, and A. A. Shiplyuk, AIAA J p. 0153 (2002).
- T. J. Horvath, S. A. Berry, B. R. Hollis, and B. A. Chang, AIAA J 2743, 0153 (2002).
- E. Akervik, U. Ehrenstein, F. Gallaire, and D. S. Henningson, Eur. J. Mech. B/Fluids 27, 501 (2008).
- F. Alizard and J.-C. Robinet, Physics of Fluid 19, 3229 (2007).
- U. Ehrenstein and F. Gallaire, J. Fluid. Mech. 536, 209 (2005).
- V. Theofilis, Progress in Aerospace Sciences 39, 3229 (2003).
- H. Fasel, U. Rist, and K. U., AIAA J 28(1), 29 (1990).
- G. Swaminathan, K. Sahu, A. Sameen, and R. Govindarajan, Theor. Comput. Fluid Dyn. 64, 25 (2011).
- M. Malik, Journal of Computational Physics 86(2), 376 (1990).
- B. Costa, D. W., and A. Simas, Proc. SCPDE: Recent Progress in Scientific Computing pp. 179–188 (2007).
- V. Theofilis, Annu. Rev. Fluid Mech. 43, 319 (2011).
- P. J. Roache, J. Fluid. Mech. 19, 405â413 (1994).
- S. J. Garret and N. Peak, Eur. J. Mech. B/Fluids 26, 344 (2007).
- R. Bhoraniya and N. Vinod, Theor. Comput. Fluid Dyn. submitted for publication (2016).