Majda Kermadi
email: kermadimajda@gmail.com

Saïd Moussaoui
email: said.moussaoui@irccyn.ec-nantes.fr

Abdelhalim Taieb Brahimi
email: abdelhalim.taiebbrahimi@univ-usto.dz

Mouloud Feliachi
email: mouloud.feliachi@univ-nantes.fr

Application of Kernel Change Detection Method in Eddy Current Non-destructive Testing (EC-NDT)

Keywords: eddy current testing, finite element method, inverse problem, kernel change detection

This paper presents an adaptation of kernel change detection (KCD) method to solve the problem of defect localization in non-destructive testing (NDT). The main outcome is the improvement of the computing time of the data processing algorithm by reducing the number of cost function evaluations. The KCD algorithm is adapted and used to detect the flawless area in a conductive plate using eddy-current testing signals, more precisely impedance measurements. The forward model linking measurements to electromagnetic properties of the plate is solved using a finite element method (FEM). The use of a defect localization step leads to reducing the number of measurement data that will be processed in order to solve the inverse problem of defect characterization. Simulation results show the efficiency of the proposed adaptation of the KCD algorithm in terms of defect localization and computing time reduction for the estimation of its geometrical parameters.

INTRODUCTION

Non-destructive testing (NDT) aims to inspect materials integrity without destroying the functionality of the component or the system [START_REF] Rivenez | Les contrôles non destructifs[END_REF]. There are several sensing NDT techniques based on different measuring modalities and signal processing methods. The focus of this paper is on Eddy Current testing (ECT) that allows treating surface flaws or relatively close to the surface of the material [START_REF] Boeing | Eddy Current Familiarization Course[END_REF]. The detection in this technique is based on the fact that the electromagnetic parameters of the material are disrupted by the presence of the defect. To solve the defect characterization problem, two main problems should be treated: the resolution of the forward problem, which consists in calculating the electromagnetic field knowing all geometric and electromagnetic parameters of the diagnosed surface, and the inverse problem which consists in using the measured impedance signal to deduce the geometric parameters (location, size and shape) and/or electromagnetic parameters (mainly permittivity and conductivity) of the defect.

Non-destructive testing techniques are generally carried out in three steps: detection, location, and characterization. Detection aims to recognize the presence of a defect, location permits to determine the damaged area, characterization step aims at finding geometrical dimensions of the defect. Actually, the localization step is very useful for the reduction of the computation time by avoiding the treatment of the whole collected data for the characterization step. Actually, each data point necessitates the solving of the forward problem and using the whole data set will induce a prohibitive computing time. The resolution of the inverse problem of defect localization and characterization in Eddy-current testing is a longstanding challenge and several methods have been proposed in a number of papers [START_REF] Bernieri | An SVM Approach to Crack Shape Reconstruction in Eddy Current Testing[END_REF][START_REF] Bernieri | Crack Shape Reconstruction in Eddy Current Testing Using Machine Learning Systems for Regression[END_REF][START_REF] Douvenot | Bases de données séquentielles et méthodes inverses appliquées au CND par courants de Foucault" Réunion plénière " Interférences d'Ondes[END_REF][START_REF] Lee | An Inverse Analysis for Crack Identification in Eddy Current NDT of Tubes[END_REF].

Explain the principle of these methods, their drawbacks and limitations and why there is a need to develop additional methods

The most issue encountered in the inversion process is the expensive computing time due to the electromagnetic model simulation in the forward problem. In this paper our objective is to reduce the calculation time by using a kernel detection method KCD [START_REF] Desobry | An Online Kernel Change Detection Algorithm[END_REF], to locate the damaged area and to apply the defect characterization algorithm to the data acquired in the neighborhood of the detected defect. In other words, the sensing data are selected so as it will be easier to estimate defects characteristics using a global optimization method such as particle swarm optimization [START_REF] Hamel | Imperialist Competitive Algorithm and Particle Swarm Optimization Comparison for Eddy Current Non-destructive Evaluation[END_REF][START_REF] Cacciola | Swarm Optimization for Imaging of Corrosion by Impedance Measurements in Eddy Current Test[END_REF].

However, the application of the KCD algorithm for the defect localization is not straightforward since the impedance variation does not present abrupt changes. So, an improvement of the detection index is presented and an optimization criterion is adopted to maximize the sparsity of this detection function.

The rest of this paper is organized as follows. We firstly present, in Section 2, the forward problem, its resolution using a Finite Element Method (FEM) and its numerical simulation for impedance calculation. Then the inverse problem of defect localization and the proposed application of the KCD algorithm for the defect localization in the context of EC-NDT are discussed in Section 3, where an adaptation is proposed in order to improve its efficiency. Section 4 gives some numerical results illustrating the usefulness of the proposed approach. Finally conclusion is drawn in Section 5.

THE FORWARD PROBLEM

