Stability of the Plane Bingham–Poiseuille Flow in an Inclined Channel

We study the stability of laminar Bingham–Poiseuille flows in a sheet of fluid (open channel) down an incline with constant slope angle β∈(0,π/2). This problem has geophysical applications to the evolution of landslides. In this article, we apply to this problem recent results of Falsaperla et al. for laminar Couette and Poiseuille flows of Newtonian fluids in inclined channels. The stability of the basic motion of the generalised Navier–Stokes system for a Bingham fluid in a horizontal channel against linear perturbations has been studied. In this article, we study the flows of a Bingham fluid when the channel is oblique and we prove a stabilizing effect of the Bingham parameter B. We also study the stability of the linear system with an energy method (Lyapunov functions) and prove that the streamwise perturbations are always stable, while the spanwise perturbations are energy-stable if the Reynolds number Re is less than the critical Reynolds number Rc obtained solving a generalised Orr equation of a maximum variational problem.


Introduction
The stability of laminar flows in an inclined channel is important in many geophysical and industrial applications. Recently, Falsaperla et al. [1,2] described the stability of laminar flows of the steady solutions of a sheet of a Newtonian fluid (open channel) both in the horizontal case [1], and down an incline with constant slope angle β [2]. In a reference frame O, i, j, k with i parallel to the incline and k orthogonal to the plane of the channel, let z be the coordinate of the axis orthogonal to the incline, the basic steady motion is given by the velocity field U(z)i and the pressure p(z). For a laminar Newtonian fluid in a horizontal layer, U(z) is linear for Couette flows and parabolic for Poiseuille flows, while, for an inclined channel, U(z) is a parabolic function which vanishes at the bottom of the channel and whose derivative with respect to z vanish at the top.
Choosing an appropriately weighted L 2 -energy equivalent to the classical energy norm, we prove in [1] that the plane Couette and Poiseuille flows are nonlinearly stable with respect to streamwise perturbations for any Reynolds number (see also [3]). We also prove in [4] that the plane Couette and Poiseuille flows are nonlinearly stable (and linearly stable in the energy norm) if the Reynolds number is less then the critical value Re Orr / sin ϕ, when the perturbation is a tilted perturbation in the x direction, with x a coordinate associated to a new reference frame with i = i cos ϕ + j sin ϕ, j = j cos ϕ − i sin ϕ, and k = k. The number Re Orr is the Orr [5] critical Reynolds number for spanwise perturbations which, for the Couette flow, is Re Orr = 177.22 and for the Poiseuille flow is Re Orr = 175.31. The results in [1] improve those landslides. We also plan to investigate, as in [1], the stability with respect to oblique perturbations and its relation with the experiments. The layout of the paper is the following. In Section 2, following [8], we introduce the governing equations for a laminar open flow of a Bingham fluid down an inclined channel. In Section 3, following [11][12][13][14][15][16], we write the linearised equations for a perturbation to the basic motion U(z) with the appropriate boundary conditions, and we solve the generalised Sommerfeld equations for the eigenvalue problem with the Chebyshev collocation method. In Section 4 we study the stability of the linearised perturbations with the Lyapunov method by using an energy norm, we find sufficient conditions for stability in energy for the streamwise and the spanwise perturbations and, for the latter, we write the Euler-Lagrange equations for the maximum problem and we solve the related generalised Orr equation. In Section 5, we draw some conclusions.

