Fully Developed Laminar Fluid Flow in a Rotating Isothermal Isosceles Triangular Channel *

A numerical study is conducted on laminar ̄ow of an incompressible viscous ̄uid through a triangular channel in the hydrodynamically fully developed ̄ow region. The triangular channel is subjected to a span-wise rotation, and secondary ̄ow motions are introduced by the Coriolis forces. A pair of counter-rotating, longitudinal vortices appears at low rotation speed (small ReReV). For moderate to large rotation speed, the vortices split into two pairs of counterrotating vortices. In certain regions, the counter-rotating vortices become highly unstable, oscillating between two solutions, i.e., one pair vortices and two pairs vortices. Isosceles triangular channels with three aspect ratios (g5 0.287, 0.5, 0.866) have been studied. A large range of ReReV and Ro (Rossby) numbers was calculated. The critical Reynolds number with respect to Rossby number, where the ̄ow instability occurs, is calculated and the dual solution region presented. The results obtained from the computations cover a broad range of parameters, especially from low rotational speed to high rotational speed. The friction factor, velocity pro®le, and streamline patterns are presented. A comparison of the numerical results with the available theoretical data is also presented.

Fluid ¯ow in rotating channels is encountered in many engineering applications.Examples include rotors of electrical machinery, gas turbines, and many other rotating machinery.Research on heat transfer in rotating pipes has intensi®ed with the aim of raising the operating 1 temperature of gas turbines.Many investigators have studied ¯ow and heat transfer characteristics in rotating ducts.A detailed literature review for rotating ducts is given by Jen (1987), and Hwang and Jen (1990).Some relevant papers are reviewed here.During 1950±60, a large number of theoretical investigations were done for ¯ow in channels with small rotation speeds.For example, Baura (1954), by using a perturbation method, obtained an approximate series solution in a rotating pipe, and showed that two counter-rotating vortices of equal size appeared.The friction factor was then computed from the series solution.These results were valid only in the case of laminar ¯ow with a small angular speed.Hart (1971) employed linear stability analysis for the laminar ¯ow region, and the result was in qualitative agreement with the experiment.He found that there existed a weak double-vortex secondary ¯ow at small rotation speed, and a secondary instability existed in the form of longitudinal rolls of non-dimensional wave number 5 at moderately high rotation speed.Hart also showed the analogy between the ¯ow instability and thermal instability of a temperature-strati®ed ¯uid.This phenomenon is also observed in the present study.Lezius and Johnston (1976) studied the Taylor-type roll cell instabilities for laminar and turbulent ¯ow in a rotating channel.They used linear stability analysis to predict the critical Reynolds number with respect to Rossby number where the instability occurs.They found the critical value for laminar plane.Poiseuille ¯ow occurred at Reynolds number of 88.53 and a Rossby number of 0.5.At low range of Rossby number (i.e.0.01±0.5), the critical Reynolds number decreased smoothly.Once it passed through the critical Reynolds number 88.53, the critical Reynolds number increased extremely rapid as Rossby number increased.They also showed that the stable limit for all modes is about Rossby number equal to 3.0.Trefethen (1957) obtained the friction factor of a circular duct experimentally in both laminar and turbulent ¯ow.Mori and Nakayama (1968) studied laminar ¯ow in rotating circular ducts by assuming velocity boundary layers along the pipe wall.The friction factor was obtained for large parameter values.Subsequently, the turbulent ¯ow in a rotating radial pipe was analyzed by Mori et al. (1971) with the same technique.Ito and Nanbu (1971) studied extensively the friction factor of a fully developed ¯ow in smooth wall straight pipes of circular cross section rotating at a constant angular speed about an axis perpendicular to its own for Reynolds number ranging from 20 to 60,000.Empirical equations for friction factors of small values of Re /Re were presented for both laminar and turbulent ¯ow.Skiadaressis and Spalding (1977) predicted the ¯ow and heat transfer characteristics for turbulent steady ¯ow in rectangular ducts.Cliord et al. (1983Cliord et al. ( , 1984) ) studied the low range of turbulent convective heat transfer in an isosceles triangular channel in the orthogonal mode experimentally.The mean and local Nusselt number variations with respect to rotational Reynolds number and rotational Rayleigh number were reported.They also showed evidence that the rotational buoyancy eect was important when the Reynolds number was less than 30,000.Speziale (1982Speziale ( , 1986) ) and Speziale and Thangan (1983) presented extensive studies of Eckaus instabilities of laminar ¯ow in rectangular ducts.For constant axial pressure gradient, they were able to study the qualitative nature of the ¯ow patterns.Hwang and Jen (1987) studied the laminar ¯ow and heat transfer problems in a square duct by using an upwind dierence scheme.The friction factor and Nusselt number ratios were plotted with respect to ReRe and PrReRe .They also showed that the eect of axial conduction should be considered for ¯uid with small Peclet number.Hwang and Jen (1990), by using the power law scheme (Patankar, 1980), extended the analysis to rectangular ducts of various aspect ratios.Hydrodynamical and thermal instabilities were also found in their study, and the discontinuity (i.e., jump) of the friction factor and Nusselt number due to the hydrodynamical and thermal instabilities were also presented.Very recently, Jen et al. (1992) and Jen and Lavine (1992) calculated the temperature and ¯ow ®elds in the entrance region of rotating isothermal rectangular ducts using a vorticity± velocity method.Their results show that the secondary ¯ow become unstable when a critical Reynolds±Rossby number combination is exceeded.In limiting cases, their results are in good agreement with the linear stability study presented by Lezius and Johnston (1976).Their results also compared favorably with experimental data.
This paper presents a theoretical analysis of forced laminar ¯ow in the hydrodynamically fully developed ¯ow region in a rotating isothermal isosceles triangular duct.The full non-linear Navier±Stokes equations are solved by the Power Law scheme (Patankar, 1980), and the ¯uid ¯ow for three dierent aspect ratios ( 0.287, 0.5 and 0.866) were analyzed.It is worth noting that most of the early studies concentrated on circular or rectangular channels.No previous ¯uid ¯ow calculations on triangular ducts have been published.In the gas turbines cooling system, triangular type channels can be found located in the trailing edge region of the blades (Cliord et al., 1984).Even though the aspect ratios of the triangles demonstrated here are slightly out of the typical range of turbine blade cooling channels, the present study still provide detailed physical inside to the ¯uid ¯ow distribution in rotating triangular channels.The results obtained from the computations cover a broad range of parameters, especially from low rotational speed to high rotational speed.The friction factor, velocity pro®le, and streamlines patterns are presented.A comparison of the numerical results with the available theoretical data is also presented.

