Stephane Brisset
email: stephane.brisset@centralelille.fr

Maxime Ogier
email: maxime.ogier@centralelille.fr

COLLABORATIVE AND MULTILEVEL OPTIMIZATIONS OF A HYBRID RAILWAY POWER SUBSTATION

Keywords: Hybrid railway power substation, multi-source systems, energy storage, collaborative optimization, Benders decomposition

Three approaches are proposed for the design of a hybrid railway power substation taking into account the control of the storage system over one year. The first one is based on a piecewise linearization of the cost function and solved with linear programming. The others decompose the whole problem in two levels with lower dimensionality. Collaborative optimization solves it with a double-loop scheme while Benders decomposition allows a sequential resolution. This last gives the same solution as linear programming in about 5 times more computing time and can be used with a non-linear cost function.

INTRODUCTION

Because of the rail network development, power demand of railway power substation (RPS) becomes increasingly important [START_REF] Aeberhard | Railway traction power supply from the state of the art to future trends[END_REF]. The integration of renewable energies could compensate this need and it reduces pollution. Nevertheless the known intermittency factor of these energies makes it difficult to estimate the amount of produced power [START_REF] Clastres | Ancillary services and optimal household energy management with photovoltaic production[END_REF]. By introducing a storage device, the hybrid RPS (HRPS) could ensure a good power quality and it reduces total cost [START_REF] Aeberhard | Railway traction power supply from the state of the art to future trends[END_REF]. The storage device should be accurately used according to consumed power, produced power and evolution of prices fixed by the power supplier [START_REF] Rigo-Mariani | Optimal sizing integrating power management for a microgrid with storage[END_REF]. To enhance the design of HRPS, energy management strategy should be suitable in the long term. This requirement leads to a very high number of decision variables that is 52560 for a time span of one year and measurement every ten minutes.

This work falls within an emerging and challenging topic in electric engineering where the main issue is to design complex systems and networks for which the dimensional and control parameters have to be optimized all at once because of their strong interaction [START_REF] Clastres | Ancillary services and optimal household energy management with photovoltaic production[END_REF]- [START_REF] Aubry | Sizing optimization methodology of a surface permanent magnet machine-converter system over a torque-speed operating profile : Application to wave energy converter[END_REF]. Renewable energies like photovoltaic panels [START_REF] Clastres | Ancillary services and optimal household energy management with photovoltaic production[END_REF], wind turbines [START_REF] Rigo-Mariani | Optimal sizing integrating power management for a microgrid with storage[END_REF], and wave energy converters [START_REF] Aubry | Sizing optimization methodology of a surface permanent magnet machine-converter system over a torque-speed operating profile : Application to wave energy converter[END_REF] are typically intermittent with variations during the day and the seasons. The use of storage systems can balance climatic irregularity but their control along a large operating profile, that is typically one year or more, introduces numerous parameters leading to large-scale optimization problems.

Three approaches could be used for the design by optimization of HRPS. The first one is to linearize the cost function and solve the energy formulation problem by using linear programming (LP) in order to find an approximate solution. If the problem scale is high and the functions are non-linear then decomposition methods can be applied to divide the complete problem into sub-problems with lower dimensionality. In this paper, collaborative optimization (CO) and Benders decomposition (BD) are investigated to find an exact solution of the whole problem. However, the difficulty is in the coordination of sub-optimizations and the convergence to the solution of the whole problem.

In the first and second parts, the model and the optimization problem of HRPS are presented. The third part details the assumptions made to linearize the problem and compares the memory required by power and energy formulations with linear programming. Collaborative and multilevel optimizations are explained and applied to the HRPS problem in parts four and five, respectively. The three approaches are compared in terms of precision and computing time in part six. Finally, some conclusions and prospects are given.

MODEL OF HYBRID RAILWAY POWER SUBSTATION

The studied HRPS is presented in Fig. 1, it is composed of photovoltaic panels, wind turbines and a storage system connected to the substation through an AC bus. Renewable energy productions are illustrated by one-way power flows on Fig. 1. Their powers vary over time and depend on the photovoltaic panels surface 𝑆 PV , and the swept surface of wind turbines 𝑆 W :

(1)

