Karol L Cascavita

Jeremy Bleyer

Xavier Chateau

Jérémy Bleyer

Alexandre Ern

Hybrid discretization methods with adaptive yield surface detection for Bingham pipe flows *

come

Hybrid discretization methods with adaptive yield surface detection for Bingham pipe flows

Introduction

The Bingham model of viscoplastic fluids, which can be traced back to Bingham [START_REF] Bingham | Fluidity and platicity[END_REF] and Oldroyd [START_REF] Oldroyd | A rational formulation of the equations of plastic flow for a Bingham solid[END_REF], describes materials which behave like a rigid solid for stresses below a critical yield stress and flow like a viscous fluid otherwise. A von Mises criterion based on a threshold on the norm of the deviatoric part of the stress tensor is employed. Such a relatively simple model turns out to be sufficient to reproduce relevant behaviors in many situations encountered in civil engineering and in other applications such as the food, petroleum, and biological industries. For recent reviews on yield stress fluids and their applications, we refer the reader to [START_REF] Balmforth | Yielding to stress: recent developments in viscoplastic fluid mechanics[END_REF][START_REF] Coussot | Bingham's heritage[END_REF]. The mathematical formulation of Bingham flow models rests on variational inequalities, and their analysis has been developed in [START_REF] Mosolov | Variational methods in the theory of the fluidity of a viscous-plastic medium[END_REF][START_REF] Mosolov | On stagnant flow regions of a viscous-plastic medium in pipes[END_REF][START_REF] Mosolov | On qualitative singularities of the flow of a viscoplastic medium in pipes[END_REF] for pipe flows and more generally in [START_REF] Duvaut | Inequalities in mechanics and physics[END_REF]; see also the book [START_REF] Fuchs | Variational methods for problems from plasticity theory and for generalized Newtonian fluids[END_REF] and the references therein for the mathematical analysis and the book [START_REF] Glowinski | Numerical analysis of variational inequalities[END_REF] for the numerical analysis of variational inequalities. The expected regularity of the solution is still an open question, but it seems unreasonable to hope for a high regularity.

The numerical simulation of viscoplastic fluids is extremely challenging because of the lack of smoothness in the constitutive relation between the stress and the strain rate tensors and because the stress field is not unique in the solid regions. Moreover, the yield surface, which separates the solid region (or unyielded region) from the flowing region (or yielded region), is not known a priori. The numerical methods used for these problems over the past three decades can be classified into two groups. The first approach hinges on introducing a small artificial parameter in the constitutive relation, thereby replacing solid zones by flowing zones with a very high viscosity [START_REF] Bercovier | A finite element method for incompressible non-Newtonian flows[END_REF][START_REF] Papanastasiou | Flows of materials with yield[END_REF]. The advantage is that the regularized equations become differentiable and are suitable for Newtonian fluid solvers. Nonetheless, this benefit comes at the expense of difficulties in accurately capturing the yield surface. The second approach is based on introducing an augmented Lagrangian and using a steepest descent method of Uzawa-type to solve the problem. Augmented Lagrangian methods have been introduced in [START_REF] Hestenes | Multiplier and gradient methods[END_REF][START_REF] Powell | A method for nonlinear constraints in minimization problems[END_REF] for nonlinear constrained minimization problems and have been successfully used in the context of Bingham flow models and nonlinear mechanics in [START_REF] Fortin | Augmented Lagrangian methods[END_REF][START_REF] Glowinski | Augmented Lagrangian and operatorsplitting methods in nonlinear mechanics[END_REF]. The work of [START_REF] Roquet | An adaptive finite element method for Bingham fluid flows around a cylinder[END_REF][START_REF] Saramito | An adaptive finite element method for viscoplastic fluid flows in pipes[END_REF] demonstrated the effectivity of the approach, combined with adaptive finite element techniques, to accurately capture the yield surface in various settings, see also [START_REF] Wang | Tracking of the yield surfaces of laminar viscoplastic flows in noncircular ducts using an adaptive scheme[END_REF] and more recently [START_REF] Zhang | An augmented Lagrangian approach to simulating yield stress fluid flows around a spherical gas bubble[END_REF]. Despite the need for introducing two additional tensor fields (a proxy for the strain rate tensor and the corresponding Lagrange multiplier), Augmented Lagrangian methods have progressively emerged over the last decade as the method of choice to simulate viscoplastic flows. For a recent review, we refer the reader to [START_REF] Saramito | Progress in numerical simulation of yield stress fluid flows[END_REF].

The goal of the present work is to combine an Augmented Lagrangian method with a hybrid discretization method of the velocity field to accurately track the yield surface using adaptive mesh refinement. Hybrid discretization methods are formulated in terms of discrete unknowns attached to mesh faces. Cell-based unknowns are also introduced, but they are eliminated locally by a Schur complement technique known as static condensation. Salient examples of hybrid discretization methods include in the loworder case Hybrid Finite Volumes (HFV) [START_REF] Eymard | Discretization of heterogeneous and anisotropic diffusion problems on general nonconforming meshes. SUSHI: a scheme using stabilization and hybrid interfaces[END_REF] and Mimetic Finite Differences (MFD) [START_REF] Brezzi | Convergence of the mimetic finite difference method for diffusion problems on polyhedral meshes[END_REF] (a unifying viewpoint between these two methods is developed in [START_REF] Droniou | A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods[END_REF]), and in the higher-order case, Hybridizable Discontinuous Galerkin (HDG) methods [START_REF] Cockburn | Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems[END_REF] and Hybrid High-Order (HHO) methods [START_REF] Di Pietro | A hybrid high-order locking-free method for linear elasticity on general meshes[END_REF][START_REF] Di Pietro | An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators[END_REF]. HHO and HDG methods have been bridged in [START_REF] Cockburn | Bridging the hybrid highorder and hybridizable discontinuous Galerkin methods[END_REF], whereas in the lowest-order case, HHO methods are closely related to HFV methods, up to an equiv-alent choice of stabilization (see [START_REF] Di Pietro | An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators[END_REF]Section 2.5] for further insight). In this work, we consider a hybrid low-order discretization method belonging to the broad family of Hybrid-Mixed-Mimetic methods [START_REF] Droniou | A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods[END_REF]. The present method can also be recovered by setting the polynomial degree k = 0 in the HHO method. Specifically, the velocity unknown is approximated by a pair of discrete unknowns, consisting of a collection of unknowns attached to the mesh faces and a collection of unknowns attached to the mesh cells. In addition, we consider piecewise constant fields for the proxy of the strain rate tensor and the associated Lagrange multiplier. We do not consider higher-order approximations owing to the moderate regularity expected for the exact solution. The resulting discrete variational inequalities are solved using the Alternating Direction Method of Multipliers (ADMM), originally introduced in [START_REF] Gabay | A dual algorithm for the solution of nonlinear variational problems via finite element approximations[END_REF][START_REF] Glowinski | Sur l'approximation, par éléments finis d'ordre un, et la résolution, par pénalisation-dualité, d'une classe de problèmes de Dirichlet non linéaires[END_REF]. For a review on the ADMM and other steepestdescent methods for Augmented Lagrangian methods, we refer the reader to [START_REF] Glowinski | Numerical methods for nonlinear variational problems[END_REF]. The convergence rate to be expected with the ADMM is of the order of O(1/n), where n is the iteration number [START_REF] He | On non-ergodic convergence rate of Douglas-Rachford alternating direction method of multipliers[END_REF]. Recent advances leading to faster convergence rates include the accelerated ADMM from [START_REF] Treskatis | An accelerated dual proximal gradient method for applications in viscoplasticity[END_REF] and the ADMM with variable step sizes from [START_REF] Bartels | Alternating direction method of multipliers with variable step sizes[END_REF]. We also mention the recent interior-point methods combined with second-order cone programming [START_REF] Bleyer | Advances in the simulation of viscoplastic fluid flows using interior-point methods[END_REF][START_REF] Bleyer | Efficient numerical computations of yield stress fluid flows using second-order cone programming[END_REF]. In what follows, we will resort to the classical ADMM.

