# Stability of a Regularized Casson Flow down an Incline: Comparison with the Bingham Case

Benedetta Calusi, Angiolo Farina, Lorenzo Fusi, Liviu-Iulian Palade

► To cite this version:

Benedetta Calusi, Angiolo Farina, Lorenzo Fusi, Liviu-Iulian Palade. Stability of a Regularized Casson Flow down an Incline: Comparison with the Bingham Case. Fluids, 2022, 7 (12), pp.380. 10.3390/fluids7120380. [https://doi.org/10.3390/fluids7120380](https://doi.org/10.3390/fluids7120380)

HAL Id: hal-04106226

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

Submitted on 11 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Stability of a Regularized Casson Flow down an Incline: Comparison with the Bingham Case

Benedetta Calusi <fnref n="1" />,* , Angiolo Farina <fnref n="1" />, , Lorenzo Fusi <fnref n="1" /> and Liviu Iulian Palade <fnref n="2" />

**Abstract:** In this paper, we study the two-dimensional linear stability of a regularized Casson fluid (i.e., a fluid whose constitutive equation is a regularization of the Casson obtained through the introduction of a smoothing parameter) flowing down an incline. The stability analysis has been performed theoretically by using the long-wave approximation method. The critical Reynolds number at which the instability arises depends on the material parameters, on the tilt angle as well as on the prescribed inlet discharge. In particular, the results show that the regularized Casson flow has stability characteristics different from the regularized Bingham. Indeed, for the regularized Casson flow an increase in the yield stress of the fluid induces a stabilizing effect, while for the Bingham case an increase in the yield stress entails flow destabilization.

**Keywords:** regularized Casson fluid; regularized Bingham fluid; linear stability analysis; long-wave approximation

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

check for
updates

**Citation:** Calusi, B.; Farina, A.; Fusi, L.; Palade, L.I. Stability of a Regularized Casson Flow down an Incline: Comparison with the Bingham Case. *Fluids* **2022**, *7*, 380. https://doi.org/10.3390/fluids7120380

Academic Editor: Mehrdad Massoudi

Received: 4 November 2022

Accepted: 6 December 2022

Published: 9 December 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

The rheological behaviour of materials such as suspensions, dispersion, and polymer solutions, is distinctly different from that of Newtonian fluids. In particular, such materials often exhibit flow properties characterized by a critical value of stress (i.e., yield stress, usually denoted as $\tau_0^*$), below which the materials do not deform, and above which they flow accordingly to their rheological properties. They are usually referred as viscoplastic materials, which include, e.g., the Bingham [1], the Herschel-Bulkley [2], and the Casson model [3].

The flow stability analysis of these models can have useful application in several industrial processes (e.g., food and pharmaceutical industries) and environmental phenomena (e.g., debris and lava flow). In general, flows are unstable when the corresponding Reynolds is larger than a critical threshold usually referred to as critical Reynolds number and denoted as $Re_c$. The pioneering works on stability of Newtonian flow down an incline has been reported in [4,5]. In these papers the authors provide a proportionality relation between the so-called critical Reynolds number, $Re_c$, and the tilt angle $\theta$ and later experimentally validated in [6]. Then, the interest to properly describe fluids with complex rheological behaviour led to an increase in theoretical, numerical, and experimental studies, see e.g., [7–32].

Recently, the onset of instability for viscoplastic fluids, flowing down an incline, has been investigated in [10,16]. In particular, a stability analysis has been performed numerically by using a spectral method in [16] and the long-wave approximation in [10], through a regularization of the Bingham law. The Bingham law describes a material characterized by the presence of a yield stress below which the continuum behaves like a rigid body and above which it flows as a linear viscous fluid. In this paper, we theoretically investigate the flow stability through the long-wave approximation technique following the approach reported in [10] and in [28]. In particular, we focus on a fluid modelled as

<footnote n="1">¹ Dipartimento di Matematica e Informatica “U. Dini”, Università degli Studi di Firenze, Viale Morgagni 67/A, 50134 Firenze, Italy</footnote>

<footnote n="2">² CNRS, Institute “Camille Jordan” UMR 5208, INSA-Lyon, Université de Lyon, 69621 Villeurbanne, France</footnote>

<footnote>* Correspondence: benedetta.calusi@unifi.it</footnote>
a regularized Casson, since it has the advantage of being easy to handle analytically. The Casson constitutive law is widely used to model blood flow [33]. Recently, studies regarding the peristaltic Casson flow (important to understand artery and vein physiology [34,35]), has been developed in [36–41].

In the “ideal” Casson model the stress is undetermined at zero strain rate. The presence of a yield stress has been widely discussed [42–46] and it is still an open debate. The use of a regularized model allows to avoid the problems due to this singularity and so to avoid several analytical and numerical issues [47,48]. Indeed, the singularity at zero strain rate can be smoothed out and the exact model can be recovered through introduction of a positive parameter, chosen quite arbitrarily, which accounts for the accuracy of the approximation [10,16,49,50].

To the best of the authors’ knowledge, the analysis of the onset of instability of a flow down an incline when the fluid is modelled as a regularized Casson material has not been presented in the literature before, and this motivates our investigations. Actually, the aim and novelty of this paper is two-fold. First, we study the stability properties of the regularized Casson flow down an incline. Then, we compare the obtained results with the one illustrated in [10] regarding the flow of a regularized Bingham. In particular, our findings highlight that the regularized Bingham fluid and the regularized Casson fluid have stability properties dramatically different. Indeed, although the two models belong to the same “viscoplastic family”, they show an opposite stability behaviour as the yield stress increases.

The paper is organized as follows: in Section 2 and 3 we formulate the mathematical problem and the main characteristics of a regularized Casson flow down an incline, respectively. In Section 4, following [10,16,51], we briefly recall linear stability analysis by using the long-wave approximation method. Then, in Section 5 and 6, we report results and some final remarks.

## 2. Mathematical Model

We proceed similarly to [10,16] briefly reporting the main theoretical background. Throughout the paper the “*” represents a dimensional quantity. Let us consider a reference framework $x^*Oy^*$ as the one depicted in Figure 1. We denote the tilt angle as $\theta \in (0, \pi/2)$ and suppose that the flow domain of the flow is given by

$$ \mathcal{D} = \{(x^*, y^*) \in \mathbb{R}^2 | 0 \le x^* \le L^*, 0 \le y^* \le h^*(x^*, t^*)\}, $$

where $L^*$ is the length of the domain and $y^* = h^*(x^*, t^*)$ is the upper free surface (not a priori known) and $H^* = \max\{h^*\}$.

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

We denote by $T^*$ the Cauchy stress tensor and set

$$ T^* = -p^*I + \tau^*, \quad (1) $$

where $\tau^*$ is the deviatoric part.
The governing equations for the two-dimensional incompressible flow, $\mathbf{v}^* = u^*\mathbf{i} + v^*\mathbf{j}$, are

$$ egin{cases} 
\rho^*(u_{t*}^* + u^*u_{x*}^* + v^*u_{y*}^*) = -p_{x*}^* + \tau_{11,x*}^* + \tau_{12,y*}^* + \rho^*g^* \sin \theta, \\ 
\rho^*(v_{t*}^* + u^*v_{x*}^* + v^*v_{y*}^*) = -p_{y*}^* + \tau_{12,x*}^* + \tau_{22,y*}^* - \rho^*g^* \cos \theta, \\ 
u_{x*}^* + v_{y*}^* = 0, 
\end{cases} 
\quad (2) $$

where $g^*$ is gravity and $\rho^*$ is the constant material density and, to take the notation as light as possible, we denote $(\cdot)_{t*} = \frac{\partial(\cdot)}{\partial t^*}$, $(\cdot)_{x*} = \frac{\partial(\cdot)}{\partial x^*}$, $(\cdot)_{y*} = \frac{\partial(\cdot)}{\partial y^*}$. We consider the non-slip and impermeability conditions $u^* = v^* = 0$ on $y^* = 0$ and the kinematical-dynamical conditions on $h^*$, namely

$$ egin{cases} 
h_{t*}^* + u^*h_{x*}^* = v^*, & y^* = h^*, \\ 
\mathbf{T}^*\mathbf{n} = 0, & y^* = h^*, 
\end{cases} 
\quad (3) $$

where $\mathbf{n}$ is the outer normal (see Figure 1). Exploiting (2)3, we rewrite (3)1 as

$$ h_{t*}^* + \left( \int_0^{h^*} u^* dy^* \right)_{x*} = 0. \quad (4) $$

We introduce the characteristic quantities

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

and the strain-rate $\dot{\gamma}^* = 1/2(\nabla^*\mathbf{v}^* + \nabla^*\mathbf{v}^{*T})$. Next, we consider the following dimensionless variables

$$ \mathbf{x} = \frac{\mathbf{x}^*}{H^*}, \quad \mathbf{v} = \frac{\mathbf{v}^*}{U^*}, \quad t = \frac{U^*}{H^*}t^*, \quad (6) $$

$$ h = \frac{h^*}{H^*}, \quad p = \frac{p^*}{p_c^*}, \quad \tau = \frac{\tau^*}{\tau_c^*}, \quad \dot{\gamma} = \frac{H^*}{U^*}\dot{\gamma}^*, \quad (7) $$

where $U^*$ denotes the reference velocity which will be selected to normalize the dimensionless longitudinal velocity. Exploiting (5)–(7), the system (2) becomes

$$ egin{cases} 
\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{cases} 
\quad (8) $$

where

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

and

$$ \mathrm{Re} = \frac{\rho^* U^* H^*}{\mu^*}, \quad Fr^2 = \frac{U^{*2}}{g^* H^*}, \quad (10) $$

are the Reynolds number and Freude number, respectively.
Finally, we recall that the flow is driven prescribing the inlet discharge which we assume to be constant in time. So, denoting as $Q^*$ the flow rate per unit fluid layer width, we have

$$Q^* = U^*H^* \int_0^1 u(y)dy, \qquad (11)$$

provided that the layer is flat and its thickness is $H^*$ (which does not vary in time). From (9) and (10) we have

$$H^* = \left( \frac{\xi \mu^{*2}}{\rho^{*2} g^* \sin \theta} \mathrm{Re} \right)^{1/3}, \qquad (12)$$

$$U^* = \frac{\mu^*}{\rho^* H^*} \mathrm{Re} = \left( \frac{g^* \mu^* \sin \theta}{\xi \rho^*} \right)^{1/3} \mathrm{Re}^{2/3}.$$

Therefore, the flow rate $Q^*$ can be rewritten in terms of $\mathrm{Re}$ and $\xi$ for given physical ($\tau_0^*, \mu^*, \rho^*$) and geometrical ($\theta$) parameters. In the sequel, we shall see that $Q^*$ can be also expressed only on terms of the Reynolds number through the normalization of the dimensionless longitudinal velocity.

### 3. Regularized Casson

Similarly to [10,16,52], we introduce a dimensional regularization parameter $\epsilon^*$ into the dimensionless Casson model, setting

$$\tau = r(|\dot{\gamma}|)\dot{\gamma}, \quad r(|\dot{\gamma}|) = \left( \sqrt{2} + \frac{\sqrt{B}}{\sqrt{|\dot{\gamma}|} + \sqrt{\epsilon}} \right)^2, \qquad (13)$$

where $|\dot{\gamma}| = \sqrt{\mathrm{tr}(\dot{\gamma}^2)/2}$, $B = \frac{\tau_0^* H^*}{\mu^* U^*}$ is the Bingham number, which represents the ratio between the yield stress and the characteristic viscous stress, with $\tau_0^*$ the yield stress, and $\epsilon = 2\epsilon^*H^*/U^*$, so that when $\epsilon^* \to 0$ we formally retrieve the Casson constitutive law [38]. We recall that for a regularized Bingham fluid the function $r$, see [10,16], is given by

$$r(|\dot{\gamma}|) = 2 + \frac{B}{|\dot{\gamma}| + \epsilon}. \qquad (14)$$

The plot of $r$, by using (13) and (14), with respect to $|\dot{\gamma}|$ for selected values of $B$ and $\epsilon$ is shown in Figure 2. In particular, the function $r(|\dot{\gamma}|)$ gives the regularized relation between the shear stress $\tau$ and shear strain rate $\dot{\gamma}$.

By using (9) and (10) the Bingham number and the regularization parameter can be rewritten in terms of $\mathrm{Re}$, namely

$$B = X \frac{\xi^{2/3}}{\mathrm{Re}^{1/3}}, \qquad (15)$$

with

$$X = \frac{\lambda_1}{(\sin \theta)^{2/3}}, \quad \lambda_1 = \frac{\tau_0^*}{(g^* \mu^*)^{2/3} \rho^{*1/3}}, \qquad (16)$$

and

$$\epsilon = \frac{\lambda_2}{(\sin \theta)^{2/3}} \frac{\xi^{2/3}}{\mathrm{Re}^{1/3}}, \qquad (17)$$

with

$$\lambda_2 = \frac{\epsilon^* \mu^{*1/3}}{\rho^{*1/3} g^{*2/3}}, \qquad (18)$$

respectively. The parameters $\lambda_1$ and $\lambda_2$ do not depend on the flow and on the tilt angle $\theta$, so that the parameter $X$ is constant once the fluid and the tilt angle have been selected, i.e., it depends only on the “material” and geometrical properties.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Plot of r(|$\dot{\gamma}$|) given by (13) and (14) for B = 1.1 and $\epsilon$ = 0.1. Although the regularized Bingham and the regularized Casson belong to the same viscoplastic family, they have different behaviour.</figcaption></figure>

We look for a solution in the form $v = u(y)i, h = 1$, thus system (8) is reduced to

$$ egin{cases} 0 = -p_x + \tau_{12,y} + \zeta, \\ 0 = -p_y - \zeta \cot \theta, \end{cases} \tag{19} $$

which leads to

$$ p = \zeta \cot \theta (1 - y). \tag{20} $$

Thus, we have $\tau_{12} = \zeta(1 - y)$, and, from (13)

$$ u_y \left( 1 + \frac{\sqrt{B}}{\sqrt{u_y} + \sqrt{\epsilon}} \right)^2 = \zeta(1 - y), $$

which, by integrating with respect to y with $u(0) = 0$, leads to

$$ u(y) = -\frac{\zeta}{4} y^2 - \frac{N_1(y)}{12\zeta^{\frac{3}{2}}} y - \frac{1}{12\zeta^{\frac{3}{2}}} (N_2(y) + N_3(y) + N_4(y)), \tag{21} $$

where

$$ N_1(y) = M_1(y)M_2(y) - 6\zeta^{\frac{5}{2}} + 8a(y) B^{\frac{1}{2}}\zeta^{\frac{3}{2}} - 6c, $$

$$ N_2(y) = -72\epsilon\zeta^{\frac{1}{2}}\left(\frac{2\epsilon^{\frac{1}{2}}B^{\frac{1}{2}}}{3} + B\right)\log\left(\epsilon^{\frac{1}{2}} - B^{\frac{1}{2}} + a(y) + M_2(y)\right), $$

$$ N_3(y) = M_2(y) \left\{ a(y) \left[ 3\zeta^{\frac{3}{2}} + \zeta^{\frac{1}{2}}(-3\epsilon + 8\epsilon^{\frac{1}{2}}B^{\frac{1}{2}} + B) \right] - b + \zeta^{\frac{1}{2}} \left[ B^{\frac{3}{2}} + \epsilon^{\frac{1}{2}}(3\epsilon - 47\epsilon^{\frac{1}{2}}B^{\frac{1}{2}} + 11B) \right] \right\}, $$

$$ N_4(y) = -8a(y) \zeta^{\frac{3}{2}}B^{\frac{1}{2}} + 8\zeta^2 B^{\frac{1}{2}} + M_3 + M_4, $$

$$ M_1(y) = -3a(y) \zeta^{\frac{3}{2}} + b, $$

$$ M_2(y) = \left[ 2a(y) \left( \epsilon^{\frac{1}{2}} - B^{\frac{1}{2}} \right) + c + a^2(y) \right]^{\frac{1}{2}}, \tag{22} $$

$$ M_3 = (3\epsilon^{\frac{1}{2}} + 5B^{\frac{1}{2}})\zeta^{\frac{3}{2}} - B^{\frac{3}{2}}\zeta^{\frac{1}{2}} + \epsilon^{\frac{1}{2}}\zeta^{\frac{1}{2}}(-3\epsilon + 47\epsilon^{\frac{1}{2}}B^{\frac{1}{2}} - 11B) - M_5\zeta(-3\epsilon + 8\epsilon^{\frac{1}{2}}B^{\frac{1}{2}} + B + 3\zeta), $$

$$ M_4 = 72\epsilon\zeta^{\frac{1}{2}}\left(-\frac{2\epsilon^{\frac{1}{2}}B^{\frac{1}{2}}}{3} + B\right)\log\left(\epsilon^{\frac{1}{2}} - B^{\frac{1}{2}} + \zeta^{\frac{1}{2}} + M_5\right), $$

$$ M_5 = \left[ B + 2B^{\frac{1}{2}} \left( \epsilon^{\frac{1}{2}} - \zeta^{\frac{1}{2}} \right) + \left( \epsilon^{\frac{1}{2}} + \zeta^{\frac{1}{2}} \right)^2 \right]^{\frac{1}{2}}, $$

$$ a(y) = \zeta^{\frac{1}{2}}(1 - y)^{\frac{1}{2}}, \quad b = 3\zeta^{\frac{3}{2}}\left(\epsilon^{\frac{1}{2}} + \frac{5}{3}B^{\frac{1}{2}}\right), \quad c = \left(\epsilon^{\frac{1}{2}} + B^{\frac{1}{2}}\right)^2. $$
Now, we normalize the velocity of the free surface so that $u(1) = 1$, obtaining the implicit relation between Re and $\xi$, i.e.,

$$ \mathcal{F}(\text{Re}, \xi) = 1, \qquad (23) $$

where

$$ \mathcal{F}(\text{Re}, \xi) = -\frac{4}{\xi^{\frac{3}{2}} \text{Re}^{\frac{1}{3}}} (D_1(\text{Re}, \xi) \log(D_2(\text{Re}, \xi)) + (D_3(\text{Re}, \xi))), \quad (24) $$

and

$$
\begin{align}
D_1(\text{Re}, \xi) &= -X^{\frac{1}{2}} \epsilon^{\frac{3}{2}} \xi^{\frac{5}{6}} \text{Re}^{\frac{1}{6}} + \frac{3}{2} X \epsilon \xi^{\frac{7}{6}}, \nonumber \\
D_2(\text{Re}, \xi) &= \epsilon^{\frac{1}{2}} + \xi^{\frac{1}{2}} - \frac{X^{\frac{1}{2}} \xi^{\frac{1}{3}}}{\text{Re}^{\frac{1}{6}}} + A(\text{Re}, \xi), \nonumber \\
D_3(\text{Re}, \xi) &= D_4(\text{Re}, \xi) + D_5(\text{B}, \xi) + D_6(\text{Re}, \xi), \nonumber \\
D_4(\text{Re}, \xi) &= \frac{A(\text{Re}, \xi)}{16} \left[ \frac{1}{3} \frac{X^{\frac{1}{2}} \xi^{\frac{1}{3}}}{\text{Re}^{\frac{1}{6}}} P_1(\text{Re}, \xi) + P_2(\text{Re}, \xi) \right], \nonumber \\
D_5(\text{Re}, \xi) &= P_3(\text{Re}, \xi) \log(2\epsilon^{\frac{1}{2}})\epsilon, \nonumber \\
D_6(\text{Re}, \xi) &= \frac{P_5(\text{Re}, \xi)}{16} \left( \frac{X^{\frac{1}{2}} \xi^{\frac{1}{3}}}{\text{Re}^{\frac{1}{6}}} P_4(\text{Re}, \xi) + \epsilon^{\frac{3}{2}} \text{Re}^{\frac{1}{3}} \xi^{\frac{1}{2}} + \frac{11}{3} \epsilon^{\frac{1}{2}} X \xi^{\frac{7}{6}} \right) + P_6(\text{Re}, \xi), \nonumber \\
P_1(\text{Re}, \xi) &= -X \xi^{\frac{7}{6}} + \text{Re}^{\frac{1}{3}} \left[ 5 \xi^{\frac{3}{2}} + \epsilon^{\frac{1}{2}} (47 \epsilon^{\frac{1}{2}} \xi^{\frac{1}{2}} - 8 \xi) \right], \tag{25} \\
P_2(\text{Re}, \xi) &= -\frac{11}{3} \epsilon^{\frac{1}{2}} X \xi^{\frac{7}{6}} + \epsilon^{\frac{1}{2}} \text{Re}^{\frac{1}{3}} \xi^{\frac{3}{2}} - \frac{1}{3} X \xi^{\frac{5}{3}} - \text{Re}^{\frac{1}{3}} (\epsilon^{\frac{3}{2}} \xi^{\frac{1}{2}} - \epsilon \xi + \xi^2), \nonumber \\
P_3(\text{Re}, \xi) &= X^{\frac{1}{2}} \epsilon^{\frac{1}{2}} \xi^{\frac{5}{6}} \text{Re}^{\frac{1}{6}} - \frac{3}{2} X \xi^{\frac{7}{6}}, \nonumber \\
P_4(\text{Re}, \xi) &= \frac{1}{3} (-47 \epsilon \xi^{\frac{1}{2}} \text{Re}^{\frac{1}{3}} + X \xi^{\frac{7}{6}}), \nonumber \\
P_5(\text{Re}, \xi) &= \epsilon^{\frac{1}{2}} + \frac{X^{\frac{1}{2}} \xi^{\frac{1}{3}}}{\text{Re}^{\frac{1}{6}}}, \nonumber \\
P_6(\text{Re}, \xi) &= -\frac{X^{\frac{1}{2}} \xi^{\frac{3}{5}} \text{Re}^{\frac{1}{6}}}{4} \left( \xi^{\frac{3}{2}} \epsilon^{\frac{1}{2}} - \frac{2}{3} \xi^2 \right) - \frac{1}{8} \epsilon \xi^{\frac{3}{2}} \text{Re}^{\frac{1}{3}} - \frac{1}{16} \xi^5 \text{Re}^{\frac{1}{3}} - \frac{1}{8} X \xi^{\frac{13}{6}}, \nonumber \\
A(\text{Re}, \xi) &= \left[ 2 \frac{X^{\frac{1}{2}} \xi^{\frac{1}{3}}}{\text{Re}^{\frac{1}{6}}} (\epsilon^{\frac{1}{2}} - \xi^{\frac{1}{2}}) + \frac{X \xi^{\frac{2}{3}}}{\text{Re}^{\frac{1}{3}}} + (\epsilon^{\frac{1}{2}} + \xi^{\frac{1}{2}})^2 \right]^{\frac{1}{2}}. \nonumber
\end{align}
$$

As expected, we obtain a one-to-one relation between Re and $\xi$ and we denote by $\xi$ the unique solution to (23) such that $u(1) = 1$. The plot $\mathcal{F}(\text{Re}, \xi)$ is displayed in Figure 3. Moreover, we recall that B is expressed in terms of Re through (15), where now $\xi$ is the solution of (23). Relation (23) defines as $\mathcal{F}(\text{Re}, \xi) = 1$ which, as expected, is a one-to-one relation between Re and $\xi$. The plot of $\mathcal{F}(\text{Re}, \xi)$ is displayed in Figure 4, which highlights that, for given $\lambda_1, \lambda_2$, and $\theta$, there exists a unique $(\text{Re}, \xi)$ fulfilling (23). Consequently, recalling that Equation (23) derives from the normalization of $u$, for any Re we obtain a unique value of $\xi$, which we denoted as $\xi$, such that $u(1) = 1$.
<figure><img src="image_6.png" /><figcaption>Figure 3. Plot of F(Re, ξ) for λ₁ = 0.1 and θ = 5° for the regularized Casson model (given by (23)).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 4. Plot of u(y) given by (21) (empty circles) and (26) (solid line) with λ₁ = 0.1 and θ = 1° for different values of ε. It is worth noting that the “exact” profile of u(y), given by (26), can be retrieved from the regularized one, given by (21). In fact, the two profiles become very similar when ε < 10⁻².</figcaption></figure>

We remark that for ε → 0 we retrieve the Casson flow whose normalized velocity field is given by

$$
u(y) = \begin{cases} -\frac{\zeta y^2}{2} + (\zeta\sigma + 2B)y + \frac{4\sqrt{B}[(\zeta\sigma+B)^{3/2}+(\zeta\sigma-y\zeta+B)^{3/2}]}{3\zeta}, & \text{if } 0 \le y \le \sigma, \\ 1, & \text{if } \sigma \le y \le 1, \end{cases} \quad (26)
$$
with

$$ \sigma = 1 - \frac{\tau_0^*}{\rho^* g^* H^* \sin \theta} = 1 - \frac{B}{\xi}, \qquad (27) $$

is the flat yield surface [10,16]. The velocity field is normalized so that $u(\sigma) = 1$, i.e.,

$$ -\frac{\zeta\sigma^2}{2} + (\zeta\sigma + 2B)\sigma + \frac{4\left[-\sqrt{B}(\zeta\sigma + B)^{3/2} + B^2\right]}{3\zeta} = 1, \qquad (28) $$

that leads, using (15) and (27), to the following implicit relation between $\zeta$ and Re

$$ \mathcal{F}(\text{Re}, \zeta) = 1, \qquad (29) $$

where

$$ \mathcal{F}(\text{Re}, \zeta) = -\frac{X^2}{6} \zeta^{\frac{1}{3}} \text{Re}^{-\frac{2}{3}} - \frac{4}{3} \sqrt{\zeta} \sqrt{X \zeta^{\frac{2}{3}} \text{Re}^{-\frac{1}{3}}} + X \zeta^{\frac{2}{3}} \text{Re}^{-\frac{1}{3}} + \frac{\zeta}{2}. \qquad (30) $$

The regularized and the “exact” profiles of $u(y)$ are given by (21) and (26), respectively. They become very similar when $\epsilon < 10^{-2}$ as shown in Figure 4.

It is worth noting that in case $B = 0$, i.e., Newtonian flow, Equation (23), for $\epsilon \to 0$, and Equation (29) simply reduce to $\mathcal{F} = \zeta - 2 = 0$, whose trivial solution is $\hat{\zeta} = 2$.

The equations governing the regularized and exact Bingham flow have been reported in [10]. In this paper, we have adopted the same notations as in [10], thus the comparison between the two models can be performed easily.

## 4. Linear Stability Furthermore, Long-Wave Approximation

In this section, we briefly recall the main characteristic of the linear stability analysis as reported in [10,16,51] and we refer the readers to [10,16] for more details on the derivation of the formulas here summarized.

We consider the basic flow consisting of $h(x, t) = h_b$, with $h_b = 1$, $v_b = u_b(y)\mathbf{i}$ with $u_b$ given by (21), and, $p = p_b(y)$ where, recalling (20), $p_b(y) = \hat{\zeta} \cot \theta(1 - y)$. Then, we perturb the basic flow superimposing small disturbances, 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} \qquad (31) $$

