Raza Ur-Rehman
email: ur-rehman@irccyn.ec-nantes.fr

Stéphane Caro
email: caro@irccyn.ec-nantes.fr

Damien Chablat
email: chablat@irccyn.ec-nantes.fr

Philippe Wenger
email: wenger@irccyn.ec-nantes.fr

Multiobjective Path Placement Optimization of Parallel Kinematics Machines Based on Energy Consumption, Shaking Forces and Maximum Actuators Torques: Application to the Orthoglide

Keywords: Multiobjective Optimization, Path Planning, Parallel Kinematics Machine

This paper deals with the multi-objective path placement optimization for Parallel Kinematics Machines (PKMs) based on energy consumption, actuators torques and shaking forces. It aims at determining the optimal location of a given test path within the workspace of a PKM in order to minimize the electric energy used by the actuators, their maximal torque and the shaking forces subject to the kinematic, dynamic and geometric constraints. The proposed methodology is applied to the Orthoglide, a three-degree-of-freedom translational PKM, as an illustrative example.

Introduction

Optimal trajectory planning has been a relevant area for roboticists for many years. Several authors have worked on trajectory planning based on different optimization objectives. A review of trajectory planning techniques is given in [START_REF] Ata | Optimal trajectory planning of manipulators: A review[END_REF]. Trajectory planning usually aims at minimizing the travel distance [START_REF] Tian | Motion planning for redundant manipulators using a floating point genetic algorithm[END_REF], travel or machining time [START_REF] Chan | Genetic-based minimum time trajectory planning of articulated manipulators with torque constraints[END_REF][START_REF] Pateloup | Corner Optimization for pocket machining[END_REF] and/or the energy consumed [START_REF] Field | Iterative dynamic programming, an approach to minimum energy trajectory planning for robotic manipulators[END_REF][START_REF] Hirakawa | Trajectory planning of redundant manipulators for minimum energy consumption without matrix inversion[END_REF], while satisfying several geometric, kinematic and dynamic constraints. The trajectory planning deals with the determination of the path and velocity/acceleration profiles (or the time history of the robot's joints), the start and end points of the trajectory being predefined and fixed in the workspace.

Another less explored aspect of trajectory planning is the placement of a given path within the workspace. It aims at determining the optimum location of a predefined path to be followed by the end-effector of a Parallel Kinematics Machine (PKM) or a robot, within its workspace with respect to one or many given objective(s) and constraint(s). This path can be the shape of a component to be machined, a welded profile or an artistic/decorative profile etc. In such situations, the trajectory planner cannot alter the shape of the path but he/she can only play with the location of that path within the workspace in order to optimize one or several criterion(a). Such an approach can be very interesting in many robotic applications. For example, in machining applications, the location of the workpiece within the workspace may affect the electric energy used by its actuators.

The path placement problem has not been extensively studied in the past. Nevertheless, some researchers proposed to solve it with respect to various optimization objectives. Several performance criteria for path location problems can be considered simultaneously (multiobjective) or individually, such as travel time, different kinetostatic performance indices (such as, manipulability or the conditioning number of the normalized kinematic Jacobian matrix), kinematic performance (velocity, acceleration), collisions, wear and vibration reduction, energy consumption etc.

Nelson and Donath [START_REF] Nelson | Optimizing the location of assembly tasks in a manipulators workspace[END_REF] proposed an algorithm for the optimum location of an assembly task R. Ur-Rehman, S. Caro, D. Chablat and P. Wenger, submitted to Mech. and Mach. Theory 3 in a manipulator workspace while taking the manipulability measure as the optimization criterion. The location of the assembly task that results in the highest manipulability is considered as optimal position. Aspragathos [START_REF] Aspragathos | Optimal location of path following tasks in the workspace of a manipulator using genetic algorithm[END_REF][START_REF] Aspragathos | Optimal location of a robot path when considering velocity performance[END_REF] used the concept of the orientation of the manipulability ellipsoid relative to the desired path and introduced a criterion to characterize the best velocity performance of the robot end-effector with the path location. Fardanesh et al. [START_REF] Fardanesh | Minimum cycle time location of a task in the workspace of a robot arm[END_REF] proposed an approach for optimal positioning of a prescribed task in the workspace of a 2R-robot arm to minimize the cycle time. In another study, Feddema [START_REF] Feddema | Kinematically optimal robot placement for minimum time coordinated motion[END_REF] formulated and solved a problem of robot base placement for a minimum joint motion time within a work cell. The proposed algorithm considers only the kinematics and the maximum acceleration of each joint in order to obtain a 25% cycle time improvement for a typical example. Hemmerle [START_REF] Hemmerle | Optimal path placement for kinematically redundant manipulators[END_REF] presented an algorithm for optimum path placement of a redundant manipulator by defining a cost function related to robot joints motion and limits. The proposed approach did not consider the path as a whole but discrete points of that path. Pamanes and Zeghloul [START_REF] Pamanes | Optimal placement of robotic manipulators using multiple kinematic criteria[END_REF] considered several kinematic indices to find the optimal placement of a robot by specifying the path with a number of points and then assigning an optimization criterion to each point. The objective was to find the path location in order to have optimal values of all the criteria assigned to the path points.

In [START_REF] Pamanes | On the optimal placement and task compatibility of manipulators[END_REF], the problem of optimal placement with joint-limits and obstacle avoidance is addressed.

Lately, a general formulation was presented to determine the optimal location of a path for a redundant robotic manipulator while dealing with mono-and multi-objective problems [START_REF] Pamanes | Single and multiobjective optimization of path placement for redundant robotic manipulators[END_REF].

In another relevant work, Wang et al. [START_REF] Wang | A study on workspace, boundary workspace analysis and workpiece positioning for parallel machine tools[END_REF] proposed a mathematical model to determine the right position and orientation of a workpiece in the workspace of a Stewart platform based NC-miller machine tool. The idea of their study was to analyze the effects of geometric and nongeometric constraints on the workspace and subsequently, using these information to determine the right pose of the workpiece in the workspace a parallel machine tool. However, the only objective that they considered was the location of the workpiece for the accessibility of a specific machine tool.

With a general literature survey, it comes out that although several performance indices are R. Ur-Rehman, S. Caro, D. Chablat and P. Wenger, submitted to Mech. and Mach. Theory 4 introduced or considered, there is no emphasis on the dynamic aspects reflecting the energy consumption by the PKM actuators, while executing a specific task. Another less explored but important criterion of optimal path placement can be the minimization of the shaking forces and moments experienced by the base of the PKM. Shaking forces and moments can affect the performance of a PKM in terms of excessive loads, accuracy, wear, fatigue, etc. Accordingly, we introduce two indices characterizing the variations and the maximum value of the shaking forces.

The maximum actuators torque is considered as another optimization objective. Keeping in view these aspects, the authors recently presented an approach for the optimal path placement for a manipulator based on energy consumption [START_REF] Ur-Rehman | Path placement optimization of manipulators based on energy consumption: application to the orthoglide 3-axis[END_REF]. In the same vein, it is pertinent to follow a multiobjective approach to optimize the location of a path to be followed by the end-effector of a PKM within its workspace. The objectives are to minimize the energy consumed by its actuators, minimize the shaking forces and/or moments and minimize the maximum actuators torques. Finally, the proposed approach is applied to the Orthoglide: a three-degree-of-freedom translational PKM.

The paper is organized as follows: in section 2, a multi-objective path placement optimization problem is formulated, by means of the design variables, the objective functions and the constraints. Moreover, an algorithm is proposed to solve such a problem. In section 3, the optimum locations of a test path within the workspace of the Orthoglide are highlighted.

Path Placement Optimization

The problem aims at determining the optimal location of a predefined path to optimize some objective functions. The entire path is supposed to be known within the framework of this research work. The path location can be defined in a similar way as to define the location of a workpiece with respect to a PKM reference point. The optimization problem is composed of three sets, namely, the set of design variables, the set of objective functions and the set of design constraints. Accordingly, the optimization problem aims at determining the design variables, characterizing the path location, in order to minimize the objective functions subject R. Ur-Rehman, S. Caro, D. Chablat and P. Wenger, submitted to Mech. and Mach. Theory 5 to the design constraints.

Decision Variables: Path Location Parameters

In order to formulate and describe the problem, two frames are defined: i) the path frame F p and ii) the base frame F b , as shown in Fig. 1(a). The path frame F p , is attached to the given/required path at a suitable point such as the geometric center of the path. As F p is attached to the path, the end-effector trajectory parameters remain constant in this reference frame, no matter where it is located, i.e., the path is fully defined and constant in F p . It can also be named as workpiece frame since it characterizes the position and the orientation of the workpiece within the workspace. The base frame F b can also be called global or PKM frame. It is attached to the PKM base and is used to locate a workpiece (or F p ) with respect to the PKM coordinate system. The location and orientation of F p with respect to F b can be defined in such a way that the whole path lies within the workspace. The position of F p

         (2)
