Thanh-Trung Vo

Patrick Mutabaruka

Saeid Nezamabadi

Jean-Yves Delenne

Edouard Izard

Roland Pellenq

Farhang Radjai
email: franck.radjai@umontpellier.fr

Mechanical strength of wet particle agglomerates

Keywords:

de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The agglomeration or granulation of solid particles is used in many sectors including powder metallurgy, chemical engineering, pharmaceutical industry, and iron-making processes to produce agglomerates or granules from small particles. The agglomerates are used to improve flow properties, enhance permeability for the interstitial gas or reduce segregation in the presence of several types of particles [START_REF] Iveson | Nucleation, growth and breakage phenomena in agitated wet granulation processes: a review[END_REF][START_REF] Herminghaus | Dynamics of wet granular matter[END_REF][START_REF] Litster | The science and engineering of granulation processes[END_REF][START_REF] Liao | Effect of adding a small amount of liquid on density-induced wet granular segregation in a rotating drum[END_REF] . The binding material is generally a liquid, which is mixed with the primary particles in a granulator [START_REF] Degrève | Spray-agglomeration of npk-fertilizer in a rotating drum granulator[END_REF][START_REF] Osborne | Binder addition methods and binder distribution in high shear and fluidised bed granulation[END_REF][START_REF] Pawar | An experimental study of droplet-particle collisions[END_REF][START_REF] Pashminehazar | Three dimensional characterization of morphology and internal structure of soft material agglomerates produced in spray fluidized bed by X-ray tomography[END_REF] . Hence, the mechanical strength of the 'raw granules' produced by the agglomeration of primary particles is ensured by the action of capillary and viscous forces due to the presence of liquid clusters in the pore space [START_REF] Iveson | Growth regime map for liquid-bound granules[END_REF][START_REF] Delenne | Liquid clustering and capillary pressure in granular media[END_REF][START_REF] Lian | A theoretical study of the liquid bridge forces between two rigid spherical bodies[END_REF][START_REF] Willett | Capillary bridges between two spherical bodies[END_REF][START_REF] Št Ěpánek | Distribution and accessibility of binder in wet granules[END_REF][START_REF] Scheel | Liquid distribution and cohesion in wet granular assemblies beyond the capillary bridge regime[END_REF] .

The raw granules should be in a densely and homogeneously packed state in order to support the loads to which they are subjected during subsequent operations and give rise to strong solid granules upon sintering [START_REF] Iveson | The dynamic strength of partially saturated powder compacts: the effect of liquid properties[END_REF][START_REF] Thornton | A theoretical model for the stick/bounce behaviour of adhesive, elastic-plastic spheres[END_REF][START_REF] Liu | Coalescence of deformable granules in wet granulation processes[END_REF][START_REF] Ghadiri | Editorial: special issue -agglomeration[END_REF][START_REF] Rojek | Modeling of a sintering process at various scales[END_REF] . The granule strength is controlled by two types of parameters: 1) process parameters such as the method of mixing particles with the liquid (depending on the granulator) [START_REF] Rahmanian | Effects of process parameters on granules properties produced in a high shear granulator[END_REF][START_REF] Ileleji | Experimental investigations towards understanding important parameters in wet drum granulation of corn stover biomass[END_REF] and 2) material parameters such as the nature of the liquid and the size distribution of the particles [START_REF] Osborne | Binder addition methods and binder distribution in high shear and fluidised bed granulation[END_REF][START_REF] Ileleji | Experimental investigations towards understanding important parameters in wet drum granulation of corn stover biomass[END_REF][START_REF] Ramachandran | Experimental studies on distributions of granule size, binder content and porosity in batch drum granulation: inferences on process modelling requirements and process sensitivities[END_REF] . The physical processes governing the growth of granules are complex due to the dynamic nature of granulation involving the collisions of particles and transport of the binding liquid inside a partially wet granular material [START_REF] Iveson | Growth regime map for liquid-bound granules[END_REF][START_REF] Iveson | Fundamental studies of granule consolidation part 1: effects of binder content and binder viscosity[END_REF][START_REF] Iveson | Growth regime map for liquid-bound granules: further development and experimental validation[END_REF][START_REF] Mishra | A preliminary numerical investigation of agglomeration in a rotary drum[END_REF][START_REF] Wauters | Liquid distribution as a means to describing the granule growth mechanism[END_REF][START_REF] Št Ěpánek | Computer simulation of granule microstructure formation[END_REF] . The agglomeration process in a fluidized bed can be correctly modeled by assuming binary collisions between wet particles [START_REF] Pawar | An experimental study of droplet-particle collisions[END_REF][START_REF] Pashminehazar | Three dimensional characterization of morphology and internal structure of soft material agglomerates produced in spray fluidized bed by X-ray tomography[END_REF][START_REF] Mikami | Numerical simulation of cohesive powder behavior in a fluidized bed[END_REF] whereas in a rotating drum [START_REF] Aguado | Best Available Techniques (BAT) Reference Document for Iron and Steel Production[END_REF] , for example, the agglomeration occurs inside a dense cohesive granular flow whose rheology has only recently been studied by careful experiments [START_REF] Lefebvre | Erosion dynamics of a wet granular medium[END_REF][START_REF] Sastry | Investigation of the layering mechanism of agglomerate growth during drum pelletization[END_REF][START_REF] Wauters | Growth and compaction behaviour of copper concentrate granules in a rotating drum[END_REF][START_REF] Govender | Granular flows in rotating drums: a rheological perspective[END_REF] and simulations [START_REF] Rahmanian | Effects of process parameters on granules properties produced in a high shear granulator[END_REF][START_REF] Berger | Scaling behaviour of cohesive granular flows[END_REF] .

The effect of material parameters on the granule strength reflects both the strength of cohesive bonds between particles and the granular texture, i.e. the organization of the primary particles inside the granule. For example, the strength of pharmaceutical tablets, measured by quasi-static compression between two platens, declines with porosity, which is a function of the consolidation pressure used to manufacture the tablet [START_REF] Song | Effects of drying technique on extrusion-spheronisation granules and tablet properties[END_REF][START_REF] Fichtner | Particle size distribution and evolution in tablet structure during and after compaction[END_REF] . In the same way, the compressive strength of cohesive powder mixtures is an increasing function of the relative density [START_REF] Wu | A simple predictive model for the tensile strength of binary tablets[END_REF] . The diametrical compression test is a simple way of measuring the tensile strength of powder compacts as the compressive strength (the stress at incipient failure of a granule) is proportional to the tensile strength of the granule [START_REF] Tsoungui | Numerical model of crushing of grains inside two-dimensional granular materials[END_REF] . This method has, however, been mostly used to study the fracture stress of brittle materials composed of particles glued via solid bonds [START_REF] Fu | An experimental study of the impact breakage of wet granules[END_REF][START_REF] Nguyen | Bonded-cell model for particle fracture[END_REF][START_REF] Cantor | Three-dimensional bonded-cell model for grain fragmentation[END_REF] . In contrast, wet granules have been much less investigated. While their cohesive strength is mainly controlled by the Laplace pressure and surface energy of the liquid phase, the effects of granular texture resulting from the granulation process and material parameters are still poorly understood.

In this paper, we are interested in the influence of material parameters on the strength of wet spherical agglomerates in which the liquid is assumed to be distributed as binary bridges joining eligible particle pairs. We use the Discrete Element Method (DEM) with a capillary cohesion law in which the attraction force is an explicit function of the gap between particles and liquid-vapor surface tension, and the amount of liquid is mainly accounted for by a debonding distance. By simulating the diametrical compression of spherical agglomerates, we find that, due to particle rearrangements, they show a plastic behavior with a threshold that we analyze as a function of friction coefficient and size span of primary particles. We also introduce a model for the compressive granule strength that accounts for particle size distribution and we discuss the role of the class of fine particles for the plastic threshold.

In Section 2 , we introduce the numerical model and procedures used to prepare and simulate spherical agglomerates. In Section 3 , we discuss the evolution of the granule strength as a function of axial strain and the effects of particle size span and liquid volume. In Section 4 , we introduce an analytical model of granule strength in terms of texture parameters. We conclude in Section 5 with a short summary of salient results and routes to further research.

Numerical method and procedures

