R Énald Chalayer

Laurent Chupin

Thierry Dubois

A BI-PROJECTION METHOD FOR INCOMPRESSIBLE BINGHAM FLOWS WITH VARIABLE DENSITY, VISCOSITY AND YIELD STRESS

Keywords: May 30, 2018 Viscoplastic medium, Variable density, Bingham flows, Projection method, Fractional time-stepping, Stability and error analysis. 1

A bi-projection method for incompressible Bingham flows with variable density, viscosity and yield stress

Rénald Chalayer, Laurent Chupin, Thierry

Introduction

Viscoplastic materials are common in many industrial processes, as in food industry with dairy products, chocolate confection, pulp suspensions, in the petroleum industry with drilling mud, cement, waxy crude oil, and in various geophysical phenomena such flows of slurries, debris and lava. In many cases, mixtures of fluids with different properties occur as sliding of mud in water, dam break of granular matter in air or water. In food industry, the cleaning of a viscoplastic material from conduits is achieved through the displacement induced by water steam or air. Also in the process of oil recovery, a flow of gas, steam or even water is used to push and transport heavy and highly viscous fluids. Viscoplastic fluids flow only if the stress exceeds a threshold, namely the yield stress, otherwise they do not deform and behave like solids. Typical viscoplastic fluid flows exhibit both yielded and unyielded regions. The rigid structures, where the deformation rates vanish, may either be plugs moving as solid bodies or frozen zones, i.e. where the flow is at rest (no velocity). The most commonly used yield stress model to account for this peculiar behaviour is the Bingham model. In the constitutive law, the stress tensor is not prescribed until the yield stress is reached while, above the threshold, it is proportional to the strain-rate tensor. This singular change of rheological properties induces mathematical difficulties. From a computational point of view, reproducing accurately the yield surfaces separating yielded and unyielded regions is challenging (see [START_REF] Saramito | Progress in numerical simulation of yield stress fluid flows[END_REF] for a recent review on numerical simulation of viscoplastic flows). Two approaches have been mainly used both in the fields of theoretical studies and numerical simulations, namely regularization methods (see for instance [START_REF] Bercovier | A finite element method for incompressible non-Newtonian flows[END_REF] and [START_REF] Papanastasiou | Flows of materials with yield[END_REF]) and optimization techniques relying on the theory of variational inequalities (see [START_REF] Duvaut | Les inéquations en mécanique et en physique[END_REF] and [START_REF] Glowinsky | Numerical Analysis of Variational Inequalities[END_REF]). Note that theoretical analyses concerning the questions of existence, uniqueness and regularity of solutions have been investigated in [START_REF] Duvaut | Les inéquations en mécanique et en physique[END_REF] and more recently in [START_REF] Fuchs | Regularity results for the quasi-static Bingham variational inequality in dimensions two and three[END_REF][START_REF] Fuchs | Variational methods for problems from plasticity theory and for generalized Newtonian fluids[END_REF]. Due to their simplicity regularization methods are appealing : the singular Bingham law is replaced by a regularized form making the flow viscous everywhere with a large but finite viscosity in the unyielded regions. The main drawback of this approach is the difficulty to accurately locate the yield surfaces which are not accounted for by the regularized model. A review for this approach can be found in [START_REF] Frigaard | On the usage of viscosity regularisation methods for visco-plastic fluid flow computation[END_REF]. The second approach used to overcome the difficulty due to the non-differentiability of the definition of the stress tensor is based on the variational formulation of the Navier-Stokes equations which leads to a saddle-point problem. This optimization problem can be solved by using the Uzawa-like method or the augmented Lagrangian method (see [START_REF] Dean | On the numerical simulation of Bingham viscoplastic flow: old and new results[END_REF] and the references therein). Both algorithms allow for accurate but rather expansive numerical simulations when a fine description of the plug regions and the yield surfaces is required (see [START_REF] Muravleva | Uzawa-like methods for numerical modeling of unsteady viscoplastic Bingham medium flows[END_REF][START_REF] Marly | Augmented Lagrangian simulations study of yield-stress fluid flows in expansioncontraction and comparisons with physical experiments[END_REF]). Efforts have been recently done to accelerate the convergence of the Augmented Lagrangian algorithm (see [START_REF] Treskatis | An accelerated dual proximal gradient method for applications in viscoplasticity[END_REF][START_REF] Saramito | A damped Newton algorithm for computing viscoplastic fluid flows[END_REF]). In [START_REF] Saramito | A damped Newton algorithm for computing viscoplastic fluid flows[END_REF], a modified Newton method achieving superlinear convergence is proposed. The Uzawa-like method (see [START_REF] Dean | On the numerical simulation of Bingham viscoplastic flow: old and new results[END_REF][START_REF] Muravleva | Uzawa-like methods for numerical modeling of unsteady viscoplastic Bingham medium flows[END_REF]) relies on a (pointwise) projection operator for the computation of the extra (plastic part) stress tensor. In [START_REF] Dean | On the numerical simulation of Bingham viscoplastic flow: old and new results[END_REF], the authors suggested to add a pseudo-time relaxation term in the Bingham projection operator in order to increase the convergence rate of the fixed-point algorithm used to compute the plastic tensor. In [START_REF] Chupin | A bi-projection method for Bingham type flows[END_REF], this approach was coupled with a projection scheme [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF][START_REF] Temam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires[END_REF] for the time discretization of the Navier-Stokes equations. The resulting bi-projection scheme was analyzed in terms of stability study and error estimates. More precisely, it was shown in [START_REF] Chupin | A bi-projection method for Bingham type flows[END_REF] that the pseudo-relaxation term in the Bingham projection does not deteriorate the first-order accuracy of the time discretization scheme. Indeed, the error is of the order δt(1 + θ) which is of first order if the prescribed relaxation parameter θ is of the order of the time step δt. Through various numerical simulations in standard configurations, the bi-projection scheme was shown in [START_REF] Chupin | A bi-projection method for Bingham type flows[END_REF] to be able to reproduce the characteristic property of Bingham flows to return to rest in finite time, to accurately predict the plug regions, and to be efficient at both large Bingham and high Reynolds numbers. The main contribution of the present paper is to extend the numerical scheme proposed in [START_REF] Chupin | A bi-projection method for Bingham type flows[END_REF] to incompressible Bingham flows with variable (in space and time) density, plastic viscosity and yield stress. The key feature of the projection method for the time discretization of the Navier-Stokes equations is to decouple the computations of the velocity field and the pressure. In a first step, a predicted velocity which is not divergence free is computed from which the gradient of a scalar function (a pseudo-pressure) is then subtracted resulting in a solenoidal velocity field. The Helmholtz decomposition is invoked in the second step so that the (pseudo-)pressure is solution of a Poisson equation. While introduced in the pioneering works of Chorin [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF] and Temam [START_REF] Temam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires[END_REF] at the end of the 60's in the case of homogeneous fluid flows, i.e. with constant density and viscosity, projection methods have first been analyzed in terms of error estimates by Shen in [START_REF] Shen | On error estimates of the projection methods for the Navier-Stokes equations: first-order schemes[END_REF] in the early 90's. Let us also mention [START_REF] Guermond | On the approximation of the unsteady Navier-Stokes equations by finite element projection methods[END_REF] for an analysis of a finite element fully discrete version. Fractional step schemes in the spirit of projection methods have been later on extended to incompressible Navier-Stokes equations with variable density and used in numerical simulations, see for example [START_REF] Bell | A second-order projection method for variable-density flows[END_REF][START_REF] Almgren | Conservative adaptive projection method for the variable density incompressible Navier-Stokes equations[END_REF]. Concerning the mathematical analysis of projection schemes for variable-density flows, stability results have been proved in [START_REF] Guermond | A projection FEM for variable density incompressible flows[END_REF][START_REF] Pyo | Gauge-Uzawa methods for incompressible flows with variable density[END_REF] and error estimates have been derived in [START_REF] Guermond | Error analysis of a fractional time-stepping technique for incompressible flows with variable density[END_REF] for a fully discrete version based on the finite element method for the spatial discretization. For density-variable flows, if the projection step is achieved as in the homogeneous case, i.e. by using the Helmholtz decomposition, a variable coefficient elliptic equation has to be solved at each time step (see [START_REF] Guermond | A splitting method for incompressible flows with variable density based on a pressure Poisson equation[END_REF] for instance) in order to compute the pressure, inducing both mathematical and computational difficulties. In [START_REF] Guermond | A splitting method for incompressible flows with variable density based on a pressure Poisson equation[END_REF][START_REF] Guermond | Error analysis of a fractional time-stepping technique for incompressible flows with variable density[END_REF] an alternative approach based on the interpretation of projection methods as penalty methods is used. As a result, the fractional time-stepping methods proposed in [START_REF] Guermond | A splitting method for incompressible flows with variable density based on a pressure Poisson equation[END_REF][START_REF] Guermond | Error analysis of a fractional time-stepping technique for incompressible flows with variable density[END_REF] necessitate the resolution of only one Poisson equation per time step and are therefore more efficient from a computational point of view. The aim of this paper is to make use of these recent developments to propose a new time discretization scheme for solving the equations describing the motion of Bingham fluids with variable density, plastic viscosity and yield stress. Unlike in [START_REF] Guermond | A splitting method for incompressible flows with variable density based on a pressure Poisson equation[END_REF][START_REF] Guermond | Error analysis of a fractional time-stepping technique for incompressible flows with variable density[END_REF], we consider the temporal semi-discrete equations. Also, as a main difference, the divergence-free velocity field is used as convective velocity in the discrete version of the mass conservation equation allowing to derive lower and upper bounds as well as error estimates on the density. Concerning the treatment of the plastic tensor, a fixed-point algorithm is invoked to solve the Bingham projection, in the spirit of [START_REF] Muravleva | Uzawa-like methods for numerical modeling of unsteady viscoplastic Bingham medium flows[END_REF]. As in [START_REF] Chupin | A bi-projection method for Bingham type flows[END_REF], a pseudo-time relaxation term is added in the computation of the plastic tensor through the Bingham projection in order to provide a geometric convergence of this algorithm. The objective of this paper is to perform stability study and error analyses of the proposed bi-projection scheme. To the best of our knowledge, no stability and convergence analysis of projection schemes for incompressible variable density Bingham flows have been done yet. As in the case of homogeneous Bingham flows, the bi-projection scheme is proved to be stable and first order accurate as long as the relaxation parameter is taken of the order of the time step. The paper is organized as follows. In Section 2, the mathematical formulation for a Bingham model with variable coefficients is introduced. A projection formulation for the plastic part of the stress tensor which will be suitable for the construction of our scheme is provided. In Section 3, some additional notations are introduced and preliminary results, that will be useful in the following sections, are given. In Section 4, the bi-projection scheme, as a time approximation of the continuous model, is shown to be well-posed and bounds on the density, and the dependent coefficients, that are the viscosity and yield stress, are obtained. Sections 5 and 6 are respectively devoted to the stability and error (convergence) analyses of the scheme.

