A LUBRICATION EQUATION FOR A SIMPLIFIED MODEL OF SHEAR-THINNING FLUID

Abstract. A lubrication equation is obtained for a simplified shear-thinning fluid. The simplified rheology consists of a piecewise linear stress tensor, resulting in a two-viscosity model. This can be interpreted as a modified Bingham fluid, which can be recovered in a specific limit. The lubrication equation is obtained in two steps. First two scalings are performed on the incompressible Navier-Stokes equations, namely the long-wave scaling and the slow motion scaling. Second, the resulting equations are averaged along the vertical direction. Numerical illustrations are provided, bringing to light the different possible behaviours.


Introduction
The lubrication equation is quite a classical simplification of the incompressible Navier-Stokes system. It is obtained for thin films of fluid, when viscous effects balance the pressure force. This occurs for instance for thin films of oil, hence the name of the equation. The study of this approximation goes back to Reynolds in 1886 [9]. Several scalings are involved to obtain this model. First the aspect ratio between the thickness of the film and the characteristic length of the substrate must be small, say δ. Simultaneously, the time scale has to be of order 1/δ. This is the so-called long wave regime, and is classically used in the shallow-water approximation. The lubrication equation requires another assumption of balance between the viscous effects and the pressure effects, which amounts to neglect all kinematic effects. This simplified flow is known as the Stokes flow. The lubrication equation itself is then obtained by integration over the fluid thickness.
We are interested here in the lubrication model for a class of non Newtonian fluids. Several fluids are known to depart from the usual Newtonian rheology, where the deviatoric stress tensor is a linear function of the strain rate tensor, thus defining the dynamical viscosity of the fluid. The lubrication equation for Newtonian fluids has been studied for instance by Huppert [5]. Non Newtonian fluids arise in several applications in engineering, biology, geophysics... In particular, viscoplastic or pseudoplastic fluids are involved in various geological problems, for instance lava flows, mudslides and avalanches. We refer to [1] for a review on the subject. A model which is widely used is the so-called Bingham-plastic model. This model involves a yield stress, namely a threshold on strain rate: for values of the strain rate above this threshold the fluid behaves like a viscous fluid, for values below, it looks like a solid. This can be thought of as an infinite viscosity fluid. We refer to the papers by Liu and Mei [7] and Balmforth et al. [2] for the study of such fluids in the lubrication approximation. Both papers contain also a complete bibliography. Liu and Mei also introduced in [8] a perturbed Bingham model, which is actually a two viscosities model, with a high viscosity for small deformations. When this viscosity goes to ∞ the Bingham model is recovered, thus giving a fluid mechanics interpretation of this solid behaviour. This is precisely the two viscosities model we investigate here. First we describe the mathematical model we use, namely the incompressible Navier-Sokes equations in a time-dependent domain, since we consider a free-boundary problem. In particular we explain in some details all the scalings involved. Next, we turn to the lubrication equation itself, which is a one-dimensional equation, obtained by averaging the previous ones along the vertical thickness. Finally we provide a few numerical illustrations based on a finite volume scheme.

Mathematical model
In this section we set up the model. The starting point is the incompressible Navier-Stokes system. We limit ourselves in this paper to the two-dimensional case, thus aiming at a one-dimensional lubrication equation. Similar computations can be performed in three space dimensions. The domain we consider is Ω t defined by f b (x) < z < ϕ(t, x), for t > 0 and x ∈ (−∞, +∞), where f b is given topography, and ϕ is a free surface. The notations we use are gathered in Figure 1. Figure 1. Notations for the two viscosities fluid: ϕ is the free surface; f b is the topography of the substrate; z * is the ordinate which separates "small deformations" (white zone) from "large deformations" (green zone), see Section 3 below. We introduce the thicknesses h The incompressible Navier-Stokes equations are where ρ is the density of the fluid, U = (u, w) is the velocity field, and the stress tensor σ is written as the sum of a volumetric stress tensor, involving the pressure p, and a deviatoric stress tensor τ : where I 2 is the identity matrix in dimension 2. The density ρ is assumed to be constant here, and the tensor σ will be defined in Section 2.1 below. Boundary conditions are: z = ϕ: fluid-atmosphere interface. We have continuity of the stress tensor at the free surface, together with a kinematic boundary condition. Since the atmosphere can be viewed as an ideal fluid, the stress tensor can be taken equal to zero above ϕ. Hence we get z = f b : interface between the fluid and the substrate, which has a fixed shape. This is a material interface, on which we impose a no-slip boundary condition Here (u b , w b ) is the so-called basal velocity. The classical no-slip condition in fluid mechanics corresponds to a zero basal velocity. A non zero velocity is motivated by an application in geophysics, where f b is the interface between the earth's crust and mantle. In a crude modelling, f b has a fixed shape, and the basal velocity is the trace of the convection movements of the mantle.