THEORETICAL ANALYSIS
Consider a forced laminar ¯ow in an isothermal isosceles channel rotating at a constant angular speed about an axis normal to the channel longitudinal direction.The physical con®guration and coordinate system is shown in Figure 1.
The channel is long enough so that the end eect of the duct is negligible.A secondary ¯ow motion is introduced by the Coriolis force acting in a direction normal to the main ¯ow and the z-directions.A pair of vortices is formed by the interaction of the main ¯ow and the angular velocity.By introducing a coordinate transformation and assuming a steady, incompressible and hydrodynamically fully developed ¯ow, the continuity and momentum equations become Continuity equation The 2W term in the Y-momentum equation is the Coriolis force that drives the ¯ow to the negative Ydirection.The Z 2 term in the Z-momentum equation represents the centrifugal force that acts in the Z-direction.This centrifugal force is balanced by the hydrostatic pressure distribution in the Z-direction.The term 2V in the Eq.[5] is also the Coriolis force coming from the interaction between the Y-direction velocity and the angular velocity.Note that this Coriolis force acts on the Z-direction non-uniformly.
In order to obtain governing parameters relevant to the present problem, the following dimensionless transformations for the dependent and independent variables are introduced where D e is the hydraulic diameter, L is the heating length along the channel, U c is the characteristic velocity in both X-and Y-directions to be de®ned, W is the average axial velocity, P z (Z) is the hydrodynamically fully developed pressure distribution and 1 2 Z 2 2 is the hydrostatic pressure variation due to the centrifugal force, and P c is the characteristic pressure in cross-sectional directions to be de®ned.In this problem we consider case L/D e 1.