Let x = [x Op y Op z Op φ θ ψ] T define the path location within the workspace in F b .

Then, the components of x are the decision variables of the optimization problem.

In the context of a general machining process like milling operation, the feature to be machined in the workpiece is defined with respect to its frame, namely F p . Likewise, the machining operation conditions such as machining velocity and acceleration are fully defined in F p . Finally the part location in the workspace is defined in F b .

Optimization Objectives

The path placement optimization is performed with the aim of minimizing four objective functions, namely, (i) electric energy consumption, (ii) variations in the shaking forces, (iii) maximum shaking forces and (iv) maximum actuators torques. In the following sections, the mathematical expressions of these objective functions are given.

The energy used by the motors depends on their velocity and torque. As the electric current drawn varies with motors required velocities and torques, the motor's self-inductance phenomenon appears. The current I drawn by the motors and the motor electromotive potential V e can be calculated as a function of the required torque τ and the angular velocity ω of the actuators, namely,

I = τ /K t (3)
V e = K e ω (4)
K t being the torque sensitivity factor expressed in [Nm/A] and K e the back electromotive force constant expressed in [V.(rad/sec) -1 ]. The total electric power P T is composed of resistive power loss, P J ; inductive power loss, P L ; and the power used to produce electromotive force, P EM [START_REF] Lacroux | Les actionneurs lectriques pour la robotique et les asservissements[END_REF], i.e.,

P T = P J + P L + P EM (5)
where

P J = RI 2 , P L = LI dI dt , P EM = V e I (6)
R being the motor winding resistance expressed in Ohm [Ω] and L the motor inductance coefficient expressed in Henry [H]. Finally, the energy E consumed by a motor can be evaluated by integrating P T over the total trajectory time T, namely,

E = T 0 P T dt (7)
P T being the instantaneous electric power at instantaneous time t, defined in Eq. ( 5).

It should be noted that Eq. ( 3) allows us to consider the energy used by the actuators while they do not move but still produce a torque to keep the PKM at a certain stationary configuration (with respect to that particular direction or actuator). Finally, the total energy E t consumed by n actuators can be written as:

E t = n i=1 E i (8)
the first objective of the path placement optimization problem can be written as:

f 1 (x) = E t → min (9)

Shaking Forces

Shaking or dynamic forces and moments are the inertial forces and moments exerted on the base of a PKM due to the uneven mass distribution. These forces/moments may deteriorate system performance by introducing excessive vibration, noise, wear and fatigue. Furthermore required input torques and forces may also increase to account for these shaking effects. Hence, in order to overcome these drawbacks and to improve system performance in terms of accuracy, precision, fatigue life, vibration reduction, motion planning and control, the study of the shaking forces and moments is of prime importance.

Dynamic balancing has been an area of research for some decades and several authors have contributed to this domain by formulating and analyzing the problems for either some particular applications or for more general prospects [START_REF] Chiou | Optimum balancing design of the drag-link drive of mechanical presses for precision cutting[END_REF][START_REF] Wu | Synthesis of reactionless spatial 3-dof and 6-dof mechanisms without separate counter-rotations[END_REF][START_REF] Moore | Determination of the complete set of shaking force and shaking moment balanced planar four-bar linkages[END_REF]. Here, we will use this concept in order to optimize the path location within the workspace of the PKM to minimize the shaking forces effects.

For the purpose of simplicity, shaking forces are exclusively considered in this research work while shaking moments are left for future works.

In order to assess the effect of shaking forces, two indices are proposed, namely, 1. Shaking force variation index, I δf , considers the maximum variation in the shaking forces along the trajectory, i.e.,

I δf = max (dF sh ) (10)
where dF sh = n i=1 dF i sh is the sum of the variations in the n shaking forces, n being the number of limbs. The variation in the shaking forces for the i th limb is given by,

dF i sh = [max(F i x ) -min(F i x )] 2 + max(F i y ) -min(F i y ) 2 + [max(F i z ) -min(F i z )] 2 (11)
where F i x , F i y and F i z are 1×N t dimensional vectors, representing the respective force profile along trajectory, N t being the number of time steps. The variations in the shaking forces will influence the base platform vibrations. Hence, the second objective function of the path optimization problem can be written as,

f 2 (x) = I δf → min (12)
2. Maximum shaking force index, I f , characterizes the magnitude of the maximum shaking forces experienced by the base platform along the trajectory, i.e.,

I f = max (F sh ) ( 13
)
where F sh = n i=1 F i sh is 1 × N t dimensional vector representing the sum of the shaking forces experienced by the n limbs with N t time steps. The maximum shaking force gives an idea of the maximum extra loads experienced by the actuators resulting from the shaking forces, for a particular trajectory. Hence the knowledge of the maximum shaking forces can help the designer to select the appropriate actuators to cater for these excessive loads on the system. Accordingly, the third objective function of the path optimization problem can be written as,

f 3 (x) = I f → min (14)

Maximum Torque

In order to reduce the actuators loads, the magnitude of the maximum torque τ max experienced by the PKM actuators is considered. For a PKM of n actuators, τ max is defined as,

τ max = max (τ 1max , . . . , τ imax , . . . , τ nmax ) (15)
τ imax being the magnitude of the maximum torque experienced by the i th actuator along the trajectory. Consequently, the fourth objective function of the path optimization problem can be written as,

f 4 (x) = τ max → min (16)

Constraints

The path placement optimization problem is subject to geometric, kinematic and dynamic constraints. Geometric constraints include joint limits and the boundaries of the workspace.

Kinematic constraints deal with the maximum actuators velocities whereas dynamic constraints are related to actuators wrenches, namely,

q il ≤ q i ≤ q iu , qi ≤ qiu , τ i ≤ τ iu (17)
q i , qi , τ i are respectively the i th actuator displacement, rate and torque. q il is the lower bound and q iu (resp. qiu and τ iu ) is the upper bound of i th actuator displacement (resp. rate and torque). For a given path placement vector x, these constraints can be evaluated by means of the PKM kinematic, velocity and dynamic models.

It is noteworthy that the PKM geometric constraints guarantee that the whole path lies inside the prescribed workspace. Similarly, the bounds on actuators rate ( qiu ) and torque (τ iu ) ensure that the PKM will not go through any singular configuration while following the path.

Problem Statement and Solution

The goal of this research work is to help the path planner find the best location of the path to be followed by the PKM in order to minimize the four objective functions defined in the foregoing subsections. The multi-objective path placement optimization problem can be formulated as:

"For a predefined path in F p , find the optimum location and orientation of F p with respect to F b , defined by the decision variables x, in order to minimize the objective functions f 1 , f 2 , f 3 and f 4 while respecting the kinematic, velocity and dynamic constraints of the PKM "

Mathematically, the problem can be formulated as follows:

min x (f 1 (x), f 2 (x), f 3 (x), f 4 (x)) subject to:            q il ≤ q i ≤ q iu qi ≤ qiu τ i ≤ τ iu (18)
R. Ur-Rehman, S. Caro, D. Chablat and P. Wenger, submitted to Mech. and Mach. Theory11

To solve the problem, a general optimization approach based on Multi-Objective Genetic Algorithm (MOGA) is proposed as illustrated in Fig. 3 Case Study: Application to the Orthoglide

Description of the Orthoglide

The Orthoglide is a Delta-type PKM [START_REF] Clavel | Delta, a fast robot with parallel geometry[END_REF] dedicated to 3-axis rapid machining applications developed in IRCCyN [START_REF] Wenger | Kinematic analysis of a new parallel machine tool: The Orthoglide[END_REF]. It gathers the advantages of both serial and parallel kinematic architectures such as regular workspace, homogeneous performances, good dynamic performances and stiffness. The Orthoglide is composed of three identical legs, as shown in Fig. 3 In the scope of this study, L workspace is equal to 0.200m. Accordingly, the coordinates of points Q + , Q -and C for the said workspace are given in Table 1.

