Alexander Kolpakov

Jun Murakami

Combinatorial decompositions, Kirillov-Reshetikhin invariants and the Volume Conjecture for hyperbolic polyhedra

Keywords: Combinatorial Decompositions, Kirillov-Reshetikhin Invariants, and the Volume Conjecture for Hyperbolic Polyhedra hyperbolic polyhedron, quantum invariant, I-H move, Whitehead move, volume

come   L'archive ouverte pluridisciplinaire

Introduction

In this paper we suggest a method of computing volume for a polytope P in three-dimensional hyperbolic space H 3 . We require the polytope P to be simple, so that the computation of volumes may combine two processes of initially different nature. First, a reduction of P (as a trivalent graph) to the tetrahedron by a sequence of I -H (or Whitehead) and so-called "capping" moves. And second, a decomposition of P into a number of generalised tetrahedra T i , i = 1, . . . , n, such that T i partition P into geometric parts, and Vol P = n i=1 Vol T i .

1

With such a decomposition we associate a function Φ(ℓ 1 , . . . , ℓ m ), m ≥ 1, that depends on some additional geometric parameters of the decomposition (in fact. lengths of some common perpendiculars to the faces of P ). Under some conditions, we show that Vol P can be expressed through the value of Φ(ℓ * 1 , . . . , ℓ * m ) at a critical point (ℓ * 1 , . . . , ℓ * m ). Finally, we use our method in a large number of examples, by implementing it in Wolfram Mathematica R [?]. We also check our computations, where possible, with the Orb software [?]. This part of our work falls in line with the recent study on the Volume Conjecture for generalised tetrahedra by Costantino and Murakami [?] and for general hyperbolic polytopes by Costantino, Guéritaud and van der Veen [?]. Next, given a polytope P , we consider its 1-skeleton Γ as a trivalent graph, compute the Kirillov-Reshetikhin invariants of Γ with an appropriate colouring of its edges determined by the corresponding dihedral angles of P and study their asymptotic behaviour. We conjecture that the sequence of combinatorial moves that is used in order to compute the Kirillov-Reshetikhin invariants of the corresponding trivalent graph Γ is associated with the sequence of moves used in the combinatorial reduction of P . Then the geometric decomposition of P associated with the respective sequence of I -H and capping moves might be used in order to establish a link between the volume of P and the asymptotic behaviour of the Kirillov-Reshetikhin invariants of Γ. Such a link has been colloquially named the Volume Conjecture by various authors, and is first established in [?]. A similar study of quantum invariants in the case of manifolds with boundary has been recently performed by Chen and Yang [?]. The paper is structured as follows: after recalling some preliminary results on the volume formulas for generalised hyperbolic tetrahedra (cf. [?, ?] and [?]), we proceed to the description of our method and formulate the main statement of the paper. We illustrate our method in two main working examples: computing the volume of a hyperbolic prism Π, and that of a hyperbolic "pleated" prism P. Then, we provide a more computationally complicated example of a dodecahedron with various Coxeter dihedral angles. However, the method of present paper is not restricted to Coxeter polytopes. We expect it to work in a wide variety of dihedral angles, and expect that it can be generalised to allow computing volumes of knotted trivalent graphs. Finally, we state a number of conjectures relating the volume of a hyperbolic polyhedron P to the Kirillov-Reshetikhin invariants of its 1-skeleton Γ, viewed as a trivalent graph in the topological 3-sphere. Various numeric experiments are described, that corroborate our conjectures.

Acknowledgements. A.K. was supported by the Swiss National Science Foundation (SNSF project no. P300P2-151316) and the Japan Society for the Promotion of Science (Invitation Programs for Research project no. S-14021). A.K. is thankful to Waseda University for hospitality during his visit. J.M. was supported by Waseda University (Grant for Special Research Projects no. 2014A-345) and the Japan Society for the Promotion of Science (Grant-in-Aid projects no. 25287014, no. 25610022).

Preliminaries

Below, we describe a method to compute the volume of a simple polyhedron P ⊂ H 3 that admits a decomposition into generalised hyperbolic tetrahedra. We recall, that a convex polytope P is simple if the valence of each vertex equals the dimension of P . Here, it means that each vertex of P is trivalent. Below we define a generalised tetrahedron, and the decomposition of a given polyhedron P into such generalised tetrahedra following a simple combinatorial procedure. We shall consider a polyhedron P both as a combinatorial object, i.e. a trivalent planar graph, and a geometric object. First of all, we deal with combinatorial polyhedra and the associated combinatorial operations or "moves" that are applied to them as to graphs. Namely, we describe two combinatorial operations on the one-skeleton of P , associated with a decomposition into generalised tetrahedra. The first is the classical I -H move, and the second is capping a triangular face with a tetrahedron: Figure 2: A "capping" move These moves have been introduced in [?], and later used in [?, ?] for the purposes of defining a statistical sum over a spin network associated with a hyperbolic polyhedron. By following the general argument for trivalent knotted graphs in [?, §7], it's not hard to see that the following statement holds.

Figure 3: A pentagonal prism Π Lemma 1 Suppose that P is a simple polyhedron, which is not a tetrahedron. Then there exists a sequence of I -H moves and capping moves transforming P into a tetrahedron. From the matrix G we produce the dual (generalised) Gram matrix G * of T , that corresponds to its "vertices" v k , although some of them are outer normals to the respective polar planes P v k . We refer the reader to [?] for the standard definitions of the Gram matrices G and G * , as well as relations between their entries. If v k is ultra-ideal, we associate with it a plane H v k = P v k , otherwise we put H v k = H 3 . Given a plane H n with the outer normal n, let H + n denote the half-space which n points to and let H - n denote the half-space which -n points to. Then, the generalised hyperbolic tetrahedron T is defined as

T = 4 k=1 H - n k ∩ 4 k=1 H - v k .
Thus, each generalised tetrahedron T still can be characterised by six parameters a 1 , a 2 , . . . , a 6 , each of which corresponds to an edge of the "standard" tetrahedron. However, we have to distinguish two cases:

1. if the planes H n i and H n j intersect along an edge e ij with parameter a k , then we set a k = e iα k , where α k is the internal dihedral angle along e ij ;

2. if the planes H n i and H n j are ultra-parallel, then the vertices v k and v l are ultra-ideal and the planes H v k and H v l intersect along the common perpendicular p ij to H n i and H n j : thus we set the respective parameter a k = e -ℓ k , where ℓ k is the length of p ij .

In Fig. ?? the case of a prism-truncated tetrahedron (see [?, ?] and [?]) is illustrated in detail. This type of generalised tetrahedron will appear quite often in the examples below, as well as several other types that we enumerate below.

G =      1 -a 1 +1/a 1 2 -a 2 +1/a 2 2 -a 6 +1/a 6 2 -a 1 +1/a 1 2 1 -a 3 +1/a 3 2 -a 5 +1/a 5 2 -a 2 +1/a 2 2 -a 3 +1/a 3 2 1 -a 4 +1/a 4 2 -a 6 +1/a 6 2 -a 5 +1/a 5 2 -a 4 +1/a 4 2 1      , (1)
as if T were a "usual" tetrahedron: it is still sufficient to determine T uniquely, up to an isometry (see [?, ?] and [?]).

Given a generalised hyperbolic tetrahedron T with Gram matrix (??), we define the following auxiliary quantities in order to state a formula for Vol T . Let U = U (a 1 , a 2 , a 3 , a 4 , a 5 , a 6 , z) denote the function

U = Li 2 (z) + Li 2 (a 1 a 2 a 4 a 5 z) + Li 2 (a 1 a 3 a 4 a 6 z) + Li 2 (a 2 a 3 a 5 a 6 z) (2)
-Li 2 (-a 1 a 2 a 3 z) -Li 2 (-a 1 a 5 a 6 z) -Li 2 (-a 2 a 4 a 6 z) -Li 2 (-a 3 a 4 a 5 z)
depending on seven complex variables a k , k = 1, 6 and z, where Li 2 (•) is the dilogarithm function.

Let z -and z + be two solutions to the equation e z ∂U ∂z = 1 in the variable z. According to [?], these are

z -= -q 1 -q 2 1 -4q 0 q 2 2q 2 and z + = -q 1 + q 2 1 -4q 0 q 2 2q 2 , (3)
where

q 0 = 1+ a 1 a 2 a 3 + a 1 a 5 a 6 + a 2 a 4 a 6 + a 3 a 4 a 5 + a 1 a 2 a 4 a 5 + a 1 a 3 a 4 a 6 + a 2 a 3 a 5 a 6 , q 1 = -a 1 a 2 a 3 a 4 a 5 a 6 a 1 - 1 a 1 a 4 - 1 a 4 + a 2 - 1 a 2 a 5 - 1 a 5 + a 3 - 1 a 3 a 6 - 1 a 6 , (4)
q 2 = a 1 a 2 a 3 a 4 a 5 a 6 (a 1 a 4 + a 2 a 5 + a 3 a 6 + a 1 a 2 a 6 + a 1 a 3 a 5 + a 2 a 3 a 4 + a 4 a 5 a 6 + a 1 a 2 a 3 a 4 a 5 a 6 ).
Given a function f (x, y, . . . , z), let f (x, y, . . . , z) | z=z - z=z + denote the difference f (x, y, . . . , z -) -f (x, y, . . . , z + ). Now we define the following function V = V (a 1 , a 2 , a 3 , a 4 , a 5 , a 6 , z) by means of the equality

V = i 4 U (a 1 , a 2 , a 3 , a 4 , a 5 , a 6 , z) -z ∂U ∂z log z z=z - z=z + . ( 5
)
Depending on the type t|a . . . |p . . . of the tetrahedron T , we consider the following cases.