and

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

where $\alpha \in \mathbb{R}$ is the wave number, $c \in \mathbb{C}$ is the complex wave speed and the notation $\hat{\cdot}$ represents the infinitesimal disturbance. We write the velocity field in terms of the stream function, i.e.,

$$ \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)}, \qquad (33) $$

where, here and in the sequel, $(\cdot)'$ denotes the differentiation with regard to $y$. Defining by $\Re(c)$ and $\Im(c)$ the real and imaginary part of $c$, we recall that $\Im(c)$ gives the growth/attenuation factor of the $\alpha^{th}$ mode. Hence, the basic flow $h_b$, $v_b$, $p_b(y)$ is unstable when the parameters involved in the problem, namely Re, $\lambda_1$, $\lambda_2$ and $\theta$, are selected so that $\Im(c) > 0$. The transition between the two regimes is identified by the so-called marginal or neutral curve, i.e., the set of Re, $\lambda_1$, $\lambda_2$ and $\theta$ at which $\Im(c) = 0$.
Then, we consider disturbances of long wavelength $2\pi/\alpha \gg 1$, i.e., $\alpha \ll 1$, expanding $\phi$ and $c$ in powers of $\alpha$ up to the first order in $\alpha$, namely

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

where $(\phi_0, c_0)$ and $(\phi_1, c_1)$ solve

