<figure><img src="image_1.png" /></figure>

HAL
open science

Linear stability of viscoplastic flows down an incline

Benedetta Calusi, Angiolo Farina, Lorenzo Fusi, Liviu Iulian Palade, Fabio Rosso

<figure><img src="image_2.png" /></figure>

► To cite this version:

Benedetta Calusi, Angiolo Farina, Lorenzo Fusi, Liviu Iulian Palade, Fabio Rosso. Linear stability of viscoplastic flows down an incline. Journal of Physics: Conference Series, 2024, 2701 (1), pp.012071. 10.1088/1742-6596/2701/1/012071. [hal-04769880](https://hal.science/hal-04769880v1)

HAL Id: hal-04769880

https://hal.science/hal-04769880v1

Submitted on 24 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

Linear stability of viscoplastic flows down an
incline

To cite this article: Benedetta Calusi et al 2024 J. Phys.: Conf. Ser. 2701 012071

View the [article online](#) for updates and enhancements.

You may also like

- A novel test technology for a standardized stability control system and its test platform
  Biao Xu, Liqiang Xu, Yuxin Xu et al.

- (Invited) Electrochemistry at Metal-Oxide Interfaces: At the Crossroads of Activity and Stability
  Nenad M Markovic

- Study on Gradation and Performance of Cold Patch Asphalt Mixtures
  Li Qiang, Liu Feiyi, Luo Xiao et al.

<figure><img src="image_4.png" /></figure>
# Linear stability of viscoplastic flows down an incline

Benedetta Calusi¹, Angiolo Farina¹, Lorenzo Fusi¹, Liviu Iulian Palade² and Fabio Rosso¹

¹Università degli Studi di Firenze, Dipartimento di Matematica e Informatica “Ulisse Dini”, Italy

²CNRS, Institute “Camille Jordan” UMR 5208, INSA-Lyon, Université de Lyon, 69621 Villeurbanne, France

E-mail: benedetta.calusi@unifi.it

**Abstract.** 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.

## 1. 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, $\tau_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 viscoplastic 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.

The pioneering works reported in [3, 4] paved the way for the analysis of stability in the Newtonian case, by providing a proportionality relation, experimentally validated in [5], between $Re_c$ and the tilt angle $\theta$. In fact, the number of studies in the theoretical, numerical, and experimental fields with the aim to properly describe fluids with complex rheological behaviour are increasing, see e.g. [6, 7, 8, 9, 12, 10, 11, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24].

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

<figure><img src="image_5.png" /></figure>
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.

## 2. 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 $\theta \in (0, \pi/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 + \tau$ the Cauchy stress tensor and as $\tau$ the deviatoric part.

We recall that for the two-dimensional incompressible flow, $v = u\mathbf{e}_x + v\mathbf{e}_y$, we have in dimensional form

$$ \begin{cases} \rho \dot{\mathbf{v}} = -\nabla p + \text{div}(\boldsymbol{\tau}) + \mathbf{g}, \\ \text{div}(\mathbf{v}) = 0, \end{cases} \qquad (1) $$

where $\mathbf{g} = \rho g \sin \theta \mathbf{e}_x - \rho g \cos \theta \mathbf{e}_y$, $g$ is gravity and $\rho$ 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 $\mathbf{n}$ the outer normal (see Fig. 1), the kinematical-dynamical conditions is

$$ \begin{cases} h_t + uh_x = v, & y = h, \\ \mathbb{T}\mathbf{n} = 0, & y = h. \end{cases} \qquad (2) $$

<figure><img src="image_6.png" /><figcaption>Figure 1: Reference framework.</figcaption></figure>

We proceed by introducing the adimensionalization

$$ \mathbf{x}^* = \frac{\mathbf{x}}{H}, \quad \mathbf{v}^* = \frac{\mathbf{v}}{U}, \quad t^* = \frac{U}{H}t, \qquad (3) $$
$$h^* = \frac{h}{H}, \quad p^* = \frac{p}{p_c}, \quad \tau^* = \frac{\tau}{\tau_c}, \quad \dot{\gamma}^* = \frac{H}{U}\dot{\gamma}, \qquad (4)$$

where $\dot{\gamma} = 1/2 (\nabla v + \nabla v^T)$ is the strain-rate,

$$\tau_c = \frac{\mu U}{H}, \quad p_c = \frac{\mu U}{H}, \qquad (5)$$

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

$$\xi = \frac{\mathrm{Re}}{\mathrm{Fr}^2} \sin \theta = \frac{\rho g H^2}{\mu U} \sin \theta, \qquad (6)$$

where the Reynolds number and Freude number are

$$\mathrm{Re} = \frac{\rho U H}{\mu}, \quad \mathrm{Fr}^2 = \frac{U^2}{gH}, \qquad (7)$$

respectively, we have

$$\left\{ \begin{aligned} \mathrm{Re} (u_{t*}^* + u^* u_{x*}^* + v^* u_{y*}^*) &= -p_{x*}^* + \tau_{11,x*}^* + \tau_{12,y*}^* + \xi, \\ \mathrm{Re} (v_{t*}^* + u^* v_{x*}^* + v^* v_{y*}^*) &= -p_{y*}^* + \tau_{12,x*}^* + \tau_{22,y*}^* - \xi \cot \theta, \\ u_{x*}^* + v_{y*}^* &= 0. \end{aligned} \right. \qquad (8)$$

From (6) and (7) we have

$$\begin{aligned} H &= \left( \frac{\xi \mu^2}{\rho^2 g \sin \theta} \mathrm{Re} \right)^{1/3}, \\ U &= \frac{\mu}{\rho H} \mathrm{Re} = \left( \frac{g \mu \sin \theta}{\xi \rho} \right)^{1/3} \mathrm{Re}^{2/3}. \end{aligned} \qquad (9)$$

Now, we consider $h = 1$ and the flow as

$$v = u(y)e_x, \qquad (10)$$

thus, we can rewrite (8) as

$$\left\{ \begin{aligned} 0 &= -p_{x*}^* + \tau_{12,y*}^* + \xi, \\ 0 &= -p_{y*}^* - \xi \cot \theta, \end{aligned} \right. \qquad (11)$$

and $p = \tau_{12} = 0$ on $y^* = 1$, since $\tau_{12}^*$ is the only non-vanishing component of $\tau^*$.

### 3. “Regularized” Models

Here, we summarize the “regularized” constitutive laws in dimensionless form for a Bingham and Casson fluids through a regularization parameter, $\epsilon$. The corresponding “exact” models and
more details are in [32, 33, 34].

We consider the chosen dimensionless regularization constitutive models

$$ \tau^* = r(|\dot{\gamma}^*|) \dot{\gamma}^*, $$

$$ r(|\dot{\gamma}^*|) = \begin{cases} \frac{B}{|\dot{\gamma}^*| + \epsilon}, & (\text{Regularized Bingham - like}) \\ 2 + \left\{ \begin{aligned} & B \left( \frac{1 - e^{-\frac{2|\dot{\gamma}^*|}{\epsilon}}}{|\dot{\gamma}^*|} \right), & (\text{Papanastasiou model}) \\ & \left( \sqrt{2} + \frac{\sqrt{B}}{\sqrt{|\dot{\gamma}^*|} + \sqrt{\epsilon}} \right)^2, & (\text{Regularized Casson - like}) \end{aligned} \right. \end{cases} \quad (12) $$

where $|\dot{\gamma}^*| = \sqrt{tr(\dot{\gamma}^{*2})}/2$, $B = \frac{\tau_0 H}{\mu U}$ is the Bingham number where $\tau_0$ is the yield stress, and $\epsilon^* = 2\epsilon H/U$, so that $\epsilon \to 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

$$ B = \Gamma \frac{\xi^{2/3}}{\mathrm{Re}^{1/3}}, \quad (13) $$

with

$$ \Gamma = \frac{\lambda}{(\sin \theta)^{2/3}}, \quad \lambda = \frac{\tau_0^*}{(g^* \mu^*)^{2/3} \rho^{*1/3}}. \quad (14) $$

The parameter $\Gamma$ 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

$$ p = \xi \cot \theta (1 - y), \quad (15) $$

and, noting that $2|\dot{\gamma}| = u_y(y)$ and by using (12)

$$ \tau_{12} = r(u_y(y)) u_y(y), $$

$$ r(u_y(y)) = \begin{cases} \frac{B}{u_y(y) + \epsilon}, & (\text{Regularized Bingham - like}) \\ 1 + \left\{ \begin{aligned} & B \left( \frac{1 - e^{-\frac{u_y(y)}{\epsilon}}}{u_y(y)} \right), & (\text{Papanastasiou model}) \\ & \left( 1 + \frac{\sqrt{B}}{\sqrt{u_y(y)} + \sqrt{\epsilon}} \right)^2, & (\text{Regularized Casson - like}) \end{aligned} \right. \end{cases} \quad (16) $$
thus, we have

$$ \xi(1-y) = \begin{cases} \frac{B}{u_y(y) + \epsilon}, & (\text{Regularized Bingham-like}) \\ 1 + \left\{ \begin{aligned} & B \left( \frac{1 - e^{-\frac{u_y(y)}{\epsilon}}}{u_y(y)} \right), & (\text{Papanastasiou model}) \\ & \left( 1 + \frac{\sqrt{B}}{\sqrt{u_y(y)} + \sqrt{\epsilon}} \right)^2. & (\text{Regularized Casson-like}) \end{aligned} \right. & \end{cases} \quad (17) $$

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 $\mathbb{R}$ and $\xi$ for each “regularized” model, namely

$$ \mathcal{F}(\mathbb{R}, \xi) = 1, \quad (18) $$

meaning that i.e., for every positive $\mathbb{R}$, there exists $\xi = \xi(\mathbb{R})$ fulfilling (18) for given $\lambda, \epsilon, \theta$. The explicit expression of $u(y)$ and $\mathcal{F}(\mathbb{R}, \xi)$ are reported in [32, 33, 34] for each “regularized” model.

## 4. 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)\mathbf{i}$ with $u_b$ (obtained through (17) and reported in [32, 33, 34]), and, $p = p_b(y)$ where $p_b(y) = \xi \cot \theta(1 - y)$ (see (15)). First, we consider a small perturbation of the basic flow, in the form of travelling waves, so that

$$ \begin{aligned} h &= 1 + \hat{h}(y)e^{i\alpha(x-ct)}, & u &= u_b + \hat{u}(y)e^{i\alpha(x-ct)}, \\ v &= \hat{v}(y)e^{i\alpha(x-ct)}, & p &= p_b + \hat{p}(y)e^{i\alpha(x-ct)}, \end{aligned} \quad (19) $$

thus,

$$ \gamma = \gamma_b + \hat{\gamma}, \quad \tau = \tau_b + \hat{\tau}, \quad (20) $$

where we denote the wave number as $\alpha \in \mathbb{R}$, the complex wave speed as $c \in \mathbb{R}$ and the infinitesimal disturbance with the notation $\hat{\cdot}$. Then, the velocity field can be rewritten as stream function, namely

$$ \hat{\psi}(x, y, t) = \phi(y)e^{i\alpha(x-ct)}, $$

as

$$ \hat{u} = \hat{\psi}_y = \phi'(y) e^{i\alpha(x-ct)}, \quad \hat{v} = -\hat{\psi}_x = -i\alpha\phi(y)e^{i\alpha(x-ct)}, \quad (21) $$

where, here and in the sequel, $\hat{\cdot}(\cdot)'$ represents the differentiation with respect to $y$.

Recalling that the growth/attenuation factor of the $\alpha^{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 $(\mathbb{R}, \lambda, \epsilon$ and $\theta$). 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 $\mathbb{R}, \lambda, \epsilon$ and $\theta$.
By using (19)-(21), the corresponding perturbed problem reduces to

$$
\left\{
\begin{array}{l}
i\alpha \operatorname{Re} [(u_b - c)(\phi'' + \alpha^2\phi) - u_b''\phi] = -2\alpha^2 [q\phi']' + \\
\qquad + [s(\phi'' + \alpha^2\phi)]'' + \alpha^2 [s(\phi'' + \alpha^2\phi)], \\[2ex]
\hat{h}(1-c) = -\phi'(1), \\[2ex]
\hat{h}p_b'(1) + \hat{p}'(1) - \hat{\tau}_{22}(1) = 0, \\[2ex]
\hat{h}\tau_{b12}'(1) + \hat{\tau}_{12}(1) = 0, \\[2ex]
\hat{p}(1) = 2i\alpha\phi'(1) + \frac{\phi'''(1) + \alpha\phi'(1)}{i\alpha} - \operatorname{Re}\phi'(1)(1-c),
\end{array}
\right.
\qquad (22)
$$

which is a boundary value problem and

$$
q(y) = r \left( \frac{u'_{b}(y)}{2} \right), \quad s(y) = \frac{1}{2} \left[ r \left( \frac{u'_{b}}{2} \right) + \frac{u'_{b}}{2} \frac{dr}{dz} \left( \frac{u'_{b}}{2} \right) \right], \quad (23)
$$

for each “regularized” model given by (16). It is worth noting that if $\mathbf{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 \in \mathbb{C}$ and the
corresponding eigenfunctions) allow one to determine $\mathbf{Re}_c$ so that $\mathrm{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]).

4.1. Analytical approach: Long-wave approximation

In this section we briefly report the main characteristics of the long-wave approximation
approach. Such method consists in disturbances of long wavelength $2\pi/\alpha \gg 1$ ($\alpha \ll 1$) through
the expansion of both $\phi$ and $c$ in powers of $\alpha$ (up to the first order in $\alpha$)

$$
\begin{equation}
\begin{aligned}
\phi(y) &= \phi_0(y) + \alpha\phi_1(y), \\
c &= c_0 + \alpha c_1,
\end{aligned}
\tag{24}
\end{equation}
$$

where $(\phi_0, c_0)$ and $(\phi_1, c_1)$ fulfil system

$$
\left\{
\begin{array}{l}
(s(y)\phi_0(y)')'' = 0, \\
\\
\phi_0(0) = \phi'_0(0) = 0, \\
\\
\phi_0''(1) - \phi_0(1) \frac{\hat{\xi}}{s(1)(c_0-1)} = 0, \\
\\
s(1)\phi_0'''(1) + \phi_0(1) \frac{\hat{\xi}s'(1)}{s(1)(c_0-1)} = 0,
\end{array}
\right.
\qquad (25)
$$
and system

$$ 
\left\{
\begin{array}{l}
(s(y)\phi_1''(y))'' = i\operatorname{Re}[\phi_0''(y)(u_b(y) - c_0) - u_b''(y)\phi_0(y)], \\
\phi_1(0) = \phi_1'(0) = 0, \\
\phi_1(1) = 0, \\
\phi_1''(1) + \phi_0(1)\frac{\hat{\xi}c_1}{s(1)(c_0-1)^2} = 0, \\
-is(1)\phi_1'''(1) + \phi_0'(1)\operatorname{Re}(c_0-1) - \phi_0(1)\frac{\hat{\xi}}{s(1)}\left(\frac{s(1)\cot\theta}{c_0-1} - i\frac{c_1s'(1)}{(c_0-1)^2}\right) = 0.
\end{array}
\right.
\qquad (26) $$

Since $\operatorname{Im}(c_0) = 0$ and $\operatorname{Re}(c_1) = 0$ (see [32, 33]), the value of $\operatorname{Re}_c$ is such that

$$ \operatorname{Im}(c) = \operatorname{Im}(c_1(\operatorname{Re}_c, \lambda, \epsilon, \theta)) = 0, \qquad (27) $$

through the prescription of the material characteristics and the tilt angle (i.e., $\lambda$, $\epsilon$, and $\theta$). Hence, the $\alpha^{th}$ mode is stable when $\operatorname{Re} < \operatorname{Re}_c$, while if $\operatorname{Re} > \operatorname{Re}_c$ instability occurs. Indeed, $\operatorname{Im}(c_1) < 0$ when $\operatorname{Re} < \operatorname{Re}_c$ and vice versa.

We compute the value of $\operatorname{Re}_c$ through the solution of the system given by (18) and (27) by using MATLAB® 2022a by means of function FSOLVE.

## 4.2. 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.,

$$ y = \frac{z+1}{2}, \quad z = -1+2y \qquad \frac{d^n}{dy^n} = 2^n \frac{d^n}{dz^n}. $$

and by introducing

$$ \phi(z) = \hat{\psi}(y), \quad U(z) = u_b(y), \quad q(z) = q(y), \quad s(z) = s(y), $$

problem (22) can be rewritten as

$$ i\alpha\operatorname{Re}\left[(U-c)\left(4\frac{d^2\phi}{dz^2} - \alpha^2\phi\right) - 4\frac{d^2U}{dz^2}\phi\right] = -8\alpha^2\frac{d}{dz}\left[q\frac{d\phi}{dz}\right] + \\ + 4\frac{d^2}{dz^2}\left[s\left(4\frac{d^2\phi}{dz^2} + \alpha^2\phi\right)\right] + \alpha^2\left[s\left(4\frac{d^2\phi}{dz^2} + \alpha^2\phi\right)\right], \quad (28) $$
$$ 
\left\{
\begin{array}{ll}
\phi = \frac{d\phi}{dz} = 0, & y = -1, \\
\\
4\frac{d^2\phi}{dz^2} + \phi\left[\alpha^2 - \frac{\xi}{s(c-1)}\right] = 0, & y = 1, \\
\\
-8is\frac{d^3\phi}{dz^3} + 2\alpha\frac{d\phi}{dz}\left[i\alpha(2q-s) + \text{Re}(c-1)\right] + \\
\qquad -\phi\left[(\alpha s \cot\theta + 2i\frac{ds}{dz})\frac{\xi}{s(c-1)}\right] = 0, & y = 1.
\end{array}
\right.
\qquad (29) $$

We consider equations (28)-(29) collocated at Chebyshev-Gauss-Lobatto points

$$ z_k = \cos\left(\frac{k\pi}{N}\right), \quad k = 0, \dots, N, \qquad (30) $$

and $\phi$ expanded as Chebyshev polynomials truncated series

$$ \phi(z) = \sum_{k=0}^{N} a_k T_k(z). $$

The problem (28)-(29) is rewritten as

$$ 
\left\{
\begin{array}{l}
\mathcal{L}_1\phi = c\mathcal{L}_2\phi, \\
\\
\mathcal{A}_1\phi = 0, \\
\\
\mathcal{B}_1\phi = 0, \\
\\
\mathcal{M}_1\phi = c\mathcal{M}_2\phi, \\
\\
\mathcal{N}_1\phi = c\mathcal{N}_2\phi + c^2\mathcal{N}_3\phi,
\end{array}
\right.
\qquad (31) $$

where

$$ \mathcal{L}_1 = (16s) \frac{d^4}{dz^4} + \left(32 \frac{ds}{dz}\right) \frac{d^3}{dz^3} - \left(8\alpha^2(q-s) - 16 \frac{d^2s}{dz^2} + 4i\alpha\text{Re}U\right) \frac{d^2}{dz^2} + \\ + \left(8\alpha^2 \frac{ds}{dz} - 8\alpha^2 \frac{dq}{dz}\right) \frac{d}{dz} + \left(4\alpha^2 \frac{d^2s}{dz^2} + \alpha^4 s + i\alpha^3\text{Re}U + 4i\alpha\text{Re}\frac{d^2U}{dz^2}\right), \qquad (32) $$

$$ \mathcal{L}_2 = (-4i\alpha\text{Re}) \frac{d^2}{dz^2} + i\alpha^3\text{Re}, \qquad (33) $$

$$ \mathcal{A}_1 = 1, \qquad \mathcal{B}_1 = \frac{d}{dz}, \qquad (34) $$

$$ \mathcal{M}_1 = (4s) \frac{d^2}{dz^2} + \alpha^2 s + \xi, \qquad \mathcal{M}_2 = (4s) \frac{d^2}{dz^2} + \alpha^2 s, \qquad (35) $$

$$ \mathcal{N}_1 = (8s^2) \frac{d^3}{dz^3} - \left[2\alpha^2(2qs - s^2) + 2i\alpha s \text{Re}\right] \frac{d}{dz} + \left(i\alpha s \xi \cot\theta - 2 \frac{ds}{dz}\xi\right), \qquad (36) $$
$$ \mathcal{N}_2 = (8s^2) \frac{d^3}{dz^3} - [2\alpha^2(2qs - s^2) + 4i\alpha s \operatorname{Re}] \frac{d}{dz}, \qquad (37) $$

$$ \mathcal{N}_3 = (2i\alpha s \operatorname{Re}) \frac{d}{dz}. \qquad (38) $$

We discretize $\phi$ 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

$\mathbf{L}_1, \mathbf{L}_2, \mathbf{A}_1, \mathbf{B}_1, \mathbf{M}_1, \mathbf{M}_2, \mathbf{N}_1, \mathbf{N}_2, \mathbf{N}_3$.

the discretized operators (32)-(38), i.e. $(N+1) \times (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

$$ \mathbf{P}\phi = c\mathbf{R}\phi + c^2\mathbf{W}\phi. \qquad (39) $$

where

$$ \mathbf{P} = \begin{bmatrix} (\mathbf{N}_1)_{0,0} & \cdots & (\mathbf{N}_1)_{0,N} \\ (\mathbf{M}_1)_{0,0} & \cdots & (\mathbf{M}_1)_{0,N} \\ (\mathbf{L}_1)_{2,0} & \cdots & (\mathbf{L}_1)_{2,N} \\ \vdots & \ddots & \vdots \\ (\mathbf{L}_1)_{N-2,0} & \cdots & (\mathbf{L}_1)_{N-2,N} \\ (\mathbf{A}_1)_{0,0} & \cdots & (\mathbf{A}_1)_{0,N} \\ (\mathbf{B}_1)_{0,0} & \cdots & (\mathbf{B}_1)_{0,N} \end{bmatrix}, $$

$$ \mathbf{R} = \begin{bmatrix} (\mathbf{N}_2)_{0,0} & \cdots & (\mathbf{N}_2)_{0,N} \\ (\mathbf{M}_2)_{0,0} & \cdots & (\mathbf{M}_2)_{0,N} \\ (\mathbf{L}_2)_{2,0} & \cdots & (\mathbf{L}_2)_{2,N} \\ \vdots & \ddots & \vdots \\ (\mathbf{L}_2)_{N-2,0} & \cdots & (\mathbf{L}_2)_{N-2,N} \\ 0 & \cdots & 0 \\ 0 & \cdots & 0 \end{bmatrix}, $$

$$ \mathbf{W} = \begin{bmatrix} (\mathbf{N}_3)_{0,0} & \cdots & (\mathbf{N}_3)_{0,N} \\ 0 & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & 0 \end{bmatrix}. $$

Then, we transform problem (39) into the linear generalized eigenvalue problem

$$ \begin{bmatrix} \mathbf{P} & \mathbf{0} \\ \mathbf{0} & \mathbf{I} \end{bmatrix} \begin{bmatrix} \phi \\ c\phi \end{bmatrix} = c \begin{bmatrix} \mathbf{R} & \mathbf{W} \\ \mathbf{I} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \phi \\ c\phi \end{bmatrix}, \qquad (40) $$

where we denote as $\mathbf{0}$ the $(N+1) \times (N+1)$ null and as $\mathbf{I}$ the $(N+1) \times (N+1)$ identity matrix. 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]).
**5. 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 $\theta$ for different values of $\lambda$ 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 ($\lambda = 0, 0.01, 0.1, 0.5$ and $\epsilon = 0.01$, analytical approach through long-wave approximation) and Papanastasiou model ($\lambda = 0, 0.4, 0.8, 1.5, 3$ and $\epsilon = 0.1$, numerical approach through spectral method) as depicted in Figure 2(a,b). In particular, we recall that $\epsilon$ represents the level of approximation for each “regularized” model and the choice of its value depends on the typical values of $\lambda \in [10^{-2}, 1]$ [34]. In particular, for Bingham-type materials, recalling the relation of proportion (14) of $\lambda$ and $\tau_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 $\theta$ 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 $\lambda = 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 comparing 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 $\alpha$ and $Re$ for selected values of $\Gamma$, 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 $\Gamma, \alpha, 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.

**6. 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, $\epsilon$. The benchmark, i.e. the classical relation between $Re_c$ and $\theta$ ([3, 4]), has been recovered when $\lambda = 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 $\tau_0$ increases (i.e. through the relation (14) with the “material” parameter $\lambda$) with a stabilizing effect. While for a Regularized Bingham-like fluid $Re_c$ decreases when $\tau_0$ increases, i.e. there is a destabilizing effect [32, 34].

Therefore, our analytical and numerical analysis (2D linear stability analysis) highlights that, although in a “regularized” formulation, the stability characteristics of models of the same class of viscoplastic fluids can be completely different [33]. Moreover, our findings can pave the way to deepen complex fluids with further more exhaustive stability analysis through the synergy of
<figure><img src="image_7.png" /><figcaption>(a) Papanastasiou model.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) Regularized Bingham-like model.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c) Regularized Casson-like model.</figcaption></figure>

Figure 2: Critical Reynolds values as a function of $\theta$ for given value of $\lambda$ and $\epsilon$ when the fluids is modelled as (a) Papanastasiou model; (b) Regularized Bingham-like model; (c) Regularized Casson-like model.
<figure><img src="image_10.png" /><figcaption>Figure 3: Evolution of Im(c) for any value of Re, $\Gamma = \lambda/(\sin\theta)^{2/3}$ and $\alpha$ when the fluids is modelled as the “exact” Bingham model.</figcaption></figure>

theoretical and numerical studies coupled with experiments.

## 7. Acknowledgements

This research was partially supported by GNFM of Italian INDAM and by National Recovery and Resilience Plan, Mission 4 Component 2 - Investment 1.4 - NATIONAL CENTER FOR HPC, BIG DATA AND QUANTUM COMPUTING - funded by the European Union - NextGenerationEU - CUP B83C22002830001).

## References

[1] Farina, A., Fusi, L.: Viscoplastic Fluids: Mathematical Modeling and Applications, Non-Newtonian Fluid Mechanics and Complex Flows: Levico Terme, Italy 2016, Springer International Publishing, 2018, 229–298.

[2] Huilgol, R.R.: Fluid Mechanics of Viscoplasticity, Springer, 2015.

[3] Benjamin, T.B.: Wave formation in laminar flow down an inclined plane, J. Fluid Mech., 1957, 2, 554.

[4] Yih, C.-S.: Stability of liquid flow down an inclined plane, Phys. Fluids, 1963, 6, 321.

[5] Liu, J., Paul, J. D., Gollub, J. P.: Measurements of the primary instabilities of film flows, J. Fluid Mech., 1993, 250, 69-101.

[6] Allouche, M.H., Millet, S., Botton, V., Henry, D., Hadid, H.B., Rousset, F.: Stability of a flow down an incline with respect to two-dimensional and three-dimensional disturbances for Newtonian and non-Newtonian fluids. Phys. Rev. E, 2015, 92, 063010.

[7] Allouche, M.H., Botton, V., Millet, S., Henry, D., Dagois-Bohy, S., Güzel, B., Hadid, H. B.: Primary instability of a shear-thinning film flow down an incline: Experimental study, J. Fluid Mech., 2017, 821.

[8] Calusi, B., Fusi, L., Farina, A.: On a free boundary problem arising in snow avalanche dynamics, ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik, 2016, 453–465.

[9] Chakraborty, S., Sheu, T. W.-H., Ghosh, S.: Dynamics and stability of a power-law film flowing down a slippery slope, Phys. Fluids, 2019, 31, 013102.

[10] Forterre, Y., Pouliquen, O.: Long-surface-wave instability in dense granular flows, J. Fluid Mech., 2003, 486, 21–50.

[11] Fusi, L.: Channel flow of viscoplastic fluids with pressure-dependent rheological parameters, Phys. Fluids, 2018, 30, 073102.
[12] Falsaperla, P., Giacobbe, A., Mulone, G.: Stability of the plane Bingham-Poiseuille flow in an inclined channel, Fluids, 2020, 5, 141.

[13] Hu, J., Millet, S., Botton, V., Hadid, H. B., Henry, D.: Inertialess temporal and spatio-temporal stability analysis of the two-layer film flow with density stratification, Phys. Fluids, 2006, 18, 104101.

[14] Hu, J., Yin, X.Y., Hadid, H.B., Henry, D.: Linear temporal and spatiotemporal stability analysis of two-layer falling films with density stratification, Phys. Rev. E, 2008, 77, 026302.

[15] Hu, J., Hadid, H.B., Henry, D., Mojtabi, A.: Linear temporal and spatiotemporal stability analysis of a binary liquid film flowing down an inclined uniformly heated plate, J. Fluid Mech., 2008, 599, 269-298.

[16] Hu, T., Fu, Q.-fei, Xing, Y., Yang, L.-jun, Xie, L.: Stability of a thin viscoelastic film falling down an inclined plane, Phys. Rev. Fluids, 2021, 6, 083902.

[17] Millet, S., Botton, V., Rousset, F., Hadid, H.B.: Wave celerity on a shear thinning fluid film flowing down an incline, Phys. Fluids, 2008, 20, 031701.

[18] Millet, S., Botton, V., Hadid, H. B., Henry, D., Rousset, F.: Stability of two-layer shear-thinning film flows, Phys. Rev. E, 2013, 88, 043004.

[19] Millet, S., Usha, R., Botton, V., Rousset, F.: The mechanism of long-wave instability in a shear-thinning film flow on a porous substrate, Acta Mech., 2019, 230, 2201-2220.

[20] Mogilevskiy, E.: Stability of a non-Newtonian falling film due to three-dimensional disturbances, Phys. Fluids, 2020, 32, 073101.

[21] Ng, C.-O., Mei, C.C.: Roll waves on a shallow layer of mud modelled as a power-law fluid, J. Fluid Mech., 1994, 263, 151-184.

[22] Nsom, B., Ramifidisoa, L., Latrache, N., Ghaemizadeh, F.: Linear stability of shear-thinning fluid down an inclined plane, J. Mol. Liquids, 2019, 277, 1036-1046.

[23] Rousset, F., Millet, S., Botton, V., Hadid, H.B.: Temporal stability of carreau fluid flow down an incline, J. Fluids Eng., 2007, 129, 913-920.

[24] Ruyer-Quil, C., Chakraborty, S., Dandapat, B.S.: Wavy regime of a power law film flow, J. Fluid Mech., 2012, 692, 220-256.

[25] Astarita, G.: Letter to the editor: The engineering reality of the yield stress, J. Rheol., 1990, 34, 275-277.

[26] Barnes, H.A., Walters, K.: The yield stress myth?, Rheologica Acta, 1985, 24, 323-326.

[27] Barnes, H.A.: The yield stress – a review or ‘$\pi\alpha\nu\tau\alpha \rho\epsilon\iota$’ – everything flows?, J. Non-Newtonian Fluid Mech., 1999, 81, 133-178.

[28] Frigaard, I. and, Nouar, C.: On the usage of viscosity regularisation methods for visco-plastic fluid flow computation, J. Non-Newtonian Fluid Mech., 2005, 127, 1-26.

[29] Frigaard, I.A., Paso, K.G., de Souza Mendes, P.R.: Bingham's model in the oil and gas industry, Rheologica Acta, 2017, 56, 259-282.

[30] Bercovier, M., Engleman, M.: A finite-element method for incompressible non-Newtonian flows, J. Comput. Phys., 1980, 36, 313-326.

[31] Fusi, L., Farina, A., Rajagopal, K.R., Vergori, L.: Channel flows of shear-thinning fluids that mimic the mechanical response of a Bingham fluid, International Journal of Non-Linear Mechanics, 2022.

[32] Calusi, B., Farina, A., Fusi, L., Rosso, F.: Long-wave instability of a regularized Bingham flow down an incline. Physics of Fluids, 2022, 34 (5), 054111.

[33] Calusi, B., Farina, A., Fusi, L.: Stability of a Regularized Casson Flow down an Incline: Comparison with the Bingham Case. Fluids, 2022, 7, 380.

[34] Fusi, L., Calusi, B., Farina, A., Rosso, F.: Stability of laminar viscoplastic flows down an inclined open channel. European Journal of Mechanics-B/Fluid, 2022, 95, 137-147.

[35] Guadagli, S., Palade, L.I., Fusi, L., Farina, A.: On a Casson Fluid Motion: Nonuniform Width Symmetric Channel and Peristaltic Flows, Fluids, 2021, 6, 356.

[36] Pascal, J. P.: Linear stability of fluid flow down a porous inclined plane, Journal of Physics D: Applied Physics, 1999, 32 (4).

[37] Canuto, C., Hussaini, M. Y., Quarteroni, A., Zang, T.A.: Spectral methods in fluid dynamics. Springer-Verlag, New York, 1987.

[38] Horng, A.T.L., Chebyshev collocation method on solving stability of a liquid layer flowing down an inclined plane, The Mechanics of Thin Film Coatings, 1996, 95-106.

[39] Güttel, S., Tisseur F., The nonlinear eigenvalue problem, Acta Numerica, 2017, 26, 1-94.

[40] Mounkaila Noma, D., Dagois-Bohy, S., Millet, S., Botton, V., Henry, D., Ben Hadid, H.: Primary instability of a visco-plastic film down an inclined plane: Experimental study, J. Fluid Mech., 2021, 922, R2.