Viscoplastic flows in channels with small aspect ratio: Bingham versus regularised models

We investigate the two-dimensional flows of a viscoplastic fluid in symmetric channels with impermeable walls under no-slip boundary conditions. To characterise the mechanical response of the viscoplastic fluid we consider both the celebrated Bingham model and a very general class of its regularisations. In order to make the problem amenable to analysis, we assume that the aspect ratio of the channel is small so that the lubrication approximation can be used. This allows us to obtain analytical solutions, perform an asymptotic analysis of the regularised solutions and compare the results predicted by the Bingham model and its regularisations. We find that in the limit as the regularisation parameter tends to zero, the regularised flow tends to those predicted by the Bingham model only in plane channels. In channels with curved walls, the results are instead markedly different.


Introduction
Several materials like mineral slurries, paints, suspensions, crude oils, and extrusions in 3D printing are modelled as viscoplastic fluids.The main constitutive feature of these models is represented by the so-called yield stress, namely a threshold to be exceeded for the fluid to start flowing.Based on the pioneering studies by Bingham (1922), Casson (1959) and Herschel and Bulkley (1926), several constitutive models for the mechanical response of viscoplastic fluids have been proposed as attested in the detailed reviews by Farina and Fusi (2018) and Huilgol and Georgiou (2022).
However, the existence of a real yield stress has been widely discussed in the literature (Astarita, 1990;Barnes, 1999;Barnes & Walters, 1985;Frigaard & Nouar, 2005;Frigaard, Paso, & de Souza Mendes, 2017) and, still nowadays, is the subject of an open debate.Indeed, the existence of a threshold for the stress above which the fluid starts flowing might be just a consequence of the experimental techniques and apparatuses used to investigate the rheology of the materials.Whether or not a fluid exhibits a yield stress might be due to the time and length scales of the observations and of the applied forces.For instance, in the famous pith drop experiment performed by Edgeworth, Dalton, and Parnell (1984), it took almost nine years for a drop of pitch to fall.This experiment clearly shows that it is inappropriate to introduce a yield stress for pitch as claimed by Bhatia, Aggarwal, Chari, and Jain (1977), and suggests that pitch can be regarded as a fluid whose effective viscosity at small rates of shear is so large that the shear rate is too small to be measured by any instrument.This might be the case of many other fluids that seem to exhibit a yield stress.
On the other hand, the numerical integration of the equations governing the flows of fluids that seem to exhibit a yield stress is very challenging.The constitutive equation of a yield stress fluid relates the rate of strain tensor to the stress tensor (Rajagopal, 2011).Since such a relation is, in general, not invertible, one cannot obtain a partial differential equation for the velocity field as in the case of Newtonian fluids.Instead, one has to solve simultaneously the equation of balance of linear momentum and the constitutive relation: a tough task.In addition, when the shear stress attains the critical limit, the effective viscosity diverges and no straightforward computation is possible.To overcome these difficulties, regularisation methods have been proposed.They consist in replacing the singular constitutive model with a smooth approximation that depends on a suitable positive regularisation parameter (say ).In the limit as  → 0 + , the regularised model tends pointwise to the exact viscoplastic model.For an exhaustive discussion on regularisation methods, we refer the readers to Frigaard and Nouar (2005) and references therein.
Introducing a regularised model is natural from a mathematical point of view as smoothing a non-differentiable minimisation problem allows the use of advantageous variational techniques that can be implemented in highly performing numerical schemes.For instance, many commercial computational fluid dynamics (CFD) codes do not solve the equations governing the flows of a Bingham fluid, but proceed with internal regularisation routines.Moreover, for simple test problems (e.g. the flow in a flat channel or in a cylindrical tube) the velocity obtained by using a regularised constitutive model converges uniformly to the one obtained by the exact (non-smooth) model as  → 0 + .Finally, the introduction of a regularised model is justified also from a physical point of view as the regularised effective viscosity approaches a large but finite value as the shear rate tends to zero.Since no one is able to establish with absolute certainty whether or not a fluid actually exhibits a yield stress, a regularised model may be more physically plausible than an exact viscoplastic one.
Particular attention has to be paid to the convergence of the flow predicted by the regularised model to the one predicted by the exact model as  → 0 + .Limiting to Bingham fluids, the convergence in mean of the regularised flows to the exact flows has been addressed by Duvaut and Lions (1976), Lions, Glowinski, and Tremolieres (1976) and Bercovier and Engelman (1980).In particular, for steady Stokes flows in a bounded two-dimensional domain with Lipschitz boundary subjected to Dirichlet boundary conditions, the difference between the exact Bingham and regularised flows is estimated in the  1 -norm as   , with  ∈ (0, 1) depending essentially on the particular regularised model used.However, the location of the yield surface (which is undoubtedly the main problem in viscoplastic flows) lies outside the convergence in mean of the regularised velocity field, since it is linked to the convergence of the extra stress.To show that the pseudo-yield surface converges to the yield surface predicted by the exact non-smooth model one needs results on the uniform convergence of the regularised stress to the exact stress which, except for very few simple cases (e.g.flow in a flat channel), has not yet been proved.
The fact that there may be significant discrepancies in the predictions obtained by using the Bingham model and its regularisations has been recently highlighted in studies concerning the stability of simple steady flows.For instance, Falsaperla, Giacobbe, and Mulone (2020) have shown that the Bingham flow down an inclined plane is unconditionally stable, whereas, in perfect agreement with the experimental observations of Noma (2021), Calusi, Farina, Fusi, andRosso (2022) and Fusi, Calusi, Farina, and Rosso (2022) have proven that instabilities may occur above a critical tilt angle when using a regularisation of the Bingham model which (because of its simple functional form) is known as the simple model (Allouche, Frigaard, & Sona, 2000), or the Papanastasiou (1987) model. Furthermore, Fusi, Farina, Rajagopal, andVergori (2022) have shown that the critical linear instability thresholds for Poiseuille flows of a viscoplastic fluid obtained by using the simple and Papanastasiou models are different from the ones predicted by the exact Bingham model.
The main purpose of this paper is to investigate to what extent the flows predicted by the exact Bingham model and its regularisations differ.In particular, we focus on the possible differences between the pseudo-yield surfaces and the 'exact' yield surface.To this aim, moved by the intention to look for analytical solutions, we consider the flow in symmetric channels with small aspect ratio so that we can introduce the lubrication approximation and make the governing equations simpler to solve analytically.We solve the approximated governing equations by considering the exact Bingham model and a very general class of regularised models.Such a class contains the most widely used regularisations of the Bingham model available in the literature such as the simple, Bercovier and Engelman (1980) and Papanastasiou models.We then study the asymptotic behaviour of the regularised solutions and compare the results with those predicted by the exact Bingham model.We see that, irrespective of the regularisation considered, in the limit as  → 0 + , the regularised flow and pseudo-yield surface tend to expressions that depend exclusively on the profile of the channel walls and the value of the Bingham number.The resulting asymptotic pseudo-yield surface coincides with the yield surface predicted by the exact Bingham model when the fluid flows in a plane channel.If the flow takes place in a non-plane channel, the two surfaces are markedly different.This is due to the fact that only in plane channels regularised models envisage, in the limit as  → 0 + , a rigid core like when considering the exact Bingham model.
To realise whether or not the lubrication approximation may cause the differences in the predictions, we solve numerically the set of governing equations without neglecting the small terms in the aspect ratio, by considering the simple, Bercovier & Engelman and Papanastasiou models, and for plane, sinusoidal, divergent and convergent straight-walled channels.To this end, we set up a numerical model with Comsol Multiphysics ® and find that, irrespective of the three regularised models used, at small enough values of , the numerical pseudo-yield surface practically coincides with the asymptotic one.We then conclude that also the numerical solution of the full set of governing equations confirms the significant differences in the yield and pseudo-yield surfaces predicted, respectively, by the Bingham and regularised models when the channel is not plane, whereas the numerical and asymptotic pseudo-yield surfaces coincide with the exact yield surface for plane channels.
The structure of the paper is as follows.In Section 2, we recall the basic features of the Bingham model and introduce a very general class of regularised models.In Section 3, we introduce the equations governing the two-dimensional flows of a viscoplastic fluid in a symmetric channel and the appropriate set of boundary conditions.Next, assuming that the channel has small aspect ratio we introduce the lubrication approximation (Section 4) and solve the approximated boundary value problem by considering the Bingham model (Section 5) and its regularisations (Section 6).We then study the behaviour of the regularised flow in the limit as the regularisation parameter tends to zero (Section 7) and compare the results of the asymptotic analysis with those predicted by the Bingham model (Section 8).The paper terminates with some remarks and discussion on the results obtained (Section 9).