The model of a Bingham viscoplastic flow with variable density

2.1. The mathematical model. Let T > 0 be a positive real number and Ω a bounded domain of R 3 ; we denote by Γ the boundary of Ω. The isothermal flow of an incompressible viscoplastic medium with variable density is modeled by the following equations, satisfied by the velocity u, the density ρ and the pressure p,

     ∂ t ρ + u • ∇ρ = 0, ρ ∂ t u + u • ∇u + ∇p = div τ , div u = 0, in (0, T ) × Ω. (1)
The deviatoric stress tensor τ is defined by the relation

(2) τ = 2µ(ρ)Du + α(ρ)Σ,
with Du = 1 2 ∇u + T ∇u the strain-rate tensor and Σ is the extra (plastic) part of the stress tensor. The plastic viscosity µ and the yield stress α are functions of the density ρ of class C 1 and are assumed to be respectively positive and non-negative. We denote by

Q 0 = {(t, x) ∈ (0, T ) × Ω; Du(t, x) = 0}
the sub-domain of (0, T ) × Ω where the strain-rate tensor vanishes. The plastic tensor Σ ∈ R 3×3 is defined by (see [START_REF] Duvaut | Les inéquations en mécanique et en physique[END_REF] for instance)

(3)    Σ(t, x) = Du(t, x) |Du(t, x)| for (t, x) ∈ (0, T ) × Ω \ Q 0 , |Σ(t, x)| ≤ 1, T Σ(t, x) = Σ(t, x), tr Σ(t, x) = 0 for (t, x) ∈ Q 0 ,
where, for all λ ∈ R 3×3 , we denote by ( 4)

|λ| 2 = 1 2 tr T λλ
its second invariant (see [START_REF] Balmforth | Yielding to stress: recent developments in viscoplastic fluid mechanics[END_REF]). The system (1)-( 2)-( 3) is supplemented with the following initial and boundary conditions ( 5)

ρ| t=0 = ρ 0 , u| t=0 = u 0 , u| Γ = 0.
Finally, we assume there exists ρ 1 , ρ 2 > 0,

ρ 1 ≤ ρ 0 ≤ ρ 2 , a.e. (6)
in Ω, and u 0 ∈ H 1 0 (Ω) with div u 0 = 0. If the velocity field u solution of (1) is sufficiently regular (see [START_REF] Boyer | Éléments d'analyse pour l'étude de quelques modèles d'écoulements de fluides visqueux incompressibles[END_REF]), as ρ is solution of a transport equation, u is divergence free and the initial density satisfies (6), we have [START_REF] Chupin | A bi-projection method for Bingham type flows[END_REF] ρ 1 ≤ ρ(t, x) ≤ ρ 2 a.e. in (0, T ) × Ω.

This ensures that, almost everywhere in (0, T ) × Ω, we have

(8) ∃ µ 1 , µ 2 > 0, µ 1 ≤ µ(ρ(t, x)) ≤ µ 2 , ∃ α 1 , α 2 ≥ 0, α 1 ≤ α(ρ(t, x)) ≤ α 2 .
In the Newtonian case, such regularity results can be obtained if the data are small enough (see for instance [START_REF] Huang | Global strong solution of 3D inhomogeneous Navier-Stokes equations with densitydependent viscosity[END_REF] and the references therein). For yield stress fluids, we assume that similar regularity results can be obtained. In the rest of the paper, we assume that solutions of (1) satisfying [START_REF] Chupin | A bi-projection method for Bingham type flows[END_REF] exists so that (8) is satisfied.

2.2.

A projection formulation. We now introduce a projection formulation for the plastic tensor Σ, which is more suited for the construction of our numerical scheme.

Proposition 1. For all ℓ > 0, the condition (3) is equivalent to the relation

(9) Σ(t, x) = P Σ(t, x) + ℓDu(t, x) ,
where P is the projection operator on the closed convex set defined by

Λ = {λ ∈ R 3×3 ; |λ| ≤ 1, T λ = λ, tr λ = 0}.
Proof. If Du(t, x) = 0, the equivalence is obvious as we have in both cases Σ(t, x) ∈ Λ. We now assume that Du(t, x) = 0. An explicit expression of the projection is given by ( 10)

P(λ) = λ if |λ| ≤ 1, λ |λ| if |λ| > 1.
Let ℓ > 0 be a positive real number. If (3) is satisfied, we have

Σ(t, x) = Du(t, x) |Du(t, x)| . Therefore, Σ(t, x) + ℓDu(t, x) = 1 + ℓ Du(t, x) > 1,
so that, according to (10), we have

P Σ(t, x) + ℓDu(t, x) = Σ(t, x) + ℓDu(t, x) |Σ(t, x) + ℓDu(t, x)| = Du(t, x) |Du(t, x)| = Σ(t, x),
hence, (9) follows.

Let us now assume that ( 9) is fulfilled. We have |Σ(t, x) + ℓDu(t, x)| > 1 (otherwise, ℓ = 0 due to [START_REF] Evans | On solving certain nonlinear partial differential equations by accretive operator methods[END_REF]). Therefore, according to [START_REF] Evans | On solving certain nonlinear partial differential equations by accretive operator methods[END_REF] we have

Σ(t, x) = Σ(t, x) + ℓDu(t, x) |Σ(t, x) + ℓDu(t, x)| , so that |Σ(t, x)| = 1.
We also deduce that

Σ(t, x) |Σ(t, x) + ℓDu(t, x)| -1 = ℓDu(t, x), so that |Du(t, x)| = |Σ(t, x) + ℓDu(t, x)| -1 ℓ . Then Σ(t, x) = Du(t, x) |Du(t, x)| which concludes the proof.
Remark 1. When Du is only Lebesgue-integrable in space, we define Σ as above with a representative of Du. Then, we define Σ as the class of Σ.

Introducing Λ = {f ∈ L 2 (Ω) 3×3 ; |f | ≤ 1, T f = f , tr f = 0, a.e. in Ω}, we have Σ(t) ∈ Λ.
With the help of Proposition 1, we can rewrite the mathematical model for incompressible visco-plastic flows with variable density, viscosity and yield stress as, for any ℓ > 0,

     ∂ t ρ + u • ∇ρ = 0, ρ ∂ t u + u • ∇u + ∇p -2 div µ(ρ)Du = div α(ρ)Σ , Σ = P(Σ + ℓDu), div u = 0, (11)
which is supplemented with the initial and boundary conditions [START_REF] Boyer | Éléments d'analyse pour l'étude de quelques modèles d'écoulements de fluides visqueux incompressibles[END_REF]. Note that for constant density flows, the system [START_REF] Frigaard | On the usage of viscosity regularisation methods for visco-plastic fluid flow computation[END_REF] reduces to the mathematical model studied in [START_REF] Chupin | A bi-projection method for Bingham type flows[END_REF].

Preliminaries

We now introduce some of the notations used in the sequel. For two vectors u and

v in R 3 , we denote by u • v their inner product, u • v = 1≤i≤3 u i v i ,
A : B = 1≤i,j≤3 A ij B ij .
Let us note that the tensorial norm | • | defined by ( 4) is not induced by the above inner product. We will make use of the standard notations L p (Ω), H k (Ω) and H k 0 (Ω) to denote the usual Lebesque and Sobolev spaces over Ω. We denote by L 2 0 (Ω) the subspace of L 2 (Ω) of functions with vanishing mean value. The norm corresponding to H k (Ω) will be denoted by • k . In particular, we will use • to denote the norm in L 2 (Ω) and (•, •) to denote the scalar product in L 2 (Ω). For each space above, we will use bold letters to denote their vectorial counterpart: for example, L 2 (Ω)

3 will be denoted L 2 (Ω). Also other functional spaces may be considered in the sequel. In such cases, they are always indicated as subscripts: for instance • L ∞ (Ω) denotes the norm associated to the space L ∞ (Ω).

The incompressibility constraint leads us to consider the following space

H = {u ∈ L 2 (Ω); div u = 0, u • n| Γ = 0},
where n is the unit outward normal to Γ. For any sequence (a n ) n∈N , we introduce

δa n = a n -a n-1 and δ 2 a n = δ(δa n ) = a n -2a n-1 + a n-2 .
The following lemma of Gronwall type will be used to derive the error estimates in Section 6.

Note that a proof can be found in [START_REF] Heywood | Finite-element approximation of the nonstationary Navier-Stokes problem. IV. Error analysis for second-order time discretization[END_REF].

Lemma 1. (Discrete Gronwall lemma). Let M ∈ N * , τ , B and C be non-negative parameters, (y n ) n , (h n ) n , (g n ) n and (f n ) n be non-negative sequences satisfying, for all m such that 0 ≤ m ≤ M , y m + τ m n=0 h n ≤ B + τ m n=0 (g n y n + f n ), with τ M n=0 g n ≤ C.
Assume τ g n < 1 for all n such that 0 ≤ n ≤ M and let

σ n = (1 -τ g n ) -1 , σ = max 0≤n≤M σ n ,
then, for all m such that 0 ≤ m ≤ M , we have

y m + τ m n=0 h n ≤ exp(σC) B + τ m n=0 f n .
Finally, we define the bilinear form B by

B(u, v) = (u • ∇)v + v 2 div u,
and the associated trilinear form: b(u, v, w) = ( B(u, v), w). It can be easily shown that

(12) b(u, v, v) = 0, ∀u ∈ H 1 (Ω) with u • n| Γ = 0, ∀v ∈ H 1 (Ω) .
Another orthogonality property will be also useful in the sequel, namely we have ( 13)

Ω ϕ u • ∇ϕ dx = 0, ∀u ∈ L 3 (Ω) ∩ H, ∀ϕ ∈ H 1 (Ω).
4. The bi-projection scheme 4.1. The temporal discretization. In this section, we propose a scheme to discretize with respect to the time variable the model [START_REF] Frigaard | On the usage of viscosity regularisation methods for visco-plastic fluid flow computation[END_REF]. Let us introduce some additional notations. Let N > 0 be an integer and (t n ) n∈{0,...,N } a sequence of discrete time levels in [0, T ]. For the sake of simplicity, we consider a uniform discretization, that is

t n = nδt with δt = T N .
We introduce two additional numerical parameters r > 0 and θ > 0.

We start with ρ 0 , u 0 , û0 , q 0 , p 0 and Σ 0 (initialization step). We choose ρ 0 = ρ 0 , u 0 = û0 = u 0 , q 0 = 0, and arbitrary Σ 0 and p 0 . For n ≥ 0, assuming that ρ n , u n , ûn , q n , p n and Σ n are known, we first compute ρ n+1 by solving

(14) ρ n+1 -ρ n δt + ûn • ∇ρ n+1 = 0,
and, with the help of ( 8), we define µ n+1 and α n+1 by the following relations ( 15)

µ n+1 = µ(ρ n+1 ), α n+1 = α(ρ n+1 ).
Next, u n+1 and Σ n+1 are solutions of the system ( 16)

               1 δt 1 2 (ρ n+1 + ρ n )u n+1 -ρ n u n + B(ρ n+1 ûn , u n+1 ) -div(2µ n+1 Du n+1 ) + ∇(p n + q n ) = div(α n+1 Σ n+1 ), Σ n+1 = P Σ n+1 + rα n+1 Du n+1 + θ(Σ n -Σ n+1 ) , u n+1 | Γ = 0.
Finally, we compute q n+1 by solving the Poisson equation ( 17)

∆q n+1 = ρ 1 δt div(u n+1 ), ∂ n q n+1 | Γ = 0,
and we write

p n+1 = p n + q n+1 , (18)
ûn+1 = u n+1 - δt ρ 1 ∇q n+1 . ( 19
)
Remark 2. Using (17), ( 18) and [START_REF] Guermond | Error analysis of a fractional time-stepping technique for incompressible flows with variable density[END_REF], we observe that, for all 0 ≤ n ≤ N, ûn satisfies ûn • n| Γ = 0 and div ûn = 0.

Note that a fractional time-stepping method in the spirit of [START_REF] Guermond | A splitting method for incompressible flows with variable density based on a pressure Poisson equation[END_REF], [START_REF] Guermond | Error analysis of a fractional time-stepping technique for incompressible flows with variable density[END_REF] is employed to decouple the velocity and the pressure. A major difference between [START_REF] Guermond | A splitting method for incompressible flows with variable density based on a pressure Poisson equation[END_REF], [START_REF] Guermond | Error analysis of a fractional time-stepping technique for incompressible flows with variable density[END_REF] and the above algorithm ( 14)- [START_REF] Guermond | Error analysis of a fractional time-stepping technique for incompressible flows with variable density[END_REF] resides in the treatment of the mass conservation equation, namely the divergencefree velocity ûn is used as the convective velocity in the equation of the density [START_REF] Glowinsky | Numerical Analysis of Variational Inequalities[END_REF]. With this approach, lower and upper bounds on the density (ρ n ) n≥0 can be derived. More precisely, we will prove that, for all 0

≤ n ≤ N , ρ 1 ≤ ρ n ≤ ρ 2 .
In order to resolve the coupling in ( 16) between the velocity u n+1 and the plastic tensor Σ n+1 , we advocate a fixed-point algorithm, in the spirit of [START_REF] Muravleva | Uzawa-like methods for numerical modeling of unsteady viscoplastic Bingham medium flows[END_REF], which is detailed in Section 4.3. As in [START_REF] Chupin | A bi-projection method for Bingham type flows[END_REF], a pseudo-time relaxation term θ(Σ n -Σ n+1 ) is added in the projection operator defining Σ n+1 in ( 16). This additional term guarantees a geometric convergence of this algorithm. System ( 17) is a classical Poisson equation. As an immediate consequence, we have q n+1 ∈ H 2 (Ω) and so ûn+1 ∈ H 1 (Ω) exist and are uniquely defined. Finally, let us mention that the coefficient rα n+1 is used, instead of r only, in the Bingham projection defining Σ n+1 . We will see in Section 4.3 that the use of rα n+1 is required to derive proper estimates and convergence of the fixed-point algorithm.

4.2.

Discrete maximum principle on the density. The use of the divergence free velocity in the mass conservation equation ( 14) allows us to derive positive upper and lower bounds on the density. This property will be repeatedly used in the rest of the paper. Let us first establish this result. Let ǫ > 0 be a parameter. In order to properly define solutions of ( 14), we introduce the following regularized form (20)

   ρ n+1 ǫ -ρ n δt + ûn • ∇ρ n+1 ǫ -ǫ∆ρ n+1 ǫ = 0, ∇ρ n+1 ǫ • n| Γ = 0.
The limit of ρ n+1 ǫ when ǫ tends to 0 is called the viscosity solution of ( 14) (see [START_REF] Evans | On solving certain nonlinear partial differential equations by accretive operator methods[END_REF] and [START_REF] Grandall | Viscosity solutions of Hamilton-Jacobi equations[END_REF]). In the remaining part of the paper, ρ n+1 will always be used to denote the viscosity solution of [START_REF] Glowinsky | Numerical Analysis of Variational Inequalities[END_REF].

Proposition 2. If ûn ∈ H ∩ L 3 (Ω) and ρ 1 ≤ ρ n ≤ ρ 2 almost
Find ρ n+1 ǫ ∈ H 1 (Ω) so that, for all ϕ ∈ H 1 (Ω), we have

Ω ρ n+1 ǫ ϕ -δt Ω ρ n+1 ǫ ûn • ∇ϕ + δtǫ Ω ∇ρ n+1 ǫ • ∇ϕ = Ω ρ n ϕ.
By virtue of the Lax-Milgram theorem (using ûn ∈ H ∩ L 3 (Ω) and the property ( 13)), there exists a unique ρ n+1 ǫ ∈ H 1 (Ω) solution of [START_REF] Huang | Global strong solution of 3D inhomogeneous Navier-Stokes equations with densitydependent viscosity[END_REF], which satisfies [START_REF] Heywood | Finite-element approximation of the nonstationary Navier-Stokes problem. IV. Error analysis for second-order time discretization[END_REF] almost everywhere in Ω and Γ.

Step 2: Bound of the regularized solution. We want to prove that ρ 1 ≤ ρ n+1 ǫ ≤ ρ 2 almost everywhere in Ω. By choosing ϕ = (ρ n+1 ǫ -ρ 2 ) + in [START_REF] Huang | Global strong solution of 3D inhomogeneous Navier-Stokes equations with densitydependent viscosity[END_REF], where the + superscript denotes the positive part (similarly a -superscript denotes the negative part), using again [START_REF] Fuchs | Variational methods for problems from plasticity theory and for generalized Newtonian fluids[END_REF], and defining the sign function by sign(x) = 1 if x ≥ 0 and sign(x) = -1 otherwise, we obtain

Ω (ρ n+1 ǫ -ρ 2 )(ρ n+1 ǫ -ρ 2 ) + + δtǫ 2 Ω |∇ρ n+1 ǫ | 2 1+ sign(ρ n+1 ǫ -ρ 2 ) = Ω (ρ n -ρ 2 )(ρ n+1 ǫ -ρ 2 ) + .
Observing that the right-hand side is negative and the second term of the left-hand side is non-negative, we have

Ω (ρ n+1 ǫ -ρ 2 )(ρ n+1 ǫ -ρ 2 ) + ≤ 0.
On the other hand, xx + is non-negative for any x ∈ R. As a consequence, 22)

