D Grasso

D Borgogno

E Tassi

Numerical investigation of a three-dimensional four field model for collisionless magnetic reconnection

Keywords: Nonlinear dynamics, plasma instabilities, fluid instabilities, numerical simulations

In this paper we present the numerical investigation of a three-dimensional four field model for magnetic reconnection in collisionless regimes. The model describes the evolution of the magnetic flux and vorticity together with the perturbations of the parallel magnetic and velocity fields. We explored the different behavior of vorticity and current density structures in low and high β regimes, β being the ratio between the plasma and magnetic pressure. A detailed analysis of the velocity field advecting the relevant physical quantities is presented. We show that, as the reconnection process evolves, velocity layers develop and become more and more localized. The shear of these layers increases with time ending up with the occurrence of secondary instabilities of the Kelvin-Helmholtz type. We also show how the β parameter influences the different evolution of the current density structures, that preserve for longer time a laminar behavior at smaller β values. A qualitative explanation of the structures formation on the different z-sections is also presented.

Introduction

Magnetic reconnection is a fundamental process in highly conductive fluids and plasmas [START_REF] Priest | Magnetic Reconnection[END_REF][START_REF] Biskamp | Magnetic Reconnection in Plasmas[END_REF]. It can be defined as a change in the topology of the magnetic field lines, which decouple their motion from the fluid one. It is associated with a release of magnetic energy into heat, plasma kinetic energy and fast particle energy and it is characterized by the formation of current density sheets in the reconnection region along with strong velocity layers. The range of phenomena involving magnetic reconnection is very wide. It includes solar flares, geomagnetic substorms, interaction of the solar wind with the magnetopause, sawtooth oscillations and disruptions in laboratory plasma, such as in Tokamaks. One of the problems in magnetic reconnection is to identify the appropriate generalized Ohm's law and the physical mechanisms which cause the diffusion of the magnetic field through the plasma. In weakly collisional plasmas, such as the high temperature ones in Tokamaks, the inverse of the electron-ion collision frequency is larger than the relaxation time of internal sawtooth oscillations. This consideration made collisionless reconnection to become a frontier subject in the early '90 [START_REF] Aydemir | [END_REF]4]. In this regime the typical length scale of the reconnection process is given by the collisionless skin depth, d e . Although a kinetic approach should be invoked in order to treat such collisionless regimes, in the presence of a strong guide field, fluid models, which offer a computational advantage, are often used. Many studies have been done in the framework of two-fluid models. In this context a description of the plasma behavior can be made considering a simpler two-field description [5] or a more sophisticated four-field description [6]. In particular, the two-field model in [5] has been extensively analyzed in the last decade both in two-dimensional and three-dimensional configurations [7,8,9,10]. In this model the evolution of the magnetic flux and plasma stream function is followed assuming that variations of the magnetic field and of the plasma velocity, along the direction parallel to the guide field, are negligible. The fingerprint of this approach is the coupling between the evolution of the current density and vorticity fields, which evolve in ordered or turbulent structures depending on the value of the electron temperature, which enters the equation through the the ion sound Larmor radius ̺ s [7,9], which is a further characteristic scale length of the phenomenon. For values of ̺ s ≥ d e the velocity layers,advecting the current density and vorticity, evolve in ordered coherent structures aligned with the separatrix of the magnetic island. On the other hand, for ̺ s << d e the velocity layers tend to be aligned with the neutral line giving rise to strong shears that lead to the onset of secondary Kelvin-Helmholtz instabilities. This twofold picture has been also confirmed in three-dimensional configurations [10,11]. When we allow the system to develop also magnetic and plasma velocity perturbations along the direction parallel to the guide field, by considering the four-field model in [6], the picture becomes richer and a new scenario appears, where the evolutions of the current density and vorticity fields are no longer coupled. The crucial role in this change is played by the β parameter, expressing the ratio between the plasma and the magnetic field pressures. Indeed, when the low β limit is abandoned, the effects of the strong velocity shears that develop in the reconnection region, on one hand lead to a turbulent vorticity, but on the other hand they get suppressed in the evolution of the current density [11,12].

Recently, an extension of this four-field model to three dimensions has also been derived [13]. The model belongs to the class of fluid models for plasma physics, for which a noncanonical Hamiltonian formulation is known. The origin of this class can be traced back to the seminal work of Morrison and Greene [14] on ideal magnetohydrodynamics. The class was subsequently enlarged to a great extent by Phil Morrison and co-workers over the years, with the discovery of Hamiltonian structures for several reduced fluid models.

In this article we analyze this model by performing numerical simulations with the aim of investigating the evolution, and its dependence on β, of velocity and current density fields in a full 3D setting. In particular we address the question of understanding in what fields the secondary Kelvin-Helmholtz instability dominates and in what fields it is suppressed, and how such instability extends along the z direction. We also recall that the absence of an ignorable coordinate raises the problem of interpreting a dynamics which is much less constrained than the two-dimensional one [15,10,11]. In particular, the system no longer possesses the infinity of Casimir invariants which constrain the 2D case and in terms of which it is possible to explain the different structures observed in the current density and vorticity fields [5,8,16,12]. Such explanation was indeed possible, due to the presence of advected scalar fields whose existence was related to the presence of an infinite number of Casimirs. The paper is organized as follows: in Sec. 2 we present the model equations; in Sec. 3 the simulation results are shown and discussed; in Sec. 4 we focus on the energy partition; conclusions are drawn in Sec. 5.

Model equations

Our investigation is based on the 3D four-field model for collisionless reconnection, described in Ref. [13]. Considering a Cartesian coordinate system (x, y, z), with the constant magnetic guide field directed along z, the model equations are given by

∂(ψ -d 2 e ∇ 2 ⊥ ψ) ∂t + [ϕ, ψ -d 2 e ∇ 2 ⊥ ψ] -d β [ψ, Z] + ∂ϕ ∂z + d β ∂Z ∂z = 0, (1)
∂Z ∂t + [ϕ, Z] -c β [v, ψ] -d β [∇ 2 ⊥ ψ, ψ] -c β ∂v ∂z -d β ∂∇ 2 ⊥ ψ ∂z = 0, (2)
∂∇ 2 ⊥ ϕ ∂t + [ϕ, ∇ 2 ⊥ ϕ] + [∇ 2 ⊥ ψ, ψ] + ∂∇ 2 ⊥ ψ ∂z = 0, (3)
∂v ∂t + [ϕ, v] -c β [Z, ψ] -c β ∂Z ∂z = 0, ( 4
)
where ψ is the poloidal magnetic flux function, ϕ is the E × B stream function, Z is proportional to the parallel magnetic perturbation and v is the parallel plasma velocity. The parameter d e represents the electron skin depth, whereas the other two parameters c β and d β are defined by c β ≡ β/(1 + β), with β indicating the ratio between the plasma pressure and the toroidal magnetic pressure, and by d β ≡ d i c β , where d i is the ion skin depth. The Poisson bracket is defined by [f, g] = ẑ • (∇f × ∇g) for generic fields f and g, whereas the symbol ∇ ⊥ refers to the gradient perpendicular to ẑ. Equations ( 1)-( 4) are written in a dimensionless form, according to which, the time is normalized with respect to a characteristic Alfvèn time t A , lengths are normalized with respect to a characteristic scale L and magnetic fields with respect to a characteristic value B. If such value is taken to be the amplitude B 0 of the toroidal guide field, then the following ordering applies:

ψ ∼ ϕ ∼ Z ∼ v ∼ ∂ ∂t ∼ ∂ ∂z ∼ ǫ ≡ B p B 0 ≪ 1, (5)
∂ ∂x ∼ ∂ ∂y ∼ 1, (6)
where B p indicates a characteristic value for the poloidal magnetic field (remark that this ordering differs from that of Ref. [13], in which a normalization based on B p , instead of B 0 , was adopted). Note that the Eqs. ( 1)-( 4) are exact only at the order ǫ 2 . Indeed, in their derivation, only the contributions up to order ǫ of the normalized magnetic field

B = ẑ + ∇ψ × ẑ + c β Z ẑ + O(ǫ 2 ), (7)
have been used. A term -c β x 0 dx ′ ∂Z/∂z x , which is of order ǫ 2 , is required in order to have a divergence-free magnetic field, but it produces only higher order corrections to (1)- (4). By construction, the system (1)-( 4) is Hamiltonian and conserves the total energy

H = 1 2 D d 3 x (d 2 e J 2 + |∇ ⊥ ϕ| 2 + v 2 + |∇ ⊥ ψ| 2 + Z 2 ), (8)
where J = -∇ 2 ⊥ ψ is the parallel current density and D is the domain of interest. The corresponding Poisson bracket has been provided in Ref. [13].

Numerical results

The equations ( 1)-( 4) are integrated numerically by splitting all the fields in two parts: an equilibrium, independent on time, and an evolving perturbation. The perturbed component is advanced in time by an explicit, fourth order Adam-Bashforth scheme. The equations are solved in a three dimensional slab, with periodic boundary conditions along the y and z directions. Dirichlet conditions are applied at the edges of the x axis, imposing that all the perturbed fields go to zero. A compact finite difference [17] algorithm, suitable for non-equispaced grid, is used for the spatial operations along the x direction, while pseudo-spectral methods are adopted for the periodic directions. Numerical filters are introduced in order to control the numerical error propagation, which is a relevant issue in the absence of any dissipative effect, as in the problem we are considering. As described in [17], these filters smooth out the small spatial scales below a chosen cutoff, while leaving unchanged the large scale dynamics on all the time evolution of the process. Finally, in order to address the three-dimensional problem, the code has been parallelized adopting the MPI libraries. We set up a numerical experiment of spontaneous, 3D, collisionless reconnection process in a static equilibrium configuration with ψ eq = -log cosh(x) and ϕ eq = v eq = Z eq = 0.

(

) 9
The integration domain is defined by -L x < x < L x , -L y < y < L y and -L z < z < L z , where L x = 11.32, L y = 4π, L z = 64π. The choice of L x avoids any influence of boundary conditions imposed along the x direction on the reconnection dynamics. The value of L y allow to consider the magnetic reconnection instability induced by highly unstable modes, that in the twodimensional limit fall in the so-called large ∆ ′ regime. It is accepted [START_REF] Aydemir | [END_REF]4,18] that in this regime magnetic reconnection can develop on a relatively fast time scale. Finally, the choice of L z >> L y guarantees that k z = 2πn/L z < k y = 2πm/L y for all the pairs of mode numbers (m, n) that develop during the nonlinear evolution of the process, as imposed by the conditions ( 5)- (6).

In order to take into account the 3D effects on the magnetic reconnection we perturbed the equilibrium ( 9) by introducing the following double helicity perturbation on the current density field