1. One pair of ultra-parallel faces. This is a prism truncated tetrahedron T , that has the symbol t|aijklm|pn for i, j, k, l, m, n ∈ {1, . . . , 6}. In Fig. ??, we illustrate the case t|a12356|p4. The volume of T can be expressed by the formula

Vol T = ℜ -V + a 4 ∂V ∂a 4 log a 4 , (6)
where V = V (a 1 , . . . , a 6 ) is the "volume function", as described in [?].

2. Two pairs of ultra-parallel faces. In Fig. ??, illustrate the case t|a2356|p14.

Here, all the vertices of T are ultra-ideal. In complete analogy to the proof of [?, Proposition 3], we can show that the volume of T in this case equals The volume of T is

Vol T = ℜ   -V + k∈{1,4} a k ∂V ∂a k log a k   . (7
Vol T = ℜ   -V + k∈{5,6} a k ∂V ∂a k log a k   . (8)
T = ℜ   -V + k∈{4,5,6} a k ∂V ∂a k log a k   . (9)
5. Four pairs of ultra-parallel faces. We illustrate the case t|a14|p2356 in Fig. ??. Here, T has only ultra-ideal vertices and its faces are arranged in

Vol T = ℜ   -V + k∈{2,3,5,6} a k ∂V ∂a k log a k   . (10)
Lemma 2 The volume of a generalised hyperbolic tetrahedron can be expressed by one of the formulae (??)-(??).

The proof of the lemma above is analogous to the proof of [?, Proposition 3] (each of the above formulae are simplified versions of [?, Theorem 1]). We also refer the reader to [?] for a geometric interpretation of the extra terms involving derivatives of V in the volume formulae above.

Let us now reconsider the previous examples of a pentagonal prism Π and a pleated prism P. We shall derive a geometric decomposition for each of them from the associated sequence of combinatorial transformations.

Remark. The geometric decomposition may look quite different from the combinatorial picture. Although it is directed by a sequence of I H and "capping" moves, each such a move result in "chopping off" a generalised tetrahedron from the polyhedron P . The geometric shapes of polyhedra in each step of such a decomposition will vary from what is present in the combinatorial picture where only the 1-skeleton of P is modified as a trivalent graph.

Example. Let us consider the first transformation ih 1 of the prism Π, as depicted in Fig. ??. Let us suppose that the common perpendicular p 67 to the plane 6 and the plane 7 is situated entirely inside Π. Then we draw two planes S 1 and S 2 orthogonal to the sides 1 and 2 of Π, and passing through p 67 . Suppose these planes land on the sides of Π in its interior. Then we cut along Π along the planes S 1 and S 2 , and thus "chop off" a polyhedron T 1 bounded by the planes of the faces 1, 2, 6 and 7 intersecting along the respective edges (or, more precisely, parts of edges) of the initial prism Π, and the planes S 1 and S 2 intersecting along p 67 . Let us note that T 1 is nothing but a generalised tetrahedron of type t|a12345|p6.

The capping transformations cap 1 , cap 2 and cap 3 in Fig. ?? correspond to further "chopping off" of generalised tetrahedra (of type t|a12345|p6):

-the tetrahedron T 2 , formed by the planes 2, 3, 6, 7, and the planes S 2 , S 3 (orthogonal to the faces 2 and 3, respectively),

-the tetrahedron the tetrahedron T 3 , formed by the planes 1, 5, 6, 7, and the planes S 3 , S 5 (orthogonal to the faces 1 and 5, respectively),

-the tetrahedron the tetrahedron T 4 , formed by the planes 3, 4, 6, 7, and the planes S 3 , S 4 (orthogonal to the faces 3 and 4, respectively).

Figure 13: A geometric decomposition of Π into tetrahedra T i , i = 1, . . . , 5, of type t|a12345|p6

In each case we assume that the plane S k lands orthogonally on the k-th face of Π in its interior. Thus, the prism Π is divided into generalised tetrahedra. The last tetrahedron in Fig. ?? corresponds to the generalised tetrahedron T 5 formed by the planes 4, 5, 6, 7 and S 4 , S 5 .

The geometric decomposition of the prism Π following the combinatorial moves in Fig. ?? is depicted in Fig. ??, with all the tetrahedra T i marked. ♦

Example. Let us consider the first transformation ih 1 of the pleated prism P. Let us suppose that the common perpendiculars to any pairs of faces lie inside P, if they exist. Let us take the common perpendicular p 67 to the plane 6 and the plane 7. Then we perform the transformation ih 1 on the edge e 12 (the intersection of the planes 1 and 2). This results in detaching a generalised tetrahedron T 1 , corresponding to the planes 1, 2, 6 and 7 (the normals to these planes determine T 1 completely through the corresponding Gram matrix, however there are two more polar planes bounding it), whose edges consist of the edges involved in the ih 1 move, and the edges coming from the intersection of the planes S 1 and S 2 , orthogonal to the faces 1 and 2, and passing through the common perpendicular p 67 to the plane 6 and the plane 7. The planes S 1 and S 2 intersect along p 67 and land on the faces 1 and 2, respectively, as well as they intersect the planes 6 and 7, yielding seven new edges in total.

Here and below we suppose that each plane S k , orthogonal to the k-th face of P, lands on it orthogonally and inside P.

The transformations cap 1 and cap 2 , which result in capping the triangular faces in Fig. ?? mean further detaching of the generalised tetrahedra T 2 and T 3 , each defined by the outer normals to the four planes 2, 3, 6, 7 and 1, 5, 6, 7, respectively. Now we consider the transformation ih 2 performed on the edge e 48 (the intersection of planes 4 and 8). We take the common perpendicular p 35 to the planes 3 and 5, and again suppose that it lies inside the pleated prism P. We draw two planes passing through p 35 : S 4 orthogonal to the face 4, and S 8 orthogonal to the face 8. Now, we detach one more generalised tetrahedron T 4 defined by the planes 3, 4, 5 and 8.

The transformation c 3 and c 4 correspond to detaching two more tetrahedra: T 5 defined by the planes 3, 5, 6 and 8, and T 6 defined by the planes 3, 5, 4 and 7.

Here we supposed that all possible common perpendiculars are situated inside P, the planes coming through the common perpendiculars orthogonally to the faces of P land in the interior of the faces, and the generalised tetrahedra T k do not overlap. In other words, we suppose that the tetrahedra T k form a decomposition of P: we shall write P = ⊔ k T k . ♦ Below we give examples where the above decomposition is realised for a concrete polyhedron, having the combinatorial type of a pentagonal prism or a pleated prism and prescribed dihedral angles.

Figure 14: A geometric realisation of the prism Π with given dihedral angles, and its decomposition into generalised tetrahedra

Example. Let us consider the pentagonal prism Π depicted in Fig. ??. Each edge of Π is labelled with an label q/p, corresponding to a dihedral angle of pπ/q. In this case, the common perpendicular p 67 to the planes 6 and 7 lies entirely inside the prism Π, and the planes S k , k = 1, . . . , 5 passing through p 67 and orthogonal to the sides of Π subdivide the prism into five "prism truncated" tetrahedra T k (i.e. tetrahedra of type t|a12345|p6), k = 1, 5, which are mutually isometric. In this decomposition all the generalised tetrahedra T k have disjoint interiors, and Π = ⊔ 5 k=1 T k . The volume of each T k , as well as the volume of the prism Π is computed in [?]: Vol T k = 0.52639, Vol Π = 5 • Vol T k ≈ 2.63200. Here we use formula (??) to compute the volumes of all T k , k = 1, . . . , 5. As described in [?], the volume computation proceeds by defining the function Φ(ℓ) = This equation can be easily solved, which allows us to determine the necessary parameter ℓ. ♦ Figure 15: A geometric realisation of the prism P with given dihedral angles, and its decomposition into generalised tetrahedra

Example. Now we consider the pleated prism P depicted in Fig. ?? together with the labels for its dihedral angles. The common perpendiculars p 1 = p 67 (to the planes 6 and 7) and p 2 = p 35 (to the planes 3 and 5) in this case lie inside P, and the planes S k , k = 1, . . . , 8, passing through each of them, respectively, land orthogonally into the interior of the k-th face. Then the pleated prism P is decomposed into seven generalised tetrahedra T k : six are of type t|a12345|p6 (the tetrahedra 1-3 and 5-7), and one is of type t|a2356|p14 (the tetrahedron 4 in the middle). We compute the volume of P. In order to do so, let ℓ k , k = 1, 2, denote the length of p k . Then we have

Vol P = 7 k=1 Vol T k ,
where we use formula (??) for the tetrahedra T k , k ∈ {1, 2, 3, 5, 6, 7} (of type t|a12345|p6) and formula (??) for the tetrahedron T 4 (of type t|a2356|p14).

The tetrahedra T k form the pleated prism P when the angle sum around each p l equals 2π. The dihedral angle along the edge playing the role of p l in each generalised tetrahedron equals

α p l (T k ) = 2a 6 ∂V ∂a 6 (T k ) mod π,
with l = 1 for the tetrahedra T k , k ∈ {1, 2, 3} and l = 2 for the tetrahedra T k , k ∈ {5, 6, 7}, while

α p 1 (T 4 ) = 2a 1 ∂V ∂a 1 (T 4 ) mod π, α p 2 (T 4 ) = 2a 4 ∂V ∂a 4 (T 4 ) mod π. Let us put Φ(ℓ 1 , ℓ 2 ) = 7 k=1 Vol T k + π k∈{1,2} ℓ k .
Then, from the above equations, we deduce that ℓ k is a solution to the system of equations

∂Φ ∂ℓ k = 0, k = 1, 2.
We can transform the above equations into a system of algebraic equations in the variables x k = e -ℓ k :

e ∂Φ ∂ℓ k = 1, k = 1, 2.
This equation is as a polynomial equation in the variables x k , and can be easily solved by computer. Moreover, the functions ∂Φ ∂ℓ k are monotone in each variable ℓ k , k = 1, 2, while keeping the other variable fixed (in the domain where the function Φ is defined). This follows from the fact that the respective dihedral angles of the tetrahedra T k , k = 1, . . . , 7, are monotone functions of ℓ k , k = 1, 2 (cf. Solving the above equations about ℓ k , we obtain that they have a unique common solution (ℓ 1 , ℓ 2 ) ≈ (0.383438, 1.06239), and the volume of P is consequently equal to ∼ 2.34308. Here we use formula (??) to compute the volumes of T k , k ∈ {1, 2, 3, 5, 6, 7}, and formula (??) to compute the volume of T 4 . ♦

Volume of a polyhedron

We prove the following theorem, that allows us for computing the volume of a finite-volume simple polyhedron P , provided it has a decomposition into generalised tetrahedra, e.g. produced by a sequence of I --H and "capping" moves. Here, by a decomposition we mean a subdivision of P into non-overlapping generalised tetrahedra. An analogous synthesis of a combinatorial decomposition followed by a geometric volume computation may be found in [?, ?] as applied to Euclidean polyhedra.

Theorem 1 Let P ⊂ H 3 be a finite-volume simple polyhedron that admits a decomposition into generalised tetrahedra T k . Then the volume of P can be expressed as

Vol P = k Vol T k (a k1 , . . . , a k6 ), (11)
where a kl are the six angle/length parameters corresponding to each generalised tetrahedron, the angle parameters a kl = e iα kl come from the original dihedral angles of the polyhedron, and the edge parameters a kl = e -ℓ kl are determined from a solution to the system of algebraic equations, in which each equation is associated with a common perpendicular p kl to some faces F k and F l of P serving as an edge for a number of generalised tetrahedra

T m , e ∂Φ kl (ℓ kl ) ∂ℓ kl = 1, ( 12
)
with Φ kl (ℓ kl ) = p kl ∩Tm =∅ Vol T m (a m1 , . . . , a m6 ) + π ℓ kl , ( 13
)
where each monomial in the variables a kl has degree at most four.

Proof. We have that

∂Φ kl (ℓ kl ) ∂ℓ kl = - p kl ∩Tm =∅ α kl 2 + π, (14)
where α kl is the dihedral angle along an edge of a generalised tetrahedron T m , corresponding to a common perpendicular p kl to some faces F k and F l of P .

Since the polyhedron P admits a decomposition into generalised tetrahedra T k , then for some length parameters ℓ kl (which will be the same for each tetrahedron T m whose edge ℓ kl coincides with the common perpendicular p kl ) we have

p kl ∩Tm =∅ α kl = 2π. ( 15
)
Thus, (??) implies (??) for these length ℓ kl .

Looking at the expression for the function Vol T k , we observe that ∂Φ kl (ℓ kl ) ∂ℓ kl consists of expressions having the form log(1 + m t ), where m t is a monomial in the parameters a kl . The longest monomial m has the form

a k 1 l 1 a k 2 l 2 a k 3 l 3 a k 4 l 4 .
Thus, if all a kmkn happen to be edge parameters (i.e. a kmkn = e -ℓ kmkn ) then we have an equation of degree four, at most. Let us call a finite-volume polyhedron P ⊂ H 3 "width uniform" if all the common perpendiculars to its facets lie inside P . We also recall that P ⊂ H 3 is called a Coxeter polyhedron if all its dihedral angles are integer submultiples of π, i.e. have the form π/m, for an integer m ≥ 2. These polyhedra form an important class widely used in further constructions of hyperbolic orbifolds and manifolds [?]. The following fact obviously holds as a consequence of geodesic convexity.

Lemma 3 Every Coxeter polytope (or, more generally, every acute-angled polyhedron) is width uniform.

A width uniform polyhedron has a geometric decomposition into generalised tetrahedra. We also notice that each polytope P may be decomposed into a number of generalised tetrahedra T k following a sequence of I -H and "capping" moves, although Theorem ?? may not apply directly. E.g. some of the T k 's may overlap, and thus we shall need to use the "inclusion-exclusion" formula in order to express the volume of P through the volumes of T k . As well, some of the polar hyperplanes bounding T k may land on the faces of P entirely or partially outside P forming "butterfly" generalised tetrahedra, some of which are depicted in Fig. ??,c.f. [?].

The volume formula of [?] works in either case, however, the glueing equations (??)-(??) will not have a unique solution. Still, the value of the volume of P can be found by using some solution among (usually few) solutions of the glueing equations.

Figure 17: Some examples of "butterfly" generalised tetrahedra

Further computational examples

Below we give more numeric examples of volume computation for several polytopes. We start with the pleated prism P, and then proceed to a dodecahedron D. We shall use Coxeter dihedral angles (i.e. dihedral angles of the form π/n, for n a natural number ≥ 2), so we can "double-check" the volume of each polytope P with the following procedure: the volume of P equals the volume of the reflection orbifold obtained from P by mirroring all its faces. We create a "double" of the respective reflection orbifold, which is an orientable orbifold O(P ) with underlying space the three-sphere S 3 , and singular set combinatorially isomorphic to the 1-skeleton of P . Each singular edge e carries an angle of 2π/n, corresponding to the angle of π/n in the original polytope P . The volume of O(P ) can be computed by the "Orb" software [?], and Vol O(P ) = 2 VolP . First, we pick the pleated prism P from our previous example, and reduce it to a tetrahedron by a different sequence of I -H and capping moves depicted in Fig. ??. Thus, the position of the common perpendiculars in the new subdivision will be different.

As shown in Fig. ??, we have again 7 tetrahedra to consider, five of which are of the type t|a12356|p4 (the tetrahedra T i , i ∈ {1, 2, 3, 6, 7}), and two are of the type t|a1234|p56 (the tetrahedra T i , i ∈ {4, 5}). Here, our computations find that the length of the common perpendicular p 1 = p 67 to the planes 6 and 7 is ℓ 1 ≈ 0.383438, and the length of the common perpendicular p 2 = p 47 to the planes 4 and 7 is ℓ 2 ≈ 0.626516. The volume of the prism P equals, as before, Vol P ≈ 2.34308. By using various sequences of I -H and capping moves, and therefore where Γ, c (r) is the unitary spin network (the Kirillov-Reshetikhin invariant with framing c (r) , r = 3, 5, . . . ) associated with Γ, c.f. [?, ?] and the root of unity is specified at q = exp 4πi r . First we recall some necessary basic definition that are used in our computation below. Let r ≥ 3 be a odd integer, q = exp(4πi/r) and the quantum number {n} and factorial {n}! be defined as {n} = q n/2 -q -n/2 , {n, k} = {n} {n-1} • • • {n-k+1}, {n}! = {n, n}.

z {z + 1}! {1} {z -A 1 }! {z -A 2 }! {z -A 3 }! {z -A 4 }! {B 1 -z}! {B 2 -z}! {B 3 -z}! ,
where

A 1 = a + b + e, A 2 = a + c + f, A 3 = b + d + f, A 4 = c + d + e, B 1 = a + b + c + d, B 2 = a + d + e + f, B 3 = b + c + e + f, ∆(x, y, z) = {-x + y + z}! {x -y + z}! {x + y -z}! {1} {x + y + z + 1}! .
and z runs from max(A 1 , A Example. Let Π be the pentagonal prism depicted in Fig. ??. Let Γ be the planar graph given by the 1-skeleton of Π, with edges e 1 , . . . , e 15 , and the corresponding dihedral angles α 1 , . . . , α 15 along them. Here, we have that α i ∈ π 2 , π 3 , 2π 5 . Let c (r) be the sequence of colourings s (r)

x + y + z ∈ Z, x + y + z ≤ r -2, 0 ≤ x + y -z, x -y + z, -x + y + z.
i , i = 1, . . . , 15, r = 3, 5, 7, . . . of the edges e i of Γ such that the conditions of Conjecture ?? are satisfied. This means that we have s

(r) i = ⌊ r-3 4 ⌋ for each α i = π 2 , s (r) i = ⌊ r-3
3 ⌋ for each α i = π 3 , and s (r) i = ⌊ 3(r-3) 10 ⌋ for each α i = 2π 5 . Then the unitary spin network associated with Γ evaluates into

Γ, c (r) = k (-1) k {k + 1} {1}   ⌊ r-3 3 ⌋ ⌊ r-3 3 ⌋ ⌊ 3(r-3) 10 ⌋ ⌊ r-3 4 ⌋ ⌊ r-3 4 ⌋ k RW q=exp( 4πi r )   5 .
In Table ?? we collect the results of numerical computations of the quantity Analogous, however way more tedious computations may be carried out in the case of the hyperbolic dodecahedra D i , i = 1, 2, 3, described in Section ??. ♦ Each time an I -H move or a capping move is performed on Γ, we have a generalised tetrahedron detached from the polyhedron P . As a step towards verification of Conjecture ?? we suggest investigating if the following statements hold.

V (r) = 2π r log Γ, c (
Conjecture 2 Let α, β, γ, δ, ǫ, φ be the dihedral angles of a generalised hyperbolic tetrahedron T , which may have ideal and ultra-ideal vertices, though no truncating polar planes determined by its ultra-ideal vertices intersect. Let a

(r) , b (r) , • • • , f (r) , r = 3, 5, 7, • • • be sequences of integers satisfying 4 π lim r→∞ a (r) r = π -α, 4 π lim r→∞ b (r) r = π -β, 4 π lim r→∞ c (r) r = π -γ, 4 π lim r→∞ d (r) r = π -δ, 4 π lim r→∞ e (r) r = π -ǫ, 4 π lim r→∞ f (r) r = π -φ,
where r is always an odd integer. Then *

2 π lim r→∞ 1 r log a (r) b (r) e (r) d (r) c (r) f (r) RW q=exp( 4πi r ) = Vol(T ). (16)
In Fig. ?? we illustrate the values of V (α) = Vol T α , for a regular tetrahedron with all dihedral angles equal to α. The tetrahedron T α is spherical for α > arccos(1/3), Euclidean for α = arccos(1/3), and hyperbolic for π/3 < α < arccos(1/3). It is ideal for α = π/3, truncated (i.e. has ultra-ideal vertices) hyperbolic for 0 < α < π/3, and it deforms into a regular ideal octahedron O for α = 0. The octahedron O has all dihedral angles equal to π/2 and its volume is Vol O ≈ 3.663.... Let us notice that the tetrahedron considered in the above conjecture is supposed to have only mild vertex truncations, c.f. [?]. Another conjecture considers a generalised hyperbolic tetrahedron T with some vertices * We suppose that the sequences a (r) , . . . , f (r) are such that the corresponding 6j-symbol is well-defined. intensely truncated. We specify this conjecture for the case of a prism truncated tetrahedron, considered in [?, ?]. Although, one may state an analogous conjecture for a generalised tetrahedron of any type considered in Section ??. {(j + 1) (2 f (r) + 1)} {j} = Vol(T ). ( 17) † However, only two truncating polar planes determined by its ultra-ideal vertices intersect.

‡ We suppose that the sequences a (r) , . . . , f (r) are such that the corresponding 6j-symbol is well-defined.

Let T α,β be the doubly truncated tetrahedron with dihedral angle α at its usual five edges and β at the edge arising from intense truncation. Let U (k, l) = Since the volume of P is a sum of volumes of the generalised tetrahedra T i composing it, we hope that there might be a correspondence between the asymptotic behaviour of the sums (??)-(??) associated with each tetrahedron T i , and the asymptotic behaviour of the Kirillov-Reshetikhin invariant Γ, c (r) of the whole polyhedron P .

Figure 1 :

 1 Figure 1: An I -H move

Example.

  Let us take a pentagonal prism Π depicted in Fig. ?? and apply a sequence of combinatorial transformations to its one-skeleton, as shown in Fig. ??.

Figure 4 :Figure 5 :Figure 6 :

 456 Figure 4: A sequence of I -H and capping moves applied to the prism Π Thus, having applied a short sequence of combinatorial moves we obtain a tetrahedron. ♦

Figure 7 :

 7 Figure 7: A prism-truncated tetrahedron

Figure 8 :

 8 Figure 8: A tetrahedron of type t|a12356|p4

Figure 9 :

 9 Figure 9: A tetrahedron of type t|a2356|p14

) 3 .

 3 Two pairs of ultra-parallel faces, another configuration. A tetrahedron of type t|a1234|p56 is depicted in Fig.??. Here, three vertices v 1 , v 2 and v 3 are ultra-ideal and the vertex v 4 is proper.

Figure 10 :

 10 Figure 10: A tetrahedron of type t|a1234|p56

Figure 11 : A tetrahedron of type t|a123|p456 4 .

 114 Figure 11: A tetrahedron of type t|a123|p456 4. One face that is ultra-parallel to three others. In Fig. ??, a tetrahedron of type t|a123|p456 is shown. Here, three faces are adjacent around the proper vertex v 4 , while three other vertices v 1 , v 2 and v 3 are ultra-ideal. The volume of T is

Figure 12 :

 12 Figure 12: A tetrahedron of type t|a14|p2356

5 k=1

 5 Vol T k + πℓ, where ℓ is the length of the common perpendicular p 56 , and solving the equation ∂Φ(ℓ) ∂ℓ = 0.We can transform the last equation into an algebraic equation about x = e -ℓ :

  Fig. ??).

Figure 16 :

 16 Figure 16: The dihedral angles of the tetrahedron T 4 as fucntions of edge lengths ℓ 1 and ℓ 2 (the graph is symmetric under exchanging ℓ 1 and ℓ 2 )

Figure 18 :

 18 Figure 18: Another sequence of I -H and capping moves reducing the pleated prism P to a tetrahedron

Figure 19 :Figure 21 :Figure 23 :

 192123 Figure19: A geometric realisation of the pleated prism P with given dihedral angles, and its decomposition into generalised tetrahedra

.

  The Racah-Wigner quantum 6j-symbol is defined for six half-integers a, b, c, d, e, f in the set {0, 1/2, 1, 3/2, • • • , (r -2)/2} and the actual formula is given by [?] as follows: a, b, e) ∆(a, c, f ) ∆(b, d, f ) ∆(c, d, e)× z (-1)

  The actual computation of the Kirillov-Reshetikhin invariant of a coloured trivalent graph Γ with colouring c follows the skein relations (4.1)-(4.13) and (5.4), (5.6) described in [?].

  Figure 25: Values of the quantum 6j-symbol and V (α).

Conjecture 3

 3 Let α, β, γ, δ, ǫ, φ be the dihedral angles of a prism truncated hyperbolic tetrahedron T , which may have ideal or ultra-ideal vertices. † Let the parameter φ represent the dihedral angle at the edge arising from intense truncation. Let a (r) , b (r) , • • • , f (r) (r = 3, 5, 7, • • • ) be sequences of integers satisfying4 π lim r→∞ a (r) r = π -α, 4 π lim r→∞ b (r) r = π -β, 4 π lim r→∞ c (r) r = π -γ, 4 π lim r→∞ d (r) r = π -δ, 4 π lim r→∞ e (r) r = π -ǫ, 4 π lim r→∞ f (r) r = π -φ,where r is always an odd integer. Then ‡ ) b (r) e (r) d (r) c (r)

  j + 1)(2 l + 1)} {2j + 1} .In Fig.?? the values of U ((r -3)/3, l) and V (β) = Vol T π/3,β are shown.