The Discrete Element Method (DEM) has been extensively used for the simulation of granular materials [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF][START_REF] Thornton | Quasi-static shear deformation of a soft particle system[END_REF][START_REF] Herrmann | Modeling granular media with the computer[END_REF][START_REF] Radjai | Discrete-element Modeling of Granular Materials[END_REF] . It is based on the step-wise integration of the equations of motion for all particles by taking into account the particle interactions. In advanced applications of the DEM, it is now possible to implement also the presence of an interstitial fluid or a solid binding matrix [START_REF] Lian | A theoretical study of the liquid bridge forces between two rigid spherical bodies[END_REF][START_REF] Radjai | Force Transmission in Cohesive Granular Media[END_REF] . However, such applications require substantially more computation power and memory in order to discretize the degrees of freedom associated with the interstitial phase. For this reason, in DEM simulations of granular processes, it is necessary to set up a modeling strategy by making appropriate choices that allow for a balance between computational efficiency (large number of particles) and physical realism. In the case of unsaturated wet granular materials, it is found that the fluid phase can be correctly represented by its cohesive and viscous effect in the particle-particle interactions [START_REF] Richefeu | Shear Strength of Unsaturated Soils: Experiments, Dem Simulations, and Micromechanical Analysis[END_REF] . Hence, we rely on this approach to model the binding liquid in the granulation process.

On the other hand, the granule strength depends on its internal structure, which is controlled to some extent by the granulation device. Here, we are interested in the simpler case of 'ideal' granules of spherical shape where the primary particles and binding liquid are homogeneously distributed. This simplification allows us to investigate the effect of basic parameters such as the particle size distribution and friction between particles on the strength in the absence of specific granulation process parameters. Furthermore, in association with texture analysis, the results of this investigation can provide a reference behavior against which the effects of process parameters can be quantified in the next step. We use the DEM with a capillary force law and a simple algorithm for the construction of spherical granules. It is worth mentioning here that this approach has been used for the simulation of an assembly of wet agglomerates in the pendular state in application to powder processes such as the coalescence of granules upon collision [START_REF] Lian | Discrete particle simulation of agglomerate impact coalescence[END_REF] and impact breakage of crystalline agglomerates [START_REF] Kafui | Numerical simulations of impact breakage of a spherical crystalline agglomerate[END_REF][START_REF] Mishra | Impact breakage of particle agglomerates[END_REF][START_REF] Thornton | Numerical simulations of diametrical compression tests on agglomerates[END_REF][START_REF] Liu | Impact breakage of spherical, cuboidal and cylindrical agglomerates[END_REF] . In this section, we describe both the numerical method and the procedure that we used to create our 'ideal' granules.

Numerical method

In the DEM, the particles are modeled as rigid particles interacting via visco-elastic force laws relating the contact force to the local strain expressed from the relative particle displacements. The simulation of rigid particles requires a stiff repulsive potential and high time resolution. The motion of each particle i with radius R i is governed by Newton's second law:

m i d 2 s i dt 2 = f n n + f t t , ( 1
)
where f n and f t are the normal and tangential components, respectively, and n and t are the corresponding unit vectors pointing in the normal and tangential directions. m i and s i are the mass and position vector of particle i , respectively. In our simulations, we used the velocity-Verlet time-stepping scheme [START_REF] Radjai | Discrete-element Modeling of Granular Materials[END_REF][START_REF] Duran | Powders, and Grains: An Introduction to the Physics of Granular Materials, Partially ordered systems[END_REF] . The force laws involve normal repulsion, normal damping, capillary cohesion and Coulomb friction. The normal force f n has three different sources:

f n = f e n + f d n + f c n . ( 2
)
The first term in this equation is the normal repulsive contact force. The normal repulsion force is a linear function of the normal elastic deflection δ n approximated by the overlap between two particles [START_REF] Richefeu | Shear strength properties of wet granular materials[END_REF][START_REF] Vo | Strength of wet agglomerates of spherical particles: effects of friction and size distribution[END_REF] :

f e n = k n δ n , ( 3
)
where k n is the normal stiffness constant. The normal damping force f d n is assumed to be a viscous force proportional to the relative normal velocity ˙

δ n : f d n = γ n ˙ δ n
, where γ n is the damping coefficient. Both these forces disappear when there is no overlap, i.e.

for δ n > 0.

The last term in the Eq. ( 2) is the capillary cohesion force f c n due to the liquid bond between two particles. It depends on the gap, which we denote by δ n as for overlaps but with negative values (positive values representing a contact deflection), liquid volume V b , surface tension γ s , and the particle-liquid-gas contact angle θ ; see Fig. 1 [11,47,54] . The capillary cohesion force can be determined by integrating the Laplace-Young equations. Various solutions have been proposed for this equation including recent new analytical solutions [START_REF] Gagneux | Analytic calculation of capillary bridge properties deduced as an inverse problem from experimental data[END_REF][START_REF] Lian | The capillary bridge between two spheres: new closed-form equations in a two century old problem[END_REF] . We use the following explicit expression of the capillary force, which is in good agreement with experiments [START_REF] Mikami | Numerical simulation of cohesive powder behavior in a fluidized bed[END_REF][START_REF] Richefeu | Stress transmission in wet granular materials[END_REF] :

f c n = -κ R, for δ n < 0 , -κ R e -δn /λ , for 0 ≤ δ n ≤ d rupt , 0 , for δ n > d rupt , ( 4
)
where R = R i R j is the geometrical mean of the particle radii R i and R j , and [START_REF] Richefeu | Stress transmission in wet granular materials[END_REF] κ = 2 πγ s cos θ [START_REF] Degrève | Spray-agglomeration of npk-fertilizer in a rotating drum granulator[END_REF] where d rupt is the debonding distance, defined as the distance beyond which the capillary bridge is unstable and the bond breaks. It is related to liquid volume by [START_REF] Richefeu | Shear strength properties of wet granular materials[END_REF][START_REF] Vo | Strength of wet agglomerates of spherical particles: effects of friction and size distribution[END_REF] :

d rupt = 1 + θ 2 V 1 / 3 b ( 6
)
The length λ is the factor that controls the exponential falloff of the capillary attraction force in Eq. ( 4) . This factor is a function of liquid volume V b , the harmonic mean radius (

R = 2 R i R j / (R i + R j )
) and the size ratio between two particles in contact

r = max { R i /R j ; R j /R i } : λ = c h (r) V b R 1 / 2 (7)
This form fits well the capillary force obtained from direct integration of the Laplace-Young equation by setting h (r) = r -1 / 2 and c 0.9 [START_REF] Richefeu | Shear strength properties of wet granular materials[END_REF][START_REF] Richefeu | Stress transmission in wet granular materials[END_REF][START_REF] Radjaï | Bond anisotropy and cohesion of wet granular materials[END_REF]] . The capillary force declines in absolute value as a function of the gap δ n up to the debonding distance d rupt .

For the tangential force f t between particles in contact, we use a combination of a tangential elastic force k t δ t , where k t is the tangential stiffness constant, a Coulomb friction force threshold

μ( f n -f c n ) ,
f t = -min (k t δ t + γ t ˙ δ t ) , μ( f n -f c n ) . ( 8
)

Ideal granules

In order to create homogeneous agglomerates of particles of spherical shape, we first prepared large samples by means of isotropic compaction inside a box. The primary particles are spheres with their diameters defined in a range [ d min , d max ] with a given size ratio α = d max /d min . The size distribution is assumed to be uniform by particle volume fractions, i.e. with all size classes having the same volume. As the total volume of particles in each size class d i is proportional to d 3 i , the uniform distribution by volume fractions is defined by the condition that the product n i d 3 i is a constant and i P i = 1 , where P i is the numerical fraction of particles in class i . These conditions lead to the following distribution P of particles of diameters d :

P (d) = 2 d 2 max d 2 min d 2 max -d 2 min d -3 (9)
This distribution has the advantage of allowing the particles belonging to each size class to be correctly represented by their volume, i.e. a large number of small particles and a small number of large particles. This distribution leads to a dense packing as the pore space between large particles is filled by smaller particles [START_REF] Voivret | Space-filling properties of polydisperse granular media[END_REF][START_REF] Voivret | Multiscale force networks in highly polydisperse granular media[END_REF] . We constructed different samples with five different values of the size ratio α= 1, 2, 3, 4 and 5. Note that the mean particle diameter d is a function of d min and d max : The largest particle size was kept to a constant value d max = 10 μm, and d min was decreased from 10 μm to 2 μm. Eq. [START_REF] Delenne | Liquid clustering and capillary pressure in granular media[END_REF] shows that, since d max is fixed, the average diameter d declines by a factor 3 when α is increased from 1 to 5.