Models for viscoplastic fluids
The Cauchy stress tensor  of an incompressible viscous fluid is usually split into the sum of an arbitrary spherical part −, where  is the Lagrange multiplier associated with the constraint of incompressibility, and a traceless tensor field  usually called the deviatoric or extra stress.The indeterminate scalar field  equals then the opposite of the mean normal stress and for this reason it is often called the ''mechanical pressure''.This terminology is however inappropriate as  is not given by a constitutive relation.For a detailed discussion on the use and misuse of the term 'pressure' we refer the reader to Rajagopal (2015).As the extra stress  is concerned, it is usually modelled as a function of the rate of strain tensor  = [∇ + (∇)  ]∕2, where  denotes the velocity field.In the literature, there are available a plethora of mathematical models for the extra stress tensors of various classes of viscous fluids.Based on the pioneer works of Bingham (1917Bingham ( , 1922)), one of most used constitutive relation for the mechanical response of an incompressible viscoplastic fluid reads where the positive constants  and  0 represent, respectively, the plastic viscosity and the so-called yield stress, and, given a symmetric second-order tensor  with components   (, Observe that if  is also traceless (like  and ), then || equals the square root of the opposite of the second principal scalar invariant of .
The constitutive Eqs.(1), known as the Bingham model, tell us that the deviatoric stress tensor can be related to the rate of strain tensor providing that || exceeds the yield stress.Otherwise,  is indeterminate as the motion of the viscoplastic fluid is rigid.Consequently, provided || >  0 , the effective viscosity   is well defined and (1) can be implemented in numerical codes.However, as || →  0 , the effective viscosity blows up.The non-smoothness of the effective viscosity makes impossible any straightforward numerical computation (Frigaard & Nouar, 2005) and may lead to mathematical paradoxes (Fusi, Farina, & Rosso, 2014).This explains why several regularised viscosity models (called also regularised models or regularisations) have been proposed to replace the Bingham model (1).
By a regularised model it is usually meant a constitutive relations for the deviatoric stress in terms of the rate of strain tensor in the form where the regularisation parameter  is positive and has the same physical dimensions as the rate of strain, the effective viscosity for all  ∈ R, and at small enough values of  it tends to the effective viscosity in the Bingham model in the sense that Unlike the Bingham model ( 1), constitutive relations in the form (2) are not based on the notion of yield stress.In any regularised model, the regularisation parameter  is taken so that (0 <)∕ 0 ≪ 1, whence, in view of (3), the effective viscosity is very large but finite also at infinitesimally small strain rates.Only in the limit as ∕ 0 → 0 + the effective viscosity blows up as the strain rate tends to zero.In fact, in this limit the constitutive relation ( 2) tends (pointwise) to the Bingham model (1) (c.f. ( 5)).Like in the Bingham model, at large enough strain rates the mechanical behaviour of a viscoplastic fluid modelled with (2) resembles that of a Newtonian fluid with viscosity  (c.f. ( 3)).Finally, the requirement (4) imposes that the shear stress-shear strain rate relationship prescribed by a regularised model differs from the one for a Newtonian fluid with viscosity  by a monotonically increasing function of the rate of strain.

Basic equations for steady two-dimensional flows in a symmetric channel
Henceforth, we shall focus on the steady two-dimensional flow of a viscoplastic fluid in a symmetric channel where ℎ() ≠ 0 for all  ∈ [0, ] (Fig. 1).
The equations governing the steady two-dimensional flows of a viscoplastic fluid in  are provided by the equation of the constraint of incompressibility and the components of the equation of balance of linear momentum which, assumed that the velocity field is in the form (, ) = (, ) 1 + (, ) 2 , in the absence of body forces read where  is the (constant) mass density, and the components of the deviatoric stress and rate of strain tensors satisfy the constitutive relations (1) or (2).
A. Farina et al.Clearly, Eqs. ( 11) have to be solved under appropriate boundary conditions.Here, we assume that the walls bounding the channel are impermeable and the fluid does not slip at the boundaries.Thus, we require that (, ±ℎ()) ⋅ (, ±ℎ()) = 0 (impermeability boundary conditions), (12a) and (, ±ℎ()) ⋅ (, ±ℎ()) = 0 (no-slip boundary conditions), where represent, respectively, the unit tangent and outward normal vector fields at the boundaries  = ±ℎ().Note that inserting (13) into (12) leads to the homogeneous boundary conditions on the components of the velocity field As a preliminary result, we observe that the discharge is constant.In fact, in view of the constraint of incompressibility (c.f. the first equation in ( 11)) and the boundary conditions ( 17), we have In what follows the discharge is a datum of the problem and is assumed to be positive.
In view of the homogeneous boundary conditions ( 14) and the symmetry of the channel , it is reasonable to expect that symmetric flows (more precisely, motions for which the velocity field is such that (, ) = (, −) and (, ) = −(, −) for all (, ) ∈ ) occur.This allows us to limit our analysis to the upper part of the channel , providing that we consider the boundary conditions Note that the conditions on  and  at  = 0 are necessary for the symmetry of the velocity field about the -axis.

Lubrication approximation
The boundary value problem ( 11) and ( 17) for a viscoplastic fluid is quite difficult to solve analytically.Approximated solutions can be found when the aspect ratio  of the channel is small by introducing the lubrication approximation.In this regard, setting  = max ∈[0,] ℎ(), we assume that For the subsequent analysis it is convenient to non-dimensionalise the governing Eqs. ( 11), the boundary conditions ( 17) and the constitutive Eqs. ( 1) and ( 2) by taking  = ∕ as the characteristic velocity and introducing the following dimensionless quantities Inserting ( 19) into ( 10), ( 1), ( 2), ( 11) and ( 15) and omitting the tildes for simplicity of notation yield: • the dimensionless parameterisation of the channel where 0 < ℎ() ≤ 1 for all  ∈ [0, 1]; • the dimensionless Bingham model where is the Bingham number and the dimensionless rate of strain tensor reads whence • the dimensionless regularised model 2 • the dimensionless governing equation where  =  ∕ is the Reynolds number; • the dimensionless equation for the half discharge Inserting the dimensionless quantities ( 19) into (17) and, again, omitting the tildes produce dimensionless boundary conditions that are formally identical to their dimensional version.
In what follows we shall limit our analysis to flow regimes for which the Bingham number is of order of unity and the Reynolds number is of order of unity or less.Then, collecting terms of order (1) in ( 23)-( 26) and using the subscript notation for differentiation with respect to the spatial variables  and  we obtain this set of approximate governing equations where, when using a regularised model, the shear stress  12 and the rate of shear strain   are related by When using the Bingham model, collecting terms of like order in (21) reveals that  12 is expressible in terms of   only if the shear stress exceeds in modulus the Bingham number, otherwise the motion is rigid (see Fusi, Farina, Rosso, and Roscani (2015) for a detailed explanation of this point).In formulae, we have From ( 28) we deduce that the Lagrange multiplier  depends only on , and, since the fluid flows from the inlet to the outlet (as the discharge is positive),   is negative.
A. Farina et al.

Approximated flow of a Bingham fluid
Within the lubrication approximation, we determine the flow of a viscoplastic fluid by considering the constitutive model (30).In this case Eqs.(28) do not hold in the unyielded plug region, i.e. the domain   ⊆  where | 12 | ≤  and the motion is rigid.Since the stress in the unyielded plug region is indeterminate, the appropriate equation that governs the motion in   is the equation of balance of linear momentum written in its integral version (Fusi et al., 2015;Safronchik, 1959).To make our analysis as simple as possible and for simplicity of presentation, we assume that the unyielded plug region is connected and extends from the inlet to the outlet of the channel so that   can be parameterised as where  = ±() are the yield surfaces, i.e. the loci at which | 12 | = .Thus, by following similar arguments as in Fusi et al. (2015), using the symmetry of the flow and in view of ( 18), the equation of balance of linear momentum of the unyielded plug region lying in  + approximates to Eqs. ( 26) hold instead for the yielded region, namely the domain  −   in which the shear stress exceeds the Bingham number.Since   is negative in  −   , from the second equation in ( 28) and ( 30) we deduce that also   is negative and thus the shear stress-shear strain rate relationship reduces to the simple relation Again, in view of the symmetry of the yielded region and the boundary conditions, we can limit our analysis to the region . Since   vanishes at the yield surface  = (), from (33) it follows that  12 = − at  = (), whence integrating (49) yields Next, from ( 33) and ( 34) we obtain the differential equation which integrated subject to the first boundary condition in (17) gives Combining this result with the continuity of the velocity field and the fact that the motion in the unyielded plug region is rigid, we deduce that the component of the velocity of the rigid core in the -direction is constant (say  =   ) and given by From ( 27), ( 36) and (37) we deduce that the yet unknown plug-speed   and yield surface  = () are related by the equation whence We now observe that since 0 ≤ () ≤ ℎ() for all  ∈ [0, 1], the plug-speed must satisfy the inequalities Then, we conclude that a necessary condition for the existence of a connected unyield phase extending from the inlet to the outlet is that the profile of the channel is such that In the subsequent analysis we assume that this requirement is satisfied.The restriction (41) is a special case of the necessary condition (85) found by Panaseti, Georgiou, and Ioannou (2019) for the existence of a connected yield phase in a Herschel-Bulkley fluid.This is consistent with the fact that the Herschel-Bulkley model reduces to the Bingham one when the power-law exponent  is equal to unity.
To complete the solution of the problem and determine the flow of the Bingham fluid we have to compute   .To this end, we consider (32), take into account that   is given by the second equation in (39) and  12 (, ()) = , and obtain the equation in the unknown It can be easily proven that the function  is a continuous monotonically decreasing function over the domain [1∕ℎ min , 3∕2].As an immediate consequence, Eq. ( 42) admits a unique solution in the interval [2∕3, ℎ  ] if and only if This means that, within the lubrication approximation, flows of a Bingham fluid in which the yielded region is connected and extends from the inlet to the outlet are possible only in channels with boundaries that satisfy (41) and at Bingham numbers in the range (43).Under these restrictions on the boundaries of the channel and on the Bingham number, from the first equation in ( 26), ( 36), ( 37), ( 39) and the boundary conditions ( 17) we deduce that the components of the velocity fields are as follows and where   =  −1 ().
As an example, in the simple case in which the thickness of the channel is constant, there is no restriction on the Bingham number as, for ℎ ≡ 1, the interval in (43) coincides with [0, +∞[, the unique root (belonging to the interval [1, 3∕2]) of Eq. ( 42) is whence the yield surface is the plane and the components of the velocity field ( 44) and ( 45) read and  ≡ 0. (48)

Approximated flow of a viscoplastic fluid modelled with a regularised stress-rate of strain relationship
Within the lubrication approximation, we now determine the flow in  of a viscoplastic fluid whose mechanical response is characterised by a regularised model in the form (2).
Firstly, we rewrite the shear stress-shear strain relation (29) as As illustrations, when using the generalised Bercovier-Engelman model the dimensionless effective viscosity in (49) reads while for the Papanastasiou model η is given by From ( 3)-( 5), we have that, for any positive , η(⋅; ) ∶ [0, +∞[→]1, 1 + ∕] is a smooth monotonically decreasing function η(⋅; ) such that A. Farina et al. for any  ∈ R, and (54) It is worth pointing out that, in view of the properties of the effective viscosity, for any regularised model and any positive regularising parameter , the shear stress is well defined in terms of the shear strain rate by (49).Only in the limit as  → 0 + ,  12 is indeterminate in the regions where the shear strain rate vanishes (c.f. ( 54)).
As consequences of the properties of the dimensionless effective viscosity η, for any positive  the function  (⋅; ) in ( 49) is defined over the whole real axis, odd, monotonically increasing and such that lim →±∞  (; ) = ±∞.Its inverse  −1 (⋅; ) is such that is differentiable and, as a direct consequence of ( 53), its first derivative with the prime denoting differentiation with respect to ||, is positive.
Moving to the solution of ( 28), since   is negative, from ( 17) and ( 49) we have that, for any  ∈ [0, 1],   (, ) is negative for all 0 <  ≤ ℎ() and vanishes at  = 0.These deductions on the sign of the shear strain rate and the properties of  (⋅; ) allow us to determine the component of the velocity field along  1 ,

𝑢(𝑥, 𝑦
The gradient of the Lagrange multiplier in ( 57) is however as yet unknown.To determine it, we combine the dimensionless equation for the discharge ( 27) with ( 57) and obtain the equation in the unknown Eq. ( 58) admits a unique solution.To prove this assertion, for any  ∈ [0, 1] we introduce the function From ( 55) we deduce that lim →0   (; ) = 0 and lim   (⋅; ) is differentiable and, as a consequence of (53), its first derivative is negative.In summary, for any  ∈ [0, 1],   (⋅; ) is continuous, monotonically decreasing and satisfies (60).Thus, for any  ∈ [0, 1] there exists a unique p (; ) that satisfies the equation   (; ) = 1.The existence and uniqueness of the solution of Eq. ( 58) is thus proven.Next, the differentiability of  −1 (⋅; ) implies the differentiability of p (⋅; ), and implicit differentiation of (58) yields that In view of the properties of  (⋅; ) and (56),  (⋅; ) is positive in [0, 1].Thus, we deduce that sign( p (; )) = sign(ℎ  ()), and p (; ) vanishes if and only if ℎ  () vanishes.In particular, p (⋅; ) is identically constant if and only if ℎ ≡ 1, that is for plane channels.Finally, inserting p (⋅; ) into (57) gives the solution for the component  of the velocity field, which, combined with the incompressibility constraint (c.f. the first equation in ( 26)), yields the solution for the component of the velocity field along  2

𝑣(𝑥, 𝑦
We now define the pseudo-yield surface (that lies in the upper part of the channel)  = σ(; ) as the locus at which the modulus of the shear stress equals the Bingham number, i.e | 12 (, σ(; ))| = .From ( 49) and ( 57), the pseudo-yield surface is found to have equation If the fluid flows in a plane channel (for which ℎ ≡ 1), the pseudo-yield surface ( 64) is actually a pseudo-yield plane (c.f. ( 62)).Recalling ( 46) and ( 48), the yield surface in a Bingham flow in a flat channel is the plane (47).As it is reasonable to expect, the regularised and exact yield planes differ because Eq. ( 64) does depend on the specific regularising model used.However, as we shall prove in the next section, irrespective of the regularised model used, these two planes coincides in the limit as  → 0 + (see Fig. 2).On the contrary, if the viscoplastic fluid flows in a non-plane channel, the pseudo-yield surface (64) differs significantly from the yield surface  = 3∕  − 2ℎ() predicted by the Bingham model.Suffice to observe that sign(  ) = −sign(ℎ  ) in [0, 1], whereas, in view of ( 62) and ( 64), sign( σ (; )) = sign(ℎ  ()) for all  ∈ [0, 1].In simple terms, the exact yield surface  = 3∕  − 2ℎ() has the opposite monotonicity of the profile of the upper wall of the channel, whereas the pseudo-yield surface (64) has the same monotonicity as the boundary  = ℎ() (see Fig. 3).

Asymptotic results for regularised models
The analysis carried out in the previous section reveals that, within the lubrication approximation and by using a regularised model in the form (2), the approximate equations governing the two-dimensional flows in a symmetric channel admit a unique solution for any profile of the channel walls and at any Bingham number.On the other hand, since any regularised model satisfies the limit (54), it is natural to investigate the behaviour of the solutions ( 57) and ( 63) in the limit as  → 0 + .
To perform such an asymptotic analysis, we first observe that, in view of ( 53), it can be proven that • in the limit as  → 0 + ,  −1 (⋅; ) converges uniformly to a piecewise continuously differentiable function, as  → 0 + ; • for all  > 0, sup ∈R d −1 d (; ) < 1; (66) • in the limit as  → 0 + , one has (67) Next, we define the new variable () = ∕|  ()|, substitute |  ()| by ∕() into (58), and from (65) deduce that in the limit as  → 0 + , Eq. ( 58) reduces to For any Bingham number and any  ∈ [0, 1], Eq. ( 68) admits as unique solution in the interval [0, ℎ()].Again in light of (65), we deduce that Combining ( 70) with ( 57), ( 63), ( 65), ( 67), ( 66) and ( 70), and by means of the dominated convergence theorem, we conclude that, in the limit as  → 0 + , (⋅, ⋅; ) converges uniformly to Here, we have assumed that the channel is plane and  = 3.For  = 0.01 the pseudo-yield surfaces predicted by the two regularised models practically coincide with the asymptotic surface (47) which, in turn, is identical to the yield surface predicted by the Bingham model.and (⋅, ⋅; ) converges uniformly to The locus  =  0 () and the asymptotic solutions  0 and  0 depend exclusively on the profile of the boundaries of the channel as they do not keep track of the particular regularised model used.These asymptotic results are then universal for the class of regularised models in the form (2). Incidentally, the asymptotic flows in a channel obtained by Housiadas and Georgiou (2023) when using the Papanastasiou model coincide with ( 71) and ( 72).
We finally observe that the shear stress is indeterminate in the region as the shear strain rate  0 vanishes there.However, since the velocity field in  0 is constant only if ℎ ≡ 1, the motion of the part of the fluid that occupies  0 is rigid only if the channel is plane.In this case,  0 coincides with the unyielded plug region   predicted by the Bingham model.If the viscoplastic fluid flows in a non-plane channel,  0 and  0 are not constant in  0 as they depend on the spatial variable  (c.f. ( 71) and ( 72)), whence, on using the same terminology as that adopted by Frigaard and Ryan (2004),  0 is a pseudo-plug region as it is not truly unyielded, and  =  0 () is a pseudo-yield surface.

Comparison between the results predicted by the Bingham and regularised models
We are now in position to compare the results predicted by using the Bingham and regularised models within the lubrication approximation.
When considering the Bingham model, two-dimensional flows of a viscoplastic fluid with a connected unyielded plug region extending from the inlet to the outlet of the channel are admissible providing that the boundaries of  satisfy the restrictions (41) and the Bingham number lies in the interval (43).On the contrary, when using a regularised model, the boundary value problem ( 17) and ( 26) is uniquely solvable for any profile of  and any value of the Bingham number, no part of the fluid moves rigidly, and the region bounded by the pseudo-yield surfaces  = ± σ( ∶ ) (c.f. ( 64)) is connected and extends from the inlet to the outlet.
On the other hand, in the limit as the regularising parameter  tends to zero, despite the convergence of any regularising model to the Bingham model, the asymptotic analysis carried out in the previous section confirms the absence of parts of the fluid that move rigidly if the channel is non-flat (c.f. ( 71) and ( 72)).This means that, no regularised model in the form (2) is appropriate to reproduce the steady two-dimensional flows of a Bingham fluid in a non-plane channel.On the contrary, for plane channels (for which ℎ ≡ 1) the asymptotic surface  =  0 (), with  0 as in (69), and the yield surface predicted by the Bingham model coincide and have Eq.( 47).Moreover, the asymptotic solutions ( 71) and ( 72) coincide with (48).We may then conclude that, within the lubrication approximation and in the limit as  → 0 + , any regularised model in the form (2) reproduces precisely the two-dimensional flows of a Bingham fluid in a plane channel.It might sound curious that the accuracy of a regularised model to reproduce the two-dimensional flows of a Bingham fluid does depend on the shape of the channel.However, this is a consequence of the fact that, in the limit as  → 0 + , a viscoplastic fluid modelled with a regularised stress-rate of strain relationship presents a rigid core like a Bingham fluid only in plane channels.
To analyse further the differences between the results predicted by using the Bingham and regularised models we have solved numerically the system of governing Eqs.(26) without introducing the lubrication approximation.The numerical simulations have been performed by using the CFD model built in Comsol Multiphysics ® .Specifically, we have considered a channel with aspect ratio  = 0.01 (c.f. ( 18)) and taken the half discharge  and the constitutive parameters  and  0 so that  = 0.1 and  = 3. Taking advantage of the symmetry of the channel and the flow, we have discretised  + with 10000 4-node quadrilateral elements.With this scheme, we have performed numerical simulations for plane, sinusoidal, divergent and convergent straight-walled channels by using the Papanastasiou (7), simple (8) and Bercovier-Engelman (9) models, and considering different values of the (dimensionless) regularisation parameter .In each simulation we have determined the surface  =   (; ) at which the shear stress  12 equals in modulus the Bingham number and recorded the uniform norm ‖  (⋅; ) −  0 ‖ ∞ ( Table 1).What emerges is that for any type of the channels considered ‖  (⋅; ) −  0 ‖ ∞ decreases as  decreases, and for small enough values of  the two surfaces  =   (; ) and  =  0 () almost coincide (see Fig. 4).This means that, for a given profile of the upper boundary of  + , in the limit as  → 0 + also the numerical approach based on the solution of the full system of governing Eqs. ( 26) and boundary conditions (17) leads to asymptotic surfaces that depend only on the shape of the boundary of the channel, and are independent of the specific regularised model used.This means that for sinusoidal, divergent and convergent channels the surfaces  =   (; ) differ significantly from the yield surface predicted by the Bingham model.Only in a plane channel the surface  =   (; ), with  small enough, practically coincides with the yield plane (47) predicted by the Bingham model (see Fig. 4(d)).

Conclusions
We have studied the two-dimensional flows of a viscoplastic fluid in a symmetric channel .As constitutive functions for the deviatoric stress tensor, we have considered both the classical Bingham model (1) and the very general class of smooth models in the form (2) that, at small values of the dimensionless regularisation parameter  (c.f. ( 25)), approximate the material response of those fluids that, under specific experimental conditions, seem to exhibit a yield stress.For simplicity, we have assumed that the aspect ratio of the channel (defined as the ratio between the maximum half width and the length of ) is very small so that it makes sense to introduce the lubrication approximation.In so doing, we have been able to solve analytically the approximated boundary value problem that governs symmetric flows in  both when the mechanical response of the viscoplastic fluid is modelled by the Bingham model and when the stress-rate of strain relationship is given by a smoothed model.In both cases, we have determined the pseudo-yield surfaces, i.e. the non-material surface at which the modulus of the shear stress equals the yield stress, and observed that at small  the pseudo-yield surfaces predicted by regularised Bingham models are practically identical to the exact yield surface if the channel is plane, whereas they are markedly different when the walls of  are not flat.Specifically, in non-plane channels, the yield surface (that lies in the upper part of the channel) predicted by the Bingham model has the opposite monotonicity of the profile  = ℎ() of the upper boundary of the channel (this is in perfect agreement with the results by Frigaard and Ryan (2004), Muravleva andMuravleva (2011) andFusi et al. (2015)), whereas the pseudo-yield surface predicted by the regularised model has the same monotonicity as the upper wall of .Such a discrepancy has been confirmed also by the asymptotic analysis carried out in Section 7 and the numerical investigation in Section 8 performed without introducing the lubrication approximation.This result is very important since it shows that the regularised models are not able to reliably reproduce the exact yield surface in complex geometries.Specifically, in non-flat channels, the loci at which the shear stress equals the Bingham number determined numerically by using CFD models that use regularisations of the Bingham model do not represent the exact yield surfaces.
Finally, the asymptotic analysis has revealed that, irrespective of the particular regularised model adopted, in the limit as  → 0 + the asymptotic pseudo-yield surface and velocity field depend exclusively on the profiles of the boundaries of the channel and the value of the Bingham number.These universal results for the class of regularised models have not been obtained before.The numerical simulations in Section 8 performed by using the simple, Bercovier-Engelman and Papanastasiou models show that, for  = 0.01, the pseudo-yield surfaces almost coincide with the asymptotic pseudo-yield surface  =  0 (), with  0 as in ( 69).This

Fig. 1 .
Fig. 1.Sketch of the channel in which the viscoplastic fluid flows.

Fig. 2 .
Fig. 2.Pseudo-yield surfaces predicted by the simple and Papanastasiou models at different values of the regularisation parameter .Here, we have assumed that the channel is plane and  = 3.For  = 0.01 the pseudo-yield surfaces predicted by the two regularised models practically coincide with the asymptotic surface (47) which, in turn, is identical to the yield surface predicted by the Bingham model.

Fig. 3 .
Fig.3.Pseudo-yield surfaces predicted by the simple and Papanastasiou models at different values of the regularisation parameter .Here, we have assumed that the upper boundary of the channel has equation  = 0.9 + 0.1 sin(2) and  = 3.For  = 0.01 the yield surfaces predicted by the two regularised models coincide with the asymptotic yield surface (47).The yield surface predicted by the Bingham model differs markedly from the asymptotic pseudo-yield surface, and has the opposite monotonicity of the profile of the upper wall.

Fig. 4 .
Fig. 4. Yield surfaces predicted by the Bingham model for sinusoidal (a), divergent straight-walled (b), convergent straight-walled (c) and plane channels.Except for the case in which the thickness of the channel is constant (i.e.ℎ ≡ 1), the yield surfaces differ significantly from the asymptotic surface  =  0 (), with  0 as in (39), to which the pseudo-yield surface predicted by any regularised model tend in the limit as  → 0 + , and from the surface  =   (; 0.01) obtained numerically by using the Papanastasiou model.The latter two surfaces are almost identical in all the four types of channel (the solid black and red curves coincide).The contour plots in the background of the four figures refer to the modulus of the shear stress evaluated numerically.