In Bingham pipe flow simulations (i.e., the antiplanar configuration where the velocity field is unidirectional), the velocity field is usually discretized using conforming Lagrange finite elements [START_REF] Saramito | An adaptive finite element method for viscoplastic fluid flows in pipes[END_REF], and in the more general case of vector velocity fields with an incompressibility constraint, one uses either infsup stable mixed finite elements [START_REF] Roquet | An adaptive finite element method for Bingham fluid flows around a cylinder[END_REF] or stabilized equal-order finite element pairs [START_REF] Latché | Analysis of the Brezzi-Pitkäranta stabilized Galerkin scheme for creeping flows of Bingham fluids[END_REF]. The advantage of considering hybrid discretization methods, as proposed in the present work, is twofold. First, hybrid discretization methods naturally lead to local conservation properties at the cell level, whereas this conservation property is somewhat less local using nodal-based finite elements (typically, one needs to consider cell patches around mesh vertices, see, e.g., [START_REF] Ern | Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations[END_REF]). Second, and more importantly in the present context, hybrid discretization methods support polygonal/polyhedral meshes including cells with hanging nodes. The use of such meshes can be exploited when performing local mesh adaptation, including local refinement near the yield surface and possibly agglomeration-based coarsening in unyielded regions. One of the main goals of this work is to highlight this possibility, offered by hybrid discretization methods. For earlier work on agglomeration-based coarsening, we refer, e.g., to [START_REF] Bassi | On the flexibility of agglomeration based physical space discontinuous Galerkin discretizations[END_REF][START_REF] Bassi | Agglomeration based discontinuous Galerkin discretization of the Euler and Navier-Stokes equations[END_REF] and references therein. For simplicity, we focus on Bingham pipe flows, but we mention that incompressible flows can be approximated by hybrid discretization methods, as shown, e.g., in [START_REF] Cockburn | Devising HDG methods for Stokes flow: An overview[END_REF][START_REF] Nguyen | A hybridizable discontinuous Galerkin methods for Stokes flow[END_REF] for HDG methods and in [START_REF] Di Pietro | A discontinuous skeletal method for the viscosity-dependent Stokes problem[END_REF] for HHO methods.

This paper is organized as follows. In Section 2, we present the model problem and the Augmented Lagrangian setting. Section 3 is devoted to the hybrid discretization method and the discrete iterative algorithm. In Section 4, we outline the adaptive algorithm for the detection of the yield surface, and we present numerical results on pipes with circular cross-section and with an eccentric annular cross-section [START_REF] Szabo | Flow of viscoplastic fluids in eccentric annular geometries[END_REF][START_REF] Wachs | Numerical simulaiton of steady Bingham flow through an eccentric annular cross-section by distributed Lagrange multiplier/fictitious domain and augmented Lagrangian methods[END_REF]. Finally, conclusions are drawn in Section 5.

Continuous setting

In this section, we briefly outline the Bingham flow model for a viscoplastic fluid and its simplified version for pipe flows. We also recall the Augmented Lagrangian setting.

Let Ω ⊂ R d , d ≥ 2, denote an open, bounded, connected subset of R d with a Lipschitz boundary. Given an external force field f : Ω → R d and considering, for simplicity, homogeneous Dirichlet boundary conditions, the Bingham flow model consists of looking for the total symmetric stress tensor field σ tot : Ω → R d×d s and the velocity field u : Ω → R d such that

∇•σ tot + f = 0 in Ω, ∇•u = 0 in Ω, u = 0 on ∂Ω. (1)
The constitutive equation for a viscoplastic fluid involves a von Mises criterion, i.e., a threshold on the norm of the deviatoric part of the stress tensor σ = σ tot -1 d tr(σ tot )I d with I d the identity tensor of R d×d . Specifically, σ is related to the symmetric velocity gradient ∇ s u = 1 2 (∇u + ∇u T ) as follows:

   σ = 2µ∇ s u + √ 2σ 0 ∇ s u |∇ s u| 2 for |σ| 2 > √ 2σ 0 , ∇ s u = 0 for |σ| 2 ≤ √ 2σ 0 , (2)
with µ > 0 and σ 0 ≥ 0 denoting, respectively, the shear viscosity and the yield stress, and |τ | 2 = √ τ : τ denoting the Frobenius norm for any tensor τ ∈ R d×d . The factor √ 2 is only a convenient convention for simple shear flows in pipes to avoid counting twice the partial derivatives from the symmetric velocity gradient. The region where |σ| 2 > √ 2σ 0 is called the yielded zone and corresponds to liquid behavior, whereas the region where |σ| 2 ≤ √ 2σ 0 is called the unyielded zone and corresponds to solid behavior. The yield surface where |σ| 2 = √ 2σ 0 separates the two regions. The yield surface is not known a priori, and its prediction is an important aspect of viscoplastic flow simulation.