(ρ n+1 ǫ -ρ 2 )(ρ n+1 ǫ -ρ 2 ) + = 0 a.e.
ρ 1 ≤ ρ n+1 ε ≤ ρ 2 , a.e. in Ω.
Step 3: Passage to the limit on ǫ and viscosity solution. By choosing ϕ = ρ n+1 ǫ in (21), we have

Ω |ρ n+1 ǫ | 2 + δtǫ Ω |∇ρ n+1 ǫ | 2 = Ω ρ n ρ n+1 ǫ .
By invoking the Cauchy-Schwarz inequality, it follows

ρ n+1 ǫ 2 + δt √ ǫ∇ρ n+1 ǫ 2 ≤ ρ n ρ n+1 ǫ .
Then, with the help of Young's inequality, we obtain

ρ n+1 ǫ 2 + 2δt √ ǫ∇ρ n+1 ǫ 2 ≤ ρ n 2 .
By using this last inequality and ( 22), we deduce that (ρ n+1 ǫ ) ǫ , and ( √ ǫ∇ρ n+1 ǫ ) ǫ , are bounded sequences in respectively L ∞ (Ω) and L 2 (Ω). Hence, there exists weakly convergent subsequences

(23) ρ n+1 ǫ * ⇀ ǫ→0 ρ n+1 in L ∞ (Ω), √ ǫ∇ρ n+1 ǫ ⇀ ǫ→0 ψ in L 2 (Ω) .
As a consequence of the weak convergence and by using [START_REF] Lewicka | On the optimal constants in Korn's and geometric rigidity estimates, in bounded and unbounded domains, under Neumann boundary conditions[END_REF],

ρ 1 ≤ ρ n+1 ≤ ρ 2 , a.e. in Ω.
Finally, for all ϕ ∈ D(Ω), we have

Ω ρ n+1 ǫ ϕ -→ ǫ→0 Ω ρ n+1 ϕ, -δt Ω ρ n+1 ǫ ûn • ∇ϕ -→ ǫ→0 -δt Ω ρ n+1 ûn • ∇ϕ, δtǫ Ω ∇ρ n+1 ǫ • ∇ϕ -→ ǫ→0 0, so that Ω ρ n+1 ϕ -δt Ω ρ n+1 ûn • ∇ϕ = Ω ρ n ϕ.
We conclude that ρ n+1 satisfies the transport equation ( 14) in the distributions sense.

We are now able to establish the following result:

Theorem 1. If u 0 and ρ 0 satisfy (6), the sequences (û n , ρ n ) 0≤n≤N solutions of ( 14)-( 19) satisfy, for all n such that 0

≤ n ≤ N , ûn ∈ H ∩ H 1 (Ω), and
ρ 1 ≤ ρ n ≤ ρ 2 , a.e. in Ω.
Proof. We proceed by induction.

Case n = 0. From (6), using ρ 0 = ρ 0 and û0 = u 0 , we infer

ρ 1 ≤ ρ 0 ≤ ρ 2 , a.e.
in Ω, and û0 ∈ H ∩ H 1 (Ω) .

Case n ≥ 1. We assume

ρ 1 ≤ ρ n ≤ ρ 2 , a.e.
in Ω, and ûn ∈ H ∩ H 1 (Ω) .

