Linear stability of viscoplastic flows down an incline

The stability analysis of viscoplastic flows down an inclined plane is done by comparing results obtained through theoretical and numerical studies of “regularized” models. The theoretical analysis is performed for Regularized Bingham and Casson-like fluids via the long-wave approximation method. In particular, the Bingham and the Casson flow have different stability characteristics, for Bingham-type materials an increase in yield stress leads to flow destabilization, while Casson-type materials show the opposite behaviour. The numerical study is performed by using the Papanastasiou and the “exact” Bingham model via a spectral method. The comparison between theoretical and numerical results shows excellent agreement. Our findings highlight that “regularized” and “exact” flow and can have stability characteristics, although they are “practically indistinguishable”. We validate our approach with the Regularized Bingham-like model, which is in rather satisfactory agreement with the experimental data.


Introduction
Materials with complex and different rheological behaviour from Newtonian fluids are usually involved in natural phenomena and industrial processes.Often their flow properties are characterized by a critical value of stress (yield stress, τ 0 ), i.e. above the yield stress the materials flow accordingly to their rheological properties, otherwise there is no material deformation [1,2].They are usually named as viscoplastic materials.The Bingham and Casson models are an example of viscoplatic materials.Properly describing fluids with complex rheological behaviour and the detection of critical conditions for the onset of instabilities may be useful to optimize industrial processes and to predict critical environmental situations.
In general, the critical Reynolds number, Re c , is the Reynolds number threshold above which flow instabilities can occur.
Since there is still an open debate regarding the yield stress existence [25,26,27,28,29], viscosity regularization methods have mainly been used in order to avoid the issues due to the inherent singularity in the study of viscoplastic flows."Regularized" constitutive equations are characterized by the introduction of a parameter, which is positive and provides the approximation level, and through which the "exact" model can be retrieved [30,34].
Our study focuses on the analysis of theoretical and numerical works regarding linear stability of viscoplastic flow down an inclined plane by using "regularized" models [32,33,34].In particular, the theoretical study has been performed through the long-wave approximation approach in the case of Regularized Bingham and Casson-like materials [32,33].The numerical study has been performed by using the Papanastasiou model (a Bingham law regularization) through a spectral collocation method that makes use of Chebyshev polynomials [34].
The theoretical background on the mathematical formulation of the problem and the main features of a "regularized" fluid flowing down an inclined plane in Section 2 and Section 3, respectively.Then, Section 4, similarly to [32,34,36], summarizes the analysis of linear stability through both analytical and numerical approaches, i.e. the long-wave approximation and spectral collocation method, respectively.Section 5 and 6, provides our findings and conclusions.

Mathematical Background and Problem Formulation
Here, we summarize the theoretical formulation of the problem, details are given in [32,33,34].Figure 1 shows the flow domain, where we denote the tilt angle as θ ∈ (0, π/2), the length of the domain as L and the upper free surface (not a priori known) as y = h(x, t) with H = max{h}.
We denote as T = −pI + τ the Cauchy stress tensor and as τ the deviatoric part.We recall that for the two-dimensional incompressible flow, v = ue x + ve y , we have in dimensional form where g = ρg sin θe x −ρg cos θe y , g is gravity and ρ is material density (constant).The boundary conditions consists in the non-slip and impermeability conditions u = v = 0 on y = 0 and, by denoting as n the outer normal (see Fig. 1), the kinematical-dynamical conditions is Figure 1: Reference framework.
We proceed by introducing the adimensionalization IC-MSQUARE-2023 Journal of Physics: Conference Series 2701 (2024) 012071 where γ = 1/2 ∇v + ∇v T is the strain-rate, and U denotes the characteristic velocity to be chosen in the sequel so that u(h) = 1.By inserting ( 5)-( 4) in (1) and by denoting as where the Reynolds number and Freude number are respectively, we have From ( 6) and (7) we have Now, we consider h = 1 and the flow as thus, we can rewrite (8) as and p = τ 12 = 0 on y * = 1, since τ * 12 is the only non-vanishing component of τ * .