In the specific situation of pipe flows, the above model problem can be simplified by assuming that the velocity field is a divariate uni-directional field, that is, Ω ⊂ R 2 now denotes the cross-section of the pipe, and letting x = (x, y) be the position vector in Ω, we have

u(x) = (0, 0, u(x)) T , 2∇ s u =   0 0 ∂ x u 0 0 ∂ y u ∂ x u ∂ y u 0   , (3)
so that the constitutive equation (2) becomes

   σ = µ∇u + σ 0 ∇u |∇u| 2 for |σ| 2 > σ 0 , ∇u = 0 for |σ| 2 ≤ σ 0 , (4)
where σ is now R 2 -valued and |τ | 2 now denotes the Euclidean norm of any R 2 -valued field τ . Besides, the governing equations (1) become

∇•σ + f = 0 in Ω, u = 0 on ∂Ω, (5)
where the external force field is now scalar-valued. Let L 2 (Ω) denote the space of square-integrable functions over Ω and let us set L 2 (Ω) := L 2 (Ω; R 2 ). The inner product in both spaces is denoted (•, •) Ω and the associated norm • Ω . Let H 1 0 (Ω) denote the space of functions from the Sobolev space H 1 (Ω) (that is, functions from L 2 (Ω) whose weak gradient is in L 2 (Ω)) with null trace on ∂Ω. Assuming f ∈ L 2 (Ω), it is well-known that the velocity field solving (4)-( 5) is the unique minimizer in H 1 0 (Ω) of the following energy functional:

u = arg min v∈H 1 0 (Ω) (H(∇v), 1) Ω -(f, v) Ω , (6)
with the dissipation energy density H : R 2 → R such that

H(g) = µ 2 |g| 2 2 + σ 0 |g| 2 . (7)
The first term in the right-hand side of (7) evaluated with g = ∇v represents the viscous energy and the second term the plastic energy. The weak form of (5) together with ( 4) is (σ, ∇v

) Ω = (f, v) Ω for all v ∈ H 1 0 (Ω) where the stress vector σ is in the sub-differential ∂H(∇u) = {τ ∈ R 2 | H(j) - H(∇u) ≥ τ •(j -∇u), ∀j ∈ R 2 }.
The minimization problem ( 6) can be tackled by decomposition-coordination methods. To decouple the nonlinearity from the velocity field, one introduces an auxiliary field γ ∈ L 2 (Ω), and one enforces that γ = ∇u by means of the Lagrange multiplier σ ∈ L 2 (Ω) (which turns out to be indeed the stress vector) together with a least-squares penalty term. Let us set

X(Ω) := H 1 0 (Ω) × L 2 (Ω) × L 2 (Ω). ( 8
)
The augmented Lagrangian is defined as L :

X(Ω) → R such that L(u, γ, σ) := (H(γ), 1) Ω + (σ, ∇u -γ) Ω + α 2 ∇u -γ 2 Ω -(f, u) Ω , (9)
where α > 0 is the augmentation parameter. The triple (u, γ, σ) ∈ X(Ω) is a saddle-point of the Lagrangian L, that is to say, L(u, γ, τ ) ≤ L(u, γ, σ) ≤ L(v, δ, σ) for all (v, δ, τ ) ∈ X(Ω), if and only if γ = ∇u and the pair (σ, u) solves ( 4)-( 5).

The Augmented Lagrangian method is an iterative method that, at each step, minimizes first the augmented Lagrangian L jointly with respect to the pair (u, γ) and then updates the Lagrange multiplier σ. The idea in the Alternating Direction Method of Multipliers (ADMM) [START_REF] Gabay | A dual algorithm for the solution of nonlinear variational problems via finite element approximations[END_REF][START_REF] Glowinski | Sur l'approximation, par éléments finis d'ordre un, et la résolution, par pénalisation-dualité, d'une classe de problèmes de Dirichlet non linéaires[END_REF] is that the joint minimization is replaced by a successive minimization with respect to γ and u. Let the superscript n ≥ 0 denote the iteration index in the ADMM. Given the initial values u 0 ∈ H 1 0 (Ω) and σ 0 ∈ L 2 (Ω), we perform, for all n ≥ 0, the following three steps:

• Step 1 : Given u n ∈ H 1 0 (Ω) and σ n ∈ L 2 (Ω)
, we define the field γ n+1 ∈ L 2 (Ω) pointwise (a.e.) in Ω as follows:

γ n+1 :=    1 (α + µ) (|θ n | 2 -σ 0 ) θ n |θ n | 2 if |θ n | 2 > σ 0 , 0 if |θ n | 2 ≤ σ 0 , (10)
where θ n := σ n + α∇u n . A straightforward computation shows that we are actually enforcing ∂ γ L(u n , γ n+1 , σ n ) = 0.

• Step 2 : We seek u n+1 ∈ H 1 0 (Ω) solving the following variational problem where stress terms are treated as explicit force terms: For all v ∈ H 1 0 (Ω),

α(∇u n+1 , ∇v) Ω = (f, v) Ω -(σ n -αγ n+1 , ∇v) Ω , (11)
which amounts to enforcing ∂ u L(u n+1 , γ n+1 , σ n ) = 0.

• Step 3 : Finally, we update the Lagrange multiplier

σ n+1 ∈ L 2 (Ω) by setting σ n+1 = σ n + α(∇u n+1 -γ n+1 ). ( 12
)
As the ADMM approaches convergence, we have (∇u

n+1 -γ n+1 ) → 0 so that ∂ σ L(u n+1 , γ n+1 , σ n+1 ) → 0.
The stopping criterion for the ADMM is R n ≤ with > 0 a prescribed stopping tolerance, and where the residual is defined as follows:

R n := σ n+1 -σ n 2 Ω + α 2 ∇u n+1 -∇u n 2 Ω 1 2 . ( 13
)
Remark 2.1 (Variants for decomposition-coordination). It is possible to consider other decompositions for the nonlinear convex minimization problem (6), e.g., one where the auxiliary variable is only used for the plastic energy. It is observed in [START_REF] Dean | On the numerical simulation of Bingham visco-plastic flow: Old and new results[END_REF] that both formulations lead to fairly comparable computational costs. Therefore, we only consider in what follows the method where the auxiliary variable is used for the whole dissipation energy.