From Proposition 2, we have ρ 1 ≤ ρ n+1 ≤ ρ 2 almost everywhere. The existence and uniqueness for the solution (u n+1 , Σ n+1 ) ∈ H 1 0 (Ω) ×Λ to the system ( 16) is a consequence of the theory of variational inequalities. The precise result for this system can be found in [START_REF] Dean | On the numerical simulation of Bingham viscoplastic flow: old and new results[END_REF] (page 47). Then, by solving [START_REF] Guermond | A projection FEM for variable density incompressible flows[END_REF], we obtain the existence and uniqueness of q n+1 ∈ H 2 (Ω). As a consequence, ûn+1 given by [START_REF] Guermond | Error analysis of a fractional time-stepping technique for incompressible flows with variable density[END_REF] lies in H ∩ H 1 (Ω), which concludes the proof. Remark 3. Using (8) and the above result, we observe that, for all 0 ≤ n ≤ N ,

µ 1 ≤ µ n ≤ µ 2 , α 1 ≤ α n ≤ α 2 , a.e. in Ω.

4.3.

Practical implementation of the Bingham projection. The equations ( 16) involve a coupling between u n+1 and Σ n+1 . In order to solve this system in practice and following [START_REF] Chupin | A bi-projection method for Bingham type flows[END_REF], we employ a fixed-point iteration procedure and proceed as it follows. Let us start with Σ n,0 = Σ n . For k ≥ 0 we assume that Σ n,k ∈ Λ is known and we compute u n,k ∈ H 1 0 (Ω) by solving the following elliptic problem ( 24)

           1 δt 1 2 (ρ n+1 + ρ n )u n,k -ρ n u n + B(ρ n+1 ûn , u n,k ) -div(2µ n+1 Du n,k ) + ∇(p n + q n ) = div(α n+1 Σ n,k ), u n,k | Γ = 0.
Next, a projection is used to explicitly deduce the extra-stress tensor Σ n,k+1 ∈ Λ, namely

Σ n,k+1 = P Σ n,k + rα n+1 Du n,k + θ(Σ n -Σ n,k ) . (25)
The following result establishes that the sequence (u n,k , Σ n,k ) k converges to the desired solution (u n+1 , Σ n+1 ) of the equation ( 16).

Theorem 2. If 4θ + r (α 2 ) 2 µ 1 ≤ 4, then, for all 0 ≤ n ≤ N -1, the sequence (u n,k , Σ n,k ) k tends to (u n+1 , Σ n+1 ) when k tends to infinity. Moreover the convergence is geometric with common ratio 1 -θ.

Proof. We denote u k = u n,ku n+1 and Σ k = Σ n,k -Σ n+1 . By subtracting ( 24)-( 25) to ( 16), we obtain

                 1 2δt (ρ n+1 + ρ n )u k + B(ρ n+1 ûn , u k ) -div(2µ n+1 Du k ) = div(α n+1 Σ k ), u k | Γ = 0, Σ k+1 = P Σ n,k + rα n+1 Du n,k + θ(Σ n -Σ n,k ) -P Σ n+1 + rα n+1 Du n+1 + θ(Σ n -Σ n+1 ) . (26)
We now take the inner product of the first equation in [START_REF] Pyo | Gauge-Uzawa methods for incompressible flows with variable density[END_REF] with u k in L 2 (Ω). Using [START_REF] Fuchs | Regularity results for the quasi-static Bingham variational inequality in dimensions two and three[END_REF], we deduce

(27) 1 δt ρ n+1 + ρ n 2 u k 2 + 4 µ n+1 Du k 2 = -(α n+1 Σ k , Du k ).
From the second equation of ( 26), since P is a projection, we derive

|Σ k+1 | ≤ |(1 -θ)Σ k + rα n+1 Du k |.
By taking the L 2 (Ω)-norm, it follows (28)

Σ k+1 2 ≤ (1 -θ) 2 Σ k 2 + r 2 α n+1 Du k 2 + r(1 -θ)(α n+1 Σ k , Du k ).
Combining ( 27), [START_REF] Saramito | Progress in numerical simulation of yield stress fluid flows[END_REF] and making use of α n+1 ≤ α 2 , µ n+1 ≥ µ 1 and ρ n+1 +ρ n 2 ≥ ρ 1 , we obtain

Σ k+1 2 + rρ 1 (1 -θ) δt u k 2 + r 4µ 1 (1 -θ) -rα 2 2 Du k 2 ≤ (1 -θ) 2 Σ k 2 .
From the assumption 4θ + r

α 2 2 µ 1 ≤ 4, we deduce Σ k ≤ (1 -θ) k Σ 0 and u k ≤ δt(1 -θ) rρ 1 (1 -θ) k Σ 0 which concludes the proof of Theorem 2.
Note that, if 4θ + r α 2 2 µ 1 < 4, we also have convergence in H 1 (Ω) for the velocity, namely

Du k ≤ 1 -θ r(4µ 1 (1 -θ) -rα 2 2 ) (1 -θ) k Σ 0 .

Stability analysis

In Theorem 1, we have proved that the density satisfies a discrete maximum principle which ensures, as it is noted in Remark 3, that the plastic viscosity and the yield stress have lower and upper bounds. By using this preliminary stability result, we can now establish a stability result for all sequences computed with the scheme ( 14)-( 19).

Theorem 3. If u 0 and ρ 0 satisfy (6), if r

α 2 2 µ 1 ≤ 3 2 and θ ≤ 1 2 , then, for all 1 ≤ n ≤ N , we have √ ρ n u n 2 + δt n-1 k=0 2µ n+1 Du k+1 2 + δt 2 ρ 1 ∇p n 2 + 2θ r δt Σ n 2 ≤ ρ 0 u 0 2 + δt 2 ρ 1 ∇p 0 2 + 2θ r δt Σ 0 2 .
Proof. We take the inner product of the first equation in ( 16) with 2δtu n+1 in L 2 (Ω). By observing that

2u n+1 , 1 2 (ρ n+1 + ρ n )u n+1 -ρ n u n = σ n+1 u n+1 2 -σ n u n 2 + σ n (u n+1 -u n ) 2 ,
where we denote σ n = √ ρ n , and recalling equation ( 12), we deduce ( 29)

σ n+1 u n+1 2 -σ n u n 2 + σ n (u n+1 -u n ) 2 + 4δt 2µ n+1 Du n+1 2 = -2δt(∇(p n + q n ), u n+1 ) -2δt(α n+1 Σ n+1 , Du n+1 ).
Now, we have to control the two terms in the right-hand side of ( 29) in order to obtain the expected inequality.

Term -2δt(∇(p n + q n ), u n+1 ). Recalling that p n -p n-1 = q n (see ( 18)), we write

-2δt(∇(2p n -p n-1 ), u n+1 ) = 2δt(∇(p n+1 -2p n + p n-1 ), u n+1 ) -2δt(∇p n+1 , u n+1 ).
Taking the inner product of the first equation in ( 17)

with 2δt 2 ρ 1 (p n+1 -p n ) -(p n -p n-1 ) in L 2 (Ω), we obtain δt 2 ρ 1 ∇(p n+1 -p n ) 2 -∇(p n -p n-1 ) 2 + ∇(p n+1 -2p n + p n-1 ) 2 = 2δt(u n+1 , ∇(p n+1 -2p n + p n-1 )).
Taking the inner product of the first equation in [START_REF] Guermond | A projection FEM for variable density incompressible flows[END_REF] with -2δt 2 ρ 1 p n+1 in L 2 (Ω), we obtain

- δt 2 ρ 1 ∇p n+1 2 -∇p n 2 + ∇(p n+1 -p n ) 2 = -2δt(u n+1 , ∇p n+1 ).
Using the last three equalities, we deduce

(30) -2δt(∇(2p n -p n-1 ), u n+1 ) = δt 2 ρ 1 ∇(p n+1 -2p n + p n-1 ) 2 + δt 2 ρ 1 -∇p n+1 2 + ∇p n 2 -∇(p n -p n-1 ) 2 .
Using ( 17), we have

(31) δt 2 ρ 1 ∇(p n+1 -2p n + p n-1 ) 2 ≤ σ n (u n+1 -u n ) 2 .
Combining ( 29), ( 30) and ( 31), we finally deduce (32)

σ n+1 u n+1 2 -σ n u n 2 + 4δt 2µ n+1 Du n+1 2 + δt 2 ρ 1 ∇p n+1 2 -∇p n 2 + ∇(p n -p n-1 ) 2 ≤ -2δt(α n+1 Σ n+1 , Du n+1 ).
Term -2δt(α n+1 Σ n+1 , Du n+1 ). Using the second equation of ( 16), we have

|Σ n+1 | ≤ Σ n+1 + rα n+1 Du n+1 + θ(Σ n -Σ n+1 ) .
By taking the square of this inequality, using the identity 2a(a -b) = a 2 -b 2 + (a -b) 2 , the Cauchy-Schwarz and Young's inequalities and integrating over Ω, we obtain

(33) -2δt(α n+1 Σ n+1 , Du n+1 ) ≤ 2r α 2 2 µ 1 δt 2µ n+1 Du n+1 2 + 2θ(2θ -1) r δt Σ n -Σ n+1 2 - 2θ r δt Σ n+1 2 + 2θ r δt Σ n 2 .
Combining (32) and (33), we obtain

σ n+1 u n+1 2 -σ n u n 2 + 2δt 2 -r α 2 2 µ 1 2µ n+1 Du n+1 2 + δt 2 ρ 1 ∇p n+1 2 -∇p n 2 + ∇(p n -p n-1 ) 2 + 2θ(1 -2θ) r δt Σ n -Σ n+1 2 + 2θ r δt Σ n+1 2 - 2θ r δt Σ n 2 ≤ 0.
We conclude the proof by using r α 2 2 µ 1 ≤ 3 2 and θ ≤ 1 2 , and summing from 0 up to n -1 the last inequality.

Note that the above estimate does not provide a bound on the pressure independent on δt. We will derive such a stability result by using the final error estimate (47) in the proof of Theorem 5 (see Remark 5). Moreover, using the second equation in ( 16), Σ n ∈ Λ. Hence, for all 0 ≤ n ≤ N , we have

Σ n L ∞ (Ω) ≤ 1.

Error analysis

As usual in error estimates, we have to assume that the solutions of ( 11) have a sufficient regularity. We assume in this section that the solution (u, ρ, p, Σ) of ( 11) satisfies