d = d max d min d P (d) δd = d min 2 α 1 + α = d max 2 1 + α (10)
For isotropic compaction, the particles were introduced in a box and equal compressive stresses σ 0 were applied to the box walls without gravity until a packing in static equilibrium was achieved. During this step, the capillary force was set to zero and the friction coefficient to 0.1 at all contacts between particles in order to obtain a dense sample. Once all particles reached a state of static equilibrium, a spherical probe was placed in the center of the box and its radius was increased until exactly 50 0 0 particles were inside the probe. These particles were then extracted and allowed to relax with the capillary force law activated. The common numerical parameters in these simulations are θ = 0 , γ n = γ t = 5 × 10 -5 Ns/m and γ s = 0 . 072 N/m (water).

We subjected the granules prepared by the above procedure to axial compression between two platens, as illustrated in Fig. 2 (a). The bottom platen is fixed and a downward motion is applied to the top platen with a constant velocity v 0 = 0 . 1 ms -1 . Hence, with time step δt = 10 -9 s used in our simulations, the total downward displacement during one time step is v 0 δt = 10 -10 m, which is 10 -5 times the size of the primary particles. This means that the diametrical deformation applied to the granules is slow enough to allow for a quasi-static compression test.

To see how quasi-static is the compression, one may also compare the average elastic force increment δf e between particles with the cohesive force f c = πγ s d. The incremental force between particles is simply given by the contact normal stiffness k n multiplied by the average normal displacement δ n at the contact points between primary particles. The latter is given by the mean diametral deformation v 0 δt / D g times the primary particle size d . We thus can define the dimensionless index:

I e = δ f e f c = k n δt πγ s v 0 D g . ( 11
)
For monodisperse particles, from our parameter values we get I e 2 . 5 × 10 -3 , which means that the force increments are generally small compared to the cohesion force. Finally, the time associated with the dissipation of kinetic energy is given by the ra- In Fig. 3 , we display a snapshot of a granule for α = 1 at the end of diametrical compression. At the beginning of the test, at most three primary particles of the granule are in contact with either of the two platens. As the compression proceeds, the granule spreads without breaking between the two platens and the number of contacts between the granule and each platen increases. The mechanical response of the granule under diametrical compression requires a measure of the axial force and vertical deformation. The vertical deformation in compression is given by

ε = ln 1 + h D g ≈ h D g ( 12
)
where h = v 0 t is the total downward displacement of the top platen (the lower platen being fixed). We also measure the vertical force component F between the granule and the top platen by summing up the normal forces between the primary particles and the platen. Let us note that, as the granule is in static equilibrium, all forces at each horizontal layer of the granule are balanced so that the vertical force acting between all horizontal layers is equal to F . By dividing this force by the sectional area π a 2 of the granule, where a is the radius of the central section of the granule perpendicular to the compression axis, we get the average vertical stress σ zz = F / (π a 2 ) in the center of the granule. The value of a is estimated from the positions of the particles located at the boundary of the actual central section.

The stress can also be obtained from the values of normal forces and branch vectors (vectors joining particle centers) using [START_REF] Rothenburg | A Micromechanical Definition of the Cauchy Stress Tensor for Particulate Media[END_REF][START_REF] Christoffersen | A micromechanical description of granular material behavior[END_REF][START_REF] Mehrabadi | On statistical description of stress and fabric in granular materials[END_REF][START_REF] Moreau | New Computation Methods in Granular Dynamics[END_REF][START_REF] Staron | Multi-scale analysis of the stress state in a granular slope in transition to failure[END_REF]

σ zz = 1 V g N b k =1 f k z k z = n b f k z k z k , ( 13
)
where V g is the volume of the granule, N b is the number of bonds, n b = N b /V g is the number density of bonds, f k z and k z are the zcomponents of the bond force vector and branch vector, respectively, at the contact k , including all internal contacts as well as the contacts with the platens. The symbol k denotes averaging over all contacts k in the volume. We find that σ zz σ zz during vertical compression. Note that we can define a characteristic cohesive stress σ c from the capillary force and the mean particle diameter d :

σ c = κ d ( 14
)
This stress depends on d and, by virtue of Eq. ( 10) , varies linearly with α as Its value increases in our simulations by a factor 3 as α is increased from 1 to 5 by reducing d min .

σ c = κ d max 1 + α 2 (15)

Granule strength

Fig. 4 shows the evolution of the mean vertical stress σ zz normalized by the reference cohesive stress σ c as a function of the vertical strain ε for different values of granule parameters. In all cases, σ zz first increases with strain and reaches a plateau more or less fast depending on the values of parameters. It then declines smoothly as a result of the gradual loss of cohesive contacts. The stress plateau is a signature of plastic deformation due to particle rearrangements. However, this plateau does not persist as the cohesive contacts break apart and do not heal. This irreversible character of cohesive contacts is assumed to reflect the fact that the liquid contained in a capillary bridge between two particles is not physically available (e.g. due to evaporation or drainage) once the bridge disappears. This may not always be the case as liquid drops may survive at the surface of the particles and migrate through the vapor phase or by diffusion at the surface of the particles to the newly-formed contacts during a continuous deformation of the granule.

In this way, the observed behavior can globally be qualified as ductile with a well-defined plastic plateau before the beginning of a progressive loss of cohesion at strains above 0.2. This means that the debonding events between primary particles do not lead to spontaneous formation of a fracture surface. We also see that the plastic stress threshold σ p is of the order of 0.4 σ c for monodisperse particles and 0.3 σ c for α = 5 . On the plateau, the particles are well-connected with one another, and the loss of one or two cohesive bonds of a particle does not lead to macroscopic rupture. The values of parameters affect not only the plastic threshold but also the initial build-up and later fall-off of the stress.

The effect of size span α on the plastic strength is displayed in Fig. 5 for several values of the debonding distance d rupt . We see that the ratio σ p / σ c declines by nearly the same amount in all cases as α increases from 1 to 5. This is, however, only a small relative loss of strength with respect to σ c given that, according to Eq. ( 15) , σ c is multiplied by 3 when α increases from 1 to 5.

Hence, in absolute value, the cohesive strength of the granule increases by nearly a factor 3 as shown in the inset to Fig. 5 . Fig. 6 displays the cohesive (plastic) strength as a function of d rupt for different values of α. We see that this dependence is linear and quite weak for all values of α. As we shall see in the next section, this increase of cohesive strength reflects that of the connectivity of primary particles by liquid bonds as the debonding distance increases.

Analytical model

The granule strength under diametrical compression reflects the microstructure of the granule, which depends in our simulations on the effects of size span and debonding distance (related to the amount of liquid). The granular microstructure can be described in terms of various scalar and tensorial variables such as the coordination number, packing fraction and fabric tensor [START_REF] Radjai | Modeling granular materials: century-long research across scales[END_REF] . Fig. 7 shows that the initial value of the wet coordination number Z 0 is an increasing function of both α and the debonding distance d rupt .

The wet coordination number is defined as the number of capillary bonds per particle. This is slightly above the "dry" coordination number, which accounts for only the geometrically touching particles. We see that the increase of Z 0 is more significant with d rupt than with α. Its increase is nonlinear with d rupt and it levels off around Z 0 = 12 .

In order to get a more clear understanding of the relation between the connectivity of particles and the plastic strength, we may use the expression [START_REF] Št Ěpánek | Distribution and accessibility of binder in wet granules[END_REF] of the stress tensor applied in our case to the whole volume V g of the granule. To obtain an analytical expression of the cohesive stress, we consider the cohesive forces f i j c = πγ s d i d j between particles of diameters d i and d j . Most of the cohesive strength is carried by this cohesive force acting between particles that are in contact. But many particles are connected by capillary bridges with a nonzero gap, where the cohesive force is below f i j c [START_REF] Radjaï | Bond anisotropy and cohesion of wet granular materials[END_REF] . At plastic threshold, a large number of bonds along the directions perpendicular to the compression axis are tensile whereas many others along the compression axis are compressive.

