Flow past a heated slippery corner

Discussed in this paper is the unsteady, laminar, two-dimensional, viscous flow around a corner subject to Navier-slip conditions characterized by a slip length, γ, with 0 ≤ γ ≤ 1. Comparisons with no-slip (γ = 0) boundary conditions are made to assess the impact of a slippery surface. Both isothermal and heated flows are considered including forced and mixed convection. A numerical solution procedure utilizing boundary-layer coordinates is outlined. Numerous simulations have been compiled for the Reynolds number Re = 100 and Rayleigh numbers Ra = 0, 1, 10. Various results and comparisons are presented and discussed.


Introduction
Although the no-slip condition has become the default boundary condition applied at a solidliquid interface, the validity of this condition was challenged at length during the nineteenth and early twentieth century.The notion of a slip condition was originally proposed by Navier [1].Cases in which the no-slip condition are known to fail include the flows of rarified gases or flows within microfluidic/nanofluidic devices [2] where the mean free path of the fluid approaches the characteristic length scale.For situations involving hydrophobic solids, a general Navier-slip condition yields results that are in better agreement with experimental observations than the traditional no-slip condition [3].Beavers and Joseph [4] proposed a semi-empirical slip condition for a liquid-porous medium interface which was later extended by Saffman [5].The surface roughness scale distinguishes hydrophobic surfaces from superhydrophobic surfaces.Hydrophobic surfaces can be characterized as having micron-sized protrusions, while superhydrophobic surfaces typically have nanometer-sized ridges or posts.The importance of superhydrophobic surfaces lies in the observation that these surfaces have been shown to produce a significant reduction in drag for both laminar and turbulent flows.Although the boundary condition experienced by a fluid over a superhydrophobic surface is no-slip, the air-fluid interfaces existing between the surface features are essentially shear free and hence responsible for reducing drag.The engineering of superhydrophobic surfaces has impacted various technologies ranging from microfluidic/nanofluidic devices to marine vessels.In this study we investigate the impact of slip conditions on corner flows.
Corner flow occurs when fluid is trapped between two intersecting planes.There are several different problems related to isothermal corner flow utilizing no-slip conditions.For example, for two-dimensional steady creeping flow near a corner there is the so-called paint-scraper problem solved by Taylor [6] where one plane moves parallel to itself with constant velocity, and there is also the hinged-plate problem solved by Moffatt [7] where the planes are rotating toward or away from each other about their line of intersection with constant angular velocity.Further, if the planes are stationary and the motion is driven by far-field stirring, Moffatt showed that there exists an infinite sequence of reversing eddies of decreasing size provided the corner angle is less than 73 • .In addition, there is also the three-dimensional problem of flow along a rightangle corner.This problem was orginally solved within the framework of laminar boundary-layer theory by Carrier [8] and later by Rubin and Grossman [9].Another well-known closely related problem is that of steady two-dimensional stagnation flow first studied by Hiemenz [10].Here, a uniform flow along the y axis impinges on a flat wall located at y = 0.The flow takes place in the upper x − y plane and the stagnation point, located at the origin, divides the flow into two streams along the wall which then depart in oppposite directions: one in the positive x direction and the other in the negative x direction.The case of unsteady flow was later investigated by Rott [11], and a class of similarity solutions pertaining to stagnation point flows having a far-field accelerating or decelerating incoming flow was discovered by Kolomenskiy and Moffatt [12].Lastly, D'Alessio [13] recently solved the corner flow problem by implenting boundary-layer coordinates with no-slip conditions.The focus of that study was on the formulation and use of boundary-layer coordinates with corner flow used as one illustrative example.
The corner flow problems subjected to Navier-slip conditions considered in this study include the isothermal, unsteady, two-dimensional flow around a right-angle corner, and the forced and mixed convection problems whereby the walls are maintained at a variable temperature that exceeds that of the surrounding fluid.In both cases we will solve the fully nonlinear Navier-Stokes equations instead of the simplified creeping flow or boundary-layer equations.Past studies regarding flow and heat transfer from a corner including free convection from heated walls and plates (vertical and horizontal) are well documented in [14], [15].These invesigations focus on aspects such as the corner angle, heating, and turbulence, and all involve no-slip boundary conditions.Thus, the current work adds new results to supplement those already existing in the literature.Lastly, some related thermally-driven laminar flow problems in bounded cavities include the research reported in [16] - [20].
The main focus of the present study is two-fold: to extend the work of D'Alessio [13] to include Navier-slip conditions, and to present new results for the corresponding forced and mixed convection problems.The paper is structured as follows.In the next section we formulate the problems considered in terms of boundary-layer coordinates.These problems include isothermal and heated flows.Following that, in section 3 a numerical solution procedure utilizing finite differences is proposed to solve the governing equations.Numerical results are then presented and discussed in section 4, and detailed comparisons are made with the case of no-slip boundary conditions.Finally, a brief summary is provided in section 5.