𝑃 PV (𝑡) = 𝑆 PV 𝑃 irr (𝑡) ∀𝑡 ∈ [0, 𝑇] (2) 𝑃 W (𝑡) = 𝑘 p 𝑆 W 𝑣 adj 3 (𝑡) with 𝑣 adj (𝑡) = { 𝑣 W (𝑡) 0 ≤ 𝑣 W (𝑡) < 𝑣 rated 𝑣 rated 𝑣 rated ≤ 𝑣 W (𝑡) ≤ 𝑣 off 0 𝑣 off < 𝑣 W (𝑡) ∀𝑡 ∈ [0, 𝑇]
where 𝑇 is the end of time range, 𝑣 W (𝑡) is the wind speed in m/s, 𝑃 irr (𝑡) is the solar irradiance in W/m 2 , 𝑣 off is the turbine cut-off wind speed in m/s, 𝑣 rated is the turbine rated wind speed in m/s, and 𝑘 p = 𝜌 𝑎𝑖𝑟 𝐶 𝑝 2 ⁄ where 𝜌 𝑎𝑖𝑟 is the air density and 𝐶 𝑝 the power coefficient of the wind turbine, assumed constant.

All power flows of this multi-source system have to satisfy the power balance:

(3) 𝑃 load (𝑡) -𝑃 sto (𝑡) -𝑃 PV (𝑡) -𝑃 W (𝑡) -𝑃 grid (𝑡) = 0 ∀𝑡 ∈ [0, 𝑇]
where 𝑃 grid (𝑡) is the power from the grid, 𝑃 sto (𝑡) is the power from the storage system, and 𝑃 load (𝑡) is the power consumption of trains. Given the maximum amount of stored energy and power exchanged between the storage system and other substation components, lower and upper bounds apply for 𝑃 sto (𝑡) and 𝐸 sto (𝑡) as expressed in (4)-( 5):

(4) -𝑃 max ≤ 𝑃 sto (𝑡) ≤ 𝑃 max ∀𝑡 ∈ [0, 𝑇] (5) 0 ≤ 𝐸 sto (𝑡) ≤ 𝐸 max ∀𝑡 ∈ [0, 𝑇]
where the maximum amounts of power 𝑃 max in W and energy 𝐸 max in J define the storage capacity. At the end of the time range 𝑇, the amount of energy stored is supposed to return to its initial value 𝐸 ini for periodicity:

(6) 𝐸 sto (𝑇) = 𝐸 ini
The power 𝑃 sto (𝑡) and the energy 𝐸 sto (𝑡) can both express the temporal fluctuation of energy in the storage device. The stored energy is computed from the history of 𝑃 sto (𝑡) and the initial amount of energy in the storage system 𝐸 ini assuming that there is no loss in the storage system: [START_REF] Costa | A survey on Benders decomposition applied to fixed-charge network design problems[END_REF] 𝐸 sto (𝑡) = 𝐸 ini -∫ 𝑃 sto (𝜏) 𝑑𝜏 𝑡 τ=0

∀𝑡 ∈ [0, 𝑇]

Reciprocally, the power 𝑃 sto (𝑡) can be computed from the derivative of the stored energy 𝐸 sto (𝑡):

(8) 𝑃 sto (𝑡) = -𝑑𝐸 sto (𝑡) 𝑑𝑡 ⁄ ∀𝑡 ∈ [0, 𝑇]
OPTIMIZATION PROBLEMS In equations ( 1)-( 8), the values of the wind speed 𝑣 W (𝑡), the solar irradiance 𝑃 irr (𝑡), the turbine cut-off wind speed 𝑣 off , the turbine rated wind speed 𝑣 rated , and the power consumption of trains 𝑃 load (𝑡) are known. In order to solve the optimization problem, the time is discretized with a time step Δ𝑡 resulting in 𝑛 t = 𝑇 Δ𝑡 ⁄ time steps. Measurement of the mean power consumption of trains is made every 10 minutes therefore, the time step is taken equal to Δ𝑡 = 600 seconds. The wind speed and solar irradiance are measured every hour and interpolations are used to compute values every 600 seconds. The time range is considered to be one year as the slowest variation is the one of the solar irradiance, resulting in 𝑛 t = 52560 time steps. Now, it is possible to write the discretized forms of equations ( 7) and ( 8):

(7d) 𝐸 sto (𝑡 = 𝑖 𝛥𝑡) = 𝐸 sto,i = 𝐸 ini -∆𝑡 ∑ 𝑃 sto,j ∀𝑖 ∈ ⟦1, 𝑛 t ⟧ i j=1 (8d) 𝑃 sto (𝑡 = 𝑖 𝛥𝑡) = 𝑃 sto,i = [𝐸 sto,i-1 -𝐸 sto,i ]/∆𝑡 ∀𝑖 ∈ ⟦1, 𝑛 t ⟧ with 𝐸 sto,0 = 𝐸 ini
The objective function is expressed in [START_REF] Geoffrion | Generalized benders decomposition[END_REF] for 𝑛 y = 20 years operating mode, that is the lifespan of equipments, by considering investments cost 𝐶 inv of all devices and also the non-linear energy cost for purchasing or selling electricity. Bold letters are used for vectors. [START_REF] Geoffrion | Generalized benders decomposition[END_REF] OF(𝑆 PV , 𝑆 W , 𝑃 max , 𝐸 max , 𝑷 𝐠𝐫𝐢𝐝 ) = 𝐶 inv (𝑆 PV , 𝑆 W , 𝑃 max , 𝐸 max ) + 𝑛 y ∑ 𝑐𝑜𝑠𝑡(𝑃 grid,i )

n t i=1
The decision variables are the photovoltaic panels surface 𝑆 PV , the swept surface of wind turbines 𝑆 W , the maximum amount of power 𝑃 max , the maximum amount of energy in the storage system 𝐸 max , the initial amount of energy in the storage system 𝐸 ini , and a high dimension decision variable for the control of the storage system that can be the power 𝑷 𝐬𝐭𝐨 or the energy 𝑬 𝐬𝐭𝐨 . Therefore, two optimization problems can be expressed.

In the first formulation, the decision variables for the control of the storage system are the powers 𝑃 sto (𝑡) and (7d) is used to rewrite the discretized forms of equations ( 4)-( 6):

(4d) -𝑃 max ≤ 𝑃 sto,i ≤ 𝑃 max ∀𝑖 ∈ ⟦1, 𝑛 t ⟧ (10) 0 ≤ 𝐸 ini -∆𝑡 ∑ 𝑃 sto,j i j=1 ≤ 𝐸 max ∀𝑖 ∈ ⟦1, 𝑛 t ⟧ (11) ∑ 𝑃 sto,j n t j=1 = 0
Equations ( 1)-( 3) are merged to express the power from the grid:

(12) 𝑃 grid,i = 𝑃 load,i -𝑃 sto,i -𝑆 PV 𝑃 irr,i -𝑘 p 𝑆 W 𝑣 adj,i 3 ∀𝑖 ∈ ⟦1, 𝑛 t ⟧
Therefore, the optimization problem with the power formulation is stated as:

(13)

min 𝑆 PV ,𝑆 W ,𝑃 max ,𝐸 max ,𝐸 ini ,𝑷 𝐬𝐭𝐨 OF ( 𝑆 PV , 𝑆 W , 𝑃 max , 𝐸 max , 𝑷 𝐥𝐨𝐚𝐝 -𝑷 𝐬𝐭𝐨 -𝑆 PV 𝑷 𝐢𝐫𝐫 -𝑘 p 𝑆 𝑊 𝒗 𝐚𝐝𝐣 𝟑 ) 𝑠. 𝑡. (4d), (10), (11) 𝑆 PV ∈ 𝑅 + , 𝑆 W ∈ 𝑅 + , 𝑃 max ∈ 𝑅 + , 𝐸 max ∈ 𝑅 + , 𝐸 ini ∈ 𝑅 + , 𝑷 𝐬𝐭𝐨 ∈ 𝑅 n t
In the second formulation, the decision variables for the control of the storage system are the energies 𝐸 sto (𝑡) and ( 8) is used to rewrite equations ( 4)-( 6) and ( 12): To summarize, the optimization problems for both formulations include 𝑛 t + 5 = 52565 design variables, one non-linear objective function, one linear equality constraint, and 4 𝑛 t = 210240 linear inequality constraints. It is to mention that the left inequality in (5d) for the optimization problem with the energy formulation that is 0 ≤ 𝐸 sto,i can be ignored since the lower bound of energy is zero. This removes 𝑛 t = 52560 inequality constraints.

(5d) 0 ≤ 𝐸 sto,i ≤ 𝐸 max ∀𝑖 ∈ ⟦1, 𝑛 t ⟧ (6d) 𝐸 sto,n t = 𝐸 ini (14) -𝑃 max ∆𝑡 ≤ 𝐸 sto,i-1 -𝐸

LINEARIZATION OF COSTS

Several pricing systems exist in the electricity market. For instance, the energy purchase price may vary according to the hour and the month while the selling price of renewable energies is often constant. In this case, the energy cost is a piecewise linear function of the power consumption or production according to the hour, date, and direction of the power flow. On the other hand, the delivery cost of electricity is a non-linear function that is not investigated in this paper but constitutes a challenging prospect. Linear programming is used to solve the optimization problem for both formulations. In Matlab programming, specific matrices are required to express the equality and inequality constraints. Due to a time step of 10 minutes, each temporal variable induces 52560 unknowns for one year cycle thus sparse matrices are employed to reduce the memory requirements. In Table 2, both formulations are compared according to the size of matrices required to express the equality and inequality constraints such as the amount of memory required to store them by using full or sparse matrices. Fig. 2.a. and Fig. 2.b. present respectively the nonzero elements of equality and inequality matrices for power and energy formulations during one day, i.e. 144 time steps. On Fig. 2.a., it can be seen that the power formulation induces a higher number of nonzero elements which makes difficult the storage of matrices for a high number of decision variables. By employing energy formulation, the memory required is greatly reduced and the optimization of HRPS is achievable in about one minute on a laptop computer for a complete year. The implementation is made with Matlab 2015b and the dual-simplex algorithm in linprog function from the Optimization toolbox. The tolerance has to be reduced to 10 -9 to achieve a good precision for a complete year. Other options take the default values. The linearization of costs appears to be a good way to get an approximate solution of the HRPS optimization problem but the prospect is to optimize with the real costs that are non-linear functions. Therefore, new strategies are required for the non-linear optimization with a high number of decision variables.

COLLABORATIVE OPTIMIZATION

The principle of collaborative optimization is to decompose the complete problem into sub-problems, each of them having a lower number of decision variables [START_REF] Zadeh | Metamodel-based collaborative optimization framework[END_REF]. Sequential quadratic programming (SQP) and active-set, among others, are algorithms able to accurately solve non-linear optimization problems with numerous active constraints at optimum in a reduced number of evaluations. Unfortunately, the number of variables is limited to hundreds. Interior-point method (IPM) can extend the number of variables to thousands but the precision of the solution may decrease because of the barrier method used to handle constraints. Attempts to solve problems (13) and ( 16) by using SQP and IPM have been carried out with an increasing number of time steps (Fig. 3). It appears that the convergence of SQP is not achieved in an affordable computing time beyond 432 time steps (3 days) while the limit of IPM is 1296 time steps (9 days). For an increasing number of time steps, IPM shows better solution and lower computing time compared to SQP. However, only the solution found by SQP for a time range of one day is as good as LP. Therefore, the decomposition of the HRPS problem is guided by the number of variables in each subproblem. It is proposed to decompose in two levels: year and day. This is shown in Fig. 4 where the year level is at the top and the day level at bottom. For the sake of simplicity, the year is reduced to 7 days in this figure. At the year level, the design variables are the four dimensional parameters since their values are not time-dependent. At the day level, the idea is to decompose the complete year in 𝑛 d = 365 days with optimizations made in parallel. This is enabled by introducing additional linking variables at the year level that are the amount of energy in the storage system at the beginning of each day. Additional constraints are added at the day level to impose the amount of energy stored at the beginning and the end of each day. The number of design variables for one day is 𝑛 t 𝑛 d ⁄ = 144, that is also the number of time steps for this period, i.e. 24 hours multiplied by 6 time steps per hour. At the year level, the number of design variables is 𝑛 d + 4 = 369 corresponding to 𝑛 d = 365 linking variables and 4 dimensional parameters.

The minimum energy costs (24) at the day level given the values of the dimensional and linking variables are found by 𝑛 d optimizations in parallel and added as a second member in the objective function ( 25) at the year level.

(24) As the objective function at the year level includes the results of optimizations at day level, those lasts are nested in the year level model as shown in Fig. 5. This leads to a double-loop approach with high computing time and slow convergence. Two cases are considered for the inner-loop. In the first case (CO 1st case), the energy cost is a non-linear function and problems (24) are solved with non-linear programming (NLP) and distributed computing. In the second case (CO 2nd case), the energy cost is approximated by a piecewise linear function and problems (24) are solved with linear programming. For both cases, the outer-loop that minimizes (25) requires NLP because the energy cost for all days 𝐶 ener,d are implicit functions. The derivatives of 𝐶 ener,d relative to the variables 𝑆 PV , 𝑆 W , 𝑃 max , 𝐸 max , 𝑬 𝐢𝐧𝐢 are computed by using the finite difference technique. There is no guaranty of convergence for the outer-loop due to the fact that 𝐶 ener,d may not be convex. However, all the constraints are linear and the feasible space is thus convex. This double-loop approach appears to have unaffordable computing time even for a small time range in the first case while the same solution as LP is found in a reasonable computing time for few days in the second case (Fig. 6). The collaborative optimization with linear programming in the innerloop (2 nd case) will be investigated for a complete year in the results part.

𝐶 ener,d = min 𝑬 𝐬𝐭𝐨 ∑ 𝑐𝑜𝑠𝑡(𝑃 load,i -𝑃 PV,i -𝑃 W,i -[𝐸 sto,i-1 -𝐸 sto,i ]/∆𝑡) d n t n d ⁄ 𝑖=1+(d-

MULTILEVEL OPTIMIZATION

The advantage of multilevel optimization is to avoid the nested optimization by using an approach where the optimizations are solved sequentially. At the year level, the optimization is performed by using substitute models 𝐶 ener,d ~ of the result of optimizations at the day level 𝐶 ener,d for each day. This way, the costly optimizations at the day level are replaced by fast evaluations of functions at the year level that are usually linear. Then, the optimizations are made at the day level and the substitute models are refined. The loop between levels stops when the discrepancies between the substitute models and the minimum energy costs are negligible. Figure 7 gives the flowchart of this approach and the main steps are detailed below. 24) is an implicit function of the design variables at the day level, the difficulty is to build a substitute model 𝐶 ̃ener,d = 𝑠 𝑑 (𝐸 max , 𝑃 max , 𝑆 PV , 𝑆 W , 𝐸 ini,d , 𝐸 ini,d+1 ) that is an explicit function of the design variables at the year level.

Based on the theory of duality, the Benders decomposition [START_REF] Costa | A survey on Benders decomposition applied to fixed-charge network design problems[END_REF], [START_REF] Benders | Partitioning Procedures for Solving Mixed-Variables Programming Problems[END_REF] takes advantage of the information on the Lagrange multipliers to define a lower bound for the objective function of the sub-problems, as detailed in Appendix. When the objective function and the constraints of the sub-problems are linear, and the objective function and the constraints of the master problem (at the year level) are continuous, then Benders decomposition converges in a finite number of steps [START_REF] Benders | Partitioning Procedures for Solving Mixed-Variables Programming Problems[END_REF]. If the whole problem has a finite optimum solution, Benders decomposition converges to an optimum solution. Hence, the application of Benders decomposition on problem (22) or (23) gives an optimal solution. If the objective function or the constraints of the sub-problems are no longer linear, the problem can be solved using a generalized Benders decomposition [START_REF] Geoffrion | Generalized benders decomposition[END_REF]. This method is proven to be convergent when the sub-problems consist in minimizing a convex function over a convex constraint set [START_REF] Geoffrion | Generalized benders decomposition[END_REF].

In the HRPS optimization problem, the objective functions are the costs of energy for all days. In the subproblems, the decision variables 𝐸 max , 𝑃 max , 𝑆 PV , 𝑆 W , 𝑬 𝐢𝐧𝐢 at the year level are constant. The constraints (5d), (6d), (14), and (15) are rewritten such as these variables appear in their RHS: The algorithm starts with an initial value of the variables at the year level. After solving the sub-problems, the Lagrange multipliers are available for all constraints: (𝜆 ineq1,i ) The algorithm stops if the gap between the sums of energy costs found at both levels is small enough: where 𝜀 is the precision chosen equal to 1€ what is about 10 -7 relative to the minimum value of the objective function.

(

RESULTS AND COMPARISON

As already mentioned, the energy cost is assumed piecewise linear in order to compare the results of collaborative optimization and Benders decomposition to the solving of the complete optimization problem with linear programming. As expected, this last gives the lower total cost in the shortest time and is considered as the reference for comparison.

Figures 8 show the value of the objective function and the computing time for the three approaches when the number of days considered in the problem increases. LP is linear programming, CO is collaborative optimization with linear programming for the inner-loop, and BD is the Benders decomposition. The cost is normalized with respect to LP value. It can be seen that BD gives the same solution as LP with about 5 times more computing time but this ratio decreases as the number of days considered increases. CO fails to find the solution when considering more than 120 days. The time required by CO is about one thousand times more than LP.

CONCLUSIONS

Benders decomposition appears to be very efficient for the optimization of HRPS on a large time range taking into account the control of the storage system. For a piecewise linear cost function, the design is performed in a few minutes on a laptop computer for a complete year. The solution is the same as the one found with linear programming and the computing time is about 5 times higher only. Collaborative optimization fails to find the same solution and its computing time is at least one hundred times higher.

As the energy cost function is not linear in reality, the linear programming cannot be used. Therefore, one prospect is to apply Benders decomposition to design the HRPS with a non-linear energy cost function as it is not possible to solve this optimization problem as a whole.

Figure 1 .

 1 Figure 1. Power flows in hybrid railway power substation and dimensional decision variables.

Figure 2 .

 2 Figure 2. Comparison of matrices for one day. Each column is a variable and each row is a constraint.

Figure 3 .

 3 Figure 3. Comparison of algorithms to solve HRPS complete problem: cost (left) and time (right). The cost is normalized with respect to LP value.

Figure 4 .

 4 Figure 4. Decomposition of the HRPS in two levels.

Figure 5 . 1 𝐶 1 Figure 6 .

 5116 Figure 5. Collaborative optimization (double-loop) structure.

Figure 7 .

 7 Figure 7. Multilevel optimization (sequential) structure. As the objective function 𝐶 ener,d = 𝑓(𝐸 sto,1+(d-1) n t n d ⁄ , ⋯ , 𝐸 sto,d n t n d ⁄ ) expressed in (24) is an implicit function of the design variables at the day level, the difficulty is to build a substitute model 𝐶 ̃ener,d = 𝑠 𝑑 (𝐸 max , 𝑃 max , 𝑆 PV , 𝑆 W , 𝐸 ini,d , 𝐸 ini,d+1 ) that is an explicit function of the design variables at the year level.

Figure 8 .

 8 Figure 8. Comparison of approaches to solve HRPS problem: cost (left) and time (right).

  𝑆 PV , 𝑆 W , 𝑃 max , 𝐸 max , 𝑷 𝐥𝐨𝐚𝐝 -∆𝑬 𝐬𝐭𝐨 /∆𝑡 -𝑆 PV 𝑷 𝐢𝐫𝐫 -𝑘 p 𝑆 W 𝒗 𝐚𝐝𝐣 𝟑 ) 𝑠. 𝑡. (5d), (6d), (14) 𝑆 PV ∈ 𝑅 + , 𝑆 W ∈ 𝑅 + , 𝑃 max ∈ 𝑅 + , 𝐸 max ∈ 𝑅 + , 𝐸 ini ∈ 𝑅 + , 𝑬 𝐬𝐭𝐨 ∈ 𝑅 +n t

sto,i ≤ 𝑃 max ∆𝑡 ∀𝑖 ∈ ⟦1, 𝑛 t ⟧ with 𝐸 sto,0 = 𝐸 ini (15) 𝑃 grid,i = 𝑃 load,i -∆𝐸 sto,i /∆𝑡 -𝑆 PV 𝑃 irr,i -𝑘 p 𝑆 W 𝑣 adj,i 3 ∀𝑖 ∈ ⟦1, 𝑛 t ⟧ with ∆𝐸 sto,i = 𝐸 sto,i-1 -𝐸 sto,i Therefore, the optimization problem with the energy formulation is stated as: (16) min 𝑆 PV ,𝑆 W ,𝑃 max ,𝐸 max ,𝐸 ini ,𝐸 sto OF (

Table 1 . Description of both formulations for linearized optimization problem Formulation power energy

 1 The power from the grid 𝑷 𝐠𝐫𝐢𝐝 is decomposed in two parts, the power received 𝑷 𝐠𝐫𝐢𝐝 𝑆 PV ∈ 𝑅 + , 𝑆 W ∈ 𝑅 + , 𝑃 max ∈ 𝑅 + , 𝐸 max ∈ 𝑅 + , 𝐸 ini ∈ 𝑅 + 𝑬 𝐬𝐭𝐨 ∈ 𝑅 +n t , 𝑷 𝐠𝐫𝐢𝐝The linearized optimization problems for both formulations include 3 𝑛 t + 5 = 157685 decision variables, one linear objective function, 𝑛 t + 1 = 52561 linear equality constraints, and 4 𝑛 t = 210240 linear inequality constraints. Table1summarizes the decision variables and constraints for both formulations.

	+	and the power sent

-∀𝑖 ∈ ⟦1, 𝑛 t ⟧ This allows changing the non-linear energy cost function by a piecewise linear function: } 𝑠. 𝑡. (5d), (6d), (14), (21) + ∈ 𝑅 +n t , 𝑷 𝐠𝐫𝐢𝐝 -∈ 𝑅 -n t

Table 2 . Comparison of memory requirements for both formulations (𝒏 𝐭

 2 = 𝟓𝟐𝟓𝟔𝟎)

	Formulation	Expression	power	Memory	Expression	energy	Memory
	Size of matrices Number of nonzero	(5 𝑛 t + 1)(3 𝑛 t + 5) 14 𝑛 t + 𝑛 t 2	309 GB (full) 62 GB (sparse)	(4 𝑛 t + 1)(3 𝑛 t + 5) 14 𝑛 t + 4	247 GB (full) 17 MB (sparse)

  𝑆 PV , 𝑆 W , 𝑃 max , 𝐸 max , 𝑬 𝐢𝐧𝐢 𝐸 ini,1 𝐸 ini,2 𝐸 ini,d 𝐸 ini,d+1 𝐸 ini,n d 𝐸 ini,1

		initialize
	solve		solve		solve
	(24) for day	…	(24) for day	…	(24) for day
	#1		#d		#𝑛 d
	C ener,1	C ener,d	C ener,n d
	𝜆 eq,1		𝜆 eq,d	𝜆 eq,n d
	𝜆 ineq,1	𝜆 ineq,d	𝜆 ineq,n d
		solve master problem (35)
			𝑪 𝐞𝐧𝐞𝐫 ~
	𝑆 PV , 𝑆 W , 𝑃 max , 𝐸 max , 𝑬 𝐢𝐧𝐢		(36)?
		no

  𝐸 sto,i-1 -𝐸 sto,i + 𝑃 grid,i ∆𝑡 = [𝑃 load,i -𝑆 PV 𝑃 irr,i -𝑘 p 𝑆 W 𝑣 adj,i

	26)	𝐸 sto,i ≤ 𝐸 max	∀ 𝑖 ∈ ⟦1 + (𝑑 -1) 𝑛 t 𝑛 d ⁄ , 𝑑 𝑛 t 𝑛 d ⁄ ⟧
	(27)	-𝐸 sto,1+(d-1)n t n d ⁄	≤ 𝑃 max ∆𝑡 -𝐸 ini,d
	(28)	𝐸 sto,i-1 -𝐸 sto,i ≤ 𝑃 max ∆𝑡	∀ 𝑖 ∈ ⟦2 + (𝑑 -1) 𝑛 t 𝑛 d ⁄ , 𝑑 𝑛 t 𝑛 d ⁄ ⟧
	(29)	𝐸 sto,1+(d-1)n t n d ⁄	≤ 𝑃 max ∆𝑡 + 𝐸 ini,d
	(30)	𝐸 sto,i -𝐸 sto,i-1 ≤ 𝑃 max ∆𝑡	∀ 𝑖 ∈ ⟦2 + (𝑑 -1) 𝑛 t 𝑛 d ⁄ , 𝑑 𝑛 t 𝑛 d ⁄ ⟧
	(31)	-𝐸 sto,1+(d-1)n t n d ⁄ -𝑘 p 𝑆 W 𝑣 adj,1+(d-1)n t n d + 𝑃 grid,1+(d-1)n t n d ⁄ ∆𝑡 = 𝑃 load,1+(d-1)n t n d ⁄ ⁄ 3 ∆𝑡 -𝐸 ini,d	⁄ ∆𝑡 -𝑆 PV 𝑃 irr,1+(d-1)n t n d	∆𝑡
	(32)				3 ] ∆𝑡
			∀ 𝑖 ∈ ⟦2 + (𝑑 -1) 𝑛 t 𝑛 d ⁄ , 𝑑 𝑛 t 𝑛 d ⁄ ⟧
	(33)	𝐸 sto,n t n d ⁄	= 𝐸 ini,d+1

  𝜆 eq2 for (33). The multiplier is a measurement of the expected decrease of the energy cost for an increase of the RHS of a constraint in the sub-problem. As all RHS are expressed explicitly as linear functions of the variables at the year level, the lower bounds for all energy costs 𝑪 ̃𝐞𝐧𝐞𝐫 are deduced and introduced as additional constraints (34) in the optimization problem at the year level as explained in the Appendix.The master problem (35) at the year level is then solved with these additional constraints, one for each day, and new values for the design variables are deduced such as the values of the expected energy costs for all days, 𝐶 ̃ener,d . At each iteration, 𝑛 d new constraints on the energy costs are added. 𝑆 W ,𝑃 max ,𝐸 max ,𝑬 𝐢𝐧𝐢 ,𝑪 ̃𝐞𝐧𝐞𝐫 𝐶 inv (𝑆 PV , 𝑆 W , 𝑃 max , 𝐸 max ) + 𝑛 y ∑ 𝐶 ener,d ñd d=1 𝑠. 𝑡. (34) ∀ 𝑑 ∈ ⟦1, 𝑛 d ⟧ 𝑆 PV ∈ 𝑅 + , 𝑆 W ∈ 𝑅 + , 𝑃 max ∈ 𝑅 + , 𝐸 max ∈ 𝑅 + , 𝑬 𝐢𝐧𝐢 ∈ 𝑅 +n d , 𝑪 𝐞𝐧𝐞𝐫 ~∈ 𝑅 n d

		min
	(35)	𝑆 PV ,
					1+(d-1)n t n d ⁄ ≤ i ≤ d n t n d ⁄	for the inequality
	constraint (26), 𝜆 ineq2,1+(d-1) n t n d ⁄	for (27), (𝜆 ineq2,i ) 2+(d-1)n t n d ⁄ ≤ i ≤ d n t n d ⁄	for (28), 𝜆 ineq3,1+(d-1) n t n d ⁄	for
	(29), (𝜆 ineq3,i ) 2+(d-1)n t n d ⁄ ≤ i ≤ d n t n d ⁄	for (30), 𝜆 eq1,1+(d-1)n t n d ⁄	for the equality constraint (31),
	(𝜆 eq1,i ) 2+(d-1)n t n d ⁄ ≤ i ≤ d n t n d ⁄ for (32), and (34) -𝐶 ̃ener,d -∑ 𝜆 ineq1,i d n t n d ⁄ i=1+(d-1)n t n d ⁄ 𝐸 max -∑	d n t n d ⁄ i=1+(d-1)n t n d ⁄	(𝜆 ineq2,i + 𝜆 ineq3,i )	𝑃 max ∆𝑡 +
		(𝜆 ineq2,1+(d-1)n t n d ⁄	-𝜆 ineq3,1+(d-1)n t n d ⁄	+ 𝜆 eq1,1+(d-1)n t n d ⁄ ) 𝐸 ini,d -𝜆 eq2 𝐸 ini,d+1
		+∆𝑡 [𝜆 eq1,1+(d-1)n t n d ⁄ ≤ [𝜆 eq1,1+(d-1)n t n d ⁄	⋯ 𝜆 eq1,d n t n d ⁄ ] {[ 𝑃 load,1+(d-1)n t n d 𝑃 irr,1+(d-1)n t n d ⁄ ⋮ 𝑃 irr,d n t n d ⁄ ⁄ ⁄ 𝑃 load,d n t n d ⋯ 𝜆 eq1,d n t n d ⁄ ] [ ⋮ ] ∆𝑡 ] 𝑆 PV + 𝑘 p [	⁄ 𝑣 adj,1+(d-1)n t n d 3 ⋮ ⁄ 3 𝑣 adj,d n t n d	] 𝑆 W }

