Effect of slip on the linear stability of the rotating disk boundary layer

The linear stability of the rotating disk boundary layer with surface roughness is investigated. Surface roughness is modeled using slip boundary conditions [M. Miklavcˇicˇ and C. Y. Wang, Z. Angew. Math. Phys. 55 , 235–246 (2004)], which establish concentric grooves, radial grooves, and isotropic roughness. The effect on the stationary crossﬂow and Coriolis instabilities is analyzed by applying slip conditions to the undisturbed ﬂow and linear disturbances. This analysis builds on the work of Cooper et al. [Phys. Fluids 27 , 014107 (2015)], who modeled slip effects on the base ﬂow but applied the no-slip condition to the linear perturbations. Neutral stability curves and critical parameter settings for linearly unstable behavior are computed for several radial and azimuthal slip length settings. The application of slip on the linear disturbances has a signiﬁcant impact on the ﬂow stability. In particular, the Coriolis instability undergoes considerable destabilization in the instance of concentric grooves (i.e., radial slip) and radial grooves with sufﬁciently large azimuthal


I. INTRODUCTION
The rotating disk boundary layer is an archetype example for investigating instability mechanisms that trigger laminar-turbulent transition in wall-bounded external flows. When the disk rotates about its axis with a constant angular velocity, a fluid is drawn toward the disk before being thrown radially outward. The resulting flow admits the similarity solution first derived by von K arm an 3 and has since been the study of numerous theoretical and experimental investigations. 4,5 The influential work by Gregory et al. 6 demonstrated that the rotating disk is susceptible to an inviscid crossflow instability, much like that found in the boundary layers on swept wing bodies. The crossflow or type-I instability is convective and takes the form of stationary spiral vortices. Following the initial observations of Gregory et al., 6 crossflow disturbances have been observed by many experimentally [7][8][9][10][11][12] and emerge for Reynolds numbers Re < 300, with transition to turbulence developing on the interval 500 < Re < 560. [A formal definition for the Reynolds number, Re, is given below in Eq. (5).] Malik 13 employed linear stability theory to predict the onset of the crossflow instability, with the critical Reynolds number Re c % 286.
In addition to the type-I instability, the rotating disk is susceptible to a viscous type-II mode brought about by Coriolis effects 14 and a spatially damped type-III mode. 15 Malik 13 and Balakumar and Malik 16 showed that the stationary type-II instability emerges for Reynolds numbers Re տ 451. Further investigations of the rotating disk concern the development of absolute instability 11,17 and global instability mechanisms. [18][19][20][21][22][23] Due to its similarities with the flow on a swept wing, the rotating disk is often used to test strategies for controlling laminar-turbulent transition processes, including mass transfer, [24][25][26] compliant surfaces, 27,28 temperature-dependent viscosity, 29 axial flow, 30 and disk modulation. 31,32 In addition, small-scale surface roughness, as a method for delaying the onset of crossflow, has been the focus of several recent investigations. 2,[33][34][35][36][37][38][39] Roughness on the surfaces of airplane wings and rotating bodies, such as the rotating disk, was thought to enhance drag and trigger the premature onset of laminar-turbulent transition. However, recent research suggests that the right-sort-of roughness can establish favorable control benefits. [40][41][42] Currently, two distinct models for small-scale roughness on the disk surface exist in the literature, which permit the study of concentric, radial, and isotropic roughness configurations. Concentric roughness refers to surface grooves along the azimuthal direction, e.g., a phonograph record or laser-etched disk. Radial roughness refers to radial grooves that form a star pattern on a disk. The first approach, by Yoon et al. 43 (YHP), implements concentric roughness by prescribing the surface geometry as a function of the radial position. The linear stability of the YHP model was undertaken by Garrett et al. 33 The second method, by Miklavčič and Wang 1 (MW), is the focus of the current investigation, where a so-called partial-slip condition replaces the usual no-slip boundary condition on the disk surface.
The MW approach for modeling surface roughness makes use of the slip boundary condition 44 for a velocity field, u, slip length, k, and wall-normal direction, z. According to this condition, the fluid's tangential velocity near the disk surface, is directly proportional to the shear rate at the surface. (In the instance k ¼ 0, the no-slip condition is recovered.) In this way, Miklavčič and Wang 1 modeled radial and azimuthal slip, leading to concentric and radial grooves, respectively. Moreover, coupling nonzero radial and azimuthal slip lengths establish isotropic roughness arrangements.
Cooper et al. 2 undertook a linear stability analysis of the rotating disk with small-scale surface roughness imposed using the MW partial-slip conditions. Slip applied along the radial direction (matched to concentric grooves) stabilized both stationary type-I crossflow and type-II Coriolis instabilities, raising the critical Reynolds number to significantly larger values than in the case of no-slip. Indeed, the latter mode disappeared from the neutral stability curves for relatively small radial slip lengths. On the other hand, slip along the azimuthal direction (matched to radial grooves) raised the critical Reynolds number of the stationary crossflow instability but destabilized the Coriolis instability to the extent that this mode soon became the dominant form of instability.
Following the initial study by Cooper et al., 2 the MW modeling was extended to include the asymptotic structure of the inviscid crossflow mode 35 and local and global stability mechanisms. 37 The BEK family of rotating flows (that includes the B€ odewadt, Ekman, and von K arm an boundary layers), 34 rotating cones, 39 and non-Newtonian flow effects 36 were also included. In these earlier investigations, the surface slip was applied only to the undisturbed flow to model flow over a rough surface. Crucially, these studies neglected slip effects in the linear stability analysis. Instead, no-slip conditions were imposed on the linear perturbations. As a result, fundamental stability characteristics were overlooked.
In the context of plane Poiseuille flow in a two-dimensional channel, several studies [45][46][47] have shown that applying slip to the undisturbed flow and linear perturbations has significant stability implications. In particular, Ceccacci et al. 47 showed that the slip length required to make the flow linearly stable depended on whether the linear stability analysis was subject to no-slip or slip conditions. The current study aims to extend the earlier investigation of Cooper et al. 2 by applying slip to both the undisturbed flow and linear perturbations. Slip configurations matched to concentric and radial grooves and isotropic roughness are considered, enabling us to determine the correct effect of surface roughness on the stability of the rotating disk boundary layer.
The remainder of this paper is outlined as follows. In the following section, the formulation is described, including the model, undisturbed flow, and linearized Navier-Stokes (LNS) equations. In Secs. III and IV, linear stability and energy analyses are undertaken on the rotating disk with slip, including conditions matched to concentric grooves, radial grooves, and isotropic roughness. Finally, conclusions are presented in Sec. V.

