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

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.


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 [31]. They are characterized by two numerical quantities: a yield stress τ Y that must be exceeded for strain to appear, and a fluid viscosity µ f that describes its linear behaviour once it starts to flow (see figure 1). 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 (6)), 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 [27,28], 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 [19,22,18], where non-homogeneous coefficients and different boundary conditions arise. Two-fluid anti-plane shear flows that arise in oilfield cementing are studied in [16,17]. Settling of particles under gravity, not necessarily in anti-plane configurations is also considered in [23,30]. Finally, the previous work [15] 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]. On the numerical aspects, there are several methods available in the literature for the computation of limit loads [7] and Cheeger sets [8,9,6], and both of these problems are closely related to ours, as we shall see below.
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 [15], 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 [12] and an upwind finite difference discretization [11]. 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.

The model
The constitutive law for an incompressible Bingham fluid in three dimensions is given by the von Mises criterion where v is its velocity (for which incompressibility implies div v = 0), and Ev = (∇v + ∇v )/2 is the linearized strain, ∇v ∈ R 3×3 being the Jacobian matrix of the vector v. We denote by σ D the deviatoric part of the Cauchy stress tensor σ(x, y, z) ∈ R 3×3 sym , that is where p is the pressure and tr σ D = 0. These equations state that as long as a certain stress is not reached, there is no response of the fluid (see Figure 1). The geometry we consider consists of a Bingham fluid filling a vertical cylindrical domain Ω × R ⊂ R 3 and a solid inclusionΩ s × R ⊂Ω × R, wherê so thatΩ s is composed of disconnected particles that do not touch the boundary of the domain. We denote byΩ f =Ω \Ω s the portion of the domain occupied by the fluid, and by ρ s , ρ f the corresponding constant densities. We focus on a vertical stationary flow, meaning that the velocity is of the form v =ω(0, 0, 1) T and constant in time. Moreover, all quantities are invariant along the vertical direction, so we can directly consider a scalar velocityω :Ω → R (ω is the velocity of the fluid onΩ f and of the solid inΩ s ), see Figure 2. For the rest of the article, the differential operators denoted by ∇ and div are the two-dimensional ones.
Additionally to incompressibility, we consider the stronger condition of an exchange flow problem, meaning that we require that the total flux across the horizontal slice is zero,Ωω = 0. ( A word on this condition is required. If the cylindrical domain was closed by a bottom fluid reservoir on which no-slip boundaries are assumed, one could use incompressibility, the divergence theorem and the boundary conditions to obtain (3) in any horizontal plane. In our case, while not strictly consistent with an infinite cylinder, it is added as a modelling assumption, reflecting that the region of interest is far away from the bottom of the 3D domain. The same approximation has been used in previous works treating models of drilling and cementing of oil wells [14,16] and justified experimentally in [20] with applications to magma in volcanic conduits.
In the anti-plane case, the Bingham constitutive law (1) can be written in terms of the vector of shear stressesτ = (σ xz , σ yz ) to obtain Since the material occupying the regionΩ s is perfectly rigid, the corresponding constitutive law is ∇ω = 0 onΩ s .
Noting the decomposition of the stress tensor (2), the balance laws for the fluid and the solid particles then write with p z the pressure gradient along the vertical direction. The second equation in (6) expresses that for a steady fall motion, the gravity and buoyancy forces should be in equilibrium with the shear forces exterted by the fluid on each particle [32]. The buoyancy forces b i on each solid particle should be understood as resulting from Archimedes' principle and originating outside the region of interest, being exerted by the bottom reservoir of fluid. This interpretation implies that these forces are proportional to the volume of the solids and the vertical difference of pressure, a fact that we obtain as a consequence of the exchange flow condition in (9). In this equation, n f is the exterior unit normal to ∂Ω f , which at ∂Ω f ∩ ∂Ω s is the interior unit normal to ∂Ω s . These equations are complemented by the following boundary conditions: we assume that on the boundaries ofΩ, we have a no-slip boundary condition and similarly we assume thatω is continuous across the interface ∂Ω s , [ω] ∂Ωs = 0. Figure 2: Anti-plane situation: A falling cylinder, with gravity along its axis of symmetry.

Eigenvalue problems
We assume thatΩ andΩ s are bounded and strongly Lipschitz,Ω s ⊂Ω, that ∂Ω s ∩ ∂Ω = ∅ and thatΩ s has finitely many connected components. Following [15,30], we introduce the functional with the set of admissible velocitieŝ where the argument m is a scalar multiplier for the exchange flow condition (3). Writing the Euler-Lagrange equations in theω argument at an optimal pair for the saddle point problem, we obtain a solution of our constitutive and balance equations (4) and (6), with Notice that since we work inĤ , the no-slip boundary condition (7) and solid constitutive law (5) are automatically satisfied, and adequate testing directions are constant on connected components ofΩ s , which leads to the force balance condition in the second part of (6). Condition (8) is implied (in an appropriate weak form) by the fact thatω ∈ H 1 (Ω). SinceF is convex in its first argument and concave on the second, we can introduce the integral constraint in the space, and focus on the equivalent formulation of finding minimizers ofĜ We proceed to simplify the dimensions in the above functional, so that we can work with just one parameter. Assuming a given length scaleL, we define the buoyancy number Y and a velocity scaleω 0 by so that defining the rescaled velocity ω and corresponding domains by we end up with the functional to be minimized over By the direct method it is easy to prove (see for instance [15]) that G Y has a unique minimizer, which we denote by ω Y and that corresponds to the weak solution of (3), (4), (5), (6), (7), and (8) in physical dimensions through the scaling in (10). Now, noticing that u → Y´Ω f |∇u|−´Ω s u is convex, and that the Gâteaux derivative of u →´Ω f |∇u| 2 at the point ω Y in direction h iś Ω f ∇ω Y · ∇h, differentiating in the direction v − ω Y , as done in [13,Section I.3.5.4] shows that for every v ∈ H , As in [15], one can introduce Y c := sup ω∈H ´Ω s ώ Ω |∇ω| (13) and test inequality (12) with v = 0 and v = 2ω Y to obtain From this, and using the definition of Y c in (13) it follows that The last inequality implies, thanks to the homogeneous boundary conditions on ω, that ω Y = 0 in Ω f as soon as Y Y c .

Relaxed problem and physical meaning
We determine the critical yield stress Y c , defined in (13) and properties of the associated eigenfunction. The optimization problem (13) is equivalent to computing minimizers of the functional Because E might not attain a minimizer in H , we consider a relaxed formulation on a subset of functions of bounded variation.

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.
is said to be of bounded variation if its distributional gradient ∇v is a Radon measure with finite mass, which we denote by TV(v).
In particular, if ∇v ∈ L 1 (A), then TV(v) =´A |∇v|. Similarly, for a set B with finite Lebesgue measure |B| < +∞ we define its perimeter to be the total variation of its characteristic function 1 B , that is, Per(B) = TV(1 B ). In addition, for any sequence (w n ) that converges to some w in L 1 , TV(w) lim inf TV(w n ).
We frequently use the coarea and layer cake formulas: Lemma 1. Let u ∈ BV(R 2 ) with compact support, then the coarea formula [2, Theorem 3.40] (16) holds. If u ∈ L 1 (R 2 ) is non-negative, then we also have the layer cake formula [26, Theorem An important role in characterizing constrained minimizers of the TV functional is played by Cheeger sets, which we now define.

Definition 2.
A set is called Cheeger set of A ⊆ R 2 if it minimizes the ratio Per(·)/| · | among the subsets of A.
The following result is well known and has been stated for instance 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 TV(·)/ · L 1 (A) . One example is the set Ω of Figure 6 below.

Generalized minimizers of E
Using the compactness Theorem 2, it follows that the relaxed quotient Ωs ω of (15) attains a minimizer in the space Note that the quotient E is invariant with respect to scalar multiplication, and we can therefore add the constraint to B without changing the minimal value of the functional E. Thus, the problem of minimizing E over B is equivalent to the following problem: Find a minimizer of TV over the set By using standard compactness and lower semicontinuity results in BV(R 2 ), it is easy to see [15] 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 L 1 topology.
Remark 2. Notice that BV is larger than the optimization space (28) used in [15] , where it has been assumed that v = const. in Ω s . See also Section 5.3.

The critical yield limit
We investigate the limit of ω Y (the minimizer of G Y , defined in (11)) when Y → Y c . For this purpose we first prove and summing, we get which implies the assertion.
We are now ready to investigate the convergence of ω Y and its rate.
Moreover, the sequence of rescaled profiles converges in the sense of Theorem 2, up to possibly taking a sequence, to a solution of Problem 1.
Proof. The first part of the proof is already presented in [13, Section VI 8.3, Equation (8.20)] but we reproduce it here for convenience. As before, let Y c Y 1 > Y 2 0. We use (12) for Y 1 and v = ω Y 2 as well as the same inequality for Y 2 and v = ω Y 1 and sum the inequalities obtained to getˆΩ On the other hand, the Cauchy-Schwarz inequality giveŝ Putting these two inequalities together, we obtain which leads to (19). Now, the associated functions v Y , defined in (20), have total variation 1 and zero mean. From Theorem 2 it follows that v Y converges in L 1 to some v c . Now, it follows directly from (21) and (14) that and therefore, using the L 1 convergence of v Y , its definition (20) and that´Ω |∇ω y | =´Ω f |∇ω y |, Recalling that TV(v Y ) = 1, the semi-continuity of the total variation with respect to L 1 convergence implies TV(v c ) 1, which yields From the above result, we see that a minimizer of the quotient´Ω |∇v| Ωs v 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.

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

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 [15] in the simpler situation when the velocities were considered uniformly constant on the whole Ω s . For the proof of Theorem 5, we proceed in two steps: 1. We prove the existence of a minimizer for Problem 1 which attains only finitely many values. This is accomplished by convexity arguments reminiscent of slicing by the coarea (16) and layer cake (17) 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.
The core of the proof of Theorem 5 is the following lemma, that states that a simplified version of the minimization problem can be solved with finitely many values.

Lemma 2.
Let Ω 1 ⊂ Ω 0 be two bounded measurable sets, ν ∈ R. Then, there exists a minimizer of TV on the set where the range consists of at most five values, one of them being zero.
In turn our proof of Lemma 2 is based on the following minimizing property of level sets, which we believe could be of interest in itself.
The proofs of these two lemmas are located after the proof of Theorem 5.
Proof of Theorem 5.
Step 1. A minimizer with finite range.
To begin the proof, we assume that we are given a minimizer u of the total variation in BV , that is, a solution of Problem 1. We represent Ω s by its connected components Ω i s , i = 1, . . . , N , Since u belongs to BV , u is constant on every Ω i s , and we introduce the constants γ i such that We can assume that γ i γ i+1 . Note that the constraint (18) reads Defining Notice that each u i minimizes the total variation among functions with fixed integral´Ω u i , and satisfying the boundary conditions u = γ i on {u γ i } and u = γ i+1 on {u γ i+1 }.
As a result, the function v i : shows that v i can be replaced by a five level-set functionṽ i which has total variation smaller or equal to TV(v i ). Hence u i can be replaced by the five level-set functioñ u i := γ i +ṽ i (γ i+1 − γ i ) without increasing the total variation. Therefore, the finitely-valued functioñ is again a solution of Problem 1 (the functions u andũ coincide on Ω s , so the constraint ffl Ωsũ = 1 is satisfied).
Step 2. Construction of a three-valued minimizer.
Step 1 provides a solutionũ of Problem 1 that reaches a finite number (denoted as p + 1) of values. We denote its range (listed in increasing order) by Let us now define, for i < 0, where E i ⊂ E j whenever i < j < 0 or i > j > 0. We also have as well as α i < 0 for i < 0 and α i > 0 for i > 0. The constraint on the sign of the α i is made such that the formula (24) holds. Indeed, if the α i change signs, the right hand side of (24) is only an upper bound for TV(ũ).
Introducing the vectors a = (Per(E p − ), · · · , Per(E p + )), b = (|E s p − |, · · · , |E s p + |), c = (|E p − |, · · · , |E p + |), x = (α p − , · · · , α p + ) , minimizing (24) forũ of the form (23) and with the constrained mentioned above is reformulated into finding a minimizer of (a, x) → a T |x| 1 , Denoting by σ ∈ {−1, 1} p ⊆ R p indexed by i ∈ {p − , · · · , p + } with σ i = −1 for i < 0 and σ i = 1 for i > 0, this minimization problem can be rewritten as where σ : x := (x 1 σ 1 , · · · , x p σ p , x p+1 σ p+1 ). The space of constraints is then a (possibly empty) polyhedron given by the intersection of the quadrant σ : x 0 with the two hyperplanes c T x = 0 and b T x = |Ω s |. Now for a point of a polyhedron in R p to be a vertex, we must have that at least p constraints are active at it. Therefore, at least p − 2 of these constraints should be of those defining the quadrant σ : x 0, meaning that at a vertex, at least p − 2 coefficients of x are zero. This polyhedron could be unbounded, but since a 0 and σ : x 0 componentwise, the minimization of a T (σ : x) must have at least one solution in it. Moreover, since it is contained in a quadrant (σ : x 0), it clearly does not contain any line, so it must have at least one vertex ([5, Theorem 2.6]). Since the function to minimize is linear in x, it has a minimum at one such vertex ([5, Theorem 2.7]). That proves the existence of a minimizer of (25) with at least p − 2 of the (α i ) being zero. This corresponds to a minimizer for Problem 1 which has only two level-sets with nonzero values, finishing the proof of Theorem 5.
With this splitting, w − can be seen to be a minimizer of TV over By Theorem 3, almost every level set of w − is a Cheeger set of Ω 0 \{w > 0}, the complement where C 0 is one such Cheeger set, the total variation doesn't increase. Therefore, there exists a minimizerw − of TV on A − that reaches only one non-zero value.
With an analogous argumentation we see that, because w 1+ minimizes TV on the set there exists a minimizerw 1+ that writesw where C 1 is a Cheeger set of {w 1} and ζ 0 is a constant. Moreover, defining µ :=ˆΩ w (0,1) , w (0,1) minimizes TV on the set Lebesgue differentiation theorem implies that E s = E s and E (0) s = Ω \ E s a.e. Now, since the level-sets are nested, the function s → |E s | is nonincreasing. Therefore, there exists s µ such that for s > s µ , |E s | µ, and for s < s µ , |E s | µ.
Let us now define We then have the following fact, to be proved below: Claim. If E ± is not empty, 1 E ± minimizes total variation in A (0,1) To finish the proof of Lemma 2, we distinguish two alternatives. Either E + or E − has mass µ, in which case the claim above implies Lemma 2, or E ± are both nonempty and |E + | < µ and |E − | > µ.
The same steps lead to the same (27). Finally, one just write (we use (27), the coarea and the layer-cake formulas) As a result, one can replace w (0,1) in the decomposition (26) by a three valued minimizer w (0,1) of TV in A sn → E + in L 1 , one has |E + | = lim |E (1) sn | = lim |E sn | and the semicontinuity for the perimeter gives Per(E + ) lim inf Per(E (1) sn ).
In fact, the sequence Per(E sn ) is bounded. To see this, we fix a valueŝ < s µ and since E s 1 ⊂ E (1) sn ⊂ Eŝ we can write for some t n ∈ (0, 1)
Now, let us assume that there exists v ∈ BV(Ω) with´v = |E + | and TV(v) < Per(E + ) − ε. By the above, for every δ > 0 we can find n such that |E + | |E sn | |E + | − δ and Per(E + ) Per(E (1) sn ) + δ. Now, if δ < ε/10 is small enough, we can find a ball B n ⊂ Ω such that´Ω v · 1 Ω\Bn = |E sn | and v ∞ Per(B n ) ε/10, so we get and therefore we get a contradiction with the TV-minimality of E sn . Selecting an increasing sequences n s µ and such that Ω \ E |Es n | , we obtain similarly that 1 E − minimizes TV in A (0,1)

Proof of Lemma 3
Proof of Lemma 3. Since the arguments Ω 0 , Ω 1 are fixed for the course of this proof, we will denote the sets A τ (Ω 0 , Ω 1 ) by A τ for each τ > 0. First, note that for every s 1 < s 2 , the function where (u · 1 u<s 1 − s 1 ) + (v(s 2 − s 1 ) + s 1 ) + (u · 1 u>s 2 − s 2 ) ∈ A ν , which is a contradiction with the minimality of u.
Letting s 0 as in the assumptions, we have just seen that for every h > 0, Finally, let us assume that 1 Es 0 does not minimize total variation in A |Es 0 | . Then, there would exist ε > 0 and u 0 ∈ A |Es 0 | such that Since s 0 is a Lebesgue point, one can find δ > 0 such that for every h δ, Let h δ and B be a ball such that Per(B) ε 4 u 0 ∞ . There exists α such that the function u 0 + α1 B satisfiesˆu Reducing h if needed, one can enforce that |α| 2 u 0 ∞ . Then, which contradicts the minimality of u [s 0 −h,s 0 +h] and proves the claim.

Minimizers with connected level-sets
In this subsection, we refine our analysis slightly, and show the existence of three-valued minimizers for Problem 1 with additional properties. We start with the following definition: This notion is in fact a natural measure-theoretic sense of connectedness for sets for finite perimeter, for more information about it see [1].
Remark 3. By computing the Fenchel dual of Problem 1, it can be seen that the non-zero level-sets of any solution are minimizers of the functional This optimality property in turn implies lower bounds only depending on k for the perimeter and mass of E, and in case it can be decomposed in the sense of Definition 3, the same lower bounds also hold for each set in such a decomposition. In consequence, E can only be decomposed in at most a finite number of sets. The proof of these statements relies heavily on the results of [1], and is presented in [10] for the unconstrained case, and [21] for the case with Dirichlet constraints, as used here.
Assuming these results, one can simplify the level sets of solutions further: Theorem 6. There exists a minimizer for Problem 1 attaining exactly three values for which all non-zero level-sets are indecomposable.
Proof. First, we consider the positive level-set and assume that it is decomposable in two sets Ω 1 , Ω 2 as in Definition 3. Then the corresponding minimizer u can be written as where α, β > 0. Consider a perturbation of u of the form with |h| α, |k| α, and |l| β. Then, since Ω 1 ∩ Ω 2 = ∅, u h ∈ BV if and only if where Ω s i := Ω i ∩ Ω s . These two equations lead to Under our assumptions on h, k, l, Ω 1 and Ω 2 , and since 1 Ω 1 + 1 Ω 2 = 1 Ω 1 ∪Ω 2 , the total variation of the perturbed function u h can be written as Then, because u is a minimizer of TV, it follows that h Per(Ω 1 ) + k Per(Ω 2 ) + l Per(Ω − ) 0.
Since the left hand side and k, l are linear in h, one can replace h by −h and obtain h Per(Ω 1 ) + k Per(Ω 2 ) + l Per(Ω − ) = 0 which shows that u h is also a minimizer. Now since we have one can choose h such that h = −α or k = −α without violating |l| ≤ β, and therefore produce a minimizer whose positive part is either Ω 2 or Ω 1 , respectively. We proceed similarly for the negative part and therefore obtain an indecomposable negative level-set.
Remark 4. In the above proof, through an adequate choice of components for deletion, one can even obtain simply connected level sets. The measure-theoretic notion corresponding to simple connectedness is defined in [1] to be boundedness of the connected components of the complement of the set, these connected components having been defined through indecomposability. For example, assuming that Ω 2 is fully enclosed in Ω − (that is if ∂Ω 2 ∩ ∂Ω − = ∂Ω 2 ), then the variation of u h can also be written which is linear in h as long as k −α−β−l. The equality case in this last constraint corresponds to joining Ω 2 to Ω − , and avoiding creating a "hole" in Ω − by the procedure mentioned above (which replaces α1 Ω 2 by zero). Clearly, this procedure can also be performed for the positive level set, and in fact the "holes" to be deleted could also be connected components of the zero level set. Therefore, a solution in which both the positive and negative level set are simply connected can be obtained.
Remark 5. The intuition behind these last results is that, like in the proof of Theorem 5, the constraints of the problem are linear with respect to the values, and the total variation is also linear as long as the signs of the differences of values at the interfaces do not change. In particular, the points at which the topology of the level sets changes are situations in which these signs change (that is, the values of two adjacent level sets are equal).

Numerical scheme and results
We now turn our attention to the numerical computation of solutions to the eigenvalue for Problem 1. At first, for simplicity, we limit ourselves to the case (considered in [15]) in which the velocities are assumed constant on the whole Ω s . That is, the problem considered is minimization of the total variation in the space where the constraint u ≡ 1 in Ω s corresponds to (18) under this simplification. This restriction corresponds to the case in which either Ω s is connected, so that there is only one solid particle, or all the particles are constrained to move with the same velocity. In Section 5.4 we point out the required modifications for the multi-particle case and present a variety of computed examples.
To compute a minimizer of TV in BV ,1 , we use a standard primal dual algorithm [12]. The constraint´Ω v = 0 is enforced through a scalar Lagrange multiplier q, whereas the conditions v = 0 on ∂Ω and v = 1 on Ω s are encoded as indicator functions. Our discretization of choice is finite differences on a rectangular grid {1, . . . , m} × {1, . . . n}, where in this whole section, for simplicity, we assume that n = m and Ω (0, 1) 2 . This leads to a saddle point problem of the form min v∈X max p∈X 4 q∈R Here, X = R n 2 denotes the space of real-valued discrete functions on the square grid G n = {1, . . . , n}×{1, . . . n}. Since we use Dirichlet boundary conditions, the grid encloses the physical domain. The corresponding constraint set is then where Ω n and Ω n s denote the parts of the grid corresponding to Ω and Ω s respectively (note that to correctly account for perimeter at the boundary we must have Ω n ⊂ {2, . . . , n − 1} × {2, . . . n − 1}). The indicator function (in the convex analysis sense) of a set A is denoted by χ A , so that χ A (x) = 0 if x ∈ A, and +∞ otherwise. ∇ stands for a suitable discrete gradient, whose choice we now discuss.

Discretization
We discretize the problem using the "upwind" scheme of [11] which has the advantage of carrying a high degree of isotropy. The discrete velocity is denoted by v ij , and we use the signed gradient (∇v) ij introduced in [11], containing separate components for forward and backward differences with opposite signs: therefore, at each grid point (i, j) ∈ G n the signed gradient and its corresponding multiplier variables ∇v ij , p ij ∈ (R 2 ) 2 . We note that to compute the gradient when any of the indices is 1 or n one needs to extends the functions outside the grid, but for the problem at hand any choice will do, since Ω n never touches the boundary of the grid.
For us it is important to use a discretization that takes into account derivatives in all coordinate directions equally, since we aim to resolve sharp geometric interfaces that are not induced by a regularization data term. Figure 3 contains a comparison with the results obtained when using forward differences. In that case, the geometry of the interfaces is distorted according to their orientations, a phenomenon which is minimized in the upwind scheme. Using centered differences is also not adequate, since the centered difference operator has a nontrivial kernel and our solutions are constant in large parts of the domain.

Convergence of the discretization
It is well-known that the standard finite difference discretizations of the total variation converge, in the sense of Γ-convergence with respect to the L 1 topology [11], where the discrete functionals are appropriately defined for piecewise constant functions. We now aim to demonstrate that the chosen discretization and penalization scheme still converges and correctly accounts for the boundary conditions in the limit. We introduce, for each (i, j) ∈ {1, . . . , n − 1} 2 , First, we need to decide which constraint to use in the discrete setting. We denote by Our choice is to take Ω n s := such that the discrete constraints are less restrictive than the continuous ones (see Figure 4) and We define TV n as in [11], when the function is piecewise constant on the R n ij and +∞ otherwise.
with ∇v ij ∨ 0 denotes the positive components of ∇v ij , which was defined in (31), therefore picking only the 'upwind' variations. The norm is computed using the inner product in R 2×2 . We first prove the following lemma, which states that the continuous total variation may be computed with multipliers with positive components, mimicking the discrete definition.
We can now prove Gamma-convergence of the discrete problems, implying convergence of the corresponding minimizers.
Proof. First, we study the Γ-liminf and assume that v n → v in L 1 . Notice that we can write TV n (v n ) as a dual formulation where div n p ∈ R 2 is the signed divergence corresponding to (31), and defined by (div n (p)) ij :=(p 1, This is obtained easily by a (discrete) integration by parts in the expression Now, we note that every p : G → (R 2 ) 2 can be viewed as the discretization of some smooth function p : [0, 1] 2 → (R 2 ) 2 , for example stating As a result, one can write TV n (v n ) = sup v n · div n (p) p ∈ C 1 0 [0, 1] 2 , (R 2 ) 2 , |p| 1, p 0 .
Therefore, using Lemma 4 we get For χ C n , let us first assume χ C (v) = +∞, that is either v ≡ 0 on [0, 1] 2 \Ω or v ≡ 1 on Ω s . If the latter holds, then for ε small enough, Ω s ∩ ({v > 1 + 2ε} ∪ {v < 1 − 2ε}) has positive measure and thanks to the L 1 convergence of v n , must have a positive measure for n big enough. That implies χ C n (v n ) = +∞ and the Γ-liminf inequality is trivially true. If χ C (v) < ∞, then χ C (v) = 0 and the inequality is also true since C n ⊂ C. Let now v ∈ BV((0, 1) 2 ). For the Γ-limsup inequality we want to construct a sequence v n → v such that If v / ∈ C, any v n → v gives the inequality.
If v ∈ C, then we first introduce where ψ δ is a convolution kernel with width δ.
that satisfies χ C n (v δ,n ) = 0, and compute Then since v δ ∈ C 1 , it is clear that the right hand side converges to |∂ x v δ |. Note that in the 'upwind' gradient of a smooth function, only one term by direction can be active, then it is also true for v δ,n if n is large enough and therefore TV n (v δ,n ) → TV(v δ ). By a diagonal argument on δ and n, we conclude.

Single particle results
In this section, we again restrict ourselves to the case in which there is either only one particle, or the particles are constrained to move with the same velocity.
In [15], it is shown analytically that the minimizers of TV over the set BV ,1 defined in (28) have level-sets that minimize some geometrical quantities. In particular, Theorem 4.10 shows that there exists a minimizer of the form where Ω − is the maximal Cheeger set of Ω \ Ω s , and Ω 1 is a minimizer of Unfortunately, determining Cheeger sets analytically is only possible in a very narrow range of sets, which makes useful the numerical computation of minimizers. We present two examples of the output of the numerical method for (29) with the constraint (30). First, we consider the "Pacman" shaped Ω s within again a square Ω; see Figure 5 (left). This example induces both asymmetry (left-right) and non-convexity of Ω s which is showed in [15] to influence the geometry of the minimizer. The solution is shown in the central panel of Figure 5 and the right-hand panel shows a histogram of the solution. Figure 5: Numerical results with a single particle, illustrating the results of [15]. On the left of the two subfigures, the free part Ω n \ Ω n s of the computational domain G n is in gray, and the particles Ω n s are white. The minimizers are on the right, where the blue colour represents the negative values and the red colour, the positive ones. More precisely, the left result has nonzero values {−2.38, 7.41} whereas the right result has {−2.35, 7.41}. The corresponding computed critical yield numbers are is Y c = 0.0576 and Y c = 0.0596, with Ω having side length 1.
The second example concerns the geometry depicted in Figure 6 (top panel), in which Ω s denotes the two L-shaped regions in the white dumbbell-shaped domain Ω. By giving a close look, it is clear that there is a Cheeger set of Ω \ Ω s in each half of the domain, which implies the non uniqueness of the minimizer. The question is which solution the computations will converge to. Figure 6 (lower, left and right) show that different minimizers are selected numerically, in this case by using different numerical resolution.

Several particles
We now extend the numerical scheme of to optimize also over the velocities γ i on each component Ω i s . The corresponding problem is again the minimization (29), but with the new constraint set Here, (Ω n s ) i denotes the i-th component of the discrete domain, corresponding to Ω i s . The set C n is the discrete counterpart to the set BV used in sections 3 and 4.
We give several examples that illustrate the behavior of TV-minimizers in BV with a disconnected Ω s . Figure 7 shows the influence of the positions of particles with respect to each other and to the boundary, which might lump up in different configurations. Figure 8 shows two generic situations: 8(a), the flowing part is concentrated around one connected component of Ω s whereas on 8(b), it is concentrated around the whole Ω s .
We also give an example where uniqueness of the minimizer is not expected. In Figure 9, we consider a grid of circular particles in a square. It is easy to see analytically that any subset of the particles can be chosen as positive part of the minimizer. We present two computations at different numerical resolutions that pick two different subsets.
Since the solutions we compute correspond to limit profiles of the original flows (Theorem 7), the results presented both here and in Section 4.1 mean that near the stopping regime Y → Y c the transition between yielded and unyielded regions of the fluid typically happens closer and closer to the particle boundaries and the domain boundaries. This is consistent with the Cheeger set interpretation of the buoyancy case (which was already present in [15]) and the many previous works on non-buoyancy cases ( [28,19], for example).

A random distribution of small particles
We also present two examples of random distribution of square particles in a bigger square. Figure 10 shows the same number of particles distributed in two different ways and the corresponding minimizers. This example shows that the yield number depends strongly on the geometry of the problem, not only on the ratio solid/fluid. An interesting problem would be to investigate the optimal distribution to maximize/minimize this yield number.