APPENDIX: BENDERS DECOMPOSITION

We consider a linear optimization problem with two sets of variables. The complete formulation is: (A1) min

x,y 𝑐 t 𝑥 + 𝑑 t 𝑦 (A2) s.t.

𝐴 𝑥 + 𝐵 𝑦 ≤ 𝑏 𝐴 eq 𝑥 + 𝐵 eq 𝑦 = 𝑏 eq 𝐷 𝑥 ≤ 𝑒

In order to perform the Benders decomposition, the variables 𝑥 are assumed to belong to the master problem while the variables 𝑦 are in the sub-problem. Classically, the variables in the master problem are integer and the ones in the sub-problem are continuous but Benders decomposition can also be applied to problems with only continuous variables.

At the beginning, an initial value 𝑥 0 for the variables in the master problem is given. The first iteration (𝑘 = 1) starts by solving the sub-problem expressed as:

After solving the sub-problem (A4)-(A6), a marginal cost or Lagrange multiplier is available for each constraint: 𝜆 ineq for inequality constraint, 𝜆 eq for equality constraint, 𝜆 ub for upper bound, and 𝜆 lb for lower bound. The multiplier is a measurement of the expected decrease of the objective for an increase of the RHS of the constraint. The multiplier for an inequality constraint is positive if the constraint is active else it is equal to zero. For equality constraint, the multiplier is positive or negative.