Hybrid discretization

In this section, we describe the hybrid discretization method for our model problem, together with the resulting discrete augmented Lagrangian and ADMM.

Discrete minimization problem on the velocity

We employ a hybrid low-order discretization method [START_REF] Di Pietro | An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators[END_REF][START_REF] Droniou | A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods[END_REF][START_REF] Eymard | Discretization of heterogeneous and anisotropic diffusion problems on general nonconforming meshes. SUSHI: a scheme using stabilization and hybrid interfaces[END_REF], with scalar velocity u approximated by a pair of discrete unknowns, consisting of a collection of scalar unknowns attached to the mesh faces and a collection of scalar unknowns attached to the mesh cells. The auxiliary variable γ and the Lagrange multiplier (or stress) σ are both discretized by a collection of vector unknowns attached to the mesh cells. One key advantage of considering a hybrid low-order discretization method is that all the discrete vector fields are piecewise constant on the mesh, which substantially facilitates the enforcement at the discrete level of Step 1 above. Moreover, there is not much expected benefit of using a high-order method, at least near the yield surface, owing to the lack of smoothness of the velocity. We consider a mesh sequence (T ) ∈N where, for all ∈ N, T is a finite collection of nonempty, disjoint, open cells whose closures either cover Ω exactly or cover a close approximation thereof if Ω has curved boundaries. The subscript typically refers to a step in an adaptive mesh procedure which can involve local mesh refinement (and possibly coarsening). A generic mesh cell is denoted T , its boundary ∂T , its diameter h T , and its unit outward normal n T . The mesh cells can have a polygonal shape, and the meshes can possess hanging nodes. We set h := min T ∈T h T (we consider the minimum value instead of the more usual maximum value since the meshes are locally refined). The mesh sequence is assumed to be shaperegular in the sense, e.g., of [START_REF] Di Pietro | A hybrid high-order locking-free method for linear elasticity on general meshes[END_REF], meaning that every mesh T admits a simplicial matching refinement belonging to a shape-regular mesh sequence in the sense of Ciarlet. This assumption can be relaxed to allow for polygonal meshes with some face degeneration as in [START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF]; this is useful for instance when considering agglomeration-based mesh coarsening procedures. We say that the one-dimensional subset F ⊂ Ω is a mesh interface (resp., mesh boundary face) if there are two distinct mesh cells T F,1 , T F,2 ∈ T and an affine hyperplane H F so that F = ∂T F,1 ∩ ∂T F,2 ∩ H F (resp., there is a mesh cell T F ∈ T h and an affine hyperplane H F so that F = ∂T F ∩ ∂Ω ∩ H F , where Ω := int( T ∈T T )). Mesh interfaces are collected in the set F i , mesh boundary faces in F b , and we let F := F i ∪ F b be the collection of the mesh faces. The diameter of a mesh face F ∈ F is denoted by h F . For all T ∈ T , we let F ∂T := {F ∈ F | F ⊂ ∂T } be the collection of the mesh faces composing ∂T . For all F ∈ F ∂T , n T,F is defined as the unit normal to F pointing outward T ; notice that, by definition, n T,F is a constant vector over F . Let A ⊂ Ω be a mesh cell, its boundary, or one of its faces; the inner product in L 2 (A) and in L 2 (A) := L 2 (A; R 2 ) is denoted (•, •) A and the associated norm • A . For a finite set S, |S| denotes its cardinal number.

Let T ∈ T be a mesh cell. Recall from [START_REF] Di Pietro | A hybrid high-order locking-free method for linear elasticity on general meshes[END_REF][START_REF] Di Pietro | An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators[END_REF] that HHO methods are devised from a local reconstruction operator and a local stabilization operator. In the present lowest-order setting where k = 0, both operators can be written explicitly, see also [START_REF] Droniou | A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods[END_REF][START_REF] Eymard | Discretization of heterogeneous and anisotropic diffusion problems on general nonconforming meshes. SUSHI: a scheme using stabilization and hybrid interfaces[END_REF]. The local space of discrete velocity unknowns is defined as

U T := R × R |F ∂T | ∀T ∈ T , (14)
where the first component of the pair is related to the cell and the second component to the faces composing its boundary. We use the notation

(v T , v ∂T ) = (v T , (v F ) F ∈F ∂T ) for a generic element in U T . Let P 0 (T ; R d ) be composed of constant R d -valued functions in T . The local gradient recon- struction operator G T : U T → P 0 (T ; R d ) is such that, for all (v T , v ∂T ) ∈ U T , G T (v T , v ∂T ) = F ∈F ∂T |F | d-1 |T | d (v F -v T )n T,F , (15)
where |F | d-1 is the measure of the face F and |T | d that of the cell T . Let P 0 (F ∂T ) denote the space composed of piecewise constant functions over ∂T . The stabilization operator S T : U T → P 0 (F ∂T ) is such that, for all (v T , v ∂T ) ∈ U T , its restriction to a face F ∈ F ∂T is defined as

S T (v T , v ∂T )| F := v F -v T + G T (0, v ∂T -v T ) • (x F -x T ), (16)
where x F and x T are, respectively, the barycenter of F and of T . The above local gradient reconstruction and stabilization operators are to be used jointly to discretize the viscous energy. The main motivation is the following stability property (see, e.g., [START_REF] Di Pietro | An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators[END_REF]Lemma 4]): letting η ∂T be the piecewise constant function on ∂T such that η ∂T | F = h -1 F , there is a real number η > 0, independent of , such that, for all (v T , v ∂T ) ∈ U T and all T ∈ T ,

η|(v T , v ∂T )| 2 T ≤ G T (v T , v ∂T ) 2 T + η 1 2 ∂T S T (v T , v ∂T ) 2 ∂T ≤ η -1 |(v T , v ∂T )| 2 T , (17)
with the local H 1 -like seminorm on U T such that

|(v T , v ∂T )| 2 T := F ∈F ∂T h -1 F v F -v T 2 F . (18)
This result shows that the gradient reconstruction operator together with the stabilization operator are able to act as an H 1 -like seminorm on U T .

We can now assemble the convex nonlinear minimization problem solved by the discrete velocity. Let us set