"Regularized" Models
Here, we summarize the "regularized" constitutive laws in dimensionless form for a Bingham and Casson fluids through a regularization parameter, .The corresponding "exact" models and more details are in [32,33,34].We consider the chosen dimensionless regularization constitutive models where is the Bingham number where τ 0 is the yield stress, and * = 2 H/U , so that → 0 provides the corresponding "exact" model [32,33,34,35].In the sequel, the " * " are omitted to maintain the notation as light as possible.By using (9) the Bingham number can be rewritten as The parameter Γ depends only on the "material" and geometrical features, thus it is constant when the fluid and the tilt angle have been chosen.
From system (11) we get and, noting that 2| γ| = u y (y) and by using ( 12) Equations ( 17) provide the expressions of u(y), i.e. the basic profile u b (y), and by requiring that u(1) = 1, we get a bijective function of Re and ξ for each "regularized" model, namely meaning that i.e., for every positive Re , there exists ξ = ξ(Re) fulfilling ( 18) for given λ, , θ.

Linear Stability
We report the main features of the linear stability analysis and the details on the formulas computation and derivation here summarized are given in [32,33,34].The basic flow is given by h(x, t) = h b with h b = 1, v b = u b (y)i with u b (obtained through (17) and reported in [32,33,34]), and, p = p b (y) where p b (y) = ξ cot θ(1 − y) (see (15)).First, we consider a small perturbation of the basic flow, in the form of travelling waves, so that where we denote the wave number as α ∈ R, the complex wave speed as c ∈ R and the infinitesimal disturbance with the notation (•).Then, the velocity field can be rewritten as stream function, namely ψ(x, y, t) = φ(y)e iα(x−ct) , where, here and in the sequel, (•) represents the differentiation with respect to y.
Recalling that the growth/attenuation factor of the α th mode is the imaginary part of c, Im (c), the basic flow h b , v b , p b (y) is unstable if Im (c) > 0 for selected (Re, λ, and θ).The transition between the unstable (Im (c) > 0) and stable (Im (c) < 0) regimes is called marginal or neutral curve, i.e. when Im (c) = 0 for selected Re, λ, and θ.
By using ( 19)-( 21), the corresponding perturbed problem reduces to which is a boundary value problem and for each "regularized" model given by (16).It is worth noting that if B = 0 (Newtonian case) we obtain q = 2, s = 1 and from equation ( 22) 1 we retrieve the Orr-Sommerfeld equation [32,33,34].Therefore, non-trivial solution to the problem (22) (the eigenvalues c ∈ C and the corresponding eigenfunctions) allow one to determine Re c so that Im (c) = 0.In the next Subsections the main features of the analytical and numerical approaches used to solve (22) are briefly reported (details are in [32,33,34]).
Hence, the α th mode is stable when Re < Re c , while if Re > Re c instability occurs.Indeed, Im (c 1 ) < 0 when Re < Re c and vice versa.We compute the value of Re c through the solution of the system given by ( 18) and ( 27) by using MATLAB ® 2022a by means of function fsolve.

Numerical approach: Spectral method
We briefly report the developed numerical procedure for solving the boundary value problem (22) based spectral collocation method by means of Chebyshev polynomials of the first kind, so that the boundary value problem (22) can be expressed through discretized differential operators.Indeed, Chebyshev polynomials are appropriate to describe boundary layer phenomena which are commonly encountered in flow stability problems.By mapping the physical domain [0, 1] into the computational domain [−1, 1], i.e., and by introducing can be rewritten as We consider equations ( 28)-( 29) collocated at Chebyshev-Gauss-Lobatto points and φ expanded as Chebyshev polynomials truncated series The problem ( 28)-( 29) is rewritten as where We discretize φ at the collocation points ( 30) and so functions s, q, U .Then, we approximate the derivative d/dz through the Chebyshev differentiation matrix D N (see [37] for the definition).
We denote as the discretized operators ( 32)-( 38), i.e. (N + 1) × (N + 1) matrices with complex entries.We recall that for the implementation of the boundary conditions we follow [38].The discretized problem is a quadratic matrix eigenvalue problem of the type where  .
It is worth highlighting that the problem (40) (equivalent to (39)) is twice as large as (39).To solve the system we use MATLAB ® environment through QZ decomposition in the case of Papanastasiou and "exact" Bingham model.One must be careful to choose the proper N , since the round-off errors begin to accumulate for "large" N .We maintain high accuracy by using (39) with a smaller value of N to determine a good approximation of c, then to improve such a value we use the Newton-QR inverse iteration method (described in [39]).