δJ(x, y, z) = Ĵ1 (x) exp(ik y1 y + ik z1 z) + Ĵ2 (x) exp(ik y2 y -ik z2 z). ( 10
)
Here we considered the pairs of wave numbers (1, 1) and (1, -1), that correspond to the wave vectors components k y1 = 1/4, k z1 = 1/64 and k y2 = 1/4, k z2 = -1/64. These two helicities, that are both linearly unstable, have resonant surfaces at x s1 = 0.06258 and x s2 = -0.06258. Ĵ1 (x) and Ĵ2 (x) are functions localized within a width of the order d e , around x s1 and x s2 , respectively. Their typical amplitude is of order 10 -4 . In order to properly treat the small scale structures that typically generate in the collisionless regimes a mesh of nx = 801, ny = nz = 512 grid points has been adopted.

In this paper we analyzed the influence of the 3D effects on the dynamics of the magnetic reconnection instability by considering two sets of the physical parameters appearing in eqs. ( 1)-( 4). In both cases we have used d e = 0.24 and d β = 0.96, while c β is varied, assuming the values c β = 0.4 and c β = 0.8. This choice allows to address two regimes with different values of the ion skin depth d i and of the ratio, β, between the plasma and the magnetic field pressure. In particular, when c β = 0.4, β = 0.19, while c β = 0.8 gives β = 1.78, which we refer to as a small and a high β regimes.

Qualitative aspects of the dynamics along the z direction

As shown in the framework of a 3D, two-field, dissipationless model ( [15], [10], [11]), the nonlinear interaction between two initially imposed linearly unstable helicities strongly affects the spatial distribution of the small scale structures where current density and vorticity fields are located. The shape of these patterns is the consequence of the interaction between the current layers associated to the island chains corresponding to the imposed helicities. According to the 2D results [7], in the collisionless cases when ̺ s = 0 the current density field exhibits a cross-shape structure distributed along the separatrices of the magnetic islands, and is strongly peaked at the corresponding X points. For a single helicity with (k y , k z ) wave vector, the value of the y coordinate of the magnetic island chain X points, varies along z following the lines with slope k z /k y , while the value of the x coordinate does not depend on z. This remains valid also in presence of more than one helicity, at least at the very early stages of the nonlinear evolution. In the early phase, indeed, each helicity is characterized by a given number, depending on k y , of thin current density channels with a characteristic spatial orientation. Similarly to electric wires, the current density channels of various helicities interact among each other by means of attracting forces whose strength depends on their geometry. In particular, the maximum of the attraction is localized on the z = const sections where the centers of the current density channels, localized at the X-points of the corresponding magnetic islands, have the same y coordinates. On these sections the interaction leads to the merging of the current channels. An opposite scenario characterizes the z = const planes where the X and O-points of the magnetic islands with different helicities face each other. Due to the fact that negative minimum current density values lie on the O points, very small attracting or even repulsive interactions act on the maximum current density peaks, as it is the case when two electric wires with currents flowing in opposite directions, face each other. Fig. 1 sketches this behavior for a double helicity case with poloidal and toroidal mode number m = 1, n = 1 and m = 1, n = -1. The figure shows the current density distribution on four z = const sections at the times t 1 and t 2 , that correspond to the very beginning and the advanced phases of the nonlinear evolution of the process. In the planes z = -L z and z = 0, where the X points of the structures have the same y coordinate, the current layers experience the maximum attraction and merge by generating a single X point pattern at t = t 2 . On such planes, the dynamics is quite similar to that observed in the single-helicity case. On the contrary, on the planes z = L z /2 and z = -L z /2, where the maximum current density of an island chain faces the minimum of the second chain, current layers merge by generating a new spatial distribution where two X-points are still present. The above described situations refer to the two extreme cases that can take place. Through the z = const planes, lying between those described above, the dynamics varies continuously with z, with a force between the current layers, which depends on the relative shift between the positions of the X-points along the y direction.

Interpretation of the numerical simulations

Due to the high value of d β , the early stages of the nonlinear evolution of the reconnection process are similar in both the cases c β = 0.4 and c β = 0.8.

We define here the end of the linear phase by the time in which the growth of the two initially imposed modes, deviates from the exponential one, which they would have when considered separately. In this early nonlinear phase of the evolution, the contribution of the terms proportional to c β is almost negligible with respect to the Poisson brackets involving d β . Since, in this limit, the parallel magnetic field perturbation Z becomes proportional to the vorticity U , the equations for F = ψ -d 2 e ∇ 2 ⊥ ψ and for U decouple from the equations for Z and v and reduce to the two field limit of the model ( 1)-( 4), described in [15], where ρ s = d β [16]. According to the two field model results [11], also in this case the current density and vorticity field concentrate in analogous structures characterized by thin and sharp layers on all the z = const sections. Due to the proportionality with U , similar patterns develop also in Z, while the parallel velocity component v is characterized by large scale cell structures similarly to v ⊥ = |∇ ⊥ ϕ|, the field which it is mainly transported by. This behavior is illustrated in figure 2, where the contour plots of U , Z, J and v are plotted on the particular plane z = 0 at t = 90τ A for the c β = 0.8 case. When the influence of c β effects becomes important, at the later stages of the process, the vorticity layers tend to broaden, starting from the z = const sections where the X points of the linear magnetic islands associated with the different helicities have similar y coordinate. This leads to the formation of new vorticity structures, localized around the x=0 axis for the particular choice of symmetric initial conditions we assumed, and, as a consequence, of highly localized, bar shaped, perpendicular velocity layers similar to jets as shown in figure 3. As already pointed out in Ref. [13], when the influence of a finite β is not negligible, the structures of U and Z eventually differ. On the other hand, from Eqs. ( 2) and (3), one sees that in the limit c β → 0 (but with finite d β ) the relation Z = d β U satisfies the system, if the initial conditions on Z and U are the same, which is the case here. This explains the similarity in the structures of Z and U for negligible β. The amplitude of the perpendicular velocity tends to become stronger as the magnetic reconnection develops. In the advanced nonlinear phase it provides the dominant contribution both in the equations of the vorticity and the parallel velocity fields. U and v are simply transported along the perpendicular velocity and as a consequence they assume the same spatial distribution of v ⊥ . On the other hand, due to the value of d β we considered, the time evolution of J and Z continues to be dominated by the nonlinear terms involving the perpendicular component of the magnetic field,

d β [ψ, Z] and d β [∇ 2 ⊥ ψ, ψ]
respectively. In figure 4 the intense jets appearing on the z = 0 section for the perpendicular and parallel velocities are shown together with the more ordered structures characterizing Z and J.

The fronts of both the parallel and perpendicular velocity jets tend to move towards the center of the computational box starting from the boundary surfaces y = ±4π. Approaching the y = 0 plane, the tops of the jets roll up and assume a "mushroom cap" shape which makes the velocity layers very similar to the "finger" like structures typical of the hydrodynamic Rayleigh-Taylor instability. This suggests the fluid like behavior of the system in these structures. Despite the magnetic field line stochasticity which characterizes the 3D settings we treat here, in fact, magnetic field components parallel to the jets are practically zero, which makes the contribution of the Poisson brackets [∇ 2 ⊥ ψ, ψ] and [Z, ψ] locally negligible and reduces the U and v equations to purely fluid advection equations. This justifies the occurrence, at later times, when the width of the jets is of the order of d e and the velocity is comparable to v A , of a Kelvin-Helmholtz like secondary instability which is responsible for the break of the velocity patterns. At this stage c β effects appear also on J and Z fields producing a turbulent redistribution of the patterns in the regions where v ⊥ and v are picked, which coexists with the laminar filaments concentrated at the boundaries of the reconnection area. This behavior is clearly illustrated in fig. 5, where the perpendicular velocity and current density are shown on the z = 8π section at the latest stage of the process for the case c β = 0.8. On the x -y planes with a larger distance between the y coordinates of the linear magnetic island X-points the contribution of c β remains rather small all along the reconnection process. On these sections, in fact, the nonlinear coupling between the current density structures do not allow the formation of the intense and extremely localized velocity layers described above, as illustrated in fig. 6. Here a comparison of the merging of the current density layers is shown on two different planes. On the left is the contour plot as it appears at z = -64π and on the right as it appears at z = 32π plane. These planes correspond to the first and second sketch of fig. 1 respectively. We can see that the current density layers tend to broaden and their intensity is reduced as far as we consider z-planes where the distance between the y coordinates of the X-points, corresponding to the two initially imposed helicities, increases. Because of the dominant terms depending on d β , the U , J and Z fields exhibit analogous patterns, while v ⊥ and v are characterized by spatial structures larger than d e and amplitude much smaller than v A .

It is also interesting to comment on the difference between a small and a high β regimes, which we consider as represented by the c β = 0.4 and c β = 0.8 cases, respectively. In particular, in the high-β regime, one has d β ≈ d i and this parameter measures the strength of the Hall term in the equations governing the magnetic dynamics. In the 2D limit, the model becomes then identical to the Hall MHD model investigated in Ref. [19], apart from the difference in the mechanism breaking the frozen-in condition. The comparison between the c β = 0.4 and c β = 0.8 cases confirms that also in 3D geometry larger values of the β parameter accelerate the formation of the fluid jets, whose occurrence happens at earlier stages of the reconnection process. Parallel and perpendicular velocities inside the layers are typically smaller in the lower β case, for comparable extensions of the reconnection region on a given z = const section. This difference is highlighted in fig. 7, where the profiles of the velocity fields v ⊥ and v are shown on the section z = 0 for the two different values of c β mentioned above. In particular, on the first row v ⊥ (left frame) and v (right frame) are plotted as function of x, in the limited range -4 < x < 4 and -2 < x < 2 respectively, at the y coordinate corresponding to their maximum. We see that, when considered at times which correspond to the same phase of the evolution of the reconnection process in the two cases, the velocity fields are smaller by a factor of order two in the lower c β limit, drawn in red. This characteristic is seen also in the second row of fig. 7, where the velocity fields are plotted as function of y at x = 0, where the fluids jets are located. Here, we can also appreciate that the extension of the fluid jets, which start at the boundaries y = ±4π for both cases, is shorter when c β = 0.4. In this case (red lines) the jets end at y = ±5 while for c β = 0.8 (black lines) they end at y = ±1 approximately. This is due to the fact that the jets formation started earlier for the c β = 0.8. As it will be seen in the next section, greater outflow speeds in the parallel direction, at high β, can be ascribed to a transfer of energy, mediated by the Lorentz force, from the magnetic to the kinetic form. The amount of transferred energy is indeed proportional to c β . Therefore one can expect a greater acceleration at higher β. We remark also that, because the Hamiltonian ( 8) is independent on β, the total energy available in the small and high β case is the same. Therefore, increase of the parallel kinetic energy with β is a pure consequence of the dynamical transfer and is not influenced by differences in the total energy of the system. Finally, we notice that, as a consequence of having less steep velocity layers for small β, in the latter case the Kelvin-Helmholtz destabilization of the fluid jets is delayed, making its influence on the density current and parallel magnetic field perturbation rather smooth.

Energy considerations