U := R |T | × R |F | , (19)
where the first component of the pair is related to the mesh cells and the second component to the mesh faces. We use the notation (v T , v F ) = ((v T ) T ∈T , (v F ) F ∈F ) for a generic element in U , and for any mesh cell T ∈ T and any (v T , v F ) ∈ U , we let (v T , v ∂T ) ∈ U T be composed of the components of (v T , v F ) attached to the cell T and its faces. We enforce the homogeneous Dirichlet condition on the discrete velocity explicitly by introducing the subspace

U ,0 := {(v T , v F ) ∈ U | v F = 0, ∀F ∈ F b }. (20)
The discrete counterpart of ( 6) is to seek the discrete velocity field (u T , u F ) ∈ U ,0 such that (u T , u F ) = arg min

(v T ,v F )∈U ,0 T ∈T (H(G T (v T , v ∂T )), 1) T + α 2 η 1 2 ∂T S T (v T , v ∂T ) 2 ∂T -(f, v T ) T , (21)
with the dissipation energy density H : R 2 → R defined in [START_REF] Bleyer | Advances in the simulation of viscoplastic fluid flows using interior-point methods[END_REF] and where α is the augmentation parameter of the Lagrangian. The stabilization term is needed to ensure stability in the discrete system of step 2 of the ADMM (see below). Our numerical experiments in section 4.1 show that this term does not have an influence on the velocity error if the parameter α is not too small (e.g., of the order of the shear viscosity or larger).

Discrete augmented Lagrangian and ADMM

The derivation of the discrete augmented Lagrangian is identical to the continuous setting in that we introduce, for all T ∈ T , the auxiliary field γ T ∈ R 2 , and we enforce that γ T = G T (u T , u ∂T ) by means of the Lagrange multiplier σ T ∈ R 2 . Then, setting

X(T ) := U ,0 × R 2|T | × R 2|T | , (22)
the discrete augmented Lagrangian is defined as L : X(T ) → R such that

L ((u T , u F ), γ T , σ T ) := T ∈T (H(γ T ), 1) T + α 2 η 1 2 ∂T S T (v T , v ∂T ) 2 ∂T -(f, u T ) T + (σ T , G T (u T , u ∂T ) -γ T ) T + α 2 G T (u T , u ∂T ) -γ T 2 T , (23)
where we used the notation γ T = (γ T ) T ∈T and

σ T = (σ T ) T ∈T for all γ T , σ T ∈ R 2|T | . If ((u T , u F ), γ T , σ T ) is a saddle-point of the discrete Lagrangian L , then G T (u T , u ∂T ) = γ
T for all T ∈ T and (u T , u F ) is the minimizer of ( 21) with β = α.

The ADMM applied to the discrete augmented Lagrangian reads as follows: Given the initial values (u 0 T , u 0 F ) ∈ U ,0 and σ 0 T ∈ R 2|T | , we perform, for all n ≥ 0, the following three steps:

• Step 1: Given (u n T , u n F ) ∈ U ,0 and σ n T ∈ R 2|T | , we define the field γ n+1 T ∈ R 2|T | such that, for all T ∈ T , γ n+1 T :=    1 (α + µ) (|θ n T | 2 -σ 0 ) θ n T |θ n T | 2 if |θ n T | 2 > σ 0 , 0 if |θ n T | 2 ≤ σ 0 , (24)
where θ n T := σ n T + αG T (u n T , u n ∂T ).

• Step 2: We seek (u n+1 T , u n+1 F ) ∈ U ,0 solving the following discrete variational problem where the stress terms are treated as explicit force terms: For all (v T , v F ) ∈ U ,0 ,

T ∈T α(G T (u n+1 T , u n+1 ∂T ), G T (v T , v ∂T )) T + α(η ∂T S T (u n+1 T , u n+1 ∂T ), S T (v T , v ∂T )) ∂T = T ∈T (f, v T ) T -(σ n T -αγ n+1 T , G T (v T , v ∂T )) T . ( 25
)
As it is customary with hybrid discretization methods, the cell unknowns can be locally eliminated in each mesh cell by using a Schur complement technique (also known as static condensation) whereby they are expressed locally in terms of the face unknowns and the force terms. After this elimination is performed, (25) reduces to a linear system in terms of the face unknowns whose stencil couples neighboring faces sharing a mesh cell. As shown, e.g., in [START_REF] Cockburn | Bridging the hybrid highorder and hybridizable discontinuous Galerkin methods[END_REF], this linear system can be interpreted as a global transmission problem expressing the law of action and reaction for equilibrated face-based tractions. We refer the reader to [START_REF] Cicuttin | Implementation of Discontinuous Skeletal methods on arbitrary-dimensional, polytopal meshes using generic programming[END_REF] for further aspects of the implementation of HHO methods.

• Step 3: Finally, we update the Lagrange multiplier σ n+1 T ∈ R 2|T | by setting, for all T ∈ T ,

σ n+1 T = σ n T + α(G T (u n+1 T , u n+1 ∂T ) -γ n+1 T ). ( 26
)
The stopping criterion for the ADMM is R n ≤ with > 0 a prescribed stopping tolerance, and where the discrete residual is defined as follows:

R n := T ∈T σ n+1 T -σ n T 2 T + α 2 T ∈T G T (u n+1 T , u n+1 ∂T ) -G T (u n T , u n ∂T ) 2 T 1 2 . ( 27
)
Remark 3.1 (Jump-based plastic energy). In each mesh cell T ∈ T , the stabilization operator S T defined in (16) essentially depends on the difference between the face unknowns and the trace of cell unknown, i.e., it represents a jump term. We notice that in [START_REF] Duvaut | Inequalities in mechanics and physics[END_REF] and in [START_REF] Eymard | Discretization of heterogeneous and anisotropic diffusion problems on general nonconforming meshes. SUSHI: a scheme using stabilization and hybrid interfaces[END_REF] this jump term only contributes to the viscous energy. This is rather natural since the jump term tends to zero as the mesh is refined (see, e.g., [18, Eq. ( 45)]). One can consider an additional jump-based contribution to the plastic energy. This presents, however, the drawback of introducing additional auxiliary variables and Lagrange multipliers, thus significantly increasing the computational costs. Moreover, we have observed in our numerical experiments that these additional auxiliary variables and Lagrange multipliers do not improve the accuracy of the simulations, even on relatively coarse meshes.

Mesh adaptation for yield surface detection