The ductility of the particle at failure is induced by the effect of tensile bonds that prevent the primary particles from sudden rupture. We will account below for these bonds by a prefactor estimated from simulations. The vertical component at plastic threshold is given by

σ p = σ zz = n b f (i j) c (i j) (i j) cos 2 θ (i j) z (i j) , (16)
where the summation runs over all contacts ( ij ) and we have

i j = 1 2 (d i + d j ) (17)
The angle θ (i j) z is the angle between the contact normal and the vertical axis. This angle is assumed to be uncorrelated with f (i j) c (i j) . For a nearly isotropic distribution of contact orientations, we have cos 2 θ z = 1 / 3 , so that

σ p = κ 12 n b (d i + d j ) d i d j (i j) ( 18
)
The number density (number per unit volume) n b of bonds can be estimated as the number of bonds per particle Z /2 divided by the free volume V f occupied by each particle. This volume is simply the particle volume divided by the mean packing fraction of the granule such that the sum of all particle free volumes is equal to the granule volume V g . Hence,

V f = π 6 d 3 (19) so that n b = Z 2 V f = 3 Z π d 3 (20)
Introducing this expression in Eq. ( 18) , we get

σ p = Z κ 4 π (d i + d j ) d i d j d 3 (21)
The geometric factor (d i + d j ) d i d j can be evaluated as [START_REF] Richefeu | Shear strength properties of wet granular materials[END_REF] (

d i + d j ) d i d j = d 3 / 2 i d 1 / 2 j + d 3 / 2 j d 1 / 2 i = 2 d 3 / 2 d 1 / 2 (22)
with the underlying assumption that there is no size segregation so that d i and d j are not correlated. Hence, we finally get

σ p = η sZ 2 π κ d = η sZ 2 π σ c (23)
The prefactor η is introduced here to account for the bonds with a nonzero gap where the cohesive force is below f i j c . Since Z denotes the coordination number for capillary bonds, Eq. ( 23) with η = 1 can be considered as an upper bound for the plastic strength. In the same way, for Z equal to the coordination number only for contacts, Eq. ( 23) provides a lower bound of the plastic strength.

As Z varies between 8 and 12 (see Fig. 7 ), we thus expect that η is generally below but close to 1.

Fig. 8 displays the values of η computed from our simulations for all values of α and d rupt . We see that η declines slightly with α and d rupt . For small values of d rupt , its value is 0.7 for monodisperse granules and 0.6 for polydisperse granules. With this prefactor, Eq. ( 23) predicts that the dependence of the normalized strength σ p / σ c with respect to α and d rupt is mediated by that of sZ . Fig. 9 displays sZ as a function of d rupt /d max for different values of α. Interestingly, up to insignificant statistical fluctuations, sZ is independent of α and a nearly linear function of d rupt /d max as that of σ p in Fig. 6 .

We see that the analytical model presented in this section correctly links the microstructure to the overall strength of the agglomerate. The trends are well predicted by the model up to the prefactor η, which appears to be weakly dependent on the material parameters. The physical interpretation of its value may be related to the presence of the large proportion of capillary bonds with nonzero gap, as briefly discussed previously. However, it may also be a consequence of the inhomogeneous stress transmission inside the agglomerate.

Conclusions

In this paper, we used a 3D particle dynamics algorithm together with a capillary force law to analyze the cohesive strength and microstructure of spherical agglomerates. The agglomerates were constructed by extraction of spherical samples from a granular bed prepared by compaction for different values of particle size span, and subjected to diametrical compression between two platens for different values of debonding distance, which accounts for the amount of liquid in the capillary bonds. Despite the irreversible nature of cohesive bonds (i.e. no new cohesive bond), we observe a plastic plateau before the onset of failure. We showed that the plastic strength is proportional to the characteristic capillary stress γ s / d with a multiplicative factor that is a linear function of the debonding distance, increases with the size ratio and is nearly independent of the friction coefficient (not shown here). We also introduced an analytical expression of the cohesive strength in terms of the packing fraction, wet coordination number, size polydispersity and debonding distance. This model is in excellent agreement with the observed trends up to a prefactor that we estimated from the numerical data, and which is weakly dependent on the material parameters.

As previously discussed in this paper, our results provide the behavior of an 'ideal' granule in the sense that the granules were not created by an agglomeration process. We presently work on the simulation of the granulation process in a rotating drum that will allow us to investigate the important complementary problem of predicting the agglomerate microstructure from the process parameters such as rotating speed and filling rate. The ideal granule can then be used as a reference system with which we will compare part of our results. Given the broad applications of the agglomeration process, it is also desirable to validate the simulations by comparison with experimental observations. Systematic diametrical compression tests are presently underway in order to determine the effects of material parameters on the granule strength allowing for the validation of numerical results. Our experiments are in good agreement with the order of magnitude of the cohesive strength of agglomerates for η 0.4. In these experiments, particles of the same typical size as in simulations were used. They were mixed with water and agglomerated into granules in a rotating drum. Preliminary experimental results were presented in [START_REF] Contreras | Cohesive strength of iron ore granules[END_REF] .

Fig. 1 .

 1 Fig. 1. Schematic drawing of the forces acting on particle i by a contacting particle j and by a non-contacting neighboring particle k .

Fig. 2 .

 2 Fig. 2. Schematic representation of diametrical compression test (a), and snapshots of granules composed from 50 0 0 primary spherical particles before diametrical compression for (b) α = 1 ; (c) α = 2 , and (d) α = 5 .

Fig. 3 .

 3 Fig. 3. Snapshot of a granule under diametrical compression for α = 1 . The lines are force chains joining the particle centers.

Fig. 4 .

 4 Fig. 4. Vertical strength σ zz normalized by the reference stress σ c for different val- ues of friction coefficient μ, size ratio α and the debonding distance d rupt , as a func- tion of cumulative vertical strain ε.

Fig. 5 .

 5 Fig. 5. Normalized plastic strength of the granule for several values of the debonding distance as a function of the size ratio α. The inset shows the non-normalized value of the strength as a function of α.

Fig. 6 .

 6 Fig. 6. Normalized plastic strength σ p / σ c of the granules for different values of the size ratio as function of debonding distance d rupt .

Fig. 7 .

 7 Fig. 7. The initial wet coordination number Z 0 of the granules for different values of debonding distance as a function of the size ratio α.

Fig. 8 .

 8 Fig. 8. The prefactor η in Eq. (23) as a function of the debonding distance d rupt for different values of α.

Fig. 9 .

 9 Fig. 9. The Z s in Eq. (23) as a function of the debonding distance d rupt for different values of α.


Kaixiang Zhang
email: kaixiangzhang@zju.edu.cn

Jian Chen
email: jchen@zju.edu.cn

Franc ¸ois Chaumette
email: francois.chaumette@inria.fr.

F Chaumette

Visual Servoing With Trifocal Tensor

In this paper, a trifocal tensor-based approach is developed for 6 degrees-of-freedom visual servoing. The trifocal tensor model among the current, desired, and initial views is introduced to describe the geometric relationship. Then, the tensor elements are refined to construct the visual feedback without resorting to explicit estimation of the camera pose. Based on the extracted tensor features, an adaptive controller is designed to drive the camera to a desired pose and compensate for the unknown distance scale factor. Moreover, Lyapunovbased techniques are exploited to analyze the system stability and convergence domain. Simulation results are provided to demonstrate the effectiveness of the developed approach.

I. INTRODUCTION

Since closing the control loop with vision sensors can increase the flexibility and accuracy of a robotic system, many efforts have been devoted to visual servoing in the past few decades [START_REF] Hutchinson | A tutorial on visual servo control[END_REF]- [START_REF] Dame | Mutual information-based visual servoing[END_REF]. To develop vision-based control strategies, a good choice is to use the geometric correlation among multiple views.