Figure 26 :

 26 Figure 26: Values of U ( r-3 3 , l) and V (β).

  2 , A 3 , A 4 ) to min(B 1 , B 2 , B 3 ). Above we suppose that the triples (a, b, e), (a, c, f ), (b, d, f ) and(c, d, e) are all admissible. Here, a triple of half integers (x, y, z) in {0, 1/2, • • • , (r-2)/2} is said to be admissible if x, y, z satisfy

Table 1 :

 1 r) for various odd values of r. The values of Kirillov-Reshetikhin invariants for Π We also perform numerical computations for the pleated prism P depicted in Fig.??. The numeric values of its Kirillov-Reshetikhin invariants are given in Table?? for various odd values of r.

	r	483	963	1923	3843	Vol Π
	V (r) 2.27094 2.42388 2.51421 2.56627 ≈ 2.63200

Table 2 :

 2 The values of Kirillov-Reshetikhin invariants for P


Pierre Saramito
email: pierre.saramito@imag.fr

A new constitutive equation for elastoviscoplastic fluid flows

Keywords: non-Newtonian fluid, viscoelasticity, viscoplasticity, constitutive equation

From thermodynamic theory, a new three-dimensional model for elastoviscoplastic fluid flows is presented. It extends both the Bingham viscoplastic and the Oldroyd viscoelastic models. Fundamental flows are studied: simple shear flow, uniaxial elongation and large amplitude oscillatory shear. The complex moduli (G ′ , G ′′ ) are found to be in qualitative agreement with experimental data for materials that present microscopic network structures and large scale rearrangements. Various fluids of practical interest, such as liquid foams, droplet emulsions or blood, present such elastoviscoplastic behavior : at low stress, the material behaves as a viscoelastic solid, whereas at stresses above a yield stress, the material behaves as a fluid.