DIMENSIONAL ANALYSIS
By substituting the dimensionless transformation [6] into Eqs.
[2]±[5], and considering that the order of magnitude of the Coriolis force 2W is equal to that of the viscous terms and the pressure term in the Eq.[4], we ®nd 0W 0 @ 2 V @X 2 , 0W 0 @P 3 @Y : 7 We obtain the characteristic quantities U c and P c where Re WD e / is the Reynolds number and Re D e / is the rotational Reynolds number.De®ne the Rossby number Ro Re/2Re .The continuity and momentum equations become where C ÿ@P z =@Z 1 D 2 e =W is a constant determined by the global continuity condition i.e., w 1 (see Eq. [21]).
It is seen from the continuity equation (i.e., Eq. [10]) that the ¯ow is essentially two-dimensional given a sucient channel length.The parameter ReRe , indicates a combi- nation eect of the axial velocity W and the angular speed, , i.e., the combination eect of Coriolis force and viscous force.The Coriolis forces acting in the y-and z-directions are represented by the term ÿ2w in Eq.
[12] and ÿReRe v/Ro in Eq. [13], respectively.If the Rossby number is large (for approximately Ro > 10), the rotational speed is small, and the ÿReRe v/Ro term negligible.In the present study we consider the Rossby number from 1 (i.e., small rotation speed) to 0.05 (i.e.high rotation speed).
[14] and [15] are: @w @x 0 along the center line x 0, 16 w @ @n 0 on the channel wall: 17 Note that only symmetrical vortices will be obtained with the present boundary conditions [16], set for half of the channel.In the present study, it can be seen from Eqs. [13] and [14] that there are two parameters ReRe and Ro that determine the ¯ow characteristics of the rotating duct ¯ow system.Note also that other governing parameters can be chosen, such as Re and Re .However, the use of ReRe and Ro bears the advantage that the axial Coriolis force term, v/Ro, in Eq. [13] is negligible when the rotational speed is small (i.e., Ro ! 1).Thus, only one parameter, ReRe , is needed to describe the ¯ow behavior in a rotating system (Hwang and Jen, 1990).

FLUID FLOW CHARACTERISTICS
The ¯ow in a channel can be characterized by the friction factor.Following the conventional de®nition, the friction factor is written as: On the other hand, the friction factor can also be derived from the overall force balance, the result being f 1 Re C=2: 20 It was found from the computations that both formulations (i.e., Eqs. [19] and [20]) are in very good agreement (within 1% error).For simplicity, the latter formula is used throughout this study.

METHOD OF SOLUTION
The solutions for unknown variables u, v, w, , and in Eqs.
[13]±[15] with unknown constants C and satisfying boundary conditions [16] and [17] is a matter of considerable mathematical diculty.A numerical ®nite-dierence scheme is employed in the present paper to obtain the solution of the set Eqs.
[13]±[15], the power law ®nite-dierence approximation (Patankar, 1980) is used for the non-linear terms in these equations.The numerical procedure is (1) Assign initial values for unknown u, v, w, , and , and for parameters ReRe , and Ro.
(2) Give an initial guess for constant C and solved Eq.where n is the nth number of iteration.
(4) Then the vorticity transport Eq. [14] is solved for , and the associated boundary vorticity is obtained from Roache (1972), the formula is Grid function convergence tests are carried out for all dierent aspect ratios to ensure the grid independence in the numerical solutions.Table I shows fRe computed by using mesh sizes 40 2 40, 80 2 80, and 120 2 120 at ReRe equal to 0, 10 3 , 10 4 .It is seen that the dierence of fRe for grid sizes 80 2 80 and 120 2 120 are generally less than 1%.
Grid function convergence tests indicate that the numerical result becomes grid independent for relatively coarse grids, thus suggests that the error due to the Power law scheme (Patankar, 1980) becomes small enough with modest grid re®nement.The grid size of 80 2 80 was ®nally selected throughout the entire computation of this study.