The two-view geometry, such as homography [START_REF] Malis | 2-1/2-D visual servoing[END_REF]- [START_REF] Silveira | Direct visual servoing: Vision-based estimation and control using only nonmetric information[END_REF] and epipolar geometry [START_REF] Mariottini | Image-based visual servoing for nonholonomic mobile robots using epipolar geometry[END_REF], [START_REF] Becerra | A sliding-modecontrol law for mobile robots based on epipolar visual servoing from three views[END_REF], has been widely used in visual servoing. More precisely, both homography and epipolar based methods construct the geometric relationship between the current and desired views to facilitate the control development. The geometric relationship is formulated by homography matrix or fundamental matrix, which can be calculated through the corresponding feature points in different views. However, both homography and epipolar based methods have drawbacks. The decomposition of homography matrix requires some knowledge about the desired pose to determine the unique solution, while the epipolar geometry becomes ill-conditioned with short baseline and with planar scenes [START_REF] López-Nicolás | Visual control through the trifocal tensor for nonholonomic robots[END_REF].

Different from homography and epipolar geometry, trifocal tensor encapsulates the intrinsic geometric correlation among three views and is independent of the observed scene. Due to this fact, the trifocal tensor has great potential in addressing visual servoing [START_REF] Andreff | Laser steering using virtual trifocal visual servoing[END_REF]. Most of the existing trifocal tensor based methods focus on controlling a nonholonomic mobile robot to achieve different tasks, mainly including regulation [START_REF] López-Nicolás | Visual control through the trifocal tensor for nonholonomic robots[END_REF], [START_REF] Becerra | Omnidirectional visual control of mobile robots based on the 1d trifocal tensor[END_REF], [START_REF] Becerra | Exploiting the trifocal tensor in dynamic pose estimation for visual control[END_REF], path following [START_REF] Sabatta | Vision-based path following using the 1D trifocal tensor[END_REF], and trajectory tracking [START_REF] Chen | Trifocal tensor-based adaptive visual trajectory tracking control of mobile robots[END_REF]. In [START_REF] López-Nicolás | Visual control through the trifocal tensor for nonholonomic robots[END_REF], by considering the planar motion constraint, a subset of the trifocal tensor is utilized to regulate a mobile robot toward a desired pose. In [START_REF] Becerra | Omnidirectional visual control of mobile robots based on the 1d trifocal tensor[END_REF], a twostep switching strategy is proposed with 1D trifocal tensor. Moreover, in [START_REF] Becerra | Exploiting the trifocal tensor in dynamic pose estimation for visual control[END_REF], the measurements of the trifocal tensor are exploited to estimate the pose of a mobile robot, and then the regulation task is accomplished with the estimated pose. Except for nonholonomic mobile robots, the trifocal tensor is rarely extended to address the 6 degrees-of-freedom (DOF) visual servoing. In [START_REF] Shademan | Three-view uncalibrated visual servoing[END_REF], all elements of the trifocal tensor are used as visual features to design an uncalibrated control scheme for manipulators. The redundant feature information is exploited to estimate the interaction matrix numerically, and thus it is difficult to ensure the system stability theoretically.

In this paper, a 6 DOF visual servoing approach is presented to regulate a camera to a desired pose. The geometric relationship among the current, desired, and initial views is described by the scene-independent trifocal tensor. To obtain a set of visual features with satisfactory decoupling properties, an auxiliary tensor variable is introduced. After that, 3 elements of the trifocal tensor and 6 elements of the auxiliary tensor variable are chosen based on the geometric connotation to define the visual features. The analytical form of the interaction matrix is derived to relate the control inputs with the variations of the tensor features. Furthermore, considering the unknown scale factor, an adaptive controller is developed via Lyapunov-based techniques. This paper is closely related to the works [START_REF] López-Nicolás | Visual control through the trifocal tensor for nonholonomic robots[END_REF], [START_REF] Chen | Trifocal tensor-based adaptive visual trajectory tracking control of mobile robots[END_REF], [START_REF] Shademan | Three-view uncalibrated visual servoing[END_REF] with distinct differences. First, this paper focuses on the 6 DOF eye-inhand visual servoing, while [START_REF] López-Nicolás | Visual control through the trifocal tensor for nonholonomic robots[END_REF], [START_REF] Chen | Trifocal tensor-based adaptive visual trajectory tracking control of mobile robots[END_REF] consider the 3 DOF vision-based control of mobile robots. Second, in our work, the tensor elements are selected based on the geometric relationship to facilitate the controller design and stability analysis. However, in [START_REF] López-Nicolás | Visual control through the trifocal tensor for nonholonomic robots[END_REF], the tensor elements are chosen experimentally. In [START_REF] Shademan | Three-view uncalibrated visual servoing[END_REF], all the trifocal tensor elements are used in the control development resulting in that the error system is cumbersome. Third, instead of estimating the interaction matrix on-line as [START_REF] Shademan | Three-view uncalibrated visual servoing[END_REF], the analytical form of the interaction matrix is derived in the paper, and the theoretical analysis of the system stability is presented.

The remainder of this paper is organized as follows. In Section II, the vision system is modeled using trifocal tensorbased techniques. The adaptive controller is designed in Section III, and the stability analysis is developed in Section IV. Furthermore, simulation results and conclusions are given in Section V and Section VI, respectively.

II. VISION SYSTEM MODEL

A. Problem Statement and Notations

As illustrated in Fig. 1, F c , F d , and F i denote the current, desired, and initial coordinate frames of the camera, respectively. Given the desired image captured in F d , the objective is to develop an adaptive controller to ensure that the current camera frame F c asymptotically converges from F i to F d using trifocal tensor-based techniques.

Some notations are introduced to improve the readability of this paper. Denote 0 n×n ∈ R n×n and 0 n ∈ R n as the nby-n zero matrix and n-by-1 zero vector, respectively. [•] × ∈ R 3×3 is the skew symmetric matrix associated to a 3-by-1 vector, and [•] ×(j) is the j-th column of [•] × . Given a vector c ∈ R n , c (j) ∈ R denotes the j-th element of c. Given a matrix C ∈ R n×n , C (j) ∈ R n is the j-th column of C and C (kj) ∈ R is the element on the k-th row, j-th column of C. For a trifocal tensor variable C ∈ R 3×3×3 , it can be seen as a collection of three matrices C (1) , C (2) , C (3) ∈ R 3×3 . Denote C (j) ∈ R 3×3 as the j-th matrix of C. Then, C (jl) ∈ R 3 is the l-th column of C (j) and C (jkl) ∈ R is the element on the k-th row, l-th column of C (j) . Moreover, a trifocal tensor variable, or matrix, or vector accompanied with a bracket (t) implies that its value varies with time.

B. Trifocal Tensor Model

As also shown in Fig. 1, to relate the camera frames, let c R i (t) ∈ SO 3 and c t i (t) ∈ R 3 be the rotation and translation between F c and F i expressed in F c . Likewise, the relative rotation and translation between F d and F i are denoted as d R i ∈ SO 3 and d t i ∈ R 3 , which are expressed in F d . Consider a static feature point O in the scene, its corresponding normalized Cartesian coordinates in the views F c , F d , and F i are denoted as m c (t), m d , m i ∈ R 3 , respectively. Let T (t) ∈ R 3×3×3 be the trifocal tensor among the current, desired, and initial views. Then the geometric relationship of the point correspondences m c (t), m d , and m i can be described as follows [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]:

[m c ] ×   3 j=1 m i(j) T (j)   [m d ] × = 0 3×3 . (1)
By using the relative pose information among three views, the trifocal tensor can be formulated into the following form [START_REF] López-Nicolás | Visual control through the trifocal tensor for nonholonomic robots[END_REF], [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]:

T (j) = c R i(j) d t T i -c t i d R T i(j) . (2)
Then, the expression of T (jkl) (t) ∈ R (j, k, l = 1, 2, 3) can be derived from [START_REF] Chaumette | Visual servo control part I: Basic approaches[END_REF], which is given by

T (jkl) = c R i(kj) d t i(l) -c t i(k) d R i(lj) . (3)
From ( 2) and (3), it can be seen that if T (t) is available, then the camera pose can be extracted from the trifocal tensor by using singular value decomposition (SVD). Although using the explicit pose information as feedback signals can simplify the controller design, SVD-based pose extraction is complicated and sensitive to image noises. To avoid the aforementioned problem, in this paper, we focus on utilizing the elements selected from the tensor variables to define the visual feedback. An intuitive idea is to define the feature signals with the elements of T (t). However, based on (3), it can be found that the time-varying pose information c R i (t) and c t i (t) are coupled into the expression of T (jkl) (t). Therefore, if the trifocal tensor elements are directly chosen as the visual features, the derived interaction matrix will not present satisfactory decoupling characteristics, which increases the complexity of controller design and stability analysis. Considering this issue, an auxiliary tensor variable