Similarly, the prismatic actuators bounds, ρ min and ρ max , can be calculated [START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF]. Table 2 shows the lower and upper bounds of the prismatic joints displacements and their maximum allowable velocity and torque for the Orthoglide. The geometric, kinematic and dynamic parameters of the Orthoglide are defined in [START_REF] Wenger | Kinematic analysis of a new parallel machine tool: The Orthoglide[END_REF][START_REF] Pashkevich | Stiffness analysis of overconstrained parallel manipulators[END_REF][START_REF] Chablat | Architecture Optimization, of a 3-DOF Parallel Mechanism for Machining Applications, The Orthoglide[END_REF][START_REF] Guégan | Identification of the Dynamic Parameters of the Orthoglide[END_REF]. The electric energy E i used by each actuator is calculated by means of Eqs. ( 3)- [START_REF] Nelson | Optimizing the location of assembly tasks in a manipulators workspace[END_REF]. As the Orthoglide has three 3-phase Sanyo Denki synchronous servo motors (ref erence : P 30B0604D), Eq. ( 5) is multiplied by 3 to cater for the power consumed by the each phase of the motor in order to calculate the electric power P T i used by each actuator, i.e.,

P T i = 3(RI 2 + LI dI dt + V e I) (19)

Shaking forces formulation for the Orthoglide

Shaking forces f sh at the base platform of a PKM depend on the mass and the acceleration of the center of mass of each moving element. For a system of w masses, f sh can be expressed as:

f sh = w j=1 m j cj ( 20
)
where m j is the mass and cj is the acceleration of the center of mass of the j th element. In order to calculate the shaking forces at the Orthoglide base-frame, three reference frames, corresponding to each leg, are defined at points A 1 , A 2 and A 3 , as shown in Fig. 4. The shaking forces for each leg are calculated independently in the respective reference frame. Each leg is supposed to be composed of six components, namely, Besides these leg elements, there are two other moving masses: the mass of the moving part of the prismatic actuator, M a , and the one of the end-effector, M p . The mass of each member M j is denoted by m j and the center of mass with respect to the point A i (i = 1, 2, 3) of each element 0.932 kg is denoted by C ji , (j = 1 . . . 6, a, p), as shown in Fig. 5. The mass of each element and the geometric parameters of the Orthoglide are presented in Table 3. Assuming that the material of the leg components is homogeneous with center of mass at their geometric center, vectors c ji of the j th element of the i th leg can be expressed in terms of the geometric (L b , L f , d, e, λ) and the configuration (q 1i , q 2i , q 3i ) parameters of the Orthoglide, as shown in Fig. 5. Therefore, the acceleration cji of each element can be calculated. These relations are given in the Appendix.

M 1 : foot of length L f ,
F 0 X 0 Y 0 Z 0 A 1 X 1 Y 1 Z 1 A 2 X 2 Y 2 Z 2 A 3 X 3 Y 3 Z 3 F b O b X b Y b Z b
Knowing the mass and acceleration of each element of the i th leg, the shaking forces at the base point of the leg (point A i ) can be calculated by using Eq. ( 20), i.e.,

f i sh = (m 1 + m a ) c1i + m 2 c2i + m 3 c3i + m 4 c4i + m 5 c5i + m 6 c6i + m p cpi (21)
The total shaking force at the base frame of the Orthoglide can be obtained by summing up the forces experienced at points A 1 , A 2 and A 3 . With the definitions of reference frames at points A i , the total shaking force in the reference frame F 0 , of origin point A 1 , as shown in Fig. 4, can be written as:

A i A i X i X i Y i Y i Z i Z i P P M 1 M 2 M 3 M 4 M 5 M 6 C ai C ai C 1i C 1i C 2i C 2i C 3i C 3i C 4i C 4i C 5i C 5i C 6i C 6i C pi C pi d λ L f L b L b c o s q3 i c o s q2 i L b sin q 3i L b cos q 3i sin q 2i q 1i q 1i q 2i q 2i q 3i q 3i
f 0 sh =       F 1 x + F 2 z + F 3 y F 1 y + F 2 x + F 3 z F 1 z + F 2 y + F 3 x       (22)
The magnitude of the shaking force experienced by the i th actuator at point A i , can be written as:

F i sh = (F i x ) 2 + F i y 2 + (F i z ) 2 (23)
Accordingly, two shaking force indices, defined in Sec. 2.2.2, can be calculated for the given discrete time steps.

Maximum Torque for the Orthoglide

The magnitude of the maximum torque τ max experienced by the Orthoglide actuators can be written as:

τ max = max (τ 1max , τ 2max , τ 3max ) ( 24
)
with τ imax = max(τ i ), τ i being the 1×N t dimensional vector of the i th actuator for N t trajectory points. The inverse dynamic model of the Orthoglide described in [START_REF] Guegan | Dynamic modeling of the Orthoglide[END_REF][START_REF] Guégan | Identification of the Dynamic Parameters of the Orthoglide[END_REF] is used to evaluate the actuators torques.

Trajectory Planning and External Forces

In order to apply the methodology proposed for path placement optimization, a rectangular test path is proposed. The test path is defined by the length L and the width W of the rectangle, as shown in Fig. 6 are considered to be always parallel. Accordingly, there are four path placement variables, i.e.,

x = [x Op y Op z Op φ] T , as illustrated in Fig. 6(a).

The magnitude of the end-effector velocity is supposed to be constant along the path. Hence, [ ẋP p ẏP p żP p ] T in the path frame can be evaluated as a function of time. Figure 7 shows the position and velocity profiles in F p for a 0.05 m×0.10 m rectangular path and for a constant end-effector velocity of 1.0 m.s -1 . Position and velocity vectors defined in F p can be expressed in F b by means of the transformation matrix defined in Eq. ( 2), namely,

Ob X b Y b p F b F b X p Y p O p P p Fp
         x P b y P b z P b 1          =          Cφ -Sφ 0 x Op Sφ Cφ 0 y Op 0 0 1 z Op 0 0 0 1                   x Pp y Pp z Pp 1          ,          ẋP b ẏP b żP b 1          =          Cφ -Sφ 0 0 Sφ Cφ 0 0 0 0 1 0 0 0 0 1                   ẋPp ẏPp żPp 1         
where C and S respectively stands for cosine and sine functions. For the matter of simplicity and not to deal with tangent and curvature discontinuities, the path is considered to be composed of four independent line segments. Therefore the discontinuities between the segments are neglected.

In order to analyze the effect of external cutting/machining forces in the generation of a given R. Ur-Rehman, S. Caro, D. Chablat and P. Wenger, submitted to Mech. and Mach. Theory19

X b Y b Z b F a F r F f V p

Results and Discussions

Multiobjective optimum path placement problems are solved for the Orthoglide for the test path described in the previous section. The constraints of the optimization problem are the geometric, kinematic and dynamic ones. The geometric constraints are the upper and lower bounds of the prismatic joints variables. The kinematic constraints are the maximum velocities that the prismatic actuators can produce whereas the dynamic constraints are the limits on the actuators torques. Kinematic and dynamic constraints are specified by the actuators designer and are obtained from the catalogue. As already mentioned, the optimization variables are the coordinates of the origin of path reference frame F p and the orientation angle of F p with respect to F b . The numerical values of the constraints used for the Orthoglide are given in Table 2.

Multi-objective path placement optimization problem for the Orthoglide can be expressed as,

min x (f 1 : E t , f 2 : I δf , f 3 : I f , f 4 : τ max )
subject to:

           ρ min ≤ ρ x,y,z ≤ ρ max v x,y,z ≤ v max τ x,y,z ≤ τ max (25)
where x = [x Op y Op z Op φ] T . The subscripts x, y and z are used for three prismatic actuators or three Cartesian directions.

The problem is modeled with MATLAB and modeFRONTIER [29]. A multiobjective ge- Table 5 shows that the variations in the shaking forces (I δf ) can be significantly reduced (almost to zero) with an appropriate path placement. As a matter of fact, up to 60% of the energy consumption can be saved with a proper path placement. Maximum shaking forces I f and maximum actutators torque τ max can be reduced to 17% and 25%, respectively.

Figure 12 shows the shaking forces experienced by the three prismatic actuators of Orthoglide for minimum energy consumption and minimum I δf path locations. It can be seen that although the maximum values of the shaking forces for both Pareto-points are almost the same, their variations can be reduced considerably with proper path location. The smoother the shaking force variations, the lower the vibrations in the mechanism. From Table 5 and Fig. 11(a) it can be noted that optimum points with respect to E t , I δf and I f lie in the neighbourhood of the isotropic configuration of the Orthoglide (x p = y p = z p ≈ 0) whereas for τ max , z p attain their minimum value, i.e., at the base of the workspace. Similarly with respect to orientation, φ, E t is minimum for φ ≈ 0 • and the other three objectives are minimum for φ ≈ 45 • .

Figure 13 summarizes the variational trends as well as the inter-dependency between the objective functions and design variables by means of a scatter matrix. The lower triangular part of The correlation coefficients vary from -1 to 1. Two variables are strongly dependent when their correlation coefficient is close to -1 or 1 and independent when the latter is null. From Fig. 13,