The range constraint (A6) is transformed in two inequality constraints to standardize the notation:

The second step in the iteration is to solve the master problem expressed as: The iterations are continued until a stopping criterion is fulfilled. This stopping criterion is usually:

where 𝜀 is the requested precision. A new constraint on 𝑓 ~ is added to the set of inequality constraint at each iteration.


John Dodds

Fabienne Espitalier
email: fabienne.espitalier@enstimac.fr

Olivier Louisnard
email: olivier.louisnard@enstimac.fr

Romain Grossier

René David

Myriam Hassoun

Fabien Baillon

Cendrine Gatumel

RAPSODEE Nathalie Lyczko

The Effect of Ultrasound on Crystallisation-Precipitation Processes: Some Examples and a New Segregation Model

Keywords: crystallization, nucleation, precipitation, segregation model, ultrasound

This paper discusses the effects of ultrasound on the production of particles by precipitation and crystallization. Examples are given from the formation of crystals of BaSO 4 , K 2 SO 4 , TiO 2 and sucrose. It is shown that ultrasound reduces the induction time, narrows the width of the metastable zone and leads to the production of more, finer, and more uniform crystals in some cases. The reasons for these effects of ultrasound on the nucleation of crystals are discussed and a possible mechanism is presented.

Introduction

Ultrasound has been used for over 20 years to control the formation of solids in crystallization processes. Experiments have been reported where supersaturation has been generated by cooling (a solution or a melt), by addition of a solvent or by reaction. A large variety of products have been studied: inorganic (titanium dioxide, ammonium sulfate, aluminium alloy, zeolite...), organic (sugars, pharmaceutical drug or excipients, cocoa butter...) and ice. All this work shows that ultrasound has an effect on nucleation by shifting the size distribution towards small particles and modifying the morphology or the polymorphs synthetized. [START_REF] Ruecroft | Sonocrystallization: the use of ultraound for improved industrial crystallization[END_REF] presented two large scale industrial processes based on a flow cell with multiple transducers for a high insonation and a device using a combination of a vortex-mixing with a transducer for preparation of micron sized particles. However, the mechanism or the link between nucleation and ultrasonic cavitation remains unclear and uncertain [START_REF] Ruecroft | Sonocrystallization: the use of ultraound for improved industrial crystallization[END_REF] and only a few papers have tried to quantify the observed phenomena. For example in a recent paper, [START_REF] Virone | Ultrasound for reproducible nucleation in batch crystallization[END_REF] attempted to correlate the collapse pressure of cavitation bubble with supersaturation and therefore with the nucleation rate. In their paper, they presented results for crystallization induction time determined in a system with a well-defined pressure profile avoiding uncontrolled reflection. It was shown that the induction times calculated from the nucleation as a function of local pressure and growth rate were much lower than those determined experimentally: 1 s instead of 30 min. In the first part of this paper some characteristic results will be presented and in the second part, the modeling of molecular segregation by pressure gradients around a single cavitation bubble will be described. This is proposed as a mechanism for the effect of ultrasound on particle generation and some with practical consequences will be presented.