(34) u ∈ L ∞ (0, T ; L ∞ (Ω)), ∇u ∈ L ∞ (0, T ; L ∞ (Ω)), ∂ t u ∈ L ∞ (0, T ; L 3 (Ω)), ∂ tt u ∈ L 2 (0, T ; L 6/5 (Ω)), ∇ρ ∈ L ∞ (0, T ; L ∞ (Ω)), ∂ tt ρ ∈ L 2 (0, T ; L 2 (Ω)), ∂ t p ∈ L 2 (0, T, H 1 (Ω)), ∂ t Σ ∈ L 2 (0, T ; L 2 (Ω)).
Under this regularity assumptions on u (see [START_REF] Boyer | Éléments d'analyse pour l'étude de quelques modèles d'écoulements de fluides visqueux incompressibles[END_REF]), the solution ρ of the mass conservation equation ( 1) satisfies the maximum principle, that is almost everywhere in (0, T ) × Ω, we have

ρ 1 ≤ ρ ≤ ρ 2 .
We introduce the following errors

ε n = ρ(t n ) -ρ n , e n = u(t n ) -u n , ên = u(t n ) -ûn , r n = p(t n ) -p n , s n = Σ(t n ) -Σ n .
We

µ(ρ(t n )) -µ n ≤ µ ′ ∞ ε n and α(ρ(t n )) -α n ≤ α ′ ∞ ε n .
In the rest of the paper, C denotes a constant independent of the time step δt. Since p| t=0 is not prescribed in the initial conditions, we can not choose p 0 = p| t=0 . Therefore r 0 = 0 so that we assume (35)

r 0 2 + ∇r 0 2 ≤ Cδt 2 .
We will also assume (36) δt < 1/6.

Remark 4. Due to [START_REF] Guermond | Error analysis of a fractional time-stepping technique for incompressible flows with variable density[END_REF] and the above definitions, we have e n = ênδt ρ 1 ∇q n . Recalling that div ên = 0 and ên • n| Γ = 0 (see Remark 2), we deduce

ên 2 ≤ ên 2 + δt 2 ρ 2 1 ∇q n 2 = e n 2 .
6.1. Preliminary error estimates for the density. As for the stability analysis, we first focus on the sequence (ρ n ) {0≤n≤N } governed by equation ( 14).

Theorem 4. Under the regularity assumptions (34), if u 0 and ρ 0 satisfy (6), then, for all n such that 0 ≤ n ≤ N -1, we have

(37) δt n k=0 ρ(t k+1 ) -ρ k+1 2 ≤ Cδt 2 + Cδt n k=0 u(t k+1 ) -u k+1 2 .
Proof. The L ∞ (Ω)-regularity of ρ n+1 is not sufficient to derive the expected result. For this reason, the regularized solution ρ n ǫ ∈ H 1 (Ω), solution of ( 20), is considered. We denote

ε n ǫ = ρ(t n ) -ρ n ǫ .
By subtracting the first equation of (1) taken at time t n+1 from (20), we obtain (38)

ε n+1 ǫ -ε n δt = ρ(t n+1 ) -ρ(t n ) δt -∂ t ρ(t n+1 ) -ûn • ∇ε n+1 ǫ -ên • ∇ρ(t n+1 ) -δu(t n+1 ) • ∇ρ(t n+1 ) -ǫ∆ρ n+1 ǫ ,
where by invoking Taylor's formulae, we have

ρ(t n+1 ) -ρ(t n ) δt -∂ t ρ(t n+1 ) = 1 δt t n+1 tn (t n -t)∂ tt ρ(t)dt.
Now, by taking the L 2 (Ω)-inner product of the error equation ( 38) with 2δtε n+1 ǫ and using the above equality, we obtain

ε n+1 ǫ 2 -ε n 2 + ε n+1 ǫ -ε n 2 = 2 t n+1 tn (t n -t)∂ tt ρ(t)dt, ε n+1 ǫ -2δt Ω ε n+1 ǫ ên • ∇ρ(t n+1 ) + δu(t n+1 ) • ∇ρ(t n+1 ) + 2δtǫ(∇ρ n+1 ǫ , ∇ε n+1 ǫ ). By writing 2δtǫ(∇ρ n+1 ǫ , ∇ε n+1 ǫ ) = 2δtǫ(∇ρ n+1 ǫ , ∇ρ(t n+1 )) -2δt √ ǫ∇ρ n+1 ǫ 2 ,
we obtain (39)

ε n+1 ǫ 2 -ε n 2 + ε n+1 ǫ -ε n 2 + 2δt √ ǫ∇ρ n+1 ǫ 2 = 2 t n+1 tn (t n -t)∂ tt ρ(t)dt, ε n+1 ǫ -2δt Ω ên • ∇ρ(t n+1 )ε n+1 ǫ -2δt Ω δu(t n+1 ) • ∇ρ(t n+1 )ε n+1 ǫ + 2δtǫ(∇ρ n+1 ǫ , ∇ρ(t n+1 )).
Making ǫ tends to 0 in (39), using [START_REF] Marly | Augmented Lagrangian simulations study of yield-stress fluid flows in expansioncontraction and comparisons with physical experiments[END_REF] and the inequality

ε n+1 ≤ lim inf ǫ ε n+1 ǫ , we derive (40) ε n+1 2 -ε n 2 ≤ 2 t n+1 tn (t n -t)∂ tt ρ(t)dt, ε n+1 -2δt Ω ên • ∇ρ(t n+1 )ε n+1 -2δt Ω δu(t n+1 ) • ∇ρ(t n+1 )ε n+1 .
We now estimate the terms in the right-hand side of (40)

2 Ω t n+1 tn (t n -t)∂ tt dtρ(t)dt ε n+1 ≤ δt 2 t n+1 tn ∂ tt ρ(t) 2 dt + δt ε n+1 2 , -2δt Ω ên • ∇ρ(t n+1 )ε n+1 ≤ δt ∇ρ 2 L ∞ (0,T,L ∞ (Ω)) ên 2 + δt ε n+1 2 ,
and

-2δt Ω δu(t n+1 ) • ∇ρ(t n+1 )ε n+1 ≤ 2δt 3 2 t n+1 tn ∂ t u(t) 2 L 4 (Ω) dt 1 2 ∇ρ L ∞ (0,T,L 4 (Ω)) ε n+1 ≤ δt 2 ∇ρ 2 L ∞ (0,T,L 4 (Ω)) t n+1 tn ∂ t u(t) 2 L 4 (Ω) dt + δt ε n+1 2 .
Reporting the above inequalities in (40) and summing over n from 0 to j, with 0 ≤ j ≤ N -1, we have

ε j+1 2 ≤ δt 2 ∂ tt ρ 2 L 2 (0,T,L 2 (Ω)) + ∂ t u 2 L 2 (0,T,L 4 (Ω)) ∇ρ 2 L ∞ (0,T,L 4 (Ω)) + 3δt j n=0 ε n+1 2 + δt ∇ρ 2 L ∞ (0,T,L ∞ (Ω)) j n=0 ên 2 .
Using the discrete Gronwall Lemma 1, the assumption (36), and Remark 4 we derive

ε j+1 2 ≤δt 2 exp(6T ) ∂ tt ρ 2 L 2 (0,T,L 2 (Ω)) + ∂ t u 2 L 2 (0,T,L 4 (Ω)) ∇ρ 2 L ∞ (0,T,L 4 (Ω)) + δt exp(6T ) ∇ρ 2 L ∞ (0,T,L ∞ (Ω)) j n=0 e n 2 .
We conclude the proof of Theorem 4 by multiplying by δt and summing over j from 0 to k, with k ≤ N -1, the above inequality.

6.2. Final error estimates.

Theorem 5. Under the regularity assumption (34), if the initial conditions satisfy (6) and (35),

if θ < 1 2 , δt < 1/6, rα 2 2 ≤ µ 1 12(1+C(Ω) 2 )
where C(Ω) is the Poincaré's constant, then, for all n such that 1 ≤ n ≤ N , the discrete velocity solution of the scheme (14)-( 19) satisfies the following error estimate

u(t n ) -u n 2 + δt N -1 k=0 u(t k+1 ) -u k+1 2 1 ≤ Cδt(δt + θ).
Proof. We first derive the error equation which is obtained by subtracting the two first equations in [START_REF] Frigaard | On the usage of viscosity regularisation methods for visco-plastic fluid flow computation[END_REF] to the first equation in ( 16), namely we find (41)

1 δt 1 2 (ρ n+1 + ρ n )e n+1 -ρ n e n -div 2µ n+1 De n+1 + ∇ p(t n+1 ) -(p n + q n ) = div α(ρ(t n+1 ))Σ(t n+1 ) -α n+1 Σ n+1 + div 2 µ(ρ(t n+1 )) -µ n+1 Du(t n+1 ) + R n+1 1 + R n+1 2 + R n+1 3 + R n+1 4 ,
where

R n+1 1 = ρ n u(t n+1 ) -u(t n ) δt -ρ(t n+1 )∂ t u(t n+1 ), R n+1 2 = 1 2 u(t n+1 ) ρ n+1 -ρ n δt -∂ t ρ(t n+1 ) , R n+1 3 = ρ n+1 ûn • ∇u n+1 -ρ(t n+1 )u(t n+1 ) • ∇u(t n+1 ), R n+1 4 = 1 2 u n+1 div(ρ n+1 ûn ) - 1 2 u(t n+1 ) div u(t n+1 )ρ(t n+1 ) .
By virtue of the Korn's inequality (see [START_REF] Lewicka | On the optimal constants in Korn's and geometric rigidity estimates, in bounded and unbounded domains, under Neumann boundary conditions[END_REF]), since e n+1 ∈ H 1 0 (Ω), we have ∇e n+1 2 ≤ 2 De n+1 2 . Invoking the Poincaré's inequality ( e n+1 ≤ C(Ω) ∇e n+1 ) and denoting

C P = (1 + C(Ω) 2 ) -1 , we obtain 4µ 1 C P δt e n+1 2 1 ≤ 4δt 2µ n+1 De n+1 2 .
We now take the inner product of equation ( 41) with 2δte n+1 in L 2 (Ω). Using the last inequality we obtain (42)

σ n+1 e n+1 2 + σ n (e n+1 -e n ) 2 -σ n e n 2 + 4µ 1 C P δt e n+1 2 1 ≤ -2δt ∇ p(t n+1 ) -(p n + q n ) , e n+1 -2δt α(ρ(t n+1 ))Σ(t n+1 ) -α n+1 Σ n+1 , De n+1 -4δt µ(ρ(t n+1 )) -µ n+1 Du(t n+1 ), De n+1 + 2δt R n+1 1 , e n+1 + 2δt R n+1 2 , e n+1 + 2δt R n+1 3 , e n+1 + 2δt R n+1 4 , e n+1 .
In the following steps, we derive estimates of the terms occurring in the right-hand side of (42). They are treated in order of appearance in (42).