X b [m] Y b [m] Z b [m] τ min I f min E min δI f min Y b Z b O b Q - Q + -0.
X b [m] Y b [m] Z b [m] δI f max E max I f max τ max Q + Q - O b Z b Y b X b -0.
F sh [N] Shaking forces for minimum E t Time [sec] F sh [N]
Shaking forces for minimum I δf • the correlation of I f with E t , I δf and τ max is very low as the corresponding correlation coefficients are equal to -0.113, 0.185 and -0.065, respectively.

F 1 F 2 F 3 F 1 F 2 F 3 0
• φ has strong and direct correlation with E t (0.907) whereas it has strong and inverse correlation with I δf and τ max (-0.966, -0.828, respectively);

• x p , y p , z p have very week or unpredictable relations with respect to all objectives and parameters.

As I δf and τ max are linearly related, the lower I δf , the lower τ max , no matter the scale of variation of both functions. E t and τ max are antagonistic. Likewise, E t and I δf are antagonistic.

Regarding the design variables, the path orientation φ is the most influential for the considered test path. Finally, the foregoing results are only valid for the Orthoglide and the given test path. However, the methodology illustrated in this section is appropriate for any PKM and test path. The proposed methodology is applied to the Orthoglide, a 3-DOF translating PKM with a regular cubic workspace. A rectangular shaped test path is considered as an illustrative example.

Such a path can be used to realize pocketing operations.

The use of the electric energy instead of the mechanical energy as an optimization criterion is pertinent. Although actuator electric energy consumption depends on the mechanical energy requirements, the electric energy evaluation is more comprehensive than its mechanical counterpart. The general approach used to calculate the mechanical energy with the help of velocity and dynamic models, i.e., by using actuators torques and velocities, may lead to an under estimation of the energy requirements in the case where actuators are experiencing torques with zero velocities. Besides, usual mechanical energy calculations do not consider the resistive energy loss in the motors windings as well as the energy loss due to the variations in the actuators velocities. Those variations affect the current requirements and hence induce electromotive forces in the actuators. Accordingly, the electric energy formulation takes into account all these energy losses.

The energy required to perform a given task depends on the position and the orientation of the R. Ur-Rehman, S. Caro, D. Chablat and P. Wenger, submitted to Mech. and Mach. Theory26 task within the workspace. Accordingly, some electric energy can be saved by properly selecting the position and the orientation of the task. Indeed, a misplaced task can cause excessive energy consumption and can force the actuators to go over their performance limits. For the Orthoglide, the optimum path location is found to be in the neighbourhood of the isotropic configuration but there is no general rule to predict the exact optimal position and orientation of a task particularly for a complicated three dimensional task or for an irregular workspace.

However, a detailed analysis of the energy variation within the workspace for a given task can lead to the optimal position/orientation of that particular task. Numerical optimization algorithms are useful for such a comprehensive analysis in which all the problem constraints and performance measures can be considered simultaneously.

Variations in the shaking forces experienced by the base platform of a PKM also vary with the path location. However, the lower the shaking forces, the higher the energy consumption.

The magnitude of the maximum shaking forces can also be reduced with path location but to a lower extent compared to the reduction of the variations in the shaking forces. Similarly, the magnitude of the maximum actuators torques can be reduced with appropriate path location.

However, the energy consumption would increase with the minimization of the shaking forces effects and maximum torques. Hence, a trade-off has to be made by the user.

In the future work, some multi-objective problems will be treated with more objectives and constraints and with complex tasks and irregular workspaces. Furthermore, along with shaking forces, shaking moments will also be considered.

Figure 1 :Fp ( 1 )sin θ cos θ sin ψ cos θ cos ψ

 11 Figure 1: Path placement characterization, F b and F p being the base and path frames

2 .

 2 Starting with the input variables, the inverse kinematic model (IKM), the inverse velocity model (IVM) and the inverse dynamic model (IDM) of the PKM are used to determine the objectives functions and constraints for particular designs variables. MOGA is used to generate the initial population, evaluate the objective functions and constraints in order to generate new populations by carrying out the reproduction, cross over and mutation operations. The Pareto optimal solutions are obtained from the final feasible population.

Figure 2 :

 2 Figure 2: Flowchart of the path location optimization process

Figure 3 :

 3 Figure 3: A snap and workspace of the Orthoglide

M 2 :M 6 :

 26 the small side of the parallelogram joint attached to the foot. Its length is equal to d, M 3 , M 4 : the longer sides of the parallelogram joint. Their length is equal to L b , M 5 : the small side of the parallelogram joint attached to the end-effector, the link between M 5 and the end-effector of length e.

Figure 4 :

 4 Figure 4: Definition of the reference frames

Figure 5 :

 5 Figure 5: Orthoglide i th leg

  (a). Path reference frame F p is located at the geometric center of the rectangle. This type of path can be the example of the generation of a rectangular pocket like that of Fig. 6(b). The position of F p in the base frame F b is defined with the Cartesian coordinates of the origin of F p , O p (x Op , y Op , z Op ) and the orientation of F p with respect to F b is given by Euler's angles, as depicted in Fig. 1(b). For the sake of simplicity, only one of the three rotation angles is considered i.e, rotation about Z b -axis while X b Y b and X p Y p planes R. Ur-Rehman, S. Caro, D. Chablat and P. Wenger, submitted to Mech. and Mach. Theory18

Figure 6 :

 6 Figure 6: Test path characterization

Figure 7 :

 7 Figure 7: Test trajectory for a rectangular path of size 0.05 × 0.10 m

Figure 8 :

 8 Figure 8: Cutting forces

  Figure 9: modeFRONTIER model

Figure 10 :

 10 Figure 10: Pareto frontier for the Orthoglide path placement problem

Figure 11 :

 11 Figure 11: Path locations in the Orthoglide workspace for minimum and maximum objective functions

Figure 12 :

 12 Figure 12: Shaking forces experienced by three actuators for E min and I δf min

Figure 13 :

 13 Figure 13: Scatter matrix for objective functions and design variables

Table 2 :

 2 Orthoglide

		actuators parameters (i = x, y, z)
	ρ i min	0.126 m
	ρ imax	0.383 m
	v imax	1.00 m.s -1
	τ imax	1.274 Nm

Table 3

 3

		: Orthoglide parameters
	Parameter	Value	Parameter Value
	m a	0.300 kg	L b	0.310 m
	m 1	0.248 kg	d	0.080 m
	m 2 , m 5	0.095 kg	L f	0.150 m
	m 3 , m 4	0.117 kg	e	0.031 m
	m 6	0.010 kg	λ	45 •
	m p