The main effects of ultrasound on the production of particles by precipitation and crystallization will be presented with examples from previous papers [START_REF] Bird | Transport phenomena[END_REF][START_REF] Gatumel | Nucleation control in precipitation processes by ultrasound[END_REF][START_REF] Grossier | Mixture segregation by an inertial cavitation bubble[END_REF] involving BaSO 4 , K 2 SO 4 , TiO 2 and sucrose [START_REF] Gatumel | Nucleation control in precipitation processes by ultrasound[END_REF][START_REF] Lyczko | Effect of ultrasound on the induction time and the metastable zone widths of potassium sulphate[END_REF][START_REF] Baillon | Procédé de synthèse du dioxyde de titane: analyse et modélisation des solutions Titane-Sulfate; influence des ultrasons sur la précipitation[END_REF], Hassoun et al. 2003). In these experiments, the dissipated ultrasound power measured by calorimetry ranged from 0.01 to 0.15 W/g solution, except in the case of continuous precipitation of barium sulfate (1.6 W/g sol.). The frequency of ultrasound was 20-23 kHz. The results obtained from all these different studies, involving different products and using different ultrasonic transducers showed that ultrasound makes crystallization more reproducible. and that the major effect of ultrasound is on primary or of more reproducible and secondary nucleation.

Barium Sulfate

In the study of barium sulfate precipitation, two modes of precipitation were used: semibatch in a stirred vessel of 1.5 L and in continuous mode [START_REF] Gatumel | Nucleation control in precipitation processes by ultrasound[END_REF]. The continuous mode was used to separate the nucleation step from the crystal growth step. The experimental equipment comprised two vessels in series: a tube cell (0.050 l) and a stirred vessel (0.785 l). Ultrasound was emitted respectively by a tip transducer in the small cell and from a large bottom surface in the 0.785 l vessel (Figure 1).