Laminar Open Flow of a Bingham Fluid Down an Inclined Channel
Consider the stationary laminar motion of a mud layer of thickness D along an inclined plane of infinite length with constant slope angle β, 0 < β < π/2 [8] (p. 100) in a reference frame O, i, j, k with associated coordinates x, y, z. The motion takes place in a channel D 1 = R 2 x,y × [0, D] z of thickness D. The x-axis is taken along the slope direction while the z-axis is perpendicular to the slope. The y-axis is also parallel to the slope and orthogonal to the slope direction x. The channel extends indefinitely in the x, y directions and has a finite depth D in the z direction (see Figure 1). pl ug la ye r sh ea r la ye r so lid la ye r g x z D S D β Figure 1. Inclined layer of an angle β. The light grey layer is the shear layer which is moving with a parabolic law. The wall corresponding to z = 0 is the rigid terrain. The superior layer, the plug region D S ≤ z ≤ D, is rigidly moving at a constant velocity U(D S ). D is the depth of the layer, D S is the depth of the shear layer and D − D S is the depth of the plug region. g is the gravity.
We consider the flow of an incompressible Bingham fluid with yield stress τ 0 and plastic viscosity (Bingham viscosity) µ B . The wall placed in z = 0 is rigid (the velocity of the flow vanishes) while that placed at z = D is stress-free (the shear stress is zero). The governing equations are (cf. [11][12][13][14][15][16]): where U is the velocity field, ρ is the constant density,p is the pressure, τ is the deviatoric extra-stress tensor, ∇ is the gradient operator, and g is the gravity, g = g sin β i − g cos β k. The velocity vector is of the form U = Ui + Vj + Wk, where U, V, W are the velocity components, and i, j, k are unit vectors in the streamwise (x) and spanwise (y) directions, and normal to the wall (z), respectively. As in [11][12][13][14][15][16], we assume the constitutive relations between the stress and the strain rate for a Bingham fluid whereγ and τ are the strain rate and deviatoric stress tensors respectively. The yield stress is denoted by τ 0 and η denotes the effective viscosity, defined by where µ B is the Bingham viscosity (see De Blasio [8]), also called limiting viscosity (cf. [11]), τ 0 is the yield stress,γ and τ are respectively the second invariants of the strain rate and the deviatoric stress tensors. They are given byγ = γ ijγij /2, τ = τ ij τ ij /2,γ ij = U i,j + U j,i , and the Einstein convention on summation of repeated indices has been assumed. Following De Blasio [8], we assume that "the thickness of the mud flow does not change along the x coordinate, hence, U andp are independent of x". So, we consider laminar shear flows U = U(z)i,p =p(z) which are solutions to the stationary equations Solving these equations, we have    τ xz = ρg sin βz +Ā, p = −ρg cos βz +B, whereĀ andB are integration constants that can be determined by imposing the chosen boundary conditions, that, in our case, are vanishing shear stress at the upper surface of the mud, i.e., ∂ z U(D) = 0, and pressure equals to the atmospheric pressure p atm . Applying Bingham Equation (1) and denoting Following De Blasio [8] (p. 101), the first of two, Equation (3), can now be integrated from the base up to the level where the right-hand side becomes zero. This defines a layer known as the shear layer (the yielded zone) whose thickness D S can be found from Equation (3) imposing τ = τ 0 . We find

From Equation (3) and the definition of D S we have
We note that for z ≤ D S the velocity U(z) increases with height. For z ≥ D S , as De Blasio observes, the weight of the overlying material is not sufficient to shear the material, as it is obvious from the second of Equation (4). As a consequence, the material moves rigidly, U(z) = c 0 (c 0 a constant) like a plug. The thickness of this plug layer (unyielded zone) is Integrating (4) and adopting the boundary conditions We have the basic Bingham-Poiesuille flow,p = p atm + ρg cos β(D − z), and These equations are non-dimensionalized using a length scale D, a velocity scale U 0 = U(D S ) = (ρg sin βD 2 S )/(2µ B ), and a pressure-stress scale ρU 2 0 . Introducing the non-dimensional Reynolds and Bingham numbers where z S = D S /D is a positive number that represents the depth of the yield surface. Such depth can be computed directly, using the velocity U 0 and the Bingham number B, and solving the equation , using the momentum conservation equation in the plug region, where we have τ 0 = |τ w |, with τ w the wall shear stress |τ w | = τ(z = 1) in the yielded region. In non-dimensional form, we have We observe that in the limit when B tends to 0, the fluid becomes a Newtonian fluid, and z S tends to 1. In such a limit case there is no plug region.