Rheology
For a fluid, the deviatoric stress tensor τ is usually a function of the strain rate tensoṙ A Newtonian fluid is characterized by a linear relation, defining the viscosity of the fluid, which is assumed here to be isotropic and constant. Therefore we introduce the dynamical viscosity coefficient µ, and define the Newtonian stress tensor by τ N = 2µε = 2ρνε, where ν = µ/ρ is the kinematic viscosity.
Fluids that do not follow this kind of constitutive law are non-Newtonian. In the general case, the material invariance principle implies that the stress tensor depends only on the similarity invariants of the strain rate tensor, in particular the coefficients of its characteristic polynomial. In dimension 2 there are only two such coefficientsε I andε II . Namelyε I is the trace of the matrix andε II its determinant. For an incompressible fluid, the trace is zero, and moreover we havė This allows to define the strain rateγ aṡ In a similar way we can check that the Frobenius norm ofε, that is ε 2 = i,j (ε ij ) 2 , satisfies A very sketchy illustration of the possible behaviours of non-Newtonian fluids is given in Figure 2. We will be mostly interested in this work in the so-called pseudoplastic case, that is the red curve in Figure 2, for which experimental evidence can be given, see [4]. This kind of models are also used in geophysics, see [3,10,12]   We wish to give a simplified model for this pseudo-plastic fluid, that allows to handle explicit computations. The main feature of this kind of fluids is a nonlinear viscosity, decreasing with the strain rate. Mimicking the Bingham model, which is based on a threshold on the shear stress, we consider a model with a threshold on the strain rate: the viscosity is equal to some large µ B for small deformations, that isγ < γ c , where γ c > 0 is a given constant, and to another value µ for large deformations,γ > γ c . Such models were introduced by Liu and Mei [8], and the limit case ν B → ∞, which leads to a Bingham fluid, is studied in [7] and [2]. Notice that using (8) the threshold γ c onγ can be replaced by a threshold γ c = γ c / √ 2 on ε . A multidimensional formulation for these simplified pseudo-plastic fluids is therefore A particular limit case is ν B → ∞, which leads to a Bingham type fluid. To view this, it is convenient to define the following quantities (see Figure 3 below for an illustration in dimension 1) so that definition (9) can be rewritten It is clear on this formulation that the relevant limit is ν B → +∞, together with γ c → 0, keeping ν B γ c = τ c . In doing so, we recover the classical Bingham stress tensor, with threshold τ c : Finally, notice that in the pseudo-plastic (or shear thinning) context, we consider 0 < ν < ν B , but similar computations can be performed in any case. It is convenient for the scalings below to rewrite expression (11) using an equivalent kinematic viscosity ν eq , which satisfies ν ν eq ν B :