Figure 2 shows the supersaturation ratio S at the outlet of the growth vessel as function of the supersaturation ratio at the outlet of the nucleation cell with or without ultrasound. The supersaturation ratio S is the ratio between the concentration and the equilibrium concentration. The concentrations are calculated from conductivity measurements. S was always lower at the outlet of the growth vessel than at the outlet of the nucleation cell due to the consumption of solute by growth or secondary nucleation in the growth vessel. Three zones can be identified: zone 1 corresponds to experiments with ultrasound, zone 2 to experiments with ultrasound in the growth vessel and zone 3 with ultrasound in the nucleation cell and eventually in the growth vessel. In all cases, the application of ultrasound brings about a decrease in S and in the mean surface diameter (Figure 3). In this figure, the numbers of crystals at the exit of two vessels are given. In presence of ultrasound the number of crystals increases ten fold. This increase is directly linked with the decrease of the mean diameter. In addition on Figure 2, shows the SEM pictures of crystals obtained at the exit of the cell with and without ultrasound. This indicates that for a dissipated power greater than 1.5 W/g sol, ultrasound modifies the morphology of crystals to make more uniform crystals than without ultrasound. Figure 4 presents the surface size distribution obtained at the exit of the nucleation cell (a) with and without ultrasound and at the exit of the growth vessel with or without ultrasound in the nucleation cell (b).

The application of ultrasound narrows the particle size distribution and gives a monomodal distribution.

Potassium Sulfate

Induction times have been measured at 15 °C in a vessel of 200 ml and 1 L for potassium sulfate [START_REF] Lyczko | Effect of ultrasound on the induction time and the metastable zone widths of potassium sulphate[END_REF]. The experimental equipment (Figure 5) are thermostated vessels (0.2 L and 1L) respectively with a magnetic stirrer at a constant rotation speed of about 500 rpm and with a three-blade propeller (Mixel TT) at a constant rotation speed about 600 rpm. Ultrasound is respectively applied at the top of the vessel by a stainless ultrasound source with a titanium tip and at the base of the vessel with a cup-horn source (from Sinaptec Co.). Two ultrasonic powers levels have been measured: 0.05 and 0.12 W/g solution for the tip and 0.03 W/g solution for a cup horn source. Saturated aqueous solutions were prepared from potassium sulfate crystals (synthesis grade 99,9%) and distilled water. The solution is rapidly cooled from the saturation temperature plus 5 °C to a final temperature of about 15 ± 0.4 °C. When this final temperature is reached, the conductivity solution is monitored in order to detect the appearance of the solid. The conductivity and the temperature of the solution were measured with an instrumental accuracy of about ± 0.1 mS/cm and ± 0.1 °C, respectively. Ultrasound is applied at the end of the cooling phase. In the 0.2 L vessel, the initial saturation temperatures ranged from 22 to 25 °C and from 19.5 to 23 °C for experiments without and with ultrasound, respectively. In the 1 L vessel, the initial saturation temperatures ranged from 24 to 30 °C for experiments without ultrasound and from 19 to 25 °C for experiments with ultrasound. Figure 6 presents the induction times at 15 °C for the two vessels as function of the absolute supersaturations, and it is seen that the induction time is greatly reduced in the presence of ultrasound [START_REF] Lyczko | Effect of ultrasound on the induction time and the metastable zone widths of potassium sulphate[END_REF]. For instance, in the case of 0.2 L vessel, for absolute supersaturations higher than 0.018 g de solid/g water, ultrasound has no influence on the induction time. For absolute supersaturations lower than 0.012 g de solid/g water, the induction time decreases as the ultrasonic power increases: the effect of ultrasound seems to be greater at low absolute supersaturation. The experimental procedure for crystal recovery is summarized on Figure 7. For experiments carried out in the first apparatus, the magnetic stirrer is replaced by a mechanical stirrer when the crystals start to appear. This change is necessary to maintain the potassium sulfate crystals in suspension in the solution. For experiments without ultrasound, crystals are left growing for only 10 minutes, as they become too large to be analyzed after that duration. For experiments with ultrasound, the ultrasound is stopped after the nucleation step and crystals are left growing for one hour. For all experiments in the second apparatus, crystals grow for one hour from the moment they are detected,. The suspension is then filtered with a 0.45 lm pore diameter filter and the crystals are dried in an oven at 100 °C during 24 hours. The crystals numbers have been measured by a particle counter (Coulter Multisizer II). Figure 8 shows, the number of crystals per unit volume against the absolute supersaturation for experiments with and without ultrasound, in the two experimental apparatus. In each vessel more crystals are formed in presence of ultrasound. For an equivalent absolute supersaturation (0.015 g K 2 SO 4 /g water), the number of crystals is about 3.10 8 #/m 3 without ultrasound and higher than 2.5 10 10 #/m 3 with ultrasound. Thus, ultrasound enhances the production of crystals. For experiments without ultrasound, more crystals are formed in the 0.2 L vessel than in the 1 L vessel and the same observation can be made for experiments with ultrasound. Therefore, the tip transducer seems more efficient that the "cup-horn".

For the two following studies (titanium dioxide and sucrose crystallization), the ultrasonic power reported throughout the paper is the active power indicated by the sonotrode generator and therefore is only an indication of the power really transmitted to the solution.

Titanium Dioxide

In the case of titanium dioxide, the induction time is about 30-35 min without ultrasound, is less than 25 min for a power of 60 W and less than 10 min for a power of 100 W at 100 °C [START_REF] Baillon | Procédé de synthèse du dioxyde de titane: analyse et modélisation des solutions Titane-Sulfate; influence des ultrasons sur la précipitation[END_REF]. For this product, the morphology is unchanged (Figure 9). In order to identify the crystals obtained, the X-ray diffraction patterns for titanium dioxide powders collected at the end of precipitation have been acquired in a range of 2 l from 20°to 75°by steps of 0.020°. For every point, the intensity data have been collected over 5 s. Figure 10 and Figure 11 show patterns obtained for powders produced by the hydrolysis of moderately and highly supersaturated solutions (respectively around 3000 and 8000), with and without ultrasound. The patterns have been deliberately shifted vertically for easier reading. Without ultrasound, it is seen that the pattern is characteristic of a pure anatase TiO 2 powder. On the other hand, with ultrasound, the XRD pattern also reveals the presence of a certain quantity of rutile within a matrix of anatase at low and high supersaturation. In order to clarify the influence of supersaturation, Figure 12 shows the refined XRD patterns (fine acquisition on the restricted 2 l range between 20°and 32°by steps of 0.005°) for two experi-ments with ultrasound, respectively at low and high supersaturation. It appears that in the former case, the product is mainly anatase, as shown previously, while in the latter, the pattern is clearly biased toward the rutile reference peak. Thus, for strong supersaturation, the application of ultrasound orients markedly the precipitation towards the rutile form. Deconvolution of the peak leads to an estimation of 60 % for the mass-fraction of rutile in the final powder.

Sucrose

Sucrose crystallization has been studied in batch and continuous mode [START_REF] Hassoun | Criblage des paramètres influant sur la cristallisation par refroidissement d'un produit organique assistée par ultrasons en cristallisoir discontinu et continu[END_REF], Hassoun et al. 2003).

Only results concerning continuous mode are reported in this paper. The crystallization was performed at 40 °C with a feed concentration of 76 % (w/w) and a feed flowrate of 0.15 l/min. Under these operating conditions, no crystals appeared in the 1.5 L crystallizer and in this particular case, ultrasound can be used to induce nucleation by a brief insonation or applied continuously in order to maintain a low size. The latter mode was chosen. The size distributions obtained in the steady state have been measured by acoustic attenuation (Ultrasizer, MALVERN) for different ultrasound powers. The operating conditions are given in Table 1. On Figure 13, the mean volume diameter d [START_REF] Gatumel | Nucleation control in precipitation processes by ultrasound[END_REF][START_REF] Bird | Transport phenomena[END_REF] and the variation of concentration between the entrance and the exit of the crystallizer DC E-S are reported. The d [START_REF] Gatumel | Nucleation control in precipitation processes by ultrasound[END_REF][START_REF] Bird | Transport phenomena[END_REF] decreases for a ultrasonic power between 25 and 35 W and stays constant for ultrasonic power higher than 35 W. DC E-S is highest at 25 W and therefore the absolute supersaturation is lower. At this power ultrasound would have only not act on nucleation, but growth will also be favored: crystals would therefore be bigger. For ultrasound powers higher than 25 W, the phenomenon of cavitation seems to increase with ultrasound power, leading to more nucleation, and as a result the creation of a bigger number of small crystals. For ultrasound power higher than 30 W, the size of the particles does not vary any more. This phenomenon was already noted during the study of the crystallization of barium sulfate under ultrasound. It can be explained by a phenomenon of autodamping of ultrasound with high ultrasonic power. When the ultrasonic power increases, the number of cavitation bubbles becomes more important and holds up the spread of the wave which then becomes less efficient. The semi-logarithmic plot of the population density according to the size of crystals for different ultrasound power (25-70 W) reveals some deviations from predictions from the model for a MSMPR type crystallizer (Mixed Suspension Mixed Product Removal) (Fig- ure 14). This phenomenon was already be noticed in the case of sucrose crystallizations [START_REF] Hartel | A kinetic study of the nucleation and growth of sucrose crystals in a continuous cooling, crystallizer[END_REF][START_REF] Hartel | Crystallization kinetics for the sucrosewater system[END_REF] and were attributed to a growth rate dependent on the size of crystals. The growth rate, G, has been described by the model of Abbeg, Stevens and Larson:

GL G 0 " 1 L G 0 s # b
where G 0 is the growth rate of small crystals and b describes the effect of size on the growth rate. With this type of law for the growth rate, the population balance for a MSMPR crystallizer can be rewritten:

nL n 0 " 1 L G 0 s # Àb exp " 1 À h 1 L=G 0 s 1Àb i 1 À b # :
The values of n 0 (the population density for the zero size), G 0 and b have been identified (Table 2). The experiment at 25 W is not taken in to account. The growth rates are all of the same order 6.42 ± 0.5 10 -10 m/s. This therefore implies that for an ultrasound power higher than 30 W the growth rate does not increase any more. We can note a slight difference on the growth rate calculated at 50 W. The calculated growth rate of 10 lm crystals is comparable with those obtained by other authors without ultrasound in the same range of supersaturation [START_REF] Hartel | A kinetic study of the nucleation and growth of sucrose crystals in a continuous cooling, crystallizer[END_REF][START_REF] Hartel | Crystallization kinetics for the sucrosewater system[END_REF]). The population densities obtained with ultrasound are 100 to 1000 times higher than those reported in the literature with higher residence time.

The density at zero n 0 increases according to applied ultrasound power. Therefore in most cases ultrasound would have an influence on nucleation rate: the higher the applied power the greater the number of created germs created.

For two experiments at 25 and 70 W, ultrasound was stopped after the steady state regime has been attained. The sucrose concentration was followed for one hour after this interruption (Figure 15 for P = 25 W): the sucrose concentration is constant. It seems that ultrasound initiates the crystallization, after this initiation ultrasound can be stopped or applied alternatively. In all the three cases (K 2 SO 4 crystallization by cooling, or BaSO 4 reactive crystallization), ultrasound brings about an increase in the number of crystals and a decrease in their size. All this at equivalent supersaturation and without modifying the crystalline structure if the solid does not present polymorphs. Other experiments with TiO 2 , not reported here also confirm these conclusions. In addition it was found that the use of ultrasound induced a change in crystalline structure if the solid presents polymorphs. Moreover, despite the sucrose solution viscosity being a 100 to 3000 times higher than the viscosity of water at the same temperature, the effects of ultrasound are found to be identical as previously mentioned: decrease of induction time and increase of the nucleation rate. Different explanations for these effects caused by ultrasound have been proposed in the literature: -the direct effect of high pressure arising around cavitation bubbles, or the effect of rapid local cooling rates-, however, none

Segregation Model

The apparent effects of ultrasonic cavitation on nucleation are the same that one would obtain by increasing supersaturation: more numerous and small crystals, nu-cleating more rapidly. A naive question would therefore be: could cavitation increase supersaturation? In classical cavitation experiments, the liquid composition is homogeneous and in mechanical equilibrium, so that supersaturation can be defined at the macroscopic scale. In the case of a liquid undergoing acoustic cavitation, microscopic bubbles are driven in radial oscillation by the sound field. For a sufficiently high excitation, the bubble motion is mainly driven by the inertia of the liquid: after a large expansion phase during the field depression phase, the bubble collapses very quickly on a time-scale of the order of several ns. This phenomenon is repeated at each acoustic period, for billions of bubbles. Under these conditions, it is easily conceivable that there is no mechanical equilibrium at the microscopic scale in the liquid, possibly resulting in local variations of supersaturation, but the precise phenomenon involved remains to be identified. Pressure gradients may be a possible source of mixture segregation. Following diffusion theory [START_REF] Hirschfelder | Molecular theory of gases and liquids[END_REF][START_REF] Bird | Transport phenomena[END_REF], when a mixture of two species is submitted to a pressure gradient, the lightest of the two is pushed toward low pressure regions. This forced diffusion process, known as pressure diffusion, remains generally weak but is commonly invoked to explain slow sedimentation in a still fluid mixture, or in ultracentrifuge applications where it is used profitably to separate species of large molecular weight from a solvent [START_REF] Archibald | The process of diffusion in a centrifugal field of force[END_REF]. Since the outward acceleration of the bubble at the end of its collapse is about twelve orders of magnitude higher than gravity, it is conceivable that the corresponding huge pressure gradient would segregate very efficiently the species present in the liquid.

To explore this issue, we have developed a perturbation method to solve the convection-diffusion mass transport equation of a species in a binary mixture around a single cavitation bubble, including the pressure diffusion effect (Louisnard & al., 2006). The transport equation reads

∂C A ∂t m:∇C A D∇:∇C A ∼ bC A ∇p 1
where C A is the mass-fraction of species A, v(r,t) and p(r,t) are respectively the radial velocity and pressure field around the oscillating bubble, and D is the diffusion coefficient of species A in the mixture. The para-Fig. 11: XRD patterns obtained from thermal hydrolysis at high supersaturation with and without ultrasound in the case of titanium dioxide.

Fig. 12: XRD patterns obtained from thermal hydrolysis at moderate and high supersaturation with ultrasound in the case of titanium dioxide. meter ∼ b is the segregation parameter linked to the difference between the apparent densities of the two species:

∼ b M A RT V A M A À V B M B ! : 2
where M i and V i are the molecular weight and partial molal volume of species i, respectively. It must be noted that the theory remains valid if one of the species, say A, is replaced by sufficiently small particles (typically at the nanoscopic scale):in this case V A =M A represents 1/ q A the inverse of the density of the particle material and the molecular weigth is calculated by M A = N am q A V A , where N am and V A is respectively the Avogado's number and the particle volume. The approximate analytical solution obtained yields the concentration of species A at the bubble wall in the following form:

C A t; 0 C A0 exp bI " 1 À b D xR 2 0 ! 1=2 Gt # 3
where C A0 is the mass-fraction of A in the undisturbed mixture, I is a positive constant, G(t) an oscillating function, and b is the dimensionless version of the segregation parameter ∼ b. Both quantities I and G(t) depend on the bubble dynamics and can be easily calculated once the time-dependant bubble radius is known. In the case of inertial cavitation, where the bubbles collapse violently at each acoustic period, G(t) is mainly a large positive short pulse repeating periodically at the end of the bubble collapse, lingering on a time-scale of several ns. In order to obtain a global picture of the segregation phenomenon, it is interesting to calculate the quantities:

X m expbI 4 for a given mixture around a typical cavitation bubble driven at increasing power levels. The quantity X m is the segregation rate at the bubble wall averaged over one bubble cycle and DX m is the peak value of the variation of the segregation rate, reached at the end of each bubble collapse. When b = 0, pressure diffusion does not act and we recover an unsegregated mixture. When b < 0 (A is in this case the heaviest species), X m is lower than 1, which means that on average species A is depleted at the bubble wall, and DX m is positive, so that the mixture is overconcentrated periodically at each bubble collapse.

As an example, we consider a mixture of water with spherical copper nanoparticles, with radii ranging from 1 nm to 10 nm. The diffusion coefficient D is calculated from Stokes-Einstein theory. The bubble is assumed to be filled of air, has a radius of 4 lm at rest, and is driven by an acoustic field of frequency 20 kHz, which amplitude ranges from 0.8 to 1.6 bar, which are typical values in cavitation experiments.