II. FORMULATION A. Undisturbed flow
Consider a disk of infinite radius that rotates in an infinite body of incompressible fluid at a constant angular velocity X Ã ¼ ð0; 0; X Ã Þ. The modeling is carried out in cylindrical polar coordinates r Ã ¼ ðr Ã ; h; z Ã Þ, where r Ã is the radial distance from the vertical axis of rotation, h is the azimuthal angle, and z Ã is the wall-normal direction. (Here, an asterisk denotes dimensional quantities.) In what follows, all quantities are defined in the frame of reference that rotates with the disk.
The Navier-Stokes equation and the continuity equation are considered in cylindrical polar coordinates and are defined as for the kinematic viscosity, Ã , and fluid density, q Ã . Here, t Ã denotes dimensional time, U Ã ¼ ðU Ã ; V Ã ; W Ã Þ is the dimensional velocity, and P Ã is the dimensional pressure. The von K arm an 3 similarity solution is found by setting and where F, G, and H represent the non-dimensional velocity profiles along the radial, azimuthal, and wall-normal directions, respectively. The non-dimensional pressure field is denoted P and the scaled wall- is the constant boundary-layer thickness that is used to non-dimensionalise units of length. Additionally, the non-dimensional radius r ¼ r Ã =d Ã .
On substituting Eq. (2) into Eq. (1), the following system of ordinary differential equations for F, G, H, and P is derived: which is solved subject to the boundary conditions and Here, a prime denotes differentiation with respect to z, and parameters k and g represent the slip lengths acting along the radial r-and azimuthal h-directions. In reference to those earlier studies that used slip to model surface roughness, 1,2,34-39 setting k > 0 and g ¼ 0 establishes concentric grooves, g > 0 and k ¼ 0 radial grooves, and k ¼ g > 0 isotropic roughness. (It is worth noting that in the earlier investigations undertaken by Cooper et al. 2 and others, [34][35][36][37] the relationship between the slip length and the form of roughness was misinterpreted. There, they assigned the radial slip length, k, to radial or star roughness and the azimuthal slip length, g, to concentric roughness. The reverse of that modeled herein. However, as pointed out by Thomas (personal communication, 2021), the interpretation by Cooper and coworkers was incorrect, as the radial slip length, k, imposes roughness along the radial r-direction and establishes concentric roughness when the disk is viewed from above. Similarly, the azimuthal slip length, g, establishes roughness along the azimuthal h-direction, leading to the radial or star roughness when viewed from above. Hence, the differences in terminology are presented here and that given in the previous investigations. The system of ordinary differential equations (3) was solved using the boundary value routines in MATLAB, namely, bvp4c. Figure 1 depicts steady velocity profiles F, G, and H as functions of the wall-normal z-direction, for different combinations of k and g.
On non-dimensionalizing the velocity and pressure fields on the respective scales r Ã o X Ã and q Ã r Ã2 o X Ã2 , the undisturbed base flow can be written as where the Reynolds number of the flow is defined as for a non-dimensional reference radius r o .

B. Linearized Navier-Stokes equations
The undisturbed flow (4) is perturbed by infinitesimally small disturbances, with the total velocity U and pressure P given as where e ( 1 and u ¼ ðu; v; wÞ: Substituting Eqs. (6) into Eq. (1) and linearizing with respect to e gives the linearized Navier-Stokes (LNS) equations and the perturbation continuity equation where Boundary conditions on the disk surface are modeled in two distinct ways. The first approach, method A, follows the modeling implemented by Cooper et al. 2 among others, whereby perturbations at the disk wall satisfy the no-slip condition However, this particular approach neglects the effect of slip on the perturbations [and is only applied to the undisturbed flow (4)].
Consequently, a second approach, method B, is devised, in which perturbations satisfy the slip condition at the disk wall Here, k and g are the same valued slip lengths as that specified in the boundary conditions (3e) and (3f) for the undisturbed flow (4). In the far-field, perturbations of both methods satisfy The motivation for method B is based on those studies undertaken on plane Poiseuille flow in a two-dimensional channel, [45][46][47] where slip was applied to the total velocity field, U, i.e., both the undisturbed base flow, U B , and the linear perturbation, u. Applying slip conditions to the total velocity field accounts for the potential slip effects brought about by both the undisturbed base flow and the perturbation and accurately captures the stability of the boundary layer.
The commonly used parallel flow approximation is implemented, whereby the radius, r, is replaced with the Reynolds number, Re. Additionally, as Re ) 1, terms of order Re À2 and smaller are neglected, as per Malik. 13 Replacing r by Re in Eqs. (7) establishes a set of equations that are separable with respect to r, h, and t, which allow perturbations to be decomposed into the normal mode form fu; pgðr; h; z; tÞ ¼ fûðzÞ;pðzÞg exp i ar þ bReh À xt ð Þ ð Þ ; (11) where a 2 C denotes the radial wavenumber, x 2 R is the frequency, and b 2 R represents the azimuthal wavenumber. (For a spatial linear stability analysis, a i < 0 corresponds to unstable behavior.) Given the above-mentioned assumptions, Eqs. (7) becomes where j 2 ¼ a 2 þ b 2 and a ¼ a À i=Re.

III. LINEAR STABILITY
The linear stability [Eqs. (12)] was numerically solved using the Chebyshev collocation method developed by Trefethen. 48 Derivatives along the wall-normal z-direction were replaced with Chebyshev matrix approximations, with N Chebyshev mesh points mapped from the semi-infinite physical domain z 2 ½0; 1Þ onto a finite computational interval n 2 ½À1; 1, via the coordinate transformation for a stretching parameter, L. In addition, the system of Eq. (12) was recast as a quadratic eigenvalue problem of the following form: forq ¼ ðû;v;ŵ;pÞ T and square matrices A j (for j ¼ 0, 1, 2). The eigenvalue problem was then solved using the eig command in MATLAB. Table I compares the radial wavenumber a computed for the slip lengths k ¼ g ¼ 0:5 and method B. In each instance, the Reynolds number Re ¼ 500, azimuthal wavenumber b ¼ 0:08, and frequency x ¼ 0. The number of Chebyshev mesh points increases from N ¼ 25 to N ¼ 200 for stretching parameters L ¼ 1, 2, and 3. Results are identical for all three values of L modeled when N ! 50. For the remainder of this study, N ¼ 100 and L ¼ 2, were sufficient to undertake an accurate linear stability investigation. A spatial linear stability analysis of the stationary type-I crossflow and type-II Coriolis instabilities was undertaken by setting the frequency x ¼ 0. Neutral points ðRe; b; aÞ for linear instability were then determined at azimuthal wavenumber increments Db ¼ 10 À4 , for various combinations of the two slip lengths, k and g. In the subsequent discussion, the hat notation introduced in Eqs. (11), (12), and (14) has been dropped.
In the instance k ¼ g ¼ 0 (i.e., the von K arm an flow without slip), the neutral stability curve is characterized by two distinct branches, with the upper branch matched to the type-I crossflow instability and the lower branch to the type-II Coriolis instability. In addition, critical conditions for linear instability are consistent with those results obtained by Malik 13 and Balakumar and Malik 16 and others: the critical Reynolds number Re c % 286:1 for the crossflow instability, whereas Re c % 451:4 for the Coriolis instability.
Neutral stability curves and critical conditions obtained for k > 0 and method A are consistent with those results presented by Cooper et al. 2 A strong stabilizing effect is established for both the type-I and type-II modes, with neutral curves moving to the right. The critical Reynolds number, Re c , associated with the type-I mode increases rapidly with increasing k, with Re c more than 500 for k ¼ 0:5, and the lower branch of the neutral stability curve (corresponding to the type-II mode) disappears for k > 0:05.  Physics of Fluids ARTICLE pubs.aip.org/aip/pof The method B results presented in Fig. 2(b) behave very differently from those calculations obtained for method A. The application of slip on the linear perturbations [in addition to the undisturbed flow (4)] reduces the critical Reynolds number for linearly unstable behavior, with neutral curves moving to the left and upward in the ðRe; bÞplane. Consequently, the effect of slip applied along the radial r-direction (i.e., k > 0) is to destabilize the rotating disk boundary layer. Figure 2(c) displays neutral stability curves for method B over a reduced ðRe; bÞ-parameter range and radial slip lengths k ¼ 0:05 through k ¼ 0:25 at step intervals Dk ¼ 0:05. The plot demonstrates that the destabilizing effect is primarily due to the type-II Coriolis instability; the lower branch moves to the left, and smaller Reynolds numbers, Re, at a considerably faster rate than the upper branch. Moreover, the Coriolis instability becomes the most unstable mode for a radial slip length on the interval 0:15 < k < 0:2.
Critical conditions Re c , b c , and a c , for linearly unstable behavior, are plotted in Fig. 3 as a function of the radial slip length, k, for step sizes Dk ¼ 10 À2 . Blue solid curves depict the critical conditions obtained for the type-I crossflow instability using method A. Red dashed and yellow chain lines illustrate the equivalent solutions obtained via method B for the type-I and type-II modes, respectively. (Curves matching the method A type-II Coriolis instability are not included in Fig. 3 as the critical conditions for linear instability could not be determined for k ! 0:05.) The strong stabilizing effect associated with the method A approach is evident in Fig. 3(a), with the critical Reynolds number, Re c , increasing over the range of k shown. On the other hand, the destabilizing effect brought about by method B is demonstrated, with Re c decreasing for both forms of linear instability. On the interval 0 k 0:18, the critical Reynolds number, Re c , for the crossflow instability is the smallest. Thus, the type-I crossflow mode is the most unstable for this range of k. However, near k % 0:18 (and Re c % 270), the dominant mode (i.e., most linearly unstable) switches to the Coriolis instability. In addition, for all slip lengths k > 0:18, critical conditions for the crossflow instability could not be determined. Hence, the abrupt end in the red dashed lines at k ¼ 0:18. The critical Reynolds number, Re c , for the Coriolis instability continues to decrease for larger k, with a minimum value of Re c % 179 attained near k % 0:7; for k > 0:7, Re c marginally increases. Furthermore, for g ¼ 0, the critical azimuthal and radial wavenumbers, b c and a c , grow in size as k increases.
Critical Reynolds numbers, Re c , azimuthal wavenumbers, b c , and radial wavenumbers, a c , are presented in Tables II and III for the respective methods A and B, for different combinations of the slip lengths k 2 ½0; 1 and g 2 ½0; 1. The results of method A, presented in Table II, are in excellent agreement with those given in Cooper et al. 2 and Alvero glu et al. 34 (see the Appendix for direct comparisons). Figure 4 compares the absolute value of the radial u-, azimuthal v-, and wall-normal w-velocity perturbation fields established for methods A (blue solid lines) and B (red dashed). The two sets of solutions correspond to the slip lengths k ¼ 0:2 and g ¼ 0, and flow conditions Re ¼ 400 and b ¼ 0:1. In each instance, solutions correspond to the type-I crossflow instability and are normalized on the maximum amplitude of the v-velocity field, i.e., maxjvj. The resulting radial wavenumbers for methods A and B are a ¼ 0:3835 À 0:0107i and a ¼ 0:3588 À 0:0530i, respectively, i.e., both modes are linearly unstable. However, the latter approach establishes stronger radial growth. (Recall that linearly unstable behavior corresponds to a i < 0.) The general structure of all three velocity fields is qualitatively the same for the two methods. The notable exception is the behavior near the disk wall. The radial u-velocity is non-zero at z ¼ 0 for method B due to the slip boundary condition (9). In addition, there are subtle differences in the character of the azimuthal v-and wall-normal wvelocity fields near the disk surface. Most notably, the w-velocity displays a non-zero gradient at z ¼ 0, which is a direct consequence of the slip condition (9) and the requirement that the continuity [Eq. (12d)] be satisfied at z ¼ 0, i.e., w 0 ð0Þ ¼ Àiaku 0 ð0Þ for method B.
[For method A, based on the application of the no-slip condition (8) w 0 ¼ 0 at z ¼ 0, as shown in Fig. 4(c)].
B. Radial grooves: g ! 0 and k 5 0 The above-mentioned analysis is extended to include the modeling of slip lengths g ! 0 and k ¼ 0. This particular combination of ðk; gÞ establishes a roughness model for radial grooves. 2  Comparing the solutions of the two approaches, the upper branch (matched to the type-I crossflow mode) appears to be relatively unchanged in most cases. The noticeable exception is the appearance of a new lobe located in the upper right-hand corner of Fig. 5(b) that protrudes from the upper branch of the neutral curve for the case g ¼ 0:5 (turquoise dashed line).
The appearance of this new lobe is surprising, but expanding the ðRe; bÞ-plane [see Fig. 5(c)], this particular feature of the neutral curve is found for all five non-zero g cases modeled. As the slip length, g, decreases, the lobe moves toward larger Reynolds numbers, Re, and larger azimuthal wavenumbers, b, enveloping a greater region of the ðRe; bÞ-plane. For sufficiently large Re, the lobe reattaches to the upper branch of the neutral curve. It is worth noting that the lobe was neither found for the rotating disk boundary layer without slip (that is, k ¼ g ¼ 0) nor was it found for any configuration modeled using the method A approach. Moreover, the lobe was not found, at least for the range of Re and b considered, for radial slip lengths k > 0. Another distinguishing difference between those neutral stability curves plotted in Fig. 5(a) for method A and those results displayed in Fig. 5(b) for method B is the behavior of the type-II Coriolis instability along the lower branch. For method A, as g increases, the lower branch of the neutral curve progressively shifts to the left and smaller Reynolds numbers Re. However, for method B, the lower branch moves to the right and larger Re for g ¼ 0:1, but for g > 0:1, neutral stability curves shift to the left and smaller Re. Figure 6 displays critical conditions Re c , b c , and a c for linear stability, as a function of the azimuthal slip length, g, with Dg ¼ 10 À2 . Line types are the same as those given in Fig. 3, with purple dotted lines representing solutions obtained via method A for the type-II Coriolis instability. The two solutions obtained for the type-I mode are almost identical over the interval g 2 ½0; 0:55; the critical Reynolds number, Re c , increases for both methods A and B. (For g > 0:55, Re c could not be determined for the method B crossflow instability due to the significant destabilization of the Coriolis instability.) On the other hand, there are considerable differences between the two sets of results for the type-II mode; Re c decreases continuously for method A, whereas for method B, Re c initially grows and achieves a peak value of Re c % 591 at g % 0:13 before decreasing rapidly for larger g. Indeed, about g % 0:46, the Coriolis instability becomes more unstable than the crossflow instability and achieves a value of Re c % 184 at g ¼ 1.
(Note this is comparable with the behavior observed for the method A approach, in which the Coriolis instability is the most unstable mode  for g տ 0:4.) Consequently, to optimize the stabilization of the rotating disk boundary layer via radial roughness, the azimuthal slip length should not exceed g % 0:46. The new lobe found along the upper branch of those neutral curves plotted in Figs. 5(b) and 5(c) raises the question of the accuracy and robustness of the numerical scheme and the stability characteristics associated with this feature. Regarding the numerical accuracy, the number of Chebyshev points, N, was raised to N ¼ 200. In addition, stretching parameters L ¼ 1 and L ¼ 3 were modeled. However, results were not affected by either of these parameter changes. Moreover, a second numerical formulation, based on a finite wall-normal z-direction (i.e., z 2 ½0; z max for z max ¼ 20, rather than the semi-infinite domain modeled herein) and Gauss-Lobatto grid spacing, 34 verified the existence of this lobe structure for non-zero g and k ¼ 0. Hence, confirming that the lobe is not a numerical artifact of the numerical scheme. Concerning the physical properties of the lobe, additional stability characteristics are presented in Figs. 7 and 8 to determine the nature of this aspect of the stability parameter space.
Contours of the radial growth rate a i 2 ½À0:09; 0 are plotted in the ðRe; bÞ-plane in Fig. 7(a). The azimuthal slip length g ¼ 0:5 (and k ¼ 0). [The dotted straight line at b ¼ 0:12 marks the location where results are analyzed in Fig. 7(b)]. The outer blue solid line depicts the critical conditions for linear instability (i.e., a i ¼ 0), while the innermost purple chain line corresponds to a i ¼ À0:09. Interestingly, within the ðRe; bÞ-region bounded by the lobe, a i > À5 Â 10 À3 . Thus, the radial growth rate within this feature of the neutral curve is small. Indeed, this is significantly smaller than the maximum growth rate realized by the crossflow instability at equivalent Reynolds numbers Re; the type-I mode attains a maximum growth rate a i % À0:08 at Re % 1000. Furthermore, the maximum growth rate found within the lobe is less than that found along the lower branch; on inspection of the contours, the Coriolis instability achieves growth rates a i տ À 0:02 (which is also considerably smaller than the radial growth observed for the crossflow instability).
The small radial growth rate, a i , found within the lobe is illustrated in Fig. 7(b). Real (blue solid lines) and imaginary (red dashed) components of the radial wavenumber, a, are plotted against the Reynolds number, Re, for b ¼ 0:12, i.e., an azimuthal wavenumber, b, which passes through the center of the lobe. Linearly unstable behavior is realized on the interval 870 < Re < 1410, with a maximum growth rate a i % À1:4 Â 10 À3 . Figure 8 displays the real and imaginary parts of the radial wavenumber, a, as a function of the azimuthal wavenumber, b, in the instance the Reynolds number Re ¼ 1000. Slip settings are the same as those in Fig. 7: g ¼ 0:5 and k ¼ 0. Blue solid and red dashed lines depict the solutions of the respective type-I crossflow and type-II Coriolis instabilities. On 0:011 < b < 0:036, the type-II mode is linearly unstable, while the type-I mode is unstable on the interval 0:022 < b < 0:084; the crossflow instability achieves radial growth rates, a i , four times greater than the Coriolis instability. In addition, both forms of instability are linearly unstable on the interval 0:022 < b < 0:036, leading to algebraic radial growth, 17 i.e., modal coalescence of the type-I and type-II modes. Interestingly, the type-II mode is again unstable on the interval 0:087 < b < 0:149, with a i > À10 À3 . Hence, the lobe along the upper branch of the neutral curve is due to the Coriolis instability. We speculate that this new feature of the neutral stability curve is an artifact of correctly applying  slip along the azimuthal h-direction, i.e., to both the undisturbed flow (4) and linear perturbation (11) Critical stability conditions for g 2 ½0; 1 and k ¼ 0 are presented in Table III for method B. In all cases modeled, the critical Reynolds number, Re c , was located along either the upper branch (i.e., the type-I mode) or the lower branch (i.e., the type-II mode). In no instance were critical conditions found along the new lobe.
The absolute value of the three components of velocity, u, are plotted in Fig. 9 for slip lengths g ¼ 0:5 and k ¼ 0, and parameter settings Re ¼ 500 and b ¼ 0:025 (left-hand side of Fig. 9) and Re ¼ 1000 and b ¼ 0:12 (right-hand side of Fig. 9). As before, blue solid lines and red dashed lines represent results obtained for method A and method B, respectively. Both parameter sets are matched to the type-II Coriolis mode and linearly unstable behavior (i.e., a i < 0), except for the latter parameter set and method A where a i > 0. Similar to those velocity fields presented in Fig. 4 for the concentric roughness configuration, method A and method B solutions are qualitatively the same except near the disk wall. As a result of the slip condition (9) for method B, the v-velocity field is non-zero at the disk wall. Consequently, the wall-normal w-velocity field displays a non-zero gradient at the disk surface, i.e., for the continuity [Eq. (12d)] to be satisfied at z ¼ 0, w 0 ð0Þ ¼ Àibgv 0 ð0Þ is required.
C. Isotropic roughness: k 5 g ! 0 The final model investigated corresponds to isotropic roughness, where the two slip lengths are non-zero and of equal magnitude, i.e., k ¼ g > 0. Figures 10(a) and 10(b) display neutral stability curves for the respective methods A and B, for slip lengths k ¼ g ¼ 0 through to k ¼ g ¼ 0:5. Method A solutions are consistent with those results presented in earlier studies, 2 with both forms of linear instability stabilized by slip. However, solutions obtained for method B suggest the correct effect of slip on the type-I and type-II modes is less straightforward.  Critical conditions for linear instability are plotted in Fig. 11 as a function of the slip length k ¼ g. Line types are the same as those used in Fig. 3. (Similar to Fig. 3, Re c could not be determined for the method A Coriolis instability when k ¼ g > 0:05.) In both instances, the critical Reynolds number, Re c , for the crossflow instability increases with the slip length. However, Re c grows at a reduced rate for the method B approach. For instance, when k ¼ g ¼ 1; Re c % 624 for method A, while Re c % 412 for method B. Hence, a weaker stabilizing effect is realized when slip is applied to the undisturbed flow (4) and linear perturbation (11) The effect of isotropic roughness on the type-II Coriolis instability is comparable with that presented in Fig. 6 for radial grooves, i.e., g > 0 and k ¼ 0. Initially, the critical Reynolds number, Re c , increases as the slip length increases, with a maximum attained at Re c % 615 for k ¼ g % 0:4. However, for slip lengths k ¼ g > 0:4, Re c decreases, with Re c % 486 for k ¼ g ¼ 1. Indeed, given the trend of the results plotted in Fig. 11(a), the type-II instability is expected to become more unstable than the type-I instability for a slip length k ¼ g > 1. Thus, there is a limitation on the stabilizing effect realized by this form of roughness.
Linear stability characteristics presented for the isotropic roughness configuration are closer in appearance to those results established for the radial roughness geometry than the concentric roughness arrangement. This suggests that the azimuthal slip length, g, has a greater impact on stability properties than the radial slip length, k, especially when the two slip lengths are of equal size. Indeed, we might have anticipated this particular observation given the behavior of the radial and azimuthal steady velocity fields, F and G, plotted in Fig. 1. Variations in the radial velocity, F, are considerably smaller for the isotropic roughness [see Fig. 1(c)] when compared with the concentric roughness [see Fig. 1(a)]. On the other hand, the azimuthal velocity, G, experiences a marginally larger variation for the isotropic roughness [see Fig. 1(f)] than the radial roughness [see Fig. 1(e)].

IV. ENERGY ANALYSIS
An energy balance analysis is now undertaken, similar to that presented by several related studies, 2,27,29,33,36 to further our understanding of the underlying physical mechanisms brought about by correctly applying slip to both the undisturbed flow (4) and linear perturbation (11) The energy equation is derived by first multiplying the linearized momentum Eqs. (7a)-(7c) by the velocity perturbation fields u, v, and w and summing together to give the following kinetic energy equation for the linear disturbances: where the viscous stress terms where an overbar denotes a period-averaged quantity and the w subscript in (III) and (IV) denotes a quantity computed at the disk wall.  On the left-hand side of Eq. (16), (A) represents the average kinetic energy convected by the radial component of the undisturbed flow, (B) the work done by the pressure perturbation field, and (C) the work done by the viscous stresses across the boundary layer. On the righthand side, (I) represents the Reynolds-stress energy production, (II) the viscous dissipation energy removal, (III) the pressure work, (IV) the contributions from work done on the wall by viscous stresses, and (V) the streamline curvature effects and the three-dimensionality of the undisturbed flow. The energy balance Eq. (16) is normalized by the integrated mechanical energy flux to give where a i denotes the imaginary component of the radial wavenumber, a. Positive terms indicate the production of energy, while negative terms indicate the dissipation of energy. When energy production exceeds energy dissipation in the system, the linear perturbation grows and is unstable, i.e., a i < 0. Previous studies 2,27,31,33,36 showed that the Reynolds stress, P 2 , and viscous dissipation, D, dominate the energy production and energy removal in the rotating disk boundary layer without slip. Whereas terms PW 2 and S i are zero as a consequence of the boundary conditions. Moreover, all other terms are negligible. Similar conclusions are found here for method B in the case slip is applied to both the undisturbed flow (4) and linear perturbation (11) Figure 12 illustrates the radial growth rate, Àa i , as a function of the azimuthal wavenumber, b, for the Reynolds number Re ¼ 400. Solutions correspond to method B, and slip lengths k 2 ½0; 0:5 and g 2 ½0; 0:5 that establish concentric, radial, and isotropic roughnesses. There is a marked reduction in the radial growth rate, Àa i , as the azimuthal slip length, g, increases for the radial and isotropic roughness arrangements [see Figs. 12(b) and 12(c)]. For these two configurations, the type-I crossflow instability is the dominant mode, and for the range of g considered, there is a comparable reduction in the radial growth rate. (There is, however, a noticeable variation in the interval of unstable azimuthal wavenumbers, b, with solutions moving to the left and smaller b in the instance of radial roughness.) Moreover, short intervals of small amplification type-II Coriolis instability emerge in the case of radial roughness with g ¼ 0:4 and g ¼ 0:5 [see the green solid and turquoise dashed lines in Fig. 12(b)], which is consistent with the behavior shown in Fig.  5(b). On the other hand, for the concentric roughness configuration [see Fig. 12(a)], the radial growth rate, Àa i , initially grows as the radial slip length, k, increases. A peak radial growth rate, Àa i , is then realized for k ¼ 0:2, with Àa i decreasing for larger radial slip lengths, k. In this instance, the crossflow instability is dominant for k < 0:2, with the Coriolis instability responsible for stronger amplification rates for k > 0:2, which is again consistent with those neutral stability curves plotted in Fig. 2(b). Thus, while larger radial slip lengths, k, destabilize the rotating disk boundary layer (and, in particular, the Coriolis instability) and bring about smaller critical Reynolds numbers, Re c [recall Fig. 2(b) and Table III], the radial  growth rate, Àa i , is weaker than that found at smaller k. This trend continues for larger k.
The energy balance analysis, outlined above in Eqs. (15) through (17), is applied to those modes that achieve the maximum growth rate, Àa i , in Fig. 12 for Re ¼ 400, with results plotted in Fig. 13. The respective combinations of the radial wavenumber, a, and azimuthal wavenumber, b, are given in Table IV. Similar to the aforementioned investigations on the rotating disk, 2,27,31,33,36 the Reynolds stress, P 2 , and viscous dissipation, D, are again the most significant energy production and energy removal quantities, respectively. Moreover, both quantities are affected by the application of surface slip, leading to considerable variations in the overall energy change. Figures 13(b) and 13(c) demonstrate that the stabilizing effect obtained for the radial and isotropic roughnesses is a direct consequence of the significant reduction in the Reynolds stress, P 2 . In addition, for the isotropic roughness arrangement, the absolute value of the viscous dissipation, D, increases as the slip length increases, leading to further energy reductions. Figure 13(a) shows that for concentric roughness, an increasing radial slip length, k, destabilizes the flow up to a threshold value of k % 0:2. Destabilization is brought about by an increase in the Reynolds stress, P 2 , and a marginal reduction in the viscous dissipation, D. For k > 0:2, the energy production due to the Reynolds stress, P 2 , decreases, while the energy removal due to the viscous dissipation, D, increases. Thus, the rotating disk with concentric roughness experiences a stabilizing effect at the Reynolds number Re ¼ 400, providing the radial slip length is sufficiently large, i.e., k > 0:2.

V. CONCLUSIONS
A linear stability investigation has been undertaken on the rotating disk boundary layer with surface roughness, which was realized by imposing slip conditions on the disk wall. 1 Slip conditions were applied to the undisturbed flow (4) and linear perturbations (11) extending the earlier analysis undertaken by Cooper et al., 2 who modeled slip effects on the base flow only; slip was not applied to the linear perturbations in their study. The stability of the crossflow and Coriolis instabilities was computed for a range of slip conditions matched to concentric grooves (k > 0 and g ¼ 0), radial grooves (g > 0 and k ¼ 0), and isotropic roughness (k ¼ g > 0).
Modeling slip on both the undisturbed flow and linear perturbations has a considerable impact on the behavior of neutral curves in the stability plane and the critical conditions for linearly unstable behavior. In particular, the concentric roughness configuration destabilizes both forms of linear instability. Moreover, for large enough radial slip lengths, k, the type-II Coriolis instability is more unstable than the type-I crossflow instability. These observations are the direct opposite of that found by Cooper et al., 2 where both modes underwent a significant stabilization. An energy balance analysis indicates that the destabilization is due to an increase in the Reynolds stress energy production and a reduction in the viscous dissipation energy removal.
There were fewer differences between the results of the current study and those presented by Cooper et al. 2 for the radial and isotropic roughness arrangements. For the latter modeling approach, results for the crossflow instability were qualitatively similar to that obtained by Cooper and coworkers; a stabilizing effect was again realized for the isotropic roughness configuration, but the benefits were weaker than that observed in the earlier study. However, the Coriolis instability experienced a destabilization in the approach implemented herein for FIG. 13. Energy balance terms for the type-I mode, for Re ¼ 400 and method B. (a) k ! 0 and g ¼ 0, (b) g ! 0 and k ¼ 0, and (c) k ¼ g ! 0. Azimuthal and radial wavenumbers, b and a, are given in Table IV very large slip lengths. On the other hand, for the radial roughness geometry, there were negligible variations in the critical conditions for the crossflow instability, with a stabilizing effect realized in both instances. In addition, both approaches saw a destabilization of the Coriolis instability for sufficiently large azimuthal slip lengths, g, although some stabilization was realized for the current method and g Շ 0:13. Moreover, a new lobe brought about by the Coriolis instability was found along the upper branch of the neutral stability curve. However, for those parameter settings modeled herein, critical conditions for linear instability were always due to the crossflow instability or the Coriolis instability along the lower branch of the neutral curve. We now compare the results of the current study for the MW slip model with those computations obtained for the YHP 43 roughness model. The YHP model establishes concentric roughness on the rotating disk through a coordinate transformation that engineers radial oscillations in the undisturbed base flow. Garrett et al. 33 undertook a linear stability study of this model by applying a radial averaging of the flow, i.e., surface variations were neglected. For this flow configuration, the crossflow instability underwent a significant stabilizing effect. Meanwhile, the Coriolis instability was destabilized and eventually became the dominant mode (see Fig. 3 of Garrett et al. 33 ). Subsequently, Garrett et al. compared the YHP results for concentric roughness with what was assumed to be the equivalent set of results for the MW model. 2 However, as noted previously, earlier investigations concerning the MW model 2,33-37 misinterpreted the relationship between the slip length and the form of roughness. The azimuthal slip length, g, was thought to establish concentric roughness when, in fact, it brings about the radial or star roughness. Consequently, Garrett et al. 33 compared the YHP solutions with the MW results 2 for radial roughness (i.e., azimuthal slip) rather than the concentric roughness induced via the radial slip length, k.
On comparing the neutral stability curves presented in Fig. 2 for concentric roughness with those results given in Garrett et al. 33 for the YHP model (see Fig. 3 of their paper), several conclusions are drawn. In the case of method A, roughness has a comparable stabilizing effect on the type-I crossflow instability as the YHP approach, with critical conditions for linear instability shifting to larger Reynolds numbers. However, the two approaches exhibit different behavior for the type-II Coriolis instability; this form of instability disappears for the method A approach, whereas for the YHP model, the Coriolis instability undergoes a significant destabilization. On the other hand, for method B, both the Coriolis and crossflow instabilities experience a destabilizing effect. Thus, there are considerable differences between the two approaches for modeling concentric roughness, regardless of whether the YHP results are compared against the MW solutions based on method A or B.
The significant differences between the MW and YHP models for concentric roughness are, at first, concerning. However, the reasoning for the variations in behavior quickly becomes apparent when comparing the respective undisturbed base flow profiles. Comparing the radial velocity profiles, F, shown here in Fig. 1(a) for the MW concentric roughness against the equivalent YHP solutions given in Garrett et al. 33 [see their Fig. 1(a)] reveals that while the peak F-value increases in both cases, the velocity at the disk wall is vastly different (due to the application of slip in the MW approach). In addition, the gradient of the azimuthal velocity profile, G 0 , near the disk surface, increases for the MW model [see Fig. 1(d)] but decreases for the YHP model as the respective roughness levels are increased [see Fig. 1(b) of Garrett et al. 33 ]. Hence, the two approaches for modeling concentric roughness establish very different base flow profiles, leading to very different stability outcomes.
The current investigation implements the parallel flow approximation, in which the radial dependence of the undisturbed base flow was neglected by setting the radius equal to the Reynolds number. Future studies should follow the direction of Thomas and Davies 18,23,37 and include non-parallel effects. In addition, the MW modeling approach implemented in the current study assumes that surface roughness can be achieved by replacing the no-slip boundary condition with slip conditions on the disk surface. Selecting different radial and azimuthal slip lengths establishes concentric, radial, and isotropic roughnesses of different heights. However, the relationship between slip and roughness height is not yet known, and it is necessary to carry out experiments that verify the results presented herein. Even more so, given the contrasting behavior between the results for the MW slip approach and that obtained for the YHP model. 33 Some attempt toward this aim has been made by the group at the University of Warwick. [49][50][51] A series of test disks were manufactured by cutting concentric grooves into an aluminum plate. Their experimental study demonstrated that concentric surface roughness can bring about some stabilization of the crossflow instability on a rotating disk. Given the observations of the Warwick group and the initial misinterpretation of how slip relates to roughness, 2,33-37 the results presented in Figs. 2 and 3 for concentric roughness suggests that the MW slip model may not be appropriate for modeling surface roughness, at least on the rotating disk.
A theoretical study that implements the correct (physical) no-slip boundary conditions over a rough wall and accounts for the surface geometry without a radial averaging of the flow would shed light on the appropriate boundary conditions to consider for the base flow and perturbations. This will be the focus of our subsequent investigations,  the aim of which is to determine the optimal roughness distribution to stabilize crossflow instability, leading to practical applications.

AUTHOR DECLARATIONS Conflict of Interest
The authors have no conflicts to disclose.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.  34 (in italics). With the exception of a few minor discrepancies for cases matched to radial roughness (i.e., g 6 ¼ 0 and k ¼ 0), excellent agreement is achieved in all cases.