$$ \left\{ \begin{aligned} (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{aligned} \right. \tag{35} $$

and

$$ \left\{ \begin{aligned} & (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{aligned} \right. \tag{36} $$

with

$$ 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], $$

and $r$ given by (13). In particular, we have that $\Im(c_0) = 0$, while $\Re(c_1) = 0$, thus

$$ \phi(y) e^{i\alpha(x-ct)} = \phi(y) \underbrace{e^{i\alpha(x-c_0t)}}_{\text{travelling wave growth/attenuation}} \underbrace{e^{a^2\Im(c_1)t}}_{\text{attenuation}}. \tag{37} $$

In particular, we can find the critical value of Re, denoted as $Re_c$, such that

$$ \Im(c) = \Im(c_1(\operatorname{Re}_c, \lambda_1, \lambda_2, \theta)) = 0, \tag{38} $$

by prescribing the material characteristics and the tilt angle (i.e., $\lambda_1$, $\lambda_2$, and $\theta$). Hence, for $Re < Re_c$ the $a^{th}$ mode is stable, while instability arises when $Re > Re_c$, since $\Im(c_1) < 0$ for $Re < Re_c$ and vice versa. Moreover, we eventually remark that the identification of $Re_c$ means, from the practical point of view, the identification of a critical discharge, $Q_c^*$, above which the flow becomes unstable.
## 5. Results

The critical value of the Reynolds number, $Re_c$, is computed by solving the system given by the system of algebraic Equations (23) and (38) with MATLAB® 2022a, using the function FSOLVE.

Figure 5 shows the variation of $Re_c$ with respect to the tilt angle $\theta$ for different value of the material parameter $\lambda_1$, when $\epsilon = 0.01$ by considering the regularized Bingham and Casson model. Similar to [10], at a given $\theta$, $Re_c$ decreases for increasing values of $\lambda_1$ and coincides with $5/4 \cot \theta$, when $\lambda_1 = 0$, i.e., $B = 0$. However, recalling the proportionality relation (16) between $\lambda_1$ and $\tau_0^*$, we have that (see Figure 5A) the yield stress destabilizes the flow when this is modelled using a regularized Bingham flow, while (see Figure 5B) the regularized Casson flow is more stable than the Newtonian flow (i.e., the yield stress has a stabilizing effect on the flow when the material is modelled as a regularized Casson fluid). Coherently, in the case of regularized Casson fluid, $Re_c$ is an increasing function of $\lambda_1$, namely an increase in the yield stress $\tau_0^*$ leads to a flow stabilization (Figure 6C,D). For the regularized Bingham, we have an opposite behaviour. Indeed, Figure 6A,B highlight that $Re_c$ decreases as $\lambda_1$ (i.e., $\tau_0^*$) increases. Moreover, it is worth noting that, as physically expected, an increase of $\theta$ leads to flow destabilization in both cases.

In Table 1, we report the values of $Re_c$ and $c_0$ when $\theta = 5^\circ$ for various values of $\lambda_1$ when a regularized Bingham and Casson models are considered. We notice that $c_0 = 2$ when $\lambda_1 = 0$, i.e., when the flow is Newtonian, as in [10]. Again, coherently with the results obtained in [10], as $\lambda_1$ increases the superficial wave speed increases also for the regularized Casson fluid.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Evolution of the critical Reynolds number, $Re_c$, with respect to the tilt angle, $\theta$, with $\epsilon = 0.01$ for different values of $\lambda_1$ in the case of the flow modelled as a regularized Bingham (A) and regularized Casson (B) fluid. The theoretical Newtonian flow (i.e., $5/4 \cot \theta$) is given by the red circles. The continuous line is the Newtonian flow computed by our code. We emphasize that the theoretical curve and the computed one coincide.</figcaption></figure>

<table><caption>**Table 1.** Values of $Re_c$ and $c_0$ for given values of $\lambda_1$ with $\theta = 5^\circ$ and $\epsilon = 0.01$ for both the regularized Bingham and Casson models.</caption><thead><tr><th rowspan="2">λ<sub>1</sub></th><th colspan="2">Re<sub>c</sub></th><th colspan="2">c<sub>0</sub></th></tr><tr><th>Reg. Bingham</th><th>Reg. Casson</th><th>Reg. Bingham</th><th>Reg. Casson</th></tr></thead><tbody><tr><td>0</td><td>14.29</td><td>14.29</td><td>2</td><td>2</td></tr><tr><td>0.01</td><td>14.08</td><td>15.61</td><td>2.03</td><td>2.15</td></tr><tr><td>0.1</td><td>12.45</td><td>19.08</td><td>2.34</td><td>2.52</td></tr><tr><td>0.5</td><td>9.05</td><td>28.13</td><td>4.31</td><td>3.30</td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Plot of the critical Reynolds number, $Re_c$, as a function of $\lambda_1$ with $\epsilon = 0.01$ for different values of $\theta$ for a regularized Bingham (A,B) and regularized Casson (C,D) fluid. The case of Newtonian flow corresponds to $\lambda_1 = 0$, i.e., B = 0.</figcaption></figure>

**6. Conclusions**

In this paper, the stability analysis of a free surface regularized Casson flow down an incline has been theoretically investigated and compared to the one obtained for the regularized Bingham flow. In both cases the benchmark represented by the Newtonian case (i.e., when $\lambda_1 = 0$, thus B = 0) has been recovered. Our results show that for a regularized Casson fluid $Re_c$ increases with increasing values of the “material” parameter $\lambda_1$ (that is proportional to $\tau_0^*$ see (16)), while for a regularized Bingham fluid $Re_c$ decreases when $\lambda_1$ (i.e., $\tau_0^*$) increases [10,16]. Therefore, our findings (obtained within the long-wave approximation method) show that the flow of the regularized Casson fluid is stabilized by increasing yield stress $\tau_0^*$ contrary to what happens with the regularized Bingham.

The stability analysis of the exact Bingham model, investigated in [16], shows that the flow down an incline is unconditionally stable for every Reynolds number. Therefore, our results are unexpected, highlighting that, although in a regularized formulation, models belonging to the same class of viscoplastic fluids can have stability characteristics completely different. Although we are not aware of any studies on the stability of the Casson fluid flowing down an incline, we suppose that the results of [12] can also be extended to this case. Therefore, we show that (as for the regularized Bingham fluid [10]) the regularized Casson can have stability properties that are different from the classic Casson flow. It is worth remarking that the our study has been developed by applying the long-wave approximation to the flow of regularized Bingham and Casson fluids down an incline.

We have in fact shown that the flow along an incline of a Casson-type material becomes increasingly stable as the yield stress increases. Exactly the opposite behaviour occurs with the Bingham fluid. Therefore, this feature can be used, from an experimental point of view,
to highlight the difference between the two rheological models. We believe that our results
can pave the way to experimental studies on the flow down an incline.

**Author Contributions:** Conceptualization, B.C., L.I.P., L.F. and A.F.; methodology, B.C., L.I.P., L.F. and A.F.; software, B.C.; formal analysis, L.F. and A.F.; investigation, B.C., L.I.P., L.F. and A.F.; writing—original draft preparation, B.C., L.I.P., L.F. and A.F.; writing—review and editing, B.C., L.F. and A.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partially supported by GNFM of Italian INDAM.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflicts of interest.

## References

1. Bingham, E.C. *Fluidity and Plasticity*; McGraw-Hill: New York, NY, USA, 1922.

2. Herschel, W.H.; Bulkley, R. Konsistenzmessungen von gummi-benzollosungen. *Kolloid-Z.* **1926**, *39*, 291–300. [CrossRef]

3. Casson, N. A Flow Equation for Pigment-oil Suspensions of the Printing Ink Type. In *Rheology of Disperse Systems*; Mill, C.C., Ed.; Pergamon Press: Oxford, UK, 1959; pp. 84–104.

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

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

6. Liu, J.; Paul, J.D.; Gollub, J.P. Measurements of the primary instabilities of film flows. *J. Fluid Mech.* **1993**, *250*, 69–101. [CrossRef]

7. 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. [CrossRef]

8. 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*, R1. [CrossRef]

9. Balmforth, N.J.; Liu, J.J. Roll waves in mud. *J. Fluid Mech.* **2004**, *519*, 33–54. [CrossRef]

10. Calusi, B.; Farina, A.; Fusi, L.; Rosso, F. Long-wave instability of a regularized Bingham flow down an incline. *Phys. Fluids* **2022**, *34*, 054111. [CrossRef]

11. 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. [CrossRef]

12. Falsaperla, P.; Giacobbe, A.; Mulone, G. Stability of the plane Bingham-Poiseuille flow in an inclined channel. *Fluids* **2020**, *5*, 141. [CrossRef]

13. Fernandez-Nieto, E.D.; Noble, P.; Vila, J.-P. Shallow water equations for non-Newtonian fluids. *J. Non-Newtonian Fluid Mech.* **2010**, *165*, 712–732. [CrossRef]

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

15. Fusi, L. Channel flow of viscoplastic fluids with pressure-dependent rheological parameters. *Phys. Fluids* **2018**, *30*, 073102. [CrossRef]

16. Fusi, L.; Calusi, B.; Farina, A.; Rosso, F. Stability of laminar viscoplastic flows down an inclined open channel. *Eur. J. Mech.-B/Fluid* **2022**, *95*, 137–147. [CrossRef]

17. 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. [CrossRef]

18. 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. [CrossRef]

19. 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. [CrossRef]

20. Hu, T.; Fu, Q.F.; Xing, Y.; Yang, L.J.; Xie, L. Stability of a thin viscoelastic film falling down an inclined plane. *Phys. Rev. Fluids* **2021**, *6*, 083902. [CrossRef]

21. Métivier, C.; Nouar, C. Stability of a Rayleigh-Bénard Poiseuille flow for yield stress fluids—Comparison between Bingham and regularized models. *Int. J. Non-Linear Mech.* **2011**, *46*, 1205–1212. [CrossRef]

22. 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. [CrossRef]

23. 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. [CrossRef]

24. 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. [CrossRef]

25. Mogilevskiy, E. Stability of a non-Newtonian falling film due to three-dimensional disturbances. *Phys. Fluids* **2020**, *32*, 073101. [CrossRef]
26. 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.
[CrossRef]

27. Noble, P.; Vila, J.-P. Thin power-law film flow down an inclined plane: Consistent shallow-water models and stability under large-scale perturbations. *J. Fluid Mech.* **2013**, *735*, 29–60. [CrossRef]

28. 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. [CrossRef]

29. Nsom, B.; Ramifidioa, L.; Latrache, N.; Ghaemizadeh, F. Linear stability of shear-thinning fluid down an inclined plane. *J. Mol. Liquids* **2019**, *277*, 1036–1046. [CrossRef]

30. Pascal, J.P.; D'Alessio, S.J.D. Instability of power-law fluid flows down an incline subjected to wind stress. *Appl. Math. Model.* **2007**, *31*, 1229–1248. [CrossRef]

31. 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. [CrossRef]

32. Ruyer-Quil, C.; Chakraborty, S.; Dandapat, B.S. Wavy regime of a powerlaw film flow. *J. Fluid Mech.* **2012**, *692*, 220–256. [CrossRef]

33. Merrill, E.W.; Margetts, W.G.; Cokelet, G.C.; Gilliland, E.R. The Casson equation and rheology of the blood near shear zero. In *Proceedings Fourth International Congress on Rheology*; Copley, A.L., Ed.; Interscience: New York, NY, USA, 1965; Part 4, pp. 135–143.

34. Farina, A.; Fasano, A.; Rosso, F. Mathematical models for some aspects of blood microcirculation. *Symmetry* **2021**, *13*, 1020.
[CrossRef]

35. Fasano, A.; Sequeira, A. *Hemomath: The Mathematics of Blood*; Springer: Berlin/Heidelberg, Germany, 2017.

36. Charakopoulos, A.; Karakasidis, T.; Sarris, I. Analysis of magnetohydrodynamic channel flow through complex network analysis. *Chaos* **2021**, *31*, 043123. [CrossRef] [PubMed]

37. Fusi, L. Lubrication flow of a generalized Casson fluid with pressure-dependent rheological parameters. *J. Non-Newton. Fluid Mech.* **2019**, *274*, 104199. [CrossRef]

38. 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. [CrossRef]

39. Reddy, M.G.; Kumara, B.C.P.; Makinde, O.D. Cross Diffusion Impacts on Hydromagnetic Radiative Peristaltic Carreau-Casson Nanofluids Flow in an Irregular Channel. *Defect Diffus. Forum* **2017**, *377*, 62–83. [CrossRef]

40. Sofos, F.; Karakasidis, T.; Spetsiotis, D. Molecular dynamics simulations of ion separation in nano-channel water flows using an electric field. *Mol. Simul.* **2019**, *45*, 1395–1402. [CrossRef]

41. Song, L.; Huo, X.; Zhang, L.; Xie, Y.; Yang, M. Fast estimation on the pressure of detonation products of cyclotetramethylene tetranitramine through molecular dynamics simulations. *Int. J. Mod. Phys. B* **2021**, *35*, 2150106. [CrossRef]

42. Astarita, G. Letter to the editor: The engineering reality of the yield stress. *J. Rheol.* **1990**, *34*, 275–277. [CrossRef]

43. Barnes, H.A. The yield stress—A review or ἀνταρει'—Everything flows? *J. Non-Newton. Fluid Mech.* **1999**, *81*, 133–178.
[CrossRef]

44. Barnes, H.A.; Walters, K. The yield stress myth? *Rheol. Acta* **1985**, *24*, 323–326. [CrossRef]

45. Frigaard, I.; Nouar, C. On the usage of viscosity regularisation methods for visco-plastic fluid flow computation. *J. Non-Newton. Fluid Mech.* **2005**, *127*, 1–26. [CrossRef]

46. Frigaard, I.A.; Paso, K.G.; de Souza Mendes, P.R. Bingham's model in the oil and gas industry. *Rheol. Acta* **2017**, *56*, 259–282.
[CrossRef]

47. Fusi, L.; Farina, A.; Rosso, F. On the mathematical paradoxes for the flow of a viscoplastic film down an inclined surface. *Int. J. Non-Linear Mech.* **2014**, *58*, 139–150. [CrossRef]

48. Fusi, L.; Farina, A.; Rosso, F. Roscani, S. Pressure driven lubrication flow of a Bingham fluid in a channel: A novel approach. *J. Non-Newtonian Fluid Mech.* **2015**, *221*, 66–75. [CrossRef]

49. Bercovier, M.; Engleman, M. A finite-element method for incompressible non-Newtonian flows. *J. Comput. Phys.* **1980**, *36*, 313–326.
[CrossRef]

50. Papanastasiou, T.C. Flows of materials with yield. *J. Rheol.* **1987**, *31*, 385–404. [CrossRef]

51. Pascal, J.P. Linear stability of fluid flow down a porous inclined plane. *J. Phys. D Appl. Phys.* **1999**, *32*, 417. [CrossRef]

52. Allouche, M.; Frigaard, I.A.; Sona, G. Static wall layers in the displacement of two visco-plastic fluids in a plane channel. *J. Fluid Mech.* **2000**, *424*, 243–277. [CrossRef]