Mathematical formulation
In this section we formulate the problems for isothermal flow (i.e.constant temperature) and heated flows spanning forced and mixed convection.We begin with isothermal flow.

Isothermal flow
The flow configuration past a 90 degree corner is shown in Figure 1, and has walls located at x = 0 and y = 0. Here, we consider the unsteady, two-dimensional, laminar flow of a viscous, incompressible, Newtonain fluid around the corner in the first quadrant.We will assume the flow remains laminar and two dimensional for all time.
In the above, u, v denote the velocities in the x, y directions, respectively, while P is the pressure.The parameter Re refers to the familiar Reynolds number defined by Re = UL/ν where U, L are velocity and length scales, respectively, and ν is the kinematic viscosity.Equations ( 1) -( 3) are to be solved for x, y > 0 using two different sets of boundary conditions.One set of boundary conditions considered are the usual impermeability and no-slip conditions given by u = v = 0 along x = 0 and y = 0 .
The corresponding far-field conditions can be obtained by considering the equation which represents an approximation to equation (1) for large x with the understanding that v ≈ 0. The solution satisfying u = 0 at y = 0 and u → 1 as y → ∞ is given by [21] Applying the same argument for large y we obtain the similar result The second set of boundary conditions adopted are the impermeability and Navier decrease to zero, and hence become no-slip.As expected, γ = 0 recovers the no-slip condition.For a porous material, the dimensional slip length, denoted by γ * , is related to the permeability Π through the relation γ * = √ Π/Λ where Λ is the Beavers-Joseph dimensionless constant which takes on values in the range 0.1 -0.4 [4]; for common materials √ Π can be as large as 10 −2 mm [22].The far-field condition along the wall y = 0 is also expected to satisfy In this case the solution satisfying u = γ∂u/∂y at y = 0 and u → 1 as y → ∞ is given by [21] with the similar expression and plotted in Figure 2 are several far-field velocity profiles, u γ (y), for various slip lengths, γ.
We see that as γ → 0 the no-slip profile is approached.Since the flow is two dimensional it is beneficial to work in terms of a stream function (ψ)vorticity (ω) formulation.This is accomplished by defining which automatically satisfies the continuity equation, and System (1) -( 3) then simplifies to As noted in D'Alessio [13], the parameter λ is a measure of the boundary-layer thickness formed along the walls, and the far-field conditions suggest introducing boundary-layer coordinates ξ and η defined by An interpretation of working in terms of boundary-layer coordinates ξ, η is that the physical, or Cartesian, coordinates x, y become moving coordinates; that is, lines of constant ξ and η expand in time when plotted in a Cartesian coordinate system.This is ideal from a numerical point of view since the grid lines are alive and allowed to expand with the growing boundary layers to ensure adequate resolution with the passage of time.This is particularly important in the early development of the flow when the boundary layers grow more rapidly.Thus, in terms of the boundary-layer coordinates, and rescaling the stream function and vorticity according to equations ( 4) -( 5) transform to In terms of Ψ and Ω and the boundary-layer coordinates ξ and η the impermeability and no-slip conditions now take the form Ψ = 0 , ∂ψ ∂η = 0 along η = 0 and Ψ = 0 , ∂ψ ∂ξ = 0 along ξ = 0 .
Likewise, the far-field conditions become The far-field conditions are more complicated.Here, we only list those along the wall η = 0 with the understanding that similar expressions hold along the wall ξ = 0 The initial conditions for Ψ, Ω, at t = 0, denoted by Ψ 0 and Ω 0 , respectively, can be obtained by setting t = 0 in ( 7) -( 8), and thus satisfy System ( 9) -( 10) must be solved numerically as well as system ( 7) -( 8), and this will be discussed shortly.We note that for isothermal flow the problem is controlled by the parameters Re and γ with γ = 0 denoting the no-slip case.