Q(t) ∈ R 3×3×3 is constructed to separate the time-varying signals c R i (t) and c t i (t). More precisely, Q(t) is designed as Q (j) = T (j) d R i(j) × = c R i(j) d t T i d R i(j) × (4)
where [START_REF] Chaumette | Visual servo control part I: Basic approaches[END_REF]

and d R T i(j) d R i(j)

×

= 0 are used. The auxiliary tensor variable Q(t) is constructed with the aid of the trifocal tensor T (t) and the rotation information d R i between the desired and initial views. Since these two views are recorded before starting the control task, d R i can be obtained offline with high accuracy by utilizing different vision-based techniques [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF], [START_REF] Ma | An Invitation to 3-D Vision: From Images to Geometric Models[END_REF], and we will show in Section V that our controller is robust to coarse approximation of d R i .

Besides, based on ( 2) and ( 4), it can be found that if the distance between the desired frame F d and the initial frame F i is zero, i.e., d t i = 0, then the tensor variables T (t) and Q(t) will not contain any terms related to the rotation matrix c R i (t). It is clear that under this circumstance, the rotation error cannot be eliminated with the trifocal tensor among the current, desired, and initial views, and an efficient way to address this issue is to construct the vision system with homography model. To ensure that the trifocal tensor model is applicable for the visual servoing, it is assumed that d t i = 0 in the following development.

C. Tensor Normalization

Based on (1), the trifocal tensor T (t) can be estimated up to a scale from point correspondences among three views [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF], i.e., T λ (t) = λT (t), where T λ (t) ∈ R 3×3×3 is the obtained scaled trifocal tensor and λ ∈ R is the scale parameter. Then, the scaled auxiliary tensor variable

Q λ (t) ∈ R 3×3×3 can be calculated by Q λ(j) (t) = T λ(j) (t) d R i(j) × , implying that Q λ (t) = λQ(t)
. Since λ is different each time the tensor variables are estimated, a normalization method should be introduced to ensure that the tensor variables are scaled by a common factor during the control procedure [START_REF] López-Nicolás | Visual control through the trifocal tensor for nonholonomic robots[END_REF]. According to (4) and the facts that

c R T i (t) c R i (t) = c R i (t) c R T i (t) = d R T i d R i = d R i d R T i = I 3×3 , and Q λ (t) = λQ(t), it can be determined that 3 j=1 3 k=1 3 l=1 Q 2 λ(jkl) = 2λ 2 d t T i d t i = 2λ 2 d * 2 . ( 5
)
where

d * d t T i d t i ∈ R is the constant distance between F d and F i .
Owing to the relationship shown in [START_REF] Dame | Mutual information-based visual servoing[END_REF], the normalized tensor variables T (t), Q(t) ∈ R 3×3×3 can be calculated by

T(j) = T λ(j) 1 2 3 j=1 3 k=1 3 l=1 Q 2 λ(jkl) Q(j) = Q λ(j) 1 2 3 j=1 3 k=1 3 l=1 Q 2 λ(jkl) . (6)
Using ( 5), ( 6), T λ (t) = λT (t), and Q λ (t) = λQ(t), T (t) and Q(t) can be rewritten as T(j) (t) =

T (j) (t) d * and Q(j) (t) = Q (j) (t) d * .
Moreover, according to ( 2) and ( 4), the normalized tensor variables T (t) and Q(t) can be formulated in terms of pose information and distance scale factor as follows:

T(j) = 1 d * c R i(j) d t T i -c t i d R T i(j) Q(j) = 1 d * c R i(j) d t T i d R i(j) × . (7)
To define the visual errors, the desired normalized tensor variables Td , Qd ∈ R 3×3×3 corresponding to T (t) and Q(t) need to be introduced. Obviously, these two desired normalized tensor variables are defined in terms of the desired pose, initial pose, and current pose being equal to the desired pose. Hence, Td and Qd can be computed from the desired and initial images before the control task starts. Besides, according to the definition of desired normalized tensor variables, Td and Qd can be expressed with the desired pose information and distance scale factor, as follows:

Td(j) = 1 d * d R i(j) d t T i -d t i d R T i(j) Qd(j) = 1 d * d R i(j) d t T i d R i(j) × . (8)
In next section, the tensor elements are selected to design an adaptive controller.

III. CONTROL DEVELOPMENT

A. Tensor Derivation

Before presenting the control strategy, the tensor derivation needs to be developed. The motion dynamics of the pose signals shown in Fig. 1 can be expressed as [START_REF] Chaumette | Visual servo control part I: Basic approaches[END_REF] c Ṙi = -

[ω] × c R i , c ṫi = -v + [ c t i ] × ω, d Ṙi = 0 3×3 , d ṫi = 0 3 (9)
with v(t), ω(t) ∈ R 3 being the linear and angular velocities of the camera, respectively. Based on ( 7), [START_REF] Silveira | Direct visual servoing: Vision-based estimation and control using only nonmetric information[END_REF], and the fact that

[ c t i (t)] × ω(t) = -[ω(t)] × c t i (t)
, the time derivative of T(j) (t) and Q(j) (t) can be deduced as follows:

Ṫ(j) = 1 d * v d R T i(j) -[ω] × T(j) Q(j) = -[ω] × Q(j) . (10)

B. Open-Loop Error System

To accomplish the visual servoing task, 3 elements of T (t) and 6 elements of Q(t) are chosen to construct system errors. Specifically, for the normalized trifocal tensor, T(jl) (t) with d R i(lj) = 0 is selected to define the translation errors, and for the normalized auxiliary tensor variable, Q(1l) (t) and Q(2l) (t) with QT (1l) (t) Q(1l) (t) = 0 and QT (2l) (t) Q(2l) (t) = 0 1 are selected to define the rotation errors.

Without loss of generality, to clearly show why this selection criterion is applicable and to facilitate the controller design, let us consider for instance that

d R i(11) = 0 and QT ( 11
) (t) Q(11) (t), QT (21)
(t) Q(21) (t) = 0 (i.e., d t T i d R i(1) ×(1) , d t T i d R i(2) ×(1)
= 0), and hence T(11

) (t) = T(111) (t) T(121) (t) T(131) (t) T , Q(11) (t) = Q(111) (t) Q(121) (t) Q(131) (t) T ,
According to [START_REF] Malis | 2 1/2 D visual servoing with respect to unknown objects through a new estimation scheme of camera displacement[END_REF], [START_REF] López-Nicolás | Homography-based control scheme for mobile robots with nonholonomic and field-of-view constraints[END_REF], and ( 11), e T (t) and e Q (t) can be rewritten as

e T = 1 d * d t i(1) c R i(1) -d R i(1) -d R i(11) c t i -d t i e Q = 1 d *    d t T i d R i(1) ×(1) c R i(1) -d R i(1) d t T i d R i(2) ×(1) c R i(2) -d R i(2)    .
(12) From ( 12), it can be concluded that if e T (t), e Q (t) → 0, then

c t i (t) → d t i and c R i(1) (t) c R i(2) (t) → d R i(1) d R i(2) provided that d R i(11) = 0, d t T i d R i(1) ×(1)
= 0, and

d t T i d R i(2) ×(1)
= 0. Note that the third column of a 1 It can be determined from [START_REF] Malis | 2 1/2 D visual servoing with respect to unknown objects through a new estimation scheme of camera displacement[END_REF] that

QT (1l) (t) Q(1l) (t) = 1 d * d t T i d R i(1) ×(l) 2 and QT (2l) (t) Q(2l) (t) = 1 d * d t T i d R i(2) ×(l) 2 . Therefore, QT (1l) (t) Q(1l) (t) = 0 and QT (2l) (t) Q(2l) (t) = 0 indicate that d t T i d R i(1) ×(l) = 0 and d t T i d R i(2) ×(l) = 0.
rotation matrix can be represented by the cross-product of the other two, and thus

