Critical yield numbers and limiting yield surfaces of particle arrays settling in a Bingham fluid

# Critical yield numbers and limiting yield surfaces of particle arrays settling in a Bingham fluid

José A. Iglesias jose.iglesias@ricam.oeaw.ac.at Johann Radon Institute for Computational and Applied Mathematics (RICAM)
Austrian Academy of Sciences, Linz, Austria.
Gwenael Mercier gwenael.mercier@ricam.oeaw.ac.at Johann Radon Institute for Computational and Applied Mathematics (RICAM)
Austrian Academy of Sciences, Linz, Austria.
Otmar Scherzer otmar.scherzer@univie.ac.at
###### Abstract

We consider the flow of multiple particles in a Bingham fluid in an anti-plane shear flow configuration. The limiting situation in which the internal and applied forces balance and the fluid and particles stop flowing, that is, when the flow settles, is formulated as finding the optimal ratio between the total variation functional and a linear functional. The minimal value for this quotient is referred to as the critical yield number or, in analogy to Rayleigh quotients, generalized eigenvalue. This minimum value can in general only be attained by discontinuous, hence not physical, velocities. However, we prove that these generalized eigenfunctions, whose jumps we refer to as limiting yield surfaces, appear as rescaled limits of the physical velocities. Then, we show the existence of geometrically simple minimizers. Furthermore, a numerical method for the minimization is then considered. It is based on a nonlinear finite difference discretization, whose consistency is proven, and a standard primal-dual descent scheme. Finally, numerical examples show a variety of geometric solutions exhibiting the properties discussed in the theoretical sections.

## 1 Introduction

In this article, we investigate the stationary flow of particles in a Bingham fluid. Such fluids are important examples of non-Newtonian fluids, describing for instance cement, toothpaste, and crude oil [24]. They are characterized by two numerical quantities: a yield stress that must be exceeded for strain to appear, and a fluid viscosity that describes its linear behaviour once it starts to flow (see figure 1).

An important property of Bingham fluid flows is the occurrence of plugs, which are regions where the fluid moves like a rigid body. Such rigid movements occur at positions where the stress does not exceed the yield stress.

In this paper we consider anti-plane shear flow in an infinite cylinder, where an ensemble of inclusions move under their own weight inside a Bingham fluid of lower density, and in which the gravity and viscous forces are in equilibrium (cf (5)), therefore inducing a flow which is steady or stationary, that is, in which the velocity does not depend on time. For such a configuration, we are interested in determining the ratio between applied forces and the yield stress such that the Bingham fluid stops flowing completely. This ratio is called critical yield number.

#### Related work.

To our knowledge, the first mathematical studies of critical yield numbers were conducted by Mosolov & Miasnikov [20, 21], who also considered the anti-plane situation for flows inside a pipe. In particular, they discovered the geometrical nature of the problem and related the critical yield number to what in modern teminology is known as the Cheeger constant of the cross-section of the region containing the fluid. Very similar situations appear in the modelling of the onset of landslides [13, 15, 12], where non-homogeneous coefficients and different boundary conditions arise. Two-fluid anti-plane shear flows that arise in oilfield cementing are studied in [10, 11]. Settling of particles under gravity, not necessarily in anti-plane configurations is also considered in [16, 23]. Finally, the previous work [9] also focuses in the anti-plane settling problem. There, the analysis is limited to the case in which all particles move with the same velocity and where the main interest is to extract the critical yield numbers from geometric quantities. In the current work we lift this restriction and focus on the calculations of the limiting velocities, also from a numerical point of view. Various applications of the critical yield stress of suspensions are pointed out in [4, Section 4.3].

#### Structure of the paper.