The mesh sequence (T ) ∈N is generated iteratively by means of an adaptive procedure. The step ≥ 0 of this procedure consists of solving the discrete problem on T , marking some cells of T for refinement, and producing the new mesh T +1 . Here, we only consider local mesh refinement in view of yield surface detection. An agglomeration-based mesh coarsening procedure in the unyielded region is briefly outlined at the end of Section 4.1 below.

To describe our adaptive procedure, we introduce, for any mesh cell T ∈ T , the subset N v T = {T ∈ T | ∂T ∩ ∂T = ∅} which collects all the mesh cells (including T itself) sharing at least a point with T and the subset

N f T = {T ∈ T | ∂T ∩ ∂T ∈ F h }
which collects all the mesh cells sharing a face with T . Furthermore, during the adaptive procedure, any mesh cell is assigned a level ι(T ) ∈ N which is equal to the number of refinements that were necessary to obtain the cell T from its ancestor in the original mesh T 0 (and all the mesh cells in T 0 are assigned the level 0); thus, if T ∈ T , then ι(T ) ∈ {0, . . . , }.

The marking of the cells in T for refinement involves two steps. First, all the cells in T are flagged according to the value of the quantities {|θ n T | 2 } T ∈T once the ADMM has reached convergence, that is, the cell flag is set to flag(T ) = SOLID if |θ n T | 2 < σ 0 and to flag(T ) = LIQUID otherwise. Then, the first marking step increases by one the level of all the mesh cells such that the set M v T := {flag(T ) | T ∈ N v T } contains both LIQUID and SOLID labels. The aim of the second marking step is to control the level of hanging nodes produced by the local refinement procedure. Although the present hybrid discretization method can support more than one hanging node per edge, we actually limit the number of hanging nodes per edge to one at most so as to ensure a smoother transition of the mesh near the yield surface. To achieve this goal, we perform a second marking of the cells to ensure that in any subset N f T , the difference of cell levels is one at most. The whole marking procedure is summarized in Algorithm 1. Finally, the mesh refinement step consists in performing a subdivision of the marked cells in 4 congruent subcells (known in the literature as red refinement [START_REF] Verfürth | A posteriori error estimation techniques for finite element methods[END_REF]). Mark T and increase level ι(T ) = ι(T ) + 1 Build the subset of neighbors N f T 10:

for T ∈ N f T do 11: if |ι(T ) -ι(T )| > 1 then 12:
Mark and increase level of cell with lower level end for 15: end for Remark 3.2 (Comparison). Other criteria for flagging the cells near the yield surface can be considered; for instance, in [START_REF] Wang | Tracking of the yield surfaces of laminar viscoplastic flows in noncircular ducts using an adaptive scheme[END_REF][START_REF] Zhang | An augmented Lagrangian approach to simulating yield stress fluid flows around a spherical gas bubble[END_REF], the contrast in the norm of the strain rate tensor in vertex-based cell patches is used. We also notice that the present local mesh refinement procedure is more local than those considered in the context of conforming finite elements which need to build a mesh without hanging nodes. Here, there is still some non-locality because we limit the number of hanging nodes to one, but we can also consider a larger value.

Numerical results

In this section, we present our numerical results for Bingham flow problems. We consider flows in pipes with a circular cross-section and an eccentric annular cross-section. The first setting leads to one of the few Bingham pipe flow problems with known exact solution, whereas the second setting is well-documented in the literature (see, among others, [START_REF] Bleyer | Advances in the simulation of viscoplastic fluid flows using interior-point methods[END_REF][START_REF] Szabo | Flow of viscoplastic fluids in eccentric annular geometries[END_REF][START_REF] Wachs | Numerical simulaiton of steady Bingham flow through an eccentric annular cross-section by distributed Lagrange multiplier/fictitious domain and augmented Lagrangian methods[END_REF]). In all cases, the external force f , which represents the transverse pressure gradient forcing the flow, is constant over the cross-section. The (non-dimensional) Bingham number is defined as Bi = σ 0 R µV where R is a reference length (the radius of the outer boundary) and the reference velocity V is computed as

V = f R 2 2µ
, so that we have Bi = 2σ 0 f R . In our numerical experiments, we set the convergence threshold for the discrete residual defined in ( 27) to = 10 -8 . The augmentation parameter is set to α = 10 and the pipe radius to R = 1.

Circular cross-section

The cross-section Ω is a disk centered at the origin and of radius R. The exact velocity only depends on the radial coordinate r and is given by

u(r) V =        1 2 1 - r 2 R 2 -Bi 1 - r R if r R ≥ Bi, 1 2 (1 -Bi) 2 if r R ≤ Bi. (28)
The flow consists of a central solid region of radius Bi × R moving at a uniform velocity and a liquid region connecting the central rigid plug to the pipe walls. In particular, no flow occurs (u = 0) if the Bingham number is larger than 1 whereas the case Bi = 0 (no solid region) corresponds to the classical Newtonian Poiseuille flow. We use the hybrid discretization method from Section 3.1 together with the augmented Lagrangian and ADMM from Section 3.2. We actually discretize only one eighth of the unit disk and enforce symmetry conditions on the sides of the resulting angular sector.

In Figure 1, we display the velocity profile along a line of azimuthal angle π 6 for uniform triangulations with mesh-sizes h ∈ {0.5, 0.1, 0.05}. In each mesh cell, the velocity is reconstructed as an affine function using the cell-value to assign its mean-value and the gradient reconstruction operator to evaluate its gradient. We observe very close profiles between h = 0.1 and h = 0.05. Table 1: Energy difference between circumscribed and inscribed uniform triangulations with radial boundary divided by J = 2 j segments and h J the maximum mesh-size of the circumscribed mesh-triangulation. j h J Bi = 0.3 Bi = 0.9 3 1.55E-1 4.323e-4 9.667e-5 4 7.08E-2 1.076e-4 2.159e-5 5 3.57E-2 2.699e-5 5.753e-6 6 1.81E-2 9.237e-6 2.131e-6

We show in Figure 3 the velocity error u-u h L 2 (T ) vs. the mesh-size for the values of the augmentation parameter α = 10 k with k ∈ {-2, -1, 0, 1, 2}. We first observe for α ≥ 1, the velocity error is essentially independent of this parameter (and converging with first-order in the mesh-size). For lower values of α, we see a degradation of the velocity error, which is more pronounced for coarser meshes. Additionally, we found that the number of iterations to reach a residual threshold of 10 -5 for the above values of α is {5, 4, 14, 93, 913}. The choice α = 10 made below appears as a reasonable compromise between (velocity) accuracy and computational efficiency.