c R i(1) (t) c R i(2) (t) → d R i(1) d R i(2) indicates that c R i(3) (t) → d R i(3) , i.e., c R i (t) → d R i .
Based on the above analysis, it can be seen that the selection criterion of the tensor elements ensures that the constructed system errors are applicable for the visual servoing.

Taking the time derivative of ( 12) and substituting from [START_REF] Mariottini | Image-based visual servoing for nonholonomic mobile robots using epipolar geometry[END_REF], the open-loop error system can be deduced, as follows:

ėT = d R i(11) d * v -[ω] × T(11) ėQ = L Q ω ( 13
)
where

L Q (t) ∈ R 6×3 is given by L Q Q(11) × Q(21) × . (14)

C. Controller Design

Based on the structure of the open-loop error system given in [START_REF] Andreff | Laser steering using virtual trifocal visual servoing[END_REF], the control inputs v(t) and ω(t) are designed as follows:

v = 1 d R i(11) -k 1 e T + d * [ω] × T(11) ω = -k 2 L + Q e Q (15) where k 1 , k 2 ∈ R are positive constant gains, L + Q (t) L T Q (t)L Q (t) -1 L T Q (t) ∈ R 3×6 is the pseudo-inverse of L Q (t)
, and d * (t) ∈ R is the estimate of the unknown distance scale factor d * . Note that a property will be given in Section IV to show that L T Q (t)L Q (t) -1 is symmetric and positive definite, and thus the calculation of L + Q (t) is always feasible. Moreover, to compensate for the unknown distance information, the update laws for d * (t) is given by

ḋ * = -k 3 e T T [ω] × T(11) (16)
with k 3 ∈ R being a positive constant gain. After substituting (15) into ( 13), the closed-loop error system can be derived, which is given by ėT = -

1 d * k 1 e T - d * d * [ω] × T(11) ėQ = -k 2 L Q L + Q e Q (17) where d * (t) d * -d * (t) ∈ R is the estimate error of d * .

IV. STABILITY ANALYSIS

To facilitate the stability analysis for the proposed approach, a property about the interaction matrix L Q (t) is presented firstly.

Property 1:

The matrix L T Q (t)L Q (t)
-1 is symmetric and positive definite.

Proof: According to ( 7) and ( 14), L Q (t) can be rewritten as

L Q = 1 d *    d t T i d R i(1) ×(1) c R i(1) × d t T i d R i(2) ×(1) c R i(2) ×    . (18)
Then, it can be derived that

L T Q L Q = 1 d * 2 d t T i d R i(1)
×( 1)

2 c R i(1) T × c R i(1) × + d t T i d R i(2) ×(1) 2 c R i(2) T × c R i(2) × = 1 d * 2 -d t T i d R i(1)
×( 1)

2 c R i(1) 2 × -d t T i d R i(2)
×( 1)

2 c R i(2) 2 × . ( 19
) Since c R i(1) (t) and c R i(2) (t) are linearly independent, d t T i d R i(1) ×(1)
= 0, and

d t T i d R i(2) ×(1) = 0, it can be concluded from (19) that L T Q (t)L Q (t)
is symmetric and positive definite [START_REF] Tayebi | Inertial vector measurements based velocity-free attitude stabilization[END_REF]. Thus, the corresponding inverse matrix

L T Q (t)L Q (t)
-1 is also symmetric and positive definite.

Theorem 1: Consider the system (13) under the control inputs [START_REF] Becerra | Exploiting the trifocal tensor in dynamic pose estimation for visual control[END_REF] and the update law [START_REF] Sabatta | Vision-based path following using the 1D trifocal tensor[END_REF]. Then,

1) The equilibria of the closed-loop system ( 17) are given by

Ω 1 : e T e Q = 0 9 , Ω 2 : e T e Q = -2   0 3 Qd(11) Qd(21)   , Ω 3 : e T e Q = -2   0 3 Qd(11) 0 3   , Ω 4 : e T e Q = -2   0 3 0 3 Qd(21)   . (20)
2) The equilibrium point Ω 1 is asymptotically stable provided that

e T Q (0)e Q (0) < 4 min {a 1 , a 2 } ( 21
)
where the positive constants a 1 , a 2 ∈ R are defined as

a 1 1 d * d t T i d R i(1) ×(1) 2 a 2 1 d * d t T i d R i(2) ×(1) 2 . ( 22
)
Proof: To prove Theorem 1, a non-negative Lyapunov function V (t) ∈ R is defined as follows:

V 1 2 d * e T T e T + 1 2 e T Q e Q + 1 2k 3 d * 2 . ( 23
)
After taking the time derivative of ( 23) and exploiting ( 16), [START_REF] Chen | Trifocal tensor-based adaptive visual trajectory tracking control of mobile robots[END_REF], and Property 1, it can be concluded that

V = -k 1 e T T e T -k 2 e T Q L Q L + Q e Q = -k 1 e T T e T -k 2 e T Q L T Q L Q -1 e Q ≤ 0 (24)
where e Q (t) L T Q (t)e Q (t) ∈ R 3 . Based on ( 23) and ( 24), it can be concluded that e T (t), e Q (t), d * ∈ L ∞ and e T (t), e Q (t) ∈ L 2 . Then, standard signal chasing arguments can be used to obtain that ėT (t), ė Q (t) ∈ L ∞ . As e T (t), e Q (t) ∈ L 2 and ėT (t), ė Q (t) ∈ L ∞ , Barbalat's lemma [START_REF] Khalil | Nonlinear Systems[END_REF] can be exploited to infer that lim t→∞ e T (t), e Q (t) = 0.

From Property 1, it can be found that the rank of L T Q (t) is three, and there exists a null space in such a way that ∀u ∈ ker L T Q (t) , L T Q (t)u = 0. Therefore, lim t→∞ e Q (t) = lim t→∞ L T Q (t)e Q (t) = 0 does not indicate that lim t→∞ e Q (t) = 0, i.e., there may exist multiple equilibria for the closed-loop system [START_REF] Chen | Trifocal tensor-based adaptive visual trajectory tracking control of mobile robots[END_REF]. It is clear that Ω 1 is one of the equilibria. In the following, we prove that Ω 2 , Ω 3 , and Ω 4 are also equilibria.

According to [START_REF] Becerra | Omnidirectional visual control of mobile robots based on the 1d trifocal tensor[END_REF], a basis of null space of L T Q (t) is given by

u 1 = Q(11) 0 3 , u 2 = 0 3 Q(21) , u 3 = Q(21) Q(11) . (25) If e Q (t) ∈ ker L T Q (t)
, then there exist constants λ 1 , λ 2 , λ 3 ∈ R, which cannot be zero simultaneously, such that

e Q = λ 1 u 1 + λ 2 u 2 + λ 3 u 3 . (26)
Based on ( 11) and ( 25), ( 26) can be rewritten as Qd( 11)

Qd(21) = - (λ 1 -1) Q(11) + λ 3 Q(21) (λ 2 -1) Q(21) + λ 3 Q(11) . (27)
Furthermore, from ( 7), [START_REF] López-Nicolás | Homography-based control scheme for mobile robots with nonholonomic and field-of-view constraints[END_REF], and the properties of rotation matrix, the following expressions can be derived:

QT d(11) Qd(21) = 0, QT d(11) Qd(11) = a 1 , QT d(21) Qd(21) = a 2 (28) QT (11) Q(21) = 0, QT (11)
Q(11) = a 1 , QT (21)
Q(21) = a 2 (29)
where a 1 and a 2 are defined in [START_REF] Khalil | Nonlinear Systems[END_REF]. Substituting ( 27) into (28) and using (29) to collect similar terms, it can be determined that

(λ 1 -1) λ 3 a 1 + (λ 2 -1) λ 3 a 2 = 0, (λ 1 -1) 2 a 1 + λ 2 3 a 2 = a 1 , (λ 2 -1) 2 a 2 + λ 2 3 a 1 = a 2 . (30)
Utilizing the last two equations of (30) to eliminate the terms related to λ 3 , it can be obtained that