RESULTS AND DISCUSSIONS
As shown in Figure 3(a), the ¯uid in the core region is driven by the Coriolis force in the negative Y-direction normal to the axis of rotation.The core region ¯uid pushes the bottom layer ¯uid to the side wall, and thus the ¯uid near the side wall ¯ows to the positive Ydirection, and a pair of counter-rotating eddies is  LAMINAR FLUID FLOW generated.As rotation speed increases, the Coriolis force becomes stronger and pushes the ¯uid more strongly to the negative Y-direction.This means that the secondary ¯ow intensity increases as rotation speed increases.The ¯ow pattern becomes unstable when the Coriolis force acting on a ¯uid particle is not balanced by the local pressure gradient (Lezius and Johnston, 1976).Another pair of counter-rotating vortices is generated as shown in Figure 3(b).This phenomenon will be described in more detail in later sections.
In order to illustrate the eects of the Coriolis force on ¯ow characteristics, the constant axial velocity contours (isovels), and streamlines in the isosceles triangular ducts are shown.Selected typical ¯ow patterns for low ReRe to high ReRe (10 4 , 10 5 , and 7.5 2 10 5 ) with aspect ratio 0.287, 0.5, 0.866 and Ro 1.0 (i.e. the axial Coriolis force term is not negligible) are presented in Figures 4(a)± (c).For the case with 0.886 (i.e., apex angle equals to 60 ), it can be seen from Figure 4(a) (i) that the isovels in the core region are pushed towards the negative Ydirection as ReRe increases.Furthermore, the location of maximum axial velocity is also moved downward by the strong secondary ¯ow in the center core region.This means that a large axial velocity gradient is formed near the bottom horizontal wall of the triangular duct, and the core region constant axial velocity lines become more ¯at.A high shear layer, called the Ekman layer, is formed near the side wall of the duct.It is also seen that the Eckhaus instability (Guo and Finlay, 1991) exists when ReRe equals to 10 5 (Figure 4(a) (ii)).The constant axial velocity lines near the bisector of the bottom horizontal wall are pushed back to the positive Y-direction and a closed high axial velocity region is found.The onset of the Eckhaus instability in fact increase the gradient near the bottom wall, thus the friction factor is expected to increase.At larger value of ReRe (i.e., ReRe 7.5 2 10 5 ), the Eckhaus instability disappears, and the isovels become denser near the bottom wall (Figure 4(a) (iii)).This is because the rotational speed has increased substantially.
The streamlines for ReRe 10 4 , 10 5 , and 7.5 2 10 5 are also shown in the ®gures.The center of the vortex moves toward the side and the bottom walls as ReRe increases from 10 4 to 10 5 .For ReRe 10 5 , the second pair of counter-rotating vortices exists near the bisector of the bottom horizontal wall.Considering the characteristic quantities in Eq. [8], and the streamlines in Figure 4(a), we know that the intensity of the secondary ¯ow is proportional to (ReRe ) n , where n < 1.From Eq. [15] the distance between two nearby streamline is inversely proportional to intensity of the secondary ¯ow.An exponent of n 0.5 is observed among the ¯ow patterns of ReRe 10 4 , 10 5 , and 7.5 2 10 5 .It is interesting to note that the parameter (ReRe ) 0.5 is analogous to the Dean number (K ) for the curved pipe ¯ow (Dennis and Ng, 1982).
The cases with 0.5 (90 apex angle) and 0. on half domain of the duct.That means that for the small apex angle isosceles, there exist two local peak axial velocities at each side of the duct.This is due to the eect of Ekman layer, and the axial direction Coriolis force.Furthermore, if we compare the height of the unstable regions with respect to the height of the duct, it is interesting to ®nd that for all three ducts, the ratios are, 0.343 for 0.886, 0.349 for 0.5, and 0.347 for 0.287.
To see clearly the onset of the Eckhaus instability, the isovels and streamlines distributions for aspect ratio of 0.287 (i.e. 120 apex angle) with Ro 1 are shown in Figure 5.It is observed that the instability appears between ReRe from 15,000 to 18,000.This means that at a critical Reynold number and Rossby number the breakdown of one pair of eddies to two pairs of counter-rotating eddies occurs.The streamlines distribution shown in the ®gure demonstrates that at ReRe 15,000, the secondary ¯ow near the bisector of the bottom wall is almost stagnant, when ReRe is increased to 18,000, the large vortex breaks down.A second pair of counter-rotating vortices appears at the center of the bottom wall of the duct.
For the axial velocity distribution in the centerline, Figure 6(a) depicts the dimensionless axial velocity distribution in an isosceles triangular duct with 120 apex angle (i.e., 0.287) for various values of the parameter ReRe at Ro 1.As ReRe goes from 0 to 10 4 before the instability appears, the curves are pushed toward the negative Ydirection by the secondary ¯ow.This causes a large axial velocity gradient near the wall as ReRe increases.Thus, a larger friction loss should be expected near the bottom wall of the triangular duct.Once the vortex breakdown occurs, the centerline axial velocity near the bottom wall are pushed back toward the positive Y-direction due to the generation of the second pair vortices, as indicated by the curves for ReRe 2 2 10 4 and 5 2 10 4 .At ReRe 10 5 , the Eckhaus instability disappears, and the axial velocity pushes further to the bottom wall.It is also observed that the centerline axial velocity distribution is almost linear at ReRe < 10 4 .
Figure 6(b) depicts the eect of axial velocity across the duct at y 0.25b (where b is the height of the isosceles duct), for 0.287, Ro 1.It is observed that the axial velocity gradient at the slope wall increases as ReRe increases.The Coriolis force acting on the z-direction, i.e., v/Ro term in Eq. [13], pushes the center core axial ¯ow upstream, and downstream near the slope wall, depending on the sign of the v component of the secondary ¯ow.Thus, a uniform axial ¯ow distribution at the core region can be observed in the ®gure.It can also be seen that the core region axial velocity becomes smaller as ReRe increases.
This reveals that the upstream Coriolis force becomes stronger for the case of increasing ReRe .For this case, the Eckhaus instability occurs at the curve of ReRe 10 5 .A small ¯uctuation in axial velocity distribution near the core region is observed in the ®gure.It can also be seen clearly that there is a layer formed near the slope wall.This layer is called the Ekman layer.This layer has the thickness of the order of O(E 1/4 ), where E is the Ekman number de®ned below (Bennets and Hocking, 1974;Hocking, 1967): In the present study, the layer thickness is determined from the wall to the location of maximum velocity.It is observed that the order of magnitude of the thickness of the Ekman layer is in good agreement with this theoretical analysis.
In Eq. [13], the constant C is determined by satisfying the condition [21] (i.e., global continuity).The product of friction factor and Reynolds number can be easily obtained from Eq. [20].The curves for the ratio fRe/( f Re) 0 versus ReRe with various Rossby number and aspect ratios are shown in Figures 7(a)±(c).
A detailed discussion is presented for Figure 7(a) (i.e., 0.287) only since the behaviors in Figures 7(a)±(c) are quite similar.We can see from the ®gure that when the Rossby number is greater than 10 (i.e., small rotation speed), the two fRe/( fRe) 0 curves, corresponding to Ro 1 and Ro 10, are essentially the same.This means that the eect of the Rossby number is small in the large Rossby number region.Thus in Eq. [13], the Coriolis force induced by the secondary ¯ow term ÿReRe v/Ro is negligible when Rossby number is greater than 10.In this region, the only dierence is that they have dierent jump values (i.e., discontinuities) of friction factor ratio at the onset of the instabilities.This means that the onset of the LAMINAR FLUID FLOW Eckhaus instability does indeed increase the friction factor.This instability exists for large to moderate Rossby number, from 1 to 0.5 (not shown), disappears at about Ro 0.2 for all three aspect ratios with the parameter range shown.This tells us that there must exist a stable limit of the Eckhaus instability.This will be discussed later.
Another interesting result is that the jump value of the friction factor ratio where the instability occurs is analogous to the local height of the second pair vortices located at the center of the bottom wall.Note for the curve of Ro 1, the increase in friction factor ratio (with respect to fRe 0 ) is about 33% due to the onset of Eckhaus instability.This is in agreement with our earlier analysis that the second pairs of the cells are about one-third of the height of the duct.This curve shows that the instability disappears at about ReRe 75,000.This result is similar to that obtained by Cheng et al. (1976) for a curved pipe at high Dean number.It can be observed from the ®gure that the hysteresis loop appears in the range where Eckhaus instability exists.From the curve with Ro 1, the onset of the Eckhaus instability can be seen when the ReRe 16,000.Now, if we perform our numerical experiments starting from this value (with second pair of the cells) and calculate backward to smaller ReRe , we can see that the ¯ow instability still exists until ReRe 10 4 .A similar situation can be seen when the ¯ow instability disappears at ReRe 75,000.If we calculate backward using this value, it can be seen from the ®gure that the ¯ow instability does not re-occur until ReRe reaches about 40,000.This means that a dual solutions region is formed where both two and four vortices solution exist.Table II shows values of Ro 1 and 0.866 in a typical dual solution region.
Note that this is analogous to the results obtained by Dennis and Ng (1982) for a curved pipe for low to medium Dean number.The stable limit points for dierent aspect ratios of isosceles triangular ducts are obtained numerically, and plotted in Figure 8 alongside the theoretical results of Hart (1971) and Lezius and Johnston (1976), which were obtained from linear stability analysis.The numerical results from Jen and Lavine (1992) for rectangular ducts are also plotted in this ®gure for comparisons.Note that the de®nition of the Rossby number in Lezius and Johnston (1976) is D/w and is dierent from that de®ned here.The relation becomes Ro 0 1/(4Ro), where Ro 0 is the Rossby number de®ned in Lezius and Johnston (1976), and is the abscissa in this ®gure.There is also a dierence of a factor of two in the Reynolds number used in this study and that de®ned in Lezius and Johnston (1976) because the hydraulic diameter as de®ned here is two times larger.These results tend to be in qualitative agreement in trend with both Hart (1971) and Lezius and Johnston (1976).
Two curves with aspect ratios, 0.287 and 0.5, are presented in this study.It is seen from the ®gure that these  curves lie quite close to that analyzed by Lezius and Johnston (1976) when Rossby number (Ro) is large (or small Ro 0 ).These are reasonable results since as the apex angle approaches 180 degree the solution should tend to be the same as the rotating parallel plate ¯ow as analyzed in Lezius and Johnston (1976).At smaller Rossby number, the curves deviate from the linear stability curve.It is worth noting that our calculations are solved using the full Navier±Stokes equations including all the nonlinear terms in the equations, while Lezius and Johnston's results are from a Linear stability analysis.As shown in Speziale (1982) the nonlinear (i.e., convective) terms in the governing equations plays an important role when strong asymmetry occurs.Therefore, in such geometry, the full Navier± Stokes equations must be solved to ensure accurate results.It is also observed that the critical Re 0 number increases with aspect ratio.The Ro 0 number corresponding to the critical Re 0 number decreases as aspect ratio increases.The curves in the ®gure show that the slope of the curves after the critical point increases as the aspect ratio increases.The all mode stability limit for the Rossby number decreases as apex angle decreases.
Because of the lack of experimental data in the fully developed region of an isothermal isosceles triangular duct, the present numerical curves are compared with the existing rectangular data in the fully developed region (Hwang and Jen, 1990).Figure 9 shows the values of fRe/ ( fRe) 0 vs ReRe for Ro ! 1 with three dierent aspect ratios.It can be seen from the ®gure that the equilateral triangle ( 0.866) has the largest friction factor ratio, and the friction factor ratio decreases as aspect ratio decreases (or apex angle increases).For the equilateral triangular duct, it is almost coincident with the square duct case ( 1, & symbols) from Hwang and Jen (1990).For the case with 0.287 of the triangular duct, it seems that the friction factor ratio is much larger than the 0.2 case of the rectangular duct.From the data shown above, a brief conclusion can be drawn here that the triangular ducts seem to have larger friction loss than the rectangular ducts.The present analysis also agrees qualitatively with the ¯ow structure with the curved triangular channel presented by Nandakuma et al. (1993) (not shown).