We begin in Section 2 by recalling the mathematical models describing the stationary Bingham fluid flow in an anti-plane configuration, and an optimization formulation for determining the critical yield number.
Next, in Section 3 we consider a relaxed formulation of this optimization problem, which is naturally set in spaces of functions of bounded variation, and show that the limiting velocity profile as the flow stops is a minimizer of this relaxed problem.
In Section 4, as in the case of a single particle [9], we prove that there exists a minimizer that attains only two non zero velocity values.
Finally, in Section 5 we present a numerical approach to compute minimizers. This approach is based on the non-smooth convex optimization scheme of Chambolle-Pock [8] and an upwind finite difference discretization [7]. We prove the convergence of the discrete minimizers to continuous ones as the grid size decreases to zero. We then use this scheme to illustrate the theoretical results of Section 4.

## 2 The model

We consider a Bingham fluid filling a vertical cylindrical domain and a solid inclusion . We denote by the portion of the domain occupied by the fluid, and by the corresponding constant densities. We focus on a vertical stationary flow, meaning that the solid flows down at a vertical velocity which is constant in time. Moreover, all quantities are invariant along the vertical direction, so we can directly consider a scalar velocity ( is the velocity of the fluid on and of the solid in ), see Figure 2.

The constitutive equations write

 {div^τ=pz+ρfgon ^Ωf,div^τ=pz and ∇^ω=0 on ^Ωs, (1)

with the pressure gradient along the vertical direction, and where is the stress, which on is related to the velocity by the von Mises criterion

 ⎧⎪ ⎪⎨⎪ ⎪⎩^τ=(μf+τY|∇^ω|)∇^ωif ^τ⩾τY,∇^ω=0if ^τ⩽τY. (2)

These equations state that as long as a certain stress is not reached, there is no response of the fluid (see Figure 1).

We are considering the exchange flow problem, meaning that we require that the total flux across the horizontal slice is zero, that is

 ∫^Ω^ω=0. (3)

Mathematically, can be seen as a Lagrange multiplier for (3). Moreover, we consider the following boundary conditions:

1. We assume that on the boundaries of , we have a no-slip boundary condition

 ^ω=0on ∂^Ω. (4)
2. For a steady fall motion, the forces exterted by gravity and by the fluid on the solid should be in equilibrium [25]

 ∫∂^Ωs^τ⋅n=ρsg|^Ωs| (5)

where denotes the outward-pointing normal to . Using (5) and the divergence theorem in the second equation of (1), we get that the total pressure gradient on amounts to a buoyancy force

 ∫^Ωspz=(ρs−ρf)g|^Ωs|. (6)

### 2.1 Eigenvalue problems

Following [9, 23], we study minimizers of the functional

 ^G⋄(^ω):=⎧⎪⎨⎪⎩μf2∫^Ωf|∇^ω|2+τY∫^Ωf|∇^ω|−(ρs−ρf)g∫^Ωs^ωif ^ω∈^H⋄+∞else, (7)

over

 ^H⋄={v∈H10(^Ω) ∣∣∣ ∫^Ωv=0, ∇v=0 in ^Ωs}. (8)

Writing the Euler-Lagrange equations for we obtain (1) and (2). Notice that since we work in , the no-slip boundary condition (4) is automatically satisfied, and adequate testing directions are constant on connected components of , which leads to the force balance condition (5).

We proceed to simplify the dimensions in the above functional, so that we can work with just one parameter. Assuming a given length scale , we define the buoyancy number and a velocity scale by

 Y:=τY(ρs−ρf)g^L,^ω0:=(ρs−ρf)g^L2μf, (9)

so that defining the rescaled velocity and corresponding domains by

 ω(x):=^ω(^Lx)^ω0,Ω:=^Ω^L,Ωf:=^Ωf^L, % and Ωs:=^Ωs^L, (10)

we end up with the functional

 G⋄Y(ω):=⎧⎪⎨⎪⎩12∫Ωf|∇ω|2+Y∫Ωf|∇ω|−∫Ωsωif ω∈H⋄+∞else, (11)

over

 H⋄={v∈H10(Ω) ∣∣∣ ∫Ωv=0, ∇v=0 in Ωs}. (12)

By the direct method it is easy to prove (see for instance [9]) that has a unique minimizer, which corresponds to the weak solution in physical dimensions through the transformations in (10). Denoting this minimizer by we have that for every ,

 ∫Ωf∇ωY⋅∇(v−ωY)+Y∫Ωf|∇v|−Y∫Ωf|∇ωY|⩾∫Ωs(v−ωY). (13)

As in [9], one can introduce

 Yc:=supω∈H⋄∫Ωsω∫Ω|∇ω| (14)

and test inequality (13) with and to obtain

 ∫Ω|∇ωY|2=∫Ωf|∇ωY|2=∫ΩsωY−Y∫Ωf|∇ωY|. (15)

From this, and using the definition of in (14) it follows that

 ∫Ω|∇ωY|2⩽∫Ωf|∇ωY|[supω∈H⋄∫Ωsω∫Ω|∇ω|−Y]=(Yc−Y)∫Ωf|∇ωY|. (16)

The last inequality implies, thanks to the homogeneous boundary conditions on , that in as soon as

## 3 Relaxed problem and physical meaning

We determine the critical yield stress , defined in (14) and properties of the associated eigenfunction. The optimization problem (14) is equivalent to computing minimizers of the functional

 E(ω):=∫Ω|∇ω|∫Ωsω % over H⋄. (17)

Because might not attain a minimizer in , we consider a relaxed formulation on a subset of functions of bounded variation.

### 3.1 Functions of bounded variations and their properties

We recall the definition of the space of functions of bounded variation and some properties of such functions that we will use below. Proofs and further results can be found in [2], for example.

###### Definition 1.

Let be open. A function is said to be of bounded variation if its distributional gradient is a Radon measure with finite mass, which we denote by . In particular, if , then . Similarly, for a set with finite Lebesgue measure we define its perimeter to be the total variation of its characteristic function , that is, .

###### Theorem 1.

The space of functions of bounded variation on , denoted , is a Banach space when associated with the norm

 ∥v∥BV(A):=∥v∥L1(A)+TV(v).

The space of functions of bounded variation satisfies the following compactness property [2, Theorem 3.44]:

###### Theorem 2 (Compactness and lower semi-continuity in BV).

Let be a sequence of functions such that is bounded. Then there exists for which, possibly upon taking a subsequence, we have

 vnL1−→v.

In addition, for any sequence that converges to some in ,

 TV(w)⩽liminfTV(wn).

We frequently use the coarea and layer cake formulas:

###### Lemma 1.

Let with compact support, then the coarea formula [2, Theorem 3.40]

 TV(u)=∫∞−∞Per(u>t)dt=∫∞−∞Per(u

holds. If is non-negative, then we also have the layer cake formula [19, Theorem 1.13]

 ∫R2u=∫∞0|{u>t}|dt. (19)

An important role in characterizing constrained minimizers of the functional is played by Cheeger sets, which we now define.

###### Definition 2.

A set is called Cheeger set of if it minimizes the ratio among the subsets of .

The following result is well known and has been stated for instance in [18, Proposition 3.5, iii] and [22, Proposition 3.1]:

###### Theorem 3.

For every non-empty measurable set open, there exists at least one Cheeger set, and its characteristic function minimizes the quotient in . Moreover, almost every level set of every minimizer of this quotient is a Cheeger set.

###### Remark 1.

Some sets may have more than one Cheeger set, which introduces nonuniqueness in the minimizers of the quotient . One example is the set of Figure 6 below.

### 3.2 Generalized minimizers of E

Using the compactness Theorem 2, it follows that the relaxed quotient

 E(ω):=TV(ω)∫Ωsω

of (17) attains a minimizer in the space

 B:={v∈BV(R2) ∣∣∣∫Ωv=0, ∇v=0 on Ωs, v=0 on R2∖Ω}.

Note that the quotient is invariant with respect to scalar multiplication, and we can therefore add the constraint

 \fintΩsv:=1|Ωs|∫Ωsv=1 (20)

to without changing the minimal value of the functional . Thus, the problem of minimizing over is equivalent to the following problem:

###### Problem 1.

Find a minimizer of over the set

 BV⋄:={v∈BV(R2) ∣∣∣∫Ωv=0, \fintΩsv=1, ∇v=0 on Ωs, v=0 on R2∖Ω}.

By using standard compactness and lower semicontinuity results in , it is easy to see [9] that there is at least one solution to Problem 1. In particular, we emphasize that all the constraints above are closed with respect to the topology.

###### Remark 2.

Notice that is larger than the optimization space (42) used in [9] , where it has been assumed that in . See also Section 5.3.

### 3.3 The critical yield limit

We investigate the limit of (the minimizer of , defined in (11)) when . For this purpose we first prove

###### Proposition 1.

The quantity is nonincreasing with respect to . In particular, it is bounded.

###### Proof.

Let . Then, from the definition (11) of being a minimizer of it follows that

 G⋄Y2(ωY2)⩽G⋄Y2(ωY1)=G⋄Y1(ωY1)+(Y2−Y1)∫Ωf|∇ωY1|,
 G⋄Y1(ωY1)⩽G⋄Y1(ωY2)=G⋄Y2(ωY2)+(Y1−Y2)∫Ωf|∇ωY2|,

and summing, we get

 (Y1−Y2)(∫Ωf|∇ωY2|−∫Ωf|∇ωY1|)⩾0,

which implies the assertion. ∎

We are now ready to investigate the convergence of and its rate.

###### Theorem 4.

For and any , we have

 ∫Ω|∇ωY|2=o((Yc−Y)2−ϵ). (21)

Moreover, the sequence of rescaled profiles

 vY:=ωY∫Ω|∇ωY| (22)

converges in the sense of Theorem 2, up to possibly taking a sequence, to a solution of Problem 1.

###### Proof.

Using (15), the definition of (14) and since vanishes on , and is bounded by Proposition 1 we get

 0⩽∫ΩsωY−Y∫Ω|∇ωY|=∫Ω|∇ωY|2⩽(Yc−Y)∫Ω|∇ωY|⩽C(Yc−Y), (23)

where here and in the following we denote by a generic constant, which may be different in each appearance.

That is, we get in particular

 ∫Ω|∇ωY|2⩽C(Yc−Y). (24)

The rate in this estimate can be improved. To do so, we use Poincare’s inequality, which combined with (24) implies that

 ∫ΩsωY⩽C∥ωY∥L2(Ω)⩽C√∫Ω|∇ωY|2⩽C(√Yc−Y). (25)

Now, using (25) and (24) in (23) we get

 0⩽C(√Yc−Y)−Y∫Ω|∇ωY|⩽C(Yc−Y), (26)

which can only hold as if

 ∫Ω|∇ωY|=O(√Yc−Y). (27)

But then, using (27) in (23) provides us with Iterating the above substitutions using the improved estimates for instead of (24) (that is, using a bootstrapping argument) we obtain (21) for any .

Now, the associated functions , defined in (22), have total variation and zero mean. From Theorem 2 it follows that converges in to some . Now, it follows directly from (23) that

 limY→YcY∫Ω|∇ωY|∫ΩsωY=1, (28)

and therefore, using the convergence of , its definition (22) and that , (28) implies

 ∫Ωsvc=limY→Yc∫ΩsvY=limY→Yc∫ΩsωY∫Ωf|∇ωy|=Yc.

Recalling that , the semi-continuity of the total variation with respect to convergence implies , which yields

 Yc∫Ω|∇vc|−∫Ωsvc⩽0,

which can be rewritten as

 Yc⩽∫Ωsvc∫Ω|∇vc|

so is a maximizer of . ∎

From the above result, we see that a minimizer of the quotient can be obtained as a limit of rescaled physical velocities, and therefore carries information about their geometry. For this reason, we will focus on these minimizers in the following.

## 4 Piecewise constant minimizers

We prove the existence of solutions of Problem 1 with particular properties. In our previous work [9] this problem was considered under the assumption that the velocity is constant in the whole . In the situation considered here, the physical velocity is constant only on every connected component of , and the velocity of each solid particle is an unknown. Therefore, the candidates of limiting profiles over which we optimize (belonging to ) also satisfy on .

### 4.1 A minimizer with three values

###### Theorem 5.

There is a solution of Problem 1 that attains only two non-zero values.

The same result has been proved in [9] in the simpler situation when the velocities were considered uniformly constant on the whole . For the proof of Theorem 5, we proceed in two steps:

1. We prove the existence of a minimizer for Problem 1 with attains only finitely many values. This is accomplished by convexity arguments reminiscent of slicing by the coarea (18) and layer cake (19) formulas, but more involved.

2. When considered over functions with finitely many values, the minimization of the total variation with integral constraints is a simple finite-dimensional optimization problem, and standard linear programming arguments provide the result.

#### Step 1. A minimizer with finite range.

To begin the proof, we assume that we are given a minimizer of the total variation in , that is, a solution of Problem 1. We represent by its connected components , ,

 Ωs=N⋃i=1Ωis. (29)

Since belongs to , is constant on every , and we introduce the constants such that

 u∣∣∣∣Ωis=γi. (30)

Note that the constraint (20) reads

 1∑ni=1|Ωis|N∑i=1γi|Ωis|=1. (31)

Defining

 ui:=u⋅1{γi

we have

 u=N∑i=1(ui−γi).

Notice that each minimizes the total variation among functions with fixed integral , and satisfying the boundary conditions on and on .

As a result, the function minimizes the total variation with constraints , and prescribed integral. Lemma 2 below (applied with and ) shows that can be replaced by a five level-set function which has total variation smaller or equal to . Hence can be replaced by the five level-set function without increasing the total variation.

Therefore, the finitely-valued function

 ~u:=N∑i=1(~ui−γi)

is again a solution of Problem 1 (the functions and coincide on , so the constraint is satisfied).

We complete Step 1 by proving Lemma 2, which again requires the proof of an auxiliary lemma.

###### Lemma 2.

Let be two bounded measurable sets, . Then, there exists a minimizer of on the set

 (32)

where the range consists of at most five values, one of them being zero.

###### Proof.

Let be an arbitrary minimizer of in . Splitting at and we can write

 w=(w1+−1)+w(0,1)−w− (33)

with , , and the usual negative part. We see from the coarea formula that

 TV(w) =∫s⩽0Per(w⩾s)+∫0

With this splitting, can be seen to be a minimizer of over

 A−:={v∈BV(R2) ∣∣∣ v=0 % on {w>0}∪R2∖Ω0, ∫Ωv=∫Ωw−}.

By Theorem 3, almost every level set of is a Cheeger set of , the complement of . In particular, if we replace by , where is one such Cheeger set, the total variation doesn’t increase. Therefore, there exists a minimizer of on that reaches only one non-zero value.

With an analogous argumentation we see that, because minimizes on the set

 A1+:={v∈BV(R2) ∣∣∣ v=1 % on {w<1}, ∫Ωv=∫Ωw1+},

there exists a minimizer that writes

 ~w1+=1+ζ1C1

where is a Cheeger set of and is a constant.

Moreover, defining

 μ:=∫Ωw(0,1),

minimizes on the set

 A(0,1)μ:={v∈BV(R2) ∣∣∣ v=1 on {w⩾1},v=0 on {w⩽0} and ∫Ωv=μ}.

We now show that there exists a minimizer of in that attains only three values. Since is one of them, there exists some minimizer of in with values in . We denote by a generic one. In what follows, we denote by the level-sets of .

We prove in Lemma 3 that for almost every , minimizes in . That implies in particular that for a.e. , minimizes perimeter with fixed mass. We introduce the set of points of density for and the set of points of density 0 for , that is

 E(1)s:={x∈Ω∣∣∣limr→0|Es∩Br(x)||Br(x)|=1}andE(0)s:={x∈Ω∣∣∣limr→0|Es∩Br(x)||Br(x)|=0}.

Lebesgue differentiation theorem implies that and a.e.

Now, since the level-sets are nested, the function is nonincreasing. Therefore, there exists such that

 for s>sμ,|Es|⩽μ, and for s

Let us now define

 E+:=⋃s>sμE(1)sandE−:=⋂s

We then have

###### Claim.

If is not empty, minimizes total variation in , with

###### Proof of claim.

By Lemma 3, minimizes total variation in for almost every . Then, let us select a decreasing sequence such that for each , minimizes total variation in . Since in , one has and the semicontinuity for the perimeter gives

 Per(E+)⩽liminfPer(E(1)sn).

In fact, the sequence is bounded. To see this, we fix a value and since we can write for some

 |E(1)sn|=tn|E^s|+(1−tn)|Es1|.

Therefore, applying Lemma 3 again we obtain

 Per(E(1)sn)⩽TV(tn1E^s+(1−tn)1Es1)⩽Per(E^s)+Per(Es1).

Now, let us assume that there exists with and . By the above, for every we can find such that and

 Per(E+)⩽Per(E(1)sn)+δ.

Now, if is small enough, we can find a ball such that and , so we get

 TV(v⋅1Ω∖Bn) ⩽TV(v)+∥v∥∞Per(Bn)⩽Per(E+)−ε+∥v∥∞Per(Bn) (34) ⩽Per(E(1)sn)+∥v∥∞Per(Bn)+δ−ε⩽Per(E(1)sn)−ε2,

and therefore we get a contradiction with the -minimality of .

Selecting an increasing sequence and such that minimizes in , we obtain similarly that minimizes in

To finish the proof of Lemma 2, we distinguish two alternatives. Either or has mass , in which case the claim above implies Lemma 2, or are both nonempty and

 |E+|<μand|E−|>μ.

In the second case, let . Then, and there exists such that The function therefore belongs to . Since is a minimizer of in this set, one must have

 Per(E−)⩽TV(t1Es+(1−t)1E+)⩽|E−|−|E+||Es|−|E+|Per(Es)+|Es|−|E−||Es|−|E+|Per(E+).

This equation rewrites

 Per(Es)⩾|Es|−|E+||E−|−|E+|Per(E−)+|E−|−|Es||E−|−|E+|Per(E+). (35)

Similarly, if , one has and is a convex combination of The same steps lead to the same (35). Finally, one just write (we use (35), the coarea and the layer-cake formulas)

 TV(u) =∫10Per(Es)⩾∫10(|Es|−|E+|)Per(E−)+(|E−|−|Es|)Per(E+)|E−|−|E+| ⩾∫10Per(E−)−Per(E+)|E−|−|E+||Es|+|E−|Per(E+)−|E+|Per(E−)|E−|−|E+| =Per(E−)−Per(E+)|E−|−|E+|μ+|E−|Per(E+)−|E+|Per(E−)|E−|−|E+| =TV(λ1E−+(1−λ)1E+)

with .

As a result, one can replace in the decomposition (33) by a three valued minimizer of in . Therefore, combining the three modified parts we see that there exists a minimizer in

 ~w:=(~w1+−1)+~w(0,1)−~w−

which attains at most five values. ∎

###### Lemma 3.

Let be a minimizer of in with values in , and . Let be a Lebesgue point of and (these two functions are measurable, so almost every is a Lebesgue point for them). Then minimizes