Christine Frances
email: christine.frances@ensiacet.fr.

Alain Lin

Comminution Process Modeling Based on the Monovariate and Bivariate Direct Quadrature Method of Moments

Keywords: particle technology, mathematical modeling, grinding, population balance, direct quadrature method of moments

Population balance equations (PBE) applied to comminution processes are commonly based on selection and breakage functions, which allow the description of changes of particle-size distributions vs. time. However, other properties, such as the particle strength, also influence the grinding kinetics. A bivariate PBE was developed and resolved by the direct quadrature method of moments. This equation includes both the particle size and strength, the latter of which is defined as the minimum energy required for breakage. The monovariate case was first validated by comparing the predicted moments with those calculated from the size distributions given by an analytical solution of the PBE derived for specific selection and breakage functions. The bivariate model was then compared with a discretized model to evaluate its validity. Finally, the benefit of the bivariate model was proven by analyzing the sensitivity of some parameters and compar-ing the results of the monovariate and bivariate cases.

Introduction

Population balance equations (PBE) are commonly used to model particulate processes. In the framework of comminution processes, PBE tracks the change in the density distribution of the weight of particles distributed over defined size classes. Batch grinding mass balance equations usually incorporate the selection function, which expresses the breakage probability of particles as a function of their size, and the breakage function, which describes the daughter particle-size distribution. This approach permits the optimization of the size reduction process and has intensively been used in the past to model ball mills 1,2 and more recently to help determine breakage and selection functions parameters. 3,4 The kinetics parameters are then back-calculated from the PBE. A similar approach was also used to model the simultaneous breakage and reagglomeration of fines in a grinding process 5 and to identify the grinding mechanisms in a submicron milling process. 4,[START_REF] Hennart | Characterization and modeling of a sub-micron milling process limited by agglomeration phenomena[END_REF] The modeling of continuous grinding processes is also based on PBE, which combines the fragmentation kinetics with the particle residence time distribution. The batch mass balance equation is first solved, which yields the temporal change of the mass fractions in each size class. The grinding kinetics is then combined with a global transport model to obtain the size distribution of the ground product during a continuous process. [START_REF] Frances | On modeling of submicronic wet milling processes in bead mills[END_REF] Many analytical or approximated solutions of the batch grinding equation have been proposed in the literature 2,8 depending on the specific forms of the selection and breakage functions. Irrespective of the method used to solve the mass balance equation, the mass balance can describe the fragment-size distribution as a function of the grinding time.

However, size is not the only important solid characteristic in a grinding process, as the particle strength (or grindability) is another key parameter that directly affects the process kinetics. Grindability can be the cause of what Austin et al. [START_REF] Austin | Process Engineering of Size Reduction: Ball Milling[END_REF] called an abnormal breakage behavior, which is characterized by an increase in the selection function vs. time (when a fatigue phenomenon occurs) or a decreasing tendency (when small fragments are more resistant than larger ones).

The usefulness of the actual population balance models for mill design, optimization, or control is thus reduced if the particle size is the unique property. An innovative approach was recently proposed by Crespo, [START_REF] Crespo | Application of particle fracture energy distributions to ball milling kinetics[END_REF] who developed a PBE for a ball-milling process in which the weight distribution was jointly characterized by the particle size and strength. In the context of grinding processes, the particle strength can be practically defined by expressing the minimum energy required to break the particle. Such a property can be evaluated by performing crushing tests on a single particle, such as with the ultrafast load cell (UFCL). [START_REF] King | Measurement of fracture energy during single-impact fracture[END_REF][START_REF] Tavares | Single-particle fracture under impact loading[END_REF] The method used to characterize the fracture energy leads to a statistical treatment of the results to consider the variability of the particle properties. The particle strength is thus defined by an energy spectrum that depends on particle size. The PBE applied to comminution processes can finally be written for the breakage energy-size distribution.

An important quantity for grinding processes is the energy consumed by the mill, a part of which is transmitted to the particles. For example, Kwade [START_REF] Kwade | Determination of the most important grinding mechanism in stirred media mills by calculating stress intensity and stress number[END_REF] developed theoretical considerations of interactions between the fracture energy spectra induced by bead motion to characterize the material properties and impact energy for stirred media mills. The impact energy spectra depend on the process conditions and can be determined by computing the collision energy via the discrete element method, [START_REF] Tuzcu | Modeling breakage rates in mills with inpact energy spectra and ultrafast load cell data[END_REF][START_REF] Tuzcu | Modeling ball mills and semiautogeneous mills with impact energy spectra-based population balance concept[END_REF] simulating the hydrodynamics inside the grinding chamber using computation fluid dynamics (CFD) codes based on Reynolds-averaged Navier-Stokes [START_REF] Lane | CFD modeling of a stirred bead mill for fine grinding[END_REF] or direct numerical simulation. [START_REF] Gers | Numerical simulation of grinding in a stirred media mill: hydrodynamics and collision characteristics[END_REF] Thus, a first main challenge of comminution processes modeling is to simultaneously consider two variables, such the particle size and fracture energy in the PBE, which is addressed in this article. A further challenge would be to combine the PBE with a simulation tool to predict the solid property by considering the impact energy spectra.

Crespo [START_REF] Crespo | Application of particle fracture energy distributions to ball milling kinetics[END_REF] solved a discretized form of the homogeneous bivariate PBE. However, this method only considers a few size and fracture energy intervals and cannot easily be coupled to CFD code if flow or impact energy heterogeneities must be considered.

In this work, the batch grinding equation is solved by the direct quadrature method of moments (DQMOM); one of the advantages of this method is that it can solve the PBE by incorporating several properties of the solid phase. In the QMOM or DQMOM approaches, the mixed moments of the density distribution must be expressed by the quadrature formula by incorporating the weights and abscissas defined for each property. Over the last decade, QMOM and DQMOM modeling methods have been successfully used for a variety of particulate processes: aggregation and breakage processes in homogeneous or one-dimensional flows, [START_REF] Marchisio | Quadrature method of moments for aggregation-breakage processes[END_REF][START_REF] Marchisio | Implementation of the quadrature method of moments in CFD codes for aggregation-breakage problems[END_REF][START_REF] Marchisio | Solution of population balance equations using the direct quadrature method of moments[END_REF] growth and crystallization processes, [START_REF] Marchisio | Quadrature method of moments for population-balance equations[END_REF] coagulation and sintering, [START_REF] Fox | Bivariate direct quadrature method for coagulation and sintering of particle populations[END_REF] and nanoparticle processes, [START_REF] Zucca | Validation of bivariate DQMOM for nanoparticle processes simulation[END_REF] among others. Examples of the use of these methods have also been reported for other dispersed systems, such as liquid sprays undergoing droplet coalescence and evaporation, [START_REF] Fox | Numerical simulation of spray coalescence in an Eulerian framework: direct quadrature method of moments and multifluid method[END_REF] bubble coalescence and breakage in vertical gas-liquid flows, [START_REF] Deju | Study of isothermal vertical bubbly flow using direct quadrature method of moments[END_REF][START_REF] Deju | Capturing coalescence and break-up processes in vertical gas-liquid flows: assessment of population balance methods[END_REF] or in stirred tank reactors. [START_REF] Petitti | Bubble size distribution modeling in stirred gas-liquid reactors with QMOM augmented by a new correction algorithm[END_REF][START_REF] Buffo | Multidimensional population balance model for the simulation of turbulent gas-liquid systems in stirred tank reactors[END_REF][START_REF] Buffo | Multivariate quadrature-based moments methods for turbulent polydisperse gasliquid systems[END_REF] The breakage phenomena in dispersed systems have thus been intensively treated in the literature by QMOM and DQMOM models. However, the breakage kernels used in these previous studies were usually different from those commonly used in comminution processes and did not consider two properties for the solid phase, such as the particle size and strength; the second being a complex property conditioned by the first. Moreover, another advantage of the DQMOM method is that it solves the spatial and temporal transport equations directly on weights and abscissas, which allows the implementation of the population balance model in CFD codes. [START_REF] Marchisio | Investigation of soot formation in turbulent flames with a pseudobivariate population balance model[END_REF][START_REF] Pasquale | Model validation for precipitation in solvent-displacement processes[END_REF][START_REF] Silva | Implementation and analysis of numerical solution of the population balance equation in CFD packages[END_REF][START_REF] Silva | Development and implementation of a polydispersed multiphase flow model in OpenFOAM[END_REF] Thus, DQMOM model provides a promising approach for prospective studies on heterogeneous comminution processes.