Step 1. From ( 17) we obtain, for all g ∈ H 1 (Ω), (43) ∇(r n+1 -r n ), ∇g = ρ 1 δt (e n+1 , ∇g) + (∇δp(t n+1 ), ∇g).

Taking g = -2δt 2 ρ 1 δ 2 r n+1 in (43), we have

- δt 2 ρ 1 ∇(r n+1 -r n ) 2 -∇(r n -r n-1 ) 2 + ∇δ 2 r n+1 2 = -2δt(e n+1 , ∇δ 2 r n+1 ) - 2δt 2 ρ 1 ∇δp(t n+1 ), ∇δ 2 r n+1 .
Taking g = 2δt 2 ρ 1 r n+1 in (43), we have

δt 2 ρ 1 ∇r n+1 2 -∇r n 2 + ∇δr n+1 2 = 2δt(e n+1 , ∇r n+1 ) + 2δt 2 ρ 1 ∇δp(t n+1 ), ∇r n+1 .
By summing these last two equations, we deduce (44)

δt 2 ρ 1 ∇r n+1 2 -∇r n 2 + ∇δr n 2 - δt 2 ρ 1 ∇δ 2 r n+1 2
= 2δt e n+1 , ∇(2r n -r n-1 ) + 2δt 2 ρ 1 ∇δp(t n+1 ), ∇(2r n -r n-1 ) .

From (43), we have ∇δ 2 r n+1 , ∇g = ρ 1 δt (δe n+1 , ∇g) + (∇δ 2 p(t n+1 ), ∇g).

Taking g = δt √ ρ 1 δ 2 r n+1 in the above relation, we obtain (45)

δt 2 ρ 1 ∇δ 2 r n+1 2 ≤ 1 ρ 1 ρ 1 δe n+1 + δt∇δ 2 p(t n+1 ) 2 ≤ σ n δe n+1 2 + δt 2 ρ 1 ∇δ 2 p(t n+1 ) 2 + 2δt δe n+1 , ∇δ 2 p(t n+1 ) .
By summing (44) and (45), we deduce

δt 2 ρ 1 ∇r n+1 2 -∇r n 2 + ∇δr n 2 ≤ 2δt e n+1 , ∇(2r n -r n-1 ) + 2δt 2 ρ 1 δp(t n+1 ), ∇(2r n -r n-1 ) + σ n δe n+1 2 + δt 2 ρ 1 ∇δ 2 p(t n+1 ) 2 + 2δt δe n+1 , ∇δ 2 p(t n+1 ) ,
which can be rewritten as -2δt e n+1 , ∇(2r n -r n-1 ) -2δt δe n+1 , ∇δ 2 p(t n+1 ) -2δt e n , ∇δ 2 p(t n+1 )

≤ -2δt e n , ∇δ 2 p(t n+1 )

a 1 - δt 2 ρ 1 ∇r n+1 2 -∇r n 2 + ∇δr n 2 + 2δt 2 ρ 1 δp(t n+1 ), ∇(2r n -r n-1 ) a 2 + σ n δe n+1 2 + δt 2 ρ 1 ∇δ 2 p(t n+1 ) 2 a 3
.

Using the Cauchy-Schwarz and Young's inequalities, we obtain

a 1 ≤ δt e n 2 + 2δt 2 t n+1 t n-1 ∂ t p(s) 2 1 ds, a 2 = 2δt 2 ρ 1 δp(t n+1 ), ∇r n + δp(t n+1 ), ∇δr n ≤ 2δt 2 ρ 1 t n+1 tn ∂ t p(s) 2 ds + δt 3 ρ 1 ∇r n 2 + δt 3 ρ 1 ∇δr n 2 , a 3 ≤ 2δt 3 ρ 1 t n+1 t n-1 ∂ t p(s) 2 1 ds.
We finally deduce

-2δt ∇ p(t n+1 ) -(p n + q n ) , e n+1 ≤ - δt 2 ρ 1 ∇r n+1 2 -∇r n 2 + δt 3 ρ 1 ∇r n 2 + σ n δe n+1 2 + δt 2 2 + 4δt ρ 1 + 1 t n+1 t n-1 ∂ t p(s) 2 1 ds + δt e n 2 .
Step 2. We write (46)

-2δt α(ρ(t n+1 ))Σ(t n+1 ) -α n+1 Σ n+1 , De n+1 = -2δt α(ρ(t n+1 )) -α n+1 Σ(t n+1 ), De n+1 -2δt s n+1 , α n+1 De n+1 ≤ δt µ 1 C P 2 e n+1 2 1 + δt 2 α ′ 2 ∞ µ 1 C P ε n+1 2 -2δt s n+1 , α n+1 De n+1 .
In order to treat the last term in the right-hand side of (46), we use the relation ( 9) on Σ with ℓ = rα n+1 , namely

Σ(t n+1 ) = P Σ(t n+1 ) + rα n+1 Du(t n+1 ) , so that s n+1 = P Σ(t n+1 ) + rα n+1 Du(t n+1 ) -P Σ n+1 + rα n+1 Du n+1 + θ(Σ n -Σ n+1 ) .
Since P is a projection, we have

|s n+1 | ≤ |s n+1 + rα n+1 De n+1 + θ(s n -s n+1 ) + θδΣ(t n+1 )|.
Taking the square of this inequality leads to

|s n+1 | 2 ≤|s n+1 | 2 + r 2 (α n+1 ) 2 |De n+1 | 2 + θ 2 |s n -s n+1 | 2 + θ 2 |δΣ(t n+1 )| 2
+ rs n+1 : α n+1 De n+1 + θs n+1 : (s ns n+1 ) + θs n+1 : δΣ(t n+1 ) + rθα n+1 De n+1 : (s ns n+1 ) + rθα n+1 De n+1 : δΣ(t n+1 ) + θ 2 (s ns n+1 ) : δΣ(t n+1 ).

Using the Cauchy-Schwarz and the Young's inequalities, the identity 2a(a-b) = a 2 -b 2 +(a-b) 2 , integrating over Ω, and multiplying by 2δt r , we obtain

-2δt(s n+1 , α n+1 De n+1 ) ≤ 6δtrα 2 2 De n+1 2 + 2θ r (3θδt + 1) δΣ(t n+1 ) 2 - 2θ r δt(1 -δt) s n+1 2 + 2θ r δt s n 2 .
By reporting the last inequality in (46) we deduce

-2δt α(ρ(t n+1 ))Σ(t n+1 ) -α n+1 Σ n+1 , De n+1 ≤ µ 1 C P 2 + 6rα 2 2 δt e n+1 2 1 + δt 2 α ′ 2 ∞ µ 1 C P ε n+1 2 + δt 2θ(3θδt + 1) r t n+1 tn ∂ t Σ(s) 2 ds - 2δtθ r (1 -δt) s n+1 2 -s n 2 .
Step 3. The third term in the right-hand side of (42) is bounded as follows

-4δt µ(ρ(t n+1 )) -µ n+1 Du(t n+1 ), De n+1 ≤δt 8 µ ′ 2 ∞ µ 1 C P Du 2 L ∞ (0,T,L ∞ (Ω)) ε n+1 2 + δt µ 1 C P 2 e n+1 2 1 .
Step 4. Using the Sobolev embedding H 1 (Ω) ⊂ L 6 (Ω) in addition to the Cauchy-Schwarz inequality, the Young's inequality and Taylor's formulae, we deduce

2δt(R n+1 1 , e n+1 ) ≤ 6C S µ 1 C P δt ε n 2 ∂ t u 2 L ∞ (0,T,L 3 (Ω)) + δt µ 1 C P 2 e n+1 2 1 + 6C S µ 1 C P δt 2 ρ 2 2 t n+1 tn ∂ tt u(s) 2 L 6/5 (Ω) ds + ∂ t u 2 L ∞ (0,T,L 3 (Ω)) t n+1 tn ∂ t ρ(s) 2 ds ,
where C S is the continuity constant from the Sobolev embedding

H 1 (Ω) ⊂ L 6 (Ω), that is, for all f ∈ H 1 (Ω), we have f 2 L 6 (Ω) ≤ C S f 2 1 .
Step 5. Using ( 14) and the first equation of ( 11), we write

R n+1 2 = 1 2 u(t n+1 ) δu(t n+1 ) • ∇ρ(t n+1 ) + u(t n ) • ∇ε n+1 + ên • ∇ρ n+1 .
Recalling that both u and ên are divergence-free vector fields, we obtain

R n+1 2 = 1 2 u(t n+1 ) div δu(t n+1 )ρ(t n+1 ) + u(t n )ε n+1 + ên ρ n+1 .
From the above equality, performing an integration by parts and using the identity

∇(u • v) = ∇u • v + ∇v • u, for any u and v in H 1 (Ω), we derive 2δt(R n+1 2 , e n+1 ) = -δt Ω ρ n+1 ên • ∇u(t n+1 ) • e n+1 + ∇e n+1 • u(t n+1 ) + Ω ε n+1 u(t n ) • ∇u(t n+1 ) • e n+1 + ∇e n+1 • u(t n+1 ) + Ω ρ(t n+1 )δu(t n+1 ) • ∇u(t n+1 ) • e n+1 + ∇e n+1 • u(t n+1 ) .
Proceeding as in the previous steps, we deduce

2δt(R n+1 2 , e n+1 ) ≤ µ 1 C P 2 δt e n+1 2 1 + 3δt µ 1 C P u 2 L ∞ (0,T,L ∞ (Ω)) + C S ∇u 2 L ∞ (0,T,L 3 (Ω)) × ρ 2 2 ên 2 + u 2 L ∞ (0,T,L ∞ (Ω)) ε n+1 2 + 3ρ 2 2 µ 1 C P δt 2 u 2 L ∞ (0,T,L ∞ (Ω)) + C S ∇u 2 L ∞ (0,T,L 3 (Ω)) t n+1 tn ∂ t u(s) 2 ds. µ 1 C P u 2 L ∞ (0,T,L ∞ (Ω)) + 9C S ∇u 2 L ∞ (0,T,L 3 (Ω)) t n+1 tn ∂ t u(s) 2 ds.
Step 7. The final estimate is obtained by applying the discrete Gronwall Lemma 1. Let us first consider the particular case n = 0. Using hypothesis (35) and the previous error equations, we deduce e 1 2 + δt e 1 2 1 + δtθ s 1 2 + δt 2 ∇r 1 2 ≤ Cδt(θ + δt). General case (n ≥ 1). Now, using steps 1 to 7, we finally obtain Remark 5. From the inequality (47), we deduce that, for all n such that 0 ≤ n ≤ N -1,