Heated flow
We next consider the unsteady, laminar, two-dimensional forced and mixed convective flows in the first quadrant.In addition to the flow around the corner the fluid is also heated by the walls x = 0 and y = 0 with gravity acting in the vertical direction.The temperature along the walls is allowed to vary with position but not with time.The equations governing the motion and heat transfer process of a viscous, incompressible fluid are the continuity, Navier-Stokes and energy equations.The fluid is characterized by the following properties: the kinematic viscosity, ν, the thermal diffusivity, κ, the thermal expansion coefficient, α, and the thermal conductivity, k.While these fluid properties are assumed to remain constant, the fluid density, ρ, is allowed to vary with temperature in the usual linear fashion given by where T denotes the temperature while T 0 is the far-field constant ambient temperature, and ρ 0 is the corresponding density.
Expressed in terms of a stream function and vorticity as previously defined, the governing equations take the form Here, the Boussinesq approximation was used to describe the buoyancy force, and viscous dissipation was neglected.The parameters Ra and P r are known as the Rayleigh and Prandtl numbers, respectively, and are defined by where T w is the maximum wall temperature, L is the length scale, U is the velocity scale, and g is the acceleration due to gravity.The dimensionless temperature, φ, is related to the dimensional temperature, T , through the relation φ = (T − T 0 )/(T w − T 0 ).We note that the special case Ra = 0 corresponds to forced convection, and here the stream function and vorticity equations become decoupled from the energy equation.
As in the previous section, the far-field condition for φ along the wall y = 0 is expected to and the solution satisfying φ = φ ∞ at y = 0 and φ → 0 as y → ∞ is given by where φ ∞ is the far-field wall temperature.A similar expression holds along the wall x = 0.In this investigation we will consider wall temperature profiles that monotonically decrease with distance from the corner with φ ∞ = 0.In particular, in our simulations we have used the following wall temperature profiles φ(x = 0, y, t) = e −y 2 and φ(x, y = 0, t) = e −x 2 , which focusses the heating near the corner.As suggested by the above far-field condition, we switch to boundary-layer coordinates ξ, η.Further, we rescale the stream function and vorticity using Ω = λω , Ψ = ψ λ .
Equations ( 11) -( 13) then become The initial and boundary conditions for Ψ and Ω are the same as the isothermal case.The boundary conditions for φ, on the other hand, are as follows We notice that in terms of boundary-layer coordinates the wall temperature now varies with both position and time.The initial solution for φ at t = 0, denoted by φ 0 , can be obtained by solving 1 P r Since at t = 0 we have that λ = 0, it follows that the temperature along the walls will initially be constant with φ 0 (ξ = 0, η) = φ 0 (ξ, η = 0) = 1.Thus, the far-field conditions become Lastly, for heated flow the problem is completely characterized by the dimensionless parameters Re, Ra, P r and γ.

Numerical solution procedure
Numerical solutions were obtained using the finite difference method [23].The governing equations were solved on the rectangular domain [0, ξ ∞ ] × [0, η ∞ ] where ξ ∞ and η ∞ denote the outer boundaries approximating infinity.The computational domain was discretized into I equally spaced subintervals in the ξ direction and J equally spaced subintervals in the η direction forming a network of (I −1)×(J −1) interior grid points (ξ i , η j ) where ) and Δη = η ∞ /(J + 1) denoting the uniform grid spacing in the ξ, η directions, respectively.We begin by discussing the solution of the linear coupled system ( 9) -( 10) for the initial solutions Ψ 0 and Ω 0 .Replacing all derivatives by second-order, central-difference approximations leads to the discretized linear system of equations Here, we have adopted the notation f k i,j = f (ξ i , η j , t k ).The above system was solved using the Gauss-Seidel iterative procedure.The only complication that arises stems from the fact that the surface vorticity is not known.We can overcome this difficulty by deriving a second-order accurate formula for the surface vorticity as follows.
We present the details for the boundary η = 0 for the impermeability and Navier-slip conditions.Recall that along η = 0 we have and from (7) it also follows that Now, expanding Ψ in a Taylor series about η = 0 yields Using the known boundary conditions for Ψ together with the expressions for the higher derivatives yields the following expression for the surface vorticity In arriving at this result we have made use of the second-order finite difference approximation For the heated flow problem the initial energy equation ( 17) in discretized form is similar to (19) and becomes and can be solved in a similar manner.We next discuss the solution of systems ( 7) -( 8) and ( 14) -( 16).We first point out that the stream function equations ( 7) and ( 14) in discretized form are essentially the same as (18).Further, equations ( 8), ( 15) and ( 16) are similar in form, and thus, we will discuss their solution procedure simulataneously by casting them in the generic form where χ denotes either Ω or φ, and the function R(ξ, η, t) refers to all the remaining terms when brought to the right-hand side.Assuming the solution at time t is known, we can advance the solution to time t + Δt by integrating (22) and using integration by parts to obtain where Δt is the time increment.We now approximate the integrals using where β is a weight factor and F is a generic function.In general, 0 ≤ β ≤ 1 and we have used β = 1/2 which yields the well-known Crank-Nicolson scheme.With this approximation in place we obtain Upon substituting the expression for R(ξ, η, t + Δt) and replacing all spatial derivatives using second-order, central-difference approximations, equation ( 23) becomes a nonlinear system of algebraic equations which is to be solved in conjunction with (20) and the discretized stream function equation using the Gauss-Seidel iterative procedure to determine Ψ, Ω and φ (for the heated case) at time t + Δt at all the grid points.The convergence criterion adopted is that the maximum difference between successive iterates must be less than a specified tolerance .