λ 2 1 -2λ 1 a 2 1 -λ 2 2 -2λ 2 a 2 2 = 0. (31)
If λ 3 = 0, then based on the first equation of ( 30) and (31), it can be concluded that λ 1 and λ 2 which satisfy the constraint do not exist provided that a 1 = a 2 . If λ 3 = 0, then it can be obtained from (30) that (λ 1 -1) 2 a 1 = a 1 and

(λ 2 -1) 2 a 2 = a 2 .
According to these two equations, it can be concluded that (λ 1 , λ 2 , λ 3 ) can be selected as (2, 2, 0), (2, 0, 0), or (0, 2, 0). Thus, there exist three particular camera poses that will lead to the degeneration of the control inputs, i.e., there exist another three equilibria for the closed-loop system. Specifically, based on (27), the corresponding values of Q(11) (t) and Q(21) (t) for these three equilibria (poses) are given by

Ω 2 : Q(11) Q(21) = -Qd(11) Qd (21)
,

Ω 3 : Q(11) Q(21) = -Qd(11) Qd (21)
,

Ω 4 : Q(11) Q(21) = Qd(11) -Qd(21) . (32)
Then, from ( 7), [START_REF] López-Nicolás | Homography-based control scheme for mobile robots with nonholonomic and field-of-view constraints[END_REF], and (32), the corresponding value of c R i (t) can be derived as follows:

Ω 2 : c R i = -d R i(1) -d R i(2) d R i(3)
,

Ω 3 : c R i = -d R i(1) d R i(2) -d R i(3) , Ω 4 : c R i = d R i(1) -d R i(2) -d R i(3) . (33)
From a geometric point of view, (33) indicates that the rotation of the current camera frame F c at Ω 2 , Ω 3 , and Ω 4 differs from the desired frame F d by 180 • of rotation about the axes d R i(3) , d R i(2) , and d R i(1) , respectively. Moreover, according to [START_REF] Becerra | A sliding-modecontrol law for mobile robots based on epipolar visual servoing from three views[END_REF] and (32), the equilibria Ω 2 , Ω 3 , and Ω 4 given in ( 20) can be obtained. Now, the proof of the second claim of the theorem is presented. Consider the following non-negative function:

V Q (t)
e T Q (t)e Q (t) ∈ R, whose time derivative along the trajectory of ( 17) is given by

VQ = -2k 2 e T Q L Q L + Q e Q = -2k 2 e T Q L T Q L Q -1 e Q ≤ 0.
(34) Using (34), it can be concluded that ∀t ≥ 0, V Q (t) ≤ V Q (0). Moreover, based on ( 20) and (28), the corresponding values of V Q (t) at the equilibria Ω 2 , Ω 3 , and Ω 4 can be calculated, as follows:

Ω 2 : V Q = 4 (a 1 + a 2 ) , Ω 3 : V Q = 4a 1 , Ω 4 : V Q = 4a 2 .
(35) From ( 35), it can be deduced that the equilibrium in the domain

Φ = e Q ∈ R 6 |V Q = e T Q e Q < 4 min {a 1 , a 2 } is nothing else but Ω 1 . Furthermore, based on V Q (t) ≤ V Q (0),
it can be concluded that if e Q (0) ∈ Φ, then ∀t ≥ 0, e Q (t) ∈ Φ. Therefore, with the above analysis, it can be derived that Ω 1 is asymptotically stable, i.e., lim t→∞ e T (t), e Q (t) = 0, provided that e T Q (0)e Q (0) < 4 min {a 1 , a 2 }.

V. SIMULATION RESULTS

Simulation studies are performed with the aid of the opensource ViSP library [START_REF] Marchand | ViSP for visual servoing: A generic software platform with a wide class of robot control skills[END_REF] to illustrate the performance of the proposed control strategy. Specifically, a simulated freefloating camera controlled in 6 DOF is used to capture visual information, and 9 non-coplanar points are extracted from the current, desired, and initial views for the trifocal tensor estimation. The 9 point correspondences across three views are used to generate an initial solution to the trifocal tensor. Then, we use geometric minimization algorithm to provide a geometrically valid tensor [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. After obtaining the trifocal tensor, the auxiliary tensor variable can be computed with [START_REF] Collewet | Photometric visual servoing[END_REF]. By using the normalization method proposed in Section II-C, T (t) and Q(t) can be obtained. In the simulation, the desired camera pose F d is chosen as (0, 0, -0.8, 0 • , 0 • , 0 • ) expressed in the inertial coordinate frame. The gain parameters are adjusted as k 1 = 0.6, k 2 = 1, k 3 = 0.3, and the initial value of d * (t) is chosen as d * (0) = 0.4(m). In the following, three cases are presented to test the proposed approach.

• Pure Translation Case: The first simulation considers the pure translation along the three axes. The initial pose of the camera is set as (-0.5, -0.3, -1.8, 0 • , 0 • , 0 • ). Meanwhile, the corresponding simulation results are shown in Fig. 2. For this case, only the system error e T (t) needs to be eliminated. That is because the tensor elements used for the error system are selected by considering the geometric connotation, which guarantees that e T (t) and e Q (t) correspond to the translation and rotation errors, respectively. Thanks to the decoupled system errors, the camera trajectory in Cartesian space is linear without redundant rotation motion. • Large Rotation Case: For visual servoing, one of the most challenging configurations is the large rotation error around the z-axis. To further evaluate the proposed approach, a 170 • z-axis rotation is considered here. More precisely, the initial camera pose is chosen as (0.29, 0.1, -1.4, 0 • , 0 • , -170 • ). The simulation results are depicted in Fig. 3. From 3(b), it can be seen that the image trajectories of the feature points follow a spiral motion, which is exactly as expected due to the rotational motion around the z-th axis. In fact, in this case, a 1 = 0.024, a 2 = 0.77, and the initial error e T Q (t)e Q (t) is e T Q (0)e Q (0) = 3.16 which does not satisfy the condition given in [START_REF] Tayebi | Inertial vector measurements based velocity-free attitude stabilization[END_REF]. Moreover, it can be found from (20) that the initial system error e Q (t) is very close to the equilibrium Ω 2 given in Theorem 1. Nevertheless, the system errors still converge to zero as shown in Fig. 3(c) without being restricted by Ω 2 .

• General Motion Case: The initial pose of the camera in this case is set as (-0.94, 0.32, -1.21, -20 • , 30 • , -50 • ), indicating that both the translation and rotation errors along the three axes exist. As illustrated in Fig. 4, the proposed approach can regulate the camera to the desired pose effectively. Note that the auxiliary tensor variable is constructed with the aid of the rotation matrix d R i . To test the robustness of the controller with respect to coarse d R i , an error is added (5 • on each axis). Moreover, a supplementary error is also added to the camera intrinsic parameters (10%). The simulation results with coarse d R i and camera intrinsic parameters are presented in Fig. 5. Due to the introduction of rotation error, the image trajectory of the feature points given in 5(b) is quite different from the one in 4(b). However, the convergence of the camera pose and image points coordinates to their desired values demonstrates the correct realization of the task.

VI. CONCLUSION

This paper presented a trifocal tensor-based approach for 6 DOF visual servoing. Partial tensor elements were selected based on the geometric connotation of trifocal tensor model to construct the visual feedback, which can avoid explicit camera pose decomposition. Considering the unknown distance scale factor, an adaptive controller was designed to drive the camera to the desired pose. The Lyapunov-based method was exploited to analyze the stability of the control system. Moreover, the performance of the proposed approach was evaluated from simulation results.

Fig. 1 .

 1 Fig. 1. Three-view vision model.

  and Q(21) (t) = Q(211) (t) Q(221) (t) Q(231) (t) T are utilized in the following development. Let e T (t) ∈ R 3 and e Q (t) ∈ R 6 be the system errors defined by e T T(11) -Td(11) e Q Q(11) -Qd(11) Q(21) -Qd(21) .

Fig. 2 .Fig. 3 .

 23 Fig. 2. Pure translation case: (a) Camera motion in Cartesian space. (b) Image trajectory of the feature points. (c) System error convergence. (d) Control inputs.

Fig. 4 .Fig. 5 .

 45 Fig. 4. General motion case: (a) Camera motion in Cartesian space. (b) Image trajectory of the feature points. (c) System error convergence. (d) Control inputs.

This work was supported in part by the National Natural Science Foundation of China under Grant 61433013, and in part by the scholarship from Zhejiang University.