Because the model under consideration is a Hamiltonian system, the total energy (8) should be preserved during the dynamics. It is therefore important to verify that the decay of H, due to numerical dissipation, remains negligible. From Fig. 8 it is possible to see that, at the latest time of the simulation, the drop in the total energy for both the cases reported is less than 4 . We remark also that, during most of the simulations, the velocity at which the magnetic energy decreases is much larger than that at which the total energy dissipates. This reassures us of the fact that the magnetic energy loss is due to a genuine reconnection (or also ideal) process and not to numerical dissipation. In addition to remarks concerning the reliability of the simulations, the analysis of the time evolution of the different forms of energy is also suitable for physical considerations. Fig. 8 shows how, evidently, the perpendicular magnetic energy is converted into various forms during the reconnection process. In order to better understand the mechanisms of transfer of energy, it is useful to consider the following relations, which can be easily obtained from the model equations ( 1)-(4):

d dt 1 2 d 3 x|∇ ⊥ ϕ| 2 = d 3 xJ(B • ∇)ϕ + O(ǫ 4 ), (11)
d dt

1 2 d 3 xv 2 = c β d 3 xv(B • ∇)Z + O(ǫ 4 ), (12)
d dt

1 2 d 3 x(|∇ ⊥ ψ| 2 + d 2 e J 2 ) = -d 3 xJ(B • ∇)ϕ -d β d 3 xJ(B • ∇)Z + O(ǫ 4 ), (13)
d dt

1 2 d 3 xZ 2 = -c β d 3 xv(B • ∇)Z +d β d 3 xJ(B • ∇)Z + O(ǫ 4 ). (14)
In order to obtain compact and physically more meaningful expressions, we express the right-hand sides of ( 11)-( 14) in terms of the magnetic field B, instead of ψ and Z. As a consequence, by virtue of what was specified in Sec.2, the resulting expressions are exact only at the leading order and corrections of order ǫ 4 are required. Of course, giving up B in favor of ψ and Z, yields the relations whose sum gives the exact conservation of the Hamiltonian H. Note also that all the leading order terms on the right-hand sides of ( 11)-( 14) possess the common form d 3 xf (B•∇)g, for generic fields f and g. Mechanisms that produce the transfer of energy from one form to another, can then become inefficient if some fields have a weak variation along the magnetic field lines. From Eqs. ( 11) and ( 13), we see that the mechanism that acts as a source for perpendicular kinetic energy, is a sink for the energy of the poloidal magnetic field and of the parallel current. The rise of E kperp observed in Fig. 8, then implies a drop in E B + E ke . The growth of the parallel kinetic energy, on the other hand, is possible only at finite β, in which case the parallel Lorentz accelerates the fluid. This mechanism acts also as a sink for the parallel magnetic energy. Note that, in Fig. 8, for the smaller c β case, there is a competition between the conversion of magnetic energy into parallel magnetic energy and perpendicular kinetic energy, which lasts until around t = 125 i.e. well into the nonlinear phase, when E Z reaches its maximum and then drops. At the same time, an increase in E v is observed. On the other hand, for the larger c β case, the conversion into perpendicular plasma kinetic energy, is dominant from the beginning of the nonlinear phase. We notice also that the form of energies that grow the most in the last phase of the simulations, i.e. E v and E kperp , are those determined by the fields, v and |∇ ⊥ ϕ|, respectively, which exhibit the secondary instability.

Conclusions

In this paper we analyzed a three-dimensional model for collisionless reconnection, which takes into account magnetic and velocity perturbations along the guide field direction. On the basis of previous two-dimensional studies we were interested in exploring in a 3D context the different behavior of the current density and vorticity structures. To this aim, we carried out a detailed analysis of the plasma velocity fields, which play a key role in the current density and vorticity dynamics. In the presence of high values of the d β parameter, we found that, during the early stages of the nonlinear phase, the perpendicular and parallel velocities on the z = const sections are characterized by large scale cell structures. The other fields J, U and Z, on the other hand, develop thinner spatial patterns similar to the case β = 0. Following the evolution process, we observe that the perpendicular velocity tends to form highly localized patterns, aligned along the x = 0 line, on some peculiar z = const sections, which depend on the particular choice of the initial perturbation. These patterns are then recaptured also in the parallel velocity and vorticity fields, while, due to the high value of the d β parameter chosen for this analysis, the evolution of the current density and of the Z field, is dominated by the contribution of the Poisson brackets involving the magnetic flux.

At the end of the reconnection process both the perpendicular and parallel velocity layers have become so localized and intense that they develop secondary instabilities of the Kelvin-Helmholtz type. In connection with this we find that the larger c β , the sooner these instabilities develop.

The increase of the velocities fields and the occurrence of secondary instabilities influence also the energy distribution, favoring the conversion of the perpendicular magnetic energy into perpendicular and parallel kinetic energy rather than into parallel magnetic energy. The latter, in particular is transferred into parallel kinetic energy through the action of the parallel Lorentz force. and the current density J (right bottom frame) on the section z = 8π at the time t = 125τ A , for the c β = 0.8 case. Jet like patterns are present in both the velocity components, whose maximum amplitude is of the order of the Alfvèn velocity, while Z and J exhibit a filamentary structure distributed on a large fraction of the corresponding area where the magnetic topology modifications are located.

Figure 5: Contour plots of the pependicular velocity v ⊥ (left frame) and the current density J (right frame) on the section z = 8π at the time t = 145τ A , for the c β = 0.8 case. The occurrence of a secondary Kelvin-Helmholtz instability is responsible for the break of the v ⊥ field layers. At this stage of the process, a turbulent redistribution also affects the J structures, especially along the layers where v ⊥ is concentrated. Figure 6: Contour plots of the current density on the z = -64π (left frame) and on the z = -32π (right frame) at the time t = 105τ A , for the c β = 0.4 case. The current density layers tend to broaden and their intensity is reduced as far as we consider z-planes where the distance between the y coordinates of the X-points, corresponding to the two initially imposed helicities, increases.

Figure 1 :

 1 Figure1: Sketch of the dynamics of the current density structures (black solid lines) in the nonlinear evolution of the 3D reconnection process. The two columns shows the characteristic contour plots of the current density at an early, t 1 , and at a more advanced stage, t 2 , of the nonlinear phase on four different z = const sections. In each plot the gray regions represent the corresponding domains where the magnetic topology modifications produced by the reconnection instability are located.

Figure 2 :

 2 Figure2: Contour plots of the vorticity U (left top frame), the parallel magnetic field perturbation Z (right top frame), the current density J (left bottom frame) and the parallel velocity v (right bottom frame) on the section z = 0 at the time t = 90τ A , corresponding to the early nonlinear evolution, for the c β = 0.8 case.

Figure 3 :

 3 Figure 3: Contour plots of the vorticity U (left frame) and the parallel velocity v (right frame) on the section z = 8π at the time t = 105τ A , for the c β = 0.8 case. Early fluid jets are located around the x = 0 line.

Figure 4 :

 4 Figure4: Contour plots of the perpendicular velocity v ⊥ (left top frame), the parallel velocity v (right top frame), the parallel magnetic field perturbation Z (left bottom frame) and the current density J (right bottom frame) on the section z = 8π at the time t = 125τ A , for the c β = 0.8 case. Jet like patterns are present in both the velocity components, whose maximum amplitude is of the order of the Alfvèn velocity, while Z and J exhibit a filamentary structure distributed on a large fraction of the corresponding area where the magnetic topology modifications are located.

Figure 7 :

 7 Figure 7: Profiles on the z = 0 section of v ⊥ and v as function of x (top left and right frame respectively) at the y coordinate corresponding to their maximum, and as function of y (bottom left and right frame respectively) at x = 0. The curves corresponding to the c β = 0.8 case are drawn in black, while the curves corresponding to the c β = 0.4 case are drawn in red.

Figure 8 :

 8 Figure 8: Normalized deviations, of the different energy contributions from their value at t = 0, as a function of time for two different simulations. On the left is the plot corresponding to the set of parameters c β = 0.4, d i = 2.4, while on the right is the plot corresponding to the case c β = 0.8, d i = 1.2. E V refers to the parallel plasma kinetic energy (1/2) d 3 xv 2 , E Z to the parallel magnetic energy (1/2) d 3 xZ 2 , E B to the perpendicular magnetic energy (1/2) d 3 x|∇ ⊥ ψ| 2 , E ke to the kinetic energy associated to the parallel current density (d 2 e /2) d 3 xJ 2 , E kperp to the perpendicular kinetic energy (1/2) d 3 x|∇ ⊥ φ| 2 and E tot to the total energy H.

Acknowledgments

This article is dedicated to P.J. Morrison, on the occasion of his 60th birthday. The authors thank Phil for many inspiring discussions had over the years. The authors acknowledge dr. L. Comisso for providing the sketch proposed in figure 1. E. T. acknowledges fruitful discussions with the Nonlinear Dynamics group at the Centre de Physique Théorique, Luminy. This work was supported by the European Community under the contracts of Association between EURATOM and ENEA and between EURATOM, CEA, and the French Research Federation for fusion studies . The views and opinions expressed herein do not necessarily reflect those of the European Commission. Financial support was also received from the Agence Nationale de la Recherche (ANR GYPSI).


Ana Cameirao
email: cameirao@emse.fr

Amara Fezoua

Yamina Ouabbas

Jean-Michel Herri

Myriam Darbouret

Philippe Glénat

AGGLOMERATION OF GAS HYDRATE IN A WATER-IN-OIL EMULSION: EXPERIMENTAL AND MODELING STUDIES

Keywords: gas hydrates, flow loop, rheology, chord length distribution, modeling

The formation of gas hydrates in water-in-oil emulsion is investigated at the laboratory pilot scale on a flow loop. The Archimede flow loop (30 m long, and 1 cm diameter) has been built to reproduce the thermodynamic conditions encountered in deep-sea pipelines with a maximum flow rate equal to 3x10 -5 m 3 .s

Recently, a Moineau pump has been added in order to maintain the flow at a constant rate, especially after the beginning of crystallization and independently of the rheological characteristics (the flow rate can be varied in the range 2x10 -5 m 3 .s -1 to 14x10 -5 m 3 .s -1 ). This work presents some preliminary experiments to characterize the flow properties. Also some crystallization experiments were realized in order to understand the coupling between crystallization and rheology by following the size of the agglomerates and the water conversion. When water droplets crystallize to form hydrate particles, an agglomeration phenomenon is evidenced by pressure drop measurements. The agglomeration phenomenon is also detected by a Focused Beam Reflectance Measurements (FBRM) probe and is highlighted by a sharp change in the mean chord length and a spread of the Chord Length Distribution (CLD) to larger chord length.

INTRODUCTION

The crystallization of gas hydrates in pipelines is a flow assurance problem that affects the petroleum industry in respect to the increase of the operative pressure and increase of the water cut. The main reason is that the remained off-shore petroleum resources are deeper and also because longer pipelines are needed.

To prevent the formation of gas hydrates three solutions are envisaged: using additives to prevent the formation, heating the pipeline and isolating it.

Because the two last are very expensive in energy and capital, adding additives is the most used solution. These additives can be of three different types: thermodynamic, kinetic or anti-agglomerant types [2,8].