Introduction

Historical background

The development of viscoplastic rheological model based on yield stress started in 1900 when Schwedoff [START_REF] Schwedoff | La rigidité des fluides[END_REF], studying a gelatin suspension, presented a one-dimensional plastic viscoelastic version of the Maxwell model:

   ε = 0, when τ ≤ τ 0 , λ dτ dt + (τ -τ 0 ) = η m ε, when τ > τ 0 , (1)
where τ is the stress, ε the rate of deformation, η m > 0 the viscosity, τ 0 ≥ 0 the yield stress and λ ≥ 0 a relaxation time. In steady shear flow this reduces to τ = τ 0 + η m ε when τ > τ 0 . In 1922, Bingham [START_REF] Bingham | Fluidity and plasticity[END_REF] proposed the one-dimensional stress-deformation rate equation for a viscous fluid with a yield stress:

max 0, |τ | -τ 0 |τ | τ = η m ε ⇐⇒      |τ | ≤ τ 0 when | ε| = 0, τ = η m ε + τ 0 ε | ε| otherwise. (2)
Notice that this model is equivalent -up to the sign of ε, assumed positive -to the steady case of the model proposed by Schwedoff. Numerous attempts have been made to modify this simple equation to account for more complex behavior of such materials. In 1926, Herschel and Bulkley [START_REF] Herschel | Measurement of consistency as applied to rubber-benzene solutions[END_REF] proposed to model the observed shear stress dependence of the viscosity on the shear rate ε after yielding by explicitly defining the viscosity η m as a power-law function of | ε|.