The approximation of the curved boundary with the mesh triangulation introduces an error on the boundary conditions. To address this issue, we compute the energy of the discrete solution using sequences of meshes producing circumscribed and inscribed polygons to the exact circular domain. By placing the mesh nodes on the unit circle, we form an inscribed polygon, whereas by placing them on a circle with radius 1/ cos(α/2J), where J = 2 j is the number of boundary segments and j ∈ {3, 4, 5, 6}, we form a circumscribed polygon. In Table 1, we report the energy difference for the discrete solutions produced by the circumscribed and inscribed triangulations. We observe that this difference converges to zero with essentially second-order in the mesh-size.

We perform local mesh adaptation to detect the yield surface as described in Section 3.3 starting from an initial triangulation T 0 with uniform meshsize h 0 = 0.1. We run the adaptive loop five times, thereby producing the locally refined triangulations T , ∈ {1, . . . , 5}. Figure 2 shows the discrete residual R as a function of the iteration number in the ADMM for a fine uniform triangulation with mesh-size h = 0.0075 and for the adaptation levels ∈ {1, 3, 5} and for Bi = 0.3 and Bi = 0.9. We observe a similar convergence behavior for all meshes and for both Bingham numbers. Notice that we do not initialize the ADMM at some level ≥ 1 with the interpolated solution from level ( -1), but simply by zero; we indeed observed that after a few tens of iterations, the behavior of the ADMM is similar for both choices of the initialization. Moreover, we observe in this specific example that in all cases, the ADMM first reaches a convergence rate of order O(1/n), but then switches to a faster rate as the discrete residual approaches zero.

In Figures 4 and5, we show for Bi = 0.3 and Bi = 0.9, respectively, the locally refined meshes and the stress norm |σ| 2 for the adaptation levels ∈ {0, 1, 3, 5}. The exact yield surface is depicted in red. We observe the progressive capture of the yield surface as the adaptation loop progresses. Note that the locally refined meshes have a mesh-size near the yield surface that is about the same as that of the uniform triangulation (h = 0.0075), but they lead to reduction factors in the number of DOFs of 15 for Bi = 0.3 and of 6 for Bi = 0.9. To provide a more quantitative measure of accuracy in the detection of the yield surface, we report in Tables 2 and3, for Bi = 0.3 and Bi = 0.9, respectively, the error on the mean-value and variance for the distribution of the normalized radial coordinate of the mesh nodes attached to edges separating a cell marked as belonging to the unyielded region with a cell marked as belonging to the yielded region (we can view the collection of such edges as the discrete yield surface). The radial coordinate is normalized by the exact radius of the yield surface which is equal to Bi. In Tables 2 and3, the column DOFs reports the total number of face-based velocity unknowns. We observe from these tables that the error on the mean-value and the variance converge to zero, essentially with first-order in h . It is also interesting to notice that the variance ranges between one-third and one-fourth of h . Finally, we briefly comment on the possibility of agglomeration-based mesh coarsening in the unyielded region. We consider here just a simple illustration where we perform the coarsening on the locally refined mesh T 3 . We agglomerate all the cells in T 3 that are in the unyielded region and such that all their neighbors in the sense of faces are also in the unyielded region. We call the resulting mesh T * 3 . These meshes are shown in Figure 6 for Bi = 0.3 and Bi = 0.9. They lead, respectively, to 511 and 1027 DOFs (face-based velocity unknowns), which represents a relative savings of 9.2% and 58.9% with respect to the corresponding mesh without agglomeration.

Computing a new discrete solution on T * 3 , we observe that the H 1 -error on the velocity is 1.5×10 -3 whereas it is 1.1×10 -3 when computed on the mesh T 3 . Thus, performing an agglomeration-based mesh coarsening does not really hamper accuracy while at the same time reducing the computational costs, especially for Bi = 0.9 where the unyielded region is quite large.

Eccentric annular cross-section

The geometric setting for the eccentric annular cross-section is illustrated in Figure 7. The parameters are R e = 1 for the external radius, R i = 0.4 for the inner radius, and e = -0.15 for the eccentricity of the inner hole. We take advantage of the symmetry by discretizing only one half of the cross-section. The Bingham number is set to Bi = 0.2.

Figure 8 shows the discrete residual R as a function of the iteration number in the ADMM for a fine uniform triangulation with mesh-size h = 0.0075 and for the adaptation levels ∈ {2, 4, 6} with initial mesh-size h = 0.05. As before, we observe a similar convergence behavior for all meshes (with a few irregularities). However, in contrast with the previous test case, the convergence rate is O(1/n) when the discrete residual goes below 10 -6 , whereas this rate is higher when the residual is in the range [10 -6 , 10 -2 ].

The locally refined meshes at the adaptation levels ∈ {2, 6} are presented in Figure 9.

The left panel of Figure 10 shows the spatial distribution of the velocity field and of the unyielded region (notice that this region is composed of two connected components having a different shape). The right panel of Figure 10 shows the spatial distribution of the stress norm |σ| 2 . These results are obtained using a fine uniform triangulation with mesh-size h = 0.0075, and can be compared with the results presented in Figure 11 (velocity field and unyielded region) and in Figure 12 (stress norm) which are obtained using locally refined meshes. We can observe the progressive capture of the yield surface and that the predictions on the sixth locally refined mesh match well with those on the fine uniform triangulation. We notice that the locally refined mesh T 6 leads to about 10 5 DOFs as the uniform triangulation with h = 0.0075, but the local mesh-size in T 6 near the yield surface is about 10 times smaller. This leads to a much sharper resolution as shown in Figure 13, where we provide a zoom of two regions around the yield surface for the stress norm distribution.

The above results on the velocity field and the localization of the yield surface also match well with the predictions reported in [START_REF] Bleyer | Advances in the simulation of viscoplastic fluid flows using interior-point methods[END_REF]. To provide a more quantitative comparison, we present in Table 4 the value of the flux Q = (u h , 1) Ω across the cross-section as predicted using the methodology from [START_REF] Bleyer | Advances in the simulation of viscoplastic fluid flows using interior-point methods[END_REF] (recall that it is based on a conforming finite element discretization, and the numerical solver uses either the ADMM or second-order cone programming) and the present one. We consider uniformly refined meshes as the conforming finite element solver does not support locally refined meshes with hanging nodes. We consider mesh-sizes h ∈ {0.1, 0.05, 0.025}. We observe that the results obtained with both methods match quite closely.