σ n+1 e n+1 2 -
∇r n+1 2 ≤ C δt + θ δt .
Taking θ = δt, we obtain, for all n such that 0 ≤ n ≤ N -1,

∇p n+1 ≤ C,
where C is independent of δt.

Conclusion

We have proposed and analyzed in this paper a new scheme for the temporal discretization of incompressible equations describing the motion of Bingham fluids with variable density, plastic viscosity and yield stress. A fractional time-stepping algorithm is coupled with a projection formulation for the definition of the plastic tensor. This approach, also used in the Uzawa-like algorithm for solving viscoplastic flows, allows to handle the non-differentiable definition of the stress tensor in the Bingham constitutive laws. The first sub-time step of the bi-projection scheme consists in computing a non-solenoidal velocity field and is followed by a projection step so that the final velocity field is divergence free. In classical projection methods, the Helmholtz decomposition is invoked for this projection resulting in a Poisson equation satisfied by the pseudo-pressure. Extending this approach to density-variable flows will lead to a second-order elliptic equation with variable coefficients for the pressure. Another approach based on the interpretation of projection methods in terms of a penalty method, Guermond and Salgado in [START_REF] Guermond | A splitting method for incompressible flows with variable density based on a pressure Poisson equation[END_REF][START_REF] Guermond | Error analysis of a fractional time-stepping technique for incompressible flows with variable density[END_REF] derived a fractional time-stepping scheme where the computation of the pressure is achieved by solving one standard Poisson equation per time step. Unlike in [START_REF] Guermond | A splitting method for incompressible flows with variable density based on a pressure Poisson equation[END_REF][START_REF] Guermond | Error analysis of a fractional time-stepping technique for incompressible flows with variable density[END_REF], the divergence-free velocity is used as convective velocity in the mass conservation equation. This point is of major importance as it ensures a maximum principle for the discrete continuity equation. Lower and upper bounds on the density are obtained. The plastic part of the stress tensor is treated implicitly into the first sub-step (prediction step) of the fractional time-stepping method. As in the Uzawa-like method for Bingham flows, a fixed-point algorithm is used to compute the plastic tensor. Due to a pseudo-time relaxation term added in the Bingham projection operator, the convergence is geometric with common ratio (1 -θ) where θ is a prescribed relaxation parameter. The main results derived in this paper are the unconditional stability and the convergence of the bi-projection scheme. More precisely, errors committed by approximating the velocity and the density are bounded from above by a term of the order of δt(δt + θ). Hence, by choosing θ of the order of the time step ensures the error of the bi-projection scheme to be of first order.

  and by | • | the associated norm. For two tensors A and B in R 3×3 , we denote by A : B their inner product, namely

  denote µ ′ ∞ = sup x∈R |µ ′ (x)| and α ′ ∞ = sup x∈R |α ′ (x)|. Since µ and α are functions of class C 1 , we obtain

L 2 ≤ µ 1 C P 2 ,ρ 1 ∇r k+1 2 + 2θδt r s k+1 2 + δtµ 1 C P N - 1 n=0 e n+1 2 1 ≤

 22211 ∞ (0,T,L 3 (Ω))t n+1 tn ∂ t ρ(s) 2 ds , with K(u) = u 2 L ∞ (0,T,L ∞ (Ω)) + 9C S ∇u 2 L ∞ (0,T,L 3 (Ω)). Summing the last equation over n from 1 to k with 0 ≤ k ≤ N -1, using 6rα 2 Theorem 4, the discrete Gronwall Lemma 1 and the case n = 0, we deduce, for all 0 ≤ k ≤ N -1, (47)σ k+1 e k+1 2 + δt 2 Cδt(δt + θ),which concludes the proof of Theorem 5.

  everywhere in Ω, then there exists a unique viscosity solution satisfying[START_REF] Glowinsky | Numerical Analysis of Variational Inequalities[END_REF] in the distribution sense. Moreover, this solution satisfies ρ 1 ≤ ρ n+1 ≤ ρ 2 , a.e. in Ω.

Proof. Step 1: Regularized equation. Equation (20) leads us to consider the variational problem

  σ n e n 2 + δt 2 ρ 1 ∇r n+1 2 -∇r n 2

			+	2θδt r	s n+1 2 -s n 2 +	3 2	µ 1 C P -4rα 2 2 δt e n+1 2 1
	≤	δt 3 ρ 1	∇r n 2 +	2θδt 2 r	s n+1 2 + 1 +	6ρ 2 2 µ 1 C P	K(u)	δt 1 ρ 2	σ n e n 2
	+ L + 2 µ 1 C P 3 u 2 6C S µ 1 C P ∂ t u 2 L ∞ (0,T,L 3 (Ω)) δt ε n 2 +	6ρ 2 2 µ 1 C P	K(u)δt 2	tn t n+1	∂ t u(s) 2 ds
	+	2 r	(3θδt + 1)θδt	t n+1 tn	∂ t Σ(s) 2 ds + 1 +	2 ρ 1	(1 + 2δt) δt 2	t n+1 t n-1	∂ t p(s) 2 1 ds
	+	6C S µ 1 C P	2 δt 2 ρ 2		t n+1

∞ (0,T,L ∞ (Ω)) K(u) + α ′ 2 ∞ + 4 µ ′ 2 ∞ Du 2 L ∞ (0,T,L ∞ (Ω)) δt ε n+1 2 tn ∂ tt u(s) 2 L 6/5 (Ω) ds + ∂ t u 2

Acknowledgements

This research was financed by the French Government Laboratory of Excellence initiative n o ANR-10-LABX-0006, by the French National Research Agency (ANR) RAVEX project, and by the French National Center for Scientific Research (CNRS) TelluS project. This is Laboratory of Excellence ClerVolc contribution number 302.

Step 6. Using [START_REF] Fuchs | Regularity results for the quasi-static Bingham variational inequality in dimensions two and three[END_REF], we have

Proceeding as above, we obtain


Guillaume Zahnd
email: g.zahnd@erasmusmc.nl

Maciej Orkisz

Simone Balocco

André Sérusclat

Philippe Moulin

Didier Vray

Tracking arterial wall motion in a 2D+t volume

I. INTRODUCTION

Cardiovascular diseases are the primary cause of mortality worldwide [START_REF]Cardiovascular diseases (CVDs), Fact sheet number 317[END_REF]. Yet, cardiovascular risk prediction is challenging, as the screening potential of traditional imaging biomarkers, such as pulse wave velocity, intima-media thickness, or arterial calcification, remains limited [START_REF] Simon | Performance of subclinical arterial disease detection as a screening test for coronary heart disease[END_REF]. Therefore, there is a strong need for novel screening biomarkers to improve risk prediction.

The common carotid artery (CCA) has been extensively analyzed using B-mode ultrasound (US) imaging [START_REF] Golemati | Toward novel noninvasive and low-cost markers for predicting strokes in asymptomatic carotid atherosclerosis: The role of ultrasound image analysis[END_REF]. The deformation of the tissues in the direction parallel to the blood flow during the cardiac cycle has gained interest during the last decade [START_REF] Persson | A new non-invasive ultrasonic method for simultaneous measurements of longitudinal and radial arterial wall movements: First in vivo trial[END_REF], [START_REF] Golemati | Carotid artery wall motion estimated from B-mode ultrasound using region tracking and block matching[END_REF], [START_REF] Cinthio | Evaluation of an ultrasonic echo-tracking method for measurements of arterial wall movements in two dimensions[END_REF], [START_REF] Cinthio | Longitudinal movements and resulting shear strain of the arterial wall[END_REF], [START_REF] Idzenga | Estimating cyclic shear strain in the common carotid artery using radiofrequency ultrasound[END_REF], [START_REF] Soleimani | Carotid artery wall motion estimation from consecutive ultrasonic images: Comparison between block-matching and maximum-gradient algorithms[END_REF], [START_REF] Yli-Ollila | Axial and radial waveforms in common carotid artery: An advanced method for studying arterial elastic properties in ultrasound imaging[END_REF], [START_REF] Zahnd | Evaluation of a Kalman-based block matching method to assess the bidimensional motion of the carotid artery wall in B-mode ultrasound sequences[END_REF]. Such cyclic shearing motion of the intima-media complex with respect to the tunica adventitia, hereafter referred to as "longitudinal kinetics" (LOKI), has been reported to reflect vascular health [START_REF] Svedlund | Carotid artery longitudinal displacement predicts 1-year cardiovascular outcome in patients with suspected coronary artery disease[END_REF], [START_REF] Zahnd | Measurement of two-dimensional movement parameters of the carotid artery wall for early detection of arteriosclerosis: A preliminary clinical study[END_REF], [START_REF] Zahnd | Longitudinal displacement of the carotid wall and cardiovascular risk factors: associations with aging, adiposity, blood pressure and periodontal disease independent of cross-sectional distensibility and intima-media thickness[END_REF], [START_REF] Ahlgren | Longitudinal displacement and intramural shear strain of the porcine carotid artery undergo profound changes in response to catecholamines[END_REF], [START_REF] Gastounioti | Carotid artery wall motion analysis from B-mode ultrasound using adaptive block matching: in silico evaluation and in vivo application[END_REF], and may therefore constitute a valuable image-based biomarker for improved cardiovascular risk prediction.