In 1932, Prager [START_REF] Prager | Über die Ansätze der Mechanik isotroper Kontinua[END_REF], using the von Mises [START_REF] Mises | Mechanik der festen Körper im plastich-deformablen Zustand[END_REF] yielding criterion, proposed to extend the Bingham model to the three-dimensional case and Oldroyd, in 1947, in a collection of papers (see e.g. [START_REF] Oldroyd | A rational formulation of the equations of plastic flow for a Bingham solid[END_REF]) studied the Bingham three-dimensional model and its Herschel-Bulkley extension coupled with the Navier-Stokes equations for the motion of the fluid. Oldroyd also proposed a three-dimensional constitutive equation which combines the yielding criterion together with a linear Hookean elastic behavior before yielding and a viscous behavior after yielding. In the one-dimensional case, the model can be written as:

     τ = µε when τ ≤ τ 0 , |τ | -τ 0 |τ | τ = η m ε when τ > τ 0 . (3)
When compared to [START_REF] Bingham | Fluidity and plasticity[END_REF], this model is an improvement, since the material is no longer rigid before yielding.

Here, the yield stress τ 0 is related to the critical strain ε 0 = τ 0 /µ. Since the first equation in [START_REF] Doraiswamy | The Cox-Merz rule extended: a rheological model for concentrated suspensions and other materials with a yield stress[END_REF] describes stresses in term of strain ε and the second equation in [START_REF] Doraiswamy | The Cox-Merz rule extended: a rheological model for concentrated suspensions and other materials with a yield stress[END_REF] in term of strain rate ε, the stress-strain curve predicted by this model must exhibit a discontinuity at the critical strain ε = ε0 at which the stress jumps from τ = τ 0 to τ = τ 0 + η ε. This is an approximation of the true behavior of materials: the real deformation at the transition is expected to be smooth, at least continuous.

In 1950, Oldroyd [START_REF] Oldroyd | On the formulation of rheological equations of states[END_REF] developed a theory for the invariant forms of rheological equations of state and proposed a three-dimensional viscoelastic model, that can be expressed in its one-dimensional version as:

λ dτ dt + τ = η m ε (4)
where the total stress σ = η ε + τ . In this approach the stress τ is the elastic part of the total stress, from which the elastic deformation can easily be obtained. The constant η > 0 is another viscosity, often called the solvent viscosity in the context of polymer solutions.

One-dimensional presentation of the proposed model

Let us introduce an one-dimensional version of our model in order to combine the two previous models (2) and ( 4):