CONCLUDING REMARKS
(1) The equations governing the fully developed laminar ¯ow of rotating isothermal isosceles triangular ducts have been solved using vorticity±stream function method for three-dimensional ¯uid ¯ow.(2) It is found that at moderate rotational speed, the Eckhaus instability occurs, and the second pair vortices generated near the bisector of the duct.(3) It is seen from ®gures that the unstable disturbed Coriolis force is acting on the negative Y-direction towards the bottom wall.If we compare with unstable four-vortex pattern for all three dierent aspect ratios, it is found that the ratios of the height of the unstable regions with respect to the height of the duct are almost the same.(4) The numerically obtained stability boundary points for two dierent aspect ratios (i.e., 0.287 and 0.5) are plotted with the linear stability results for parallel plates.It is found that at small aspect ratio (i.e., 0.287), the boundary points follow closely at small rotational speed, and deviates at larger rotational speed.
the mean wall shear stress.Here w can be derived from the averages of the local derivatives.The friction [13]   for w simultaneously by using a point SOR method.The value C is adjusted by considering the relation The equation for and in Eq. [15] is solved for by a point iterative successive under-relaxation (SUR) scheme until the following criterion if ful®lled: the kth number of computation from Steps 2 to 5. (7) Calculate the friction factor from the Eq.[20].

FIGURE 2
FIGURE 2 Vorticity at the slope wall.

TABLE I
FIGURE 3 (a) Two cells ¯ow pattern; (b) four cells ¯ow pattern.