In previous works, we modeled the agglomeration and it effect on the viscosity of the suspension in a batch reactor [3]. In a way to better understand the coupling between formation and flowing, a flow loop has been constructed at Ecole des Mines de Saint-Etienne. Using the flow loop with the recent configuration, the effect of emulsifier on the water-in-oil emulsions can be studied. This model can be used to understand the agglomeration degree of emulsion during crystallization by varying operating parameters: water cut, concentration of additives, concentration of salt in water…

The new configuration of the Archimede flow loop allows the study of the transport of slurries with high concentration of crystals.

The experimental results obtained in this work can permit a better understanding of crystallization under flow and contribute to a further work on modeling this type of crystallization by population balances as it was done before [4].

EXPERIMENTAL

APPARATUS AND PROCEDURE

Crystallization of water-in-oil emulsions of hydrate slurries under shear stress was studied in the Archimede loop. Details on the apparatus and results are given in [4]. The crystallization of gas hydrates is studied from the rheological properties deduced from pressure drop measurements, and from the size and shape of crystal deduced from the signal provided by the FBRM probe. The Lasentec D600L FBRM probe was used in this study to measure the chords length distributions during the crystallization process. The FBRM (Focused Beam Reflectance Measurement) probe is an in-situ particle size analyzer equipped with a rotating low intensity LASER beam which is reflected when it intercepts an object, see more details on [9]. For appropriate sampling of the particles, the probe is placed at the centre of a vertical section so that all categories of particles are expected to circulate in front of the probe window. It is also necessary to avoid forming a "dead zone" in front of the measurement window.

The ARCHIMEDE flow loop reproduces the thermodynamic conditions of the deep-sea pipelines with working pressures up to 100 bar. Temperature is controlled within the range 0 ºC to 14ºC. The flow loop is equipped with a FBRM probe but also with several temperature probes, differential pressure probes and a Coriolis flowmeter. Figure 1 represents the flow loop composed by a horizontal section of 30 meters and a riser of 12 m ending in a gas-liquid separator. The Moineau pump maintains the flow rate at a constant value between 2x10 -5 m 3 .s -1 and 14x10 -5 m 3 .s -1 .

A water in kerdane emulsion stabilized with an emulsifier is circulated in the flow loop. The organic phase is a C10-C13 distillation cut which is delivered by TOTAL Fluids under the commercial denomination of Kerdane. The additive is a specific surface-active solution, IPE202, patented by IFP. The NaCl is from Sigma Aldrich.

The formation of the emulsion (total volume of 11,5L) is generated by using an Ultrat Turrax impeller at 8800 rpm during 3 min. Then, the emulsion is introduced in the flow loop and circulated until the operative temperature of 4ºC is reached. The loop is then pressured up to 80 bar and Moineau pump (non-shearing progressive capacity pump from PCM) is run. Methane dissolves in the emulsion. The loop is then pressured again to 80 bars so that final pressure after solubilization is in the range 70-75 bar.

As the nucleation is observed, the pressure becomes de decrease again because of the gas consumption due to hydrate formation. The slurry is circulated in the flow loop until the end of crystallization and agglomeration.

At the end of the experiment, once the crystallization is considered as accomplished (pressure and temperature remains constants during hours), the slurry can be studied. A rheological characterization is done by varying the flow rate and measuring the pressure loss in the horizontal pipe. Different parametric studies were done by varying the percentage of water, concentration of IPE, flow rate and NaCl concentration (see Table 1).The main results are presented below.

EXPERIMENTAL RESULTS

Before studying the crystallization, the rheology of pure water, pure kerdane and emulsion of water in kerdane (30% volume water) have been studied at different temperatures by measuring the pressure loss at different flow rate. Figure 2 shows the results of this study. For each curves two different slopes were detected. This means that there are two different regimes of flow: laminar and turbulent. The transition between both regimes is approximate equal to values between 6x10 -5 and 8x10 -5 m 3 .s -1 .

At low flow rate, the flow can be considered as laminar and pressure loss can be calculated threw the following equations [START_REF] Sinquin | Rheological and Flow Properties of Gas Hydrate Suspensions[END_REF]:

∆ ൌ ݂ ఘ௩ మ (1) ܴ݁ ൌ ఘ௩ ఓ
(2)

݂ ൌ ଵ ோ
(3)

Dissolution of methane

After the injection of methane, the pressure (in this case the pressure divided by the initial pressure) decreases quickly due to the dissolution of methane in the emulsion. Figure 3 shows the decrease of pressure during the dissolution of 6 experiments at different conditions (see Table 1):

Figure 3. Dissolution of methane in water in kerdane emulsion.

The transfer rate for each experiment is showed in

Crystallization study

The crystallization experiments were analyzed in terms of conversion and average chord length during the crystallization process.

The conversion is calculated threw the equation PV=nZRT (4) where Z is the compressibility factor of the gas defined in the case of methane by the following equation Z = -0,0196P[MPa]+0,994(5)

[4].
The conversion is then calculated assuming that the cavities are filled at 100% with

iwater crw n n = η (4)
where the number of mole of water crystallized is equal to :