Linear Stability Analysis
In order to perform the linear stability analysis, as in [11][12][13][14][15][16], we introduce ( u, p), infinitesimal perturbations to the basic flow (U,p), with << 1 and |u| and p bounded. The linearised perturbation wherever the yield stress is exceeded, the effective viscosity of the perturbed flow is expanded about the basic flow (see [14]). Denoting with µ the non dimensional effective viscosity, we have where A simple calculation shows that (7) where A depends on the number , (u, v, w) are the components of u, and the subscript indicate partial differentiation with respect to the variables. Expanding µ(U(z) + u(x, y, z, t)) in powers of , we have where U is the derivative of U with respect to z. Finally, at the first order in we have: We recall that Computingγ(U) and by observing that U ≥ 0, at the first order in , we have From this relation, we easily obtain with Collecting (5), (8) and (9), we have where we have used the convention 1 = x, 2 = y, 3 = z and e 1 = i, e 2 = j, e 3 = k. The perturbation Equation (10) in the open domain where ∆ is the Laplacian. We assume that the perturbations u, p are periodic in x and y of periods 2π/a, 2π/b, with wave number a and b in the directions x and y. We will later consider the limit cases of streamwise perturbations, i.e., perturbations which do not depend on x and are periodic only on y, and spanwise perturbations, i.e., perturbations which do not depend on y and are periodic only on x. Equation (11) must be completed with a set of boundary conditions. Since, for a yield stress fluid U = −2(z − z S )/z 2 S and 1/U is singular as z → z S , the effective viscosity goes to +∞ at the yield surface. It follows that some care is required in deriving the boundary conditions at z = z S as explained for example in [11,14]. As in [11], we assume that the yield surface z = z S is perturbed tô z = z S + h(x, y, t), with h(x, y, t) periodic in x and y of periods 2π/a, 2π/b. The boundary conditions can be obtained as in [11,[13][14][15][16] and turn out to be The boundary conditions at z = 0 are rigid, as classically found in Newtonian fluids. The Dirichlet boundary conditions u = 0 at the yield surface can be deduced from the fact that the unyielded plug zone is constrained to move as a rigid body according to the Bingham model [13]. The Neumann conditions come from linearisation of the conditionγ ij (U + u) = 0 at the perturbed yield surface, onto the unperturbed yield surface position [13]. In particular, following [13], we note that the conditions ∂ z v(x, y, z s , t) = ∂ z w(x, y, z s , t) = 0, are necessary due to the singular behaviour of 1/U on z = z S in the non-Newtonian part of (11). These boundary conditions ensure that (11) are well-defined as z → z S . We thus have that v and w are asymptotic at z S to (z − z S ) α , with α > 1. This fact will be useful in the Lyapunov investigation, in particular when using the boundary conditions in the periodicity cell Ω = [0, 2π/a] × [0, 2π/b] × [0, z S ].
In order to study the linear stability with the spectral method, we assume that the perturbations have the form (u, v, w, p, h) = [u(z), v(z), w(z), p(z), h]e i(ax+by)+Ct , where C is a complex number whose real part gives the decay rate of the perturbation. We observe that, in [13], Frigaard and Nouar prove a Squire-like (see [18]) theorem: "if the equivalent eigenvalue bounds for a Newtonian fluid yield a stability criterion, then the same stability criterion is valid for the Bingham fluid flow, but with reduced wavenumbers and Reynolds numbers". Assuming this, we consider in what follows only streamwise or only spanwise perturbations.
In the case of streamwise perturbations a = 0 and b = 0, we hence have An easy calculation proves that the streamwise perturbations are always stable (the real part of C is always less than 0). For the spanwise perturbations a = 0 and b = 0, and we have The variable v appears only in the second equation of system (12), and it can proven that v cannot destabilize the basic motion. The remaining three equations depend only on u, w and on the pressure p. Using the continuity equation, one can further eliminate the variables u, p and derive a system of one equation in the variable w. This can be obtained by subtracting ia times the z-derivative of the first equation of (12) and adding the third equation of (12) multiplied by a 2 (this amounts to writing the equation of the third component of the doube-curl of u). In this way one obtains the generalised Orr-Sommerfeld equation with boundary conditions w(x, y, 0) = w z (x, y, 0) = 0, w(x, y, z S ) = w z (x, y, z S ) = 0.
Numerical calculations show that the linear perturbations are always stable, i.e., the real part of C is negative for any Reynolds number. The Chebyshev collocation method has been used up to 100 polynomials for w. This transforms the differential equation in an algebraic generalized eigenvalue system which also takes into account the given boundary conditions. Figure 2 shows that the maximum real part of the spectrum of C is negative and decreasing with the Bingham number. Here, we have considered the significative ranges a ∈ [0, 2.6], Re ∈ [1000, 11000]. We also obtain the expected fact that as Re tends to 0, the real part of C tends to −∞.

Linear-Energy Stability of Perturbations
In this section we use the linear Lyapunov method, also typically called energy-stability method. In the periodicity cell Ω, we introduce the scalar product ( , ) and the L 2 (Ω) norm · . We have the energy identities: These identities can be simplified by applying integration by parts and using the boundary conditions. Since the term (Uu x , u) integrates to zero We have the simplified identities d dt d dt

Stability of Streamwise Perturbations
To prove that the linear streamwise perturbations are always stable, we begin observing that, for streamwise perturbations, ∂ x ≡ 0. The solenoidality of the kinetic perturbation hence reads v y + w z = 0. From (14) and (15), we have From (16), we have We observe that asz → z S we have Hence, the boundary terms vanish. Proceeding in the same way for the other terms, we easily obtain: Poincaré inequality allows us to prove that v and w decay exponentially to 0 as t → +∞. By using (13), it easy to see that also u goes to 0 as t → +∞. Therefore, the energy norm of the kinetic perturbation tends to 0 for any Reynolds and Bingham number. We conclude that the streamwise perturbations are always linearly stabilizing.

Stability of Spanwise Perturbations
To investigate the spanwise perturbations we begin observing that, in this case, ∂ y ≡ 0. It follows that the solenoidality of the kinetic perturbation reads u x + w z = 0. From (14), integrating by parts, and using the boundary conditions as before, we have: Applying the Poincaré inequality and observing that 1/U = −z 2 s /(2(z − z S )) > 0 for z < z S , we have that v goes exponentially to 0 as t → +∞. The energy equations of u and w become: Adding the two equations together we have d dt From this equation, we infer that the critical Reynolds number for spanwise perturbations is greater than the critical Reynolds number for Newtonian fluid (cf. [2]). Indeed, the term multiplied by the Bingham number is always non-positive, hence that term is always stabilizing.
The linear energy stability can be studied as in [1,2] using the classical energy In order to obtain sufficient stability condition, we use the Orr-Reynolds energy identity: whereV is the orbital time derivative (i.e., the Lagrangian derivative computed along the solutions of (11), where now ∂ y ≡ 0), and I = −(wU , u) and D = ∇u 2 + 2B We haveV and S is the space of the kinematically admissible fields In the definition of S, the space H 1 (Ω) is the Sobolev space of all functions which are in L 2 (Ω) together with their first generalized derivatives. The Euler-Lagrange equation of this maximum problem is given bȳ where λ is a Lagrange multiplier.
Since, for spanwise perturbations, v ≡ 0 and ∂ y ≡ 0, by taking the third component of the double-curl of (17) and by using the solenoidality condition u x + w z = 0, we obtain the generalised Orr equation with boundary conditions w = w z on z = 0 and z = z s . We solve this boundary value problem with the Chebyshev collocation method, and we obtain the critical energy-linear Orr Reynolds number, R c . In Figure 3 we have plotted the critical Orr-Reynolds number as a function of the Bingham number. The graphic shows the stabilizing effect of yield stress (Bingham). We note that in the limiting case as B → 0 we obtain the critical value for the Newtonian fluid (note that here the depth of the layer is D while in [2] it is 2D, therefore in the limiting case we obtain a critical value which is twice the size of the corresponding critical value of the Orr-Reynolds number obtained with the same boundary conditions).

Conclusions
In this article, we study the linear stability of laminar Bingham-Poiseuille flows for an inclined sheet of fluid (open channel) with constant slope angle β. The basic flow is parabolic in the shear region (the fluid there behaves like a Newtonian fluid) and is constant in the plug region. In the horizontal case, this problem has been investigated by many authors. In this work we follow the approach of Frigaard, Nouar and co-workers and we perform a novel investigation of the linear stability of the laminar solution which takes into account the inclination of the channel and a novel investigation of the critical value for energy stability, which possibly may suggest an experimental instability threshold.
The time evolution of linear perturbations of the basic solution is investigated both with spectral methods and with energy Lyapunov methods and assuming that the more destabilizing perturbations are two-dimensional (Squire-like Theorem). By using a Chebyshev collocation method, we prove that both the streamwise and the spanwise perturbations are stable for any Reynolds and Bingham numbers. We also show that the Bingham number has a stabilizing effect. This stability result, in the inclined channel, is similar to that of the Newtonian case [2] for the classical Couette motion in the horizontal layer.
The experiments show turbulence when the Reynolds numbers are sufficiently large (cf. also [19]), this turbulence cannot be accounted for by the spectral methods. The energy methods, instead, indicate the possible insurgence of instability for large enough Reynolds numbers (Couette paradox). In [1] we have already approached this question by introducing the energy method for a particular set of tilted perturbations, with tilt angles that appear in the experiments, and we have found a good agreement with the experimental data. We also apply to this problem the same ideas, and we obtain that the streamwise perturbations are always energy stable, while the spanwise perturbations are energy stable if the Reynolds number is less than a critical number which depends on the Bingham number. This suggests that there could be oblique perturbations that may physically destabilize the basic motion at the onset of instability.
The energy methods are also very important in the investigation of the nonlinear equations, since the nonlinear terms could be crucial in developing instabilities. For this reason we think that it is important, in a future work, to study the nonlinear stability problem and look for the most physical destabilizing perturbations and their possible relation with landslides.
Finally, we note that if we had introduced a Reynolds number independent of the inclination β, we would have obtained stability results with a different Reynolds number, inversely proportional to sin β. This means that small inclinations, as expected, are stabilizing.