Table 4 :

 4 MOGA) is used to obtain the Pareto frontiers for a rectangular test path of R. Ur-Rehman, S. Caro, D. Chablat and P. Wenger, submitted to Mech. and Mach. Theory21

	modeFRONTIER algorithm parameters
	Scheduler	MOGA-II
	Number of iterations	100
	Directional cross-over probability 0.5
	Selection probability	0.05
	Mutation probability	0.1
	DNA string mutation ratio	0.05
	DOE algorithm	Sobol
	DOE number of designs	40
	Total number of iterations	40 × 100 = 4000
	netic algorithm (

Table 5 :

 5 Design parameters that correspond to the Pareto solutions for which the objective functions are either a minimum or a maximum Vs I f Vs τ max

	Objective	x p [m]	y p [m]	z p [m] φ [deg] Value	∆	%∆
	E t [J]	E max -0.0276 -0.0348 -0.0789 E min 0.0141 -0.0044 0.0097	37.2 1.44	56.49 33.54 59.38 22.94
	I δf [N]	I δf max -0.0204 -0.0693 -0.0019 I δf min 0.0009 -0.0011 -0.0164	0.63 45.0	2.445 0.034	2.41 98.60
	I f [N]	I f max I f min	0.0389 -0.0801 -0.1114 0.0009 -0.0011 -0.0164	40.1 45.0	2.620 0.456 17.39 2.164
	τ max τ max [Nm] τ min	0.0381 -0.0056 0.0714 0.0053 0.0069 -0.1166	24.0 42.0	1.135 0.849	0.287 25.24

  Ur-Rehman, S. Caro, D. Chablat and P. Wenger, submitted to Mech. and Mach. Theory24 E t , the higher τ max and I δf .

• E t , I δf and τ max strongly dependent as the correlation coefficients between E t , I δf ; E t , τ max ; and I δf , τ max are equal to -0.822, -0.785 and 0.819, respectively. However, the lower R.

Acknowledgements

The support of the Higher Education Commission (HEC) Pakistan, under HEC overseas scholarship program, held by the first author, is dutifully acknowledged.

Appendix

Expressions of c ji , j = 1 . . . 6, a, p and i = 1, . . . , 3, used for the calculations of the shaking forces for the Orthoglide with reference to Fig. 5, (C and S respectively stand for sine and cosine functions):


Frédéric Gruy
email: gruy@emse.fr

Pr Frédéric Gruy

Population Balance for Aggregation Coupled with Morphology Changes

Keywords: population balance, bivariate, aggregation, restructuring

In the past, the kinetics of aggregation has been extensively studied. Aggregation rates were measured and calculated thanks to a population balance. Aggregate morphologies were measured or got by computer simulations. However, the link between the aggregation kinetics and the morphology changes with time is not so clear. The modelling of aggregation may be even more complex as restructuring of aggregates occurs. The aim of this paper is to propose a new formulation taking into account at once kinetics of collision and morphology change rate. We built a bivariate population balance with matter volume and porous volume as internal parameters. The population balance equation contains the standard collision term and a convective term representing the porous volume change. The latter is split into two contributions, which is due to the aggregation process itself and the other one is due to the restructuring. The expressions of the first contribution are determined for Brownian and shear aggregations.

Introduction

Aggregation of small particles often occurs during industrial processes. The first step of this phenomenon is the collision of particles, which can be due to differential sedimentation, Brownian motion or shear flow in a duct or an agitated vessel. The size of the resulting clusters or aggregates increases with time. However, aggregates become looser and looser and undergo breakage. There is a competition between the cohesive forces inside the aggregate and the shear stress due to the liquid flow. Thus, the final size of aggregates can be calculated from a balance between aggregation and breakage.

Aggregation and breakage have been investigated by a lot of researchers. Under certain operating conditions (supersaturation, mass transfer inside the aggregate …), consolidation or strengthening of aggregates happens. This takes place at the neck between adherent primary particles. The consolidation prevents the breakage of aggregates and leads to more porous clusters than the ones made without consolidation.

At the same time, interaction between the flowing liquid and the aggregate can modify the internal structure of aggregates, breakage being the consequence of such a strong interaction. This results in a relative motion of primary particles inside the aggregates and leads to denser aggregates (see, for instance, [START_REF] Di Stasio | Observation of restructuring of nanoparticle soot aggregates in a diffusion flame by static light scattering[END_REF][START_REF] Adachi | Restructuring of small flocs of polystyrene latex with polyelectrolyte[END_REF][START_REF] Yu | Flocculation of kaolin particles by two typical polyelectrolytes: A comparative study on the kinetics and floc structures[END_REF]).

The aim of this paper is to propose a new framework to take into account these phenomena. The corresponding modelling will result in a new formulation of the population balance, including aggregation coupled with changes in morphology. The second section of this paper reminds us of any fundamentals on Brownian and shear aggregation. The third section addresses the state of the art about aggregate restructuring and population balance with morphological parameters. The fourth one introduces the new formulation. The fifth one presents the relationship between previous and new formulations. The sixth one gives results for Brownian and shear aggregations and the last one is a discussion about its application.

Fundamentals on elementary aggregation processes

Several experimental works and computer simulations show that aggregates have a fractal morphology. In fact, an aggregate containing i identical primary particles of radius R 1 is characterised by: the fractal dimension D f and the outer diameter R i . These two parameters are linked by the following relations:

( ) 1/ 1 / f D i R R i S = (1)
The structure factor S is depending on D f . Its value is close to 1.

Another way to consider fractality uses a continuous variable (volume v) for matter:

1/ f D R v ∝
In case of a process with only aggregation, i.e. without nucleation and growth, the aggregates are formed from initially non aggregated particles, the mean volume value of which can be expressed as:

1 v 1 1 / v N = Φ
Φ and N 1 are respectively the volume fraction of solid and the number concentration of particles. Then, equation ( 1) can be written:

( )

1/ 1 1 / / f D i R R v Sv = (2) With 3 1 1 4 /3 v R π =
The main causes of particle collision are the Brownian motion of particles and the local shear in a laminar or a turbulent flow [START_REF] Elimelech | Particle deposition and aggregation, measurement, modelling and simulation[END_REF]. The aggregation rate between two (i and j)

aggregates is characterized by a kinetic constant or kernel .

, i j

K

The Brownian kernel can be approximated by [START_REF] Veerapaneni | Hydrodynamics of fractal aggregates with radially varying permeability[END_REF][START_REF] Gmachowski | Hydrodynamics of aggregated media[END_REF]:

( )( ) , , 0 1 1 , , , 2 3 i
j B B B B i j i j i j i j i j k T K K R R R R B α α - - = = + + μ (3a)
k B , T and µ are respectively the Boltzmann constant, the temperature and the fluid dynamic viscosity.

B , B i j
α is the collision efficiency. In the case of only attractive physical interaction between primary particles, a rough approximation of Eq.(3a) is:

, ,0 8 3 i j B B B B k T K K B α α = = μ (3b) With B
α having a constant value (≈0.5) [START_REF] Kyriakidis | A study of colloidal particle Brownian aggregation by light scattering techniques[END_REF].

The shear aggregation kernel is currently written as:

,0 3 , , , 4 ( ) 3
, s s s s i j i j i j i j i j K K R R α γ = = + α (4a)
, s i j α and γ are respectively the collision efficiency and the shear rate.

Several relations have been proposed for the collision efficiency. So, in the case of only attractive physical interaction between primary particles, Veerapaneni [START_REF] Veerapaneni | Hydrodynamics of fractal aggregates with radially varying permeability[END_REF] writes:

( )

3 , 4 3 j j s i j i i K R γ η η = + R (4b) i
η are the fluid collection efficiency by the aggregate i. From this work, the following approximation can be considered for large porous aggregates:

( ) ( ) 1.12 1 / f D k R R η - - = (5)
with k=0.55.

Then, equation (4b) can be written:

( ) ( ) ( ) 3/ 2 3 ' ' 1 , 1 / / f f s i j i j K v k v Sv v Sv γ π = + 1 (6) with (7)
' 1.56 0.5

f f = - and 1/ f f D =
. f' is smaller than f for standard f values. Brownian and shear flow kernels, that are homogeneous functions of volume, obey:

( ) ( , , h j j i i K v v K v v λ λ λ = )
with h=0 for Brownian kernel and h=3f' for shear aggregation kernel.

State of the art

Collision and aggregation lead to an increase of particle size whereas the restructuring decreases the porous volume of the aggregate. Various causes give place to densification: physical forces between primary particles, e.g. Van der Waals forces, Brownian motion of the primary particles leading to thermal restructuring, relative motion of the primary particles due to the shear flow, fragmentation-reaggregation steps, i.e fragmentation of a large and loose aggregate followed by aggregation of the two fragments. At our knowledge, the modelling of these phenomena has not yet been performed. However, changes in morphology connected with restructuring have been partially studied by several authors. Their results will be helpful to build a new formulation.

Firstly, we will describe the morphology change of a single aggregate or of a colliding two-aggregate set. Then, we will consider the relationship between morphology change and population balance.

restructuring laws

Empirical equations

As a first approximation, clusters keep the same fractal dimension along the aggregation process. However, it follows from some experimental evidence [START_REF] Di Stasio | Observation of restructuring of nanoparticle soot aggregates in a diffusion flame by static light scattering[END_REF][START_REF] Adachi | Restructuring of small flocs of polystyrene latex with polyelectrolyte[END_REF][START_REF] Yu | Flocculation of kaolin particles by two typical polyelectrolytes: A comparative study on the kinetics and floc structures[END_REF] that the fractal dimension value increases with time. So, the simplest way for describing this observed increase considers a relaxation law: Relaxation time may have a constant value [START_REF] Soos | Investigation of aggregation, breakage and restructuring kinetics of colloidal dispersions in turbulent flows by population balance modelling and static light scattering[END_REF], may be a function of the actual aggregate size [START_REF] Kostoglou | Bivariate population dynamics simulation of fractal aerosol aggregate coagulation and restructuring[END_REF] or may be a function of the mean aggregate size at a given time [START_REF] Selomulya | Evidence of shear rate dependence on restructuring and breakup of latex aggregates[END_REF][START_REF] Selomulya | Understanding the role of restructuring in flocculation: The application of a population balance model[END_REF][START_REF] Baldyga | Break up of nano-particle clusters in high-shear devices[END_REF].

( max / f dD dt D D ) / f τ = - (8)
The relaxation time [START_REF] Kostoglou | Bivariate population dynamics simulation of fractal aerosol aggregate coagulation and restructuring[END_REF][START_REF] Selomulya | Evidence of shear rate dependence on restructuring and breakup of latex aggregates[END_REF][START_REF] Selomulya | Understanding the role of restructuring in flocculation: The application of a population balance model[END_REF][START_REF] Baldyga | Break up of nano-particle clusters in high-shear devices[END_REF] is expressed by means of empirical equations containing fitted parameters. So, 2 or 3 fitted parameters are used for each physical contribution.

These empirical equations well represent the reality, but can be only applied to the considered systems.

Simulations

Starting from assumed restructuring mechanisms at the primary particle scale and being given an initial aggregate morphology, computer simulations lead to the morphology of the aggregate at a given time. As a result, the relaxation time can be deduced. So, in the case of thermal restructuring, Dalis et al. [START_REF] Dalis | A lattice chain model for the thermal restructuring of nanoparticle chain aggregates[END_REF] show that the relaxation time depends on the temperature, on the size and the number of primary particles in the aggregate:

( ) 1 , , f i R T τ ∝
The computer simulations also give the maximum value of fractal dimension.

Simulation of aggregation has been extensively studied [START_REF] Elimelech | Particle deposition and aggregation, measurement, modelling and simulation[END_REF]. We will mention the work of Gmachowski [START_REF] Gmachowski | Estimation of the dynamic size of fractal aggregates[END_REF]. He proposes a new approach of the aggregation modelling by considering an aggregation act for which the fractal dimension can vary along the process. The master equation is the following:

( ) ( ) 1/ 1/ 1/ 1 2 j i D D D F i j i j F + = + with ( ) ( ) 1/ 2 2 1.56 1.228 0.228 k k F A = - - - and 1 1 1 D A d - = - , 2 1 w d A d D = + - ( , , , w i i j d D D D D + = ) i
are respectively the space dimension, the trajectory dimension, the fractal dimension of the aggregate having i primary particles and the one of the resulting aggregate.

The author deduces that:

-cluster-cluster aggregation is such as i D D = ∀ -cluster-primary particle aggregation [START_REF] Gmachowski | Calculation of the fractal dimension of aggregates[END_REF] is such as the fractal dimension value increases from about 1.75 to 2.5 as the primary particle number increases. Then, the exponent in the mass-size relation is not equal to the standard fractal dimension when i<10 6 .

-the fractal dimension can be calculated [START_REF] Gmachowski | A method of maximum entropy modelling the aggregation kinetics[END_REF] if the asymptotical population density is assumed as a log-normal function (and by using the maximum entropy principle) -restructuring may be taken into account [START_REF] Gmachowski | Aggregate restructuring and its effect on the aggregate size distribution[END_REF] by changing the value of F 2 .

Kostoglou and Konstandopoulos [START_REF] Kostoglou | Evolution of aggregate size and fractal dimension during Brownian coagulation[END_REF] use a similar aggregation act for studying Brownian aggregation.

population balance equation and restructuring

Population balance equation (PBE) is a partial derivative equation, the solution of which is the population density. The latter is a function of the time t and internal parameters for homogeneous suspension. The phenomena acting on the population density are the nucleation, growth, agglomeration and breakage. The corresponding terms are classically included in the PBE. Two ways were proposed in order to take into account aggregate restructuring in the PBE.

1D population balance

More often, the particles are described with only one internal parameter (for instance, radius R if the particle shape is a sphere). can be also the volume v or the number of primary particles i in the aggregate if only the aggregation of primary particles is considered. So, the population density for homogeneous suspension depends on t et :

.

1 p 1 p 1 p ( )

, n p t

A simple way of considering the restructuring is as follows: if the relaxation time has a constant value or depends on the mean size at a given time, the fractal dimension will be recalculated at each time step by means of Eq. ( 8) and the new value will be introduced in the aggregation kernel [START_REF] Soos | Investigation of aggregation, breakage and restructuring kinetics of colloidal dispersions in turbulent flows by population balance modelling and static light scattering[END_REF][START_REF] Selomulya | Evidence of shear rate dependence on restructuring and breakup of latex aggregates[END_REF][START_REF] Selomulya | Understanding the role of restructuring in flocculation: The application of a population balance model[END_REF].

Another way due to Baldyga et al. [START_REF] Baldyga | Break up of nano-particle clusters in high-shear devices[END_REF] is to add a convective term or a growth term Gn R ∂ ∂ into the PBE. The restructuring law appears as:

/ i G dR dt = Then,
( )

1 ln / f i i f dD G R R R D dt = - (9)

2D population balance

The 2D PBE contains two internal parameters of the aggregate. Thus, the population density is written as

1 2 , p p ( ) 1 2
, , n p p t . These internal parameters can be:

1 2 , , , f p R v i p a D = =
a is the surface area of the aggregate.

The consolidation, i.e. the sintering, of the primary particles inside the aggregate leads to a decrease of the surface area of the aggregate. This can be taken into account (see for instance [START_REF] Tandon | Monte Carlo Simulation of Particle Aggregation and Simultaneous Restructuring[END_REF]

1 2 p v p a = = ): ( ) ( ) ( ) sin , , , aggregation tering n a v n a v n a v t t t ∂ ∂ ∂ = + ∂ ∂ ∂ (10a) with ( ) ( ) sin , , a tering n a v G n a v t a ∂ ∂ = - ∂ ∂ (10b)
with the relaxation law for the surface area:

( ) min / a G a a τ = --
It can be emphasized that the agglomeration term is such as:

i j i j i j i j v v v a a a + + = + = +
However, the collision integral contains a step function and integration limits that make it more complicated its use and solving.

Likewise, Kostoglou et al. [START_REF] Kostoglou | Bivariate population dynamics simulation of fractal aerosol aggregate coagulation and restructuring[END_REF] propose to use the internal parameters 1

2 f p v p D = =
in the case of restructuring. Then, (

, , ,

f f f aggregation restructuring n D v n D v n D v t t t ∂ ∂ ∂ = + ∂ ∂ ∂ (11a) with ( ) ( ) , , f D f f restructuring n D v G n D v t D ∂ ∂ = - ∂ ∂ f v ) (11b)
with the relaxation law:

( ) ( ) max / Df f G D D τ = -
Kostoglou and Konstandopoulos [START_REF] Kostoglou | Evolution of aggregate size and fractal dimension during Brownian coagulation[END_REF] consider a 2D PBE including an aggregation act.

The PBE is solved by means of a Monte-Carlo simulation. An asymptotical value for D f is reached after a time depending on the D f value of the initial aggregates.

New formulation

The aim of this paper is to establish a population balance equation by taking into account aggregation and restructuring. The new formalism will permit the introduction of various restructuring mechanisms in a simple way. The framework will be a bidimensional population balance. The two internal parameters will be chosen so that the additivity is respected along the aggregation process. We selected the matter volume v (or the number of primary particles) and the pore volume vp. Thus, v+vp is the volume of the aggregate. The pore volume is preferred to the porosity because the latter is not an additive property along the aggregation process. Aggregation without and with restructuring will be treated by means of the same formalism.

aggregation without restructuring

Aggregation process is split into two successive steps (Figure 1):

-Collision between two aggregates leading to a new particle composed in the two aggregates in contact. This step can be written as:

( ) ( ) ( , , , i i j j i j i j v vp v vp v v vp vp + → + + (12)
This step appears in the population balance equation as a standard collision term.

The kernel is the aggregation kernel ( ) ( 0 0 , , , , ,

i i j j i j i i j j K v vp v vp K v vp v vp = + +
) .

For instance, the corresponding expression for a Brownian aggregation is:

( ) ( ) ( ) ( ) ( ) 1/ 3 1/ 3 1/ 3 1/ 3 0 2 , , , 3 i i j j i i j j i i j j kT K v vp v vp v vp v vp v vp v vp µ - - ⎡ ⎤ ⎡ = + + + + + + ⎢ ⎥ ⎢ ⎣ ⎦ ⎣ ⎤ ⎥ ⎦ (13a)
and for a shear aggregation ( ) ( ) ( )

3 1/ 3 1/ 3 0 , , , i i j j i i j j K v vp v vp v vp v vp γ π ⎡ = + + + ⎢ ⎣ ⎤ ⎥ ⎦ (13b)
-Overlapping of the two initial aggregates and restructuring of the new particle. In the case of fractal aggregates with constant fractal dimension, the new particle will have a bigger volume than the sum of the initial ones. Thus, this stage corresponds to an expansion. This step will be modelized by using a convective or growth term:

( ) 0 , vp G n vp v vp ∂ - ∂ (13c) 0 vp
G is a function of vp and probably of v. represents the expansion rate of the two initial aggregates set. This convective term can be considered as a virtual restructuring.

0 vp G This modelling introduces no additional assumption such as fractal morphology, i.e relationship between pore volume and matter volume.

Asymptotic behaviour

The aggregation kernels (for instance, Brownian or shear aggregation, Eqs. (3b), (4b))

are homogeneous functions of the sizes of colliding particles. If the expansion rate is also an homogeneous function of vp and v, e.g. or , we can write:

0 m vp G vp ∝ ( 0 m vp G v vp ∝ + ) ) y ) ( ) ( 0 0 , , , , p x y x y K x y K x λ λ λ = ( ) ( 0 0 m vp vp G x G x λ λ = or ( ) ( ) 0 0 , , m vp vp G x y G x y λ λ λ =
It has been proved that large aggregates have a fractal-like morphology for long time:

3/ f D t vp v → ∞ ∝
Thus, we observe an asymptotical behaviour, as self-similarity:

( ) ( ) ( ) 3/ , , , f D t vp v n v vp t n v t t g vt β α ∞ → ∞ ∝ → →
The use of these relations into the PBE [START_REF] Gruy | Lois d'évolution de la distribution granulométrique de poudres d'oxydes métalliques calcinées en présence de gaz réactif[END_REF] leads to:

1 / f m D = - 3 3 for Brownian aggregation (14a) 2 / f m D = -
for shear aggregation (f=f') (14b)

aggregation with restructuring

Restructuring will be taken into account by means of the overall convective term:

0 vp vp vp G G G = + 1
) )

where is the restructuring (contraction) rate and is the expansion rate introduced by modelling of aggregation. An approximate expression for could be written:

1 vp G 0 vp G 1 vp G ( 1 / p vp R G K v vp v = - (16)
The p (p>0) exponent depends on the contraction mechanism. The kinetic constant has to consider the competition between contraction and consolidation, the latter annoying the former. As already used for drop coalescence [START_REF] Tavlarides | The analysis of interphase reactions and mass transfer in liquid-liquid dispersions[END_REF], the kinetic constant could be expressed as:

R K R K ( ,max 1 2 / 1 / R R K K τ τ = + (17) 1 2
, τ τ respectively represent the characteristic times for « free of sintering » contraction and consolidation.

Comparison between old and new formulations

Two aspects characterize the aggregation: the kinetics of the collision and the topology of the aggregates. The classical formulation considers them independent to each other.

The topology of aggregates is observed (for instance, by SEM) or got by computer simulations. Moreover, topology is obtained after long time aggregation. On the other hand, the kinetic rate is measured or modelized by means of a population balance.

The new formulation considers at the same time kinetics and topology by means of the expansion rate. Then, one searches the aggregation kernel ( ) ' , , ,

( ) ( ) ( ) ( ) ( ) ( ) 0 0 1 , 2 v n K v v v n v n v v dv n v K v v n v dv t ∞ ∂ = - - - ∂ ∫ ∫ , (18)
with ( ) ( ) 0 0 , , , ,

K v v K v vp v vp α = and
( )

vp g v = for t → ∞ for instance, 1/ 3 1/ 1 1 Df v vp v v S v ⎛ ⎞ ⎛ ⎞ + = ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠ (19a) or ( ) 3/ 1 1 Df v v vp vp g v v Sv ⎛ ⎞ + = ⎜ ⎟ ⎝ ⎠ (19b) g K K K x y K x y x y ⎛ ⎞ Δ ∂ ∂ + + ⎜ ⎟ ∂ ∂ ⎝ ⎠
x et y are the overall (pore and matter) volumes of colliding aggregates.

As a consequence, if the argument values (x and y) are identical, K'>K.

Moreover, a physical reason has to be mentioned for K'>K: as the first step leads to a very loose aggregate, the hydrodynamic resistance and the physical forces during its formation are very weak. Thus, the collision efficiency of the first step, and then the aggregation kernel, is higher than the one of the classical formulation.

As Eq.( 24) is difficult to apply, we will use the kernels (Eqs.(13a-b), i.e. The next section presents an estimate of the expansion rate (Eq. ( 23)) for Brownian and shear aggregation.

0 vp G

Calculation of the expansion rate

Brownian aggregation

The Brownian kernel is chosen as a constant:

,0 ,0 0 0

B B K K α α
The number concentration N j for aggregates with j primary particles obeys the relation [START_REF] Kruyt | Colloid Science[END_REF]:

( ) ( )

1 1 10 ' 1 ' j j j N N t t - - = + - ' / c t t t =
is the dimensionless time. and 10 N ( )

0 0 10 2 / c t K N α =
are respectively the initial number concentration in primary particles and the characteristic time for Brownian ( )

1 3 ' 1 1 / f v v β - = , 2 2 2 1 v σ β = + , and ' / c t t t = 3 ' 3/ 2 f c S t k π γ = Φ
by taking the initial conditions: β are not dependent on the initial 2 β values. However, if the exponent 3f' is higher than 1, gelation occurs, i.e. divergence of the ODE system (or the standard deviation). This asymptotical behaviour corresponds to the following values for 2 β and 1 / ' d dt β :

( )

1/ 1 3 ' 2, 2 f a β - = (35) ( ) ( )( )
' 3 5 ' / 1 3 ' / 2 3 '/2 1, / ' 1 3 ' 2 3 2 f f f f a d dt C f β - - - - ⎡ = = - + × ⎣ ⎤ ⎦ (36)
Let us apply these results to Eq.( 32) in order to obtain :

0 vp G ( ) ( ) ( )( ) ( )( ) ( ) ( ) ( 1
/1 3 ' 1 /1 3 ' 1 0 1 1 1 ' 1 f f f vp DBC t v v F f f G g v μ μ μ μ μ μ μ -+ - -+ - -+ + - + Φ Γ + = ) 6 ,' (37a) with ( ) ( ) ( ) ( ) ( )
3 1 3 3 1 0 1 , ' 1 1 1 2 1 f f f f F f f x x x x x x dx μ - ⎡ ⎤ + - -- - ⎡ ⎤ ⎣ ⎦ ⎣ ⎦ = - ∫ (37b) Moreover, ( ) ( )
( )

0 1 1 1 1 / / ' 1 0 / vp c c vp t d vp v dt t G g v dvp dt v v G - - = = = ' (38) Thus, ( ) ( )
( ) ( ) ( ) ( )( ) ( ) ( ) 6 / 3 1 /1 3 ' 0 1 1 ' / , ' / ' f f f vp v E f f g v v t G g v μ μ + -+ - = (39a) with ( ) ( ) ( ) ( ) ( )
1 /1 3 ' 1 1 3 ' 3/2 , ' 1 , ' f f E f f C S k F f f μ μ μ μ μ -+ - + - + - Γ + = (39b) Finally, ( ) ( ) ( ) (
)

0 1 ' / , ' / ex vp v f f v v p v G vp H + = 1 (40a) with ( ) ( ) ( ) 3 ' 1 / 3 ' 3 3 ' 1 3 ' 3 , ' 1 3 ' 3 3 3 ' f f f f f f H f f E f f f f μ μ μ μ μ - + ⎛ ⎞ + - + = ⎜ ⎟ - + + ⎝ ⎠ (40b)
and

1 3 ' 3 1 3 3 ' f f ex f f μ μ - + = - + (40c)
If f'=f, the asymptotic behaviour with the exponent 2 / 3

f D - is verified.
The figure 2a represents the function H(f,f') (Eq. ( 40b)). We note that the effect of f on is weaker than the one of f'. The figure 2b represents the exponent ex of the power law (Eq. (40c)) versus D 0 vp G f . Following Torres [START_REF] Torres | Simulations of coagulation in viscous flows[END_REF], simulations of shear aggregation lead to a fractal dimension (D f =1/f) equal to 1.8, while Potanin [START_REF] Potanin | On the mechanism of aggregation in the shear flow of suspensions[END_REF] suggests the value 1.98. We will take the intermediate value 1.9.

From Eqs [START_REF] Kyriakidis | A study of colloidal particle Brownian aggregation by light scattering techniques[END_REF], (35-36), (40a-c), one deduces:

( ) ( ) ( 0.97 0 1 ' / 1.5 / vp v v v p v G vp = + ) 1 (41)

Discussion

The new formulation has been applied in the case of Brownian aggregation for which experiments, calculations and computer simulations were achieved by other investigators. Monte Carlo Simulations (MCS) have been performed to solve the particle dynamics equation for coagulation of particles undergoing morphology changes (Eqs. 20 and 29). Here we use the algorithm suggested by Tandon and Rosner [START_REF] Tandon | Monte Carlo Simulation of Particle Aggregation and Simultaneous Restructuring[END_REF] when they studied aggregation coupled with surface area reduction. In our case, the sintering kinetic law was replaced by the porous volume expansion law (Eq. 29). The aggregation kernel has been taken as a constant. The choice of the time step in the MCS needs a particular consideration. As a consequence we used the method proposed by Smith and Matsoukas [START_REF] Smith | Constant-number Monte Carlo simulation of population balances[END_REF] in which the number of particles is kept constant. To start the simulations we consider a set of N (N=5000) particles, each of unit volume and porous volume equal to 0. During aggregation each aggregate is characterized by v and vp. One deduces its fractal dimension by means of Eq. 19a.

From these simulations we show that:

-At the beginning of aggregation small and very porous aggregates are present and may survive. This is due to the convective term (Eq.13c), that allows certain aggregates to increase their porous volume whereas the other aggregates collide each other. The growth of the porous volume stops when the aggregate disappears due to collision. Certain small aggregates, waiting for the collision, get a higher porosity. However, their amount is smaller than 1%.

-All the small aggregates (i>20) have the same behaviour. Let us consider aggregates with a given number of primary particles (for instance i=50). They suddenly appear and their fractal dimension value is in the range [1.9-2.0]. Their number rapidly increases, then it slowly decreases (figure 3). Whatever the time the fractal dimension distribution of these aggregates is narrow. The standard deviation of the distribution is about 0.04. During the aggregation the mean fractal dimension decreases. The value corresponding to the maximum number of aggregates (i=50) is 1.65.

-After a long time, the small aggregates disappear. The fractal dimension distribution is within the range [1.6-1.8] and the mean fractal dimension tends to the value 1.80 (figure 4). The distribution is skewed with a tail for the smallest values (Df<1.7).

These results have been compared to previous works. A controversy concerns the morphology of small aggregates. Adachi et al. [START_REF] Adachi | Applicability of fractals to the analysis of the projection of small flocs[END_REF] have studied the DLCA mechanism.

They observe an increase of the fractal dimension from 1 for very small aggregates (i<5) to 1.75-1.80 for aggregates with i>50. Yang and Biswas [START_REF] Yang | Computer simulation of the aggregation and sintering restructuring of fractal-like clusters containing limited numbers of primary particles[END_REF] did the same observation. Lattuada et al [START_REF] Lattuada | A simple model for the structure of fractal aggregates[END_REF] built aggregates by simulation and showed that the aggregates are fractal-like if i>20. Our MCS of aggregation dynamics are partially in agreement with the results of these investigators. We observe small aggregates with low Df value and a change of Df with time which is the same for the aggregates with i>20.

However, in our modelling, the morphology is a dynamic property; the Df value of aggregates is close to 1.8 during their formation (see figure 3), but it becomes smaller along their removal due to collisions.

All the aggregates contain a large number of primary particles after a long time. The fractal dimension distribution varies very slowly whereas matter volume distribution tends to an asymptotic shape [START_REF] Tandon | Monte Carlo Simulation of Particle Aggregation and Simultaneous Restructuring[END_REF]. This result was expected because of the way for calculating the expansion rate.

Summary and conclusion

The morphology changes are taken into account thanks to a 2D population balance, the internal variables of which are the volume of matter and the porous volume of the aggregate. The mechanism of aggregation is split up in two steps: the collision itself characterized by the collision kernel without collision efficiency (Eqs.(13a-b)) and the change of the porous volume of the aggregate (Eq.( 15)). The latter contains a contribution due to the porous volume increase occurring for all aggregation processes and another one due to the restructuring itself. The bivariate population balance

+

Step 1

Step 2 * E-mail address: gruy@emse.fr

  are respectively the actual fractal dimension, its maximum value and a relaxation time. The relaxation time includes the above-mentioned physical contributions.

  v vp and the expansion rate leading to the same aggregation dynamics and topology than the ones obtained by the classical formulation. 0 vp G The classical formulation uses the 1D population balance equation:

1 1α

 1 = ) corresponding to the physical characteristics of step 1.

  constant values depending on f'. The asymptotical values 2,a

Figure 1 :

 1 Figure 1: decomposition into two steps in the new formulation of aggregation

Figure 3 :

 3 Figure 2b: exponent ex of the power law (Eq. 40c) versus D f .

Figure 4 :

 4 Figure 4: fractal dimension distribution for long time (time scale following [19])

* E-mail address: gruy@emse.fr

Acknowledgments

We acknowledge the Agence Nationale de la Recherche for support of this research.

The new formulation uses the 2D population balance equation:

' , ' , , , ' ,

K v vp v v vp vp n v vp n v v vp vp d vd vp n v vp K v vp v vp n v vp dvdvp

' , , , , , ,

Thus, the problem to solve is the following:

Which are the expressions of ( ) 0 , vp G v vp and , such as ( 18) ≡ (20) and ( )? ' K ( )

The relation between the 2D and 1D population densities ( t

n v is the solution of Eq. [START_REF] Kostoglou | Evolution of aggregate size and fractal dimension during Brownian coagulation[END_REF].

By using a moment method and by applying to Eq.( 20) the operator O defined by , one obtains:

, , 0

for p=0 : 18) and ( 20) are equivalent if:

By subtracting Eq.( 21) (with p=0) at Eq.( 21) (with p=1):

Eq.( 23) calls for several remarks:

i.

The method leading to Eq. ( 23) belongs to a mean field theory.

ii. iii.

If g(v) obeys a power law as

Thus, aggregation proceeds with expansion or porous volume increase.

( )

Only zero and first order moments of the 2D population density were used in order to establish Eq.( 23) because the higher order moments ( ) are not undoubtedly bounded. By contrast, the first order moment obeys the inequality:

e, the pore volume fraction is smaller than 1).

v.

The new formulation contains two successive steps (collision and expansion): K' is the kinetic constant of the collision rate and is the expansion rate. As K, in the classical formulation, includes the two steps at the same time, the inequality K'>K should be verified, contrary to the previous result (K'=K). This possible discrepancy is due to the final aggregate volume (i.e., the arguments in the kernels K and K') considered in each formulation. The following proof can be done:

* E-mail address: gruy@emse.fr

where v is the final volume. Then,

' , , , , , ,

By taking a first-order Taylor series of F 1 (F 1 is replaced by K) and a symmetrisation for the volume:

aggregation. But, the number concentrations for Brownian aggregation have the interesting property, denoted P1:

( )

With continuous variables, the population density can be written as:

and ( ) ( ) ( ) ( )

By using Eqs 23,19b and property P1, one gets the expression:

( ) ( ) ( )

For long time ( ), this expression can be transformed into: / 1

The asymptotic behaviour with the exponent 1 / 3

The fractal dimension for Brownian aggregation without restructuring is equal to 1.8 (S=0.54). This value is coming from computer simulations. Thus, ( )

shear aggregation

The general expression of the aggregation kernel is that of Veerapaneni (Eq.( 6)).

Thus, two exponents will be introduced: one f for the fractal geometry of aggregates and another one f' for the kinetic constant:

( )

* E-mail address: gruy@emse.fr

The population density for shear aggregation has no property such as P1. Then, the exact solution of 1D population balance will be replaced by an approximate function having the property P1':

h is a homogeneous function.

The Gamma function, which is widely used for approximating experimental particle size distributions, has such a property P1' [START_REF] John | Techniques for the reconstruction of a distribution from a finite number of its moments[END_REF]:

, v σ are the mean value and the standard deviation, while

These quantities can be expressed from the distribution moments M j , that depend on the time:

A is a normalization constant:

Thus, the calculation of necessitates the three moments M 0 vp G 0 , M 1 and M 2 . In the case of shear aggregation without fragmentation, the moments obey:

Following Park and Lee [START_REF] Park | Evolution of particle size distribution due to turbulent coagulation[END_REF][START_REF] Park | Asymptotic particle size distributions attained during coagulation processes[END_REF][START_REF] Park | Change in particle size distribution of fractal agglomerates during Brownian coagulation in the free-molecule regime[END_REF] and Gruy [START_REF] Gruy | Méthode approchée pour la reconstruction de la densité de population : cas de l'agrégation-agglomération[END_REF], the three first moments, the mean value and the standard deviation can be calculated thanks to the ordinary differential equation (ODE) system:

with * E-mail address: gruy@emse.fr equation is solved with the corresponding initial conditions and boundary conditions. As the contribution of the aggregation process to the porous volume change was established in this paper (Eqs. ( 29) and ( 41)) for two important aggregation mechanisms, the one of the restructuring has to be determined as the term . This task is in progress.

1 vp G