w t w 0 M ) w (% m M m n × = =
(5) The formation of methane hydrates can be described by the following reaction:

CH 4 (g) +5.75 H 2 O (l) → CH 4 5.75 H 2 O (s)
From this reaction the number of mole of water crystallized is equal to 5.75 mole of the consumed methane. (5) Table 1 show some interesting results about conversion.

Experiments at the same conditions but different flow rate show that when the flow rate increases the conversion increases (experiments 1 and 2) because the transfer rate increases (see Table 2).

For experiments at low water cut (experiment 4 and 2 and experiment 5 and 3) the conversion is complete or higher comparing with experiments with high water percentages in which the total conversion is not attained, probably because there is not enough methane to totally convert the water.

Comparing experiments at different percentages of additive (experiment 2 and 3), when the additive decreases the conversion decreases because the dissolution rate is smaller. Finally, when the emulsion is done with salted water (experiment 6), the conversion decreases (experiment 3). Comparing experiment 1 and 2 agglomeration is more pronounced in 1 at low flow rate than in 2, because the shear rate is lower.

Experiments 4, 2, 5 and 3 have large final agglomerates.

Experiments at different percentages of additive show, has expected, that decreasing the quantity of additive increases the agglomeration.

The use of salted water also seems to increase agglomeration when comparing experiment 3 and experiment 6.

From the work of LE BA [START_REF] Hung | Formation et agglomération de particules d'hydrate de gaz dans une émulsion eau dans huile : Etude expérimentale et modélisation[END_REF] by using the model of interpretation of the CLD, we can give an estimation of the number of primary particles (N p ) and fractal dimension of the agglomerates (D f ). the results for each experiment during time are shown in Table 3.

Comparing the first times of crystallization, when the additive percentage increases, the fractal dimension decreases, probably because agglomerates are more compact.

During each crystallization, the fractal dimension and the number of primary particles increases due to agglomeration.

When the percentage of water decreases or when the quantity of emulsifier decreases agglomerates are bigger (N p is larger) and viscosity higher (table 1) this was already observed in a preceding work [START_REF] Hung | Formation et agglomération de particules d'hydrate de gaz dans une émulsion eau dans huile : Etude expérimentale et modélisation[END_REF].

The experiments with salt show large agglomerates as the experiment with low concentration of additive. Nevertheless the viscosity of the suspension with salt is lower at the same conditions of water and additive (3 and 6 in table 1).

Comparing experiments 1 and 2, at high flow rate agglomeration is lower because of the shear rate. After the crystallization a rheological study is done. This study helps in understanding the capacity to transport the slurry in a pipeline as a function of the quantity of solids in the slurry.

(a) 0.1%IPE Experiments at low quantity of emulsifier have the transition between the laminar and the turbulent regime at higher flow rate.

Comparing Figure 5 and viscosity in table 1 experiment 5 has higher viscosity than 3, so the pressure loss with flow rate is more important to 5 than to 3.

CONCLUSIONS

An experimental study has been made in order to understand the crystallization of methane hydrates in a pipeline.

Experiments were made in Archimede flow loop with different flow rates using a Moineau pump, percentage of water in the emulsion, percentage of IPE and with or without NaCl.

For experiments at high percentage of water the conversion is not complete because there is not enough gas in the system to totally convert water and also because gas diffusion is more difficult. At low water percentage, conversion is complete and agglomeration is high. At high flow rates, the conversion is higher because the methane transfer rate from gas into oil is higher.

Experiments at high water cut and low percentage of IPE show a pronounced agglomeration because the frequency of collision between droplets increases and because there is a lot of water not converted that enhances agglomeration by consolidating aggregates and increases the final viscosity, as observed in preceding works [START_REF] Hung | Formation et agglomération de particules d'hydrate de gaz dans une émulsion eau dans huile : Etude expérimentale et modélisation[END_REF][START_REF] Colombel | Cristallisation et agglomération de particules d'hydrate de fréon dans une émulsion eau dans huile : étude expérimentale et modélisation[END_REF].

Figure 1 .

 1 Figure 1. Archimède flow loop.

Figure 2 .

 2 Pressure loss with varying temperature and flow rate.

Figure 4 .

 4 Figure 4. Average chord length during crystallization.

Figure 5 .

 5 Pressure loss with varying flow rate. From Figure 5 experiments were assembled by affinity: experiments at equal percentage of IPE have the same behavior in terms of pressure loss with flowrate.

  The first work showed that water droplets start to crystallize and then agglomerate in a fractal-like structure which can be postulated and evidenced from pressure drop measurements [4].

	In 2005, a Focused Beam Reflectance
	Measurements (FBRM) probe was added to the
	flow loop in order to detect agglomeration
	phenomenon. It was concluded that agglomeration
	was highlighted by a sharp change in the mean
	chord length and a spread of the Chord Length
	Distribution (CLD) to larger chord length. In order
	to interpret the CLD measurements, a model has
	been developed [1] based on literature survey
	[6,7].

This model

[1,[START_REF] Hung | Formation et agglomération de particules d'hydrate de gaz dans une émulsion eau dans huile : Etude expérimentale et modélisation[END_REF]

allows comparing the size of aggregates from viscosity measurements and from CLD. The model allows determining the fractal dimension of aggregates and their number of constituting particles. Results obtained from the model of Leba [1, 10] were compared to those given by viscosity measurements in order to better understand the rheological properties of emulsions.

Table 2

 2

		in s -1 .
	1	2	3	4	5	6
	0.0034 0.0208 0.0087 0.0106 0.0073 0.0092

Table 2 .

 2 Transfer rate. From Table 2, it is observed that experiments with increasing flow rate have increasing dissolution rate. When the percentage of water increases the transfer rate increases and when the percentage of additive increases the transfer rate decreases. Experiments with salted water show that the transfer rate is approximately the same.

Table 1 .

 1 Methane hydrate experiments.

	Experiment	Flow rate m 3 .s -1	Volume %water	Mass % IPE	g NaCl /l water	% conv	p eq (bar)	T f ºC	Final slurry Viscosity (Pa.s)
	1	7.5x10 -5	30	0.5	-	29	45.2	5	0.0074
	2		30	0.5	-	53	39.1	4.2	-
	3		30	0.1	-	39	42.2	4.4	0.0068
	4	13.5x10 -5	10	0.5	-	100	57.1	3	0.0052
	5		50	0.1	-	32	35.4	3	0.0096
	6		30	0.1	35	13	47.7	4.3	0.0059

Table 3 .

 3 Results from the interpretation of CLD during crystallization.

			1	2	3	4	5	6
	tini cr	Np Df	176 2.0	290 2.1	338 2.2	979 2.1-2.3	275 2.2	690 1.8
	tini cr +	Np	383	385	716	1004	717	655
	100min	Df	2.3	2.1	2.7	2.5	2.4	2.3
	tini cr +	Np	383	386	1127	1004	717	1124
	200min	Df	2.3	2.3	2.3	2.5	2.7	2.5

  Sum A. K., Koh C.A., Sloan E.D., Fidel-Dufour A, Gruy F, Herri JM, Clathrate hydrates : from laboratory science to engineering practice, Ind. Eng. Chem. Res., 48, 2009, p.7457-7465 [3] Pic JS, Herri JM, Cournil M, Experimental influence of kinetic inhibitors on methane hydrate Dalmazzone D., Hamed N. , Dalmazzone C., DSC measurements and modelling of the kinetics of methane hydrate formation in water-in-oil emulsion, Chemical Engineering Science, 64, 2009, p.2020 -2026 [6] Turner D.J., Miller K.T., Sloan E.D., Direct conversion of water droplets to methane hydrate in crude oil, Chemical Engineering Science, 64, 2009, p.5066 -5072 [7] Turner D.J., Miller K.T., Sloan E.D., Methane hydrate formation an d an inward growing shell model in water-in-oil dispersions, Chemical Engineering Science, 64, 2009, p.3996 -4004

		140
		120
		100
	Pdiff (kPa)	60 80	1 4
		40
		20
		[1] Leba H, Cameirao A, Herri JM, et al., Chord length distributions measurements during crystallization and agglomeration of gas hydrate in a water-in-oil emulsion: Simulation and experimentation, CHEMICAL ENGINEERING SCIENCE, 65, 3, 2010, p.1185-1200 [2] particle size distribution during batch crystallization in water, CANADIAN JOURNAL OF CHEMICAL ENGINEERING, 79, 3, 2001, p. 374-383 [4] Fidel-Dufour A, Gruy F, Herri JM, Rheology of methane hydrate slurries during their crystallization in a water in dodecane emulsion under flowing, CHEMICAL ENGINEERING SCIENCE, 61, 2, 2006, p. 505-515 [5] 0 0 2 4 6 8 10 12 14 16 Flow rate (m 3 /s)x10 5


Etienne Parizet
email: etienne.parizet@insa-lyon.fr

Vincent Koehl

Application of free sorting tasks to sound quality experiments

Keywords: number : 43, 66, Lj

In many studies devoted to the sound quality of industrial products, a perceptual space is determined through dissimilarity judgements on pairs of stimuli. A drawback of this procedure is that it can be very time consuming if the number of stimuli is large. An alternative procedure consists in a free sorting of sounds: averaging individual results provides a set of data which are considered as indicators of dissimilarities and analyzed using a multi-dimensional scaling method. The validity of this alternative can be discussed, as the psychological processes involved in the two procedures are different.

This study compared these two approaches in a particular case (door closure sounds). In this specific case, it was observed that dissimilarities obtained from the two procedures can be different, the more so as sounds are dissimilar and these differences can lead to slightly different perceptual spaces. Nevertheless, a free sorting experiment is a reliable way of reducing the number of stimuli in a large set of sounds. It allows selecting some representative sounds and narrowing the set of sounds while keeping in the subset most of the timbre features. This provides a useful preliminary step to a paired-comparison experiment.

INTRODUCTION

In sound quality applications, the identification of important timbre features is a key factor. This knowledge allows the selection of sound metrics which can be used as input of a preference model. It also gives useful indications about the way the object should be modified in order to improve its sound quality. Numerous methods can be used to identify these attributes; a very common one consists in asking subjects to evaluate the dissimilarity between two sounds presented by pairs. After collecting the individual data, multidimensional techniques lead to a perceptual space, in which stimuli can be placed so that their relative distances are close to their perceived dissimilarities. The axes of this space can then be related to timbre features, as implicitly used by listeners when achieving the task.

These timbre features are inferred by the experimenter (e.g. by listening to sounds with very different coordinates along this axis, or by computing correlations between these coordinates and candidate sound metrics, as loudness, sharpness and so on). Such a procedure has been widely used in many applications: sounds from musical instruments (Kendal and Carterette, 1991), synthesized ones [START_REF] Mcadams | Perceptual scaling of synthesized musical timbres: common dimensions, specificities, and latent subject classes[END_REF], car sounds (Lemaître et al., 2004[START_REF] Parizet | Analysis of car door closure sound quality[END_REF], aircraft noise [START_REF] Barbot | Perceptual representation of aircraft sounds[END_REF] as examples. A major drawback of this procedure lies in the number of pairs to be presented to the listener, which is a square function of the number of sounds. This can make the experiment very long and tedious for the listener.

Another procedure consists in free sorting experiments. Subjects are asked to group stimuli together, on a similarity basis; they can make as many groups as they want to. This task is related to categorization process. But the relation between categorization and similarity evaluation is strong [START_REF] Thibaut | Similarity and categorization[END_REF], though some exceptions have been pointed out (see [START_REF] Rips | Categories and resemblance[END_REF], for an example). Individual results are co-occurrence matrices, made of 0 and 1 and indicating which sounds were grouped together by the subject. Averaging individual co-occurrence matrices provides a matrix which is considered as a dissimilarity one, and, as such, can be analyzed thanks to a multi-dimensional technique. Various examples can be given, in the field of visual perception (see [START_REF] Börner | Searching for the perfect match: A comparison of free sorting results for images by human subjects and by Latent Semantic Analysis[END_REF][START_REF] Brémaud | Perceptive free sorting and verbalization tasks with naive subjects: an alternative to descriptive mappings[END_REF], haptic perception (Picard et al., 2003 or Tiest and[START_REF] Tiest | Analysis of hatpic perception of materials by multidimensional scaling and physical measurements of roughness and compressibility[END_REF], food evaluation [START_REF] Falahee | Perceptual variation among drinking waters: the reliability of sorting and ranking data for multidimensional scaling[END_REF], or perception of everyday life situations [START_REF] Edwards | The structure of perceived qualities of situations[END_REF]. In the field of sound perception, this kind of analysis was realized in the case of environmental sounds [START_REF] Guastavino | Etude sémantique et acoustique de la perception des basses fréquences dans l'environnement sonore urbain[END_REF][START_REF] Dubois | The meaning of city noises: Investigating sound quality in Paris (France)[END_REF], everyday sounds [START_REF] Guyot | Etude de la catégorisation d'un corpus de bruits domestiques[END_REF][START_REF] Houix | Classification of everyday sounds: influence of the degree of sound source identification[END_REF], cars dashboard tapping sounds [START_REF] Montignies | Study of the perceptive space linked to dashboard tapping sounds[END_REF] or loudspeakers quality [START_REF] Lavandier | Perceptual and physical evaluation of differences among a large panel of loudspeakers[END_REF]. In spite of the advantage of such a procedure, a question should be kept in mind: as the task of listeners does not consist in evaluating dissimilarities between sounds, can averaged data be considered as dissimilarities? Will the perceptual space obtained from these data be close to the one obtained from real dissimilarities? For sound stimuli, only one published study compared results obtained from these two methods [START_REF] Bigand | Multidimensional scaling of emotional responses to music: the effect of musical expertise and of the duration of excerpts[END_REF]. This study used musical excerpts as stimuli and showed a good agreement between results. The goal of this paper is to present another example, using another kind of sounds (car door closure sounds), and to discuss how free sorting experiments can be used in typical sound quality applications. Three experiments have been conducted:

-in the first one, a set of 35 sounds was sorted by listeners; this allowed to split this set in 6 groups of similar sounds; -in the second experiment, 12 sounds extracted from this set (2 sounds for each group) were sorted, using the same procedure. It appeared that stimuli could still be clustered in 6 groups, indicating that the organization of sounds, as determined from the first experiment, was robust enough; -finally, in a third experiment, subjects were asked to evaluate dissimilarities between the same 12 sounds presented by pairs. The perceptual spaces obtained from the results of this experiment and the previous one proved to be slightly different.

The comparison of the results of these three experiments will provide some guidelines about the way a free sorting experiment can be used in such applications.

EXPERIMENT 1 : FREE SORTING OF 35 SOUND SAMPLES

Stimuli

Car door closure sounds were used in that experiment; they were part of a study presented in another paper [START_REF] Parizet | Analysis of car door closure sound quality[END_REF]. At the beginning of this study, 16 cars were used. Each car was entered in a semi-anechoic room and the driver's door closure sound was recorded with a dummy head (Bruel et Kjaer type 4133) placed outside the car, in the typical position of the driver leaving the car. Sounds were sampled at a 16-bits resolution with a 44.1 kHz sampling frequency. Some modifications of the door seals were realized on two cars, increasing the number of stimuli to 27. Finally, for 8 situations, a second recording of the closing was introduced in the set of sounds in order to check the repeatability of measurements. Thus the overall number of stimuli was 35. Their duration was of approximately one second.

Free sorting experiment

Sounds were presented to listeners through headphones in a quiet room. Subjects were informed they were listening to door closure sounds. They were asked to group together similar sounds, forming as many groups they wanted to. They could listen to sounds as many time as they felt necessary.

31 subjects participated to the experiment (students or staff member of the laboratory).

Results

For each listener, a 35 x 35 co-occurrence matrix was determined. Each cell of this matrix was 0 if the two corresponding sounds had been placed in the same group by the listener, 1 otherwise. First of all, the importance of inter-individual variability was evaluated. Using Rand Index as an indicator of agreement between two individual clusterings [START_REF] Hubert | Comparing partitions[END_REF], it was not possible to separate the panel into sub-groups of listeners.

Therefore, all individual co-occurrence matrices were averaged, which lead to a dissimilarity matrix. The values of this matrix ranged between 0 and 1. It should be noted that data thus obtained are real distances, as they fulfil the triangle inequality ( , ) ( , ) ( , )

d x z d x y d y z ≤ + .
A hierarchical clustering using Ward's method (Saporta, 2006, p. 258) was computed and lead to the dendrogram presented in figure 1.a. In this diagram, the vertical axis represents distance between stimuli: the higher the node between two stimuli, the larger the distance between them. It was decided to split up the set of sounds in six groups, represented on that figure. It can be noted that, when two recordings of the same door were included in the set of sounds, they belonged to the same group (these repetitions are labelled as (2-3), (7-8), (9-10), (11-12), (20-21), (25-26), (29-30) and (31-32)).

For the following of the study, two sounds were selected in each group. The selection rule was the following one : for each sound, the average dissimilarity between this sound and all the other ones belonging to the same group was computed. The two sounds leading to the minimum values were considered as the best representative ones of the group. These sounds are labelled as 18 and 30 (for group A), 3 and 6 (group B), 8 and 27 (group C), 13 and 16

(group D), 9 and 28 (group E), 26 and 32 (group F).

EXPERIMENT 2 : FREE SORTING OF 12 SOUND SAMPLES

Experimental procedure

The 12 sound samples selected from the previous experiment were used. 30 subjects participated to a second free sorting experiment. They had to group these 12 sounds in clusters, using the same procedure as in the first experiment (and the same playback conditions). They were members of the laboratory (staff and students). Only 8 of them had participated in the first experiment (which had taken place more than one year before that one).

Results

Individual results were collected and analyzed as in the previous experiment. Subgroups of listeners could not be built, so that the mean dissimilarity matrix was computed from all subjects. The dendrogram obtained from the hierarchical clustering analysis of this matrix is shown on figure 1.b. It can be seen that the set of sounds can be organized in the same 6 groups. Therefore, it seems that the structure underlying the set of sounds used in the first experiment was preserved by the selection of representative sounds.

: EXPERIMENT 3 : DISSIMILARITY RATINGS, 12 SOUND SAMPLES

Experimental procedure

The same set of 12 sounds was used in a dissimilarity experiment. After listening to each pair of sounds, the listener had to evaluate the dissimilarity between the two stimuli, and to give his answer by moving a cursor along a continuous scale labelled from "sounds are equal" to "sounds are extremely different". 40 people took part in this experiment. As required by the car supplier funding this part of the study, listeners did not belong to the laboratory and the jury was balanced in two ways. First of all, it was made of 19 women and 21 men. Then, in each gender group, half of the subjects were between 30 and 45 and the other ones between 46 and 60.

This experiment provided a set of individual dissimilarities, recorded as numbers varying between 0 and 1; there were 40 x 66 such values, as 66 pairs (representing the upper half of the 12x12 matrix) were presented to each listener.

Results

As in the previous experiment, no clustering of listeners could be made. Differences between groups of subjects (age and gender) were not significant, and a hierarchical clustering analysis did not allow building groups of subjects from individual results. Therefore, individual values were averaged to derive a mean dissimilarity matrix. A hierarchical clustering analysis was conducted for that matrix and the result is presented in figure 1.c.

Once again, the 6 group structure can clearly be identified: this seems to be an actual organization scheme of the sound samples and the selection of the representative sounds in experiment 1 allowed to preserve this scheme. Nevertheless, some discrepancies between dendrograms appear on figure 1. For instance, stimuli from the D group (labelled as 13 and 16) are closer to the A, E and F groups than to the C and B groups, while the reverse was true in the previous experiment. Both experiments lead to the same 6-groups clustering, but the overall structures of these groups are slightly different from one experiment to the other one.

In the following, the comparison between experiments 2 and 3 (based on the same set of sounds) will be made in more details.

COMPARISON OF EXPERIMENTS 2 AND 3

Comparison of dissimilarities

In each case, the averaged dissimilarities between sounds were computed. They are shown in figure 2. Each point represents one of the 66 pairs of sounds. The X-axis coordinate of the point is the distance between the two sounds as obtained from the free-sorting experiment and its Y-axis coordinate represents the dissimilarity computed from the individual direct pair wised estimations. A circle around a point indicates that the two values are not different in a significant way (p<0.05, Kruskal Wallis criterion). It can be seen that, when sounds are rather similar, both methods gave close results. On the contrary, in the case of rather different sounds, the distance obtained from the free-sorting procedure is larger than the value given by the direct evaluation. As an example, 17 pairs are composed of sounds which had never been grouped together by any listener. But the average distance between these sounds, as given by the direct evaluation, varied between 0.46 and 0.73. It seems as if there is a threshold: for a dissimilarity greater than a given value, sounds will be attributed to different groups by all listeners.

This can explain the discrepancies observed between the two dendrograms of figure 2 (b) and (c). Sounds were organized in the same 6 groups structures (as small dissimilarities are close between the two experiments) but the higher levels of the dendrograms were different.

Comparison of perceptual spaces

The two averaged dissimilarities matrices were analyzed using the MDSCAL multidimensional scaling method (Cox and Cox, 2001, p. 31). In both cases, keeping 6 singular values appeared to be a satisfying compromise (giving 90% of the sum of all singular values). In order to compare the results, the perceptual space obtained from free sorting data was adjusted to the one computed from the dissimilarity evaluation procedure using a procrustean transformation (rotation and dilatation of the space, Cox and Cox, 2001, p. 123). Results are shown on figure 3 (plane of axis 1 and 2): both procedures provided similar perceptual spaces, though some differences can be noted. For example, sounds 18

and 30 are clearly separated on the left diagram, which is no longer true on the right one.

Such differences could be due to the sampling of subjects (as different listeners participated to experiments 2 and 3). This hypothesis was tested using a bootstrap technique. First of all, a reference space was defined as the one obtained from the dissimilarity experiment, dissimilarities being averaged over the whole panel (40 listeners). Then, in the whole set of individual dissimilarities, 15 subjects were randomly selected, using replacement. For this sample of subjects, the averaged dissimilarity matrix was computed and a 6 dimensional perceptual space was determined from an MDSCAL analysis. This space was adjusted to the reference one using a procrustean analysis [START_REF] Cox | Multidimensional scaling[END_REF]. This operation was repeated 5000 times and the statistics of all recorded results could be computed (mean value and 95% confidence interval). Results can be seen on figure 4.a. The average position of a sound is indicated by a star. As expected, this position is close to the position in the reference space which is labelled as a number. The ellipse represents the 95% confidence location of this position.

The same computation was done using data obtained from the free sorting experiment and gave the results represented on figure 4.b. It can be seen that sounds coordinates are somewhat different from the reference ones. As an example, sound 27 is closer to sounds 3 and 6 than to sound 8, which was not true in the reference plane. For five sounds only (30, 3, 6, 8 and 28), the reference position lies inside the confidence ellipse.

From this analysis, it can be thought that the dissimilarities obtained from the free sorting experiment lead to a perceptual space slightly different from the one issued from the pairwised estimations.

DISCUSSION

In the case of this study, the averaged co-occurrence matrix did not provide a perfect alternative to direct dissimilarities. As presented in figure 2, the two kinds of data were close to each other for rather similar sounds only.

The key point is that, in a free sorting experiment, individual co-occurrence matrices only contain 0 and 1 values. Inter-individual variability is thus necessary to obtain an averaged matrix which can be considered as a dissimilarity one. Consider an extreme case in which stimuli are grouped in three families by all subjects. The averaged co-occurrence matrix will also contain 0 and 1 values only and the perceptual space will take the form of an equilateral triangle in a plane. In such a situation, direct dissimilarity estimation may give more accurate information, as it is possible that dissimilarities perceived between two of these groups will be lower than those between each of these two groups and the third one. This information will not appear in the results of the free-sorting task.

In some cases, such a small inter-individual variability can be expected: an example may be an experiment in which very different sound sources are used. It is probable that all listeners will group sounds together according to their sources. When the set of sounds is more homogeneous, as in most sound quality applications, more various individual sorting data are collected. It can then be expected that averaged data can be considered as dissimilarities, at least for the closest sounds only, as it has been shown in this study. The inter-individual variability occurring in a free-sorting experiment can be estimated from the corrected Rand Index [START_REF] Hubert | Comparing partitions[END_REF]. It could be interesting to investigate more in order to relate the accuracy of the perceptual space computed from a free sorting experiment to the range of Rand Indexes computed between each pair of subjects. That could allow fixing some rules, e.g. a maximum value of averaged Rand Index over which it cannot be recommended to derive a perceptual space from the co-occurrence matrix.

Nevertheless, this limitation in the use of free-sorting experiments should not hide a very useful application of this procedure. As confirmed by the first experiment of this study, it can be used to reduce the number of stimuli, in cases where this number is too high to conduct pair-wise dissimilarity ratings.

CONCLUSION

This study showed that a free sorting experiment can be very useful to reduce the number of sound samples while keeping the main sounds features in the reduced set. That can be of a great practical importance in sound quality applications.

It is also possible to directly build a perceptual space from the results of a free-sorting experiment. But, as noted here, this space can be slightly different from the one computed from direct evaluations of dissimilarities. It can be expected that the larger the dissimilarities between sounds, the larger these differences will be, as inter-individual variability is necessary to obtain almost continuous values when averaging free-sorting individual data.

This limitation should not hide one of the main advantages of a free sorting experiment: it can involve a much larger number of samples than a paired comparison experiment. This made it very suitable for industrial applications of studies related to sound quality.

For sound quality applications, a practical recommendation would be to include a large number of sounds in the first set of stimuli presented to a jury. Thanks to a free-sorting experiment, this set of sounds can be reduced without any loss of timbre features and the reduced set of sounds can be used in a pair-wise dissimilarity evaluation experiment in order to build a reliable perceptual space.

Collected Figure Captions

Figure 1 :

 1 Figure 1: dendrograms obtained from the three experiments. (a) : free sorting, 35 sounds.

Figure 2 :

 2 Figure 2: comparison of dissimilarities obtained from the free sorting experiment and the

Figure 3 :

 3 Figure 3: MDS mappings of stimuli. Bold letters: from dissimilarity ratings. Italic letters:

Figure 4 :Figure 1 :Figure 2 :

 412 Figure 4: MDS mappings of stimuli and application of the bootstrap technique. (a):

Figure 3 :

 3 Figure 3: MDS mappings of stimuli (axes 1 and 2). Bold letters: from dissimilarity ratings.

Figure 4 :

 4 Figure 4: MDS mappings of stimuli and application of the bootstrap technique. (a) :


Aurélie Galfré

Amara Fezoua

Yamina Ouabbas

Ana Cameirao

Jean Michel Herri
email: herri@emse.fr

CARBON DIOXIDE HYDRATES CRYSTALLISATION IN EMULSION

Keywords: gas hydrates, cyclopentane, emulsion, kinetic, thermodynamic

Greenhouse gases emissions, like carbon dioxide, have been identified as major sources responsible for global warming. To reduce carbon dioxide emissions, capture process can be envisaged to extract CO 2 from flue gases. In this work, the technology of CO 2 separation and capture by gas hydrates crystallization is investigated. The principal barriers in developing this process are the necessary high conditions of formation (P, T), low hydrates formation rate, and separation efficiency. A promoter additive, the cyclopentane (CP), is proposed to improve the quality of the capture both at the thermodynamic and kinetic levels. An experimental study looking at the influence of the CP content on the kinetic of gas hydrates formation is presented. CP has been studied as a dispersed phase (emulsion) or as a secondary and clearly distinguished second phase. The rate of crystallization (gas consumption rate) reveals to be very different according to the CP dispersion. The system is oriented to form pure gas hydrates in one case, or mix CP-Gas hydrates in the other case. Finally, several equilibrium points of CP-CO 2 hydrates are obtained. The results demonstrate that cyclopentane is a good additive which can remarkably decrease the gas hydrates formation pressure and enhances the formation rate.

INTRODUCTION

Greenhouse gases emissions, like carbon dioxide, have been identified as major sources responsible for global warming. Carbon dioxide is essentially produced by combustion of fossil energy, like in power generation plant. One possible way to reduce these emissions is to install an industrial process that can be integrated in the existing installations to extract CO 2 from flue gases. Several capture processes exist. We can mention absorption by MEA, membrane or physical adsorption.

As the flue gases are at low concentration in CO 2 (typically 5%-15 % in power plant) and at high flow rates (several cubic meters per second of CO 2 ), the capture technology needs to run at severe conditions and to minimize the cost of the process Innovative technology capture by gas hydrates crystallization can be an alternative process. For high concentrated mixtures of CO 2 in nitrogen (steel making plant, at atmospheric pressure), Duc and al. [START_REF] Duc | CO 2 capture by hydrates crystallization -A potential solution for gas emission of steelmaking industry[END_REF] showed that hydrates capture can be a competitive process. CO 2 trapping from gases mixture (CO 2 /N 2 ) can be effective in a multistage separation and is then economically possible [START_REF] Duc | CO 2 capture by hydrates crystallization -A potential solution for gas emission of steelmaking industry[END_REF].

Gas hydrates are ice-like crystalline compounds formed from a combination of water and suitable guest molecules (like CO 2 , N 2 ), under low temperature and high pressure [START_REF] Dendy | Clathrate Hydrates of Natural Gases. 3e éd[END_REF]. Pure gas hydrates can store a large quantity of gas (typically 160 volumes (STP) of methane per volume of hydrates [START_REF] Thiam | L'étude des conditions thermodynamiques et cinétiques du procédé de captage de CO 2 par formation d'hydrates de gaz Application au mélange CO 2 -CH 4[END_REF]).

To reduce the high pressure necessary to form gas hydrates, a small amount of additive can be used. In literature, cyclopentane (CP) is described as an excellent promoter [START_REF] Liang | Hydrate Equilibrium Conditions for Cyclopentane and a Quaternary Cyclopentane-Rich Mixture[END_REF][START_REF] Li | Clathrate Hydrate Capture of CO2 from Simulated Flue Gas with Cyclopentane/Water Emulsion[END_REF][START_REF] Tohidi | Equilibrium data and thermodynamic modelling of cyclopentane and neopentane hydrates[END_REF][START_REF] Zhang | Enhanced Kinetics of CO 2 Hydrate Formation under Static Conditions[END_REF][START_REF] Zhang | Equilibrium of Hydrogen + Cyclopentane and Carbon Dioxide + Cyclopentane Binary Hydrates †[END_REF][START_REF] Zhang | Thermodynamic analysis of hydrate-based precombustion capture of CO 2[END_REF]. It forms structures II hydrates at a temperature near 280 K and at atmospheric pressure [START_REF] Liang | Hydrate Equilibrium Conditions for Cyclopentane and a Quaternary Cyclopentane-Rich Mixture[END_REF]. Among the thermodynamic promoters that can be used, such as TBAB [START_REF] Thiam | L'étude des conditions thermodynamiques et cinétiques du procédé de captage de CO 2 par formation d'hydrates de gaz Application au mélange CO 2 -CH 4[END_REF], [START_REF] Herri | Gas hydrates equilibria for CO 2 -N 2 and CO 2 -CH 4 gas mixtures-Experimental studies and thermodynamic modelling[END_REF] (tetra-n-butyl ammonium bromide) or THF [START_REF] Linga | The clathrate hydrate process for post and pre-combustion capture of carbon dioxide[END_REF], [START_REF] Linga | The clathrate hydrate process for post and pre-combustion capture of carbon dioxide[END_REF][START_REF] Linga | Gas hydrate formation from hydrogen/carbon dioxide and nitrogen/carbon dioxide gas mixtures[END_REF][START_REF] Linga | A new apparatus to enhance the rate of gas hydrate formation: Application to capture of carbon dioxide[END_REF][START_REF] Kumar | Separation of carbon dioxide from flue gas ("precombustion capture") via hydrate crystallization[END_REF][START_REF] Kang | Hydrate phase equilibria of the guest mixtures containing CO2, N2 and tetrahydrofuran[END_REF] (tetrahydrofurane), CP presents the constraint to be immiscible with water. So, the equilibrium pressure of CO 2 /CP hydrates is independent of CP quantity if CP is in excess. From a kinetic point of view, CP has revealed to be a promoter (Zhang and al. [START_REF] Zhang | Enhanced Kinetics of CO 2 Hydrate Formation under Static Conditions[END_REF]).

In this work, the crystallization of CO 2 /CP hydrates is studied. The kinetic effect of CP on the crystallization step is evaluated and the thermodynamic behavior is investigated. Two systems (cyclopentane dispersed in a emulsion, or cyclopentane in a clear second phase) are studied to compare their effect on the hydrates formation rate [START_REF] Li | Clathrate Hydrate Capture of CO2 from Simulated Flue Gas with Cyclopentane/Water Emulsion[END_REF]. Several ratio of CP (1.625% to 15% v CP/v water) are tested. Also, the benefice of CP has been investigated in a (Pressure, Temperature) domain in which CP/CO 2 binary hydrates and CO 2 gas hydrates can occur simultaneously. Another (Pressure, Temperature) domain has been investigated in which only CP/CO 2 binary hydrates can crystallize.

EXPERIMENTAL PROCEDURE

Experimental set up

The experimental apparatus (figure 1) allows to investigate the thermodynamic equilibrium conditions of gas hydrates (pressure and temperature) and to study the kinetic effect of promoters. The experiments are performed in a stirred batch reactor (autoclave, 2.46 L) cooled with a double jacket. Two Rushton turbines mix both gaseous and liquid phases, so that liquid and gaseous phases are considered as homogeneous phases. Two windows (12x2cm, polycarbonate), mounted on both sides of the reactor, allow to detect the occurrence of hydrates phase by direct visual observation. The solution to be tested is poured in a pyrex cell, located in the autoclave, in which the pressure can be raised up to 100 bar. Two Pt100 sensors (Prosensor instrument, precision 0.1K) give the temperatures of the liquid and the gaseous phases. The pressure is measured with a pressure transducer (range: 0-10 MPa, Keller instrument, precision of 0.1 bar). The data acquisition unit (T, P) collects the (P,T) evolution during hydrates crystallization. The system is cooled or heated by a thermostat HUBER CC3-K6 with an accuracy of 0.1 K. The liquid can be injected into the pressurized reactor by using a HPLC pump PU-1587 from JASCO.

The gas phase is supplied using gas cylinders and transferred to the reactor by a set of valves. The composition of the gas phase is determined in line by using a gas chromatograph, after sampling by a ROLSI instrument. This tool collects a controlled volume of gas (1-5 µm 3 ) which is directly injected into the loop of the gas chromatograph (VARIAN model 450 GC). This system allows to control the purity of the gas (CO 2 ). This work presents a set of experimental data on CO 2 -CP hydrates system. Then, a volume of cyclopentane is added.

Preparation of cyclopentane in water emulsion

CP dispersed as an emulsion or CP and water in separated phases

Following the previous work of Li and Fan [START_REF] Li | Clathrate Hydrate Capture of CO2 from Simulated Flue Gas with Cyclopentane/Water Emulsion[END_REF], we have investigated the influence of CP as a dispersed phase, or as a clear and separated phase from water. After pouring IPE 202, water and CP in the Pyrex cell, the emulsion is obtained by using an Ultra-turax (IKA T50) at stirring rate of 8800 tr/min. Fine CP droplets of some micrometers are obtained (figure 2).

Figure 2 Microscope photograph of CP (2.50 % vol.) in emulsion.

Crystallization procedure

The Pyrex cell (containing the solution or the emulsion) is placed in the reactor and stirred at 400 rpm. The temperature of the thermostat is decreased to the crystallization temperature (T c =274.15 K) and maintained to this value during all the crystallization process. When the liquid temperature reaches the equilibrium temperature, the stirrer is stopped and the reactor is vacuumed and then flushed with CO 2 (three to four times).

Afterwards, the reactor is charged with CO 2 . The gas phase is sampled with the ROLSI© instrument and in-line analyzed by gas chromatography to check CO 2 purity. During the introduction of CO 2 in the reactor, the stirring is stopped. The gas dissolution is considered as null. Then, the stirrer is immediately switched on. The pressure decrease (figure 3) puts into evidence the dissolution of the gaseous component in liquids (water and cyclopentane).

After a few time, a sudden increase of the temperature (exothermic process) indicates the appearance of first crystals (nucleation). During the formation of the crystals, the pressure decreases due to the gas consumption to form hydrates. After a while, Pressure and Temperature reach equilibrium (end of crystallization).

Dissociation procedure

The gas hydrates dissociation is operated at constant volume and started by heating the reactor with steps of 1K/h. Each step of temperature is considered to be at equilibrium [START_REF] Zhang | Equilibrium of Hydrogen + Cyclopentane and Carbon Dioxide + Cyclopentane Binary Hydrates †[END_REF].

As the temperature increases, the pressure raises in proportion to the gas desorption. As Gas/liquid water/liquid CP/hydrates equilibrium is reached, a change in dP/dT is observed and gas hydrates dissociation starts. Pressure increases dramatically, until total gas hydrates disappearance. Then, a sudden decrease in dP/dT shows the end of dissociation. The value of the surface increases versus CP content (CP droplet size has observed to be independent on the CP quantity).

When the nucleation occurs, the large increase of the temperature conducts the system to be out of the pure CO 2 hydrates formation area. It means that the CP/CO 2 hydrates forms and heats the solution out of the CO 2 hydrates formation area. But as temperature decreases, it re-inters in the CO 2 hydrates formation area and allows a competition between the formations of the two types of hydrates.

The intensity of the exothermic peak is directly proportional to the initial quantity of CP, so we can say that the quantity of CP/CO 2 gas hydrates is also proportional to the initial quantity of CP. The respective quantities of CO 2 gas hydrates and CP/CO 2 hydrates can be evaluated during the dissociation step. Indeed, when hydrates dissociation occurs, the (Pressure,Temperature) coordinates follow the respective equilibrium curves of CO 2 hydrates and CO 2 -CP hydrates (figure 5). The respective quantities can be then calculated. . In the experiments with 10.00 % or 15.00 % vol. of CP, practically no CO 2 gas hydrates has been produced. All CO 2 has been consumed in the CP/CO 2 hydrates formation (figure 5).

It confirms the previous assumption: CP/CO 2 binary hydrates crystallization appears in the beginning of the process and after CO 2 gas hydrates crystallization can occur or not depending on the intensity of the heating. No exothermic peak for CO 2 gas hydrates is observed (the nucleation can not be located) but the coexistence of the two types of gas hydrates is confirmed by the dissociation curve.

This study is in agreement on the Li and Fan's [START_REF] Li | Clathrate Hydrate Capture of CO2 from Simulated Flue Gas with Cyclopentane/Water Emulsion[END_REF] and Zhang and al. [START_REF] Zhang | Equilibrium of Hydrogen + Cyclopentane and Carbon Dioxide + Cyclopentane Binary Hydrates †[END_REF] studies.

Figure 5 CP in emulsion Evolution of pressure and liquid temperature versus time during dissociation ◊: CP-CO 2 gas hydrates [START_REF] Liang | Hydrate Equilibrium Conditions for Cyclopentane and a Quaternary Cyclopentane-Rich Mixture[END_REF], [START_REF] Zhang | Equilibrium of Hydrogen + Cyclopentane and Carbon Dioxide + Cyclopentane Binary Hydrates †[END_REF], [START_REF] Mohammadi | Phase equilibria of clathrate hydrates of methyl cyclopentane, methyl cyclohexane, cyclopentane or cyclohexane+carbon dioxide[END_REF].

x: simulation from GasHyDyn software At a pressure of 23.4 bar (2.5% vol. CP), the supercooling for CO 2 hydrates is around 2 K whereas for CO 2 -CP gas hydrates, supercooling is equivalent to 14.5 K (figure 6). At the beginning of the crystallization, the CO 2 consumption is fast (table 3, value for r 5hours ) due to the CO 2 -CP gas hydrates crystallization, then the CO 2 consumption slows down (table 3, r moy ) for the following several causes.

First at all, the emulsion becomes a concentrated slurry during the crystallization process, and CO 2 gas transfer starts to be slowed down. Secondly, as the suspension becomes very dense, the gas hydrates tend to agglomerate at the periphery of the Pyrex cell ( [START_REF] Zhang | Enhanced Kinetics of CO 2 Hydrate Formation under Static Conditions[END_REF]) to form a crust which slows down the temperature transfer and generates a temperature gradient within the reactor: cooling becomes difficult at the center of the Pyrex cell where the remaining water could crystallize. Finally, the driving force (which is initially around 2K for CO 2 gas hydrates) decreases during the crystallization process because pressure decreases. (P,T) formation conditions join the CO 2 equilibrium curve and crystallization runs slowly (due to the low driving force) (except 15% vol. CP, where only CP-CO 2 gas hydrates were produced).

Consumption rates are calculated by the equation [START_REF] Duc | CO 2 capture by hydrates crystallization -A potential solution for gas emission of steelmaking industry[END_REF] with n CO2,n as the initial mol number of CO 2 and n CO2,5hours or n CO2,f as the final mole number of CO 2 .

t n n r f or hours CO i CO dis ∆ - = 5 , 2 , 2 (1)
Where CO 2 mole number is calculated by:

RT T P Z PV n G CO ) , ( 2 = (2)
The compressibility factor Z is calculated by the Soave-Redlich and Kwong (SRK) equation of state (EOS) with parameters from Danesh [START_REF] Danesh | PVT and phase behaviour of petroleum reservoir fluids[END_REF]. The reactor gas volume is given by:

L R G V V V - = (3)
Consumption rate (table 3) increases with the CP volume ratio in the emulsion. The crystallization strength is dependent of the CP droplets quantity. It means that the crystallization rate is not only dependent on the gas transfer rate at the gas/liquid interface, but also on the liquid/liquid transfer rate at the emulsion interface. When CP and water are in two different phases, CP has a different kinetic effect from the CP in emulsion system [START_REF] Li | Clathrate Hydrate Capture of CO2 from Simulated Flue Gas with Cyclopentane/Water Emulsion[END_REF]. Firstly, the CO 2 dissolution rate reveals to be independent on the CP volume. CP has a lower density than water. Due to the configuration of the system, CO 2 dissolves in CP at the CO 2 /CP interface and then CO 2 transfers from CP to water at the CP/water interface. The rate of dissolution becomes dependent on the interfacial resistances which are independent to the volume but only proportional to the surface.

The difference in time induction can not be explained clearly, not only because the induction time is often a stochastic phenomena, but also because of the process initial conditions. In the case of 5.00% vol CP, the CO 2 dissolution is started at a higher temperature. The achievement of CO 2 supersaturation takes then longer time.

Li and Fan [START_REF] Li | Clathrate Hydrate Capture of CO2 from Simulated Flue Gas with Cyclopentane/Water Emulsion[END_REF] have performed CO 2 /N 2 crystallization with CP in emulsion and CP in a separated phase. They have only crystallized CO 2 /N 2 /CP gas hydrates and we cannot completely compare their experiments to ours. Contrary to their study, at same or equivalent CP percentage, crystallization in the separated phase system seems faster (r 5hours , r moy ) than in the emulsion. Theses values can be explained by the crystallization behavior.

If CP-CO 2 gas hydrates can be formed, one mole of CO 2 can produce a higher quantity of solid in comparison of pure CO 2 gas hydrates. Indeed, a part of the cavities are occupied by CP and not by CO 2 .

During crystallization from CP in emulsion, a double effect is observed: because CP is well dispersed, the crystallization rate is enhanced, and because of the crystallization rate improvement, the temperature is heated (figure 6) in a domain in which pure CO 2 gas hydrates cannot be formed. CO 2 is then only consumed to form CP-CO 2 gas hydrates. The liquid water is quikly immobilized under the solid form. The solution becomes very viscous: mass and heat transfers fall down and gas consumption rate decreases.

During crystallization from CP in a clear and separated phase from water, the CP transfer is limited and the rate of crystallization is decreased. The exothermic peak is low (table 5) and the pure CO 2 gas hydrates can continue to form. One mole of CO 2 can not form a large volume of solid (in comparison to CP-CO 2 hydrates): the suspension remains well dispersed and gas transfer can continue to be efficient. In this case, the gas consumption rate maintains high to form essentially pure CO 2 gas hydrate: contrary to the previous case (CP in emulsion), the gas consumption rate decreases rapidly because it forms a huge quantity of CP-CO 2 hydrates which affects negatively the gas to liquid transfer rate.

CP thermodynamic effect

Fan and al. [START_REF] Liang | Hydrate Equilibrium Conditions for Cyclopentane and a Quaternary Cyclopentane-Rich Mixture[END_REF] are reported that the quadruple point (CP-CO 2 gas hydrates) is at 280.22 K and 19.8 kPa. Mohammadi and al. [START_REF] Mohammadi | Phase equilibria of clathrate hydrates of methyl cyclopentane, methyl cyclohexane, cyclopentane or cyclohexane+carbon dioxide[END_REF] are shown that cyclopentane is the strongest hydrates promoter with neopentane in the second place.

In this paragraph, the thermodynamic promotion effect of CP is investigated. The crystallization is only performed in the CP-CO 2 gas hydrates area (Table 6 and figure 7).

The crystallization occurrence is checked by the dissociation step. Some CP-CO 2 equilibrium points are obtained and allow to complete the literature data.

CONCLUSIONS

In this work, CP has been used to enhance the crystallization kinetic using two approaches.

When CP is dispersed into an emulsion, time induction is shorter (0.23 to 2.71 hours versus CP volume ratio) and exothermic peak (crystallization power) is higher (1.4 to 6.3K versus CP volume ratio). But gas consumption rate is affected by the strength of the process and becomes rapidly low (formation of CP-CO 2 hydrates, high concentrated slurry).

In the second case, CP is at the water surface (in different phases): time induction is around two hours, and pure CO 2 hydrates are suspected to be formed. The solid fraction is lower and don't affect the mass transfer. The consumption rate can maintain high.

In the both systems, two hydrates types (CP/CO 2 hydrates and CO 2 hydrates) have been produced. The relative ratio of the hydrates depends on the CP ratio volume, but also on the quality of the CP dispersion. A well dispersed emulsion will generate essentially CP/CO 2 hydrates, whereas pure CO 2 hydrates will be generated if the CP is badly dispersed. Finally, crystallization has been practiced only in the CP/CO 2 hydrates area. Gas hydrates have been obtained at a pressure lower than 1 bar. Equilibrium points of CP/CO 2 binary hydrates complete the equilibrium data from litterature.

Figure 1

 1 Figure 1 Experimental set up

Figure 3

 3 Figure 3 Crystallization step: Evolution of pressure and liquid temperature versus time.(CP in emulsion, CP = 2.50% vol)

Figure 3

 3 Figure 3 Dissociation step: Evolution of pressure and liquid temperature versus time (CP in emulsion, CP = 2.50% vol.)

Figure 6

 6 Figure 6 CP in emulsion Evolution of pressure and liquid temperature versus time during crystallisation ◊: CP-CO 2 gas hydrates [4], [8], [17], x: simulation from GasHyDyn software

Figure 7

 7 Figure 7 Experimental dissociation conditions for hydrates of carbon dioxide and cyclopentane. ∆: Literature data [4], [8], [17], x: our work

Table 1 reports

 1

	component	Purity	Supplier
	CO 2	>99.995%	Air Liquide
	C 5 H 10	>95%	Alfa Aesar
	IPE 202 a		IFP Energies
			nouvelles

the purity and suppliers of the different products used. All experiments are performed with purified water (Milli-Q plus 185).

Table 1 .

 1 Purities and suppliers of components IPE 202 is a surfactant and anti-agglomerate additive. It stabilizes the emulsion and limits the agglomeration of hydrates particles.

	IPE 202 (0.45% water mass) and water are
	weighed by a Mettler balance (precision d=0.01g).

Table 2 .

 2 Experimental data: CP in emulsion

	1.63	29.8	16.2	275.7
	2.50	29.2	16.6	275.8
	5.00	30.2	17.4	276.1
	10.00	29.9	16.0	275.2
	15.00	30.8	19.5	275.4

Figure 4 Liquid temperature evolutions versus time during the first five hours of the crystallization with different CP volume ratio (CP in emulsion)

CP %v/v water t ind (h)

	∆ ∆ ∆ ∆T
	(K)

r 5hours (mol/ h) x10 -2 r moy (mol/ h) x10 -2

	1.63	2.71	1.4	1.21	0.70
	2.50	1.24	4	1.90	1.05
	5.00	0.33	4.3	9.82	1.30
	10.00	0.23	5.7	13.70	1.13
	15.00	0.23	6.3	15.50	2.22

Table 3 Kinetic

 3

information: CP in emulsion Figure

4

shows the exothermic peak due to the crystallization. With a low CP concentration (1.63 % vol), crystals nucleation (which usually takes a long time and is random) occurs in a few hours (see values in table

3

) and seems dependent of the CP quantity. The intensity of the temperature peak (table

3

) seems proportional to the CP content. It suggests that the nucleation strength depends on the availability of the CP. It is also dependent of the total surface of exchange between CP and water.

Table 4

 4 Experimental data: CP/ water in different phases

	CP	t ind	∆ ∆ ∆ ∆T	r 5hours (mol/	r moy (mol/
	%v/v water	(h)	(K)	h)	h)
				x10 -2	x10 -2
	2.50	1.67	1.3	14.3	3.42
	5.00	1.98	0.3	8.08	3.00
	Table 5 Kinetic information: CP/ water in different
		phases

Table 6 .

 6 Experimental data: CP in water emulsion

	CP / % v/v water	P alim / bar	P eq /bar	T eq /K
	5.00	10.4	3.6	275.7
	20	8.4	0.5	277.8

ACKNOWLEDGEMENTS

This work has been supported by FUI in the framework of the ACCACIA project. The authors thank all the members of the GasHyDyn team for their constant support and especially the technical staff: Alain Lallemand, Fabien Chauvy, Richard Drogo and Albert Boyer.