Results of the analytical and numerical approaches
In this section we report the main results obtained by using both analytical and numerical approaches.More details are given in [32,33,34].Figure 2 shows the evolution of Re c as a function of the tilt angle θ for different values of λ for each "regularized" model.In particular, we obtain coherent analytical and numerical results by comparing the Bingham-type materials, i.e. the Regularized Bingham-like (λ = 0, 0.01, 0.1, 0.5 and = 0.01, analytical approach through long-wave approximation) and Papanastasiou model (λ = 0, 0.4, 0.8, 1.5, 3 and = 0.1, numerical approach through spectral method) as depicted in Figure 2(a,b).In particular, we recall that represents the level of approximation for each "regularized" model and the choice of its value depends on the typical values of λ ∈ 10 −2 , 1 [34].In particular, for Bingham-type materials, recalling the relation of proportion (14) of λ and τ 0 , the increase of yield stress leads to a destabilization effect on the flow.On the contrary, the results obtained in the case of Casson-type model by using the long-wave approximation method show that the increase of yield stress has a flow stabilizing effect.Moreover, for each "regularized" model the increase of θ leads to a decrease of Re c and we retrieve the benchmark (Newtonian case), i.e. the classical relation between the critical Reynolds number and the tilt angle when λ = 0.In particular, the Newtonian flow is more stable than a Bingham-type flow, while it is less stable than Casson-type flow.We validate our approach by compare the results obtained through the long-wave approximation with experimental data (see [40]) in the case of Regularized Bingham-like material with a rather satisfactory agreement.The details are given in [32].By using the numerical approach through the spectral collocation method, the evolution of Im (c) with respect to α and Re for selected values of Γ, i.e. selected geometrical and "material" properties see expression (14), is provided in Figure 3.In particular, the imaginary part of c, Im (c), is negative for any Γ, α, Re.Therefore, our findings highlight that the "exact" Bingham (unconditionally stable) and the Papanastasiou model (instability can arises above a critical Reynolds number) can have different stability characteristics.

Conclusions and final remarks
In this work, we investigate the stability of Bingham and Casson-type flow down an inclined plane theoretically (Regularized Bingham and Casson-like) and numerically (the Papanastasiou and "exact" Bingham) so that the all the issue concerning the existence of truly unyielded zones are avoided.In particular, the "exact" constitutive models has been approximated through the introduction of a regularized parameter, .The benchmark, i.e. the classical relation between Re c and θ ( [3,4]), has been recovered when λ = 0 (Newtonian case).Coherent results has been obtained by comparing analytical and numerical results (Regularized Bingham-like and Papanastasiou).In the case of Bingham-type materials, we compare our results with experimental data and with the "exact" Bingham flow to validate our approach and to investigate possible discrepancies with the "exact" one, respectively.In particular, we obtain a good comparison of the Regularized Bingham-like [32] with the experimental data [40].The numerical analysis of the stability for the "exact" Bingham model and of the Papanastasiou model shows that for every Reynolds number the first one is unconditionally stable, while for the second one instability may arise above a critical Reynolds number.Such discrepancies has been already observed in the case of plane Poiseuille flow [31].Moreover, our results highlight that in the case of Regularized Casson-like fluid Re c increases when τ 0 increases (i.e. through the relation (14) with the "material" parameter λ) with a stabilizing effect.While for a Regularized Bingham-like fluid Re c decreases when τ 0 increases, i.e. there is a destabilizing effect [32,34].

Figure 2 : 12 Γ
Figure 2: Critical Reynolds values as a function of θ for given value of λ and when the fluids is modelled as (a) Papanastasiou model; (b) Regularized Bingham-like model; (c) Regularized Casson-like model.

Figure 3 :
Figure 3: Evolution of Im(c) for any value of Re, Γ = λ/(sin θ) 2/3 and α when the fluids is modelled as the "exact" Bingham model.