Results and discussion
Numerous simulations were conducted to determine the optimal values of the computational parameters.For computing the initial solution at t = 0 the locations of the outer boundaries were taken to be ξ ∞ = η ∞ = 5 and the computational grid (I − 1) × (J − 1) = 100 × 100 was judged to be sufficient.This resulted in equal grid spacings of Δξ = Δη = 0.05.Since the unsteady computations were carried out to t = 25 the grid and outer boundaries were doubled in order to ensure that the flow field was adequately captured for those simulations.In all our computations the time step Δt = 0.005 was used and the convergence tolerance was = 10 −6 .Lastly, no convergence problems were encountered in our numerical scheme with the parameters listed.

Isothermal flow
We first present the initial solution which is independent of the Reynolds number.This solution was used as an initial condition for obtaining unsteady solutions for a given Reynolds number.Plotted in Figure 3 is the surface vorticity distribution along the walls ξ = 0 and η = 0 using impermeability and no-slip boundary conditions.The distributions along the walls are identical suggesting that the flow is initially symmetrical.This is confirmed in the streamline plot shown in Figure 4 which is also plotted in boundary-layer coordinates.We notice a region of clockwise circulating flow near the corner, and based on the magnitude of the values of the streamlines the flow is nearly stagnant as it is rotating slowly.
In stark contrast the initial surface vorticity distribution using impermeability and Navierslip (referred to as just slip henceforth) boundary conditions is zero according to the boundary conditions, and the corressponding streamline plot displays no eddies as illustrated in Figure 5.For comparison purposes we also present the streamline plot predicted by creeping flow in Figure 6.Apart from the size of the eddies, we observe that Figures 4 and 6 are similar.The creeping flow solution was obtained by solving the biharmonic equation ∇ 4 ψ = 0.As explained in [7], a solution is sought in polar coordinates having the form where λ n are complex-valued eigenvalues (not to be confused with the parameter λ in ( 6)).Near the corner the first term in the series will be the dominant term.Thus, the approximate solution is ψ ∼ A 1 r λ 1 f λ 1 (θ) where with λ 1 = 3.7396+1.1191i,and plotted in Figure 6 is the real part of ψ in Cartesian coordinates.To ensure that the flow remains laminar and two dimensional for all time we limit our results to Re = 100, although solutions for larger Re can easily be obtained.Simulations were carried out for 0 < t ≤ 25.As previously noted, with the passage of time a larger computational domain was needed to fully capture the growing eddies.Shown in Figures 7 and 8 are surface vorticity distributions for Re = 100 at t = 25 in Cartesian coordinates along the walls x = 0 and y = 0, respectively, for various values of γ.We point out that for Re = 100 and t = 25 the parameter λ = 1, and hence, ξ = x, η = y.That is, the boundary-layer coordinates ξ, η are exactly equal to the physical coordinates x, y at that instant in time.These figures clearly illustrate that the vorticity distributions along the walls are no longer the same.We observe that the variation in surface vorticity decreases as γ increases as a result of slippage.Also, as γ → 0 the profile approaches that corresponding to no-slip as expected.Streamline plots are shown in Figures 9 and 10 for γ = 0 and γ = 1, respectively, and they reveal that the eddy size increases significantly with the passage of time, and that the eddy size decreases as γ increases.Also evident in these diagrams is the asymmetry in the eddies, that is, they extend further along the y axis than along the x axis.