In this study, the DQMOM method was used to solve the PBE applied to homogeneous batch grinding processes. In the first step, a single variable, the particle size, is considered, and the results of the model were compared to an analytical solution of the mass balance equation. In the second step, a bivariate model is developed to track the change in the density distribution, which was simultaneously expressed in terms of the particle size and strength. The results of the model are compared with a discretized solution of the PBE by taking a finite number of size and energy intervals. Finally, the benefit of the bivariate approach is discussed based on the sensitivity of the results to the kinetic parameters of the grinding process.

Monovariate DQMOM Applied to Comminution Processes General equations: homogeneous breakage

For sake of simplicity, only the homogeneous equation of population balance is addressed. The population balance can be expressed as follows by introducing the particle size x and the mass particle-size distribution f (x, t) instead of the number-density function

@f x; t ð Þ @t 5S x x; t ð Þ (1)
In the case of pure breakage, the right-hand side contains only two components, a sink term expressing the breakage of particles of size x in smaller ones and a source term due to the breakage of particles larger than x

@f x; t ð Þ @t 52a x ð Þf x; t ð Þ1 ð 1 x a y ð Þb x; y ð Þf y; t ð Þdy (2)
The DQMOM method is based on the Gaussian quadrature. [START_REF] Marchisio | Solution of population balance equations using the direct quadrature method of moments[END_REF][START_REF] Fox | Bivariate direct quadrature method for coagulation and sintering of particle populations[END_REF] The continuous size distribution can then be substituted with the discrete decomposition

f ðx; tÞ5 X N a51 w a t ð Þd xðtÞ2x a ðtÞ ½ ( 3
)
where N is the number of nodes a, x a is the property of the node (its size), w a is its weight, and d is the Dirac function.

If this discrete expression of the mass particle-size distribution f (x, t) is injected into transport Eq. 2, one obtains the following

X N a51 d x2x a ð Þ @w a @t 2 X N a51 d 0 x2x a ð Þw a @x a @t 5S x x; t ð Þ (4)
Following Marchisio and Fox,[START_REF] Marchisio | Solution of population balance equations using the direct quadrature method of moments[END_REF] one can substitute the weighted abscissa, c a , for the abscissa, x a , which yields the following c a 5w a x a (5)

Equation 4 then takes the following form

X N a51 d x2x a ð Þ1d 0 x2x a ð Þx a ½ @w a @t 2 X N a51 d 0 x2x a ð Þ @c a @t 5S x x; t ð Þ (6)
The derivatives with time of both the weights and the weighted abscissa are independent of x. Thus, a set of ordinary differential equations is usually defined to follow the transient changes of these weights and weighted abscissas @w a @t 5E wa ðtÞ @c a @t 5E ca ðtÞ

The integer moments of the distribution were defined by Hulburt and Katz [START_REF] Hulburt | Some problems in particle technology[END_REF] as follows

m k ðtÞ5 ð 11 21 x k f ðx; tÞdx5 X N a51 w a x k a ( 8
)
Equation 6 can be transformed into a transport equation of moments by multiplying it by x k and integrating over x. Different terms need thus to be derived, such as

ð 11 21 x k dðx2x a Þdx5x k a ( 9
)
ð 11 21 x k d 0 ðx2x a Þdx52kx k21 a (10)
Thus, the final equation can be written as follows

ð12kÞ X N a51 x k a E wa 1k X N a51 x k21 a E ca 5 ð 11 21 x k S x x; t ð Þdx (11)
This set of Eq. 11 corresponds to an algebraic system, with 2N unknowns: E w1 , . . ., E wN , E c1 , . . ., E cN . Following Marchisio and Fox, [START_REF] Marchisio | Solution of population balance equations using the direct quadrature method of moments[END_REF] one can write this system in matrix form A E 5 D, where A is a square matrix (2N, 2N) and E and D are column vectors of size 2N. Thus, 2N moments are needed to close the system, so k must vary between 0 and 2N 2 1. One can consider two particular values of k, namely k 5 0 and k 5 1.

When k 5 0, the previous equation simplifies to the following form

X N a51 E wa 5 ð 11 21 S x x; t ð Þdx50 (12)
On the right hand side, the summation of the source and sink terms over the whole distribution of sizes must be null, irrespective of the kernel. Keeping in mind the definition of E wa , one can rewrite this expression as follows

@ @t w 1 1w 2 1w 3 ð Þ 5 @m 0 @t 50 ( 13
)
This result is obvious. When k 5 1, the previous equation reduces to

X N a51 E ca 5 ð 11 21 x S x x; t ð Þdx5S 1 (14)
Here, one can rewrite this expression as follows by again keeping in mind the definition of E ca

@ @t w 1 x 1 1w 2 x 2 1w 3 x 3 ð Þ 5 @m 1 @t 5S 1 ( 15
)
The transient change in the first moment could be predicted depending on the closure of the sink and source terms. The kernels associated with the present comminution process will be defined in the next paragraph.

Monovariate DQMOM validation and discussion

In the specific case where the selection function a(x) is a power function (Eq. 16) and the breakage distribution function B(x, y) is self-similar as in Eq. 17 [START_REF] Frances | On modeling of submicronic wet milling processes in bead mills[END_REF][START_REF] Austin | Solutions of the batch grinding equation leading to Rosin-Rammler distributions[END_REF][START_REF] Austin | Some results on the description of size reduction as a rate process in various mills[END_REF][START_REF] Kapur | Self-preserving size spectra of comminuted particles[END_REF]

a x ð Þ5A 0 x x max a 0 (16) B x; y ð Þ5 x y a0 or b x; y ð Þ5 dB x; y ð Þ dx 5 a 0 y x y a021 (17)
Then, the batch grinding equation has been shown to have an analytical solution

R x; t ð Þ5R x; 0 ð Þexp 2A 0 x x max a 0 t (18)
where R(x, t) is the retained cumulative-size distribution defined as follows

R x; t ð Þ5 ð 1 x f y; t ð Þdy (19)
The DQMOM model was developed using Matlab and validated considering a normal law for the initial particlesize distribution with a mean value (m) equal to 3 and a standard deviation (r) of 0.5. If the order N of the quadrature is chosen to be 3, the weights and abscissas can be calculated as function of the mean and standard deviation of the normal law by the following formula

x 1 0 ð Þ5m2 ffiffi ffi 3 p r52:13 w 1 0 ð Þ51=6 x 2 0 ð Þ5m53 w 2 0 ð Þ52=3 x 3 0 ð Þ5m1 ffiffi ffi 3 p r53:87 w 3 0 ð Þ51=6 (20)
It corresponds to the following expression for the mass particle-size initial distribution f(x, 0)

f ðx; 0Þ5 1 r ffiffiffiffiffi ffi 2p p e 2 x2m ð Þ 2 2r 2 (21)
One can thus explicitly define the analytical expression of the initial cumulative size distribution R(x, 0) as follows

Rðx; 0Þ5 ð 1 x f ðy; 0Þdy5 1 2 12erf x2m ffiffi ffi 2 p r ! (22)
The analytical solution of the mass particle-size distribution f(x, t) can be expressed as follows

f ðx; tÞ5 1 r ffiffiffiffiffi ffi 2p p e 2 x2m ð Þ 2 2r 2 e 2A 0 t x xmax ð Þ a 0 1A 0 t a 0 x a021 x max a 0 1 2 12erf x2m ffiffi ffi 2 p r ! e 2A0t x xmax ð Þ a 0 (23)
The model was validated by comparing the values of the first moments obtained by the DQMOM model with those derived from the analytical solution (Eq. 23). The results reported in Figures 1-3 correspond to the case where A 0 5 0.5 min 21 and a 0 5 2. Figure 1a represents the change in the first six moments as a function of time. Due to the definition of f(x, t), m 0 represents the overall solid mass and remains constant over time (see Eq. 13), which agrees with the mass conservation during the grinding process. The higher order moments decrease over time due to the size reduction. The simulated values closely approximated the analytical solution. The changes in the dimensionless moments, which were obtained by dividing the transient moments by their value at the initial time, are reported in Figure 1b on a semilog graph for the first minutes of the process.

In theory, one can return to the conclusions reached in the previous section given the selection and breakage kernels and explicitly define the transient changes in the moments

@m 1 @t 5S 1 52 A 0 a 0 11 1 x max a 0 m a 0 11 (24)
The following can be shown in a more general way

@m i @t 5S i 52i A 0 a 0 1i 1 x max a 0 m a01i (25)
An analytical expression of the moments as a function of time is not easily derived from these equations. However, an approximated expression of the first moments can be written over the first minutes of the process, as shown in Figure 1b m

k t ð Þ m k 0 ð Þ 5exp 2 t s k (26)
The time-scale of the local disintegration phenomena given by 1/A 0 can be compared with the time-scale of the process estimated from s k (for t 5 s k , the dimensionless moments are equal to 36.8%). The estimated values obtained for different values of A 0 are reported in Table 1. The time-scale of the process always significantly exceeded the time-scale of the disintegration process, irrespective of the value of A 0 . Furthermore, quasi-identical values of the parameters s k can determined as a function of the time-scale of the local disintegration phenomena. Specifically, the following relationships could be determined

s 1 % 11:2 A 0 ; s 2 % 7:2 A 0 ; s 3 % 5:9 A 0 ; s 4 % 5:0 A 0 ; s 5 % 4:4 A 0 ( 27
)

Monovariate DQMOM results and discussion

The changes in the abscissas and weights as a function time are, respectively, reported in Figures 2a,b. As expected, the abscissas decrease over time in relation to the size reduction process. Simply, if we assume that each couple (x i , w i ) characterizes a subpopulation of the entire size distribution, their change over time can be compared with the corresponding analytical particle-size distribution. Some examples are given in Figure 3. At each time considered, the abscissas belong to the size range, and the weights agree with the proportion of particles in the entire size distribution. Moreover, Figure 2b shows that the weight of the finer population first decreases during the first minutes of the process before increasing for longer times. Initially, the PSD is narrow, but it rapidly spreads over time as finer fragments are produced. These fragments eventually comprise the finer subpopulation, and the couple (x 1 , w 1 ) reproduces their change. Meanwhile, the weight of the coarser population first increases from its initial value before decreasing. This phenomenon results from the spreading of the PSD; the couple (x 3 , w 3 ) that represents the coarser population at any time expresses this trend. Initially, w 3 equals 1/6, it increases during the first minutes of the process while x 3 decreases (see e.g., Figure 3b at t 5 1 min) and then decreases as the coarse particles disintegrate. Thus, even though the weights and abscissas do not fully describe the particle-size distributions, these values closely approximate the corresponding analytical PSD.

Bivariate DQMOM Applied to Comminution Process General equations: homogeneous breakage

Two properties were considered for the solid phase of the bivariate case, the particle size and the particle strength, which was represented by the fracture energy, that is, the minimum energy needed to break the particle. The particlesize energy distribution f(x, e, t), which represents the mass fraction of particles characterized by size x and fracture energy e at instant t, was again approximated based on the N nodes of the Gaussian quadrature

f x; e; t ð Þ5 X N a51 w a t ð Þd x2x a ð Þd e2e a ð Þ (28)
Furthermore, the mixed moments were defined as follows

m kl ðtÞ5 ð 11 21
x k e l f ðx; e; tÞdx5

X N a51 w a x k a e l a (29)
In this bivariate case, the ODE system is based on 3N equations corresponding to the N nodal values of the weights, weighted sizes, and weighted energies @w a @t 5E wa (30)

@ w a x a ð Þ @t 5E wxa (31) @ w a e a ð Þ @t 5E wea (32)
The right hand sides of this ODE system are solutions of the algebraic system with 3N unknowns: E w1 , . . ., E wN , E wx1 , . . ., E wxN , E we1 , . . ., E weN . The generic equation to build the system is written as follows

X N a51 12k2l ð Þ x k a e l a E wa 1kx k21 a e l a E wxa 1lx k a e l21 a E wea  à 52 X N a51 w a x k a e l a a a 1 X N a51 a a b a kl w a (33)
The building of this algebraic system corresponds to the choice of 3N mixed moments or 3N couples (k, l), which are needed to close the system.

Application of bivariate DQMOM: case study

Recall the general mass balance equation for a batch grinding process

@f x; e; t ð Þ @t 52a
Crespo [START_REF] Crespo | Application of particle fracture energy distributions to ball milling kinetics[END_REF] recently developed both experimental and theoretical background for the comminution mechanisms related to this problem. In our bivariate approach, the selection function and the breakage function being based on both size and energy enable to introduce improved physics compared to monovariate case. The selection function a(x, e) describes the rate at which particles of size x and breakage energy e are broken. The breakage function b(x, e, x 0 , e 0 ) represents the probability to obtain fragments of properties (x, e) after a fracture event of particles characterized by size x 0 and breakage property e 0 .

Following Crespo, [START_REF] Crespo | Application of particle fracture energy distributions to ball milling kinetics[END_REF] the selection function can be considered as the product of two probabilities: the impact probability, which depends on the process conditions and the particle size, and the breakage probability after a stress event, which depends on the breakage energy. Assuming that the impact probability can be expressed by a power law function as in the previously considered monovariate case (Eq. 16), and that the breakage probability can be described by an exponential law, the overall selection function can be written as follows

aðx; eÞ5A 0 x x max a0 exp 2 e eà h i (35)
e* is a reference parameter that as the same dimensions as energy. Its value will be commented on in section Comparison of Monovariate and Bivariate DQMOM Methods.

Similarly, the breakage function can be decomposed into two contributions. Assuming that the size and energy of daughter particles do not depend on the mother energy after breakage, the following equation can be addressed

b x; e; x 0 ; e 0 ð Þ 5b 1 x; x 0 ð Þb 2 e ð jxÞ (36)
The first term breaks the usual probability into fragments of size x from mother particles of size x 0 , as expressed by Eq. 17. The second term expresses the energy spectrum of daughter particles conditioned by their size. The particle fracture energy can be evaluated from single-particle impact test, such as with the UFCL. [START_REF] King | Measurement of fracture energy during single-impact fracture[END_REF][START_REF] Tavares | Single-particle fracture under impact loading[END_REF] Such experiments conducted on individual particles lead to a large scatter of the data and are usually treated using statistical tools. The fracture probability can then be described by a log-normal distribution

P E ð jxÞ5 1 2 11erf ln E=E 50 ð Þ ffiffiffiffiffiffiffiffiffi ffi 2r E 2 p ! ( 37
)
where E is the mass-specific particle fracture energy (equal to the particle fracture energy divided by the mass of a particle, e/m p) . E 50 and r E are the median and geometric variance of the distribution, respectively. Tavares and King [START_REF] Tavares | Single-particle fracture under impact loading[END_REF] analyzed the effect of the particle size on the fracture energy for a variety of materials, and the following equation for the nominal size was proposed

E 50 5e 50 =m p 5E 1 11 x 0 x2x 00 / " # (38)
All parameters in this equation are constants that depend on the material properties. In our simulations, the following data obtained by Tavares and King [START_REF] Tavares | Single-particle fracture under impact loading[END_REF] for irregular particles of quartz were used: E 1 5 43.4 mJ/g, x 0 5 3.48 mm, x 00 5 5 mm, and / 5 1.61.

The fracture energy spectrum is then given by the following

b 2 ejx; t ð Þ5 @P ejx ð Þ @e 5 exp 2 ln e=e50 ð Þ ð Þ 2 2r E 2 h i e ffiffiffiffiffiffiffiffiffiffiffiffi 2pr E 2 p ( 39
)
With these assumptions, the fracture energy is a function of the particle size and the energy spectrum of the daughter particles does not depend on the mother energy spectrum. The second property of the solid is thus conditioned by the first one and the bivariate PBE could be solved with a conditional quadrature method of moments whose an algorithm was developed by Yuan and Fox. [START_REF] Yuan | Conditional quadrature method of moments for kinetic equations[END_REF] But the choice of DQMOM method rather than CQMOM method is driven by the prospect to derive a general approach of comminution problems that should be transferable to other kind of breakup mechanisms. Indeed, the conditions chosen for the case study are not universal in the comminution context and both properties have been treated here as if they were independent. Thus, the bivariate model could be used in more general situations, when the energy spectrum is not (or not only) governed by the particle size (fatigue phenomena, increase of particle strength, change in fracture mechanism during the process . . .) or even when another solid property, such as a morphological parameter, has to be considered.

Bivariate DQMOM validation and discussion

The PBE applied to a bidimensional comminution process does not have an analytical solution, such as the previously defined solution for the monovariate case. To validate the method, the change in certain mixed moments calculated by the bivariate DQMOM were compared to those calculated using a discrete solution of Eq. (34). After the discretization of the variables size and energy, Eq. ( 34) can be rewritten in the following form

@f i; j; t ð Þ @t 52f i; j; t ð Þa i; j ð Þ1 X i0<i X j0 f i 0 ; j 0 ; t ð Þa i 0 ; j 0 ð Þb i; j; i 0 ; j 0 ð Þ (40)
while considering discrete forms for the selection and the breakage functions such as

a i; j ð Þ5A 0 x i x max a 0 exp 2 e j eà h i (41) b i; j; i 0 ; j 0 ð Þ 5b 1 i; i 0 ð Þb 2 jji ð Þ5 a 0 x i0 x i x i0 a 0 21 b 2 jji ð Þ (42) with e j 5e 50 i ð Þexp ffiffiffiffiffiffiffiffiffi ffi 2r E 2 p erfinv 2P j 21 À Á h i (43)
e 50 i ð Þ5bqx 3 i E 1 11 x 0 x i 2x 00 / " # (44) b 2 jji ð Þ5 1 ffiffiffiffiffiffiffiffiffiffi ffi 2pr 2 E p 1 e j exp 2 log 2 e j =e 50 i ð Þ À Á 2r 2 E ! (45)
The problem was solved by considering one hundred size intervals between the minimum and maximum sizes (i.e., zero and 6 mm in the studied case) and one hundred energy intervals. The reference energy was defined to be equal to the median energy for the maximum particle size E 50 (x max 5 6 mm). The energy range was chosen to obtain a breakage probability between zero and 1 for each size class. The system of ordinary differential equations obtained after discretization was directly solved using the ode45 function of Matlab.

The discrete solution was used to validate the bivariate solution obtained with DQMOM for a number of nodes N 5 3. Indeed, following Fox and coworkers, [START_REF] Fox | Optimal moment sets for multivariate direct quadrature method of moments[END_REF][START_REF] Marchisio | Computational Models for Polydisperse Particulate and Multiphase Systems[END_REF] it is known that, for two solid properties, the minimum number of nodes to be considered to represent correctly the second-order moments is equal to 3. For N > 3, it is necessary to take into account mixed moments having higher orders which have no practical interest for the grinding process. The determination of moments with a high order on the base of experimental data could also lead to significant errors and this would constitute an important limitation for validation purpose.

If the number of nodes is equal to three, nine moments are needed to solve the system. In that case, it is known that there is no optimal moment set. [START_REF] Fox | Optimal moment sets for multivariate direct quadrature method of moments[END_REF][START_REF] Marchisio | Computational Models for Polydisperse Particulate and Multiphase Systems[END_REF] However, the following set of moments, that was chosen for the DQMOM resolution: (m 00 , m 10 , m 01 , m 20 , m 11 , m 02 , m 30 , m 03 , m 22 ), corresponds to a valid moment set, in accordance with the general rules given by Marchisio and Fox. [START_REF] Marchisio | Computational Models for Polydisperse Particulate and Multiphase Systems[END_REF] Indeed, all six moments of global order less than or equal to two were chosen (i.e., m 00 , m 10 , m 01 , m 20 , m 11 , m 02 ) and the remaining three moments were of global order three for two of them (m 30 and m 03 ) and of a global Order 4 for the last one (m 22 ). m 30 and m 03 were chosen to precisely predict the pure moments of order three, as suggested by Zucca et al. [START_REF] Zucca | Validation of bivariate DQMOM for nanoparticle processes simulation[END_REF] analysis. Pure moments of order two and three are needed to evaluate the d 32 (ratio of the pure moment of Order 3 on the pure moment of Order 2), which is a characteristic dimension often considered in particulate processes. Moreover, with the set of moments given above, both properties have been equally treated.

The comparison between certain mixed moments calculated by DQMOM with N 5 3 and those estimated by the discrete solution is illustrated in Figure 4. The initializing values were obtained by resolving the nonlinear system comprised of the 3N mixed moments using a Newton-Raphson iterative method. The following initial values were considered for the weights and abscissas: w 1 5 0.5, w 2 5 0.1, w 3 5 0.4; x 1 5 2.5, x 2 5 3.5, x 3 5 3.5; e 1 5 2.0, e 2 5 10, and e 3 5 2.7.

The models showed significant differences, especially for the mixed moments with a high order on the second variable (see Figure 4b). However, these errors were already present at the initial time. They are due to approximations in the double integral calculation within the discrete solution to determine the mixed moments. The precision of the discrete model could most likely be improved by increasing the number of size and/or energy intervals, which increases the computational time. Indeed, a standard computer obtains the DQMOM solution in a few minutes, whereas the calculation of the discrete model requires approximately two couple of weeks with a 150 3 150 grid. However, similar tendencies can be observed with both models. For example, both models calculated similar values for the Sauter diameter (Figure 5), which is defined by the ratio of the m 30 to the m 20 mixed moments.

Moreover, the DQMOM model provides access to the weights and abscissas as a function of time. These data are not directly correlated with the physical properties of the distributions but the couples (w a , x a ) and the particle-size distributions show good agreement. The PSD were determined as a function of time by using the discrete solution and calculating the mean values of the mass fractions over the energy range. Some examples are reported in Figure 6. The size abscissas belong to the particle-size range, and shifting the abscissas toward the lowest sizes is in concordance with the size reduction process, irrespective of the time. The increase in w 1 also accompanies the creation of a population constituted by fines particles.

Comparison of Monovariate and Bivariate DQMOM Methods

The monovariate and bivariate DQMOM methods can be directly compared by considering the mixed moments with a nil order for the second variable, m 10 , m 20 , m 30 , and so forth for the bivariate solution. The changes in the first six moments over time are reported in Figure 7a. As expected, the moments decrease more slowly with the bivariate solution. The second part of the selection function [exp(2e/e*)] acts as a decelerator for the breakage process. In the monovariate model, the energy limitation is implicitly contained in the kinetic prefactor A 0 . To approach identical conditions, we can assume that the breakage energy is independent of the particle size and constant over time. This constant is termed the critical energy and equal to e*. If the energy supplied during the process is higher than the critical energy, the breakage probability will equal 1. Conversely, it equals zero when the supplied energy is less than the critical energy. In this case, the energy breakage function can be expressed as follows

b ejx ð Þ5 @P ejx ð Þ @e 5d e2e à ð Þ (46)
The two models yield very similar values for the transient change of the six moments when the selection function is assumed to only depend on the particle size for both the monovariate and bivariate DQMOM models (Eq. 16) (Figure 7b). Figure 7c also shows that the monovariate and bivariate DQMOM solutions yielded similar changes in the Sauter diameter for the special case of a constant critical energy, whereas these values significantly differ for the general case.

The arbitrary grinding kinetics [A 0 5 0.5 min 21 and e* 5 e 50 (x max )] were considered on the Figure 7c. To prove the benefit of the bidimensionnal model over the 8a). Similarly, both models were compared for different values of the reference energy Eref while maintaining the kinetic parameter A 0 5 0.5 constant (Figure 8b). Decreasing the grinding parameter A 0 decreased the breakage kinetics and dampens the decrease of the Sauter diameter. The difference between the monovariate and bivariate solution was similar irrespective of the value of A 0 . In each case, the decrease in the Sauter diameter slowed when considering the energy barrier compared to the solution obtained with the monovariate model. A similar change was observed when considering the effect of the critical energy, e*, by dividing or multiplying it by 10. The value of the critical energy significantly affected the Sauter diameter. A low value of the critical energy implies that the material has a low strength (a high grindability) and can be ground by any amount of energy supplied to the mill. If the critical energy is independent of the particle size, the particle strength is not essential to the population balance. On the contrary, the benefit of the bivariate model compared to the classical one becomes significant for hard materials or if the energy needed for breakage depends on the particle size.

Conclusions

The DQMOM was used to solve the PBE applied to a homogeneous batch-grinding process. A property of the solid phase, the particle size, was first considered, and the model was validated by comparing the numerical results with an analytical solution of the PBE obtained by choosing a power law for the size selection function and a self-similarity expression for the breakage function. A bivariate model was then developed in which two properties were considered for the solid phase: the particle size and the minimum energy needed for their breakage. Expressions similar to those of the monovariate case were considered for the selection and breakage-size distribution in this case. Moreover, the energy spectra were assumed to be described by a log-normal probability function, and the breakage energy probability was assumed to be described by an exponential law. The bivariate DQMOM solution was very similar to a discrete solution of the PBE obtained by taking a finite number of size and energy intervals, but the computational time was significantly reduced. Finally, the benefit of the bivariate comminution model was proven by comparing the monovariate and bivariate solutions for a simplified case that adopted a critical minimum energy independent of size. The energy needed for breakage can thus be important for hard or heterogeneous materials for which the energy needed for breakage depends on the particle size. Generally speaking, our work shows the applicability of the DQMOM approach to solve a bidimensional PBE for a batch-grinding process. This model could be further developed to implement it in a CFD code to predict the solid properties while considering the heterogeneities of the impact energy inside the grinding chamber. Moreover, this work was focused on the particle size and energy for breakage, but other properties of the solid phase could also be considered. In particular, morphological characteristics could be considered to predict the change in the particle shape during a grinding process.

Literature Cited

Figure 1 .

 1 Figure 1. (a) Change in the first moments vs. timecomparison of the monovariate DQMOM solution (points) with the analytical solution (continuous lines) and (b) change in the dimensionless moments at the beginning of the process-comparison of the monovariate DQMOM solution (points) with an exponential law (lines).

Figure 2 .

 2 Figure 2. (a) Change in the abscissas-monovariate DQMOM and (b) change in the weightsmonovariate DQMOM.

Figure 3 .

 3 Figure 3. (a) Change in the PSD vs. time (analytical solution) and visualization of couples (w a , x a ) obtained by DQMOM monovariate-t 5 0 min, (b) t 5 1 min, (c) t 5 5 min, (d) t 5 20 min, and (e) t 5 100 min.

Figure 4 .

 4 Figure 4. (a) Change in certain mixed moments vs. time-comparison of monovariate DQMOM solution with the discrete solution and (b) change in certain mixed moments vs. timecomparison of bivariate DQMOM solution with the discrete solution.

Figure 5 .

 5 Figure 5. Change in the Sauter diameter vs. timecomparison of bivariate DQMOM solution with the discrete solution.

Figure 6 .

 6 Figure 6. (a) PSD obtained by the discrete solution taking mean values of the mass fractions over the energy range and visualization of the (w a , x a ) couples from bivariate DQMOM-t 5 1 min, (b) t 5 5 min, (c) t 5 20 min, and (d) t 5 100 min.

1.

  Lynch AJ. Mineral Crushing and Grinding Circuits. New York: Elsevier Scientific Publishing Company, 1977. 2. Prasher CL. Crushing and Grinding Process Handbook. New York: Wiley, 1987. 3. Chimwani N, Glasser D, Hildebrandt D, Metzger MJ, Mulenga FK. Determination of the milling parameters of a platinum group minerals ore to optimize product size distribution for flotation purposes. Miner Eng. 2013;43/44:67-78. 4. Hennart SLA, Wildeboer WJ, van Hee P, Meesters GMH. Identification of the grinding mechanisms and their origin in a stirred ball mill using population balances. Chem Eng Sci. 2009;69(19):4123-4130. 5. Fadda S, Cincotti A, Concas A, Pisu M, Cao G. Modelling breakage and reagglomeration during fine dry grinding in ball milling devices. Powder Technol. 2009;194(3):207-216.

Figure 7 .

 7 Figure 7. (a) Comparison of monovariate solution and bivariate DQMOM solution in the general case, (b) specific case (Constant critical energy), and (c) comparison of the Sauter diameter obtained with monovariate DQMOM solution and bivariate DQMOM solutions (general and special cases).

Figure 8 .

 8 Figure 8. (a) Comparison of the Sauter diameter obtained with monovariate and bivariate DQMOM solutions-effect of A 0 and (b) effect of Eref.

Table 1 .

 1 Comparison of the Time-Scale of the Local Break-Up Phenomena with the Time-Scale of the Grinding Process

	s 1	s 2	s 3	s 4	s 5

Kinetic Parameter

A 0 (min 21 ) Time-Scale Local Break-Up (min) Time-Scale of the Process (min)


domain: info.info-sc sdv chim

text:
Franc ¸ois Boulier

Franc ¸ois Fages

Ovidiu Radulescu

Satya S Samal

Andreas Schuppert

Werner Seiler

Thomas Sturm

Sebastian Walcher

Andreas Weber

The SYMBIONT Project: Symbolic Methods for Biological Networks

www.symbiont-project.org Overview Figure 1: SYMBIONT project overview SYMBIONT ranges from mathematics via computer science to systems biology, with a balanced team of researchers from those fields. At the present stage the project has a clear focus on fundamental research on mathematical methods and prototypes in software. Results are systematically benchmarked against models from computational biology databases. We summarize the motivation and aims for the project, and report on some existing results.

Main Objectives

Methods and Results

Computer algebra methods for solving systems of equations and inequalities For polynomial or rational vector fields, classical topics in the qualitative theory of ordinary differential equations such as the existence and stability of periodic solutions and equilibria can be reduced to real quantifier elimination problems. These questions are important for chemical reaction networks (CRN). However, they turn out difficult for medium size and large networks when using generic methods. Based on ideas from Stoichiometric Network Analysis (SNA) specifically virtual substitution and regular chain techniques allow the symbolic computation of Hopf bifurcation fixed points. Some crucial structural features of CRN are used to derive corresponding heuristics for quantifier elimination [START_REF] England | Symbolic versus numerical computation and visualization of parameter regions[END_REF]. "Complex Balancing," which is frequently encountered in CRNs, leads to parameterized binomial systems. For solving such systems, we are using polynomial time solution methods [START_REF] Grigoriev | Complexity of solving systems with few independent monomials and applications[END_REF]. Model order reduction often relies on quasi-steady state phenomena that lead to slow-fast variable decomposition. This is related to the classical singular perturbation scheme due to Tikhonov and Fenichel [START_REF] Goeke | Quasi-steady state -intuition, perturbation theory and algorithmic algebra[END_REF]. The aim is to determine critical manifolds, i.e., parameter values resulting in stationary points with positive dimension. Figure 4 illustrates typical results. In addition, we use scaling and more generally Lie symmetries for variable and parameter nondimensionalisation and identification of small parameters. Tropical geometry offers an alternative framework. For polynomial ODEs, we use tropical methods to rank monomial terms according to their orders of magnitude and to identify dominant terms.

Model order reduction and qualitative dynamics

When there is only one dominant term or when the dominant terms all have the same sign, the dynamics is fast and the system tends rapidly towards a region in phase space, where at least two dominant terms of opposite signs are equilibrated. We call a solution to the latter problem a tropical equilibration, which we can compute automatically (cf. [START_REF]A geometric method for model reduction of biochemical networks[END_REF] and Fig. 5). We show that tropical equilibrations can be algorithmically grouped into branches and put into correspondence with metastable dynamic regimes of chemical reaction networks (cf. Fig. 6a). Each tropical equilibration branch leads to a different reduced model. We aim to develop symbolic and numeric methods dedicated to the parameter estimation problem for models formulated by means of nonlinear integrodifferential equations (IDE). In particular, we address the computation of the model input-output equation and the numerical integration of IDE systems (cf. [START_REF] Boulier | Symbolic-numeric methods for nonlinear integro-differential modeling[END_REF] and Fig. 7). SYMBIONT eventually aims at investigating challenging models for systems biology and medicine applying the various methods described here. It will compute attractors, dominant subsystems, and metastable regimes as well as tipping points and transitions between such regimes (cf. Fig. 6b, Fig. 8). Our analysis could provide hints about the onset of a disease and support the mechanistic understanding of the disease process. For instance, network dynamic regimes are related to tumor cell phenotypes in the case of TGF-β signalling [START_REF]Geometric analysis of pathways dynamics[END_REF].

Applications in systems biology and medicine

Figure 2 :

 2 Figure 2: MAPK is an important signalling network. Reduced models serve as test cases.

Figure 3 :

 3 Figure 3: Schematic view of HoCoQ-an algorithm for computing Hopf bifurcations using convex coordinates and quantifier elimination.

Figure 4 :

 4 Figure 4: An instance where quasi-steady state reduction provides incorrect results (Michaelis-Menten system with small parameter k -1 ). The dashed line represents the concentration of product correctly, the dotted line represents the quasi-steady state approximation.

Figure 5 :

 5 Figure 5: Computation of tropical equilibrations using Constraint Logic Programming in BIOCHAM-4.

Figure 6 :Figure 7 :

 67 Figure 6: a) Flowchart indicating reconstruction and analysis of dynamical networks used as models for cell processes. Static networks are directed graphs; Dynamic networks are systems of polynomial ODEs; Branches B1, . . . , B10 of tropical equilibrations indicate metastable regimes; Transitions between metastable regimes form a non-deterministic automaton. b) Tropical equilibration branches provide proxys for metastable regimes of large biochemical networks.

Figure 8 :

 8 Figure 8: A representation of critical transitions in biological systems. The deviation of few "critical" parameters in the network (left) can shift the system's response resulting in the transition from "healthy" to "disease" state (right).

Challenges•

  A rigorous proof of the appearance of Hopf bifurcations, e.g., in the MAPK cascade • Using hidden toric structures of steady states in combination with linear conservation constraints • Model reduction for ODEs with sums of fractions • Additional constraints (e.g. stability constraints) to the tropical equilibration branches • Tipping points and dynamical regimes of large biological networks symbolically


Adrien Wion

Mathieu Bouet

Luigi Iannone

Vania Conan

Thales

Distributed Function Chaining with Anycast Routing

published or not. The documents may come

I. INTRODUCTION

Network services used to be built as an ordered set of physically wired hardware appliances that processed traffic for security or optimization purpose. With Network Functions Virtualization (NFV), middleboxes are more and more softwarebased running on top of virtualization-enabled equipment, thus allowing cost reduction and network flexibility. Nevertheless, with this new paradigm, new challenges have risen. Indeed, the service function chains are completely separated from the physical topology, and virtual functions are more ephemeral and dynamic in nature. Steering traffic through these sparsely located virtual entities, without compromising end-users sessions and Quality of Service (QoS), is therefore a complex challenge.

Despite the fact that Internet Service Providers critically rely on middleboxes for security and policy compliance [START_REF] Sherry | Making middleboxes someone else's problem: Network processing as a cloud service[END_REF], most of existing NFV management solutions rely on a logically centralized entity, generally named orchestrator. Such centralized approaches, as they require a holistic view of the whole network to perform service chaining, introduce control reactivity and resiliency (e.g., single point of failure) issues. Also, this may be quite costly for operators, since it requires the deployment of a whole new management and control infrastructure. In addition, the control part, which is meant to modify network configuration so to accommodate the orchestrator decisions, tend to be poorly interoperable with legacy appliances and is thus hard to deploy incrementally.

We believe that centralized orchestration for service function chaining is not necessary. The same functionalities can be provided in a distributed way by directly augmenting the network routing layer. In particular, we argue in this paper that it is possible to leverage on any Interior Gateway Protocol (IGP), anycast addressing, and any service chaining encapsulation, to construct a distributed service-aware distributed control-plane. We propose a modular architecture, showing that we do not need complex elements and we remain interoperable with legacy appliances. We also implemented such architecture, and early evaluation shows that our system successfully steers traffic through the intended service chain, distributing it over different instances according to available resources.

The reminder of this paper is organized as follows. First, we overview related work in Sec. II. Then, we introduce in Sec. III the main concept of our proposal: namely the service plane topology. We detail the system architecture in Sec. IV and the implementation in Sec. V. Early results supporting our proposal are presented in Sec. VI, while Sec. VII provides an agenda about research worth to be performed with respect to our proposal. Sec. VIII concludes the paper.

II. RELATED WORK

So far, NFV frameworks have been built on top of centralized cloud-based management system, which has simplified the use and implementation of resource allocation algorithms [START_REF] Palkar | E2: A framework for NFV applications[END_REF], [START_REF] Herrera | Resource allocation in NFV: A comprehensive survey[END_REF], [START_REF] Gember | Stratos: A network-aware orchestration layer for middleboxes in the cloud[END_REF], [START_REF] Bremler-Barr | OpenBox: a software-defined framework for developing, deploying, and managing network functions[END_REF]. For instance, Ghaznavi [START_REF] Ghaznavi | Distributed service function chaining[END_REF] proposes a centralized VNF (Virtual Network Functions) splitting and placement algorithm. Some solutions, such as Slick [START_REF] Anwer | Programming slick network functions[END_REF], go further proposing a programming language to define, on a central control point, high level policies. Such policies drive the decision being taken and enforced at runtime concerning chaining logic, flow forwarding, VNFs placement. Yet, such centralization, while simplifying VNF and path management, comes at the price of increased fragility (e.g., single point of failure, control loop delay, etc.) and high computation load on the central point.

While the Software-Defined Networking (SDN) paradigm helps in simplifying path management, as previously men-tioned, it still struggles to achieve traffic steering with finegrained forwarding rules. First, the size of forwarding state is limited by costly memory (TCAM). Second, dealing with forwarding rules installation when there are middleboxes (e.g., NAT -Network Address Translation -service) is another challenge to be tackled [START_REF] Fayazbakhsh | Flowtags: Enforcing network-wide policies in the presence of dynamic middlebox actions[END_REF], [START_REF] Gember | Stratos: A network-aware orchestration layer for middleboxes in the cloud[END_REF], [START_REF] Qazi | SIMPLEfying middlebox policy enforcement using SDN[END_REF].

Another approach consists in using encapsulation to convey traffic along a service chain. Recent work at the IETF (Internet Engineering Task Force) proposes Network Service Header (NSH) as a dedicated encapsulation header for service chaining [START_REF] Quinn | Network service header (NSH)[END_REF]. Segment Routing (SR) encapsulation has also been proposed to handle Service Function Chaining [START_REF] Abdelsalam | Implementation of virtual network function chaining through segment routing in a linux-based NFV infrastructure[END_REF]. It is based on a source routing model to steer traffic segment to segment. Recent work has also made the case for session protocol to build service overlay [START_REF] Zave | Dynamic service chaining with dysco[END_REF]. Even if the approach allows dynamic chaining, it relies on extending/modifying TCP, which, in an ossified Internet, is a hard task [START_REF] Honda | Is it still possible to extend TCP[END_REF].

Whether or not encapsulation is used, another main challenge in service function chaining is the coordination among flow path and middlebox state. Indeed, sometime VNF instances need to be created or reduced due to fluctuations in flow volume, migrated for resource optimization, or just recovered due to failure. Some solutions, like OpenNF [START_REF] Gember-Jacobson | OpenNF: Enabling innovation in network function control[END_REF] and Split/Merge [START_REF] Rajagopalan | Split/merge: System support for elastic execution in virtual middleboxes[END_REF], provide an open interface on middleboxes, so to allow a central coordinator to manage both forwarding and state. Dysco [START_REF] Zave | Dynamic service chaining with dysco[END_REF] proposes to solve these issues by consolidating forwarding and session state into middleboxes. Kablan et Al. [START_REF] Kablan | Stateless network functions: Breaking the tight coupling of state and processing[END_REF] instead try to avoid such state coordination problem by splitting each VNF into a stateless processing part and a consistent back-end data store.

All these works strongly rely on a centralized orchestration or were only used in such context. In the rest of the paper, we make the case for orchestrating service chaining in a distributed manner.

III. DISTRIBUTED ORCHESTRATION VIA NETWORK LAYER ROUTING AUGMENTATION

While so far service function chaining has relied on a holistic centralized orchestration to steer the traffic through sequences of virtual appliances, we believe that it can be done at the network layer routing in a distributed way.

Indeed, any network Interior Gateway Protocol (IGP) can be directly leveraged to convey the location, the type, and the necessary information associated to a virtual appliance and build an augmented network view.

Such a view, which we call the service plane topology, is formed by two different types of nodes. The first type is NFV nodes. NFV nodes are physical appliances that run the IGP and host virtual middleboxes (i.e., VNFs). NFV nodes can be datacenters, Points of Presence, or routers with VNF hosting capabilities. The second type of node, named VNF node, corresponds to the VNF instances themselves. NFV nodes can provide different types of service: Deep Packet Inspection (DPI), Firewalling, NAT, stream encoding etc. These virtual nodes run on top of NFV nodes. Since the NFV nodes that host them run the IGP, they can directly use control packets to share information on their VNF instances. This way, the Fig. 1: Network topology composed of 6 NFV nodes, with 3 of them hosting VNF instance (1a). The IGP views the two FW instances as a single entity, since they announce the same anycast IP address. A first flow (plain red line) is routed through the IDS and the top FW instance. A second flow (dashed blue line), arriving after the previous one, is then routed through the IDS and the bottom FW instance as the first FW instance is already loaded with the first flow.

VNF nodes are present in the IGP topology too. Consequently, each NFV node has a service plane view1 to not only take chaining decisions, but also VNF instantiation, scaling, or deletion decisions.

The main feature of VNF instances is the service they provide. We thus propose to leverage on anycast addressing to announce service functions (i.e., VNF instances) on the network. Different VNF instances, that are potentially hosted on different NFV nodes, but that provide the same service, will be announced by the same prefix. This way, each function is mapped to a prefix. The IGP cost to reach such prefix, announced by each NFV node, can be based on the VNF state, its available capacities, its load, or any other relevant information. 2 A link between two NFV nodes represents a topological distance (network cost), while a link to a VNF node describes some state of the VNF instance (VNF cost). Thus, a routing decision makes a tradeoff among these metrics and can be designed so as to balance the load, differentiate nodes or chains, etc. By applying the routing algorithm associated to the IGP to such a service plane topology, nodes can populate their routing table, which now becomes service-aware, since it includes routes toward all the available VNF-based services. Indeed, VNF instances providing the same service share the same prefix, hence they can be discriminated by the prefix's attribute(s).

Figure 1 illustrates with a toy example the approach we propose. Figure 1a represents the network topology constituted of NFV nodes. Each VNF instance of a given type is announced on the network with the same anycast address. In particular the two Firewall (FW) instances announce the same prefix: @IP f w . Flows have to be processed here by a unique chain: IDS + F W . The first flow is thus routed through the IDS instance and then through the top FW instance. Indeed, in this example, this VNF instance is at one hop from the NFV node that hosts the IDS instance. The NFV nodes that host the used VNFs advertise their neighbors with the new experienced load or any other relevant information. When the second flow arrives, the Firewall instance at the bottom is preferred, resulting in load balancing among the FW instances as well as dynamic path allocation for service function chains (Figure 1b). Notice that in Figure 1b, since the same prefix is announced but no adjacency is made, the flows that use a link to reach a service function (drawn as boxes) have to use the same link to go out of it. However, note as well that this link is only virtual, since it is the representation of the VNF instance in the IGP, but in reality is running directly on a NFV node.

Combining this augmented IGP with anycast addressing allows to fully benefit from what is already done at the network layer routing: network layer information exchange and route computation. Indeed, the NFV nodes can be considered as classic routers that compute the next hop(s) for the best path(s), depending on the metrics, to the different prefixes, which are in our case different network services.

As for any IGP, high level policies can be used to control the decision-making at each NFV node. They are common to all the nodes. Such policies include flow classification rules, to map traffic to the needed service chain. High level policies also concern routing decisions since all NFV nodes must share the same routing objectives. Based on the service plane topology, the NFV nodes can use the shortest path algorithm, or any other path computation algorithm, to choose which instance of the next VNF of the chain the flow will go through. Additionally, high level policies can define as well how to compute VNFs' IGP costs, stating which data to use and the function to translate such data in a cost.

To actually drive flows through the service chain they are associated to, we need to leverage on an encapsulation approach. In both the hop-by-hop model and the source routing model the encapsulation header should provide the information necessary to steer the flows through the correct sequence of VNFs. As such, the header should include i) part or all of the service chain identified at the classification step at the ingress of the network and ii) the next service step in this chain. For instance, in the example in Figure 1, the NFV node that hosts the IDS instance must have a mean to know that a packet belonging to a specific flow has been assigned to the service chain IDS + F W , that the next service to apply is F W , and which of the F W instances it actually has to go through.

IV. SYSTEM ARCHITECTURE

In this section, we describe the architecture of our system and design its main modules. A NFV node, as illustrated in Figure 2, is composed of a router providing underlay connectivity, a connector, which attaches the router to the different VNF instances, the VNFs themselves, providing the services, and a Distributed MANagement and Orchestration (D-MANO) component, which allows local autonomous management of the node.

Router: The router provides both underlay connectivity and participate in the network IGP. It exposes a control interface used by the D-MANO to inject or remove VNF anycast prefixes, announcing the services available on the node and the associated costs. This control interface is also used to get the IGP topology to build the service plane topology.

Connector: The connector allows dispatching traffic to the VNFs. It enforces chaining decisions as follows. It forwards incoming packets to the intended VNF instance, based on the encapsulation header. Once the packets have been processed, the VNF forwards them back to the connector, which enforces a forwarding decision toward the next VNF instance location (i.e., its connector) according to the service topology. These forwarding decisions are cached in the connector, indexed by a hash computed using flow-related information. The connector also exposes a control interface, used by the D-MANO, to populate the service-aware routing table and the mapping between service function and VNF instance unicast address. This information is used by the connector to enforce chaining decisions for outgoing traffic, and locally balance the load among the VNF instances that provide the same service (same prefix).

VNF: VNF instances process service flow packets according to the service they provide. Once a packet has been processed, the VNF instance updates the chaining encapsulation header to point to the next service. Each instance is monitored by the D-MANO.

Distributed MANO: The D-MANO controls and manages the other NFV node's components. It is configured with high level policies, which guide its autonomous orchestration decisions. It has two essential control functions (illustrated in Figure 2). The first one consists in monitoring VNF instances, deriving from them VNF costs, and then injecting such costs in the IGP, via the router. The second function consists in getting IGP information from the router to build the service plane topology, computing the service-aware routing table and then pushing it in the connector.

V. IMPLEMENTATION

We started to implement our proposed solution, which we describe in the present section. Furthermore, we include the technical choices we made for each component of the architecture described in the previous section.

A. System-Level Choices

Encapsulation Header: Our implementation is based on the Network Service Header (NSH) protocol to allow steering the traffic through the different services [START_REF] Quinn | Network service header (NSH)[END_REF]. Even if other encapsulations, such as Segment Routing v6 [START_REF] Abdelsalam | Implementation of virtual network function chaining through segment routing in a linux-based NFV infrastructure[END_REF], could have been used, our choice is motivated by the fact that NSH is an IETF standard explicitly designed for service chaining and is widely used in many opensource frameworks (e.g., [5], [4], [3],

[1]). In NSH, the Service Path Identifier (SPI) field uniquely identifies a set of abstract service functions (i.e., the Service Function Chain), while the Service Index (SI) points to the next function the packet has to be delivered to in the SPI set. NSH also provides extensible metadata fields that we leverage to convey the hash value used to consistently identify a flow along its chain. Such hash value is computed at the classification step with the 5-tuple of the original packet.

IGP: We build our implementation on top of an Open Shortest Path First (OSPF) underlay since this IGP is widely used and easily extensible thanks to opaque Link State Advertisements (LSA). Opaque LSAs are leveraged to propagate information about VNF instances and links. Even if flooding opaque LSAs increase control traffic overhead, it does not affect OSPF stability since they do not trigger shortest path algorithm computation. We thus define VNF opaque LSAs to convey 3 pieces of data: i) the anycast address of a VNF instance, ii) the associated VNF cost, and iii) the NSH endpoint IP address (i.e., the IP address of the next Connector). In our initial implementation, we choose to use a simple VNF metric: the remaining processing capacity of the VNF instance. The NFV nodes use the provided information to build a graph linking VNFs and NFV nodes, each link being weighted with the associated cost (Fig. 3). Thus, each NFV node is able to build the service plane topology based on the information shared via OSPF.

Service-aware path computation algorithm: We choose to use Weighted Cost MultiPath (WCMP) [START_REF] Zhou | WCMP: Weighted cost multipathing for improved fairness in data centers[END_REF] to compute nodes' service-aware routing table. It is particularly suited for our anycast-based approach as it allows balancing the traffic based on the VNF cost. As illustrated on Figure 3, we use network link costs and VNF costs to weight the paths to a VNF anycast prefix. In this example, we show the service topology as seen by node A. It is easy to see that the cost to reach the F W instance on node B is 30 and to reach the one on node D is 25. Such cost is used by WCMP to assign the flows on these instances. Since the VNF cost is regularly updated, WCMP adapts to the current load by distributing the traffic on the VNF instances that have the lower load (i.e., lower cost and thus higher WCMP weight). In Sec. VII we discuss more about the metrics.

B. Node-Level Choices

We build our NFV node using Linux and use network namespaces to isolate the components.

Router: In our implementation, we use FRRouting [2], an open source IP routing protocol suite, to implement our OSPF router. In particular, we use the OSPF API offered by FRRouting to mirror the Link-State Database (LSDB) in the D-MANO and to inject VNF opaque LSAs.

Connector: We implemented the connector logic in P4, a language for programming the dataplane [START_REF] Bosshart | P4: Programming protocol-independent packet processors[END_REF]. Our P4 code is run on the simple switch target [START_REF]P4 software switch[END_REF]. Its runtime CLI is exposed to the D-MANO to configure the switch and populate the WCMP table at runtime. VNF: VNFs are implemented as simple processes (using scapy [START_REF]Scapy[END_REF]) parsing incoming packets, decrement their NSH SI field, and forward them back to the connector. The focus of the initial implementation being on the different components of the proposed approach, we purposely choose simplistic VNFs for the time being. The Python psutil library enables us to monitor the resources used by the VNF processes.

D-MANO: The D-MANO has been implemented in Python. Its main loop runs as follows. First, it polls the resource use of the local VNF instances to build the related costs. The costs are then announced on the network with VNF opaque LSAs. Second, the D-MANO gets the VNF announces from its mirrored LSDB. With these data, it builds a service view (see Fig. 3b). Based on this topology, it computes WCMP weights and updates them on the connector.

VI. PRELIMINARY RESULTS

In this section we evaluate a simple scenario to show how we can achieve load balancing on different VNF instances of the same type by using the proposed solution.

We consider a network topology that looks like Figure 3b, except that we use one single generic service. Moreover, the link cost between Node A and Node B and between Node A and Node C are set to 160 and between Node A and Node D it is set to 30. All the VNF instances have the same initial capacity (i.e., VNF cost) set to 150. It is the maximum number of packets per second a VNF instance can process, normalized so to have the same range of values as the link costs. We use Mininet to emulate this topology [START_REF] Lantz | A network in a laptop: rapid prototyping for software-defined networks[END_REF].

Traffic has to be steered through one of the VNF instances and then toward the egress. We generate constant bit-rate flows on a source connected to Node A. Each flow lasts 50 seconds and consumes 2 of processing units at the VNFs. The arrival rate is of two flows per second. Our scenario evolves in 2 phases. Phase 1: only the VNFs on Node B and Node C are running. Phase 2: After 150 seconds, a third VNF is instantiated on at Node D, which leads to a redistribution of the traffic, since there are more VNF processing capacities now available in the service topology. of the experiment, the system load rises until it reaches its steady state. Note that, in this example, a measure of each VNF load is measured and advertised every 2s. We can see that, during the first phase, each VNF instance receives in average the same amount of traffic. Indeed, they do have the same network cost from the ingress point of view and the same initial VNF cost. Once Phase 2 starts, after the 50 seconds of transition, which lasts between t = 150s and t = 200s, a new steady state is reached. Now the VNFs on Node B and C, each process 40% of the traffic, while the VNF on Node D roughly processes 20%. This distribution of traffic corresponds to the WCMP weights that consider links' cost and VNFs' cost.

Figure 5 presents the mean traffic distribution on the instances for the steady state of the two phases of the scenario. They result from 20 runs of the experiment. We can observe that our solution is able to balance the load among the available VNFs. The mean and median loads are centered on the values we can compute from WCMP: 50/50% in Phase 1 and 40/40/20% in Phase 2. In addition, 50% of the loads are less than 3 points from the median value, while the max and min values are at most 10 points from it. Such limited variation shows that the system remains quite stable.

VII. RESEARCH AGENDA

Our preliminary results illustrate how service chaining can indeed be achieved by augmenting the network layer routing and applying high level policies. However, while opening interesting perspectives, it opens as well a number of questions. We overview them in this section.

Traffic Engineering Constraints: Forwarding traffic in service function chains, fulfilling both Service Level Agreements and cost minimization, is a hard task. Some long-lived flows require QoS guarantees (e.g., small delay for VoIP), while short ones may suffer from an initial latency in path computation (e.g., DNS query). We believe that best effort traffic and short lived flows are best handled by precomputed hop-by-hop routing decisions. Conversely, traffic requiring resource reservation would be best served using the source routing paradigm. Our service plane topology provides support for both approaches, and enables enforcing high-level policies. However, such hybrid scenarios and related tradeoff need further investigation.

VNF Metrics: In our approach, service-aware routing involves two different types of entity: namely network links and VNF instances. While assigning a cost to a link is straightforward and normal operation (based on bandwidth, latency etc.), evaluating the cost of a VNF instance is an open research area. On the one hand, such a cost may be based on a plethora of VNF state parameters [START_REF] Cao | Nfv-vital: A framework for characterizing the performance of virtual network functions[END_REF], [START_REF] Naik | NFVPerf: Online performance monitoring and bottleneck detection for NFV[END_REF]. On the other hand, the metric computation needs to be in the same order of magnitude of the links' metric, and, more importantly, it has to be additive, so to guarantee loop-free convergence even when taking into account multiple constraints [START_REF] Wang | Bandwidth-delay based routing algorithms[END_REF], [START_REF]Algorithms for finding paths with multiple constraints[END_REF].

VNF/Resource Management: Resource allocation is already a hard problem when using a centralized approach [START_REF] Herrera | Resource allocation in NFV: A comprehensive survey[END_REF]. Even if a distributed approach, like ours, improves the architecture resiliency and scalability, it adds coordination to the problem. In a distributed environment, each NFV node needs to take autonomous VNF provisioning decisions based on the exchanged information. Defining the needed information, their granularity, their update frequency, and the range of actions that each NFV node can take according to global resources availability remains to be explored. However, there is a great potential to use Machine Learning solutions that would make the network completely autonomous.

Service Modification: During a flow lifetime, the service chain it is associated to may be modified for numbers of reasons (e.g., a suspicious flow is redirected to a DPI). Such a service change implies a traffic redirection. In the source routing model, this modification could be easily handled since per flow state is concentrated at the edge. Conversely, in the hop-by-hop model, the chaining protocol should coordinate NFV nodes' state in order to modify the flow path (e.g., use Operations, Administration, and Maintenance for signaling). Existing work [START_REF] Zave | Dynamic service chaining with dysco[END_REF] identified challenges to keep end user sessions alive during reconfigurations. Coordination between VNF session state and routing decisions remains a challenging question to be explored.

Maintenance and Failure: To provide carrier grade network services, the impact of VNF unavailability on existing traffic should be minimized. With our approach maintenance can be easily handled through any existing loop-free graceful shutdown approach [START_REF] Francois | Disruption free topology reconfiguration in OSPF networks[END_REF]. Furthermore, some VNF state migration use-cases can be dealt locally with on NFV nodes [START_REF] Rajagopalan | Split/merge: System support for elastic execution in virtual middleboxes[END_REF], [START_REF] Kablan | Stateless network functions: Breaking the tight coupling of state and processing[END_REF], [START_REF] Woo | Elastic scaling of stateful network functions[END_REF]. However, VNF migration to a remote NFV node is more challenging. Indeed, NFV state migration has to be coordinated with service topology update. Certainly existing fail-over mechanisms and make-before-break approaches can be considered, yet, the design of such mechanism is an open research area.

Security: Distributing the service chaining decision raises some inherent security questions. To state if a VNF announce is valid or not, NFV nodes should trust each other. Trustworthiness can be solved by key distribution and initial authentication. Moreover, since the chaining protocol may convey sensitive information in its header or metadata, it may be useful to use encryption between authenticated NFV nodes. For instance, we could use IPsec [START_REF] Seo | Security Architecture for the Internet Protocol[END_REF] as a transport encapsulation between NFV nodes. In general, since we are augmenting the network routing layer, without revolutionizing it, there is quite a number of existing security solutions that can be considered to provide security in the service topology.

Multi-domain SFC: Even if multi-domain SFC would open new business opportunities, service providers are reluctant to share information related with their network. We believe that a distributed design can ease multi-domain orchestration. Defining an IGP routing logic to provide distributed SFC decisions is a first step for the design of multi-domain services. Indeed, the next step to be investigated is the use of interdomain routing, based on BGP [START_REF] Rekhter | A Border Gateway Protocol 4 (BGP-4)[END_REF] to provide chaining among different administrative entities, for instance based on the use of communities [START_REF] Li | BGP Communities Attribute[END_REF].

VIII. CONCLUSION

In this paper, we have made the case for orchestrating service chaining in a distributed manner. We proposed to augment the network layer routing by using anycast addressing for VNF so to build what we call the service topology, allowing embedding service chaining into routing. We designed an architecture based on this concept and implemented a first prototype. Early evaluation performed with our implementation shows that flows can be successfully driven through the chain of services according to available resources. Our approach sets itself apart from previous work, and as such it still needs to be thoroughly investigated. To this end we provide a research agenda highlighting the different aspects that need to be tackled. However, what comes out as well is quite promising and opens interesting perspectives.

  IGP logical view.

Fig. 2 :

 2 Fig. 2: NFV Node architecture. Doted arrows illustrate VNF routing control flow. Solid arrows show how VNFs state is monitored, transformed in a cost, which is then injected in the IGP.

  Service plane view at the Node A.

Fig. 3 :

 3 Fig. 3: Each node builds its service plane view (example at Node A on Fig. (b)) with the Network costs and VNF costs so as to compute the next hop(s).

Fig. 4 :

 4 Fig. 4: Traffic distribution over time on the VNF instances. During the first 150s only two VNFs are running. At t = 150s, a third VNF is instantiated.

Figure 4

 4 presents the traffic distribution over time on the VNF instances. Since each flow lasts 50s, during the first 50s

Fig. 5 :

 5 Fig. 5: Boxplot of the traffic distribution on the VNF instances during the two phases.

We use the terms service plane topology and service plane view interchangeably.

2 In Sec. VII, we discuss more about how to calculate such a metric in a meaningful and rigorous way, so to guarantee loop-free routing convergence.