The results are displayed in Figure 16: the absissa of both graphs is the amplitude of the acoustic field driving the bubble. Figure 16 (a) represents the average depletion X m of particles at the bubble wall, while Figure 16 (b) represents the amplitude DX m of the periodic particles over-concentration. It can be seen on the left graph that the particles progressively disappear from the bubble wall, as either the amplitude or the particle size is increased. On the right graph, it can be seen that as the particle size is increased at fixed P, the over-concentra-tion DX m first increases, and then decreases again for very large molecules. Similarly, for fixed particle size, DX m crosses a maximum as P increases (visible on the 5 nm and 10 nm curves). Combining the results of the two graphs, the following qualitative conclusions can be drawn:

• for small particles or low forcing amplitude, X m ≈ 1, DX m ≈ 0: the mixture remains almost unsegregated (see filled circles symbols on Figure 16),

• for medium particles or medium forcing amplitude, X m < 1, DX m large: particles are slightly depleted at the bubble wall, but largely over-concentrated at each bubble collapse (see filled squares symbols on Figure 16). Note that this overconcentration may reach more than 2 orders of magnitude,

• for large particles or large forcing amplitude, X m ≈ 0, DX m ≈ 0: particles are almost completely depleted at the bubble wall, and the over-concentration is small again (see open circles symbols on Figure 16), so that the particles are constantly held far from the bubble in this case. We now turn to apply these results to homogenous nucleation of crystals. Following classical nucleation theory [START_REF] Kaschiev | Nucleation, Basic theory with applications[END_REF], the lowest energy expense for a firstorder phase transition is achieved by the progressive aggregation of solute molecules to form molecular clusters. These clusters have the density of the new phase and they must reach a critical size, referred as ''nucleus'' for the transition to occur. Therefore, the nucleation time is mainly the time required for the clusters to reach this critical size, and any microscopic effect enhancing the formation of the smallest clusters may reduce the nucleation time. Clusters may grow by aggregation of solute molecules one by one with existing smaller clusters, as schematically displayed in Figure 17, where we have labeled the solute molecules by C 1 and a cluster formed by n molecules by C n . Direct aggregation between two clusters C m and C n is generally neglected since they are much less numerous in solution than solute molecules. It should be noted that clusters do have physical reality and they have been evidenced for example by [START_REF] Mullin | Evidence of molecular cluster formation in supersaturated solutions of citric acid[END_REF] by gravity-driven sedimentation. This is in fact a pressure diffusion effect, due to the densities difference between the solution and the clusters. Therefore, on the basis of above-described results on segregation, it is expected that clusters could be efficiently segregated by a cavitation bubble, in the following way:

• Solute molecules and small clusters would be remain unsegregated,

• Medium clusters would be periodically over-concentrated near the bubble wall at each bubble collapse (up to 2 orders of magnitude). This would favour not only their attachment to solute molecules, but would also enhance the generally neglected direct aggregation between two clusters (Figure 18), thus considerably increasing the overall nucleation rate.

• Large clusters would be finally pushed far from the bubble by the average effect. Again, this favours nucleation kinetics, in virtue of Le Chatelier principle since the ultimate products of the aggregation reactions will be constantly removed from the bubble wall. In conclusion, a cavitation bubble, by the drastic segregation effect it produces in the surrounding liquid, would promote nucleation by acting as a cluster attachment reactor. Other possible microscopic mechanisms remain to be explored (direct dependance of supersaturation on pressure, solvent evaporation in the bubble), but our calculations show that segregation is a potential candidate for the macroscopically-observed nucleation enhancement. Furthermore, an experimental study based on microscopic planar laser induced-fluorescence is in progress to assess the physical reality of the segregation phenomenon around a single bubble [START_REF] Grossier | Mixture segregation by an inertial cavitation bubble[END_REF].

Conclusions

Results are presented for several series of experiments on the formation of crystals of BaSO 4 , K 2 SO 4 , TiO 2 and Sucrose, both with and without application of ultrasound. These experiments show, as have those of others, that ultrasound has an effect on nucleation resulting in finer more uniform crystals. There is as yet no well established explanation of how these effects are caused.

Here we present a new hypothesis based on the segregation in a liquid mixture by pressure gradients induced by cavitation bubbles. This theory seems to be in qualitative agreement with the phenomena found in crystallisation with ultrasound and an experimental programme based on the observation single cavitation bubbles is at present underway to make a quantitative test of the hypothesis.

Fig. 1 :

 1 Fig.1: Experimental equipment for continuous precipitation of barium sulfate[START_REF] Gatumel | Nucleation control in precipitation processes by ultrasound[END_REF].

Fig. 2 :Fig. 3 :

 23 Fig.2: Supersaturation ratio S at the exit of the growth vessel as function of the supersaturation ratio at the exit of the nucleation cell, with or without ultrasound[START_REF] Gatumel | Nucleation control in precipitation processes by ultrasound[END_REF].

Fig. 4 :Fig. 5 :

 45 Fig. 4: Surface size distribution of BaSO 4 crystals at the exit of the nucleation cell (Gatumel et al., 1999).

Fig. 6 :Fig. 7 :Fig. 8 :

 678 Fig. 6: Induction time as function of absolute supersaturation for potassium sulfate at 15 °C assuming an activity coefficient ratio close to 1 (Lyczko et al., 2002).

Fig. 9 :

 9 Fig. 9: SEM photographs of titanium dioxide crystals at moderate supersaturation without and with ultrasound.

Fig. 10 :

 10 Fig. 10: XRD patterns obtained from thermal hydrolysis for titanium dioxide at moderate supersaturation with and without ultrasound.

Fig. 13 :

 13 Fig.13: d[START_REF] Gatumel | Nucleation control in precipitation processes by ultrasound[END_REF][START_REF] Bird | Transport phenomena[END_REF] at the exit of the sucrose crystallizer and DC E-S = C e -C s as function of the ultrasound power.

Fig. 14 :

 14 Fig. 14: Semi-logarithmic plot of the population density according to the size of crystals for different ultrasound power.

Fig. 15 :

 15 Fig. 15: Concentration and temperature as function of time at 25 W in the case of the Sucrose crystallization.

Fig. 16 :

 16 Fig. 16: Average (a) and periodic (b) concentration fields in a mixture of water with copper nanoparticles around a 4 lm air bubble driven by a 26.5 kHz acoustic field of various amplitudes P.

Fig. 17 :

 17 Fig. 17: Cluster growth by aggregation of a solute molecule C 1 on an existing cluster C n , following classical nucleation theory. Direct aggregation between clusters C m and C n has low probability to occur because of the low number of clusters compared to solute molecules.

Fig. 18 :

 18 Fig. 18: Cluster growth in presence of a cavitation bubble: pressure diffusion concentrates medium clusters near the wall of the bubble at each collapse, favouring aggregation between these clusters and a solute molecule, and possibly giving direct aggregation between two clusters. Large clusters are then held far from the bubble by the average effect.

Table 1 :

 1 Operating conditions for Sucrose crystallisation.

	Puissance	C E (%)	C S (%)	C E -C S (%)	DC en sortie
	(W)				(%)
	70	75.72	74.41 ± 0.20	1.31	3.91
	50	75.80	74.49 ± 0.17	1.31	3.99
	35	75.72	74.71 ± 0.19	0.98	4.24
	30	75.72	74.49 ± 0.18	1.23	3.99
	25	76.25	74.30 ± 0.21	1.96	3.8

of them has been demonstrated clearly either by modeling or by experiments

[START_REF] Ruecroft | Sonocrystallization: the use of ultraound for improved industrial crystallization[END_REF]

. In the following we present an alternative hypothesis based on the molecular segregation of a liquid mixture by cavitation bubbles.

Table 2 :

 2 Calculated values of G0, n0 and b for different ultrasound power in the case of Sucrose

	GL G 0	"	1	G 0 s L	# b
	Ultrsasound		s(s)	G 0 (m/s)	n 0 (nb/m -4 )	b	J (nbre/m 3 /s)
	power (W)
	30			502	9.6*10 -10	4.9*10 20	0.67	4.7 10 11
	35			509	7.6*10 -10	22.7*10 20	0.59	17.3 10 11
	50			570	1.3*10 -10	39.0*10 20	0.81	5.1 10 11
	70			502	7.2 10 -10	80.0 10 20	0.59	57.6 10 11