Heated flow
In these simulations the Reynolds and Prandtl numbers were fixed at Re = 100 and P r = 0.7 and simulations were carried out to t = 25.We begin with the initial isotherm plot in boundarylayer coordinates shown in Figure 11 which depends only on P r and represents the numerical solution of equation (17).The isotherms appear to be concentric curves which emulate a rounded corner.The corresponding initial surface vorticity distribution and initial streamline plot are identical to Figures 3 and 4, respectively.In order to demonstrate the impact of heating we first present the results for forced convection with Ra = 0. Here, the streamline and surface vorticity distributions are the same as those presented for the isothermal case (see .Contrasted in Figure 12 in Cartesian coordinates are isotherms for γ = 0 (no-slip) and γ = 1.We immediately observe that the slip condition inhibits the diffusion of heat.In both   cases the temperature variations are confined to the corner region where the heating takes place, and heat is transported by the eddies formed there.Since the eddy is larger for the no-slip case, this explains the greater extent of the isotherms.
We next present the local surface heat transfer coefficients, Nu, defined by Nu(x, y = 0, t) = −2 ∂φ ∂y , Nu(x = 0, y, t) = −2 ∂φ ∂x , also known as the Nusselt number.Illustrated in Figure 13 in Cartesian coordinates is the Nusselt number variation along the walls for both no-slip and slip cases.Although the heat transfer coefficient along the bottom wall remains positive, the side wall shows a reversal in sign where it starts off positive near the corner and becomes negative further away.This behaviour occurs for both the no-slip and slip cases and is also the result of the eddies.Since the corner is the hottest region and the eddies circulate in the clockwise direction, warmer fluid is forced up along the side wall and when it becomes warmer than the wall it causes the sign change in Nu.
The situation along the bottom is unaffected mainly because the eddy extends almost twice the distance along the side wall than it does along the bottom wall.Hence, there is little effect on the temperature gradient along the bottom wall.Also note that Nu = 0 at the origin because of our choice in the imposed temperature profile along the walls y = 0 and x = 0.The impact of the slip conditions becomes apparent when the average Nusselt number, denoted by Nu, is computed.This was done over the intervals 0 ≤ x ≤ 7 and 0 ≤ y ≤ 7 since beyond that the Nusselt number is essential zero.Along the bottom wall we find that the no-slip value is Nu = 0.36 while for slip we obtain Nu = 0.49.Thus, along the bottom wall the slip conditions enhances heat transfer from the wall.For the side wall we find that the no-slip value is Nu = −0.14whereas the slip value is Nu = −0.28.This indicates that the slip condition enhances heat transfer to the side wall.So, overall the slip condition enhances heat transfer either to or from the wall.We now explore mixed convective heat transfer for Rayleigh numbers in the range 0 < Ra ≤ 10.Here, the momentum and energy equations become coupled through the buoyancy term given by Ra∂φ/∂ξ in (15).This acts as a source term for the generation of vorticity.To illustrate the effects of buoyancy we compare the surface vorticity and Nusselt number distributions for Ra = 0, 1, 10. Figures 14 and 15 contrast the surface vorticity and Nusslet number distributions, respectively, for the no-slip case whereas Figures 16 and 17 contrast the slip case with γ = 1.We see that as Ra increases so do the amplitudes in the fluctuations of the surface vorticity.Also, for a fixed value of γ, as Ra increases so does the absolute value of the rate of heat transfer.This is confirmed in Table 1 which lists the average Nusselt number for various cases.Another trend is that for a fixed Ra, as γ increases so does the absolute value of the rate of heat transfer.
Lastly, we compare streamline and isotherm plots for Ra = 10 in Figures 18 and 19, respectively.For the no-slip case the streamline plot reveals the formation of secondary counter-

Summary
This study quantified the impact of slip conditions in a 90 degree corner by contrasting noslip and slip simulations for both isothermal and heated flows.For non-isothermal flows the heat was concentrated near the corner by imposing a specified temperature profile along the walls that decayed rapidly with distance from the corner.To guarantee that the flow remained two dimensional and laminar, computations were performed with Re = 100, and the Rayleigh number was limited to the range 0 ≤ Ra ≤ 10.The slip length was allowed to vary from γ = 0 (no-slip) to γ = 1.In all our simulations the Prandtl number was set to P r = 0.7, and the unsteady calculations were carried out to t = 25.The key findings are as follows.For isothermal flow the slip condition reduced the variation in surface vorticity as well as the size  of the eddy formed near the corner when compared to the no-slip case.For heated flows both forced and mixed convective cases were investigated.With forced convection the slip condition enhanced heat transfer to and from the walls but reduced the overall extent of the transfer of heat.Finally, for mixed convection in addition to the trends observed with forced convection, the slip condition also suppressed the formation of secondary vortices that occurred with the corresponding no-slip cases.The values of Ψ are: 0.01, 0.025, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9.