λ dτ dt + max 0, |τ | -τ 0 |τ | τ = η m ε. ( 5
)
where the total stress is σ = η ε + τ . When λ = η = 0 we obtain (2) while when τ 0 = 0 our model reduces to [START_REF] Gordon | Anisotropic fluid theory: a different approach to the dumbbell theory of dilute polymer solutions[END_REF]. Observe that (5) differs both from ( 1

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 µ η σ ε 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 η σ ε εe εm µ ηm (c) present model (d) creeping σ(t) = σ for the present model η σ τ0 ε εe ηm εm µ σ < τ 0 σ > τ 0 t ε(t) 0 0 Figure 1. The proposed elastoviscoplastic model. is σ = µε + η ε + τ .
As soon as the strain energy exceeds the level required by the von Mises criterion, the stress in the friction element attains the yield value and the element breaks allowing deformation of all the other elements. After yielding, the deformation of these elements describes the Oldroyd-type viscoelastic behavior.

The evolution in time of elongation ε(t) for a fixed imposed traction σ (creeping) is represented on Fig 1 .d. When σ ≤ τ 0 , the elongation for a fixed imposed traction is bounded in time, which means that such material behaves as a solid. Otherwise, σ > τ 0 , the elongation is unbounded in time, which means that the material behaves as a fluid. Notice the continuous behavior of the solution: the change of regime occurs at t = t 0 for the solution σ > τ 0 and there is neither jump nor loss of differentiability in the curve. The solution on Fig 1 .d is computed explicitly and details are reported in appendix B.

Comparison with other recent and closely related models

A number of other closely related models have appeared in the literature.

In 1991, Beris et al.. [START_REF] Doraiswamy | The Cox-Merz rule extended: a rheological model for concentrated suspensions and other materials with a yield stress[END_REF], in order to recover a continuous approximation of the solution of the elasticviscoplastic model (3) proposed by Oldroyd, introduced an ad-hoc recovery procedure. Before yielding, the material behaves as an elastic solid while after yielding it behaves as a power-law viscous non-Newtonian fluid. Despite the lack of a thermodynamical analysis of their model, these authors have been able to propose, based on their computational results, a useful Cox-Merz rule extension that was found to be in good agreement with experimental data on a suspension of silicon particles in polyethylene.

In 2003, Puzrin and Houlsby [15, p. 254] inserted friction and a spring into the viscoelastic Kelvin-Voigt model and the resulting model is represented in Fig. 2.a. Before yielding, the friction is rigid and the material behaves as an elastic solid, thanks to the spring µ. After yielding, the deformation is governed by the so-called standard model: a Kelvin-Voigt element plus a spring in series (see e.g. [10, p. 42]). The predicted elongation under a constant traction is presented on Fig. 2.b. The elongation jumps immediately at t = 0 from ε = 0 to ε = σ/µ, since the spring element has no time scale. Next, when t > 0, if σ ≤ τ 0 , then the elongation remains constant, since the spring is fully extended. Otherwise, if σ > τ 0 , the elongation grows and tends to a bounded value: the material behaves as a solid. By this way Puzrin and Houlsby were able to capture some relevant aspects of the behavior of saturated clays. These authors have proposed in a collection of papers [START_REF] Houlsby | Rate-dependent plasticity models derived from potential functions[END_REF][START_REF] Puzrin | Rate-dependent hyperplasticity with internal functions[END_REF] many variants of their model. Nevertheless, it is not applicable to material such as human blood or liquid foams that deform under a low stress and flow under a sufficient stress.

In 1990 Isayev and Fan [START_REF] Isayev | Viscoelastic plastic constitutive equation for flow of particle filled polymers[END_REF] proposed a viscoelastic plastic constitutive equation for flow of particles filled polymers. The main idea of these authors is to insert a friction and a spring into the viscoelastic Oldroyd model and the resulted model is represented on Fig. 2.c. Before yielding, the friction is rigid and the material behaves as a simple elastic solid, thanks to the spring µ. After yielding, the deformation is described by the Oldroyd viscoelastic model. These authors have added in parallel other Maxwell elements in order to replace the Oldroyd viscoelastic element by a multi-mode Leonov viscoelastic model. The predicted elongation under a constant traction is presented on Fig. 2.d. As for the Puzrin-Houlsby model, the elongation jumps immediately at t = 0 from ε = 0 to ε = σ/µ. Next, when t > 0, if σ ≤ τ 0 , then the elongation remains constant. Otherwise, if σ > τ 0 , the elongation grows and is not bounded: in that case, the material behaves as a fluid. The friction element was able to describe the stress generated in the disperse phase of the filled polymer melts and this contribution has represented an important conceptual advance in this domain. Since the effect of added solid particles in a polymer melt is usually to reduce the viscoelasticity, the need of this kind of formulation arises primarily in the case of highly viscoelastic polymers such as rubbers. Nevertheless, the model presented in Fig. 1.c and. 1.d is a definite improvement: the instantaneous jump at t = 0 followed by a constant elongation under traction of both the Puzrin-Houlsby and the one-mode simplified Isayev-Fan models is replaced by a smoother behavior.

During the last two decades, the thermodynamic framework developed and furnished some robust theoretical tools for managing efficiently rheological models: see e.g. [22, p. 26] or [START_REF] Maugin | The thermomechanics of plasticity and fracture[END_REF] for concise presentations of the thermodynamic framework for the development of combined viscoelastic and plastic models. The table 3 presents an historical summary of such models and signals their associated main features before and after yielding. The 3D column is marked when the model has been written in a general sense, e.g. with objective derivatives. The TH column is also marked when the model satisfies the second law of thermodynamics. The italic mark X in the 3D or the TH collumns means that the corresponding property has been showed later. For instance, the Bingham and the Herschel-Bulkley models was first proposed in a one-dimensional context, next extended to 3D and finally founded to satisfy the second law of thermodynamics.

The aim of the present article is to build the proposed model for the general three-dimensional case (section 2) and to study it on three simple flows (section 3): a simple shear flow, an uniaxial elongation and an oscillatory shear flow. The impatient reader -and the reader who is unfamiliar with the thermodynamic framework -could jump directly to the end of section 2 where the complete set of equations ( 11) governing such a flow is presented, before reading section 3 devoted to applications.

There are also two technical appendices. Appendix A contains some properties used in section 2 and related to multi-dimensional subdifferential calculus. Appendix B groups the explicit resolution of the one-dimensional creeping problems as presented in the introduction on Figs. 1 and 2.

The proposed model

Thermodynamic framework

The state of the system is described by using two independent variables : the total deformation tensor ε and an internal variable, the elastic deformation tensor ε e . We have ε = ε e + ε m where ε m represents the plastic deformation tensor. Following Halphen and Nguyen [START_REF] Halphen | Sur les matériaux standard généralisés[END_REF] (see e.g. [START_REF] Tallec | Numerical analysis of viscoelastic problems[END_REF] or [10, p. 97]) we say that a generalized standard material is characterized by the existence of a free energy function E and a potential of dissipation D, that are both convex functions of their arguments. The proposed model that combines viscoelasticity and viscoplasticity, is represented by a rheological scheme on Fig. 1.c and can be written as:

E (ε, ε e ) = µ |ε e | 2 , D ( ε, εe ) = ϕ ( ε ) + ϕ m ( ε -εe ) , (6)
where µ > 0 is the elasticity parameter and where |.| denotes the matrix norm, defined by a double contraction of indices : |ε e | 2 = ε e : ε e . The functions ϕ and ϕ m are expressed by :

ϕ ( ε) =    η | ε | 2 when tr ε = 0, +∞ otherwise, and ϕ m ( εm ) =    η m | εm | 2 + τ 0 | εm | when tr εm = 0, +∞ otherwise. ( 7
)
The ϕ function expresses the incompressible viscous behavior at macroscopic level and is associated to a macroscopic viscosity η > 0 while the ϕ m function expresses the viscoplastic behavior at microscopic level, by using a microscopic viscosity η m > 0, acting on continuous modification of the network links. and also a yield stress value τ 0 ≥ 0. When the stress becomes higher than this value, some topological modifications appear in the network of contacts. This model satisfies the second law of thermodynamics: in the framework of generalized standard materials [START_REF] Halphen | Sur les matériaux standard généralisés[END_REF][START_REF] Tallec | Numerical analysis of viscoelastic problems[END_REF][START_REF] Maugin | The thermomechanics of plasticity and fracture[END_REF] this property is a direct consequence of the convexity of both E and D.

The constitutive law

Let Ω be a bounded domain of R N , where N = 1, 2, 3. Since both ϕ and ϕ m are non-linear and nondifferentiable, the following manipulations involve subdifferential calculus from convex analysis. The material constitutive laws can be written as:

σ ∈ ∂E ∂ε + ∂D ∂ ε and 0 ∈ ∂E ∂ε e + ∂D ∂ εe , ( 8
)
where σ is the total Cauchy stress tensor. Using definition ( 6) of E and D, we get:

σ ∈ ∂ϕ ( ε) + ∂ϕ m ( ε -εe ) and 0 ∈ 2µε e -∂ϕ m ( ε -εe ) . (9)
The combination of the two previous relations leads to σ -2µε e ∈ ∂ϕ( ε). Then, by using expression (A.3) of ∂ϕ from the technical annex, and by introducing the pressure field p, we get the following expression of the total Cauchy stress tensor: σ = -p.I + 2η ε + 2µε e when tr( ε) = 0. Then, the second relation in ( 9) is equivalent to εεe ∈ ∂ϕ * m (2µε e ) where ϕ * m is the dual of ϕ m . Let us introduce the elastic stress tensor τ = 2µε e . The expression (A.2) of ∂ϕ * m in annex yields:

λ τ + max 0, |τ d | -τ 0 |τ d | τ = 2η m ε, (10)
where λ = η m /µ is the relaxation time and τ d = τ -1 N tr(τ ) I denotes the deviatoric part of τ .

The system of equations

Since the material is considered in large deformations, we choose to use the Eulerian mathematical framework, more suitable for fluids flows computations. We assume that ε = D(v) = ∇v + ∇v T /2 is the rate of deformation, while the material derivative τ of tensor τ in the Eulerian framework is expressed by the Gordon-Schowalter derivative [START_REF] Gordon | Anisotropic fluid theory: a different approach to the dumbbell theory of dilute polymer solutions[END_REF] :

2 τ = ∂τ ∂t + v.∇τ + τ W (v) -W (v)τ -a (τ D(v) + D(v)τ ) where W (v) = ∇v -∇v T /2
is the vorticity tensor. The material parameter a ∈ [-1, 1] is associated to the Gordon-Schowalter's derivative. When a = 0 we obtain the Jaumann derivative of tensors, while a = 1 and a = -1 are associated to the upper and the lower convected derivatives, respectively. The elastoviscoplastic fluid is then described by a set of three equations associated to three unknowns (τ, v, p): the differential equation ( 10) is completed with the conservation of momentum and mass:

             λ 2 τ + max 0, |τ d | -τ 0 |τ d | τ -2η m D(v) = 0, ρ ∂v ∂t + v.∇v -div (-pI + 2ηD(v) + τ ) = f , div v = 0,
where ρ denotes the constant density and f a known external force, such as the gravity. These equations are completed by some suitable initial and boundaries conditions in order to close the system. For instance the initial conditions τ (t = 0) = τ 0 and v(t = 0) = v 0 and the boundary condition v = v Γ on the boundary ∂Ω are convenient. The total Cauchy stress tensor can be written as:

σ = -pI + 2ηD(v) + τ.
Notice that when τ 0 = 0 the model reduces to the usual viscoelastic Oldroyd model [START_REF] Oldroyd | On the formulation of rheological equations of states[END_REF][START_REF] Saramito | Operator splitting for viscoelastic fluid with a differential constitutive law[END_REF] and when λ = 0 it reduces to the viscoplastic Bingham model [START_REF] Bingham | Fluidity and plasticity[END_REF][START_REF] Roquet | An adaptive finite element method for Bingham fluid flows around a cylinder[END_REF]. When both τ 0 = 0 and λ = 0 the fluid is Newtonian and the set of equations reduces to the classical Navier-Stokes equations with a total viscosity η 0 = η + η m . Conversely, when both τ 0 = 0 and λ = 0 the fluid is elastoviscoplastic.

Dimensionless formulation -Let us introduce some classical dimensionless numbers:

W e = λU L , Bi = τ 0 L η 0 U and Re = ρU L η 0 ,
where U and L are some characteristic velocity and length of the flow, respectively. We also use the retardation parameter α = η m /η 0 . The problem reduces to that of finding dimensionless fields, also denoted by (τ, v, p) such that:

             W e 2 τ + max 0, |τ d | -Bi |τ d | τ -2αD(v) = 0, Re ∂v ∂t + v.∇v -div (-pI + 2(1 -α)D(v) + τ ) = f , div v = 0, (11)
where f denotes some known dimensionless vector field. These equations are completed by the initial and boundaries conditions.

Examples

Uniaxial elongation

The fluid is at the rest at t = 0 and a constant elongational rate ε0 is applied: the Weissenberg number is W e = λ ε0 and the Bingham number Bi = τ 0 /(η 0 ε0 ). All quantities presented in this paragraph are dimensionless.

The flow is three-dimensional and the dimensionless velocity gradient writes ∇v = diag(1, -1/2, -1/2).

The problem reduces to find τ 11 , τ 22 and τ 33 such that with the initial condition τ (t = 0) = 0 and where κ = max (0, 1 -Bi/|τ d |). As above, since τ (0) = 0 and τ (t) is continuous, there exists t 0 > 0 such that when t ∈ [0, t 0 ] we have |τ d | ≤ Bi and thus κ = 0: this is the linear flow regime. The eigenvalues of the system are -2aW e and aW e. For t > t 0 , the case κ > 0 occurs. Since κ ≤ 1, when aW e > 1/2 the stress becomes unbounded in finite time, as shown on Fig. 4.a that plots the dimensionless first normal stress difference ψ = τ 11τ 22 . This drawback is still true when Bi = 0, i.e. for the Oldroyd viscoelastic model. In the context of viscoelastic models, some alternate constitutive equations that extend the Oldroyd model have been proposed. Let us consider the following constitutive equation:

     W e dτ 11 dt + (κ - 2aW
W e 2 τ + (1 + ξ tr τ ) max 0, |τ d | -Bi |τ d | τ -2αD(v) = 0, (12)
where ξ ≥ 0 is a new material parameter. For Bi = 0 we obtain the viscoelastic Phan-Thien and Tanner model (see e.g. [START_REF] Phan-Thien | A new constitutive equation derived from network theory[END_REF][START_REF] Saramito | Efficient simulation of nonlinear viscoelastic fluid flows[END_REF]). For ξ = 0 we obtain the previous elastoviscoplastic model [START_REF] Oldroyd | A rational formulation of the equations of plastic flow for a Bingham solid[END_REF] while for any ξ > 0 the solution remains bounded, as shown on Fig. 4.b for W e = 0.75.

Simple shear flow

The fluid is at the rest at t = 0 and a constant shear rate γ0 is applied: the Weissenberg number is W e = λ γ0 and the Bingham number Bi = τ 0 /(η 0 γ0 ).

The flow is two-dimensional and the dimensionless velocity gradient is constant: ∇v = ([0, 1]; [0, 0]). The problem reduces to finding τ 11 , τ 22 and τ 12 , such that, for all t > 0: and τ (t) is continuous, there exists t 0 > 0 such that when t ∈ [0, t 0 ] we have |τ d | ≤ Bi and thus κ = 0: this is the linear flow regime. The eigenvalues of the system are 0 and ±i √ 1a 2 . At t = t 0 , |τ d | reaches Bi. Then, for t > t 0 , the non-linear factor κ > 0 occurs: the corresponding term reduces the growth of the solution, which now remains bounded. When a = 1 (see Fig. 5.a, where τ 22 = 0) and when a = 0 and Bi (see Fig. 5.b, where τ 22 = -τ 11 ) is small enough, the solution tends to a constant. Remark the important overshoot of both τ 12 and τ 11τ 22 when a = 0 (Fig. 5.b) while when a = 1 there is only a tiny overshoot of τ 12 . When a = 0 and Bi becomes large, instabilities appear, while the solution remains bounded (Fig. 5.c). There has been many studies of the a parameter of the Gordon-Schowalter derivative (see e.g. [START_REF] Vlassopoulos | A generalized Giesekus constitutive model with retardation time and its association to the spurt effect[END_REF]). Based on comparisons with data measurements, most authors consider that only the case a = 1 -or at least a close to 1 -is relevant from a practical point of view.

W e d dt      τ 11 τ 22 τ 12      + W e      0 0 -(1 + a) 0 0 1 -a 1 -a 2 - 1 + a 2 0           τ 11 τ 22 τ 12      + κ      τ 11 τ 22 τ 12      =      0 0 α      , (13)
Notice that the dimensionless steady shear viscosity η s /η 0 coincides with the dimensionless shear stress σ 12 = 1α + τ 12 . When the solution becomes stationary, let us observe the dimensionless steady shear viscosity η s /η 0 as a function of W e on Fig. 6. The material presents a shear thinning character. For large W e, the shear viscosity tends to a plateau when ξ = 0 (Fig. 6.a) and decreases monotonically when ξ > 0 (Fig. 6.b). This shear thinning behavior is more pronounced when Bi increases: the value of Bi controls the plateau at small values of W e.

Periodic shear flow

An oscillatory shear flow is applied: the imposed shear strain is γ 0 sin(ωt) and the shear rate becomes γ 0 ω cos(ωt). The Weissenberg number is defined by W e = λω and the Bingham number by Bi = τ 0 /(η 0 γ 0 ω).

The flow is two-dimensional and the velocity gradient is periodic: ∇v(t) = ([0, cos t]; [0, 0]). The problem reduces to that of finding τ 11 , τ 22 and τ 12 , from R + to R, such that

Bi = 0 Bi = 1 Bi = 2 Bi = 4 (a) η s /η 0 W e -2
10 -1 10 0 10 1 10 2

Bi = 0

Bi = 1 Bi = 2 Bi = 4 (b) η s /η 0 W e 10 -2
10 -1 10 0 10 1 10 2 10 0

               W e d dt      τ 11 τ 22 τ 12      + F(t, τ ) =      0 0 α cos t      in ]0, 2π[, τ (0) = τ (2π), (14)
where

F(t, τ ) = W e cos(t)A a      τ 11 τ 22 τ 12      + κ      τ 11 τ 22 τ 12      and κ = max (0, 1 -Bi/|τ d |) and |τ d | 2 = (1/2) (τ 11 -τ 22 ) 2 + 2τ 2 12 .
The ξ material parameter introduced in (12) has been found to have an influence on the elongational flow, while in the case of the periodic shear flow, it has been found to have very minor influence. Thus only the case ξ = 0 is presented in this section.

The solution for Bi = 0

The case Bi = 0 is asymptotically reached when the amplitude γ 0 becomes large. The material behaves as a viscoelastic fluid and the system reduces to: together with periodic boundary condition τ (0) = τ (2π). By introducing the first normal stress difference ψ = τ 11τ 22 , the system can be rewritten in a more compact form:

           W e
                 W e dψ dt + ψ = 2W e cos(t)τ 12 , W e dτ 12 dt + τ 12 = α cos(t) - W e(1 -a 2 ) cos(t) 2 ψ, ψ (0)
= ψ(2π) and τ 12 (0) = τ 12 (2π),

τ 11 = 1 + a 2 ψ and τ 22 = - 1 -a 2 ψ.
When a = ±1 the system admits an explicit solution:

ψ(t) = αW e (1 + W e 2 )(1 + 4W e 2 ) (1 + 2W e 2 ) cos(2t) + 3W e sin(2t) + 1 + 4W e 2 , (15)
τ 12 (t) = α 1 + W e 2 (cos(t) + W e sin(t)). (16)
Otherwise, when |a| < 1, the solution has no explicit form to our knowledge.

The solution for large Bi

The case Bi → +∞ is reached when the amplitude γ 0 becomes small. The material behaves as an elastic solid. Moreover, there exists a finite value Bi c > 0 such that when Bi > Bi c then the solution satisfies |τ d | < Bi. In that case, the system reduces to:

           W e dτ 11 dt -W e(1 + a) cos(t)τ 12 = 0, W e dτ 22 dt + W e(1 -a) cos(t)τ 12 = 0, W e dτ 12 dt + W e cos(t) 2 ((1 -a)τ 11 -(1 + a)τ 22 ) = α cos(t),
together with the periodic boundary condition τ (0) = τ (2π). Let ψ = τ 11τ 22 . Combining the previous equations, we get

                 W e dψ dt = 2W e cos(t)τ 12 , W e dτ 12 dt = α cos(t) - W e(1 -a 2 ) cos(t) 2 ψ, ψ(0)
= ψ(2π) and τ 12 (0) = τ 12 (2π),

τ 11 = 1 + a 2 ψ and τ 22 = - 1 -a 2 ψ.
When a = ±1 the solution is known explicitly:

ψ(t) = α(c -cos(2t)) 2W e , (17)
τ 12 (t) = α sin(t) W e , ( 18
)
where c is an arbitrary constant. In that case |τ d (t)| is maximum for t = π/2 when c ≥ 0 and t = -π/2 when c ≤ 0. We have: max

t∈[0,2π] |τ d (t)| = α 2 √ 2 W e ((1 + |c|) 2 + 16) 1/2 . Let Bi c = 17 8 α W e .
Then, when Bi ≥ Bi c , there exists a periodic solution τ (t) for all t in [0, 2π] such that |τ d (t)| ≤ Bi in [0, 2π] and this solution is expressed by ( 17)-( 18).

The solution for

Bi ∈ ]0, Bi c [ Problem ( 14) is approximated by using a centered second order scheme:

               W e 2∆t      τ n+1 11 -τ n-1 11 τ n+1 22 -τ n-1 22 τ n+1 12 -τ n-1 12      + F(t n , τ n ) =      0 0 α cos t n      , τ 0 = τ N , (19)
with t n = 2πn/N , n = 0 . . . N , ∆t = 2π/N and τ n ≈ τ (t n ). This is a N by N non-linear system of equations that we solve by using a non-linear Gauss-Seidel iterative algorithm. We use N = 1000 points in the [0, 2π] interval, and iterations are stopped when the maximal residual term in the discrete nonlinear problem [START_REF] Saramito | Efficient simulation of nonlinear viscoelastic fluid flows[END_REF] becomes lower than 10 -12 . Fig. 7 shows the solution τ 12 and τ 11 (while τ 22 = 0) for a = 1, Bi = 3, α = 1 ξ = 0, and W e = 0.1, 0.05 and 0.025. When a = 0 the components τ 12 and τ 11 are comparable up to a slight variation, while τ 22 = -τ 11 , and the corresponding solution is not represented.

Fourier analysis

We expand the solution σ 12 = (1α) γ + τ 12 as:

σ 12 (t) = k≥1 a k sin(kt) + b k cos(kt),
where the coefficients are expressed by: When the response is linear we have a k = b k = 0 for all k ≥ 2. There are two cases when this situation occurs: when Bi = 0 while the material behaves as a viscoelastic fluid, and when Bi → +∞ while its behaves as an elastic solid. Otherwise, only odd harmonics appear. In [START_REF] Hyun | Large amplitude oscillatory shear as a way to classify the complex fluids[END_REF], the authors propose large amplitude oscillatory shear as a way to classify complex fluids. Following this approach, we observe on Fig. 9 that the third harmonic contributes of about 20% of the first harmonic, and higher harmonics decay very fast. See also [25, p. 74] for a comparable observation, based on a comparison between experimental data on a linear hight-density polyethylene melt and several theoretical models. See also Fig. 11 in [START_REF] Hyun | Large amplitude oscillatory shear as a way to classify the complex fluids[END_REF] for a similar observation based on experimental data for a Xanthan gum solution. Therefore, the moduli obtained from the first harmonic via Fourier transformation analysis are not substantially different from the moduli calculated neglecting higher harmonics, as in the linear regime. By extension to the linear regimes, let us introduce the in-phase moduli G ′ and the out-of-phase moduli G ′′ : these are related to the first Fourier coefficients by G ′ = a 1 and Conversely, when Bi vanishes, the material behaves as a viscoelastic fluid: the constitutive equation is also linear and the solution is expressed by [START_REF] Puzrin | Rate-dependent hyperplasticity with internal functions[END_REF]. The corresponding values of G ′ and G ′′ are denoted by G ′ 0 and G ′′ 0 :

a k = 1 π
G ′′ = b 1 . When Bi ≥ Bi c ,
G ′ 0 = αW e 1 + W e 2 and G ′′ 0 = 1 - αW e 2 1 + W e 2
In the non-linear case, the moduli (G ′ , G ′′ ) are computed by using the solution (τ n ) 0≤n≤N of the approximate problem [START_REF] Saramito | Efficient simulation of nonlinear viscoelastic fluid flows[END_REF]: 0.8% 1.1% 1.6% 2.7% 4.5% 7% 12% 23%

a 2 k + b 2 k a 2 1 + b 2
G ′ ≈ 2 N N n=1 τ n 12 sin(t n ) and G ′′ ≈ 1 -α + 2 N N n=1 τ n 12 cos(t n ).
The viscosity ratio parameter α = η m /(η + η m ) controls G ′′ ∞ , the G ′′ moduli at high values of Bi (i.e. for small amplitudes). The choice α < 1 leads to a non-vanishing G ′′ ∞ , that has been found to be consistent with data measurements; the numerical computations are presented for α = 8/9. classification introduced in [8] when using large amplitude oscillatory shear. Many materials that present such a behavior are listed on page 61 of [START_REF] Hyun | Large amplitude oscillatory shear as a way to classify the complex fluids[END_REF]. See also Fig. 7.b of this reference for the (G ′ , G ′′ ) plot of a Xanthan gum 4% solution. The micro-structure of this class of materials is characterized by a network of links that develops some complex structures. When an external strain is imposed, the complex network structure is destroyed by large deformation over a critical strain. In [START_REF] Sim | Large amplitude oscillatory shear behavior of complex fluids investigated by a network model: a guideline for classification[END_REF], based on a network model, the authors explain the overshoot of G ′′ in term of the balance between the formation and the destruction of the network junctions.

Conclusion

A new model for elastoviscoplastic fluid flows that is objective and satisfies the second law of thermodynamics is proposed in [START_REF] Oldroyd | A rational formulation of the equations of plastic flow for a Bingham solid[END_REF]. A variant of the constitutive equation is also introduced in ( 12) in order to enforce finite extensional properties of the material. Large amplitude oscillatory shear has been performed and complex moduli (G ′ , G ′′ ) are qualitatively in good agreement with experimental data for many materials that present microscopic complex network structures with large rearrangements. The model is a good candidate for numerical simulation of elastoviscoplastic in multi-dimensional geometries: future works will perform such computations and compare them with experimental data measured on flows in complex geometries [START_REF] Dollet | Flow of a foam past an elliptical obstacle[END_REF].

By differentiation, as in [START_REF] Hyun | Large amplitude oscillatory shear as a way to classify the complex fluids[END_REF], we obtain the constitutive law (5) with σ = η ε + τ . The imposed traction σ(t) = σ > 0 leads to ε = (στ )/η and the constitutive equation ( 5) is then completed with the initial condition:

     λ 1 dτ dt (t) + 1 + η η m max 0, |τ (t)| -τ 0 |τ (t)| τ (t) = σ ∀t > 0, τ (0) = 0,
where λ 1 = λη/η m . Since τ (0) = 0 and τ (t) is continuous there exists t 0 > 0 such that when t ∈ [0, t 0 ] we have |τ (t)| < τ 0 : this is the first regime and the differential equation becomes linear:

λ 1 dτ dt (t) + τ (t) = σ, ∀t ∈]0, t 0 [. Then τ (t) = σ (1 -exp (-t/λ 1 )) and ε(t) = (λ 1 σ/η) (1 -exp (-t/λ 1 )) , 0 ≤ t ≤ t 0 .
When σ ≤ τ 0 , we have always τ (t) < τ 0 and thus the previous expressions of τ and ε are true for all t ≥ 0.

In that case, the elongation for a fixed imposed traction is bounded in time, which means that the material behaves as a solid. When σ > τ 0 , we have τ (t 0 ) = τ 0 for a finite time t 0 = λ 1 log (σ/(στ 0 )). When t > t 0 we have τ (t) > τ 0 : this is the second regime and the differential equation becomes: and σ 2 = (η m σ + ητ 0 )/(η + η m ). The solution in the second regime is:

τ (t) = σ 2 -(σ 2 -τ 0 ) exp (-(t -t 0 )/λ 2 ) , t ≥ t 0 , ε(t) = σ -σ 2 η (t -t 0 ) + λτ 0 η m + λ 2 (σ 2 -τ 0 ) η (1 -exp (-(t -t 0 )/λ 2 )) , t ≥ t 0 .
In that case, the elongation for a fixed imposed traction is unbounded in time, which means that such material behaves as a fluid. As a conclusion, the material behaves as an elastic solid when the traction is lower than the yield value and as a viscoelastic fluid otherwise (see Fig. For an imposed traction σ = σ and ε m (0) = 0 the problem is a nonlinear differential equation. When σ ≤ τ 0 the solution is constant and otherwise the problem becomes linear. The solution for all t > 0 is:

ε(t) =    σ/µ when σ ≤ τ 0 , σ/µ + σ -τ 0 µ m (1 -exp(-t/λ)) otherwise.
As a conclusion, the Puzrin-Houlsby material behaves as an elastic solid when the traction is lower than the yield value and as a viscoelastic solid otherwise (see Fig. Let σ = σ > 0 be the imposed traction and assume that ε m,i (0) = 0, i = 1, 2 as initial condition. The computation is comparable to the two previous examples and is not detailed here. When σ ≤ τ 0 then the solution is constant: ε(t) = σ/µ, for all t > 0. When σ > τ 0 , after reduction to some linear problems, we obtain the solution for all t > 0: στ 0 µ m .

ε(t) =
As a conclusion, the Isayev-Fan simplified material behaves as an elastic solid when the traction is lower than the yield value and as a viscoelastic fluid otherwise (see Fig. 2.d).

  ) and (3). Schwedoff proposed a rigid behavior ε = 0 when |τ | ≤ τ 0 and Oldroyd proposed a brutal change of model when reaching the yield value. Our proposition assures a continuous change from a solid to a fluid behavior of the material. The mechanical model is represented in Fig. 1.c. A friction τ 0 has been inserted in the Oldroyd viscoelastic model (Fig. 1.b). At stresses below the yield stress, the friction element remains rigid. The level of the elastic strain energy required to break the friction element is determined by the von Mises yielding criterion. Consequently, before yielding, the whole system predicts only recoverable Kelvin-Voigt viscoelastic deformation due to the spring and the viscous element η: the Kelvin-Voigt viscoelastic model is described by a spring and a viscous body in parallel (Fig. 1.a). The elastic behavior τ = µε is expressed in (5) in differential form where µ = η m /λ is the elasticity of the spring in Fig. 1.c. Before yielding, the total stress (a) Kelvin-Voigt (b) Oldroyd

Figure 2 .

 2 Figure 2. Some other elastoviscoplastic models.

3 W 4 Figure 4 .

 344 Figure 4. Dimensionless first normal stress difference ψ = τ 11 -τ 22 for uniaxial elongation when Bi = 1, a = 1 and α = 1: (a) influence of W e for ξ = 0; (b) influence of ξ for W e = 0.75.

τ 12 τ 11 -τ 22 t 2 Figure 5 .

 12112225 Figure 5. Simple shear flow for W e = 1 and α = 8/9: (a) a = 1 and Bi = 1; (b) a = 0 and Bi = 1; (c) a = 0 and Bi = 2.

10 1 Figure 6 .

 1016 Figure 6. Shear viscosity for a = 1, α = 1 and (a) ξ = 0; (b) ξ = 0.015.

2π 0 σπ 2π 0 σ

 00 [START_REF] Oldroyd | On the formulation of rheological equations of states[END_REF] (t) sin(kt) dt and b k = 1 12 (t) cos(kt) dt.

Figure 7 .

 7 Figure 7. Periodic shear flow: influence of W e = 0.1, 0.05, 0.025 : (top) τ 12 /Bi; (bottom) τ 11 /W e.

Figure 8 .

 8 Figure 8. Lissajous plots: influence of W e = 0.1, 0.05, 0.025 : (left) τ 12 /Bi versus γ(t) = sin(t) ; (right) τ 12 /Bi versus γ(t) = cos(t).

Figure 9 .

 9 Figure 9. Periodic shear flow: intensities of the harmonics normalized by the intensity of the fundamental frequency: W e = 1/40, Bi = 3, a = 1 and α = 8/9.

Fig. 10 plots 3 Figure 10 .

 10310 Fig.10plots G ′ and G ′′ versus the dimensionless number Bi -1 = γ 0 η 0 ω/τ 0 . This is a dimensionless version of the moduli plots since most figures in the literature are based on the representation of the moduli versus the amplitude γ 0 . Only G ′′ shows an overshoot: this behavior matches the type III of the

  + τ (t) = σ 2 , ∀t > t 0 together with the initial condition τ (t 0 ) = τ 0 , where λ 2 = λη/(η + η m )

B. 2 .

 2 1.d). The Puzrin-Houlsby modelThe potentials associated to the one-dimensional version of model represented on Fig.2.a are:E (ε, ε m ) = µ m 2 |ε m | 2 + µ 2 |εε m | 2 and D ( ε, εm ) = f m ( εm ) , where f m (d) = τ 0 |d| + η m 2 |d| 2 , ∀d ∈ R.By differentiation, as in[START_REF] Hyun | Large amplitude oscillatory shear as a way to classify the complex fluids[END_REF], we obtain the consti-tutive laws: σ = µ(εε m ) and µ(εε m )µ m ε m ∈ ∂f m ( εm ). Since for all ξ ∈ R we have ∂f * m (d) = {(1/η m ) max(0, 1τ 0 /|d|) d}, ∀d ∈ R,these relations can be rewritten as: ε = ε m + σ/µ and εm = 1 η m max(0, 1τ 0 /|σµ m ε m |) (σµ m ε m ).

B. 3 . 2 |d| 2 ,

 322 2.b). The one-mode simplified Isayev-Fan modelThe potentials associated to the one-dimensional version of model[START_REF] Isayev | Viscoelastic plastic constitutive equation for flow of particle filled polymers[END_REF] represented on Fig.2.c are:E(ε, ε m,1 , ε m,2 ) = µ m 2 |ε m,1 | 2 + µ 2 |εε m,1ε m,2 | 2 , D( ε, εm,1 , εm,2 ) = f m ( εm,1 + εm,2 ) + η m,2 2 | εm,2 | 2 ,where f m (d) = τ 0 |d| + η m ∀d ∈ R. The constitutive laws are obtained as usual by differentiation.

  when σ ≤ τ 0 , σ/µ + bt + (cbλ)(1exp(-t/λ)) otherwise, where λ = η m η m,2 η m + η m,2 µ m , b = στ 0 η m + η m,2 and c = η m η m,2 η m + η m,2

  Figure 3. Summary table of the referenced contributions on elastoviscoplastic models.

	year contribution	before yielding	after yielding	3D TH
	1900 Schwedoff [20]	rigid solid	Maxwell viscoelastic fluid
	1922 Bingham [1]	rigid solid	Newtonian fluid	X X
	1926 Herschel-Bulkley [6]	rigid solid	power-law fluid	X X
	1947 Oldroyd [11]	elastic solid	Newtonian fluid	X X
	1990 Isayev-Fan [9]	elastic solid	Leonov viscoelastic fluid X X
	1991 Beris et al. [3]	elastic solid	power-law fluid	X
	2003 Puzrin-Houlsby [15, p. 254] elastic solid	Kelvin viscoelastic solid X X
	2007 present	Kelvin-Voigt viscoelastic solid Oldroyd viscoelastic fluid X X

Appendix A. Subdifferential calculus in the multi-dimensional case A.1. The ϕ m function -The subgradient ∂ϕ m , as introduced in [START_REF] Houlsby | Rate-dependent plasticity models derived from potential functions[END_REF], is defined for any tensor D by:

where τ d denotes the deviatoric part of τ .