We consider an axisymmetric non-destructive testing problem where a conductive plate is affected by two internal cracks having different geometries. We also assume that the defects have a rectangular shape and are characterized by their Positions (𝑃 ! = 10𝑚𝑚, 𝑃 ! = 25𝑚𝑚 ), Depths (𝐷 ! = 0.85𝑚𝑚,𝐷 ! = 0.508𝑚𝑚) and Lengths (𝐿 ! = 5𝑚𝑚 ,𝐿 ! = 3𝑚𝑚). The sensor is a differential probe moving along the plate (see Figure. Figure 1 shows the retained sensing configuration where the differential probe constituted by two coils allow measuring the impedance variation along the probed plate. The geometrical and the electromagnetic properties of the coils and the plate are summarized in Table 1. According the Eddy current sensing principle [START_REF] Bernieri | Crack Shape Reconstruction in Eddy Current Testing Using Machine Learning Systems for Regression[END_REF], the axisymmetric diffusion equation can be expressed in variational terms by the following energy function

Coil

𝐹 = ! !! ! ! ! !.! !" ! + !" !" ! + !"# ! 𝐴 ! -𝐽 ! 𝐴 ! 2𝜋𝑟𝑑𝑟𝑑𝑧 (1)
Where the three energy terms of the integrand appearing in equation ( 1) correspond respectively to the magnetic field, Eddy current and source current.

A Finite Element (FE) method permits to compute the magnetic vector potential 𝐴 which allows to calculate coil impedance.

The impedance of a circular filament of radius 𝑟 ! is calculated from the magnetic vector potential 𝐴 ! at 𝑟 ! and the value of current injected into the coil 𝐼 ! as follows:

𝑍 ! = - ! ! ! = !" !" ! ! = !".!!.! ! .! ! ! ! (2)
The total impedance expression of circular coil whose section is discretized in 𝑁 triangular elements is given by:

𝑍 = !".!!.! ! ! ! 𝑟 !" ! !!! . 𝐴 !" Δ ! = !".!!.! ! ! ! ! ( ! !!! 𝑟 !" Δ ! ). A !" (3)
Where: 𝑁 ! : Turns density [turns/m²] Δ ! : Surface of the i-th element. 𝑟 !" , 𝐴 !" : Central values of 𝑟, 𝐴 in the i-th element [START_REF] Bernieri | Crack Shape Reconstruction in Eddy Current Testing Using Machine Learning Systems for Regression[END_REF].

LOCALIZATION ALGORITHM

The differential probe is moving along the plate to measure variations of the coil impedance Δ𝑍 ! at each point 𝑘 = 1, … , 𝑛, the total impedance of differential probe can be obtained by summing the impedance of each coil.

To solve inverse problem of defect localization we use an abrupt change detection algorithm which is performed to detect damaged zones from this impedance measurements.

The aim is to locate defects from this set of measured impedance variations, which is considered as the input data of the algorithm KCD. The on-line abrupt change detection problem is described as follows:

Let 𝑡 be some time instant and assume that the samples in𝑥 !,! (resp. in𝑥 !,! ) are sampled i.i.d. according to some probability density function (pdf) 𝑝 ! (resp. 𝑝 ! ) Then, one of the two hypotheses holds:

𝐻 ! : 𝑝 ! = 𝑝 ! , (No abrupt change occurs) 𝐻 ! : 𝑝 ! ≠ 𝑝 ! ,
An abrupt change occurs .

This test cannot be applied in practice, however, since the pdfs 𝑝 ! and 𝑝 ! are not known. The standard practical approach uses some dissimilarity measure between 𝑝 ! and 𝑝 ! , estimated from the sole knowledge of the sets 𝑥 !,! and 𝑥 !,! . Let 𝐷 𝑥 !,! , 𝑥 !,! be such a dissimilarity measure, the previous problem can then be written as follows:

𝐻 ! : 𝐷 𝑥 !,! , 𝑥 !,! ≤ 𝜂, (No abrupt change occurs) 𝐻 ! : 𝐷 𝑥 !,! , 𝑥 !,! > 𝜂, An abrupt change occurs . ( 5
)
Where 𝜂, is a threshold that tunes the sensibility/robustness tradeoff, as in every detection framework.

The general framework for abrupt change detection requires the dissimilarity measure 𝐷 . , . aimed at comparing the sets of descriptors 𝑥 !,! and 𝑥 !,! . More precisely, a relevant 𝐷 . , . should give low values whenever 𝑥 !,! and 𝑥 !,! occupy the same region of the space 𝜒 and large values whenever 𝑥 !,! and 𝑥 !,! occupy the distinct regions.

The Kernel Change Detection algorithm was originally presented in [START_REF] Lee | An Inverse Analysis for Crack Identification in Eddy Current NDT of Tubes[END_REF]. It belongs to the family of abrupt change detection methods built on single-class 𝜈-Support vector machines.

Considering the analysis time 𝑡. Assume that two single-class classifiers are trained independently on the sets 𝑥 !,! and 𝑥 !,! , yielding two regions 𝑅 ! !,! ! and 𝑅 ! !,! ! or, equivalently in feature space 𝐻, two hyperplanes 𝑊 !,! and 𝑊 !,! parameterized by (𝑤 !,! , 𝜌 !,! ) and (𝑤 !,! , 𝜌 !,! ) .

In 𝐻, the vectors 𝑤 !,! and 𝑤 !,! define a 2-dimensional plane, denoted 𝑃 ! , that intersects the hypersphere 𝑆 along a circle with center 𝑂 and radius 1, as depicted in Fig. 3. Actually, in the least probable case where 𝑤 !,! and 𝑤 !,! are collinear (which is highly unlikely), there is an infinity of planes 𝑃 ! , and one can select any of them. In feature space, the plane 𝑊 !,! (resp. 𝑊 !,! ) bounds the segment of 𝑆 where most of the mapped points in 𝑥 !,! (𝑟𝑒𝑠𝑝. 𝑥 !,! ) lie. The dissimilarity measure is defined as [START_REF] Desobry | An Online Kernel Change Detection Algorithm[END_REF] 𝐷

! 𝑥 !,! , !,! = ! !,! ,! !,! ! !,! ,! !,! !! !,! ,! !,! = SIMULATION RESULTS
The KCD algorithm was applied to a simulated data according to the sensing configuration illustrated in Figure 1. The absolute value of impedance variation signal Δ(𝑍) is used as the input data for the solving of the defect localization problem.

The design parameters of the KCD algorithm are summarized in Figure 4 shows the input data and the resulting KCD index. One can see in Figure 4 (top) that changes in Δ(𝑍) signal are not abrupt; contrary they present a slow variation. Hence classical KCD algorithm could not detect real changes (see Figure 4 (bottom)).

Adaptation of the KCD algorithm

In order to account for the slow variation of the input signal, the KCD algorithm is modified so as to compare more distant frames of the input impedance variation (see figure 5). We introduce a horizon equal to 2ℎ between the two distant frames of the signal:

𝑥 !,! = 𝑥 !!!!! , … , 𝑥 !!!!! and 𝑥 !,! = 𝑥 !!!!! , … , 𝑥 !!!!! .
Such modification allows comparing different parts of the signal, which will exhibit more abrupt changes than in the case of successive frames. To validate the proposed modification of the KCD algorithm we compute a novel KCD index for both ℎ = 5 and ℎ = 10 (see Figure 6). We notice an improvement when we introduce horizon h. However, the choice of the value of the horizon h should be set to optimize the algorithm performances.

Optimization of the detection index

Since the KCD index will be used to detect the location of the abrupt changes, it should therefore be zero (or take small values) unless when the change occurs. Such property can be measured using the concept of sparsity. In that respect, the best value of ℎ should lead to the sparsest index ℎ = argmin ! 𝑆(𝐼).

Where the sparsity measure is defined by [START_REF] Hoyer | Non-negative Matrix Factorization with Sparseness Constraints[END_REF] 𝑆(𝐼)

= ! ! (!) ! !!! (! ! (!) ! !!! ) ! . ( 8
)
Figure 6 provides the evolution of the KCD index with respect to the value of h. It can also be noted that the best value is h=20. We present thereafter kernel change detection index for h=20 (figure 8). Results shown in figure 8 (bottom) prove the efficiency of the proposed method for detecting slow variations of the impedance. We can extract the beginnings and the ends of flawless zones according to the peaks that show KCD index with a threshold value 𝜂 = 20.

Optimization results

We apply particle swarm algorithm [START_REF] Desobry | An Online Kernel Change Detection Algorithm[END_REF] to defect characterization in order to estimate position, depth and length of two axisymmetric grooves. Particle swarm optimization is a heuristic global optimization method. It was developed in 1995 by Russel Eberhart and James Kennedy [START_REF] Kennedy | Particle Swarm Optimization[END_REF]. It can be used to explore the search space of any problem to find the set of parameters that minimize/maximize an objective. In PSO the position of each particle in the swarm is updated in order to perform its next movement by adjusting its actual velocity 𝑣 ! (𝑡) at t, its best performance 𝑃𝑏 ! and the best performance in the swarm 𝑃𝑔 ! according to these following equations [START_REF] Eberhart | A New Optimizer Using Particle Swarm Theory[END_REF][START_REF] Shi | A Modified Particle Swarm Optimizer[END_REF]:

𝑣 ! 𝑡 + 1 = 𝑤𝑣 ! 𝑡 + 𝑐 ! 𝑟 ! 𝑃𝑏 ! -𝑥 ! 𝑡 + 𝑐 ! 𝑟 ! 𝑃𝑔 ! -𝑥 ! 𝑡 (9)
𝑥 ! 𝑡 + 1 = 𝑥 ! + 𝑣 ! 𝑡 + 1 (10)
Where: 𝑣: particle velocity, 𝑥: particle position , 𝑟 ! , 𝑟 ! : random numbers generated from the interval [0,1], 𝑐 ! , 𝑐 ! : intensities of attraction towards 𝑃 ! and 𝑃 ! respectively, 𝑤: inertia factor (to control velocity). The particle swarm optimization algorithm parameters are the same as those used in [START_REF] Hamel | Imperialist Competitive Algorithm and Particle Swarm Optimization Comparison for Eddy Current Non-destructive Evaluation[END_REF].

PSO coefficient Value 𝑐 ! , 𝑐 ! 1.4 𝑤 !"# 0.9 𝑤 !"# 0.4 Particle swarm optimization method applied to characterization of the two defects shows good results with a very weak and acceptable uncertainty.

CONCLUSION

A novel KCD approach has been proposed for detecting slow variations in Eddy current nondestructive testing (EC-NDT). We first solved the forward problem where a plate is affected by two internal grooves with finite element method FEM and compute impedance variation that we use after as an input for the KCD algorithm, to extract the damaged zones in the plate from impedance variations signal. The application of the proposed KCD algorithm to Eddy current testing signal showed good results. A couple of numerical examples show the performances of this proposed method which is considered as a first step for inverse problem solution in order to reduce time computation and a number of function evaluations for flaw characterization in optimization step.

  1).

Figure 1 :

 1 Figure 1: Schematic drawing of the geometry.

3. 1

 1 Kernel change detection algorithm Most of on-line detection techniques are special instances of the following generic framework: Consider time 𝑡 and two descriptor subsets (the immediate past subset 𝑥 !,! = 𝑥 ! !!!!! ! ,…,!!! and the immediate future 𝑥 !,! = 𝑥 ! !!!,…,!!! ! !! , as represented in figure 3.

Figure 2 :

 2 Figure 2: General abrupt change detection framework based on the time series of descriptors 𝑥 ! , 𝑡 = 1,2, … , 𝑛 represented by circles.

Figure 3 :

 3 Figure 3: The SV single-class classifiers yield two regions regions 𝑅 𝑥 𝑡,1 𝐻 and 𝑅 𝑥 𝑡,2 𝐻 which are density support estimates in feature space. The circle represented corresponds to the intersection of the plane 𝑃 𝑡 (uniquely defi ned by 𝑤 𝑡,1 and 𝑤 𝑡,2 ) and 𝑆. The intersection of the (prolongated) vector 𝑤 𝑡,1 (resp. 𝑤 𝑡,2 ) with 𝑆 yields 𝑐 𝑡,1 (resp. 𝑐 𝑡,2 ), and the intersection of the hyperplane 𝑊 𝑡,1 ( resp 𝑊 𝑡,2 ) with 𝑆 in the plane 𝑃 𝑡 yields two points, one of which is denoted 𝜌 𝑡,1 (resp. 𝜌 𝑡,2 ). The situation plotted corresponds to an abrupt change, as both regions do not strongly overlap.

Figure 4 :

 4 Figure 4: Δ 𝑍 signal (top), Index KCD associated (bottom).

Figure 5 :

 5 Figure 5: Framework adopted for detecting ruptures Online (After modification).

Figure 6 :

 6 Figure 6: Δ 𝑍 signal (top), Index KCD associated for h=5 and h=10 (bottom).

Figure 7 :

 7 Figure 7: Sparsity criterion.

Figure 8 :

 8 Figure 8: Δ 𝑍 signal (top), Index KCD associated after modification (bottom).For a threshold 𝜂 = 20.

Table 1 :

 1 Coil and plate parameters.

		Descriptions
	Inner radius	7.75 mm
	Outer radius	8.5 mm
	Height along z	0.75 mm
	Distance to the	0.1 mm
	plate
	Number of turns	70
	Frequency	100 kHz
	Plate	Descriptions
	Width	1.27 mm
	Length	40 mm
	Conductivity	10 ! S/m
	Relative	1
	permeability

Table 2 :

 2 The dispersion parameter 𝜎 for the Gaussian kernel, the parameter of regularization 𝜈, the training set 𝑚 and the threshold 𝜂.

	Parameters	Values
	𝜎	10
	𝜈	0.5

Table 2 :

 2 parameters values of KCD algorithm.

Table 3 :

 3 PSO algorithm parameters.We present in Table4the performances of the PSO algorithm :

	Defects	Estimated position	Estimated length	Estimated depth
	Defect 1	10 mm	5mm	0.834mm
	Defect 2	25,03mm	2,99mm	0.509mm
		Table4: defect estimated parameters.


Olivier Ozenda

Pierre Saramito
email: pierre.saramito@imag.fr

Guillaume Chambon

Shear-induced migration in concentrated suspensions: particle mass conservation, contact pressure and jamming

Keywords: Suspension, Migration, Congested problem, Augmented Lagrangian method, Asymptotic expansion

A two-velocity fluid model is presented to describe particle migration in monodisperse suspensions of neutrally buoyant particles. In contrast to previous migration models, the proposed formulation allows us to impose explicit boundary conditions on particle velocity, and thereby to satisfy strict mass conservation for the particle phase. In addition, the upper bound on particle volume fraction (jamming limit) is strictly enforced through a non-smooth complementarity condition and the introduction of a particle jamming pressure. The model is applied to an axisymmetric Poiseuille flow and solved using a finite-element method. For that purpose, a specific, fully implicit algorithm based on nonsmooth optimisation tools is developed and validated. Preliminary comparisons with experimental data from the literature show promising agreement. In particular, the model properly captures the formation of an inner plug region, in which the suspension is saturated and jammed.

Introduction

Last decade saw the introduction of granular concepts to describe the rheology and flow of suspensions of rigid particles [START_REF] Boyer | Unifying suspension and granular rheology[END_REF][START_REF] Gallier | Rheology of sheared suspensions of rough frictional particles[END_REF]. While hydro-dynamical effects are prominent in dilute mixtures, inter-particle contacts and friction start to play a role as soon as particle volume fraction φ exceeds values of 0.2-0.25, typically [START_REF] Guazzelli | Rheology of dense granular suspensions[END_REF]. Contacts and friction, in particular, appear to be responsible for a number of specific rheological properties of concentrated suspensions, such as existence normal stress differences, particle pressure and micro structure anisotropy [START_REF] Guazzelli | Rheology of dense granular suspensions[END_REF][START_REF] Denn | Shear thickening in concentrated suspensions of smooth spheres in newtonian suspending fluids[END_REF]. Contacts and friction can also lead to shear-thickening effects, and influence the value of the critical particle volume fraction φ m above which the suspension is jammed and behaves as a solid [START_REF] Mari | Shear thickening, frictionless and frictional rheologies in non-brownian suspensions[END_REF][START_REF] Singh | A constitutive model for simple shear of dense frictional suspensions[END_REF]. Various constitutive models have been proposed to describe these properties, either through phenomenological expressions of the particle stresses [START_REF] Morris | Curvilinear flows of noncolloidal suspensions: the role of normal stresses[END_REF][START_REF] Dbouk | Shear-induced particle migration: Predictions from experimental evaluation of the particle stress tensor[END_REF][START_REF] Singh | A constitutive model for simple shear of dense frictional suspensions[END_REF], inclusion of auxiliary conformation tensors related to micro structure evolution [START_REF] Phan-Thien | A new constitutive model for monodispersed suspensions of spheres at high concentrations[END_REF][START_REF] Stickel | A constitutive model for microstructure and total stress in particulate suspensions[END_REF][START_REF] Goddard | A dissipative anisotropic fluid model for non-colloidal particle dispersions[END_REF][START_REF] Ozenda | A new rate-independent tensorial model for suspensions of non-colloidal rigid particles in newtonian fluids[END_REF][START_REF] Ozenda | Tensorial rheological model for concentrated non-colloidal suspensions: normal-stress differences[END_REF], or by explicitly accounting for granular processes [START_REF] Chauchat | A three-dimensional numerical model for incompressible two-phase flow of a granular bed submitted to a laminar shearing flow[END_REF][START_REF] Lecampion | Confined flow of suspensions modelled by a frictional rheology[END_REF].

Among the specific properties of concentrated suspensions, shear-induced particle migration received a lot of attention since the seminal study of Leighton and Acrivos [START_REF] Leighton | The shear-induced migration of particles in concentrated suspensions[END_REF]. This process is responsible for the spontaneous development of particle volume fraction heterogeneities in sheared suspensions and can lead, for sufficiently large values of the average volume fraction, to the formation of jammed plugs in which φ reaches φ m [START_REF] Lyon | An experimental study of the motion of concentrated suspensions in two-dimensional channel flow. Part 1. monodisperse systems[END_REF][START_REF] Oh | Pressure-driven suspension flow near jamming[END_REF]. While some recent approaches explore the direct numerical simulation of the fluid containing a discrete distribution of particles (see e.g. [START_REF] Wachs | Accuracy of finite volume/staggered grid distributed Lagrange multiplier/fictitious domain simulations of particulate flows[END_REF]), most numerical simulations of the migration process base on a continuous mathematical model for the fluid-particle mixture. Since the first phenomenological modelling attempt by Phillips et al. [START_REF] Phillips | A constitutive equation for concentrated suspensions that accounts for shearinduced particle migration[END_REF], various models based on two-phase mixture theory [START_REF] Jackson | Locally averaged equations of motion for a mixture of identical spherical particles and a Newtonian fluid[END_REF] have been proposed. Formulation of a closed system of conservation laws for a two-phase continuous medium requires closure assumptions to express the contributions of each phase to the Cauchy stress tensor of the mixture as well as the forces on the particle phase [START_REF] Lhuillier | Migration of rigid particles in non-Brownian viscous suspensions[END_REF][START_REF] Nott | The suspension balance model revisited[END_REF]. Specific attention should be paid to the contribution of particles and contacts to stresses, as it is now clear that shear-induced migration is driven by the existence of normal stress gradients in the suspension [START_REF] Guazzelli | Rheology of dense granular suspensions[END_REF][START_REF] Denn | Shear thickening in concentrated suspensions of smooth spheres in newtonian suspending fluids[END_REF]. The classical suspension balance model (SBM), which is based on empirical expressions for the particle stress and inter-phase drag, expresses, in absence of inertia, as a closed system of equations for the mixture velocity and particle volume fraction [START_REF] Nott | Pressure-driven flow of suspensions: simulation and theory[END_REF][START_REF] Morris | Curvilinear flows of noncolloidal suspensions: the role of normal stresses[END_REF]. Qualitatively, this model proved successful in capturing migration effects in different flow configurations [START_REF] Miller | Normal stress-driven migration and axial development in pressure-driven flow of concentrated suspensions[END_REF][START_REF] Dbouk | Shear-induced particle migration: Predictions from experimental evaluation of the particle stress tensor[END_REF]. However, difficulties arise when φ approaches the limit φ m and the strain rate vanishes. In practice, these issues are usually dealt with by considering unrealistically large values of φ m , and by adding an ad-hoc non-local term that effectively prevents the strain rate from vanishing [START_REF] Miller | Normal stress-driven migration and axial development in pressure-driven flow of concentrated suspensions[END_REF]. As a consequence the model cannot capture the formation of truly jammed plugs. Note however that a recent extension of SBM, implementing a process of inelastic compressibility through which the particle volume fraction can increase beyond the limit φ m in jammed regions, has been shown to overcome this limitation, and to effectively predict realistic plugs [START_REF] Lecampion | Confined flow of suspensions modelled by a frictional rheology[END_REF]. Physically, such solid compressibility of jammed regions can arise due to non-local effects and fluctuations induced by the neighbouring sheared regions [START_REF] Pouliquen | Fluctuating particle motion during shear induced granular compaction[END_REF][START_REF] Lecampion | Confined flow of suspensions modelled by a frictional rheology[END_REF][START_REF] Gaume | Microscopic origin of nonlocal rheology in dense granular materials[END_REF]. Finally, another drawback of such single velocity models is the difficulty to impose proper conditions on the particle flux at the boundaries of the domain. As a consequence, conservation of particle mass might not always be satisfied in numerical implementations [START_REF] Miller | Normal stress-driven migration and axial development in pressure-driven flow of concentrated suspensions[END_REF].

In this paper, we propose an alternative migration model that explicitly integrates the saturation of the mixture when the volume fraction reaches the limit φ m . A unilateral constraint is added to the system of conservation laws through a complementarity condition to ensure that φ remains effectively bounded by φ m . In addition, the model is based on a two-velocity formulation, and embodies an explicit computation of the difference between the mixture velocity and the particle phase velocity (i.e. the migration velocity). Boundary conditions on the particle flux can be explicitly imposed, such that strict conservation of particle mass is satisfied. Since the study is primarily devoted to the formulation and numerical solution of the system of conservation laws, particle stress is modelled using the classical phenomenological constitutive relation proposed by Morris and Boulay [START_REF] Morris | Curvilinear flows of noncolloidal suspensions: the role of normal stresses[END_REF]. Inter-particle friction or effects of micro structure are not directly taken into account at this stage, although they might be included in further versions of the model by using more sophisticated constitutive models [START_REF] Stickel | A constitutive model for microstructure and total stress in particulate suspensions[END_REF][START_REF] Lecampion | Confined flow of suspensions modelled by a frictional rheology[END_REF][START_REF] Ozenda | Tensorial rheological model for concentrated non-colloidal suspensions: normal-stress differences[END_REF]. As explained below, we do however introduce an additional pressure term that can be interpreted as a jamming pressure, to properly deal with the complementarity condition. Compared to SBM, our model also involves an additional inter-phase stress term required for the stability of the two-velocity formulation and to enforce boundary conditions on particle velocity.

From a mathematical standpoint, flow models involving a complementarity condition are usually referred to as congested problems [START_REF] Bresch | Singular limit of a Navier-Stokes system leading to a free/congested zones two-phase model[END_REF][START_REF] Bresch | Development of congestion in compressible flow with singular pressure[END_REF]. The mathematical properties of such systems were first studied by Lions and Masmoudi [START_REF] Lions | On a free boundary barotropic model[END_REF], and more recently by Bresch et al. [START_REF] Bresch | Singular limit of a Navier-Stokes system leading to a free/congested zones two-phase model[END_REF]. Hyperbolic variants of these models, without diffusive terms, were first applied to road traffic [START_REF] Berthelin | A model for the formation and evolution of traffic jams[END_REF][START_REF] Berthelin | A model for the evolution of traffic jams in multilane[END_REF] using the asymptotic preserving numerical method [START_REF] Degond | Numerical simulations of the euler system with congestion constraint[END_REF]. Later, these models were also applied to crowd motion [START_REF] Faure | Crowd motion from the granular standpoint[END_REF], granular media [START_REF] Perrin | Free/congested two-phase model from weak solutions to multi-dimensional compressible Navier-Stokes equations[END_REF][START_REF] Bresch | Compression effects in heterogeneous media[END_REF], and shallow flows in cavities [START_REF] Godlewski | Congested shallow water model: roof modeling in free surface flow[END_REF]. A general overview of hyperbolic systems with unilateral constraints was exposed by Bouchut et al. [START_REF] Bouchut | A hierarchy of models for two-phase flows[END_REF], and specific solutions were studied by Berthelin [START_REF] Berthelin | Existence and weak stability for a pressureless model with unilateral constraint[END_REF], Berthelin and Bouchut [START_REF] Berthelin | Weak solutions for a hyperbolic system with unilateral constraint and mass loss[END_REF] and recently Chen and Zhai [START_REF] Chen | Global large solutions and incompressible limit for the compressible Navier-Stokes equations[END_REF].

The paper is structured as follows. The proposed two-phase model for the mixture velocity, migration velocity and particle volume fraction is presented in section 2. Through an asymptotic analysis, a reduced system is then derived for the case of a uniform axisymmetric Poiseuille flow (flow in a circular tube). Section 3 proposes a numerical method to solve the model in the case of the Poiseuille problem, implementing a specific augmented Lagrangian approach to handle the nonlinearity associated with the complementarity condition. Spatial discretisation is performed using finite elements. Results are presented and discussed in section 4. After a careful investigation of the convergence properties of the numerical method, the physical characteristics of the solutions are described, and direct comparisons with experimental measurements of Oh et al. [START_REF] Oh | Pressure-driven suspension flow near jamming[END_REF]

Migration model

Two velocity formulation

Let r p be the radius of the rigid spherical particles, and η 0 be the viscosity of the interstitial Newtonian fluid. For the sake of simplicity, we assume here a neutral buoyancy: let ρ denote the constant mass density of both the fluid and the particles. This assumption is not fundamental, and the present theory could be extended to also include different mass densities and sedimentation effects.

The dynamics of the mixture, at a continuous macroscopic scale, is described by two independent velocities. Without loss of generality, we choose as independent variables the velocity of the mixture, denoted by u, and the velocity difference between the two phases, hereafter called migration velocity, denoted by w. Note that u + w represents the velocity of the particles phase.

The volume fraction is denoted by φ and is bounded by the maximal volume fraction φ m . When φ = φ m , the mixture is jammed and behaves as a solid.

For convenience, the reduced volume fraction ψ = φ/φ m is introduced. The constraint φ φ m , or equivalently ψ 1, can be expressed as a linear complementarity problem Cottle and Dantzig [START_REF] Cottle | Complementary pivot theory of mathematical programming[END_REF], Duvaut and Lions [START_REF] Duvaut | Les Inéquations En Mécanique Et Physique[END_REF]:

0 (1 -ψ) ⊥ p j 0 ⇔    1 -ψ 0 p j 0 (1 -ψ)p j = 0 (1)
Hence, the two quantities 1 -ψ and p j should always be positive, and the ⊥ notation expresses that their product should also be zero. When ψ < 1, we then have p j = 0, while p j can be nonzero when ψ = 1 in the jammed case. This constraint expresses, at the macroscopic scale, the microscopic non-penetration between the rigid particles when contacts occur. The quantity p j is a Lagrange multiplier that can be interpreted as a jamming pressure. Note that (1) should be satisfied locally, as some regions of the flow can be jammed while others are not.

Table 1 summarises the main notations used in this paper. Let Ω ⊂ R d denote the flow domain, where d 1 is the physical space dimension, and let T > 0 be the final time. The problem to solve writes: (P ): find u, w, p, ψ, p j , defined in ]0, T [×Ω, such that

                                                                       ρ ∂u ∂t + u.∇u -div (-pI + 2η 0 D(u) + τ p ) = f in ]0, T [×Ω ρ φ m ψ ∂(u+w) ∂t +(u+w).∇(u+w) -div (-p j I +τ p ) + η 0 s(ψ) r 2 p w -div (2η 0 s(ψ)D(w)) = 0 in ]0, T [×Ω div u = 0 in ]0, T [×Ω ∂ψ ∂t + div ((u + w)ψ) = 0 in ]0, T [×Ω 0 (1 -ψ) ⊥ p j 0 in ]0, T [×Ω τ p = 2η 0 ψ 1 -ψ 5φ m 2 + ψ 1 -ψ K s D (u) -2η 0 ψ 1 -ψ 2 K n |D (u)| Q in ]0, T [×Ω u = u Γ and w = w Γ in ]0, T [×∂Ω ψ = ψ Γ in ]0, T [×∂Ω - u(t = 0) = u 0 , w(t = 0) = w 0 and ψ(t = 0) = ψ 0 in Ω (2a) (2b) (2c) (2d) (2e) (2f) (2g) (2h) (2i)
Equations (2a) and (2c) express the momentum and mass conservation of the mixture, respectively, while Equations (2b) and (2d) express the momentum and mass conservation of the particle phase. Note that the latter has the form of an evolution equation for the reduced volume fraction ψ. Equation (2a) involves the expression of the Cauchy stress tensor of the mixture σ m = -pI + 2η 0 D(u) + τ p , where D(u) = (∇u + ∇u T )/2. The Lagrange multiplier p, that interprets as the bulk pressure, is introduced to enforce the mixture incompressibility constraint (2c). The quantity τ p represents the contribution of the particle phase to the extra-stress tensor, and will be discussed later.

On the right-hand-side, f represents any external body force. Equation (2b) involves the Cauchy stress of the particle phase -p j I + τ p , where p j is the jamming pressure related to the jamming constraint (2e), as explained above. Note that the particle pressure then expresses as p j -tr(τ p )/3. Equation (2b) also involves a net force exerted on the particle phase:

η 0 s(ψ) r 2 p w -div (2η 0 s(ψ)D(w)) .
The first term in this force corresponds to the drag force exerted by the fluid phase, where s(ψ) is a hindrance function for which we use the following expression proposed by Miller and Morris [START_REF] Miller | Normal stress-driven migration and axial development in pressure-driven flow of concentrated suspensions[END_REF]:

s(ψ) = 9 2(1 -φ m ψ) α-1 (1 -ψ) ,
with α ∈ [2, 5] a material parameter. As discussed by Miller and Morris [START_REF] Miller | Normal stress-driven migration and axial development in pressure-driven flow of concentrated suspensions[END_REF], this expression is adapted from the classical Richardson-Zaki function, with an additional (1 -ψ) factor in the denominator to enforce that migration ceases when maximum packing fraction is reached. The second term in the particlephase force, namely div (2η 0 s(ψ)D(w)) , represents a second-order correction with respect to particle radius r p . It is introduced here in order to allow for the imposition of boundary conditions on particle velocity u + w, or equivalently on w (see Equation (2g)).

Equation (2f) is the phenomenological constitutive relation for the particle stress tensor τ p initially proposed by Morris and Boulay [START_REF] Morris | Curvilinear flows of noncolloidal suspensions: the role of normal stresses[END_REF]. Constants K s and K n are material parameters that control the shear and normal viscosities, respectively. Note that the particle stress tensor τ p diverges when ψ → 1, so that this constitutive model cannot be used as such in saturated regions. In what follows, a regularisation will be used to overcome this limitation (see section 3). The tensor Q expresses as diag(1, λ 2 , λ 3 ) in the velocity, gradient, vorticity basis associated to viscosimetric flows (see e.g. [44, p. 158]), where λ 2 and λ 3 are material parameters controlling the two normal stress differences. Expressions of Q adapted to the case of non-viscosimetric flows have also been proposed [START_REF] Miller | Suspension flow modeling for general geometries[END_REF], but will not be used here.

Finally, the problem is closed by suitable initial and boundary conditions on velocities u and w and reduced volume fraction ψ, expressed by Equations (2g)-(2i). Here ∂Ω -denotes the upstream boundary domain for the particle phase:

∂Ω -= {x ∈ ∂Ω such that (u + w)(x) • n(x) < 0} ,
and u Γ , w Γ , ψ Γ , u 0 , w 0 , ψ 0 , are given boundary and initial data satisfying div u 0 = 0, ∂Ω u Γ ds = 0 and ψ 0 , ψ Γ ∈ [0, 1]. Observe that, from (2d) and since ψ(t = 0) 0, then ψ(t) 0 for any t 0.

Suspension balance model (SBM), which only involves the mixture velocity [START_REF] Miller | Normal stress-driven migration and axial development in pressure-driven flow of concentrated suspensions[END_REF], can be recovered by neglecting inertial terms in equations (2a) and (2b), and omitting the complementarity condition (2e) as well as the term div (-p j I + 2η 0 s(ψ)D(w)) in equation (2b). Note that, in SBM, the constraint ψ 1 is not strictly imposed, but effectively verified through the addition of an ad-hoc non-local term to the particle stress tensor τ p .

Appendix A further expands on the link between the present model and the general mixture theory developed by Jackson [START_REF] Jackson | The Dynamics Of Fluidized Particles[END_REF], including formal definitions of the variables in terms of fluid and solid phase averages. Let us also note that, as an alternative to the stress partition used in system (2), one could define the total Cauchy stress tensor of the particle phase as σ p = τ p + τ p,nl -p j I, with τ p,nl = 2η 0 s(ψ)D(w). Physically, the new corrective term τ p,nl could then be interpreted as a non-local contribution to the particle stress in unjammed regions, in the same spirit as non-local rheological models proposed for granular materials and other complex fluids [START_REF] Bouzid | Nonlocal rheology in dense granular flows. revisiting the concept of fluidity[END_REF][START_REF] Kamrin | Non-locality in granular flow: Phenomenology and modeling approaches[END_REF]. The total particle pressure then expresses as p p = p j -tr(τ p + τ p,nl )/3, also involving a contribution from τ p,nl .

Finally, let us comment on the mathematical structure of problem (2a)-(2i).

The pair (u, p) satisfies an incompressible Navier-Stokes-like subsystem (2a), (2c), while the triplet (w, ψ, p j ) satisfies a variable-density Navier-Stokes-like subsystem (2b), (2d), associated to condition (2e) that guarantees ψ ∈ [0, 1].

Coupling between the two subsystems is achieved by the tensor τ p that depends on u and ψ from (2f). Let us consider here the circular tube geometry represented on Fig. 1, with L the length of the tube and R its radius. Let (r, θ, z) be the associated cylindrical coordinate system. We consider axisymmetric flows independent upon θ. The tube is assumed to be sufficiently long, i.e. L → ∞, such that the flow is also considered to be independent upon z. Hence, mixture velocity writes u(t, r) = (0, 0, u z (t, r)), while migration velocity can develop a nonzero radial component and expresses as w(t, r) = (w r (t, r), 0, w z (t, r)). The average reduced particle volume fraction ψ 0 is constant for t 0, and we also consider that the mixture flow rate q is imposed and constant for t > 0.

Uniform Poiseuille flow

An asymptotic analysis for L → ∞, presented in Appendix B, shows that problem (2a)-(2i) then reduces to, in dimensionless form:

(Q): find u z , w z , w r , ψ, p j , defined in ]0, T [×]0, 1[ and f z in ]0, T [, s.t.                                                                            Re ∂ t u z - ∂ r r (r η app (ψ) ∂ r u z ) + f z = 0 in ]0, T [×]0, 1[ Re φ m ψ ∂ t (u z + w z ) + s(ψ) ε 2 w z - ∂ r r (r s(ψ) ∂ r w z + r (η app (ψ) -1) ∂ r u z ) = 0 in ]0, T [×]0, 1[ Re φ m ψ ∂ t w r + s(ψ) ε 2 w r - ∂ r r (2 r s(ψ) ∂ r w r ) + η n,θ (ψ) |∂ r u z | r - ∂ r r (r η n,r (ψ) |∂ r u z |) + ∂ r p j = 0 in ]0, T [×]0, 1[ ∂ t ψ + ∂ r r (rw r ψ) = 0 in ]0, T [×]0, 1[ 0 (1 -ψ) ⊥ p j 0 in ]0, T [×]0, 1[ 1 0 u z (t, r) r dr = 1 4 , ∀ t ∈ ]0, T [ ∂ r u z (t, r = 0) = u z (t, r = 1) = 0, ∀ t ∈ ]0, T [ ∂ r w z (t, r = 0) = w z (t, r = 1) = 0, ∀ t ∈ ]0, T [ w r (t, r = 0) = w r (t, r = 1) = 0, ∀ t ∈ ]0, T [ u z (t = 0) = w z (t = 0) = w r (t = 0) = 0, ψ(t = 0) = ψ 0 in ]0, 1[ (3a) (3b) (3c) (3d) (3e) (3f) (3g) (3h) (3i) (3j)
where ε = r p /R and Re = ρU R/η 0 ,, with U = 2q/(πR 2 ) the characteristic velocity scale. Here, p j denotes the Lagrange multiplier associated to the complementarity condition (3e). Similarly, f z is a Lagrange multiplier associated to the mixture flow rate constraint (3f). Note that f z interprets as the longitudinal gradient of mixture pressure. The rheological functions η app , η n,r and η n,θ involved are given by

η app (ψ) = 1 + 5φ m 2 ψ 1 -ψ + K s ψ 1 -ψ 2 η n,r (ψ) = -α n,r ψ 1 -ψ 2 η n,θ (ψ) = -α n,θ ψ 1 -ψ 2 with α n,r = λ 2 K n and α n,θ = λ 3 K n .
With these notations, the dimensionless particle stress tensor expresses as:

τ p =   η n,r (ψ) |∂ r u z | 0 (η app (ψ)-1) ∂ r u z 0 η n,θ (ψ) |∂ r u z | 0 (η app (ψ)-1) ∂ r u z 0 η n,z (ψ) |∂ r u z |   (4)
Finally, note that, unlike for 2D plane channel flow [START_REF] Lecampion | Confined flow of suspensions modelled by a frictional rheology[END_REF], particle normal stress differences are explicitly involved in the circular tube flow considered here.

Let us observe the mathematical structure of problem (3a)-(3j). Relations (3a) and (3f) constitute a linear constrained parabolic sub-problem for the unknown u z and the Lagrange multiplier f z , where ψ is considered as known. This first linear sub-problem is closed by boundary and initial conditions (3g) and (3j) for u z . Similarly, relation (3b) constitute a second linear parabolic sub-problem for the unknown w z , where both u z and ψ are considered as known. This second linear sub-problem is closed with boundary and initial conditions (3h) and (3j) for w z . Finally, relations (3c), (3d) and (3e) constitute a nonlinear constrained sub-problem for the unknowns w r and ψ and the Lagrange multiplier p j , where u z is considered as known. This problem will henceforth be called the congested flow sub-problem. This third nonlinear sub-problem is closed with boundary and initial conditions (3i) and (3j) for w r and ψ. The Lagrange multiplier p j that imposes the nonlinear constraint ψ 1 in (3e) acts on w r via its evolution equation (3c). Then, ψ is convected by w r in accordance with the mass conservation equation (3d). Hence, p j acts as a subtle indirect control upon ψ via w r . The numerical solution of the problem presented in the following section is mainly suggested by observation of this mathematical structure.

Numerical resolution

We present in this section a fully implicit algorithm for the numerical resolution of problem (Q) (uniform Poiseuille flow). At each time step of an outer loop, the three sub-problems outlined above are solved, and a fixed point inner loop ensures the convergence. Note that, unlike Degond et al. [START_REF] Degond | Numerical simulations of the euler system with congestion constraint[END_REF] and Degond and Tang [START_REF] Degond | All speed scheme for the low mach number limit of the isentropic euler equations[END_REF] who solved an hyperbolic congested flow problem with an explicit time scheme, we choose here an implicit time discretisation to avoid restrictions on time steps due to stability criteria. As the present problem also involves viscous and diffusion terms, such conditions on time step would be too restrictive.

Particle migration tends to be a slow process, for which the use of large time steps is required.

Presentation of the numerical algorithm is organised as follows. Subsection 3.1 describes the fixed point method that splits problem (Q) into the three associated sub-problems. Two of these sub-problems, namely (S 1 ) and (S 2 ), are linear and provide computations of (u z , f z ) and w z , respectively. The three other unknowns, w r , p j , ψ, are solution of a non smooth optimisation problem (S 3 ). Problem (S 3 ), which implements a complementarity condition between p j and ψ, cannot be solved by classical optimisation algorithms. This problem is thus approached by another problem ( S3 ), implementing a complementarity condition between p j and ∂ r w r /r. An optimisation method to solve problem ( S3 ) is presented in subsection 3.2. A corrective term is added to the cost function to apply an augmented Lagrangian method. The solution of ( S3 ) then expresses as a critical point of the augmented cost function, and is found using an Uzawa algorithm. Finally, subsection 3.3 describes the spatial discretisation used.

Implicit time discretisation

Let ∆t > 0 be the time step and t n = n∆t, n ∈ N, be the discrete times. A sequence of semi-discrete in time solutions u

(n) z , w (n) z , w (n) r , ψ (n) , f (n) z , p (n) j n 0
is defined by recurrence.

When n = 0, the solution is provided by the initial conditions, u

(0) z = w (0) z = w (0) r = p (0) j = 0, ψ (0) = ψ 0 and f (0) z = 0 i.e.
the material is at rest and homogeneous and ψ 0 ∈ [0, 1[ is the given initial reduced volume fraction.

When n 1, let us assume by recurrence that the numerical solution of system (Q) is given at time

t n-1 , i.e. u (n-1) z , w (n-1) z , w (n-1) r , ψ (n-1) , f (n-1) z , p (n-1) j is known. Then, u (n) z , w (n) z , w (n) r , ψ (n) , f (n) z , p (n) j
is computed using the following fixed point procedure.

Let k ∈ N denotes the index of the fixed point inner loop. At each step n 1, a sequence u

(n,k) z , w (n,k) z , w (n,k) r , ψ (n,k) , f (n) z , p (n) j k 0
is defined by recurrence.

When k = 0, this inner loop is initialised from the values at the previous time step, i.e.

u (n,0) z , w (n,0) z , w (n,0) r , ψ (n,0) , f (n,0) z i, , p (n,0) j = u (n-1) z , w (n-1) z , w (n-1) r , ψ (n-1) , f (n-1) z , p (n-1) j . When k 1, let us assume by recurrence that u (n,k-1) z , w (n,k-1) z , w (n,k-1) r , ψ (n,k-1) , f (n,k-1) z , p (n,k-1) j is known. Then, u (n,k) z , w (n,k) z , w (n,k) r , ψ (n,k) , f (n,k) z , p (n,k) j
is defined by splitting (Q) into three subsystems, two of them are linear and will be referred as (S1) and (S 2 ), the third one is a congested nonlinear problem, referred as (S 3 ).

By introducing two numerical parameters, k max and ε f p , a stopping criterion is defined:

ψ (n,k-1) -ψ (n,k) 2 + w (n,k-1) r -w (n,k) r 2 + w (n,k-1) z -w (n,k) z 2 + u (n,k-1) z -u (n,k) z 2 √ 2 ∆t ε f p or k > k max (5)
where . denotes the usual L 2 norm with axisymmetric weighting, defined for all function f by

f = 1 0 f (r) r dr 1 2
When stopping criterion ( 5) is satisfied, the fixed point loop is terminated and the last element of the sequence is simply denoted as

u (n) z , w (n) z , w (n) r , ψ (n) , f (n) z , p (n) j
, i.e. the second index k is omitted.

With the notations defined above, the two linear subsystems solved at each fixed point iteration write:

(S 1 ): find u (n,k) z defined in ]0, 1[ and f (n,k) z ∈ R such that                      Re ∆t u (n,k) z -u (n-1) z - ∂ r r r η app ψ (n,k-1) ∂ r u (n,k) z +f (n,k) z = 0 in ]0, 1[ 1 0 u (n,k) z r dr = 1 4 ∂ r u (n,k) z (r = 0) = u (n,k) z (r = 1) = 0 (6a) (6b) (6c) (S 2 ): find w (k,n) z defined in ]0, 1[ such that                      Re φ m ψ (n,k-1) ∆t w (n,k) z -w (n-1) z + ε -2 s ψ (n,k-1) w (n,k) z - ∂ r r r s ψ (n,k-1) ∂ r w (n,k) z = Re φ m ψ (n,k-1) ∆t u (n,k) z -u (n-1) z in ]0, 1[ ∂ r w (n,k) z (r = 0) = w (n,k) z (r = 1) = 0 (7a) (7b)
Similarly, the nonlinear congested subsystem write:

(S 3 ): find w (n,k) r , p (n,k) j and ψ (n,k) defined in ]0, 1[ such that                                              Re φ m ψ (n,k-1) ∆t w (n,k) r -w (n-1) r + ε -2 s ψ (n,k-1) w (n,k) r - ∂ r r 2 r s ψ (n,k-1) ∂ r w (n,k) r + ∂ r p (n,k) j = ∂ r r r η n,r ψ (n,k-1) ∂ r u (n,k) z - 1 r η n,θ ψ (n,k-1) ∂ r u (n,k) z in ]0, 1[ 0 ψ * -ψ (n,k) ⊥ p (n,k) j 0 in ]0, 1[ 1 ∆t ψ (n,k) -ψ (n-1) • X (n,k-1) + ∂ r r r w (n,k) r ψ (n,k) = 0 in ]0, 1[ w (n,k) r (r = 0) = w (n,k) r (r = 1) = 0 (8a) (8b) (8c) (8d)
where

X (n,k-1) (r) = r -∆t w (n,k-1) r
(r) in (8c) denotes a first order approximation of the characteristics.

Observe that the original unilateral constraint ψ 1 has been replaced in (8b) by ψ (n,k) ψ * , where ψ * < 1 is a numerical threshold close to 1. This threshold is used to ensure that the hindrance function s and the viscosities η app , η n,r and η n,θ remain bounded. This strategy constitutes an effective regularisation of the constitutive law in order to avoid divergence of the stress components and drag force in the saturated regions.

The two linear systems (S 1 ) and (S 2 ) are standard, while the solution of the congested nonlinear system (S 3 ) requires more work. Observe that (8c) can be explicitly solved in term of ψ (n,k) as:

1 + ∆t ∂ r r r w (n,k) r ψ (n,k) = ψ (n-1) • X (n,k-1) (9)
Indeed, for sufficiently small ∆t > 0, the first factor of the left-hand side of (8c) is strictly positive in ]0, 1[, and thus ψ (n,k) is well-defined by an explicit expression. This expression of ψ (n,k) can then be replaced in (8b). After rearrangements, we obtain a constraint in terms of w

(n,k) r : - 1 ∆t ψ * -ψ (n,k-1) • X (n,k-1) ∂ r r r w (n,k) r ⊥ p (n,k) j 0
This complementarity condition means that the compressibility of the particle phase is bounded negatively, depending on the value of the volume fraction.

Hence, congested subsystem interprets as an obstacle problem coupled to an advection equation :

( S3 ): find w (n,k) r and p (n,k) j , defined in ]0, 1[, such that                                      Re φ m ψ (n,k-1) ∆t w (n,k) r -w (n-1) r + ε -2 s ψ (n,k-1) w (n,k) r - ∂ r r 2 r s ψ (n,k-1) ∂ r w (n,k) r + ∂ r p (n,k) j = ∂ r r r η n,r ψ (n,k-1) ∂ r u (n,k) z - 1 r η n,θ ψ (n,k-1) ∂ r u (n,k) z in ]0, 1[ - 1 ∆t ψ * -ψ (n,k-1) • X (n,k-1) ∂ r r r w (n,k) r ⊥ p (n,k) j 0 in ]0, 1[ w (n,k) r (r = 0) = w (n,k) r (r = 1) = 0 (10a) (10b) (10c)
As soon as w

(n,k) r
is known, ψ (n,k) can be explicitly computed from [START_REF] Phan-Thien | A new constitutive model for monodispersed suspensions of spheres at high concentrations[END_REF].

We recognise in ( S3 ) the standard obstacle problem in mathematical physics (see e.g. [START_REF] Brezis | Sur la régularité de la solution d'inéquations elliptiques[END_REF][START_REF] Capriz | New Variational Techniques In Mathematical Physics[END_REF][START_REF] Kinderlehrer | An Introduction To Variational Inequalities And Their Applications[END_REF]). It is expressed here in term of w

(n,k) r
with a convex constraint (10b) on all the interior of the domain ]0, 1[, where p (n,k) j is the associated Lagrange multiplier.

Resolution of the non smooth obstacle sub-problem

Problem ( S3 ) defined above can be solved efficiently, without any regularisation, by an augmented Lagrangian method described in this subsection. Since there is no ambiguity, indices n and k are omitted here on the unknowns. The problem becomes:

(O): find w r and p j , defined in ]0, 1[, such that

           κw r - ∂ r r (r β ∂ r w r ) + ∂ r p j = f in ]0, 1[ g ∂ r r (r w r ) ⊥ p j 0 in ]0, 1[ w r (r = 0) = w r (r = 1) = 0 (11a) (11b) (11c)
where the following notations are introduced for the known data:

κ = Re φ m ψ (n,k-1) β = 2 s ψ (n,k-1) f = κw (n-1) r + ∂ r r r η n,r ψ (n,k-1) ∂ r u (n,k) z - 1 r η n,θ ψ (n,k-1) ∂ r u (n,k) z g = - 1 ∆t ψ * -ψ (n,k-1) • X (n,k-1)
Let L 2 , H 1 , H -1 , and H 1 0 denote the usual Hilbert functional spaces associated with the weight r for the cylindrical coordinates. Let us introduce the following convex subset of L 2 :

K = ξ ∈ L 2 ; ξ g
Note that K is indeed convex, since each convex combinations of any elements of K belongs to K. Let us also introduce the following bilinear and linear forms, defined for all w, v ∈ H 1 and q ∈ L 2 by:

a(w, v) = 1 0 (κ w v + β ∂ r w ∂ r v) r dr b(w, q) = 1 0 ∂ r (r w) q dr (v) = 1 0 f v r dr With κ, β ∈ L ∞ , and f ∈ H -1
, where L ∞ denotes the space of bounded functions. Hence, these forms are well-defined.

Moreover, we assume that g ∈ L 2 . Let B denotes the linear operator from H 1 to L 2 associated to the bilinear form b and defined for all v ∈ H 1 by

Bv = r -1 ∂ r (r v) .
It interprets as the divergence operator in the axisymmetric tube section. The quadratic function J is defined for all v ∈ H 1 by

J(v) = 1 2 a(v, v) -(v)
Problem (O) then expresses as a convex minimisation problem:

w r = arg min v∈H 1 0 J(v)
subject to Bv ∈ K Note that the convex set K is not a vector space. Thus, the previous optimisation problem is difficult to solve by finite element method, which is based on vector space approximations of functional spaces. For this reason, we introduce the indicator function

I K : L 2 → [0, ∞]
, defined for all ξ ∈ L 2 by:

I K (ξ) = 0 when ξ ∈ K ∞ otherwise
Observe that I K is a convex function since K is a convex set. The problem can then be rewritten as:

w r = arg min v∈H 1 0 J(v) + I K (Bv)
The problem now expresses as an unconstrained minimisation problem of a convex non-differentiable function on a vector space, which is more suitable to a finite element approximation. The main difficulty is to minimise with respect to I K (Bv), which is the non-differentiable part. A solution is to introduce an auxiliary variable δ, together with the additional constraint δ = Bw and its associated Lagrange multiplier, which shall coincide with the jamming pressure p j . We then introduce the following augmented Lagrangian (see e.g. [START_REF] Fortin | Augmented Lagrangian Methods[END_REF]):

L(v, ξ ; q) = J(v) + I K (ξ) + 1 0 (ξ -Bv) q r dr + µ 2 1 0 (ξ -Bv) 2 r dr (12)
where µ > 0 is the augmentation parameter. The problem is equivalent to finding the following saddle point:

(w r , δ ; p j ) = arg min v ∈ H 1 0 ξ ∈ L 2 max q∈L 2 L(v, ξ ; q)
Observe that the term factored by the augmentation parameter µ in ( 12) is the square of the constraint: the saddle-point of the Lagrangian L is thus independent of µ. The numerical parameter µ only influences the convergence of minimisation algorithm.

The solution is computed by an Uzawa descent method, fully described in Appendix C.

Finite element spatial discretisation

The dimensionless space interval [0, 1] is discretised by a uniform mesh whose step is denoted h > 0. Components u z , w z and w r of the velocities are approximated by continuous and piecewise quadratic functions, while the reduced volume fraction ψ, the jamming pressure p j and the divergence δ are approximated by continuous and piecewise linear functions. The algorithm described in the previous subsections is implemented using Rheolef C++ finite element library [START_REF] Saramito | Efficient C++ Finite Element Computing With Rheolef[END_REF]. Four different values of dimensionless mesh size h were investigated, namely h = 1/200, 1/400, 1/800, and 1/1600. The dimensionless time step ∆t is adapted according to the value of h as follows: ∆t = 400h. Hence, for each refinement of the mesh, the time step is divided by two.

Final computation time T is set such that steady-state flow regime is reached.

In practice, the time needed to reach this steady state depends on the value of average volume fraction φ 0 . For most of the numerical tests presented below (subsections 4.2, 4.3, 4.4), an arbitrary large value T = 4000 was chosen. When comparing results obtained for different values of φ 0 (subsections 4.4 and 4.5), the value of T was adapted according to a proper steady-state criterion, as explained later.

Results and discussion

This section is dedicated to a preliminary exploration of the predictions of the new migration model presented in this paper. The solutions, computed with the algorithm presented in the previous section, are compared with experimental results obtained by Oh et al. [START_REF] Oh | Pressure-driven suspension flow near jamming[END_REF]. This section starts with a presentation of the experimental setup and the choice of model material parameters. Then, validations of our numerical algorithm are presented. We start by a study of the convergence of the residual terms in the two inner loops of the implicit time discretisation scheme, and then turn to the convergence of the solution versus mesh refinement. Finally, the main physical features of the solution are described, together with preliminary comparisons to experimental data and a sensitivity analysis with respect to maximum volume fraction φ m . Oh et al. [START_REF] Oh | Pressure-driven suspension flow near jamming[END_REF] injected a mono-disperse particle suspension in a circular tube from a tank with an imposed flow rate q. At the inlet of the tube, the volume fraction is supposed to be uniform, equal to φ 0 . Particle volume fraction φ and axial velocity of the mixture u z were measured by MRI at a sufficiently long distance from the inlet to ensure fully developed flow. Note that this distance, at which the flow can be considered as fully developed, depends upon several parameters [START_REF] Oh | Pressure-driven suspension flow near jamming[END_REF]: the volume fraction at the inlet of the tube φ 0 , the density of both the fluid and the particles, the fluid viscosity η 0 , the particle radius r p , and the tube radius R. Inlet flow rate was varied between 0.5 and 3 ml/min. For the computations presented thereafter, we retained an intermediate value q = 1.88 ml/min.

The values of model constants chosen for the computations presented in this section are summarised in Table 2. The parameters K n , K s , λ 2 , λ 3 and α of the rheological model were identified from experimental data, as explained in Appendix D. The choice of the maximum volume fraction φ m will be discussed in the forthcoming paragraph dedicated to comparisons with experiments.

Convergence of the inner loops

This paragraph documents the convergence of the augmented Lagrangian loop and the fixed point loop. These numerical tests were performed with values of φ 0 = 0.32, φ m = 0.585, and an arbitrary final time T = 4000. The link between h and ∆t is provided by Table 4. norm of the auxiliary variable δ (m) during the augmented Lagrangian loop as a function of loop index m. The convergence is studied for k = 0, i.e. at the first iteration of the fixed point loop, and for t = 3T /8, i.e. fully developed flow is not reached yet. First, observe that the relative error decrease for all values of the augmentation parameter µ, as expected from theory (see Fortin and Glowinski [START_REF] Fortin | Augmented Lagrangian Methods[END_REF]). Next, observe that the convergence is faster for intermediate values of µ: the optimal value is near µ = 1280.

Let us now turn to the convergence of the fixed point loop. Fig. 2.right presents the relative error for the volume fraction ψ (n,k) versus total number of iterations k × m max of the inner loops. The initial relative error ψ (n,1) -ψ (n,0) is of about 10 -9 , and decreases to about 10 -12 , when rounding effects appear. Hence, the normalised relative error decreases to about 10 -3 . Recall that there is an inner augmented Lagrangian loop, whose index is m and maximal number of iterations is m max . To understand how the two loops interact, the value of m max has been varied. Observe that it is not necessary to iterate more than once in the augmented Lagrangian loop, for the fixed point loop to converge. Moreover, this strategy appears to be the most efficient in terms of overall convergence rate of the algorithm. Finally, the numerical parameters retained for the simulations of the next sections are grouped in Table 3. The values of µ, m max and k max ensure proper convergence of the inner loops. The stopping criterion (5) of the fixed point outer loop is set at ε f p = 10 -12 .

symbol

Spatial and temporal convergences

This paragraph is dedicated to the convergence of the solution versus the simultaneous refinement of the time and space steps. The four considered mesh configurations are summarised in Table 4. Similar to time step ∆t that tends to zero with mesh size h, we also considered that the regularisation parameter of the constitutive law ψ * (effective maximum volume fraction) evolves with mesh refinement, and tends to 1 when h → 0. Specifically, the following relation was considered based on a specific convergence study (not shown here):

ψ * = 1 -10.24h (Tab. 4).
As in the previous subsection, these numerical tests were performed for values of φ 0 = 0.32, φ m = 0.585, and T = 4000. In Fig. 3, an additional subscript h is added to variables, e.g. w r,h , ψ h , in order to indicate the mesh dependence. Fig. 3.left plots the L 2 norm of the radial migration velocity w r,h versus dimensionless time t for the four mesh refinements. Observe first that, on all meshes, w r,h decays exponentially with time: as expected, radial particle migration vanishes in steady-state regime. Furthermore, the slope of this decay appears to be mesh-independent.

An important feature of our model is that it should conserve the mass of the fluid and solid phases, as expressed by Eq. (2c) for the whole mixture and Eq. (2d) for the particle phase. At t = 0, the particle mass is given by 1 0 φ 0 r dr = φ 0 /2. Fig. 3.right shows the maximum particle mass error as a function of mesh size h. Obviously, this error is never exactly zero, as the finite element method only provides an approximation of the solution. However, the error clearly tends to zero and Fig. 3.right suggests that it converges as O(h), i.e. linearly, with mesh refinement. This linear convergence represents a major improvement compared to other existing migration models. As our formulation enables us to impose a non-penetration boundary condition (2g) on the particle velocity u + w, or equivalently on w, it prevents uncontrolled mass loss. Let us recall that this property is made possible by introducing the second-order differential term on w in (2b). Fig. 4 presents the radial profiles of various quantities at time t = T for the four mesh refinements. A clear, progressive convergence of all the profiles with h can be observed. Note that the amplitude of radial migration velocity w r is very small, confirming that the system has reached quasi-steady state. Accordingly, the bulk shear stress σ m,rz is observed to evolve linearly within the tube, as expected, for all mesh refinements. The peak of jamming pressure p j on the symmetry axis r = 0 is likely due to a numerical artefact in the treatment of the boundary condition, but appears to remain bounded as h → 0.

For each mesh size h, one can define a critical radius, denoted as r c,h (t), such that ψ h (t, r) reaches its maximum 1 -ψ * for all r ∈ [0, r c,h (t)] (see Fig. 4.top-right). This region corresponds to a central plug, where the radial migration velocity w r vanishes while the axial bulk velocity u z is constant and maximum. Table 5 indicates the values of r c,h at time t = T for the four mesh refinements (see also Tab. 4). The last column corresponds to the difference with the value obtained for the finest subdivision. Observe that this difference decreases regularly when h → 0, and suggests that the critical radius converges as O(h).

Main features of the solution

This subsection describes the main features of the solution predicted by the present migration model. The finest mesh, namely h = 1/1600, is retained for these computations. Fig. 5 shows the radial profiles of various physical quantities in the tube section at three different times. Note that values of φ 0 = 0.32,φ m = 0.585, and T = 4000 are still used here. It is observed that the central plug, in which the material is saturated (ψ = ψ * ), grows over time. At the transition between the plug and the outer sheared region (r = r c (t)), the reduced volume fraction ψ is continuous but not differentiable. As already mentioned, the radial migration velocity w r strongly decreases with time, and vanishes in the plug. Negative values of w r in the sheared region indicate that particles migrate from the wall to the centre of the tube, as expected. Interestingly, one also observes that the values of longitudinal migration velocity w z remain slightly negative even at steady-state, indicating that particle velocity lags behind fluid velocity. As expected from the complementarity condition (11b), the jamming pressure p j is nonzero inside this plug, while it is zero in the unsaturated sheared zone. Jamming pressure p j appears to increase towards the centre of the plug, as do the components of the particle stress tensor τ p,rz and (τ p ). Note however that stress values in the plug are likely to be strongly influenced by the considered regularisation for the constitutive law. One cannot argue that stress profiles in this region converge to those that would be obtained with a more sophisticated constitutive law avoiding the need for a regularisation. Finally, also observe that the jamming pressure displays a small discontinuity at the transition between the plug and the sheared region at intermediate times, but that this discontinuity tends to vanish in the steady-state regime.

Steady-state solutions have been computed for five different values of average reduced volume fraction ψ 0 ranging between 0.55 and 0.78. The corresponding values of φ 0 and φ m are indicated in Figure 6. The final computation time T used for each case was varied to ensure that a fully-developed flow regime is effectively reached. For that purpose, a steady-state criterion based on the exponential decay of w r (see Fig. 3.left) is used. More precisely, the time loop is stopped when w r (t) is reduced by a factor of 100 compared to its value at the end of the first iteration. Observe that the final dimensionless time T decreases as the average reduced volume fraction ψ 0 increases (see Fig. 6).

Figure 7 illustrates the influence of average volume fraction φ 0 on the steadystate solution, for a fixed value of φ m . As expected, the width of the plug region strongly increases with φ 0 , while the average bulk velocity decreases. The jamming pressure in the plug, as well as the different components of the particle stress tensor, are also observed to increase when φ 0 increases. For completeness, let us also mention that the Lagrange multiplier f z is negative and decreases with φ 0 : steady-state values are f z (T ) = -10.80 for φ 0 = 0.32 and f z (T ) = -43.84 for φ 0 = 0.50. Recall that f z interprets as the longitudinal gradient of the bulk pressure.

Comparison with experiments

This paragraph presents direct comparisons between the experimental measurements of Oh et al. [START_REF] Oh | Pressure-driven suspension flow near jamming[END_REF], performed on fully developed flows, and steady-state solutions of our migration model. Oh et al. [START_REF] Oh | Pressure-driven suspension flow near jamming[END_REF] report data for several values of average particle volume fraction φ 0 . We retained here the results corresponding to a semi-concentrated and to a concentrated case, with nominal values φ 0 = 0.35 and φ 0 = 0.52 respectively. We observed that these nominal values of volume fraction slightly differ from the values obtained by direct integration of the measured radial particle fraction profiles. This computation to effective values φ 0 = 0.32 and φ 0 = 0.5 for the semi-concentrated and concentrated cases, respectively. To avoid systematic discrepancies between experimental and numerical profiles, the model solutions were computed for these effective values of φ 0 .

Choosing the value of the maximal volume fraction φ m in the model requires care. As explained by Lecampion and Garagash [START_REF] Lecampion | Confined flow of suspensions modelled by a frictional rheology[END_REF], volume fraction φ is actually characterised by two noticeable limits in highly-concentrated mixtures. The first limit is the random close packing fraction φ rcp = 0.64, which cannot be exceeded. The second limit is the critical volume fraction φ c ≈ 0.585, above which the mixture is jammed, i.e. behaves as a solid. In regions where φ > φ c , volume fraction can still continue to increase by compaction, as the particle phase behaves as a compressible solid. As shown by Oh et al. [START_REF] Oh | Pressure-driven suspension flow near jamming[END_REF], the actual steadystate volume fraction reached in the jammed regions depends upon φ 0 . We may hypothesise that this steady-state volume fraction results from a complex balance between particle jamming pressure p j and friction. Since the constitutive law of Morris and Boulay [START_REF] Morris | Curvilinear flows of noncolloidal suspensions: the role of normal stresses[END_REF] considered in the present model does not account for friction nor for solid compressibility in the jammed regions, the apparent maximal volume fraction φ m to be considered remains a priori unknown in the range [φ c , φ rcp ]. Accordingly, for each average volume fraction φ 0 , solutions corresponding to several values of maximal volume fraction φ m were computed and compared to experimental data.

Fig. 8 presents comparisons between the steady-state solution at t = T and experimental measurements for the semi-concentrated case (φ 0 = 0.32). It is observed that the radial profile of mixture velocity u z across the tube section is very well reproduced by the model, and this for all choices of φ m in the range 0.585-0.64 (see Fig. 8.top-left). In particular, the velocity value in the plug appears to be well captured. Concerning particle volume fraction φ (see Fig. 8.top-right), discrepancies can be observed concerning the width of the plug and the sharpness of the transition with the sheared layer. Experimental data seem to indicate thinner plugs, and a smoother transition between the two zones. Note however that the evolution of φ in the vicinity of the wall of the tube is well reproduced by the model. It is also observed that, while the choice of φ m in the model has obviously a direct impact on the value of particle volume fraction reached in the plug, this parameter has only little influence on the profiles of φ in the sheared region. Experimental measurements appear to be best captured with φ m = 0.60 in this case.

Similar observations can be made for the concentrated case (φ 0 = 0.50). Here also, the radial profile of u z is well reproduced in the whole tube section for all choices of φ m (see Fig. 8.bottom-left). In this case, the value of volume fraction in the central plug is best captured with the upper bound φ m = 0.64 (see Fig. 8.bottom-right), suggesting that the maximum volume fraction reached in the plug appears to increase with the average concentration of the suspension. In contrast, it is again observed that the value of φ near the wall depends only slightly on the choice of φ m .

The discrepancies between numerical solutions and experimental measurements observed on volume fraction profiles are probably attributable to the limitation of the model pointed out above, namely that the constitutive law for the particle phase does not account for friction and solid compressibility in jammed regions. Experimental results clearly indicate that volume fraction continues to increase beyond φ c in the plug, and saturates only when reaching the limit φ = φ rcp . In contrast, in the model, a single maximum value φ m is considered. Extension of the constitutive law to account for such granular processes, as proposed by Lecampion and Garagash [START_REF] Lecampion | Confined flow of suspensions modelled by a frictional rheology[END_REF], thus represents a promising prospect to further improve the agreement with experimental data.

Conclusion and perspectives

This paper presents a new migration model for mono-disperse suspensions of neutrally buoyant particles. Unlike suspension balance model (SBM), which relies on a single velocity formulation [START_REF] Miller | Normal stress-driven migration and axial development in pressure-driven flow of concentrated suspensions[END_REF], our model involves two velocities and two pressures. The two-velocity formulation, coupled to the introduction of a diffusive term on the migration velocity w, allows us to properly impose non-penetration boundary conditions for the particles across walls, and thus to satisfy rigorous mass conservation for the particle phase. Physically, the new diffusive term can be interpreted as a non-local contribution to the particle stresses in unjammed regions. In addition, the unilateral constraint φ φ m on the particle volume fraction is strictly imposed through the introduction of a particle jamming pressure p j . This Lagrange multiplier, which takes nonzero values only in jammed regions, interprets as the contribution to the particle pressure of the collective interactions (contact chains) that develop between particles in these regions.

Through an asymptotic analysis, a reduced 1D migration model is derived in the case of an axisymmetric Poiseuille flow, and a fully implicit algorithm is proposed for computing numerical solutions of this reduced model. The originality lies in the handling of the unilateral constraint, through an augmented Lagrangian method embedded in a fixed point iteration at each time step. This algorithm is coupled to a finite element spatial discretisation, and the convergence properties of the scheme are carefully demonstrated. In particular, the inner loops involved in the augmented Lagrangian and the fixed point are shown to converge in a small number of iterations. The error on particle mass is also shown to converge linearly with mesh refinement.

The physical characteristics of the solutions are described as a function of average volume fraction, and quantitative comparisons with experimental measurements in fully-developed flow regime [START_REF] Oh | Pressure-driven suspension flow near jamming[END_REF] are presented for both a semiconcentrated and a concentrated case. Mixture velocity profiles, characterised by the formation of a central plug, appear to be very well reproduced. The model is thus able to accurately capture the effects of jamming and the effective viscoplastic behavior [START_REF] Boyer | Unifying suspension and granular rheology[END_REF] in highly-concentrated mixtures. Discrepancies between model predictions and experimental data remain nevertheless visible on the volume fraction profiles, particularly at the transition between the plug and the outer sheared region. In experiments, the transition between the two zones appears smoother than in the numerical solutions.

To overcome the current limitations of the model, future improvements shall concentrate on including more sophisticated constitutive relations accounting for specific granular processes in concentrated regions. In particular, inelastic compressibility beyond the jamming limit seems to be a necessary ingredient to better capture the evolution of particle volume fraction observed in plug regions [START_REF] Lecampion | Confined flow of suspensions modelled by a frictional rheology[END_REF]. Inclusion of a true yield stress in jammed zones through, e.g., frictional effects Saramito [see, e.g., 55], would also be necessary to avoid the use of a regularised apparent viscosity and predict more realistic particle stress values in these regions. Finally, the development of anisotropic micro-structures in both jammed and sheared regions could also be considered through specific tensorial constitutive relations [START_REF] Ozenda | Tensorial rheological model for concentrated non-colloidal suspensions: normal-stress differences[END_REF]. Once these more realistic constitutive models will be implemented and tested, future work shall also consider the extension of the numerical algorithm in order to efficiently address more complex geometries such as 3D channel flows, flows around obstacles [START_REF] Haddadi | Suspension flow past a cylinder: particle interactions with recirculating wakes[END_REF][START_REF] Oucheggou | Particle-laden flow around an obstacle in a square pipe: experiments and modeling[END_REF], or re-suspension experiments [START_REF] Blanc | Viscous resuspension of non-brownian particles: determination of the concentration profiles and particle normal stresses[END_REF]. Finally, another promising line of improvement of the numerical resolution would consist in implementing mesh adaptation along the interface between the plug and the sheared region, as done for Bingham fluids by, e.g., Roquet and Saramito [START_REF] Roquet | An adaptive finite element method for Bingham fluid flows around a cylinder[END_REF]. and is used in post-treatment to obtain a first order approximation in ξ of mixture pressure p. Similarly, wz appears only in equation (B.3b); consequently, this equation can be integrated explicitly after the resolution of the five remaining equations. The system is closed by suitable initial and boundary conditions, and is summarised in [START_REF] Guazzelli | Rheology of dense granular suspensions[END_REF]. Finally note that, for convenience, tilde notations are dropped in the main text (section 2.2), and the two-velocities ũ and w are normalised by the same characteristic velocity U .

Appendix C. Uzawa algorithm

Problem (O), as defined in subsection 3.2, is solved by minimising the cost function J defined by equation [START_REF] Ozenda | A new rate-independent tensorial model for suspensions of non-colloidal rigid particles in newtonian fluids[END_REF]. Let us introduce the dual function J * defined for all q ∈ L 2 by: J * (q) = -min

(v,ξ)∈H 1 0 ×L 2 L(v, ξ ; q)
The problem writes equivalently as a minimisation problem for this dual function:

p j = arg min q∈L 2 J * (q)
The numerical procedure used to solve this problem, based on Uzawa method, is described in subsection B.1, with two technical lemmas gathered in subsection B.2.

Appendix C.1. Numerical algorithm

The Uzawa algorithm for the augmented Lagrangian method expresses as a constant-step descent algorithm for the dual function J * :

• m = 0: let p (0) j be given.

• m 1: let p .

Here, m ∈ N denotes the descent loop index, which is implemented as an inner loop inside the fixed point loop with index k introduced in subsection 3.1. A numerical parameter m max is introduced, such that the stopping criterion for this inner loop is defined as m > m max . Note that the constant-descent step has been chosen equal to the augmentation parameter µ. Note also that J * is differentiable. Then, expanding its gradient, we obtain an equivalent formulation of the descent algorithm:

• m = 0: let p (0) j be given. [START_REF] Miller | Normal stress-driven migration and axial development in pressure-driven flow of concentrated suspensions[END_REF]. The present choice α = 3 corresponds to an intermediate value. We observed that the smaller the value of α, the faster the system reaches the steady state regime.

Figure 1 :

 1 Figure 1: Circular tube geometry for the Poiseuille flow.

Fig. 2 .Figure 2 :

 22 Figure 2: Convergence of the two inner loops at t = 3T /8 (h = 1/800). (left) Convergence of the augmented Lagrangian loop: normalised relative error on the auxiliary variable δ (m) versus iteration number m for different values of augmentation parameter µ. (right) Convergence of the fixed point loop: relative error for the volume fraction ψ (n,k) versus total iteration number k × mmax for different values of mmax.

dr 1 hFigure 3 :

 13 Figure 3: Convergence versus simultaneous space and time refinement (see Table4), with T = 4000, φ 0 = 0.32, φm = 0.585. (left) L 2 norm of the radial component w r,h of migration velocity versus dimensionless time t. (right) Maximum particle mass error versus mesh size h.

  h r c,h (T ) |r c,h -r c,h4 | 1

Figure 4 :Figure 5 :

 45 Figure4: Convergence versus mesh refinement: radial profiles at t = T = 4000 of longitudinal bulk velocity uz, reduced particle volume fraction ψ, radial and longitudinal migration velocity wr and wz, bulk shear stress σm,rz, and jamming pressure p j (φ 0 = 0.32, φm = 0.585).

Figure 7 :

 7 Figure7: Comparison of the steady-state solutions for two different values of initial volume fraction φ 0 = 0.32 and 0.50, while φm = 0.585 (h = 1/1600). Radial profiles of longitudinal bulk velocity uz, reduced particle volume fraction ψ, radial and longitudinal migration velocity wr and wz, jamming pressure p j , and particle stress shear and spherical τ p,rz and tr(τ p).

Figure 8 :

 8 Figure 8: Comparisons of the present migration model with experimental results of Oh et al.[START_REF] Oh | Pressure-driven suspension flow near jamming[END_REF] for the fully developed flow regime: steady-state radial profiles of bulk longitudinal velocity uz and particle volume fraction φ. (top) Semi-concentrated case φ 0 = 0.32. (bottom) Concentrated case φ 0 = 0.50. Error bars on the experimental data correspond to uncertainties estimated by Oh et al.[START_REF] Oh | Pressure-driven suspension flow near jamming[END_REF], based on measurements for the interstitial fluid and might thus be underestimated for concentrated mixtures.

2 η

 2 Figure D.9: Steady-state normal stress ratio N 2 as a function of reduced volume fraction ψ: comparison between the predictions the rheological model (with values indicated in table 2 and φm = 0.60) and experimental data of Dbouk et al. [61] and Couturier et al. [62].

Table 1 :

 1 are shown. Section 5 presents final discussions and conclusions. Notations used in the paper.

	notation description	notation	description
	u	bulk velocity	ρ	fluid and particle density
	w	migration velocity	η 0	fluid viscosity
	q	pipe flow rate	τ p	particle stress tensor
	p	mixture pressure	Q = diag(1, λ 2 , λ 3 ) normal stress tensor
	p j	jamming pressure	s(ψ)	hindrance function
	φ	volume fraction	α	exponent in s(ψ)
	φ m	maximal volume fraction	K n , K s	reduced viscosities
	ψ = φ/φ m reduced volume fraction	ε = r p /R	dimensionless particle radius
	Ω	flow domain	Re = 2ρq/(πR)	Reynolds number
	R L T r p	tube radius tube length final time particle radius	h ∆t µ	mesh size time step augmentation parameter

Table 2 :

 2 Values of model material parameters.

	4.1. Experimental setup
	symbol	value	unit
	ρ	1056	kg.m -3	symbol value unit
	η 0	3.6	Pa.s	K s	0.6
	r p	7 × 10 -5	m	K n	1
	R	3.15 × 10 -3	m	λ 2	0.9
	q	3.14 × 10 -8 m 3 .s -1	λ 3	0.5
	φ 0	0.32 ; 0.5		α	3
	φ m	0.585 -0.64

Table 3 :

 3 Numerical parameters of the algorithm.

Table 4 :

 4 Parameters of the numerical discretisation for dimensionless final time T = 4000.

Table 5 :

 5 Convergence of critical plug radius r c,h versus mesh refinement h.

  Figure 6: Final dimensionless computation time T , for which |wr| is reduced by a factor of 100 with respect to its initial value, for the five simulation runs discussed in the text.The corresponding values of average reduced volume fraction ψ 0 = φ 0 /φm, average volume fraction φ 0 , and maximal volume fraction φm, are indicated in the table. Other are given in Table2. Finest mesh size h = 1/1600 is used.

	1	w r (t, .) w r (0, .)		run 1 2 3 4
					5
	10 -1
	0 T 5 T 4			T 3 T 2	T 1 8000
				t
	run	1	2	3	4	5
	φ 0		0.32		0.50
	φ m	0.585 0.60 0.64	0.585	0.64
	ψ 0	0.55	0.53 0.50	0.85	0.78
	T	6973	5603 4940 1481	640

Appendix A. Link with mixture theory

The mixture theory developed by Jackson [START_REF] Jackson | The Dynamics Of Fluidized Particles[END_REF] is based on an asymptotic analysis for small particle radius r p . Ozenda [60, chap. 2] recently revisited this asymptotic analysis and obtained, at first order in r p , the following two-velocity system of conservation equations (see also Nott et al. [23]):

where σ h and σ c denote contributions to the Cauchy stress of the mixture due to hydrodynamic and contact interactions between the particles, respectively, p f is the fluid phase pressure, and f h is the hydrodynamic force exerted on the particle phase.

Present model Jackson's mixture theory

Table A.6: Correspondence between the notations used in system (A.1) and those used by Jackson [START_REF] Jackson | The Dynamics Of Fluidized Particles[END_REF]. Three different averages are introduced, namely the volume average over the suspension . , the volume average over the fluid phase . f , and the discrete particle phase average . p, with ν the number density of particles. See also Ozenda [60,chap. 2].

The correspondence between the notations used in the above system and the different phase averages introduced by Jackson [START_REF] Jackson | The Dynamics Of Fluidized Particles[END_REF] are explained in Table A.6 (see also Ozenda [60,p. 45] for more details). Identifying (2a) with (A.1a), and (2b) with (A.1b), we obtain the following expressions for the Cauchy stress tensor of the mixture and the forces in the particle phase, respectively:

The following closure relation can be considered:

Then, Equations (A.2a) and (A.2b) successively lead to

The first relation corresponds to a classical closure in mixture theory [46, chap. 2], through which the hydrodynamic and contact contributions σ h , σ c are lumped into a particle stress, which writes here -p j I + τ p . The second relation writes equivalently:

This expression identifies term by term with Equation (62) of Nott et al. [START_REF] Nott | The suspension balance model revisited[END_REF],

where the first term on the right-hand side represents the drag force, and the second term is a particle phase hydrodynamic stress. Note that the new corrective term introduced in our model identifies as the difference between the particle phase hydrodynamic stress and the hydrodynamic stress σ h involved in the mixture momentum conservation (A.1a), as discussed by Nott et al. [START_REF] Nott | The suspension balance model revisited[END_REF].

Since this corrective term is of second-order with respect to r p , it does not change the overall accuracy of the mixture model.

Appendix B. Poiseuille flow in a long tube

We consider the circular tube geometry represented on Fig. 1, with (r, θ, z) the associated cylindrical coordinate system. This appendix shows how problem (2a)-(2i) reduces asymptotically to problem (3a)-(3j) when tube length L becomes large with respect to its radius R. The flow is assumed to be axisymmetric, i.e. independent upon θ. No-slip conditions are assumed at tube wall i.e. u Γ = 0 in (2g). Moreover, we assume that the initial conditions u 0 and w 0 satisfies u 0,θ = w 0,θ = 0.

Let U be a characteristic velocity of the mixture and W be a characteristic migration velocity; R/U is then a characteristic time and η 0 U/R is a characteristic stress. Dimensionless variables are denoted by tildes, e.g. t = (U/R)t, r = r/R, z = z/L and ũ = (ũ r , ũθ , ũz ) = u/U. Finally, let ξ = R/L denote the tube aspect ratio, ε = r p /R denote the dimensionless particle radius, and ζ = W/U .

Appendix B.1. Mixture subsystem

Let us first consider the mixture subsystem associated to the unknowns ũ and p.

Mixture momentum and mass conservations (2a) and (2c) become:

where the Reynolds number is defined by Re = ρU R/η 0 . For ξ → 0, relation (B.1d) reduces to ∂ r (rũ r ) = 0. Hence, the boundary condition (2g) yields ũr ( t, r = 0) = ũr ( t, r = 1) = 0 at any time t ∈ ]0, T [. Thus, ũr = 0 at any time.

Moreover, since we assume ũθ = 0 at t = 0, from (B.1b), this identity remains true at any time. Thus ũ = (0, 0, ũz ).

Observe in (B.1a) that the term ∂ r τp,rr , responsible for the migration, appears at the same order in ξ as the pressure derivative ∂ r p. Conversely, in (B.1c), the term ξ ∂ z p should be at zeroth order in ξ; otherwise the suspension would not move in a long tube when ξ → 0. Finally, the pressure is assumed to admit the following expansion in ξ:

where fz depends only upon t, and p0 depends upon t and r but is independent of z, as inferred from (B.1a). Then, for ξ → 0, momentum conservation (B.1a) and (B.1c) reduce to

Two types of controls can be considered for the flow of the mixture: either the pressure drop fz or the flow rate can be imposed. Here, we choose to impose the flow rate, denoted by q. The characteristic mixture velocity is then defined as U = 2q/(πR 2 ), which is equal to twice the average velocity, such that the dimensionless flow rate expresses as

and fz in (B.2b) interprets as a Lagrange multiplier for the imposition of the flow rate constraint (B.2c). Note that when ψ = 0, from (2f), the particle stress τ p vanishes, the fluid is Newtonian and U coincides with the maximal value of the mixture velocity.

Appendix B.2. Congested subsystem

Let us now turn to the congested subsystem associated to the unknowns w, pj and ψ. Since wθ = 0 at t = 0, we assume that w = ( wr , 0, wz ) at any time.

Note that wr = 0 in general, since the particles are expected to migrate in the tube. The momentum and mass conservations of the particle phase (2b) and (2d) become:

By passing to the limit ξ → 0 and neglecting second-order terms in ζ, which are associated to inertia effects, the previous system becomes:

Finally, let us turn to the expression of particle stress tensor τ p given by (2f). Since ũ( t, r) = (0, 0, ũz ( t, r)), the norm of the strain rate writes |2 D(ũ)| = |∂ r ũz | and the components of τ p express as explicit relations involving ũz and ψ, as shown by ( 4).

The final system involves seven equations, namely (B.2a)-(B.2c), (B.3a)-(B.3c), and (1), and seven unknowns: ũz , p0 , fz , wr , wz , pj , and ψ. All these unknowns depend both upon time t and radial position r, except for fz that only depends upon time. Observe that Equation (B.2a) leads to an explicit computation of p0 ,

• m 1: let p (m-1) j be known. Then compute successively:

To simplify the simultaneous minimisation versus (v, ξ) of the Lagrangian, the algorithm is modified by decoupling the first step as:

• m = 0: let p (0) j and δ (0) be given.

• m 1: let p (m-1) j and δ (m-1) be known. Then compute successively:

r From ( 12), the Lagrangian L is quadratic and differentiable versus v. Thus, the first step of the above algorithm reduces to a linear sub-problem, namely,

∈ H 1 0 such that, for all v ∈ H 1 0 , we have ∂L ∂v w (m) r , δ (m-1) , p

The Lagrangian L is both nonlinear and non-differentiable versus ξ, but involves a sub-differential with respect to ξ. Hence, the optimal value δ (m) verifies 0 ∈ ∂ ξ L(w

), where ∂ ξ L(w

) is the sub-differential of L with respect to the variable ξ. The second step of the descent algorithm is then solved locally. The sub-gradient of the indicator of [g(r), ∞[, denoted ∇I [g(r),∞[ , verifies, for all ξ ∈ R:

{0} otherwise For all r ∈ [0, 1], the second step then expresses: find δ (m) (r) ∈ R such that: + µδ (m-1) Bv r dr (C.1a)

ii) compute explicitly:

This algorithm has been shown to converge for any µ > 0 [53, theorem 5.1]).

Recall that the solution is independent upon µ. At convergence, (C.1c) leads to δ = Bw r and then (C.1a) is exactly a weak formulation of (11a). Since w r ∈ H 1 0 in (C.1a), we also obtain (11c). Replacing δ = Bw r in (C.1b), we obtain Bw r = max g, Bw r -p j µ which leads to (11b), as shown in Appendix C.2 (lemma 2). Thus, the previous algorithm effectively provides a solution of (11a)-(11c).

Appendix C.2. Technical lemmas

Lemma 1. With the notations defined in section 3, we state the following equivalence relation:

The maximum condition is split between two cases that are equivalent to the right hand side disjunctive relation in lemma 1: Bw r = max g, Bw r -p j µ =⇒ Bw r g Bw r Bw r -p j µ

• The condition Bw r g is evidently satisfied.

• Recalling µ > 0, the positivity of p j is deduced :

Bw r Bw r -p j µ =⇒ p j 0

• The condition p j (Bw r -g) = 0 is finally stated:

Hence the unilateral condition g Bw r ⊥ p j 0 is stated.

Appendix D. Identification of the rheological parameters

Steady state flow profiles computed with our model show a very low sensitivity to the values of rheological parameters K n , λ 2 , λ 3 and α. Hence, only the value of K s was adjusted to fit with the migration data of Oh et al. [START_REF] Oh | Pressure-driven suspension flow near jamming[END_REF] (see Fig. 8). Values of λ 2 and λ 3 define the relative magnitude of the normal stress