Accurate in vivo quantification of LOKI in routinely acquired data is hindered by various factors, namely speckle decorrelation, out-of-plane motions, poor echogenicity, or tissue deformation. Aiming to estimate the trajectory of the intima-media complex during US image sequences (cineloops), a variety of motion estimation techniques have been proposed. A traditional block matching approach was used in [START_REF] Golemati | Carotid artery wall motion estimated from B-mode ultrasound using region tracking and block matching[END_REF], with a block of dimensions 3.2 × 2.5 mm 2 . This technique, although robust against noise, is likely to yield a precision loss due to the large dimensions of the block. Conversely, an echo tracking approach was proposed in [START_REF] Cinthio | Evaluation of an ultrasonic echo-tracking method for measurements of arterial wall movements in two dimensions[END_REF], with a block of dimensions 0.1 × 0.1 mm 2 . The wall motion could be observed with a great accuracy, however this technique is more sensitive to noise, and therefore more prone to error in clinical data. A contrast optimization scheme was introduced in [START_REF] Yli-Ollila | Axial and radial waveforms in common carotid artery: An advanced method for studying arterial elastic properties in ultrasound imaging[END_REF], to offer a better trade-off between noise robustness and tracking accuracy. Recent approaches have incorporated the temporal dimension into the estimation process: a Kalman filtering framework dedicated to avoid tracking divergence was proposed in [START_REF] Zahnd | Evaluation of a Kalman-based block matching method to assess the bidimensional motion of the carotid artery wall in B-mode ultrasound sequences[END_REF], and a finite impulse response filtering able to handle changes in the appearance of motion targets was introduced in [START_REF] Gastounioti | Carotid artery wall motion analysis from B-mode ultrasound using adaptive block matching: in silico evaluation and in vivo application[END_REF]. Nevertheless, all of these proposed techniques perform the motion estimation using the information provided by a subset of frames only (i.e. one or several previous frames, and the current one), and provide a solution corresponding to a succession of local optima. Instead, it could be more efficient to determine the optimal trajectory once all frames have been processed.

Dynamic programming [START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF], [START_REF] Cohen | Minimal paths and fast marching methods for image analysis[END_REF] is an efficient approach in combinatorial analysis, as the optimal solution among all possible combinations is determined without a costly full search. Dynamic programming has numerous applications in biomedical imaging, from tracking the motion of cells in fluorescence microscopy [START_REF] Padfield | Spatio-temporal cell cycle phase analysis using level sets and fast marching methods[END_REF], to segmenting the contours of the intima-media complex of the CCA in B-mode US [START_REF] Zahnd | Simultaneous extraction of carotid artery intima-media interfaces in ultrasound images: Assessment of wall thickness temporal variation during the cardiac cycle[END_REF].

The aim of this study is to present a tracking framework based on a dynamic programming scheme, to fully incorporate the temporal dimension into the estimation process. In this approach, the cine-loop is considered as a twodimensional plus time (2D+t) volume, in a similar fashion as previous work in the field of image registration [START_REF] Metz | Nonrigid registration of dynamic medical imaging data using nD+t B-splines and a groupwise optimization approach[END_REF]. Our rationale is to avoid tracking errors that can be locally provoked by the presence of noise or artifacts at some given time steps. The major contribution of this work is the research of a global optimum, as the resulting trajectory of the tracked point is determined a posteriori after all frames have been processed.

II. METHODS

A. Motion tracking

We introduce here our motion tracking framework. Since the cine-loop is considered as a 2D+t stack of frames, each point is described by three coordinates: temporal (k), longitudinal (x), and radial (y). Our approach is based on four main steps, as described below.

1) Initialization: First, the point p = (1, p x , p y ) to be tracked is manually selected within the intima-media complex, in the first frame of the cine-loop. The selection criteria are the following: the point should correspond to a bright echo scatterer that is well contrasted with its surroundings and remains perceptible during the entire cineloop (Fig. 1b) [START_REF] Cinthio | Evaluation of an ultrasonic echo-tracking method for measurements of arterial wall movements in two dimensions[END_REF], [START_REF] Zahnd | Evaluation of a Kalman-based block matching method to assess the bidimensional motion of the carotid artery wall in B-mode ultrasound sequences[END_REF]. Centered on the initial point, a block B of dimensions [B x , B y ] is extracted and stored as a reference pattern. After this initial operation, the tracking is automatically carried out, as described below.

2) Cost computation: A 2D+t cost function C is then generated, by computing the dissimilarity between the initial pattern B and all the patterns corresponding to each point location in each frame of the cine-loop (Fig. 1d). Dissimilarity is computed by means of the sum of squared differences. In this function C, regions likely to correspond to the actual trajectory of the tracked point have a lower cost, as their speckle pattern is similar to the reference pattern.

3) Dynamic programming: From the first frame forwards, a cumulated cost map C is generated via dynamic programming [START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF], [START_REF] Cohen | Minimal paths and fast marching methods for image analysis[END_REF] using a front propagation scheme (Fig. 1e). The values of C corresponding to the first frame (k = 1) are initialized according to Equation 1.

C(1, x, y) = 0, for [x, y] = [p x , p y ] C(1, x, y) = ∞, for [x, y] = [p x , p y ] (1)
For k > 1, the map C is iteratively generated using the relation described in Equation 2:

C(k, x, y) = min dx,dy C(k -1, x + dx, y + dy) + C(k, x, y) • 1 + κ x • dx 2 + κ y • dy 2 , (2)
with, dx ∈ {-Nx 2 , . . . 0, . . . Nx 2 }, dy ∈ {-Ny 2 , . . . 0, . . . Ny 2 }, N x and N y being the number of reachable neighbors in the longitudinal and radial directions, respectively, and κ x and κ y the smoothness factor in the longitudinal and radial directions, respectively.

4) Back-tracking: Finally, the trajectory of the tracked point is extracted. In the last frame of the cine-loop, the point p with the minimal cumulated cost is determined, and the optimal path is found by back-tracking the decreasing values in C, from p to p.

B. Data acquisition

Our method was evaluated on 15 healthy volunteers (4 males, mean age 50 ± 3 y.o.). B-mode US cine-loops of the left CCA (duration of at least 2 cardiac cycles) were acquired by an expert physician, using a clinical scanner (Antares, Siemens, Erlangen, Germany) with a linear probe of central frequency 8 MHz. The pixel dimensions were 30 × 30 µm, and the frame rate was 26 fps. The average number of frames per cine-loop was 145 ± 19.

C. Performance evaluation

Each resulting trajectory was compared to the corresponding reference trajectory, previously constructed by three experts performing the manual tracking of the same point [START_REF] Zahnd | Evaluation of a Kalman-based block matching method to assess the bidimensional motion of the carotid artery wall in B-mode ultrasound sequences[END_REF]. The tracking error was defined as the average absolute difference between the estimated trajectory and the reference trajectory, for both longitudinal and radial trajectory components.

D. Parameter settings

Our method was applied on all cine-loops with the following parameter settings: block dimension B x × B y = 1.50 × 0.30 mm 2 ; smoothness constraint [κ x , κ y ] = [START_REF]Cardiovascular diseases (CVDs), Fact sheet number 317[END_REF][START_REF]Cardiovascular diseases (CVDs), Fact sheet number 317[END_REF]; number of reachable neighbors [N x , N y ] = [33, 33], corresponding to 1 mm in both longitudinal and radial directions. The values B x , B y , N x and N y have been determined from previous work [START_REF] Zahnd | Evaluation of a Kalman-based block matching method to assess the bidimensional motion of the carotid artery wall in B-mode ultrasound sequences[END_REF]. The values κ x and κ y have been heuristically determined.

Fig. 2. Representative results of the longitudinal (first row) and radial (second row) components of the estimated trajectories (solid blue), compared with the reference trajectories (dotted black). (a,b) Successful estimations. (c) Our method fails to capture the entire amplitude of the second cycle, probably due to a large tissue motion. However, the correct trajectory is then successfully recovered thanks to the research of a global a posteriori optimum. (d) The estimated trajectory is similar to the reference, but the presence of a constant shift indicates that a distant point has been tracked instead of the initially selected one. This is probably be due to the fact that the initial pattern was not sufficiently discriminative. Fig. 3. Average tracking error for each cine-loop, in the longitudinal (a) and radial (b) directions. It is clearly visible that the errors are lower in the radial direction. This can be explained by the fact that the radial profile of the image presents stronger edges (i.e. lumen-intima and intima-media anatomical interfaces), whereas the longitudinal profile is more homogeneous.

III. RESULTS

Representative examples of the resulting trajectories are displayed in Figure 2. The average absolute tracking errors were 104 ± 127 µm and 29 ± 57 µm, in the longitudinal and radial directions, respectively. This accuracy level is similar to the one of state-of-the-art methods [START_REF] Zahnd | Evaluation of a Kalman-based block matching method to assess the bidimensional motion of the carotid artery wall in B-mode ultrasound sequences[END_REF]. The average tracking errors for each cine-loop are displayed in Figure 3.

These errors are also to be compared with the average peak-to-peak amplitude of the tissue motion, which was 609 ± 274 µm and 350 ± 320 µm, in the longitudinal and radial directions, respectively. Accordingly, the tracking error corresponded roughly to 17% and 8% of the total motion amplitude, in the longitudinal and radial directions, respectively.

Implemented in Matlab (MATLAB 7.13, The MathWorks 85 Inc., Natick, MA, USA, 2011), our method required, on average, 57 s to process a given cine-loop.

IV. DISCUSSION

A motion tracking method was proposed to assess the trajectory of a point located in the intima-media complex of the arterial wall in B-mode US imaging. The principal characteristic of this framework is the global evaluation of an optimal trajectory after all frames of the cine-loop have been processed. Under the context of a dynamic programming approach, we proposed a front-propagation scheme adapted from the fast-marching methodology [START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF], [START_REF] Cohen | Minimal paths and fast marching methods for image analysis[END_REF].

In this work, the parameter settings have been determined heuristically. As the smoothness constraint terms κ x and κ y are used to penalize abrupt displacements, it is likely that the values used here would be sub-optimal when applied to a different dataset where the tissue motion is different (i.e. either due to the frame rate of the sonograph, or to stiffer arteries in patients). A thorough evaluation of the actual impact of the parameter settings in a larger dataset, including diseased participants, is thus required in further studies.

The major drawback of this framework is the assumption that the tracked pattern remains similar during the entire length of the cine-loop. Actually, the initial pattern of the block B selected in the first frame is used to generate the 2D+t cost function C. As moving tissues may undergo a shearing deformation during the heart beat, the appearance of the corresponding pattern is likely to vary. Therefore, an adaptive reference pattern would probably be better tailored to keep track of the moving point. This issue will be tackled in future work.

In conclusion, fully incorporating the temporal dimension into the estimation process by considering the set of frames as a 2D+t volume could constitute a promising approach for improved estimation of the arterial wall motion in US imaging.

Fig. 1 .

 1 Fig. 1. Schematic representation of the motion tracking framework. (a) Common carotid artery in B-mode ultrasound imaging (1 st frame of the cine-loop). (b) Detailed region of the far wall. The tracked block B is represented by the red rectangle within the intima-media complex. (c) k th frame of the cine-loop (d) Corresponding cost function C. (e) Cumulated cost function C. (f) Schematic representation of the trajectory extracted from C after all frames have been processed.

ACKNOWLEDGMENT

This work was done within the French ANR LABEX PRIMES (ANR-11-LABX-0063) and CeLyA (ANR-10-LABX-0060) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR).