Dr Marie Debacq
email: marie.debacq@cnam.fr

Phahath Thammavong

Stéphane Vitu

Denis Ablitzer

Jean-Léon Houzelot

Fabrice Patisson

A HYDRODYNAMIC MODEL FOR FLIGHTED ROTARY KILNS USED FOR THE CONVERSION OF COHESIVE URANIUM POWDERS

Keywords: rotating cylinder, cohesive powder, lifter, transverse distribution, hydrodynamics, uranium

The axial and transverse solid distribution of a cohesive powder in flighted rotary kilns was calculated from an original method based on the lifter discharge law previously determined, geometrical calculations and supplementary measurements (in particular the fall velocity of the powder at high temperature). At a given axial position, the quantity of powder falling through the gas, the total quantity of flighted powder as well as the quantity remaining in the bulk are thus calculated; the axial bed depth profile is obtained from the Saeman-Kramers-Afacan model.

Highlights

► Transverse flow of cohesive powder in flighted rotary kiln was modelled. ► Axial bed depth profile was calculated from Saeman-Kramers-Afacan model. ► Solid distribution was calculated from geometrical calculations and measurements. ► The fall velocity of the cohesive powder was measured at high temperature. ► The 'snapshot' obtained allows to calculate thermal transfer areas and mass transfer.

INTRODUCTION

For over a century [START_REF] Meade | The Development of the Rotary Kiln and Its Application to Various Chemical and Metallurgical Processes[END_REF], rotary kilns have been widely used in the inorganic chemical industry. They are a key element in the production processes of cement [START_REF] Mujumdar | Simulation of Rotary Cement Kilns Using a One-Dimensional Model[END_REF], lime and pigments (titanium dioxide calcination) and in extractive metallurgy for the reduction of ore [START_REF] Boateng | Rotary Kilns[END_REF]. Rotary kilns are well suited for the conversion of coal and coke [START_REF] Lebas | Experimental study of residence time, particle movement and bed depth profile in rotary kilns[END_REF], for drying [START_REF] Shahhosseini | A simple dynamic model for solid transport in rotary dryers[END_REF], for the pyrolysis of solid waste [START_REF] Descoins | Solid transport in a pyrolysis pilot-scale rotary kiln: preliminary results-stationary and dynamic results[END_REF][START_REF] Marias | Modelling of a rotary kiln for the pyrolysis of aluminium waste[END_REF], and are also used for several conversions in the manufacture of nuclear fuel [START_REF] Debacq | Étude et modélisation des fours tournants de défluoration et réduction du difluorure d'uranyle[END_REF][START_REF] Thammavong | Experimental apparatus for studying heat transfer in externally heated rotary kilns[END_REF]. Rotary kilns used for the conversion of enriched or depleted uranium are equipped with lifters and are externally and electrically heated, but they differ in size and nominal operating conditions, as can be seen in Table 1: the type of rotary kiln used for depleted uranium conversion is henceforth named "kiln #1"; the type of conversion kiln used for enriched uranium is named "kiln #2". In these rotary kilns, uranyl difluoride (UO2F2) is converted into uranium oxides: U3O8 in the case of depleted uranium in kiln #1; UO2 in the case of enriched uranium in kiln #2 (U3O8 is an intermediate during the conversion from UO2F2 to UO2). In this type of rotary kilns, heat and mass transfer are strongly governed by hydrodynamics, which determines -among other characteristics-the quality and duration of gas/solid contact. In order to correctly model the operation of the conversion kilns, it is essential to be able to calculate precisely the solid distribution at any point in the rotary kiln.