Scalings
We introduce now the scaling laws, namely thin layer, or more precisely long wave approximation, and slow motion, in order to finally obtain the lubrication model. This kind of scalings is already present e.g. in [2] in the context of a visco-plastic fluid. Hence we propose the following family of scalings: we introduce a first set of characteristic scales, namely dimensions 0 and h 0 , characteristic velocities u 0 and v 0 , and a characteristic time t 0 . The quantities 0 and u 0 correspond to the horizontal direction, h 0 and v 0 to the vertical one. The aspect ratio δ = h 0 / 0 will be an important parameter, assumed to be small in the thin layer case. Dimensionless variables are then defined by First, we rewrite the incompressibility equation (1) in the rescaled variables. We obtain and following the least degeneracy principle [11], this implies u 0 / 0 = w 0 /h 0 , or equivalently 0 /h 0 = u 0 /w 0 . Thus w 0 /u 0 = δ, so that in the thin layer approximation w 0 is also small compared to u 0 . We turn now to the kinematic part of the equation. Using u 0 / 0 = w 0 /h 0 , we readily obtain Once again we apply the least degeneracy principle and obtain t 0 = 0 /u 0 = h 0 /w 0 , or, as expected, u 0 = 0 /t 0 and w 0 = h 0 /t 0 . We proceed in the same way for the momentum equation in v and finally obtain where we have emphasized the aspect factor δ = h 0 / 0 = w 0 /u 0 , and highlighted as characteristic quantities the horizontal velocity u 0 and the vertical extension h 0 . Following Balmforth [2], we rescale the pressure and the stress tensor by where we recall that ν is the kinematic viscosity for large deformations. We can write now the rescaled version of the Navier-Stokes momentum equations (2) and (3): At this stage, we introduce two classical dimensionless quantities, namely the Froude and Reynolds numbers We divide the previous two equations by u 2 0 /h 0 , and noticing that τ = ρReu 2 0τ , we obtain The idea now is to send δ to zero, thus implementing the thin layer assumption, but in a regime where the Reynolds number Re is kept of order 1, together with a balance between viscosity and gravity forces. Therefore we set F r 2 = δRe, This readily gives The first equality introduces another characteristic velocity, namely u 0 = (gh 2 0 )/ν, and shows that this is indeed a slow motion scaling, thus we meet the initial requirement. The latter relation is the scaling proposed in [2].
Inserting (21) in equations (19) and (20), and keeping only the dominant terms of order δ −1 gives first the dimensionless Stokes equation then the dimensionless hydrostatic relation for the pressure Now we computeτ from (12). We start by rewritingε in rescaled variableṡ From this we easily deduceτ We define a dimensionless equivalent viscosity byν eq = ν eq /ν, and rewrite equation (23) ∂z ν eq ∂ zū = ∂xp.
We turn now to the expression ofν eq . We first notice that, using (25) Hence the condition ε > γ c leads us to define a dimensionless thresholdγ c = ( We introduce a dimensionless viscosityν B and a dimensionless yield stress B by settinḡ so that the expression ofτ xz , which is the only part of the deviatoric stress tensor remaining in the equations, becomes (see Figure 3)τ This is the model proposed by Liu and Mei in [8].γ  Figure 3. Simplified shear-thinning model. We consider a piecewise linear approximation (in black) of the "theoretical" pseudoplastic law (in red). Parameters τ * and τ c are defined by (10). The blue curve is the Bingham limit: The green dashed line is the pure Newtonian limit ν B → ν.
As concerns the boundary conditions, we notice that the no-slip and kinematic boundary conditions remain unchanged by the scaling. In contrast, the continuity of the stress tensor across the free surface ϕ is greatly simplified. Recalling that ϕ = h 0φ , we indeed obtain Now making use of (26), we obtain Letting δ go to zero implies therefore that (4) becomes In other words, we recover separately the continuity of ∂ z u and the continuity of the pressure.

Remark 1.
It is worth to notice here that we assume implicitly in this paper that ν eq is bounded. However, if we wish to consider ν B → +∞, we should take care of the product δν eq that appears at several places in the equations. Thus another scaling arises, namely δν B should go to zero when δ goes to zero. This was pointed out by Liu and Mei [8].

Lubrication equation
The so-called lubrication equation is obtained by integrating equations (1) along the vertical direction. The long wave and slow motion assumptions imply that we obtain a single nonlinear equation on the depth ϕ. Similar computations were performed by Liu and Mei [8], for a two-viscosity model, in order to justify the Bingham case, which corresponds to ν B → ∞ in our context. For Bingham fluids, we refer to Liu and Mei [7], and more recently to Balmforth [2]. The final equation is obtained through three steps we present in detail now.
We recall the equations we obtained in the preceding section, dropping the bars for clarity. First we have the hydrostatic relation next, the dimensionless Stokes equation (23) where τ xz is the dimensionless deviatoric stress tensor defined by (31). These equations are coupled with the following boundary conditions (in these relations, t and x are hidden parameters): • on the free surface z = ϕ p(ϕ) = 0, ∂ z u(ϕ) = 0, • Concerning first the pressure, using the boundary condition on the free surface we obtain the usual hydrostatic approximation The averaged equation we look for is obtained by integrating in z the incompressibility equation, or mass conservation, ∂ x u + ∂ z w = 0. This is quite classical, see e.g. [2] in the same slow motion context, or [6] for shallow water approximation. We obtain The kinematic boundary condition on z = ϕ leads to w(t, x, ϕ) − u(t, x, ϕ)∂ x ϕ = ∂ t ϕ = ∂ t h, and for z = f b , we make use of the no-slip boundary condition (36), to obtain the following averaged equation The flux The first step towards the computation of the flux is to obtain the vertical velocity profile. The general structure of this profile is as follows. We have τ xz = F (∂ z u), where F is a continuous, one-to-one, increasing function, with F (0) = 0, see (31) and Figure 3. From (34) and (37) we are led to solve . Because F is increasing, ∂ z u is monotone as well, in particular, since ∂ z u = 0 for z = ϕ, its sign remains constant. Therefore |∂ z u| is decreasing in z (increasing with depth).
The threshold in formula (31) eventually splits the fluid in two layers. Let z * be defined by |∂ z u(z * )| = γ c . Provided z * ∈]f b , ϕ[ (see below for precise formulas), starting from the surface ϕ where |∂ z u| = 0 we first encounter a "small deformation" region, where |∂ z u(z)| < γ c , for z ∈]z * , ϕ[, because |∂ z u| is increasing with depth. Similarly for z ∈]f b , z * [ we have |∂ z u(z)| > γ c . Finally, according to (31), the velocity is ruled by the system of equations where for the second equation we have used that ∂ z u has a constant sign. These equations are complemented with the boundary conditions ∂ z u = 0, z = ϕ ; u = 0, z = f b . Notice that the curve z = z * is not a physical interface, yet we have continuity of the stress tensor, or equivalently here continuity of ∂ z u. Now the computations are quite easy. We integrate once the first equation between ϕ and z * , to obtain

This leads toγ
so that the value of z * and the thickness h * of this layer are given by These definitions ensure that z * f b and h * h, and are valid for ∂ x ϕ = 0 with the convention B/0 = ∞. Notice that z * can be equal to f b for weak slopes (small ∂ x ϕ), or small depths (small h). Conversely, z * → ϕ when |∂ x ϕ| goes to ∞.
Integrating once again between z * and ϕ, we obtain the velocity profile for ϕ z z * : where the constant K will be determined later. Notice for further use that by construction We turn now to the lower layer, z * z f b . The fluid here has dimensionless viscosity 1, and we use the boundary conditions (42) for z = z * , and no slip (36) at z = f b . First we get, using (42), next, integrating once again between f b and z * , where L is computed using (36), leading to Finally, we use the continuity of the velocity at z = z * to obtain the constant K: The velocity profile is therefore given by Notice that for ν B = 1, easy computations show that the profile is the same in the two layers, namely u = , which is as expected the usual parabolic profile for a Newtonian fluid.
On the other hand, letting ν B → +∞, and γ c → 0, keeping ν B γ c = τ c , we recover formally the Bingham fluid velocity, as in Balmforth [2]: It is now straightforward to obtain the flux in (40), since We have on the one hand Therefore the flux we are looking for is given by It is easy once again to check on this formula that we recover the usual cubic flux − ∂ x ϕ 3 h 3 for the Newtonian fluid ν B = 1. On the other hand the limit case ν B → ∞ gives back Balmforth's formula

Inserting (45) in the conservation equation (40) leads to the following advection-diffusion equation:
where and we recall the definitions of h * from (41), and h * Notice that 0 < D(h, ∂ x h) h 3 /(3ν B ).

Numerical illustrations
We turn now to numerical examples to illustrate the behaviour of the two-viscosity fluid. The point here is not to give an accurate specific scheme, which is by all means an interesting perspective since the diffusion term may degenerate, but is beyond the scope of this work. We merely apply here a simple finite volume strategy. The infinite space domain is replaced by some finite computational domain [a, b]. Since we do not want to cope with boundary conditions here, we merely impose a free flux on the boundaries, which is compatible with the examples we choose. Positive time and space steps ∆t and ∆x being given, we introduce the usual notations t n = n∆t, n 0, and x j = j∆x, 0 j J, where J = (b − a)/∆x. An approximation of the depth h is sought for in the form h n+1 is the numerical advection flux, and G n j+1/2 the numerical diffusion flux, both computed at interface x j+1/2 . In the following we denote u n j the discretized basal velocity, and f j the discrete topography, which are both given functions.
The advection flux is merely an upwind flux For the diffusive flux, we write G n j+1/2 = D n j+1/2 K n j+1/2 , where K n j+1/2 is the approximate value of the slope and D n j+1/2 is a discretization of (47). To obtain it we need to compute h * and h * at the interface x j+1/2 . Accordingly to (48), we put if not and (h * ) n j+1/2 = (h n j+1 + h n j )/2 − (h * ) n j+1/2 , so that D n j+1/2 is given by The time step ∆t is actually updated at each time step using the CFL condition The following simulations have been performed with J = 200 cells in the interval [−1, 1], together with σ = 0.9. All figures are gathered at the end of the paper. The first set of simulations concerns the collapse of a square-shaped stack on a horizontal flat bottom: h 0 (x) = 1 for x ∈] − 1/3, 1/3[, 0 elsewhere, with zero basal velocity (u b = v b = 0). We first propose a comparison between the two viscosities model and the high viscosity and low viscosity models. The small deformation viscosity is ν B = 100 (recall that ν = 1), and the yield stress is 0.1 in Figure 4, and 0.5 in Figure 5. These figures are complemented by Figure 6 where we display for four values of the yield stress B a timelapse of the evolution of both the total thickness of the fluid h (plain lines) and the thickness of the low velocity layer (dashed lines).
For B = 0.1, the fluid clearly behaves similarly as the low viscosity fluid in the early stages, then eventually it slows down, when the low viscosity layer tends to disappear, see Figure 6, top left. With a yield stress B = 0.5, the two viscosities model stays inbetween the other two, as expected, faster than the high viscosity model, slower than the low viscosity one, see Figure 5. However, one can check that the front hardly moves between t = 10 and t = 50, indicating that the fluid tends to behave as the high velocity one. This is made more explicit in Figure 6, top right, where for t = 10 and t = 50 the low viscosity layer is very small. In general, the thickness h * decreases with time, faster when B is larger. It is hardly observable for t = 50 when B = 2.5, indicating that the fluid is almost completely driven by the high viscosity.
Using the same initial data, we check the convergence of the two viscosities model towards the Bingham fluid when ν B goes to ∞. We take a yield stress B = 1.25, and ν B = 10, 100, 1000. As expected, the behaviour becomes close to the Bingham fluid, yet it departs from it for larger times, see Figure 7. This somehow justifies a posteriori that we have taken into account the scaling δν B → 0, see Remark 1.
We turn now to a different context, closer to the situation in geophysics. The flow here is no longer purely gravity driven, it is actually dragged along by a non zero basal velocity. The idea here is that our pseudo-plastic fluid is a very crude model of some planetary lithosphere, below which lies the mantle. The basal velocity is the upper trace of convection currents in the mantle, which are supposed to be the main drivers of plate tectonics. The initial thickness is constant equal to 1, and we use two basal velocities These velocities crudely correspond to the vertical motion of a magma bubble, which generates local perturbations of the velocity in the mantle. The first one corresponds to some bubble lift, with negative velocity on the left and positive on the right. It generates some kind of a valley surrounded by mountains, see Figure 8. Conversely, the descent of a bubble reverses the velocities, and produces a mountain surrounded with valleys, Figure 9. We notice in both cases that the small viscosity model has very little influence on the time evolution, and that the two viscosities model clearly seems to stop the motion, and leads to rather sharp angles in the thickness.