Conclusions and perspectives

We presented a new discretization scheme for Bingham pipe flows by means of a hybrid low-order method combined with the ADMM. We exploited the possibility of using polygonal meshes to perform local mesh adaptation and better capture the yield surface. We considered local mesh refinement around the yield surface. We also performed some numerical tests showing the agreement of our results with analytical solutions and results from the Table 4: Fluxes Q for σ 0 = 0.1 and different values of f computed with the finite element method and the present hybrid low-order method. f h 0.5 1 1.5 2 Finite element method 0.100 3.55e-3 2.38e-2 4.57e-2 6.77e-2 0.050 3.40e-3 2.34e-2 4.49e-2 6.67e-2 0.025 3.36e-3 2.33e-2 4.48e-2 6.64e-2

Hybrid low-order method 0.100 3.55e-3 2.37e-2 4.57e-2 6.77e-2 0.050 3.40e-3 2.34e-2 4.49e-2 6.67e-2 0.025 3.36e-3 2.32e-2 4.47e-2 6.63e-2 literature.

Several tracks can be considered for future research. First, the computational effectiveness of the actual numerical procedure can be enhanced by using hp-refinement techniques, a posteriori error estimation in the flowing (yielded) regions, and some of the acceleration techniques mentioned in the introduction for the ADMM [START_REF] Bartels | Alternating direction method of multipliers with variable step sizes[END_REF][START_REF] Treskatis | An accelerated dual proximal gradient method for applications in viscoplasticity[END_REF] or second-order cone programming [START_REF] Bleyer | Advances in the simulation of viscoplastic fluid flows using interior-point methods[END_REF][START_REF] Bleyer | Efficient numerical computations of yield stress fluid flows using second-order cone programming[END_REF]. Second, future work can include solving Stokes-like viscoplastic flows and studying their interaction with gas bubbles.

Algorithm 1 procedure 1 : 3 : 4 :

 1134 Marking Require: collection of cell flags {flag(T )} T ∈T 2: for T ∈ T do Build the subsets of neighbors N v T and their labels M v T if M v T = {LIQUID,SOLID} then 5:

Figure 1 :

 1 Figure 1: Velocity profiles along the line with azimuthal angle π 6 for uniform triangulations with mesh-sizes h ∈ {0.5, 0.1, 0.05}.

Figure 2 :

 2 Figure 2: Discrete residual vs. ADMM iteration number for Bi = 0.3 (left)and Bi = 0.9 (right) for a fine uniform triangulation with h = 0.0075 and for locally refined meshes at adaptation levels ∈ {1, 3, 5}.

Figure 3 :

 3 Figure 3: Velocity error ||u -u h || L 2 (T ) vs. the mesh-size for Bi = 0.3 and various values of the augmentation parameter α.

Figure 4 :

 4 Figure 4: Stress norm |σ| 2 for Bi = 0.3 and, from top to bottom and from left to right, for adaptation levels ∈ {0, 1, 3, 5}. The colormap (from dark blue to light blue) corresponds to the interval [σ 0 , σ 0 + δ] with σ 0 = 0.15 and δ = 0.005.

Figure 5 :

 5 Figure 5: Stress norm |σ| 2 for Bi = 0.9 and, from top to bottom and from left to right, for adaptation levels ∈ {0, 1, 3, 5}. The colormap (from dark blue to light blue) corresponds to the interval [σ 0 , σ 0 + δ] with σ 0 = 0.45 and δ = 0.005.

Figure 6 :

 6 Figure 6: Locally refined meshes with agglomeration-based coarsening in the unyielded region for Bi = 0.3 (left) and Bi = 0.9 (right).

Figure 7 :

 7 Figure 7: Geometric setting for the eccentric annular cross-section.

1 Figure 8 :

 18 Figure 8: Discrete residual vs. ADMM iteration number for adaptation levels ∈ {0, 2, 4, 6}.

Figure 9 : 20 Figure 10 :

 92010 Figure 9: Locally refined meshes for = 2 (left) and = 6 (right).

Figure 11 :

 11 Figure 11: Velocity field and unyielded region (gray). From top to bottom and left to right: adaptation levels ∈ {0, 2, 4, 6}.

Figure 12 :

 12 Figure 12: Stress norm with colormap (from dark gray to light gray) corresponding to the interval [σ 0 , σ 0 + δ] with σ 0 = 0.1 and δ = 0.005. From top to bottom and left to right: adaptation levels ∈ {0, 2, 4, 6}.

Figure 13 :

 13 Figure 13: Zoom of two regions around the yield surface. Top row: calculation using the uniform triangulation with h = 0.0075; bottom row: calculations using the locally refined mesh T 6 . Both meshes lead to approximately the same number of DOFs.

Table 2 :

 2 Error on the mean-value and variance for the distribution of the radial coordinate at the vertices of the discrete yield surface for Bi = 0.3.

	h	DOFs mean err.	var.
	0 8.10E-2	167	3.55E-2	1.78E-2
	1 5.00E-2	217	2.51E-2	1.67E-2
	2 2.50E-2	343	8.93E-3	7.79E-3
	3 1.25E-2	563	3.39E-3	3.80E-3
	4 6.25E-3 1011	1.26E-3	1.73E-3
	5 3.12E-3 1864	7.20E-4	9.07E-4

Table 3 :

 3 Error on the mean-value and variance for the distribution of the radial coordinate at the vertices of the discrete yield surface for Bi = 0.9.

	h	DOFs mean err.	var.
	0 8.10E-2	167	6.29e-3	1.27E-2
	1 4.62E-2	321	4.03e-4	9.70E-3
	2 2.31E-2	632	2.00e-3	5.71E-3
	3 1.15E-2 1256	1.49e-3	3.04E-3
	4 5.77E-3 2499	7.22e-4	1.99E-3
	5 2.88E-3 5058	4.24e-4	8.66E-4

* Draft version, May 22, 2018. This work has benefited from a French government grant managed by ANR within the framework of the National Program Investments for the Future ANR-11-LABX-0022-01.