While there is an abundant literature on the hydrodynamics of unflighted rotary kilns (see for example [START_REF] Henein | Experimental studies of transverse bed motion in rotary kilns[END_REF][START_REF] Mellmann | The transverse motion of solids in rotating cylinders-forms of motion and transition behaviour[END_REF], only a few studies concern flighted rotary kilns (see for example [START_REF] Van Puyvelde | Modelling the hold up of lifters in rotary dryers[END_REF]. The hydrodynamics can be studied globally: the effect of lifters on the total hold-up (see for example [START_REF] Friedman | Studies in rotary drying, part I. Holdup and dusting[END_REF], the time of passage (see for example [START_REF] Prutton | Factors influencing the performance of rotary dryers[END_REF], or the mean residence time (see for example [START_REF] Britton | A physical description of solids transport in flighted rotary dryers[END_REF]. A few studies propose calculation of the amount of powder falling through the gas (sometimes called airborne solid), usually in rotary dryers fitted with a large number of lifters (often a dozen, and in some cases several dozen), in which this amount is high. In such cases, it is not worth taking into account the variation of the bulk depth at the bottom of the rotary kiln with the axial position, and the calculated length of particle fall depends solely on the angular position (see for example [START_REF] Fernandes | Modeling and experimental study of hydrodynamic and drying characteristics of an industrial rotary dryer[END_REF]. Moreover, the total flighted powder fraction is mostly not calculated since, in the case of drying with hot gas circulating through the rotary kiln, the solid in the flights plays practically no role in the drying and heating processes. Thus, in these studies, only airborne (and sometimes dense phase hold-up) is quantified [START_REF] Schofield | Rotary Driers and Coolers for Granular Fertilizers[END_REF]Matchett andBaker, 1987, 1988;[START_REF] Sherritt | The Movement of Solids Through Flighted Rotating Drums. Part I: Model Formulation[END_REF][START_REF] Sherritt | The Movement of Solids Through Flighted Rotating Drums. Part II: Solids-gas interaction and model validation[END_REF][START_REF] Wang | A fundamental study of particle transport through rotary dryers for flight design and system optimisation[END_REF]Ajayi andSheehan, 2012a, 2012b).

In the rotary kilns studied here, the number of lifters is much lower than in dryers (4 or 6 depending on the type of conversion kiln). A substantial amount of solid may therefore remain in the bulk; the quantity of powder falling through the gas is extremely small (although it plays a very important role in the chemical reaction) and the solid in the lifters is likely to react and to play a role in the thermal exchanges with the wall which is externally heated. It has therefore been necessary to develop a novel model to quantify the various amounts of solid.

In a previous study [START_REF] Debacq | Transverse motion of cohesive powders in flighted rotary kilns: experimental study of unloading at ambient and high temperatures. accepted for publication in[END_REF], the discharge law of cohesive powders in flighted rotary kilns was obtained both experimentally and theoretically (i.e. from geometrical considerations), including high temperature measurements. Sunkara et al. (2013aSunkara et al. ( , 2013b) ) published a similar method for the calculation of the discharge characteristics of rectangular flights, but with experiments on non-cohesive particles and at room temperature.

In the present study a model is proposed to calculate the transverse solid distribution in a flighted rotary kiln; this model takes into account the change in the bulk density and therefore in the volumetric flow rate of the solid induced by the chemical reactions. The transverse solid distribution can be calculated by applying the lifter discharge law: for each axial position in the conversion kiln, a 'snapshot' can be determined, representative of the powder fraction F in the lifters, the powder fraction G falling through the gas and the powder fraction B in the bulk (Fig. 0 / Graphical abstract). This is particularly useful to study mass transfer between these different powder fractions and gaseous reagents. This 'snapshot' also makes it possible to calculate all the thermal transfer areas, such as the area between solid and wall or between solid and gas, etc.

HYDRODYNAMIC MODEL

Axial profile of bed depth

The axial profile of the bed depth (h versus x) was calculated using laws available in the literature: the bed depth at different axial positions along rotary kilns without inserts in rolling mode can be calculated from the differential equation ( 1), proposed by [START_REF] Saeman | Passage of solids through rotary kilns: factors affecting time of passage[END_REF], [START_REF] Vàhl | Transport of solids through horizontal rotary cylinders[END_REF] and [START_REF] Kramers | The passage of granular solids through inclined rotary kilns[END_REF]:

2 / 3 2 3 v R h R h 2 R 4 tan Q 3 cos tan dx dh                              (1)
From observations made in the pilot plant described previously [START_REF] Thammavong | Experimental apparatus for studying heat transfer in externally heated rotary kilns[END_REF], the presence of a small number of lifters does not modify the shape of the bed depth profile proposed by Saeman.

As suggested by [START_REF] Afacan | Analysis of granule breakage in a rotary mixing drum: Experimental study and distinct element analysis[END_REF], the numerical solution of the differential equation ( 1) can be found using a fourth-order Runge-Kutta method. These authors also suggested using a bed depth equal to the height of a lifter plus twice the diameter of the particles (hp + 2.dp) as the boundary condition at the cylinder exit. More recently Specht et al. (2001) found in smooth rotary kilns (without lifters or dam) that the solid depth at the discharge is several times higher than the particle diameter (and not only twice as high, as suggested by [START_REF] Afacan | Analysis of granule breakage in a rotary mixing drum: Experimental study and distinct element analysis[END_REF]. Since the particles studied here are very small (Table 2), their size is negligible compared to the height of a lifter (Table 1); after conducting tests, it appeared unnecessary to take their size into account in the boundary condition, thus only hp is used.

The local kiln hold-up Z(x) for each axial position x can be calculated using relation (2) from the local bed depth h(x) obtained by integrating the differential equation (1) using a fourth-order Runge-Kutta method.

                                             2 / 1 2 1 R x h R x h 2 R x h 1 R x h 1 cos x Z
(2)

In the case studied here, the volumetric flow rate Qv is not constant since the bulk density of the solid changes considerably during its conversion along the cylinder, going from UO2F2 at the rotary kiln entry to U3O8 or UO2 at the exit, depending on the type of conversion kiln considered (see densities in Table 2). The bulk density or the angle of repose considered at a given axial position is the mean -weighted by the composition of the solid phase at this axial position-of the bulk densities or the angles of repose of the different types of powders (given in Table 2). The composition of the solid phase is obtained by calculations involving a mass transfer model [START_REF] Debacq | Étude et modélisation des fours tournants de défluoration et réduction du difluorure d'uranyle[END_REF] that will be detailed in another paper. The hydrodynamic model and the mass transfer model are of course coupled, and are also coupled to the thermal model. Since it is not possible to solve these coupled models in parallel, an iterative calculation is required to achieve convergence on the conversion, local hold-up and temperature profiles [START_REF] Debacq | Étude et modélisation des fours tournants de défluoration et réduction du difluorure d'uranyle[END_REF] that will be the object of a third paper.

To summarize, the constant parameters in equations ( 1) and (2) are the kiln inclination angle , the rotational speed  and the inside kiln radius R. The dynamic angle of repose  and the volumetric flow rate Qv are calculated from the composition of the solid phase at each axial position. The bed depth profile h versus x is obtained from equation (1) with the boundary condition h = hp at the kiln exit. Finally the local kiln hold-up Z is calculated from equation (2) and the h(x) value at each axial position.

Transverse solid distribution

As mentioned in the introduction, some studies reported in the literature make it possible to calculate the amount of powder falling through the gas and sometimes the dense phase hold-up. Several calculations (Matchett andBaker, 1987, 1988;[START_REF] Sherritt | The Movement of Solids Through Flighted Rotating Drums. Part I: Model Formulation[END_REF][START_REF] Blumberg | Transverse bed motion and convective mass transfer in rotary cylinders. Part 2. With lifters[END_REF] are based on two hypotheses which are not verified in the present case: first, the absence of solid in the bulk (the calculation of the length of particle fall from a lifter is then simplified) and second, the calculation of the fall velocity as that of a single particle in free fall at a low Reynolds number; it is then necessary to be able to define the diameter of the particle. The particles considered in the present study, however, are not only particularly cohesive (Carr's index is more than 0.4 and the Hausner ratio is more than 1.4, as can be seen in Table 2) but are also formed of more or less coherent aggregates of small-sized elementary particles and therefore have a broad granular distribution (see Table 2); this means that it is not possible to define an appropriate particle size for the terminal fall velocity as is necessary for classical evaluation methods. [START_REF] Sherritt | The Movement of Solids Through Flighted Rotating Drums. Part II: Solids-gas interaction and model validation[END_REF] proposed an iterative procedure to calculate the change in bulk depth along the rotary kiln; however their method to calculate the hold-up in a lifter can only be applied to freeflowing granular materials and is thus not applicable to our case. Ajayi andSheehan (2012a, 2012b) recently published experimental studies using image analysis techniques to quantify the airborne solids, again for free-flowing and coarse particles, while in this work fine-grained cohesive powders are considered.

In order to find the transverse solid distribution (the 'snapshot' mentioned in the introduction and shown on Fig. 0 / Graphical abstract), i.e. not only the amount of powder falling through the gas but also the total flighted powder and the amount of powder in the bulk (which varies along the conversion kiln), the method of calculation described in the following sections is proposed.

Powder fraction falling through the gas

To calculate the powder fraction falling through the gas, it is necessary to determine the length of fall of the powder (see the Appendix) and therefore to know the quantity of solid present in the bulk. This quantity will be deduced (see section 2.2.3) from the local kiln hold-up calculated using the Saeman-Kramers-Afacan model (as explained in section 2.1), the flighted powder fraction (calculated by the procedure described in section 2.2.2) and the powder fraction falling through the gas which is sought to be determined here. An iterative calculation is therefore needed until convergence on the value of the three fractions is achieved; convergence between the local hold-up profile (values of Z for each axial position) and conversion of the powder is ensured in the global model.

The powder fraction G falling through the gas is determined by a calculation procedure involving the powder fraction q falling from one lifter between time t and t+t. For a given rotation speed , q can be calculated using equation (3) which also takes into account the slope  of the lifter discharge law  = .(F-) determined in a previous study (see notations on Fig. 1).

t 2 q       (3)
The parameters of the discharge law were measured in the operating conditions of uranium conversion kilns, i.e. at high temperatures [START_REF] Debacq | Transverse motion of cohesive powders in flighted rotary kilns: experimental study of unloading at ambient and high temperatures. accepted for publication in[END_REF]; in addition to the slope , these measurements provided the values of the following characteristic angles (the notations used are indicated on Fig. 1):  D: lifter position for which the powder it contains is no longer in contact with the powder remaining in the bulk;  I: start of 'effective' discharge, with reference to the mass transfer between solid and gaseous reagents;  F: lifter position at which the last powder particle falls.

For the calculations described in the present paper, the slope and characteristic angles considered at a given axial position are the means -weighted by the composition of the solid phase at this axial positionof the characteristic slopes and angles of the different types of powders (see values in section 3). The relevant angles for the calculation of G are I and F.

The powder fraction q falls during a time tfall which is equal to the length of fall hfall (see notations on Fig. 1 and the demonstration in the Appendix) divided by the fall velocity vfall and by the cosine of the kiln inclination angle , according to equation (4). It was observed that in a cross-section of the rotary kiln, the powder falls almost vertically from the lifters at the rotational speeds used, as shown on Fig. 2; it is therefore unnecessary to take any deviation of the particle path in the transverse plane into account. However, the kiln inclination itself can be taken into account via cosine  in equation ( 4). The fall velocity is determined experimentally at high temperature (see section 3). After falling, the powder fraction q therefore remains present in the gas phase between time t and t+tfall.

                                                    2 sin R cos tan h h R v cos tan v cos 2 cos R cos h sin h R t p o fall fall p o fall (4)
The angle  in equation ( 4) is related to the fraction B of powder remaining in the bulk by equation ( 5). As this relation is implicit in , and since an iteration loop would have further increased the already large number of iterative calculations in the global rotary kiln model [START_REF] Debacq | Étude et modélisation des fours tournants de défluoration et réduction du difluorure d'uranyle[END_REF], a numerical approximation is used to optimize the calculation time: equation ( 6) for B < 0.05 and equation ( 7) for B  0.05. The agreement between equation ( 5) and the approximations ( 6) and ( 7) is shown on Fig. 3.

      2 sin B (5)   1 ln 4870 . 5 4457 . 0 B      (6)   .1219 1 log 12.5965 6518 . 0 B      (7)
The value of  also makes it possible to recalculate the detachment angle D corresponding to the quantity of solid in the bulk, using equation ( 8).

                             2 2 for , 2 2 2 2 for , 2 2 D (8)
The powder fraction G falling through the gas at a given axial position is calculated by the following procedure:

1. Depending on the kiln rotational speed, the position of each lifter is calculated at each instant (a time-step t = 0.01 s, i.e. 0.12° at 2 rpm for example was chosen) during (1 + 1/n) revolutions (where n is the number of lifters inside the kiln); lifter number 1 is conventionally positioned at an angle  = 0 at time t = 0; the other lifters have an incremented number counterclockwise. 2. At a given axial position (and therefore for a local kiln hold-up and characteristics of the discharge law (, I and F) given by the composition of the solid phase), and for each instant, the lifters that are in 'effective' discharge are selected (i.e. angular position between I and F); this number is usually 2 or 3 for kiln #1 for example (this is the number of arrows in each box on Fig. 4). The quantity of powder q which falls from each lifter during the time-step t is then calculated according to equation (3), as well as the fall time needed for each 'clump' of powder to fall from each lifter in 'effective' discharge, using equations ( 4) and ( 6) or (7). 3. From this, the number of time-steps during which each 'clump' of powder is airborne is deduced.

The total powder fraction falling through the gas is then incremented at each subsequent timestep, as shown on Fig. 4. For example, for lifter 1 in the first time-step (beige arrow on colored version of Fig. 4 or gray arrow on black and white version of Fig. 4, on the left of Fig. 4), if it takes five time-steps for the powder fraction falling from lifter number 1 to reach the bulk in the bottom of the kiln, the corresponding quantity of powder (small beige or gray rectangles) is airborne for five time-steps; if it takes three time-steps for the powder fraction falling from lifter number 2 to reach the bulk in the bottom of the kiln (brown or black arrow on Fig. 4), the corresponding quantity of powder (small brown or black rectangles) is airborne for three timesteps, and so on for each lifter at each time-step.

4. Once the steady state regime has been reached (usually in fact in much less than one revolution), the instantaneous total powder fraction falling through the gas oscillates around an average value, as shown on Fig. 5, with a frequency corresponding to the time interval between the passage of two successive lifters at a given angular position  (i.e. 1/n revolution). That is why the powder fraction G falling through the gas is calculated by averaging over the instants between 1 and (1 + 1/n) revolutions (see Fig. 4). 5. Steps 2 to 4 of this procedure are repeated for different axial positions along the cylinder. Note that the length of the discretization segments of the cylinder always corresponds to a time interval at least equal to, and generally greater than, the time required for the kiln to perform 1 + 1/n revolutions [START_REF] Debacq | Étude et modélisation des fours tournants de défluoration et réduction du difluorure d'uranyle[END_REF].

To summarize, the constant parameters in equations ( 3) to ( 8) are the kiln inclination angle , the rotational speed , the inside kiln radius R, the dimensions of the lifters ho and hp and the fall velocity vfall. The dynamic angle of repose , the slope of discharge law  and the characteristic angles I and F are calculated from the composition of the solid phase at each axial position. This allows to calculate q from equation (3) and to determine the lifters that are in discharge in step 2 of the previous procedure. The duration of the fall can then be calculated from equation ( 4), in which the filling angle  is obtained from approximation ( 6) or ( 7) depending on the value of B in the previous iteration. The value of the filling angle  also allows to calculate D from equation ( 8); this value is necessary for the determination of F.

Step 3 of the procedure can then be carried out (Fig. 4). Once the time interval has been investigated in full, it is possible to calculate G for the considered axial position (Fig. 5) according to step 4 of the procedure. As mentioned in section 2.1, since the determinations of G, F and B are not independent, an iterative calculation is necessary to achieve convergence between these three values: an estimated value of B is needed to determine G; since the calculation of G uses the filling angle , which is itself a function of B (see equation ( 5)), this is what is used to test convergence.

Flighted powder fraction

The total flighted powder fraction F is determined using a calculation procedure that involves only the lifter discharge law  = .(F-) determined in a previous work [START_REF] Debacq | Transverse motion of cohesive powders in flighted rotary kilns: experimental study of unloading at ambient and high temperatures. accepted for publication in[END_REF] (see the notations on Fig. 1). The appropriate angles for the calculation of F are D and F.

The procedure is as follows:

1. Depending on the kiln rotational speed, the position of each lifter is calculated at each time-step (a time-step of 0.01 s was also chosen here) during 1/n revolutions (where n is again the number of lifters); lifter number 1 is again conventionally positioned at an angle  = 0 at time t = 0. 2. At a given axial position (and therefore for a local kiln hold-up and characteristics of the discharge law (, D and F) given by the composition of the solid phase), and for each instant, the powder fraction in each lifter is calculated using the discharge law. 3. In this procedure, no steady state settling is needed; the total flighted powder fraction oscillates as a function of angular position around the average value F, as shown on Fig. 6. 4. Knowing the average fraction F, the 'typical' lifter angular position that it corresponds to is determined, by taking the values obtained at all angular positions. The average whole number of lifters that are detached ( > D) and that contain solid ( < F) in this typical position is also determined. The angular position and number are useful for the mass transfer and thermal models [START_REF] Debacq | Étude et modélisation des fours tournants de défluoration et réduction du difluorure d'uranyle[END_REF]. 5. Steps 2 to 4 of this procedure are repeated for different axial positions along the conversion kiln.

Powder fraction in the bulk

The fraction B of solid in the bulk is obtained for each axial position by a simple subtraction (9). )

x ( ) x ( ) x ( Z ) x ( F G B       (9)

EXPERIMENTAL MEASUREMENTS

The dimensions and nominal operating conditions of industrial kilns studied here are given in Table 1: in the type of rotary kiln used for depleted uranium conversion and named "kiln #1", uranyl difluoride (UO2F2) is converted into uranium sesquioxide (U3O8); in the type of conversion kiln used for enriched uranium and named "kiln #2", uranyl difluoride (UO2F2) is converted into uranium dioxide (UO2) (U3O8 is an intermediate during the conversion from UO2F2 to UO2). These two types of kiln differ in dimensions and nominal operating conditions, but the same model is used to describe them. The properties of the different powders are given in Table 2. It can be seen that the oxides are clearly cohesive (Carr's index is more than 0.4 and the Hausner ratio is more than 1.4); UO2F2 is on the borderline of highly cohesive powders.

The experimental determination of the discharge law was described elsewhere [START_REF] Debacq | Transverse motion of cohesive powders in flighted rotary kilns: experimental study of unloading at ambient and high temperatures. accepted for publication in[END_REF]. The values found are:

 slope of discharge law:  = 0.00012°-1 for UO2F2 and 0.00035°-1 for oxides;  detachment angle: D = -18° for all powders (rough value measured for mean hold-up; D can be calculated from the filling angle through equation ( 8));  angle at the start of effective discharge ('effective' with reference to mass transfer between solid and gaseous reagents): I = -2° for UO2F2 and 80° for oxides at high temperatures;  angle at the end of discharge: F = 131° for all powders. As discussed previously [START_REF] Debacq | Transverse motion of cohesive powders in flighted rotary kilns: experimental study of unloading at ambient and high temperatures. accepted for publication in[END_REF], at room temperature the different powders behave in the same way (same values of slope and angles). At high temperatures (those of the process), however, while the behavior of UO2F2 remains similar to that under cold conditions, oxides have a very different behavior and powder discharge from the lifters is much like a sandcastle being unmolded from a bucket on a beach. The values of I and F used in the calculations presented in the present paper are those given above; the value of D is calculated at each axial position using equation (8). Fig. 2 is obtained by overlaying on the same image the powder curtains falling from the lifter tips from various angular positions. They were obtained from the images used to determine the discharge law [START_REF] Debacq | Transverse motion of cohesive powders in flighted rotary kilns: experimental study of unloading at ambient and high temperatures. accepted for publication in[END_REF]. In a cross-section of the drum, because of the rotational speed, the particles do not fall exactly vertically. Since the vertical deviation is weak, however, it will be neglected here. Several studies also report quasi-vertical particle paths [START_REF] Prutton | Factors influencing the performance of rotary dryers[END_REF][START_REF] Kelly | Dynamics of granular material in rotary dryers and coolers[END_REF][START_REF] Britton | A physical description of solids transport in flighted rotary dryers[END_REF][START_REF] Afacan | Analysis of granule breakage in a rotary mixing drum: Experimental study and distinct element analysis[END_REF]Ajayi andSheehan, 2012a, 2012b;Sunkara et al., 2013aSunkara et al., , 2013b)).

As mentioned in section 2.2, the classical methods used to calculate the fall velocity of the particle curtains cannot be applied to the powders considered in the present study, as they are particularly cohesive and have a broad granular distribution (Table 2). It therefore seemed preferable to determine the velocity experimentally. The setup designed to determine the discharge laws at high temperature [START_REF] Debacq | Transverse motion of cohesive powders in flighted rotary kilns: experimental study of unloading at ambient and high temperatures. accepted for publication in[END_REF] was reused for the present measurements in stagnant air. A graph plate was placed vertically behind the right-angled flight. The flight was rotated manually at a speed close to that encountered in industrial kilns. The recorded films were used to measure the velocity of the front ahead and behind the powder, which fell in 'clumps' through the gas, as described elsewhere [START_REF] Debacq | Transverse motion of cohesive powders in flighted rotary kilns: experimental study of unloading at ambient and high temperatures. accepted for publication in[END_REF]. The fall velocity measured -in temperature conditions around 700-800°C-was 1 m s -1 , with a range over the seven tests and five measurements per test, between 0.8 and 1.2 m s -1 . There was no significant difference between the velocity of the upstream and downstream powder fronts. This value is similar to that obtained by [START_REF] Kelly | Dynamics of granular material in rotary dryers and coolers[END_REF].

RESULTS AND DISCUSSION

Flighted powder fraction

If the time-step selected for the calculation procedure is sufficiently small, the rotational speed has no influence on F. Fig. 7a shows the variations in the F fraction (in green on the colored version of Fig. 7a / in black on the black and white version of Fig. 7a) as a function of the local kiln hold-up: at low hold-up, F increases with Z up to a maximum value (F=0.0284 for Z=3.5% in the conditions of Fig. 7a); it then decreases slightly until a minimum (F=0.0239 around Z=13.5% in the conditions of Fig. 7a), when the solid begins to collect in the bulk and D increases; beyond this value, F increases again since the lifters can be almost completely filled at small angular positions, as observed in the drums [START_REF] Debacq | Transverse motion of cohesive powders in flighted rotary kilns: experimental study of unloading at ambient and high temperatures. accepted for publication in[END_REF]. The first part of this curve is comparable to the results obtained by Ajayi and Sheehan (2012a) who investigated only relatively weak hold-up. Because the powders used are highly cohesive, the flights can be heavily overloaded. As a result, saturation for F is not yet achieved; F continues therefore to increase slowly, even for values of Z up to 50%.

Powder fraction falling through the gas

In the operating conditions of industrial kilns, the fraction G is always very small. As q is proportional to the rotational speed, G also varies proportionally to , as shown on Fig. 7b. Fig. 7a gives an example of the variations in the fraction G (in red/grey) as a function of the local kiln hold-up: at low hold-up, G increases with Z up to a maximum value identical to that for F (3.5% in the conditions of Fig. 7a); beyond this value, G decreases slightly since the length of fall decreases when the bed depth increases, due to the accumulation of solid in the bulk (B increases). The beginning of this curve is again comparable to the results obtained by Ajayi and Sheehan (2012a).

Axial profile of local kiln hold-up and transverse solid distribution

Fig. 8 shows the axial profile of the local kiln hold-up for kiln #1 under nominal operating conditions (see Table 1) and the average solid distribution at each axial position; molar flow rates of the reactant (UO2F2) and reaction product (U3O8) are also plotted in Fig. 8. The local kiln hold-up is very high in the entry zone; it decreases by a factor of five between the entry and exit of the conversion kiln, due to the volume contraction between UO2F2 at the entry and U3O8 at the exit (see bulk densities in Table 2). The powder fraction falling through the gas varies only slightly along the cylinder in comparison with the other fractions; it is approximately 0.3%. The total flighted powder fraction increases slightly from the rotary kiln entry to the 6 meter mark, then remains practically constant at around 4%; this corresponds to the gradual decrease in the fraction in the bulk B up to the 6 meter mark: as the detachment angle D decreases, the lifters have more space outside the bulk, so the flighted powder fraction F increases. Lastly, solid remains in the bulk almost up to the exit end of the conversion kiln. Fig. 9 shows the axial profile of local kiln hold-up and the average solid distribution along kiln #2 under nominal operating conditions (see Table 1); molar flow rates of the reactant (UO2F2), of the intermediate (U3O8) and of the reaction product (UO2) are also plotted in Fig. 9. Between the entry and exit of the conversion kiln, the local kiln hold-up decreases by a factor of three, again due to the volume contraction between UO2F2 at the entry and UO2 at the exit (see bulk densities in Table 2). This decrease takes place mainly in the first half of the kiln, i.e. in the zone where UO2F2 is consumed (compare its flow rate with hold-up in Fig. 9); as the bulk density difference between U3O8 and UO2 is smaller, the local kiln hold-up thus decreases slowly. The powder fraction falling through the gas is about 0.4% along the whole length of the kiln. The total flighted powder fraction is about 3.3%; it increases slightly from the entry of the conversion kiln to the 2.3 meter mark, simultaneously with the decrease of Z; then F decreases slightly until the kiln exit, parallel to the decrease of Z. Lastly, no solid remains outside the lifters -that is in the bulk-between the 2.3 meter mark and the exit, because Z is too low in this area: the powder falling from the lifters no longer supplies the bulk, but fills the flights directly. Furthermore the lifters cannot be completely filled and are less and less loaded as Z decreases, thus F continues to decline until the exit of the kiln.

The main difference between the two types of kiln is their hold-up: the flow rate in kiln #2 is about ten times lower than in kiln #1, while kiln #2 is only about half the size of kiln #1 (see dimensions and nominal operating conditions of the kilns in Table 1). This has an important impact on the axial profile of the local kiln hold-up (in particular on the hold-up value at the conversion kiln entry) and on the relative distribution of powder falling through the gas, flighted powder and bulk powder, as can be seen by comparing Figs. 8 and9. Industrial rotary kilns are not standalone devices: they are preceded by a hydrolysis reactor which feeds directly into the kiln through an endless screw, which prevents backflow by continuously pushing the powder. From the experience of the operators of kilns #1, problems with the screw periodically occur that are consistent with a high filling degree at this location; however the operators of kilns #2 never faced this type of problem, which is consistent with the substantially lower filling degree calculated at the entrance of this second type of kilns.

Figs. 8 and9 show that the volume contraction effect resulting from the reactions that take place in uranyl difluoride conversion kilns is very significant. In order to validate the model presented here, it would be necessary to measure the axial profiles of local kiln hold-up and the transverse solid distribution under reactive conditions, i.e. in the presence of uranium and hydrofluoric acid, at a temperature around 800°C. The few measurements reported in the literature [START_REF] Lebas | Experimental study of residence time, particle movement and bed depth profile in rotary kilns[END_REF][START_REF] Chatterjee | Flow of Materials in Rotary Kilns Used for Sponge Iron Manufacture: Part I. Effect of Some Operational Variables[END_REF][START_REF] Sai | Residence time distribution and material flow studies in a rotary kiln[END_REF][START_REF] Spurling | The transient response of granular flows in an inclined rotating cylinder[END_REF] were not only conducted on unflighted rotary kilns but more importantly, were completely out of the question here for reasons concerning the safety of the workers and of the installations. We had originally intended to exploit the weak radioactivity of uranium in order to measure the overall axial distribution of the powder along the kiln. Unfortunately this proved to be impossible due to the congested surroundings of the industrial kiln (lagging, various systems of measurement and monitoring, etc.): the gamma-ray signal was not usable.

CONCLUSIONS

In this study, a method to calculate the axial and transverse solid distribution of the powder in flighted rotary kilns is proposed (in the form of a series of 'snapshots', as shown on Fig. 0 / Graphical abstract), i.e. not only the quantity of powder falling through the gas, but also the total quantity of flighted powder, as well as the quantity remaining in the bulk. This method does not require defining a particle size, which would be meaningless in the case of fine-grained cohesive powders. The method is based on geometrical calculations and on experimental measurements conducted in the temperature conditions of the industrial plant being modeled.

The very substantial volume contraction in the modeled conversion kilns leads to a significant variation in the local hold-up between the entry (uranyl difluoride) and the exit (uranium dioxides) of the kilns. Whatever the type of conversion kiln considered, the powder fraction falling through the gas varies only slightly along the length of the kiln and is below 0.5%. The total flighted powder fraction varies between 2% and 4% along the kiln. The powder fraction remaining in the bulk varies considerably, up to 41% in kiln #1 and from 0 to 11% in kiln #2 for the nominal operating conditions of these conversion kilns. It should be emphasized that the results obtained here are specifically related to the highly cohesive character of the powders studied and to the fact that the conversion leads to a very strong variation in their properties (bulk density in particular).

The hydrodynamic model described in the present paper forms the basis for the global model of rotary kilns used for the conversion of uranium [START_REF] Debacq | Étude et modélisation des fours tournants de défluoration et réduction du difluorure d'uranyle[END_REF]. With the 'snapshots' that have been determined thanks to this hydrodynamic model, it is possible to calculate all the thermal exchange surfaces. The calculation of the flighted powder fraction, the powder fraction falling through the gas, and the fraction remaining in the bulk, enable the mass transfers to be assessed using the mass transfer model [START_REF] Debacq | Étude et modélisation des fours tournants de défluoration et réduction du difluorure d'uranyle[END_REF].

NOMENCLATURE

Latin characters dp [m]

particle size h [m] bed depth at a given axial position hfall [m] length of fall from a lifter ho [m] lifter dimension normal to the kiln wall hp [m] lifter dimension parallel to the kiln wall n [-] number of flights inside kiln q [-] powder fraction falling from one lifter between time t and t+t Qv [m 3 s -1 ] volumetric flow rate R [m] inside kiln radius tfall [s] time of fall from a lifter vfall [m s -1 ] fall velocity of the powder x [m] axial position in the kiln (distance from the solid feed) Z [-] local kiln hold-up (at a given axial position)

Greek symbols  [°] inclination of the kiln  [°-1 ] slope of the lifter discharge law  [°]
avalanche angle (instantaneous angle of repose made by particles in a flight with respect to horizontal)

 [°]
angular position of the lifter (with respect to the horizontal)

D [°] detachment angle F [°]
angle at the end of discharge I [°] angle at the start of effective discharge  [rad] filling angle (central angle extended by solid bed) at a given axial position B [-] fraction of the hold-up that remains in the bulk (in the bottom of the kiln) at a given axial position F [-] fraction of the hold-up that is flighted at a given axial position ^F  [-] total flighted powder fraction at a given angular and axial position G [-] fraction of the hold-up that is in gaseous phase (airborne) at a given axial position

       tan h h R BC h R OC p o o , and       cos R OD cos OC ,
The coordinates of the points in the reference frame (O,x,y) are as follows:

                cos h sin h R sin h cos h R A p o p o ;           sin h R cos h R B o o ;                 sin tan h h R cos tan h h R C p o p o                                   2 cos R 2 sin R 2 2 cos R G ;   F p o y cos tan h h R F      
F belongs to the straight line of slope tan  passing through G, hence the ordinate yF verifies:

                                       2 sin R cos tan h h R tan 2 cos R y p o F Lastly, as F A fall y y h  
, the length of fall in a straight cross-section of the rotary kiln with partial hold-up is: Fig. 9. Solid distribution in nominal operating conditions of kiln #2; molar flow rates obtained from the mass transfer model [START_REF] Debacq | Étude et modélisation des fours tournants de défluoration et réduction du difluorure d'uranyle[END_REF]. The densities are measured by the method of the test tube (250 mL; inside diameter 39 mm). They are obtained  20 kg m -3 . The particle size distributions of samples of oxides (U3O8 and UO2) were measured with a COULTER 1.31 1.8 in water. UO2F2 particle size distributions are measured in a COULTER 2.00 2.2 in chloroform (powder highly soluble in water). The values given in the table are the mean sizes corresponding to main peaks (bimodal distributions). The static angle of repose is measured using the crumbling cone method (NF T 20-221 standard): on a cone of powder formed using a calibrated funnel (92.9 mm larger diameter; 15.0 mm smaller diameter; 107 mm high); the funnel is in non-magnetic metal of surface roughness 0.16 µm. The dynamic angles of repose for the four powders were measured by the rotating cylinder technique, using drums without lifting devices; its diameter is the same as in the corresponding industrial kiln; the length is 0.3 m.. Sixteen experiments were performed, with hold-up ranging between 7.5% and 16 % and rotational speed from 2 to 10 rpm. Each experiment was filmed for about two minutes after reaching steady state and one image every five seconds was extracted for analysis. The dynamic angles of repose were measured from the images using specific software [START_REF] Debacq | Transverse motion of cohesive powders in flighted rotary kilns: experimental study of unloading at ambient and high temperatures. accepted for publication in[END_REF].

                                                2 sin R cos tan h h R tan 2 cos R cos h sin h R h p o p o fall

FIGURE (appendix)

FIGURESFig. 0 /

 0 FIGURES (main text)

Fig. 1 .

 1 Fig. 1. Geometrical symbols used.

Fig. 2 .

 2 Fig. 2. Falling curtains from lifter tips at different angular positions overlaid on the same image: (a) UO2F2 in drum #1,  = 3 rpm; (b) UO2 in drum #2,  = 5 rpm.

Fig

  Fig. 3.  versus       2 sin

Fig. 4 .

 4 Fig. 4. Steps 3 and 4 of the procedure used to calculate the powder fraction G falling through the gas at a given axial position (section 2.2.1).

Fig. 5 .

 5 Fig. 5. Instantaneous powder fraction falling through the gas versus time during 1 + 1/n revolution (UO2F2 in kiln #1,  = 2 rpm).

Fig. 6 .

 6 Fig. 6. Total flighted powder fraction versus angular position during one revolution (UO2F2 in kiln #1).

Fig. 7 .

 7 Fig. 7. (a) In green/black: flighted powder fraction versus hold-up; in red/grey: powder fraction falling through the gas versus hold-up (U3O8, dimensions of kiln #2,  = 7 rpm). (b) Powder fraction falling through the gas versus rotational speed (U3O8, dimensions of kiln #1, Z = 30,8%).

Fig. 8 .

 8 Fig. 8. Solid distribution in nominal operating conditions of kiln #1; molar flow rates obtained from the mass transfer model (Debacq, 2001).

Fig. A. 1 .

 1 Fig. A. 1. Symbols used to calculate the length of fall.

Table 1 .

 1 TABLES Nominal operating conditions of rotary kilns used. At the entry of the conversion process, UF6 is instantaneously and completely converted to UO2F2 in powder form; the "UO2F2 equivalent flow rate" is therefore calculated from the flow rate of UF6 and the respective molecular masses of UF6 and of UO2F2.

		kiln #1	kiln #2
	inside diameter	750 mm	348 mm
	length	10.24 m	5.34 m
	inclination	2.5%	2.5%
	number of lifters	6	4
	lifter dimension normal to the kiln wall	60 mm	35 mm
	lifter dimension parallel to the kiln wall	60 mm	35 mm
	UF6 flow rate	860 kg/h	80 kg/h
	UO2F2 equivalent flow rate*	753 kg/h	70 kg/h
	counter-current H2O flow rate	100 kg/h	18 kg/h
	counter-current H2 flow rate	18 Nm 3 /h	12 Nm 3 /h
	final product of reactions	U3O8	UO2
	rotational speed	2 rpm	7 rpm
	heat devices temperatures	750; 900; 900; 780 °C	700; 730; 740; 710; 650; 630 °C

*

Table 2 .

 2 Properties of the different powders.

	UO2F2	U3O8 type #1 U3O8 type #2 UO2

Franco-Belge de Fabrication de Combustible

Compagnie Générale des Matières Nucléaires

Acknowledgements

This work was supported by a grant from FBFC 2 Romans (France) and Cogema 3 Pierrelatte (France). We gratefully acknowledge André FEUGIER and Didier HARTMANN for their kind help in this work.


Walid Ben-Ameur
email: walid.benameur@telecom-sudparis.eu

Antoine Glorieux
email: antoine.glorieux@telecom-sudparis.eu

José Neto
email: jose.neto@telecom-sudparis.eu

A full description of polytopes related to the index of the lowest nonzero row of an assignment matrix

Consider a {0, 1} assignment matrix where each column contains exactly one coefficient equal to 1 and let h be the index of the lowest row that is not identically equal to the zero row. We give a full description of the convex hull of all feasible assignments appended with the extra parameter h. This polytope and some of its variants naturally appear in the context of several combinatorial optimization problems including frequency assignment, job scheduling, graph orientation, maximum clique, etc. We also show that the underlying separation problems are solvable in polynomial time and thus optimization over those polytopes can be done in polynomial time.

z 1 z 2 z 3 • • • z n k . . . h -→ ←-h 3 2 1
Fig. 1: On this grid representing a k-byn matrix, each column is a z i variable and the height of its black cell is the value assigned to it. h thus corresponds to the first non-empty row of the grid starting from the bottom (in this case, h = 4).

Introduction and motivations

Let us consider combinatorial optimization problems involving n ∈ N\{0} variables z i each of which can be assigned an integer number in 1, k and let h ≡ min n i=1 z i . A natural polytope related to these problems is given by

P = Conv                    y 1 k y 2 k • • • y n k . . . . . . . . . . . . y 1 2 y 2 2 • • • y n 2 y 1 1 y 2 1 • • • y n 1      , h      ∈ M k,n × N\{0} k ∑ l=1 y i l = 1, ∀i ∈ 1, n , h = min i∈ 1,n ∑ k l=1 ly i l ,         
, where M k,n is the set of all k-by-n matrices with coefficients in {0, 1} and the variables y i l are interpreted as follows : y i l = 1 if and only if z i = l. The matrix (y i j ) can be seen as a {0, 1} assignment matrix where each column contains exactly one coefficient equal to 1 while h denotes the index of the lowest row that is not identically equal to the zero row (cf Fig. 1).

Another variant of P is obtained by considering the index of the highest row that is not identically equal to the zero row. In this case we get the polytope

P = Conv                    x 1 k x 2 k • • • x n k . . . . . . . . . . . . x 1 2 x 2 2 • • • x n 2 x 1 1 x 2 1 • • • x n 1      , g      ∈ M k,n × N\{0} k ∑ l=1 x i l = 1, ∀i ∈ 1, n , g = max i∈ 1,n ∑ k l=1 lx i l ,         
.

Polytopes P and P naturally appear in the context of several combinatorial optimization. Let us for example consider the minimum-span frequency-assignment problem which is a variant of the NP-hard frequency-assignment problem [START_REF] Mann | Complexity of different ILP models of the frequency assignment problem[END_REF]. Given a simple graph G = (V, E) that is generally called the interference graph, the frequency assignment problem consists in assigning a frequency f from a set of available frequencies F to each vertex v ∈ V in such a way that each pair of antennas uv ∈ E that may interfere with one another are assigned different frequencies. Frequencies can be seen as ordered integer numbers. To reduce interferences, one might impose stronger constraints: a minimum separation between the frequencies assigned to u and v is required. If frequency i is assigned to u and j is assigned to v, then |i -j| ≥ s uv where s uv is a given number. The minimum-span frequency-assignment problem (or MS-FAP) consists in assigning frequencies to nodes taking into account the separation requirements and minimizing the difference between the largest assigned number (frequency) and the smallest assigned number (see, e.g., [START_REF] Koster | Frequency Assignment -Models and Algorihtms[END_REF]).

If we consider that

V = {v 1 , • • • , v n }, F = 1, k
where k is an upper bound of the minimum-span, then we obtain the following formulation for MS-FAP

     min g s.t. x i l + x j l ≤ 1, ∀(i, j, l, l ) ∈ 1, n 2 × 1, k 2 such that v i v j ∈ E, |l -l | < s v i v j (x, g) ∈ P , x ∈ M k,n .
where the interpretation of the x variable is the following: x i l = 1 if and only if the frequency l is assigned to the antenna v i .

Another example is the minimum makespan scheduling, which is a central problem in the scheduling area (see [START_REF] Vazirani | Minimum Makespan Scheduling' Approximations Algorithms[END_REF]). Given a set J of jobs, a set M of machines that can all process at most one job at a time, and the time t i, j ∈ N taken to process job j ∈ J on machine i ∈ M , the goal of the minimum makespan scheduling problem is to assign a machine p ∈ M for each job j ∈ J so as to minimize the makespan, i.e. the maximum processing time of any machine. Several approximation schemes have been developed to deal with this NP-hard problem [START_REF] Garey | Computers and Intractability; A Guide to the Theory of NP-Completeness[END_REF], e.g. [START_REF] Hochbaum | A polynomial approximation scheme for scheduling on uniform processors: Using the dual approximation approach[END_REF] and [START_REF] Horowitz | Exact and approximate algorithms for scheduling nonidentical processors[END_REF]. Since the processing times are integers, the timeline is discretized in identical units of time, e.g. days. We consider here the variant where all the machines in M are identical (or IM-MMS) and preemptions are not allowed. In other words, for any job j ∈ J, t i, j = t j , ∀i ∈ M. In this case, assigning a machine to each job is equivalent to assigning a day d to be the last day of processing this job, which also determines the first day d of processing and will therefore be processed by a machine free during the period [d, d ]. Now to make a formulation for IM-MMS with the set of jobs J = 1, n and m ∈ N\{0} identical machines, we take k = ∑ n i=1 t i and the variable x ∈ M k,n whose interpretation is the following: x i l = 1 if and only if the processing of the job i ends on the day l. Then we have the following formulation for IM-MMS

           min g s.t. ∑ k l=1 lx i l ≥ t i , ∀i ∈ 1, n , ∑ n i=1 ∑ min(l+t i -1,k) l =l x i l ≤ m, ∀l ∈ 1, k , (x, g) ∈ P , x ∈ M k,n .
For a job i ∈ 1, n , ∑ k l=1 lx i l ≥ t i ensures that its processing ends after enough time has passed for i to be processed, and for a day

l ∈ 1, k , ∑ n i=1 ∑ min(l+t i -1,k) l =l
x i l ≤ m ensures that no more than m jobs are being processed. Some additional constraints can be added to this formulation such as a necessary precedence or release time. If we want a job i ∈ 1, n to be processed before another job j ∈ 1, n \{i} starts processing, we add the constraint ∑ k l=1 lx i l ≤ ∑ k l=1 lx j lt j . If we want a job i ∈ 1, n to be processed before (resp. on, after) a day d ∈ 1, k , we add the constraint ∑ d-1 l=1 x i l = 1 (resp.

x i d = 1, ∑ k l=d+1 x i l = 1
). The objective function can also be any linear function depending on the x and g variables.

Another example is given by the problem of the most imbalanced orientation of a graph (or MAXIM) that consists in orienting the edges of a graph such that the minimum over all the vertices of the absolute difference between the outdegree and the indegree of a vertex is maximized (NP-complete) [START_REF] Ben-Ameur | On the most imbalanced orientation of a graph[END_REF]. In other words, for a simple graph

G = (V, E), MAXIM(G) = max Λ ∈ - → O (G) min v∈V |d + Λ (v) -d - Λ (v)|, where - → O (G) denotes the set of all the orientations of G and d + Λ (v) (resp. d - Λ (v)
) denotes the outdegree (resp. indegree) of v in G with respect to Λ . Now if we consider the graph G to be arbitrarily oriented and take its incidence matrix B ∈ {-1, 0, 1} |V |×|E| , we can describe an orientation of G with the variable x ∈ {-1, 1} |E| interpreted as follows. For each edge uv ∈ E, its orientation is kept from the original one if x uv = 1 and reversed otherwise. Then if we look at the product of B with an orientation vector

x ∈ {-1, 1} |E| we obtain B v x = d + x (v) -d - x (v), ∀v ∈ V
, where B v denotes the row of B corresponding to node v. In order to make a formulation of MAXIM, we consider indicator variables t v l ∈ {0, 1} with v ∈ V and l ∈ -k, k , k being the maximum degree of the vertices of G, that have the following interpretation:

t v l = 1 if and only if B v x = d + x (v) -d - x (v) = l
, and thus we obtain the formulation

                 max h s.t. h ≤ ∑ k l=-k |l|t v l , ∀v ∈ V, ∑ k l=-k t v l = 1, ∀v ∈ V, ∑ k l=-k lt v l = B v x, ∀v ∈ V, x ∈ [-1; 1] |E| , t ∈ M n,2k+1 , h ∈ R. Introducing variables y v l = t v -l + t v l , ∀(v, l) ∈ V × 1, k , it becomes            max h s.t. ∑ k l=-k lt v l = B v x, ∀v ∈ V, y v l = t v -l + t v l , ∀(v, l) ∈ V × 1, k , x ∈ [-1; 1] |E| , (y, h) ∈ P, t ∈ M n,2k+1 .
Considering the last formulation, a polyhedral analysis of the polytope P may be helpful in strengthening the linear relaxation of the original formulation within the framework of a cutting-plane algorithm (see, e.g. [START_REF] Cook | Combinatorial Optimization[END_REF] and [START_REF] Wolsey | Integer Programming[END_REF]).

Polytope P also appears in the context of the maximum clique problem. A discretized formulation is proposed in [START_REF] Martins | Extended and discretized formulations for the maximum clique problem[END_REF] where a variable w i q indicates whether the vertex i belongs to a clique of size q. These variables are of course linked to standard vertex variables (x i = 1 if i belongs to the maximum clique). The problem is then equivalent to maximizing q such that w i q = 1 for some i. This is again related to polytope P. More generally, several combinatorial optimization problems where discretization techniques are used can benefit from a description of either P or some of its variants.

The rest of the paper is organized as follows. First we give a complete linear description of P in Section 2. Then we show that the separation problem with respect to P can be solved in polynomial time in Section 3. Finally we give similar results for the polyhedron P and others in Section 4 and conclude in Section 5.

A full description of P

Let us define a set of inequalities that will prove to be an hyperplane representation of P.

P =        ∑ k l=1 y i l = 1, ∀i ∈ 1, n , ∑ k l=2 ∑ n i=1 λ i l y i l ≥ h -1, ∀λ ∈ Λ , ∑ h max -1 l=1 ∑ n i=1 (l -h max )y i l + h max ≤ h, ∀h max ∈ 1, k , y i l ≥ 0, ∀(i, l) ∈ 1, n × 1, k , h ∈ R,
where

Λ = λ = (λ i l ) (i,l)∈ 1,n × 1,k ∈ N nk λ i l+1 ≥ λ i l , ∀(i, l) ∈ 1, n × 1, k -1 ∑ n i=1 λ i l = l -1, ∀l ∈ 1, k .
Any element λ of Λ can be constructed as follows: we start with λ i 1 = 0, ∀i ∈ 1, n , choose an index i 2 ∈ 1, n and set λ i 2 2 = 1 and λ i 2 = 0, ∀i ∈ 1, n \{i 2 }. And we proceed like this for l = 2, • • • , k, we choose an index i l ∈ 1, n and set

λ i l l = λ i l l-1 + 1 and λ i l = λ i l-1 , ∀i ∈ 1, n \{i l }. Lemma 1 P ⊆ P
Proof. Since P is the convex hull of integer points, it suffices to show that each of those points satisfies all the inequalities in P. Let (y, h) ∈ P be one of those points, that is to say

(y, h) ∈ M k,n ×N\{0}, ∑ k l=1 y i l = 1, ∀i ∈ 1, n and h = min i∈ 1,n ∑ k l=1 ly i l . We firstly show that ∑ k l=2 ∑ n i=1 λ i l y i l ≥ h -1, ∀λ ∈ Λ . For all i ∈ 1, n , there exists l i ∈ h, k such that y i l i = 1. Now since for a given i ∈ 1, n , λ i l increases with l, we have k ∑ l=2 n ∑ i=1 λ i l y i l = n ∑ i=1 λ i l i ≥ n ∑ i=1 λ i h = h -1.
Now we take h max ∈ 1, k , and we show that h max -∑

h max -1 l=1 (h max -l) ∑ n i=1 y i l ≤ h. We have h max - h max -1 ∑ l=1 (h max -l) n ∑ i=1 y i l = h max + k ∑ l=1 min(l -h max , 0) n ∑ i=1 y i l = h max + k ∑ l=1 min(l -1, h max -1) n ∑ i=1 y i l -n(h max -1) = n ∑ i=1 k ∑ l=2 min(l -1, h max -1)y i l + 1 -(n -1)(h max -1)
There exists i * ∈ 1, n such that

y i * h = 1, then n ∑ i=1 k ∑ l=2 min(l -1, h max -1)y i l = n ∑ i=1 i =i * k ∑ l=2 min(l -1, h max -1)y i l + k ∑ l=2 min(l -1, h max -1)y i * l ≤ n ∑ i=1 i =i * k ∑ l=2 (h max -1)y i l + k ∑ l=2 (l -1)y i * l ≤ n ∑ i=1 i =i * (h max -1) + h -1 = (n -1)(h max -1) + h -1.
Now to prove that P coincides with P, we show that all facet-defining inequalities for P are among those defining P. Two inequalities are said to be equivalent if one can be obtained from the other by multiplying it by a non-zero scalar and adding a combination of equations of the type ∑ k l=1 y i l = 1.

Lemma 2 Let k ∑ l=1 n ∑ i=1 β i l y i l + γ ≥ 0, (1)
be a facet-defining inequality of P, with

β i l ∈ R, ∀(i, l) ∈ 1, n × 1, k , γ ∈ R. Then there exists (i, l) ∈ 1, n × 1, k such that (1) is equivalent to y i l ≥ 0.
For an extreme point (y, h) of P and ( ĩ, l, l ) ∈ 1, n × 1, k 2 , such that y ĩ l = 1 and l = l , we denote by (y l ĩ -→ l , h l ĩ -→ l ) the extreme point (y , h ) of P such that

y i l = y i l , ∀(i, l) ∈ 1, n × 1, k \{( ĩ, l), ( ĩ, l )}, y ĩ l = 0, y ĩ l = 1 and h = min i∈ 1,n ∑ k l=1 ly i l . For l ∈ 1, k , we denote by (y → l , l) the point of P such that y i l = 1, ∀i ∈ 1, n and y i l = 0, ∀(i, l) ∈ 1, n × ( 1, k \{ l}).
Proof. First, since for each i ∈ 1, n , we have ∑ k l=1 y i l = 1, we can replace y i 1 by 1 -∑ k l=2 y i l and get new coefficients β i 1 = 0 and

β i l = β i l -β i 1 , ∀l ∈ 2, k and a new γ = γ + ∑ n i=1 β i 1 .
Hence, without loss of generality, we can assume that β i 1 = 0 for all i ∈ 1, n . Suppose that the facet defined by ( 1) is not equivalent to a facet defined by an inequality of the type y i l ≥ 0. If we take ( ĩ, l) ∈ 1, n × 1, k -1 , we know that there exists an extreme point (y, h) of P saturating (1) and such that y ĩ l = 1, otherwise all the extreme points saturating (1) would saturate y ĩ l ≥ 0 thus contradicting the fact that (1) is facetdefining and not equivalent to y i l ≥ 0 for some

(i, l) ∈ 1, n × 1, k . Since (y , h ) = (y l ĩ -→ l+1 , h l ĩ -→ l+1 ) ∈ P, we have ∑ k l=1 ∑ n i=1 β i l y i l + γ = 0 and ∑ k l=1 ∑ n i=1 β i l y i l + γ ≥ 0 which yields β ĩ l+1 ≥ β ĩ l .
Similarly, taking an extreme point (y, h) of P saturating (1) and such that y ĩ l+1 = 1 and (y , h ) = (y l+1 ĩ -→ l , h l+1 ĩ -→ l ) ∈ P, we have ∑ k l=1 ∑ n i=1 β i l y i l + γ = 0 and ∑ k l=1 ∑ n i=1 β i l y i l + γ ≥ 0, yielding β ĩ l ≥ β ĩ l+1 . Hence, for all (i, l) ∈ 1, n × 1, k -1 , we have β i l = β i l+1 , in other words, β i l = 0, ∀(i, l) ∈ 1, n × 1, k , and (1) is not facetdefining.

Lemma 3 Let k ∑ l=1 n ∑ i=1 β i l y i l + γ ≥ h, (2)
be a facet-defining inequality of P, with

β i l ∈ R, ∀(i, l) ∈ 1, n × 1, k , γ ∈ R. Then there exists λ ∈ Λ such that (2) is equivalent to ∑ k l=2 ∑ n i=1 λ i l y i l ≥ h -1.
Proof. Again, without loss of generality, we assume that β i 1 = 0 for all i ∈ 1, n . For an ( ĩ, l) ∈ 1, n × 1, k -1 , there exists an extreme point (y, h) of P saturating (2) and such that

y ĩ l = 1. Since (y , h ) = (y l ĩ -→ l+1 , h l ĩ -→ l+1 ) ∈ P, we have ∑ k l=1 ∑ n i=1 β i l y i l + γ = h and ∑ k l=1 ∑ n i=1 β i l y i l + γ ≥ h ≥ h which yields β ĩ l+1 ≥ β ĩ l .
Hence for all i ∈ 1, n , (β i l ) l is increasing with l and therefore non-negative since β i 1 = 0, ∀i ∈ 1, n . If we consider the point (y, h) = (y →1 , 1), we obtain that γ ≥ 1. If we now consider an extreme point (y, h) of P saturating (2) and such that y ĩ 1 = 1 for an ĩ ∈ 1, n , we get ∑ k l=1 ∑ n i=1 β i l y i l + γ = h = 1. But since both the β i l and the y i l are non-negative, then so is

∑ k l=1 ∑ n i=1 β i l y i l . Hence γ ≤ 1, yielding γ = 1.
Considering (y →l , l) ∈ P for l ∈ 1, k , we obtain ∑ n i=1 β i l ≥ l -1. Let us show by induction on l that ∑ n i=1 β i l = l -1, ∀l ∈ 1, k . Our induction is already initialized by ∑ n i=1 β i 1 = 0. We suppose that for a l ∈ 1, k -1 , we have ∑ n i=1 β i l = l -1 and show that ∑ n i=1 β i l+1 = l. Suppose that all the extreme points (y, h) of P saturating (2) and such that y i l+1 = 1 for some i ∈ 1, n verify h ≤ l. Then for each ĩ ∈ 1, n , take one of those extreme saturating points such that y ĩ l+1 = 1 and let (y

, h ) = (y l+1 ĩ -→ l , h l+1 ĩ -→ l ) ∈ P. Since β ĩ l+1 ≥ β ĩ l , we have h -1 ≤ ∑ k l=1 ∑ n i=1 β i l y i l ≤ ∑ k l=1 ∑ n i=1 β i l y i l = h -1, and since h ≤ l, we have h = h and therefore, ∑ k l=1 ∑ n i=1 β i l y i l = ∑ k l=1 ∑ n i=1 β i l y i l = h -1, yielding β ĩ l+1 = β ĩ l . Thus ∑ n i=1 β i l+1 = ∑ n i=1 β i l = l -1 which contradicts ∑ n i=1 β i l+1 ≥ l.
So there exists ĩ ∈ 1, n and (y, h) an extreme point of P saturating (2) such that y ĩ l+1 = 1 and

h = l + 1. We have l ≤ ∑ n i=1 β i l+1 ≤ ∑ k l=1 ∑ n i=1 β i l y i l = l, hence ∑ n i=1 β i l+1 = l, which con- cludes our induction.

Now let us show by induction on l that for all

(i, l) ∈ 1, n × 1, k , β i
l is an integer. This induction is trivially initialized for β i 1 = 0, ∀i ∈ 1, n . We suppose that for a l ∈ 1, k -1 we have that the β i l for i ∈ 1, n are integers and we show that the same holds for the β i l+1 for i ∈ 1, n . We note α i = β i l+1β i l , ∀i ∈ 1, n and for each i ∈ 1, n we build a new set of inequality coefficients : β (i), j l = β j l -α j +δ i, j , ∀( j, l) ∈ 1, n × 1, k , where δ i, j equals 1 if i = j and 0 otherwise. Let (y, h) be an extreme point of P and for all j ∈ 1, n , let l j ∈ h, k be such that y j l j = 1. Then, since ∑ n i=1 α i = 1, we have for

i ∈ 1, n k ∑ l=1 n ∑ j=1 β (i), j l y j l = k ∑ l=1 n ∑ j=1 (β j l -α j +δ i, j )y j l = n ∑ j=1 (β j l j -α j +δ i, j ) = n ∑ j=1 β j l j ≥ n ∑ i=1 β i h = h-1, which means that for all i ∈ 1, n , ∑ k l=1 ∑ n j=1 β (i), j l y j l + 1 ≥ h is valid for P. Now since for ( j, l) ∈ 1, n × 1, k , n ∑ i=1 α i β (i) j l = n ∑ i=1 α i β j l -α j + δ i, j = n ∑ i=1 α i β j l - n ∑ i=1 α i α j + n ∑ i=1 α i δ i, j = β j l and α i ≥ 0, ∀i ∈ 1, n , (2)
is a convex combination of these inequalities. Moreover, if any of the β i l+1 , i ∈ 1, n was not an integer, then the convex combination would be non-trivial, which would contradict the fact that (2) is facet-defining. Concluding our induction, we obtain that for all (i, l) ∈ 1, n × 1, k , β i l is an integer. And thus that β ∈ Λ , i.e. (2) belongs to the set of inequalities defining P.

P is defined by n equalities, kn non-negativity constraints, k constraints of type ∑

h max -1 l=1 ∑ n i=1 (l - h max )y i l + h max ≤ h and n k-1 of inequalities of type ∑ k l=2 ∑ n i=1 λ i l y i l ≥ h -1.
The total number of inequalities is then exponential. However, the following holds.

Theorem 6

The separation problem which consists in deciding if a vector (y, h) ∈ R nk+1 is in P, and if not in returning a constraint of P violated by (y, h) can be solved in polynomial time.

Proof. Let (y, h) ∈ R nk+1 , first, one can verify in linear time if (y, h) ∈ [0, 1] nk × [1, k] is such that ∑ k l=1 y i l = 1, ∀i ∈ 1, n and verifies the k inequalities of type ∑ h max -1 l=1 ∑ n i=1 (l - h max )y i l + h max ≤ h. If not,
λ i l y i l ≥ h -1 and let λ ∈ Λ and (i 2 , • • • , i k ) ∈ 1, n k-1 the indices corresponding to the building of λ , i.e. for l = 2, • • • , k, λ i l l = λ i l l-1 + 1 and λ i l = λ i l-1 , ∀i ∈ 1, n \{i l }. By construction of i 2 , • • • , i k and ĩ2 , • • • , ĩk , we have k ∑ l=2 n ∑ i=1 λ i l y i l = k ∑ l=2 (y i l l +y i l l+1 +• • •+y i l k ) ≥ k ∑ l=2 (y ĩl l +y ĩl l+1 +• • •+y ĩl k ) = k ∑ l=2 n ∑ i=1 λ i l y i l ≥ h-1,
hence the inequality of P corresponding to λ is satisfied. So we can conclude that if (y, h) satisfies the inequality of P corresponding to λ , then it satisfies all the inequalities of P corresponding to an element of Λ . And since the construction of λ is done in polynomial time, the separation problem is indeed polynomial.

The previous result is very useful in the context of cutting plane algorithms where only violated inequalities are added and not all valid inequalities.

Variants

Linear programming formulations aiming to maximize (resp. minimize) the index of the lowest (resp. highest) nonzero row of an assignment matrix are related to polytope Q (resp. Q ) described below. Observe that h (resp. g) is only required to be less (resp. more) than or equal to min i∈ 1,n ∑ k l=1 ly i l (resp. max i∈ 1,n ∑ k l=1 lx i l ).

Q = Conv                    y 1 k y 2 k • • • y n k . . . . . . . . . . . . y 1 2 y 2 2 • • • y n 2 y 1 1 y 2 1 • • • y n 1      , h      ∈ M k,n × N\{0} ∑ k l=1 y i l = 1, ∀i ∈ 1, n , h ≤ min i∈ 1,n ∑ k l=1 ly i l ,          Q = Conv                    x 1 k x 2 k • • • x n k . . . . . . . . . . . . x 1 2 x 2 2 • • • x n 2 x 1 1 x 2 1 • • • x n 1      , g      ∈ M k,n × N\{0} ∑ k l=1 x i l = 1, ∀i ∈ 1, n , g ≥ max i∈ 1,n ∑ k l=1 lx i l ,          A full description of Q is given below. Theorem 7 Q =    ∑ k l=1 y i l = 1, ∀i ∈ 1, n , ∑ k l=2 ∑ n i=1 λ i l y i l ≥ h -1, ∀λ ∈ Λ , y i l ≥ 0, ∀(i, l) ∈ 1, n × 1, k , h ≥ 1.
Proof. It is a simple fact that h ≥ 1 is the only possible facet of type (3) while the positivity constraints are the only possible facets of type [START_REF] Ben-Ameur | On the most imbalanced orientation of a graph[END_REF]. Let us consider an inequality of type (2) defining a facet of Q. Any extreme point of Q saturating such a facet necessarily satisfies h = min i∈ 1,n ∑ k l=1 ly i l implying that it is also a point of P. Using this observation and the fact that Q and P have the same dimension, we deduce that any facet of Q of type ( 2) is also a facet of P. Using the description of P, we get the result.

Similarly to Theorem 6 we can deduce that the separation problem with respect to Q is solvable in polynomial time as well.

We can also derive a full description of P and Q from the previous results.

Theorem 8 P =        ∑ k l=1 x i l = 1, ∀i ∈ 1, n , ∑ k-1 l=1 ∑ n i=1 λ i l x i l ≤ g -k, ∀λ ∈ Λ , ∑ k l=g min +1 ∑ n i=1 (l -g min )x i l + g min ≥ g, ∀g min ∈ 1, k , x i l ≥ 0, ∀(i, l) ∈ 1, n × 1, k , g ∈ R, Q =    ∑ k l=1 x i l = 1, ∀i ∈ 1, n , ∑ k-1 l=1 ∑ n i=1 λ i l x i l ≤ g -k, ∀λ ∈ Λ , x i l ≥ 0, ∀(i, l) ∈ 1, n × 1, k , g ≤ k where Λ = λ = (λ i l ) (i,l)∈ 1,n × 1,k ∈ N nk λ i l+1 ≤ λ i l , ∀(i, l) ∈ 1, n × 1, k -1 ∑ n i=1 λ i l = k -l, ∀l ∈ 1, k .
Proof. Take an extreme point (y, h) of P, and let (x, g) ∈ M k,n × N\{0} such that x i l = y i k-l+1 , ∀(i, l) ∈ 1, n × 1, k and g = kh + 1, then g = max i∈ 1,n ∑ k l=1 lx i l , hence (x, g) ∈ P . Conversely, any extreme point (x, g) of P can be obtained from an extreme point of P in this manner. So P is obtained from P doing the change of variables x i l = y i k-l+1 , ∀(i, l) ∈ 1, n × 1, k and g = k -h+1. Therefore its hyperplane representation is obtained from that of P in Theorem 5. Similarly, Q is obtained from Q doing the same change of variables and its hyperplane representation is thus obtained from Theorem 7.

The previous results imply that the separation problems related to P and Q can be solved in polynomial time.

Conclusion

In this paper we exhibited a family of polyhedra emerging in very diverse combinatorial optimization problems including the most imbalanced orientation of a graph, the minimum span frequency assignment and some scheduling problems. Then a full description of these polyhedra has been derived. We also proved that the separation problems related to these polyhedra can be solved in polynomial time and thus optimization over them can be done in polynomial time. We think that many combinatorial optimization problems where discretization techniques are used can benefit from the description of the polyhedra introduced in this paper. We are currently carrying out experimentations to study the efficiency of cutting plane algorithms based on these polyhedra. Future work may be directed towards investigations on extensions of the polyhedra we considered here in order to get better approximations while still keeping the feature of computational tractability. One can, for example, study {0, 1} assignment matrices appended with both h and g (the index of the lowest (resp. highest) nonzero row of the matrix). The related polytope is included in the intersection of P or Q. Some preliminary investigations show that more inequalities are necessary to describe the polytope.

  we return a violated constraint. Otherwise, we build λ ∈ Λ as follows: λ i 1 = 0, ∀i ∈ 1, n , and for l = 2, • • • , k, let ĩl = arg min i∈ 1,n y i , ∀i ∈ 1, n \{ ĩl }. We will show that if (y, h) satisfies the inequality of P corresponding to λ , then it satisfies all the inequalities corresponding to an element of Λ . Suppose ∑ k l=2 ∑ n i=1

			l +y i l+1 +• • •+y i k
	and set	λ ĩl l =	λ ĩl l-1 +1 and λ i l = λ i l-1

be a facet-defining inequality of P, with

(h maxl) ∑ n i=1 y i l ≤ h for some h max ∈ 1, k .

Proof. Since for each i ∈ 1, n , we have ∑ k l=1 y i l = 1, we can replace β i l by β i l = β i lv i , for some v i ≥ 0 such that β i l ≤ 0, ∀(i, l) ∈ 1, n × 1, k , and γ by γ = γ + ∑ n i=1 v i and thus get new coefficients β which are non-positive without changing [START_REF] Garey | Computers and Intractability; A Guide to the Theory of NP-Completeness[END_REF]. So without loss of generality, we can assume that β is non-positive and furthermore that γ is minimal for a non-positive β .

For ( ĩ, l) ∈ 1, n × 2, k , we take an extreme point (y, h) of P saturating (3) such that

In other words, β i l is increasing with l, for all i ∈ 1, n . This implies that for all i ∈ 1, n , there exists l i ∈ 1, k such that β i l = 0, ∀l ∈ l i , k and β i l < 0 for l < l i because suppose there exists i ∈ 1, n for which β i k < 0, then we can replace β i l by β i lβ i k for all l ∈ 1, k and add β i k to γ and thus get new non-positive coefficients β with a γ < γ, which contradicts the minimality of γ.

, l < l ĩ), we know there exists an extreme point (y, h) of P saturating (3) such that y ĩ l = 1. Suppose that h < l, we take (y , h ) an extreme point of P such that

which is a contradiction. Hence h = l, so if we take the extreme point (y , h ) of P such that

This gives us that for all i ∈ 1, n ,

Consider the extreme point (y, h) of P such that for all i ∈ 1, n , y i l i = 1, we have γ ≤ min i∈ 1,n l i . We call h max the maximum value of h among the extreme points (y, h) of P saturating (3).

If we take an extreme point (y, h) of P realizing h max , i.e. saturating (3) and such that h = h max , we have

l is nonpositive, yields h max ≤ γ. Now for ĩ ∈ 1, n , there exists an extreme point (y, h) of P saturating (3) such that y ĩ l ĩ-1 = 1. Suppose that h < l ĩ -1, we take (y , h ) an extreme point of P such that y

that yields β ĩ l ĩ-1 = 0, which is a contradiction. So h = l ĩ -1, hence h max ≥ l ĩ -1. We obtain max i∈ 1,n l i -1 ≤ h max ≤ γ ≤ min i∈ 1,n l i ≤ max i∈ 1,n l i . There are two possibilities, either min i∈ 1,n l i = max i∈ 1,n l i , or max i∈ 1,n l i -1 = min i∈ 1,n l i . Suppose max i∈ 1,n l i -1 = h max = γ = min i∈ 1,n l i , then there exists ĩ ∈ 1, n such that l ĩ = 1 + h max , which implies that β ĩ h max = 0 and β ĩ h max = h maxγ = 0 which is a contradiction. So min i∈ 1,n l i = max i∈ 1,n l i =: L and

Let us assume that h max = L -1. Then we know that γ < L, otherwise if γ = L, the extreme point of P (y →L , L) would saturate (3) and thus contradict the maximality of h max . We consider the following inequality

Let us show that it is a valid inequality for P, that is to say, that every extreme point (y, h) of P verifies it. If h ≥ L, then ∑ L-1 l=1 ∑ n i=1 (l -L)y i l = 0 and we are done. If h ≤ L -1, then there exists ĩ ∈ 1, n such that y ĩ h = 1. Combining this with the validity of (3) implies

Moreover, if (y, h) is an extreme point of P saturating (3), then there exists ĩ ∈ 1, n such that y ĩ h = 1 and h ≤ h max = L -1 which yields

ĩ (lγ)y i l = 0 which implies that

Consequently, all points saturating (3) saturate (4), furthermore, (y →L , L) saturates (4) and not (3). This means that the face of the polyhedron defined by ( 3) is strictly contained in the face defined by (4) contradicting the fact that (3) is facet-defining, hence h max = γ = L and (3) becomes ∑ h max -1 l=1

∑ n i=1 (lh max )y i l + h max ≤ h with h max ∈ 1, n .

Theorem 5

P = P

Proof. With Lemma 1, we know that P ⊆ P. Take any facet-defining inequality of P ∑ k l=1 ∑ n i=1 β i l y i l + γ ≥ αh, where β i l ∈ R, ∀(i, l) ∈ 1, n × 1, k , (γ, α) ∈ R × {-1, 0, 1}. If α = 0 (resp. α = 1, α = -1), then Lemma 2 (resp. Lemma 3, Lemma 4) gives us that this inequality is equivalent to one of those defining P.