# On a casson fluid motion: Nonuniform width symmetric channel and peristaltic flows

Simone Guadagni, Liviu Iulian Palade, Lorenzo Fusi, Angiolo Farina

► To cite this version:

Simone Guadagni, Liviu Iulian Palade, Lorenzo Fusi, Angiolo Farina. On a casson fluid motion: Nonuniform width symmetric channel and peristaltic flows. Fluids, 2021, 6 (10), pp.356. 10.3390/fluids6100356. [https://doi.org/10.3390/fluids6100356](https://doi.org/10.3390/fluids6100356)

HAL Id: hal-04134242

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

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

# On a Casson Fluid Motion: Nonuniform Width Symmetric Channel and Peristaltic Flows

Simone Guadagni<fnref n="1" />, Liviu Iulian Palade<fnref n="2" />, Lorenzo Fusi<fnref n="1" />* and Angiolo Farina<fnref n="1" />,†

**Abstract:** Widely used for modeling biological fluids flows—in particular, blood vessel flows— a Casson flow is studied in a symmetric channel for which the aspect ratio enables one to use the lubrication approximation. Two flow driving conditions are prescribed: inlet–outlet pressure difference and peristaltic oscillations of the vessel walls. In both cases, starting from mass and momentum balance and using lubrication approximation, we investigate the conditions to be imposed on the driving mechanisms so that the inner plug does not come in touch with the walls. The study of the peristaltic flow is of great importance in view of its applications in physiology (including microcirculation applications).

**Keywords:** lubrication; non-Newtonian fluid; free boundary problems

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

check for
updates

**Citation:** Guadagni, S.; Palade, L.I.; Fusi, L.; Farina, A. On a Casson Fluid Motion: Nonuniform Width Symmetric Channel and Peristaltic Flows. *Fluids* **2021**, *6*, 356. https://doi.org/10.3390/fluids6100356

Academic Editors: Antonio Santamaría and Mehrdad Massoudi

Received: 28 August 2021
Accepted: 29 September 2021
Published: 8 October 2021

**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:** © 2021 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

During the past several decades, slurries, pastes, rude oil, drilling mud, mineral slurries, suspensions, paint, and materials such as meat extract, frequently encountered in industrial problems, have received increasing attention. The rheological behavior of these materials is usually considered to be viscoplastic, i.e., showing a critical value of stress below which the material does not flow. Such a critical value is usually referred to as yield stress. Viscoplastic materials are also called Bingham plastics, after Bingham [1], who was the first to describe several types of paint in this way in 1919. The models used for these materials included also the Herschel-Bulkley model [2], the Casson model [3] and the Heinz-Casson model [4]. A detailed review of viscoplastic fluids can be found in articles by Bird, Dai and Yarusso [5], Mitsoulis [6], Huilgol [7], and Frigaard and Nouar [8].

All of these models have various practical applications. The Herschel-Bulkley models are used for muds, foams, ceramics, and slurries. The Casson model is widely used to model blood flow [9,10], as it describes blood flow in low shear regions quite well since it captures the nonlinear dependence of viscosity on shear rate. Indeed, the Casson fluid undergoes no deformation until the shear stress is below the critical threshold; above such a threshold, it displays a shear thinning behavior. Merrill et al. [11] found that at shear rates in the range 0.1–1.0 s⁻<fnref n="1" />, the Casson constitutive equation fitted the experimental data well. However, we remark that the existence of a yield stress for blood and its use as a material parameter is currently still a controversial issue due to the sensitivity of yield stress measurements in connection to factors that are hardly controllable.

Most of the efforts in the theoretical analyses concern the extent and the shape of yielded/unyielded regions, which are the main feature of viscoplastic materials. A possible approach is the one introduced in [12,13], where the equation of motion of the unyielded part is written in an integral form. According to this method, originally developed by Safronchik [14], the unyielded region is treated as a rigid body of variable mass whose

<footnote n="1">¹ Dipartimento di Matematica e Informatica "U. Dini", Università degli Studi di Firenze, Viale Morgagni 67/A, 50134 Firenze, Italy; simone.guadagni@unifi.it (S.G.); lorenzo.fusi@unifi.it (L.F.)</footnote>

<footnote n="2">² CNRS, Institut Camille Jordan UMR 5208, INSA-Lyon, Université de Lyon, Pôle de Mathématiques, Bât. Leonard de Vinci No. 401, 21 Avenue Jean Capelle, 69621 Villeurbanne, France; liviu-iulian.palade@insa-lyon.fr</footnote>

<footnote>* Correspondence: angiolo.farina@unifi.it</footnote>
dynamics is governed by the cardinal equations. We remark that the yield surface can be determined using other methods, such as the ones illustrated in [8,15–17].

In this paper, we focus on the Casson model and analyze the flow in a symmetric channel with “small” aspect ratio so that the lubrication approximation can be safely used. Lubrication approximation is widely used in a large variety of contexts (see e.g., [18–22]). In particular, the geometrical setting of our study is similar to that studied in [21], where a similar problem is treated for a Newtonian fluid. We analyze the dynamics considering two conditions driving the flow: prescribed inlet–outlet pressure difference and peristaltic motion. The case of prescribed inlet discharge has already been analyzed in [23]. We also investigate the conditions ensuring that the inner plug does not come in touch with the walls. Peristaltic flow (for which walls are set in motion by a traveling velocity profile) is studied because of their great importance in understating artery and vein physiology [9,10]. The peristaltic motion of Carreau and Casson nanofluid under the action of a magnetic field and in nonisothermal conditions has been investigated in [24]. The magnetohydrodynamic flow was analyzed in [25], and the dynamics in nanochannels were analyzed in [26,27].

The main aim of this paper is to formulate a sound mathematical model for the motion of a Casson fluid in a vessel with small aspect ratio. As remarked, these flows characterize several applications that span from physiology to geology. Actually, the novelty of this paper is the rigorous application of the momentum equation to the unyielded phase. In the viscoplastic theory, the inner core is in fact a rigid body whose stress cannot be determined. Therefore, in order to get a reliable model, one has to correctly describe its motion equation. This is the key point of our approach. Indeed, the mathematical problem we have obtained is well posed, and this allows us to obtain useful information on the flow. The manuscript is organized as follows. In Section 2, we introduce the mathematical model. The dimensionless form is obtained in Section 3. The asymptotic expansions and the corresponding numerical results are given in Section 5. Finally, the peristaltic flow is investigated in Section 6. Finally, we underline that the study on the peristaltic motion of a Casson fluid can have repercussions in the medical field—in particular, in hemorheology.

## 2. Mathematical Model

We write the Cauchy stress tensor as $T = -pI + \tau$, where $\tau$ is its deviatoric traceless part. The constitutive equation for a Casson fluid is given by

$$ \left\{ \begin{array}{ll} \tau = \left( \sqrt{\mu} + \frac{\sqrt{\tau_0}}{\sqrt{|\dot{\gamma}|}} \right)^2 \dot{\gamma}, & \text{if } |\tau| \ge \tau_0, \\ \\ \dot{\gamma} = 0, & \text{if } |\tau| \le \tau_0, \end{array} \right. \qquad (1) $$

where

$$ \dot{\gamma} = \nabla \mathbf{v} + \nabla \mathbf{v}^T, $$

with **v** velocity field and where

$$ |\dot{\gamma}| = \sqrt{\frac{1}{2} \dot{\gamma} \cdot \dot{\gamma}}, \qquad |\tau| = \sqrt{\frac{1}{2} \tau \cdot \tau}. $$

We investigate the 2D flow in a channel such as the one depicted in Figure 1. The x axis coincides with the channel symmetry line and $y = \pm h(x)$ are the channel walls.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Flow in a symmetric channel with curved walls. The yielded region $|\tau| \ge \tau_o$ and the unyielded one $|\tau| < \tau_o$ are separated by a sharp interface $y = \pm\sigma(x, t)$. In this schematic drawing, the channel width is exaggerated for visualization purposes. We are indeed within the framework of lubrication theory such that the channel aspect ratio $\epsilon$, defined by (10), i.e., the ratio between the typical channel width $H$ and the channel length $L$, is small. Then, following the standard asymptotic analysis techniques, we look for a solution in which the main problem variables are expressed as power series of $\epsilon$. By equating the terms with the same power of $\epsilon$, we obtain a sequence of problems to the various approximation orders. In this study, we limit our analysis to the leading order, namely to the $\epsilon^0$ terms.</figcaption></figure>

The velocity field is given by

$$v(x, y, t) = u(x, y, t)e_1 + v(x, y, t)e_2,$$

and that the fluid is incompressible, i.e., tr $\dot{\gamma} = 0$ (isochoric or volume preserving flow). We assume that the yielded part of the flow $|\tau| \ge \tau_o$ and the unyielded part $|\tau| \le \tau_o$ are separated by the sharp interfaces $y = \pm\sigma(x, t)$, which actually are two free boundaries since their evolution is unknown. Because of symmetry, we may confine our analysis only to the upper part of the channel.

The motion of the unyielded region (schematically depicted in Figure 2)

$$\Omega_t = \{(x, y) : 0 \le x \le L, 0 \le y \le \sigma(x, t)\},$$

obeys to the momentum equation that we write in an integral form (see [28])

$$\int_{\Omega_t} \frac{\partial}{\partial t}(\rho \mathbf{v})dV + \int_{\partial \Omega_t} \rho \mathbf{v}(\mathbf{v} \cdot \mathbf{n})dS = \int_{\partial \Omega_t} \mathbf{\tau} \mathbf{n} dS, \quad (2)$$

where $\mathbf{n}$ is the normal to $\Omega_t$ pointing outward. In the yielded region $\sigma \le y \le h$, the governing equations of the system

$$\left\{
\begin{array}{l}
\displaystyle \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0, \\
\\
\displaystyle \rho \left( \frac{\partial u}{\partial t} + \frac{\partial u}{\partial x} u + \frac{\partial u}{\partial y} v \right) = - \frac{\partial p}{\partial x} + \frac{\partial \tau_{11}}{\partial x} + \frac{\partial \tau_{12}}{\partial y}, \\
\\
\displaystyle \rho \left( \frac{\partial v}{\partial t} + \frac{\partial v}{\partial x} u + \frac{\partial v}{\partial y} v \right) = - \frac{\partial p}{\partial y} + \frac{\partial \tau_{12}}{\partial x} + \frac{\partial \tau_{22}}{\partial y},
\end{array}
\right.
\quad (3)$$

where body forces have been neglected.
We remark that Equation (2) is nothing but the first cordinal equation for the rigid core. Because of symmetry, no rotations occur and so (2) is sufficient to determine the translation of the unyielded region.

In the unyielded part, because of symmetry, the motion is a pure translation and the velocity is given by

$$v = u_c(t)\mathbf{e}_1,$$

<figure><img src="image_5.png" /><figcaption>Figure 2. The unyielded domain $\Omega_t$. Note that $\Omega_t$ may vary on time.</figcaption></figure>

To write (2) componentwise, we look at Figure 2 where the boundary $\partial\Omega_t$ has been divided into four components

$$\partial\Omega_t = \mathcal{F}_1 \cup \mathcal{F}_2 \cup \mathcal{F}_3 \cup \mathcal{F}_4.$$

The first component of Equation (2) is given by

$$\frac{\partial}{\partial t}(\rho u_c) \int_{\Omega_t} dV + \rho u_c^2 \int_{\partial\Omega_t} (\mathbf{e}_1 \cdot \mathbf{n})dS = \int_{\partial\Omega_t} (\tau \mathbf{n} \cdot \mathbf{e}_1)dS, \quad (4)$$

while the second is given by

$$0 = \int_{\partial\Omega_t} (\tau \mathbf{n} \cdot \mathbf{e}_2)dS. \quad (5)$$

The stress tensor is given by

$$\tau = \begin{pmatrix} -p + \tau_{11} & \tau_{12} \\ \tau_{12} & -p + \tau_{22} \end{pmatrix},$$

on the surfaces $\mathcal{F}_1$ and $\mathcal{F}_3$, the only nonzero components of the applied stress are the normal ones, namely

$$\tau|_{x=0} = \begin{pmatrix} -p_{in} & 0 \\ 0 & -p_{in} \end{pmatrix}, \quad \tau|_{x=L} = \begin{pmatrix} -p_{out} & 0 \\ 0 & -p_{out} \end{pmatrix},$$

which implies that no torque is applied to the rigid core. Under these hypotheses, setting $\sigma_{in} = \sigma(0, t)$, $\sigma_{out} = \sigma(L, t)$ and recalling

$$dS = \sqrt{\sigma_x^2 + 1} dx,$$
we find that (5) is automatically fulfilled, while (4) becomes

$$2\frac{\partial}{\partial t}(\rho u_c) \int_{0}^{L} \sigma dx = 2(p_{in}\sigma_{in} - p_{out}\sigma_{out}) + \\ +2 \int_{0}^{L} (\sigma_x p - \sigma_x \tau_{11} + \tau_{12}) \Big|_{\sigma} dx.$$

Since $u_c$ depends only on the $t$ variable, integrating by parts, we find

$$\rho \frac{du_c}{dt} \int_{0}^{L} \sigma(x,t) dx = \int_{0}^{L} \left[ -\frac{\partial p}{\partial x}\sigma - \frac{\partial \sigma}{\partial x}\tau_{11} + \tau_{12} \right] \Big|_{\sigma} dx, \quad (6)$$

which is the motion equation of the unyielded phase.

To conclude the formulation of the model, we specify the boundary conditions, which are

$$u(x, H, t) = v(x, H, t) = 0, \quad (\text{No-slip and impermeability}), \quad (7)$$

$$u(x, \sigma^+, t) = u_c(t), \quad v(x, \sigma^+, t) = 0, \quad (\text{continuity of the velocity}). \quad (8)$$

The problem given by (3), (6)–(8) is very complex, even in a simple 2D setting. Analytical solutions can be found when the aspect ratio of the channel is sufficiently small. Indeed, in that case, we can make use of the lubrication approximation to determine a semianalytical solution.

Actually, the main difficulty in solving system (3), (6)–(8) lies in tracking the evolution of $\sigma$, where the constitutive equation exhibits a “singularity”. Indeed, from (1), we extrapolate the following relation between the norms $|\tau|$ and $|\dot{\gamma}|$, namely

$$\begin{cases} |\tau| = (\sqrt{\mu}|\dot{\gamma}| + \sqrt{\tau_0})^2 & |\dot{\gamma}| > 0, \\ |\tau| \in [0, \tau_0] & |\dot{\gamma}| = 0, \end{cases} \quad (9)$$

which is a graph exhibiting a singularity in $|\dot{\gamma}| = 0$. A typical way to bypass this difficulty is generally to smooth the constitutive relationship (see [6,29,30] for more details). The positive point of the smoothing methods is their simple implementation in numerical solvers. For example, many commercially CFD codes include inner regularization routines, which smooth the viscosity. However, the issue of the convergence of the yield surface determined using a smoothed constitutive model toward the one of the exact model is far to be settled down. Indeed, to prove that the regularized yield surface converges to the one obtained, solving the exact model requires a uniform convergence results which, at least in the authors' knowledge, has not yet been obtained.

### 3. Nondimensional Formulation

We rescale the system as

$$x = L\tilde{x}, \quad y = \epsilon L\tilde{y}, \quad t = \left(\frac{L}{U}\right)\tilde{t}, \quad u = U\tilde{u}, \quad v = \epsilon U\tilde{v}, \quad u_c = U\tilde{u}_c,$$

$$h = H\tilde{h}, \quad \sigma = H\tilde{\sigma}, \quad \tau = \left(\frac{\mu U}{H}\right)\tilde{\tau}, \quad \dot{\gamma} = \left(\frac{U}{H}\right)\dot{\tilde{\gamma}}, \quad p = \left(\frac{\mu U L}{H^2}\right)\tilde{p},$$
where L is the length of the channel, $H = \max_{[0,L]} h(x)$, and

$$ \epsilon = \frac{H}{L} \ll 1, \qquad (10) $$

is the aspect ratio. We also introduce the dimensionless half-discharge as

$$ Q = \frac{1}{UH} \int_{0}^{H} u \, dy = \int_{0}^{1} \bar{u} \, d\bar{y}. \qquad (11) $$

With this scaling, we get

$$ \dot{\gamma} = \begin{bmatrix} 2\epsilon \frac{\partial \bar{u}}{\partial \bar{x}} & \left( \frac{\partial \bar{u}}{\partial \bar{y}} + \epsilon^2 \frac{\partial \bar{v}}{\partial \bar{x}} \right) \\ \left( \frac{\partial \bar{u}}{\partial \bar{y}} + \epsilon^2 \frac{\partial \bar{v}}{\partial \bar{x}} \right) & 2\epsilon \frac{\partial \bar{v}}{\partial \bar{y}} \end{bmatrix}, $$

$$ \bar{\tau} = \left[ 1 + \frac{\sqrt{B}}{\sqrt{|\dot{\gamma}|}} \right]^2 \begin{bmatrix} 2\epsilon \frac{\partial \bar{u}}{\partial \bar{x}} & \left( \frac{\partial \bar{u}}{\partial \bar{y}} + \epsilon^2 \frac{\partial \bar{v}}{\partial \bar{x}} \right) \\ \left( \frac{\partial \bar{u}}{\partial \bar{y}} + \epsilon^2 \frac{\partial \bar{v}}{\partial \bar{x}} \right) & 2\epsilon \frac{\partial \bar{v}}{\partial \bar{y}} \end{bmatrix}, $$

where

$$ B = \left( \frac{\tau_0 H}{\mu U} \right), $$

is the Bingham number and where

$$ |\dot{\gamma}| = \sqrt{2\epsilon^2 \left[ \left( \frac{\partial \bar{u}}{\partial \bar{x}} \right)^2 + \left( \frac{\partial \bar{v}}{\partial \bar{y}} \right)^2 \right] + \left( \frac{\partial \bar{u}}{\partial \bar{y}} + \epsilon^2 \frac{\partial \bar{v}}{\partial \bar{x}} \right)^2}. \qquad (12) $$

Equation (3)1 becomes

$$ \frac{\partial \bar{u}}{\partial \bar{x}} + \frac{\partial \bar{v}}{\partial \bar{y}} = 0, $$

while Equations (3)2,3 become

$$ \left\{ \begin{array}{l} \epsilon \mathrm{Re} \left( \frac{\partial \bar{u}}{\partial \bar{t}} + \frac{\partial \bar{u}}{\partial \bar{x}} \bar{u} + \frac{\partial \bar{u}}{\partial \bar{y}} \bar{v} \right) = - \frac{\partial \bar{p}}{\partial \bar{x}} + \epsilon \frac{\partial \bar{\tau}_{11}}{\partial \bar{x}} + \frac{\partial \bar{\tau}_{12}}{\partial \bar{y}}, \\ \\ \epsilon^3 \mathrm{Re} \left( \frac{\partial \bar{v}}{\partial \bar{t}} + \frac{\partial \bar{v}}{\partial \bar{x}} \bar{u} + \frac{\partial \bar{v}}{\partial \bar{y}} \bar{v} \right) = - \frac{\partial \bar{p}}{\partial \bar{y}} + \epsilon^2 \frac{\partial \bar{\tau}_{12}}{\partial \bar{x}} + \epsilon \frac{\partial \bar{\tau}_{22}}{\partial \bar{y}}, \end{array} \right. \qquad (13) $$

where

$$ \mathrm{Re} = \frac{\rho U H}{\mu}, $$

is the Reynolds number. Finally, Equation (6) becomes

$$ \epsilon \mathrm{Re} \left[ \dot{v}_c \int_0^1 \bar{\sigma}(\bar{x}, \bar{t}) d\bar{x} \right] = \int_0^1 \left[ -\frac{\partial \bar{p}}{\partial \bar{x}} \bar{\sigma} - \epsilon \frac{\partial \bar{\sigma}}{\partial \bar{x}} \bar{\tau}_{11} + \bar{\tau}_{12} \right] \Big|_{\bar{\sigma}} d\bar{x}. $$

## 4. Leading Order Approximation

Let us now focus on the leading order approximation. In practice, we simplify the model retaining only those terms that do not contain $\epsilon$. This is the so-called lubrication approximation. Dropping the tildes to keep the notation as light as possible, we have
$$
\left\{
\begin{array}{l}
\displaystyle \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0, \\[1.5em]
\displaystyle \frac{\partial \tau_{12}}{\partial y} = \frac{\partial p}{\partial x'}, \qquad \text{(Yielded phase)} \quad y \in [\sigma, h], \\[1.5em]
\displaystyle \frac{\partial p}{\partial y} = 0, \\[1.5em]
\left\{
\begin{array}{l}
\displaystyle \int_0^1 \left[-\frac{\partial p}{\partial x}\sigma + \tau_{12}\right] \Big|_{\sigma} dx = 0, \\[1.5em]
\qquad u = u_c, \quad v = 0,
\end{array}
\right.
\qquad \text{(Unyielded phase)} \quad y \in [0, \sigma],
\end{array}
\right.
\qquad (14)
$$

where

$$
\tau_{12} = \frac{\partial u}{\partial y} \left[ 1 + \frac{\sqrt{B}}{\sqrt{\left| \frac{\partial u}{\partial y} \right|}} \right]^2 . \qquad (16)
$$

The yield condition is given by

$$
\left| \dot{\gamma} \right|_{\sigma} = \left| \frac{\partial u}{\partial y} \right|_{\sigma} = 0. \qquad (17)
$$

Recalling that in the upper yielded part $\partial u / \partial y < 0$, we find

$$
\frac{\partial u}{\partial y} = - \left| \frac{\partial u}{\partial y} \right|,
$$

so that (16) rewrites as

$$
\tau_{12} = - \left[ \sqrt{\left| \frac{\partial u}{\partial y} \right|} + \sqrt{B} \right]^2, \qquad (18)
$$

which, by virtue of (17), gives

$$
\tau_{12}|_{\sigma} = -B.
$$

From (14)₃, we see that $p = p(x, t)$ so that, integrating (14)₂ between $\sigma$ and $y$, we find

$$
\tau_{12} = -B + \frac{\partial p}{\partial x} (y - \sigma). \qquad (19)
$$

Inserting (18) into (19), we get

$$
\left[ \sqrt{\left| \frac{\partial u}{\partial y} \right|} + \sqrt{B} \right]^2 = B + \frac{\partial p}{\partial x} (\sigma - y),
$$
$$
\frac{\partial u}{\partial y} = - \left( \sqrt{B + \frac{\partial p}{\partial x} (\sigma - y)} - \sqrt{B} \right)^2 . \qquad (20)
$$

and hence
Assuming that the pressure drop $\Delta p > 0$, we expect $p_x < 0$ so that $p_x(\sigma - y) > 0$ in the yielded phase. Integrating (20) between $y$ and $h$ and exploiting (7), we find

$$u = \int_{y}^{h} \left( \sqrt{B + \frac{\partial p}{\partial x}(\sigma - \xi)} - \sqrt{B} \right)^{2} d\xi.$$

Setting

$$\mathcal{P} = \frac{p}{B}, \qquad (21)$$

we may rewrite $u$ as

$$\frac{u}{B} = \int_{y}^{h} \left( \sqrt{1 + \mathcal{P}_{x}(\sigma - \xi)} - 1 \right)^{2} d\xi,$$

which, after some algebra, gives

$$\frac{u}{B} = \frac{\mathcal{P}_x}{2} \left[ (\sigma - y)^2 - (\sigma - h)^2 \right] + 2(h - y) + \frac{4}{3\mathcal{P}_x} \left[ (\mathcal{P}_x(\sigma - h) + 1)^{3/2} - (\mathcal{P}_x(\sigma - y) + 1)^{3/2} \right]. \qquad (22)$$

Evaluating (22) on $y = \sigma$, we find

$$-\frac{u_c}{B} = \frac{\mathcal{P}_x}{2}(\sigma - h)^2 + 2(\sigma - h) - \frac{4}{3\mathcal{P}_x} \left[ (\mathcal{P}_x(\sigma - h) + 1)^{3/2} - 1 \right], \qquad (23)$$

where we recall that $u_c$ does not depend on $x$. Rearranging (22) and (23), we find

$$\frac{u}{u_c} = \left[ 1 - \frac{\frac{\mathcal{P}_x}{2}(\sigma - y)^2 + 2(\sigma - y) - \frac{4}{3\mathcal{P}_x} \left[ (\mathcal{P}_x(\sigma - y) + 1)^{3/2} - 1 \right]}{\frac{\mathcal{P}_x}{2}(\sigma - h)^2 + 2(\sigma - h) - \frac{4}{3\mathcal{P}_x} \left[ (\mathcal{P}_x(\sigma - h) + 1)^{3/2} - 1 \right]} \right], \qquad (24)$$

so that we can easily check that $u|_h = 0$ and $u|_\sigma = u_c$, i.e., (7) and (8) are fulfilled.

We now exploit the constraint of incompressibility $u_x + v_y = 0$. Recalling the conditions $v|_h = v|_\sigma = 0$, we find

$$0 = \int_{\sigma}^{h} \frac{\partial v}{\partial y} dy = - \int_{\sigma}^{h} \frac{\partial u}{\partial x} dy.$$

As a consequence,

$$\frac{\partial}{\partial x} \left( \int_{\sigma}^{h} u dy \right) = -u_c \frac{\partial \sigma}{\partial x} + \underbrace{\left( \int_{\sigma}^{h} \frac{\partial u}{\partial x} dy \right)}_{=0},$$

i.e.,

$$\frac{\partial}{\partial x} \left( \int_{\sigma}^{h} \frac{u}{u_c} dy + \sigma \right) = 0, \qquad (25)$$
which we may write also as

$$ \frac{\partial}{\partial x} \left( \frac{\int_{\sigma}^{h} u \, dy + \sigma u_c}{u_c} \right) = 0. \qquad (26) $$

Recalling the definition (11), we easily realize that

$$ Q = \int_{0}^{h} u \, dy = \int_{\sigma}^{h} u \, dy + \sigma u_c, \qquad (27) $$

and since $u_c$ does not depend on $x$, (26) is equivalent to $\frac{\partial Q}{\partial x} = 0$, i.e., the nondimensional discharge is the same at every cross-section $x$ of the channel. Inserting (23) into (25), we find

$$ \frac{\partial}{\partial x} \left( h - \int_{\sigma}^{h} \frac{\frac{P_x}{2}(\sigma - y)^2 + 2(\sigma - y) - \frac{4}{3P_x}[(P_x(\sigma - y) + 1)^{3/2} - 1]}{\frac{P_x}{2}(\sigma - h)^2 + 2(\sigma - h) - \frac{4}{3P_x}[(P_x(\sigma - h) + 1)^{3/2} - 1]} dy \right) = 0. \qquad (28) $$

Integrating the above in $y$, after some algebra, we find

$$ h + \frac{\frac{P_x}{6}(\sigma - h)^3 + (\sigma - h)^2 + \frac{4}{3P_x}(\sigma - h) - \frac{8}{15P_x^2}[(P_x(\sigma - h) + 1)^{5/2} - 1]}{\frac{P_x}{2}(\sigma - h)^2 + 2(\sigma - h) - \frac{4}{3P_x}[(P_x(\sigma - h) + 1)^{3/2} - 1]} = \frac{Q}{u_c}, \qquad (29) $$

where $u_c$ is the core velocity given by (23) and where we assume that $Q$, given by (27), is positive. Introducing the quantity

$$ \ell(x) = h(x) - \sigma(x) > 0, \qquad (30) $$

we rewrite the (29) as

$$ \frac{-\frac{P_x}{6}\ell^3 + \ell^2 - \frac{4}{3P_x}\ell - \frac{8}{15P_x^2}[(1 - P_x\ell)^{5/2} - 1]}{\frac{P_x}{2}\ell^2 - 2\ell - \frac{4}{3P_x}[(1 - P_x\ell)^{3/2} - 1]} = \frac{Q}{u_c} - h. \qquad (31) $$

Next, the integral equation for the unyielded phase (15)$_1$ that can be rewritten as

$$ - \int_{0}^{1} P_x(h - \ell)dx = 1. \qquad (32) $$

Finally, from (23), we observe that

$$ \frac{d}{dx} \left( \frac{u_c}{B} \right) = \frac{d}{dx} \left[ \frac{P_x}{2} \ell^2 - \ell - \frac{4}{3P_x} \left[ (1 - P_x \ell)^{3/2} - 1 \right] \right] = 0. $$

The above is due to the fact that $u_c$ cannot depend on $x$.
In conclusion, the problem to be solved is the following

$$ 
\left\{
\begin{aligned}
& \frac{1}{\mathcal{P}_x} \frac{-\frac{1}{6}(\mathcal{P}_x \ell)^3 + (\mathcal{P}_x \ell)^2 - \frac{4}{3}(\mathcal{P}_x \ell) - \frac{8}{15}[(1 - \mathcal{P}_x \ell)^{5/2} - 1]}{\frac{1}{2}(\mathcal{P}_x \ell)^2 - 2(\mathcal{P}_x \ell) - \frac{4}{3}[(1 - \mathcal{P}_x \ell)^{3/2} - 1]} = \frac{Q}{u_c} - h, \\
& - \int_0^1 \mathcal{P}_x(h - \ell)dx = 1, \\
& \frac{d}{dx}\left[\frac{\mathcal{P}_x}{2}\ell^2 - \ell - \frac{4}{3\mathcal{P}_x}[(1 - \mathcal{P}_x \ell)^{3/2} - 1]\right] = 0,
\end{aligned}
\right.
\quad (33) $$

to which we must add the boundary conditions $\mathcal{P}|_{x=0}$ and $\mathcal{P}|_{x=1}$ if the pressure difference is prescribed, or, alternatively, $Q$ and $\mathcal{P}|_{x=1}$ (actually, in place of $\mathcal{P}|_{x=1}$, we can prescribe $\mathcal{P}|_{x=0}$). The unknowns are: $(\mathcal{P}, \ell, Q)$ in the first case and $(\mathcal{P}, \ell, \mathcal{P}|_{x=1})$ in the second case. The problem is thus formally closed. However, the solution technique varies according to the mechanisms that drive the flow. The case in which the inlet discharge is prescribed has been analyzed in [23]. Here, we focus on the case in which the boundary conditions are $\mathcal{P}|_{x=0}$ and $\mathcal{P}|_{x=1}$, i.e., prescribed inlet-outlet pressure.

## 5. Solution to System (33) When the Pressure Difference is Prescribed

We assume that the pressure drop $\Delta p = p|_{x=0} - p|_{x=1}$, given, i.e., recalling (21)

$$ - \int_0^1 \mathcal{P}_x dx = \frac{\Delta p}{B}. \quad (34) $$

In particular, we stipulate $\mathcal{P}(0) = \frac{\Delta p}{B}$, and $\mathcal{P}(1) = 0$, and introduce the new function

$$ z(x) = -\mathcal{P}_x(x)\ell(x), \quad (35) $$

with $\ell$ given by (30). We remark that $z > 0$.

System (33) can be rewritten as

$$ 
\left\{
\begin{aligned}
& \frac{1}{\mathcal{P}_x} \frac{N(z)}{D(z)} = \frac{Q}{u_c} - h, \\
& \int_0^1 (-z - \mathcal{P}_x h) dx = 1, \\
& \frac{d}{dx} \left[ \frac{D(z)}{\mathcal{P}_x} \right] = 0
\end{aligned}
\right.
\quad (36) $$

where

$$ 
\begin{aligned}
N(z) &= \frac{z^3}{6} + z^2 + \frac{4}{3}z - \frac{8}{15}[(1 + z)^{5/2} - 1] \\
D(z) &= \frac{z^2}{2} + 2z - \frac{4}{3}[(1 + z)^{3/2} - 1],
\end{aligned}
\quad (37) $$

and, exploiting (23),

$$ \frac{u_c}{B} = - \frac{D(z)}{\mathcal{P}_x}. \quad (38) $$
From (36)$_1$, we get

$$ \frac{1}{P_x} = \frac{\left(\frac{Q}{u_c} - h\right) D(z)}{N(z)}, \qquad (39) $$

which, plugged into (36)$_3$, gives

$$ \frac{d}{dx} \left[ \frac{\left( \frac{Q}{u_c} - h \right) D(z)^2}{N(z)} \right] = 0. \qquad (40) $$

Computing the derivative in (40), we find

$$ \frac{dz}{dx} = \frac{h_x D(z) N(z)}{\left(\frac{Q}{u_c} - h\right) [2D'(z)N(z) - N'(z)D(z)]}, \qquad (41) $$

where $Q$ and $u_c$ are unknown at this stage. So, setting $K = Q/u_c$ (unknown parameter), and recalling (39), (41), we obtain the following Cauchy problem

$$ \left\{ \begin{array}{l} P_x = \frac{N(z)}{(K - h(x))D(z)}, \\ \\ z_x = \frac{h_x D(z) N(z)}{(K - h)[2D'(z)N(z) - N'(z)D(z)]}, \\ \\ \mathcal{P}(0) = \frac{\Delta p}{B}, \quad z(0) = z_0, \end{array} \right. \qquad (42) $$

where $z$ is given by (35) and $z_0$ is some initial guess (which, at this stage, plays the role of an unknown parameter as $K$). Solving (42), we find $\mathcal{P} = \mathcal{P}(x, K, z_0)$, and $z = z(x, K, z_0)$. To determine $K$ and $z_0$, we impose the second boundary $\mathcal{P}$, i.e., $\mathcal{P}(1) = 0$, and (36)$_2$, namely

$$ \left\{ \begin{array}{l} \mathcal{P}(1, K, z_0) = 0, \\ \\ \int_0^1 \left[ z(x) + h(x) \underbrace{\frac{N(z)}{(K - h(x))D(z)}}_{P_x} \right] dx = -1. \end{array} \right. \qquad (43) $$

The yield surface follows from (35), i.e.

$$ \sigma(x) = h(x) + \frac{z(x)}{P_x} = h(x) + z(x) \frac{N(z(x))}{(K - h(x))D(z(x))}. \qquad (44) $$

When $h_x = 0$, i.e., $h \equiv 1$, (42)$_2$ entails $z(x) = \bar{z}_0$, while (42)$_1$, (42)$_3$ implies

$$ P_x = \frac{N(\bar{z}_0)}{(\bar{K} - 1)D(\bar{z}_0)}, \quad \Rightarrow \quad P(x) = \frac{N(\bar{z}_0)}{(\bar{K} - 1)D(\bar{z}_0)} x + \frac{\Delta p}{B}. \qquad (45) $$

Therefore, system (43) rewrites

$$ \left\{ \begin{array}{l} \frac{N(\bar{z}_0)}{(\bar{K} - 1)D(\bar{z}_0)} + \frac{\Delta p}{B} = 0, \\ \\ \bar{z}_0 + \frac{N(\bar{z}_0)}{(\bar{K} - 1)D(\bar{z}_0)} = -1, \end{array} \right. \qquad (46) $$
which in turn entails

$$ \frac{N(\bar{z}_o)}{(\bar{K}-1)D(\bar{z}_o)} = -\frac{\Delta p}{B}, \quad \implies \quad (45) \qquad \mathcal{P}(x) = \frac{\Delta p}{B}(1-x), $$

and

$$ \bar{z}_o = \frac{\Delta p}{B} - 1, \quad \bar{K} = 1 - \frac{BN(\bar{z}_o)}{\Delta p D(\bar{z}_o)}, \qquad (47) $$

with N and D given by (37). Next, exploiting (44), we obtain the classical relation

$$ \sigma_o = 1 - \frac{\bar{z}_o}{\frac{\Delta p}{B}} = \frac{B}{\Delta p}. $$

We remark that our model converges to the classical results when the channel walls are flat.

### 5.1. Simulations

To solve problem (43) coupled with (46) we start considering $\bar{z}_o$ and $\bar{K}$ given by (47), i.e., the parameters corresponding to the flat channel $h(x) \equiv 1$. We then consider a generic wall profile $h(x) \neq 1$, and set up an iterative minimum-finding scheme. We proceed by determining a grid around the values $\bar{z}_o$ and $\bar{K}$ and, after having selected a pair $(z_o, K)$ in the grid, we solve the Cauchy problem (42), obtaining the pair $(\mathcal{P}(x; z_o, K), z(x; z_o, K))$ and compute

$$ d_1 = |\mathcal{P}(1, K, z_o)|, $$

$$ d_2 = \left| \int_0^1 \left( z(x) + h(x) \frac{N(z)}{(K - h(x))D(z)} + 1 \right) dx \right|. $$

We then repeat the procedure for all grid values and for each pair $(z_o, K)$ we compute the norm

$$ \|d\| = \sqrt{d_1^2 + d_2^2}. $$

We stop the procedure when we reach a pair $(z_o^*, K^*)$, whose corresponding $\|d\|$ is smaller than a prescribed tolerance. Once $(z_o, K)$ have been determined, $\sigma$ is given by (44).

Figure 3 displays the results of numerical simulations when $h(x) = 1 \pm \frac{1}{\delta} \sin(\pi x)$, $\delta = 20$, and $\frac{\Delta p}{B} = 1, 5, 2, 2.5$.

<figure><img src="image_6.png" /><figcaption>Figure 3. Simulations for $h(x) = 1 \pm \frac{1}{\delta} \sin(\pi x)$, with $\delta = 20$, for various values of $\frac{\Delta p}{B}$. We note the physical consistence of the results obtained, i.e., the width of the inner core increases as $\Delta p$ decreases. Moreover, the simulation highlights a well-known phenomenon already highlighted [12]: the monotony of $\sigma(x)$ is opposite to the one of $h(x)$.</figcaption></figure>
## 5.2. Approximate Solution

In this section, we illustrate a technique to find an approximate solution to (33) when the pressure difference is prescribed and peculiar conditions on the data are fulfilled. We recall the generalized form of the binomial Newton formula

$$ \sum_{k=0}^{\infty} \binom{\alpha}{k} x^k = (1+x)^k, \quad \text{for } |x| < 1, \qquad (48) $$

where, given $\alpha \in \mathbb{R}$ and $k \in \mathbb{N}$, $\binom{\alpha}{k}$ is the generalized binomial coefficient

$$ \binom{\alpha}{0} = 1, \quad \binom{\alpha}{k} = \frac{\alpha(\alpha-1)\cdots(\alpha-(k-1))}{k!}, \quad \text{for } k \ge 1. $$

So, assuming

$$ |\mathcal{P}_x(\sigma - h)| < 1, \qquad (49) $$

we exploit (48) in (22),

$$ \frac{u}{B} = -\frac{1}{12} \mathcal{P}_x^2 [(\sigma - h)^3 - (\sigma - y)^3] + \mathcal{O}(|\mathcal{P}_x(\sigma - h)|^3) \qquad (50) $$

where the third order terms have been neglected. Proceeding similarly in (23), we obtain

$$ \frac{u_c}{B} = \frac{1}{12} \mathcal{P}_x^2 (h - \sigma)^3 + \mathcal{O}(|\mathcal{P}_x(\sigma - h)|^3). \qquad (51) $$

Then, neglecting the higher order terms, the velocity field $u$ inside the channel can be rewritten as

$$ u = u_c = \frac{1}{12} B \mathcal{P}_x^2 (h - \sigma)^3, \quad 0 < y < \sigma, \qquad (52) $$

$$ u = u_c - \frac{1}{12} B \mathcal{P}_x^2 (y - \sigma)^3, \quad \sigma < y < h, \qquad (53) $$

where we remind that $p_x < 0$. We remark that the approximate expression (53) fulfills the no-slip condition.

Let us now rewrite Equation (25) as follows

$$ \frac{\partial}{\partial x} \left( \int_{\sigma}^{h} \frac{u}{B} dy \right) = - \frac{u_c}{B} \frac{\partial \sigma}{\partial x}. $$

After some algebra and using (50) and (51), we obtain

$$ \sigma(x) = -3h(x) + C, \qquad (54) $$

where $C$ is an unknown at this stage. To determine $C$, we exploit Equations (54) and (32), namely

$$ \left\{ \begin{array}{l} \displaystyle \int_{0}^{1} p_x \sigma \, dx = -B, \\ \sigma = -3h + C. \end{array} \right. $$

Imposing the boundary conditions for the pressure, i.e., $p(x = 0) = \Delta p$, and $p(x = 1) = 0$, we get this explicit expression for $C$

$$ C = \frac{B - 3 \int_{0}^{1} p_x h \, dx}{\Delta p}, $$
which, plugged in (54), gives

$$ \sigma = -3h(x) + \frac{B - 3 \int_{0}^{1} p_x h \, dx}{\Delta p}. \qquad (55) $$

At this point, we are in position to determine $p(x)$. We take (51), which, recalling (30), can be rewritten as

$$ \frac{u_c}{B} = \frac{1}{12} \mathcal{P}_x^2 \ell^3 = 0, $$

and differentiate it with respect to $x$

$$ p_{xx} + \frac{3}{2} p_x \left( \frac{\ell_x}{\ell} \right) = 0. $$

Now, exploiting the (54), we obtain the following boundary value problem based on an integrodifferential equation

$$ \left\{ \begin{array}{l} p_{xx} + \frac{6h_x}{3 \int_0^1 p_x h \, dx - B} p_x = 0, \\ 4h(x) + \frac{0}{\Delta p} = \Delta p, \\ p|_{x=0} = \Delta p, \\ p|_{x=1} = 0. \end{array} \right. \qquad (56) $$

**Remark 1.** As already stated, the approximation is meaningful when (49) is fulfilled, i.e., $|p_x|/|\sigma - h| < 1$.

Considering, for simplicity, $h \equiv 1$, $p_x = -\Delta p$ and, exploiting (55), we have $\sigma = \frac{B}{\Delta p}$. Hence,

$$ \frac{|p_x|}{B} |\sigma - 1| < 1 \implies \sigma > \frac{1}{2}. $$

At the same time, in order to prevent the flow come to a stop, $\sigma < 1$ i.e., $\frac{B}{\Delta p} < 1$. Hence, the approximation above developed is expected to hold true when

$$ \frac{1}{2} < \frac{B}{\Delta p} < 1. \qquad (57) $$

To solve the problem (56), we set up an iterative procedure. As first guess, we consider the linear function $p(x) = \Delta p(1 - x)$, and use (55), to give the first guess of the yield surface. The solution at the $n$th, $n \ge 1$, step is obtained solving this boundary value problem

$$ \left\{ \begin{array}{l} p_{xx}^{(n)} + 6h_x \left[ \displaystyle\frac{3 \int_0^1 p_x^{(n-1)} h \, dx - B}{\Delta p} \right]^{-1} p_x^{(k)} = 0, \\ p_x^{(n)} \Big|_{x=0} = \Delta p, \\ p_x^{(n)} \Big|_{x=1} = 0, \end{array} \right. $$
and the corresponding yield surface $\sigma^{(n)}$ is obtained by (55). We stop the procedure when $\left\|p^{(n)} - p^{(n-1)}\right\|_{L^2([0,1])}$ becomes smaller than a prescribed tolerance.

Figure 4 displays the yield surface when the wall profile is $h(x) = 1 \pm \frac{1}{\delta} \sin(\pi x)$, with $\delta = 20$, for various values of $B/\Delta p$.

<figure><img src="image_7.png" /><figcaption>Figure 4. Plots of the approximate yield surface $\sigma(x)$ when $h(x) = 1 \pm \frac{1}{20} \sin(\pi x)$, for three values of $\frac{\Delta p}{B}$, namely $\frac{5}{4}$, $\frac{3}{2}$, and $\frac{7}{4}$, fulfilling (57).</figcaption></figure>

The approximation here developed holds true for $|P_x||\sigma - h| < 1$. So, to validate the procedure, we have to verify if $|P_x||\sigma - h| < 1$ is fulfilled. Figure 5 shows the plot of $|P_x||\sigma - h|$ for the three cases displayed in Figure 4.

<figure><img src="image_8.png" /><figcaption>Figure 5. Plots of $|P_x||\sigma - h|$ for the cases in Figure 4.</figcaption></figure>

In the right panel of Figure 5, we notice that the condition under which the approximation is valid is not satisfied at the inlet and outlet.

In Figure 6, we report the comparison between $\sigma(x)$ given by (44) and $\sigma(x)$ obtained by the approximate model (55), (56). Figure 7 displays again the comparison between $\sigma(x)$ given by (44) and the approximated one, i.e., $\sigma(x)$ given by (55), (56), but now $\frac{B}{\Delta p} < \frac{1}{2}$, i.e., out of the range of validity of the approximate model. Figure 8 shows the comparison between $p(x)$ obtained solving (42), (43) and solving (56) when $\frac{B}{\Delta p} = 0.75$, and $\frac{B}{\Delta p} = 0.91$. Looking at Figure 7, one realizes that the two interfaces are quite similar, as in Figure 6, where $\frac{1}{2} < \frac{B}{\Delta p} < 1$. On the other hand, the pressure profiles, displayed in Figure 8, are monotonously decreasing but show a different behavior, although we are in the range of the approximate model.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Comparisons between $\sigma(x)$ given by (44), continuous curve, and $\sigma(x)$ given by (56), (55), dotted curve. In the upper panels, $\frac{B}{\Delta p} = \frac{10}{13}$. In the lower panels, $\frac{B}{\Delta p} = \frac{2}{3}$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Comparisons between $\sigma(x)$ given by (44), continuous curve, and $\sigma(x)$ given by (56), (55), dotted curve, when $\frac{B}{\Delta p} < \frac{1}{2}$, i.e., out of the validity range for the approximation model. In the upper panels, $\frac{B}{\Delta p} = \frac{1}{5}$. In the lower panels, $\frac{B}{\Delta p} = \frac{1}{10}$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Comparisons between $p(x)$ given by (44), continuous curve, and $p(x)$ given by the approximate model (56), dotted curve, for $\frac{B}{\Delta p} = 0.75$, $\frac{B}{\Delta p} = 0.91$. On the left side, the wall profile is $h(x) = 1 + \frac{1}{3}\sin(\pi x)$; on the right side, it is $h(x) = 1 - \frac{1}{3}\sin(\pi x)$—in both cases, with $\delta = 20$. We remark that the two profiles are different. However for the dynamics purpose, the important quantity is the pressure gradient and not the pressure.</figcaption></figure>

### 5.3. Comparison with the Pressure Driven Bingham Flow in a Channel

Problem (56) is similar to the one governing the Bingham flow in a channel, namely (see [12], Equation (39))

$$ \left\{ \begin{array}{l} p_{xx} + \left( \frac{6h_x \Delta p}{3h \Delta p + 2 \int_0^1 p_x h dx - B} \right) p_x = 0, \\ p|_{x=0} = \Delta p, \quad p|_{x=1} = 0, \end{array} \right. \qquad (58) $$

with the yield surface given by (see again [12], equation (34))

$$ \sigma(x) = -2h(x)) + \frac{B}{\Delta p} - \frac{2}{\Delta p} \int_{0}^{1} p_x h \, dx. \qquad (59) $$

In Figure 9, we report the yield surface given by (59), i.e., the Bingham yield surface, and $\sigma(x)$ obtained solving (42) and (43), i.e., the Casson yield surface. We have considered three cases, $\frac{B}{\Delta p} = 0.1, 0.5$ and $0.86$.

The plots show that the two models give rise to similar curves when $\frac{B}{\Delta p} = 0.1$ and $0.5$. If $\frac{B}{\Delta p} = 0.86$, the difference between the two curves is much more evident. We remark that, when $h(x) = 1 + \delta \sin(\pi x)$, left panels, the largest difference occurs in the central region. On the contrary, when $h(x) = 1 - \delta \sin(\pi x)$, right panels, the largest difference occurs in the external regions.
<figure><img src="image_12.png" /><figcaption>Figure 9. Comparisons between $\sigma(x)$ given by (43) and (44), i.e., Casson model, continuous curve, and $\sigma(x)$ given by model (58) and (59), Bingham model, dotted curve. We observe that the Bingham yield surface is close to the Casson yield surface when $\frac{B}{\Delta p}$ is “small”. The two surfaces highlight once more the fundamental feature of viscoplastic flow, i.e., the anti-monotony with respect to the vessel wall.</figcaption></figure>

## 6. Peristaltic Flow

We now consider the case where the duct walls move as traveling waves, i.e., the peristaltic flow. We indeed assume that

$$h(x,t) = 1 + \frac{1}{\delta} \sin(2\pi(x-t)), \qquad (60)$$

that is a sinusoidal profile whose wave length and speed are equal to 1. Next, the boundary conditions of adhesion on $y = h$, i.e., (7), change and become

$$u(x,h,t)=0, \quad v(x,h,t) = \frac{\partial h}{\partial t}.$$
Recalling (27) and following the same approach of Section 4, we have that

$$0 = \frac{\partial Q(x,t)}{\partial x} = \frac{\partial}{\partial x} \int_{0}^{h(x,t)} u dy = \frac{\partial}{\partial x} \int_{0}^{h(x,t)} \frac{\partial v}{\partial y} dy,$$

which yields

$$\frac{\partial Q}{\partial x} = -\frac{\partial h}{\partial t}.$$

or, because of (60),

$$\frac{\partial}{\partial x}(Q - h) = 0,$$

from which we conclude that $Q(x,t) - h(x,t)$ does not depend on $x$. Now, assuming that the inlet discharge $Q_{in}(t) = Q(0,t)$, we have

$$Q(x,t) - h(x,t) = Q_{in}(t) - h_{in}(t),$$

where $h_{in}(t) = h(0,t)$. Proceeding as in Section 5.2, we insert (53) in (27) and obtain

$$Q_{in} + h - h_{in}(t) = h u_{c} - \frac{1}{12} B \mathcal{P}_{x}^{2} \frac{(h - \sigma)^{4}}{4} = h u_{c} - u_{c} \frac{(h - \sigma)}{4} = \frac{u_{c}}{4} (3h - \sigma). \quad (61)$$

Now, setting

$$\begin{aligned} A(x,t) &= Q_{in}(t) + h(x,t) - h_{in}(t) \\ &= Q_{in}(t) + \frac{1}{\delta} [\sin(2\pi(x-t)) + \sin(2\pi t)], \end{aligned}$$

rearranging (61) yields

$$\sigma(x,t) = \frac{3h(x,t)u_c(t) - 4A}{u_c(t)}, \quad (62)$$

with $u_c(t)$ still unknown at this stage. To determine the unknowns $(u_c, \sigma)$, we exploit (32)$_2$ and (61), obtaining this system

$$\left\{ \begin{array}{l} \displaystyle \int_{0}^{1} \mathcal{P}_{x} \sigma \, dx = -1, \\ \\ \displaystyle \sigma = \frac{3hu_{c} - 4A}{u_{c}}. \end{array} \right. \quad (63)$$

To solve this problem, we exploit (53)$_1$ and (62) to determine $\mathcal{P}_x$ as a function of $u_c$. We thus have

$$\mathcal{P}_x = -\sqrt{\frac{12}{B} \frac{u_c^2}{(4A - 2hu_c)^{\frac{3}{2}}}}, \quad (64)$$

where we choose the negative solution according to our assumption about $p_x$. Now, we put (64) in (63)$_1$, getting

$$\sqrt{\frac{12}{B}} \int_{0}^{1} u_{c} \frac{3hu_{c} - 4A}{(4A - 2hu_{c})^{\frac{3}{2}}} dx = 1, \quad (65)$$

with $u_c = u_c(t)$.
## Simulations

To solve Equation (65), we have to impose some conditions on the integrand function. We require

$$4A - 2hu_c > 0 \Rightarrow \frac{u_c}{A} < \frac{2}{h'}$$

$$3hu_c - 4A > 0 \Rightarrow \frac{u_c}{A} > \frac{4}{3h'}$$

which, setting $\psi(x, t) = \frac{u_c(t)}{A(x, t)}$, we rewrite as

$$\frac{4}{3h} < \psi < \frac{2}{h'}$$

So, rearranging (65), we get

$$\sqrt{\frac{12}{B}} \int_{0}^{1} \sqrt{A\psi} \frac{3h\psi - 4}{(4 - 2h\psi)^{\frac{3}{2}}} dx = 1. \qquad (66)$$

So, dividing the time interval $[0, 1]$ in $n$ steps, we solve numerically (66) at any time step. The results of simulations are displayed in Figures 10 and 11.

<figure><img src="image_13.png" /><figcaption>Figure 10. Yield surface against the wall profile at five different times. $\delta = 20$, $B = 5$ and $Q(t) = 1$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 11. Yield surface $\sigma$ varying the parameters $Q$ and $B$. On the left panel, $Q(t) = 1$ while $B$ varies. On the right panel, $B$ is kept constant while $Q$ varies. In both cases, we set $\delta = 20$.</figcaption></figure>
## 7. Conclusions

In this paper, we have presented a mathematical model for a Casson flow in a symmetrical channel of varying amplitude whose walls can move over time as a traveling wave. The formulation of the fluid dynamics problem is obtained by imposing the mass and momentum balance. The latter, written for the central rigid core, results in an integral equation. We have thus determined an explicit expression for the velocity field and for the yield surface (which, being unknown, is a free boundary). The problem has been solved in two cases: (i) the driving force of the flow is the pressure difference applied between inlet and outlet; (ii) the inlet flow rate is imposed and the walls of the channel are animated by a traveling wave (peristaltic flow). Numerical simulations of the peristaltic flow have shown that as the Bingham number increases, and as the flow rate decreases, the yield surface tends to occupy the entire channel. Although an exhaustive parametric analysis has not been performed, the simulations carried out seem to show how the yield surface and channel walls oscillate in phase. This result (although not definitive) may have an important practical repercussion. Indeed, since $\sigma$ oscillates in phase with $h$, it is unlikely that the two surfaces touch and the flow stops. This fact, absent in the Bingham flows, could justify, if confirmed by further studies, the vasomotor action at the microcirculation level.

Regarding the analysis of the dynamics driven by the pressure gradient $\Delta p$, a comparison was made between the Bingham and the Casson flow. The results obtained seem to show a certain sensitivity to the $B/\Delta p$ parameter. Specifically, the two yield surfaces are similar when $B/\Delta p \sim 0.5$. As soon as $B/\Delta p$ tends to 0.9, the two surfaces detach significantly. In any case, we have found a well-known characteristic of viscoplastic flows in channels of variable amplitude: the yield surface and the channel wall have opposite monotonicity, that is, the plug shrinks (or widens) as the width of the channel increases (or decreases).

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

**Funding:** This work was funded in part by the Fondazione Cassa di Risparmio Firenze, "Bando Ricerca Scientifica e Tecnologica 2020", progetto "Uso di Equazioni Implicite per Descrivere il Comportamento di Materiali Soggetti ad Elevate Pressioni". L. I. Palade cordially acknowledges his thankfulness to DR CNRS Marie-Christine Baietto, Dean for Scientific Affairs at INSA-Lyon, for help with securing financial support for this work. Funding was provided via a research grant awarded to him. S. Guadagni received additional support from INdAM (Italy).

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** The data presented in this study.

**Conflicts of Interest:** The authors declare no conflict 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-Zeitschrift* **1926**, *39*, 291–300. [CrossRef]

3. Casson, N. *Rheology of Disperse Systems*; Mill, C.C., Ed.; Pergamon Press: Oxford, UK, 1959.

4. Fincke, A. Beitrage zur Losung Rheologischer Probleme in Der Schokoladentechnologie. Ph.D. Dissertation, TH Karlsruhe, Karlsruhe, Germany, 1961.

5. Bird, R.B.; Dai, G.C.; Yarusso, B.J. The rheology and flow of viscoplastic materials. *Rev. Chem. Eng.* **1983**, *1*, 1–70. [CrossRef]

6. Mitsoulis, E. Flows of viscoplastic materials: Models and computations. *Rheol. Rev.* **2007**, *13*, 135–178.

7. Huilgol, R.R. *Fluid Mechanics of Viscoplasticity*; Springer: Berlin/Heidelberg, Germany, 2015.

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

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

10. Fasano, A.; Sequeira, A. *Hemomath: The Mathematics of Blood*; Springer: Berlin/Heidelberg, Germany, 2017.
11. 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.

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

13. Fusi, L.; Farina, A.; Rosso, F. Planar squeeze flow of a Bingham fluid. *J. Non-Newton. Fluid Mech.* **2015**, *225*, 1–9. [CrossRef]

14. Safronchik, A.I. Nonstationary flow of a visco-plastic material between parallel walls. *J. Appl. Math. Mech.* **1959**, *23*, 1314–1327. [CrossRef]

15. Muravleva, L. Squeeze plane flow of viscoplastic bingham material. *J. Non-Newton. Fluid Mech.* **2015**, *220*, 148–161. [CrossRef]

16. Muravleva, L. Axisymmetric squeeze flow of a Casson medium. *J. Non-Newton. Fluid Mech.* **2019**, *267*, 35–50. [CrossRef]

17. Putz, A.; Frigaard, I.A.; Martinez, D.M. On the lubrication paradox and the use of regularization methods for lubrication flows. *J. Non-Newton. Fluid Mech.* **2009**, *163*, 62–77. [CrossRef]

18. Argyriadi, K.; Vlachogiannis, M.; Bontozoglou, V. Experimental study of inclined film flow along periodic corrugations: The effect of wall steepness. *Phys. Fluids* **2006**, *18*, 012102. [CrossRef]

19. Blyth, M.G.; Pozrikidis, C. Film flow down an inclined plane over a three-dimensional obstacle. *Phys. Fluids* **2006**, *18*, 052104. [CrossRef]

20. Hayes, M.; O'Brien, S.B.G.; Lammers, J.H. Green's function for steady flow over a small two-dimensional topography. *Phys. Fluids* **2000**, *12*, 2845–2858. [CrossRef]

21. Marner, F.; Gaskell, P.H.; Scholle, M. A complex-valued first integral of Navier-Stokes equations: Unsteady Couette flow in a corrugated channel system. *J. Math. Phys.* **2017**, *58*, 043102. [CrossRef]

22. Scholle, M.; Gaskell, P.H.; Marner, F. A Potential Field Description for Gravity-Driven Film Flow over Piece-Wise Planar Topography. *Fluids* **2019**, *4*, 82. [CrossRef]

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

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

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

26. 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]

27. 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]

28. Farina, A.; Fusi, L. Viscoplastic fluids: Mathematical modeling and applications. In *Non-Newtonian Fluid Mechanics and Complex Flows*; C.I.M.E. Summer School, Levico Terme, 2016, Lecture Notes in Mathematics 2212; Farina, A., Mikelic, A., Rosso, F., Eds.; Springer: Berlin/Heidelberg, Germany, 2018; pp. 229–298.

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

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