Antonin Chapoy
email: antonin.chapoy@pet.hw.ac.uk

Rod Burgass

Bahman Tohidi

Martha Hajiw

Christophe Coquelet

Thermophysical Properties, Hydrate and Phase Behaviour Modelling in Acid Gas-Rich Systems

come   L'archive ouverte pluridisciplinaire

Introduction

As the global demand for natural gas is forecasted to steadily grow, the demand will be met by supplies produced from "unconventional" gases. With sour gas fields worldwide accounting for some 40% of natural gas reserves, the development and production of these reserves are under very serious consideration [START_REF] Lallemand | Solutions for the treatment of highly sour gases[END_REF]. The main challenges facing companies developing these fields with high concentrations of acid gases are reservoir engineering, phase behaviour predictions, processing and their removal from hydrocarbon streams, transportation and storage.

Due to the huge quantities of these gases produced and more stringent environmental regulations these gases cannot be flared, thus one of the most viable options is to inject them back into the reservoir for storage as well as enhance oil recovery (EOR). For acid gases, the disposal alternatives are: injection of compressed acid gas into the formation; disposal of acid gases with the formation water or solubilise acid gases into a light hydrocarbon solvent and use the solvent as a miscible flood EOR (Jamaluddin et al., 1996). In fields with high concentration of H 2 S dehydration using glycols may lead to additional problems such as corrosion and the release of H 2 S when it is regenerated. In such cases the knowledge of multicomponent gas water content and the optimum injection pressure will be essential. Furthermore reinjection into new reservoirs (such as Kashagan field in Kazakhstan) will require extremely high pressures. The design of such compressors requires accurate thermophysical properties of multicomponent mixtures.

However limited data are available on the phase and hydrate behaviours of CO 2 -rich or acid gas systems to validate existing thermodynamic models. Therefore the applicability of the existing models and their uncertainties can lead to over or undersized designs. In this communication we present experimental techniques, equipment and thermodynamic modelling for investigating systems with high CO 2 or H 2 S concentrations, including gas reservoirs with high CO 2 content and/or CO 2 -rich systems from capture processes. In particular, experimental measurements of the locus of incipient hydrate-liquid water-vapour curve for a methane -H 2 S binary system and a CO 2 -rich natural gas (70 mole % of CO 2 and 30 mole % of light hydrocarbons C 1 to nC 4 ) in equilibrium with liquid water are presented at pressures up to 35 MPa. Experimental data are reported for water content in equilibrium with hydrates at about 150 bar and temperature range from 233 to 283 K. Density and viscosity of the mixture were also measured from 253 to 423 K at pressure up to 124 MPa. An example of dry ice formation in a CO 2 -rich natural gas will be also described.

The Cubic-Plus-Association (CPA-EoS) or the Soave-Redlich-Kwong (SRK) equation of state combined with the solid solution theory of van der [START_REF] Van Der Waals | Clathrate solutions[END_REF] as developed by [START_REF] Parrish | Dissociation pressures of gas hydrates formed by gas mixtures[END_REF] was employed to model the fluid and hydrate phase equilibria as previously described by [START_REF] Chapoy | On the Phase Behaviour of the Carbon Dioxide -Water Systems at Low Temperatures: Experimental and Modelling[END_REF]. The predictions of the thermodynamic model were compared with the experimentally measured properties (saturation pressure, dew point, frost points, hydrates). A corresponding state model developed to predict viscosity of the CO 2 -rich stream [START_REF] Chapoy | Effect of impurities on thermophysical properties and phase behaviour of a CO 2 -rich system in CCS[END_REF]) was used to evaluate the new viscosity data.

Experimental Setups and Procedures

The majority of the setups and procedures used in this paper were described in detail in Chapoy et al. (2005), [START_REF] Chapoy | On the Phase Behaviour of the Carbon Dioxide -Water Systems at Low Temperatures: Experimental and Modelling[END_REF], [START_REF] Chapoy | Effect of impurities on thermophysical properties and phase behaviour of a CO 2 -rich system in CCS[END_REF] and [START_REF] Hajiw | Effect of acide gases on the methane hydrate stability zone[END_REF]. A brief description of each setup is given below.

Saturation and Dew Pressure Measurements and Procedures

The equilibrium setup consisted of a piston-type variable volume (maximum effective volume of 300 ml), titanium cylindrical pressure vessel with mixing ball, mounted on a horizontal pivot with associated stand for pneumatic controlled rocking mechanism. Rocking of the cell through 180 degrees at a constant rate, and the subsequent movement of the mixing ball within it, ensured adequate mixing of the cell fluids. Cell volume, hence pressure, can be adjusted by injecting/withdrawal of liquid behind the moving piston. The rig has a working temperature range of 210 to 370 K, with a maximum operating pressure of 69 MPa. System temperature is controlled by circulating coolant from a cryostat within a jacket surrounding the cell. The equilibrium cell and pipework were thoroughly insulated to ensure constant temperature. The temperature was measured and monitored by means of a PRT (Platinum Resistance Thermometers) located within the cooling jacket of the cell (accuracy of ±0.05 K). A Quartzdyne pressure transducer with an accuracy of ± 0.05 MPa was used to monitor pressure. The bubble point was determined by changing the volume of the cell and finding the break over point in the pressure vs. volume curve as shown in Figure 1.

A typical test to determine the dew point is as follows: To obtain the dew point using the isochoric method the cell is loaded with the test sample and is set to 5 degrees above the estimated dew point temperature. The cell is cooled until the system has clearly become two phase. The cell temperature is then step heated, allowing sufficient time for equilibration, until the system has clearly become single phase again. Throughout the process the cell is rocked using a pneumatic pivoting system to ensure all of the cell constituents are thoroughly mixed and equilibrium is reached. The system pressure and temperature are recorded every minute using a logging program. The recorded data is then processed to determine the system pressure at each temperature step. This process results in two traces with very different slopes on a pressure versus temperature (P/T) plot, one in the single phase and one in the 2 phases region. The point where these two traces intersect is taken as the dew point (Figure 2).

Hydrate dissociation Measurements and Procedures

Dissociation point measurements were conducted using the isochoric step-heating method. The phase equilibrium is achieved in a cylindrical cell made of Hastelloy equipped with a pressure magnetic mixer. A detailed description of the apparatus and test procedure can be found elsewhere [START_REF] Hajiw | Effect of acide gases on the methane hydrate stability zone[END_REF][START_REF] Chapoy | Effect of impurities on thermophysical properties and phase behaviour of a CO 2 -rich system in CCS[END_REF]. The weight of the fluids (i.e., water and the multicomponent fluid) injected are recorded prior to any measurements and the overall feed composition can thus be calculated.

A typical test to determine the dissociation point is as follows: the cell is cleaned and dried. About half of the volume of the cell is initially preloaded with water, prior to applying vacuum to the system. Then, the fluids are loaded into the cell to reach the first desired pressure the temperature is then increased stepwise, slowly enough to allow equilibrium to be achieved at each temperature step. At temperatures below the point of complete dissociation, gas is released from decomposing hydrates, giving a marked rise in the cell pressure with each temperature step (Figure 3). However, once the cell temperature has passed the final hydrate dissociation point, and all hydrates have disappeared from the system, a further rise in the temperature will result only in a relatively small pressure rise due to thermal expansion. This process results in two traces with very different slopes on a pressure versus temperature (P/T) plot; one before and one after the dissociation point. The point where these two traces intersect (i.e., an abrupt change in the slope of the P/T plot) is taken was the dissociation point (see Figure 3). For a full discussion on accuracy and uncertainties of hydrate dissociations measurements the reader is invited to check the work of [START_REF] Stringari | Study of factors influencing equilibrium and uncertainty in isochoric hydrate dissociation measurements[END_REF][START_REF] Hajiw | Etude des Conditions de Dissociation des Hydrates de Gaz en Présence de Gaz Acides / Hydrate Mitigation in Sour and Acid Gases[END_REF]. The equilibrium dissociation point is determined as being the intersection between the hydrate dissociation (pressure increase as a result of gas release due to temperature increase and hydrate dissociation, as well as thermal expansion) and the linear thermal expansion (no hydrate) curves.

Water Content Measurements and Procedures

The core of the equipment for water content measurement has been originally described by [START_REF] Chapoy | On the Phase Behaviour of the Carbon Dioxide -Water Systems at Low Temperatures: Experimental and Modelling[END_REF] and Burgass et al. (2014). The setup comprises of an equilibrium cell and a device for measuring the water content of equilibrated fluids passed from the cell. The equilibrium cell is similar to the one described in the saturation pressure measurements. The moisture/water content measurement set-up consists of a heated line, a Tuneable Diode Laser Adsorption Spectroscopy (TDLAS) from Yokogawa and a flow meter. The estimated experimental accuracy of water content is ±5 ppm mole.

Viscosity and density Measurements and Procedures

All viscosity and density measurements were conducted using an in-house designed and constructed set-up; a schematic view is shown in Figure 4. This setup has been designed to have a maximum working pressure of 200 MPa (density measurements are limited to 140 MPa) and a maximum working temperature of 473.15 K. For these measurements, the set-up was located inside the chamber of an oven, manufactured by BINDER GmbH, capable of being used at temperatures from 200 to 443.15 K. The set-up is comprised of two small cylinders, with volumes of 15 cm 3 , connected to each other through a capillary tube with a measured length of 14.781 metres and a calculated internal diameter of 0.29478 mm. An oscillating U tube densitometer Anton Paar DMA-HPM is connected to the set-up.

.5 Frost Point Measurements and Procedures

The equipment and procedures used for determining the frost points of the solid CO 2 has been originally described by [START_REF] Longman | Measurement and Modeling of CO2 Frost Points in the CO2-Methane Systems[END_REF]. A stainless steel equilibrium cell, which is approximately 11 cm 3 in volume, is submerged in an ethanol bath. The temperature of the ethanol is controlled by a thermostat (LAUDA Proline RP 1290) with a working temperature range 183 to 320 K and an accuracy of ±0.01 K. The ethanol is stirred continuously to maintain a homogeneous temperature distribution in the bath. The equilibrium cell temperature is measured by a platinum resistance thermometer located inside the equilibrium cell. The accuracy of the measured temperature is ±0.05 K. The temperature probe was calibrated against a Prema 3040 precision thermometer and checked by measuring the triple point of pure CO 2 . The equilibrium cell pressure is measured by a Quartzdyne pressure transducer.

Materials

Methane and hydrogen sulphide were purchased from Air Liquide with 99.995 vol% certified purity for methane and 99.5 vol% for hydrogen sulphide. Deionised water was used in all experiments. Carbon dioxide (CO 2 ) has been purchased from BOC with a certified purity higher than 99.995 vol%. Compositions of the synthetic mixtures are given in Table 3.

Thermodynamic and Viscosity Modelling

Fluid and Hydrate Phase Equilibria Model

A detailed description of the original thermodynamic framework used in this work can be found elsewhere [START_REF] Haghighi | Experimental and thermodynamic modelling of systems containing water and ethylene glycol: Application to flow assurance and gas processing[END_REF]Chapoy et al., 2014). In summary, the thermodynamic model is based on the uniformity of fugacity of each component throughout all the phases. The CPA-EoS or the SRK-EoS (if no water is present) is used to determine the component fugacities in fluid phases. The hydrate phase is modelled using the solid solution theory of van der Waals and Platteeuw (1954) as developed by [START_REF] Parrish | Dissociation pressures of gas hydrates formed by gas mixtures[END_REF]. The CPA-EoS binary interaction parameters between components were determined using the group contribution method developed by Jaubert and co-workers. The model to calculate frost points was described by [START_REF] Longman | Measurement and Modeling of CO2 Frost Points in the CO2-Methane Systems[END_REF]. The developed model can predict accurately the distribution of water in the CO 2 or H 2 S-rich phase and solubility of CO 2 or H 2 S in the aqueous phase below and above the critical point of pure CO 2 as shown in Figures 5 and6 [START_REF] Duan | Accurate Thermodynamic Model for the Calculation of H2S Solubility in Pure Water and Brines[END_REF] In this work, the molar volume for CO 2 or a CO 2 -rich mixture, as calculated by the SRK-EoS, is corrected using the exact volume of pure CO 2 at the given T and P.

C EoS new V V V  
(1) Where EoS V is the molar volume obtained from the equation of state. The correction of molar volume in the Eq. ( 2), c V , is defined as:

C i N k i C V x V    1 (2)
x i is the composition of component i in the phase in which the volume is calculated. For CO 2 ,

c i V is defined by MBWR EoS CO Pure C CO V V V   2 2
(3) For the other components, c i V was set to 0. The carbon dioxide density is computed from the MBWR equation in the form suggested by [START_REF] Ely | Thermophysical properties for special high CO2 content mixtures[END_REF] ,will be the same.

) , (

r r r P T f   (5)
Based on the dilute gases considerations and kinetic theory, viscosity at critical point can be approximated as:

6 / 1 2 / 1 3 / 2 c c c T M P   (6)
Where, M denotes the Molecular weight. Thus, the reduced viscosity can be expressed as:

2 / 1 3 / 2 6 / 1 ) , ( ) , ( M P T P T P T c c C r       (7)
For one component as a reference component if the function f in Eq. ( 5) is known, it is possible to calculate the viscosity of any other components, such as component i, at any pressure and temperature. Thus,

                                 i c co i c co o c i c o i o c i c i P P P T T T T T M M P P , , 0 6 1 , , 2 1 3 2 , , ,   (8)
Where, 0 refers to the reference component. Methane with the viscosity data published by [START_REF] Hanley | Equation for the viscosity and thermal conductivity coefficients of methane[END_REF] was selected as the reference fluid in the original Pedersen model. In this work, CO 2 with the viscosity data published by [START_REF] Fenghour | The Viscosity of Carbon Dioxide[END_REF] has been selected as the reference fluid as CO 2 is the major constituent of the stream. The viscosity of CO 2 as a function of density at given T and P can be calculated from the following equation:

) , ( ) ( ) , ( 0 T T T         (9)
Where, η 0 (T) is the zero-density viscosity which can be obtained from the following equation:

  ) ( 00697 . 1 2 / 1 0    T T T    (10)
In this equation, the zero-density viscosity is in units of Pa.s and temperature, T, in K. The reduced effective cross section,

) (   T  
, is represented by the empirical equation:

        4 0 ) (ln ln i i i T a T   (11)
Where the reduced temperature, T * , is given by:

T * =kT/ε (12)
And the energy scaling parameter, ε/k =251.196 K. The coefficients a i are listed in Table 1. The second contribution in Eq. ( 9) is the excess viscosity,

             (13)
Where, the temperature is in Kelvin, the density in kg/m 3 and the excess viscosity in Pa.s.

The coefficients d ij are shown in Table 2. The corresponding states principle expressed in Eq. ( 8) for the viscosity of pure components works well for mixtures. Pedersen et al. (1984) have used the following expression to calculate the viscosity of mixtures at any pressure and temperature.

                             (14)
Where

The critical temperature and pressure for mixtures, according to recommended mixing rules by [START_REF] Murad | Corresponding states correlation for thermal conductivity of dense fluids[END_REF], can be found from:

                                                     N i N j j c j c i c i c j i N i j c i c N j j c j c i c i c j i mix c P T P T z z T T P T P T z z T 3 3 / 1 , , 3 / 1 , , , , 3 3 / 1 , , 3 / 1 , , , (16)
2 3 3 / 1 , , 3 / 1 , , , , 3 3 / 1 , , 3 / 1 , , , 8                                                                N i N j j c j c i c i c j i N i j c i c N j j c j c i c i c j i mix c P T P T z z T T P T P T z z P (17)
The mixture molecular weight is found from

  n n w mix M M M M      303 . 2

 

 N i i i N i i i w M z M z M 2 (19)   N i i i n M z M (20)
The parameter  for mixtures in Eq. ( 14) can be found from: Also,  for the reference fluid can be found from Eq. ( 21) by replacing the molecular weight of the mixture with that of the reference fluid, carbon dioxide. The reduced density, ρ r , is defined as: The procedure below should be followed to calculate the viscosity of CO 2 systems with impurities by the proposed corresponding state principle model:

1. Calculate the T cmix , P cmix and M mix from Eq. ( 16), ( 17) and ( 18), respectively. 3. The mixture parameter,  mix , and  0 should be calculated from Eq. ( 21). 4. The reference pressure and temperature, P 0 and T 0 , should be calculated from Equation ( 15). 5. Calculate the CO 2 reference fluid, ) , ( 0 0 0 T P



in Eq. ( 14) from Eq. ( 9). 6. Calculate the mixture viscosity from Eq. ( 14).

Results and Discussions

All results were compared where possible with experimental values for pure methane, pure carbon dioxide, a synthetic CO 2 -rich fluid (CO 2 : 89.83 mole%; O 2 : 5.05 mole%; Ar: 2.05 mole%; N 2 : 3.07 mole% from [START_REF] Chapoy | Effect of impurities on thermophysical properties and phase behaviour of a CO 2 -rich system in CCS[END_REF] and a typical North Sea natural gas.

For a model to predict accurately the hydrate phase behaviour or transport properties, it is essential that the phase behaviour is correctly predicted, i.e. the phase region, bubble and dew lines. For example for hydrate calculation, the hydrate stability has very sharp temperature dependency above the bubble point, an error in estimating the saturation pressure will lead to high deviations in the hydrate phase behaviour. Viscosity models are also dependant on good density predictions, if a vapour-liquid behaviour is predicted instead of a saturated liquid it will also lead to very high deviations in viscosities. As shown in Figure 7, the SRK-EoS model combined with the group contribution for k ij can predict the phase envelope of the multicomponent systems with good accuracy. The predictions are of greater accuracy for the system containing less carbon dioxide. 1964).

Methane, carbon dioxide and hydrogen sulphide are well known structure I hydrate formers. Hydrate phase equilibria of these systems have been extensively investigated and can be predicted with very high accuracy as seen in Figure 8. Multicomponent systems containing hydrogen sulphide are far scarcer. [START_REF] Hajiw | Effect of acide gases on the methane hydrate stability zone[END_REF] measured the hydrate dissociation conditions for a mixture of methane and hydrogen sulphide. Composition of the fluid is given in Table 1. As the solubility difference between methane and hydrogen sulphide is of several order of magnitude, the hydrate stability zone of this mixture is highly dependent on the fluid to water ratio as seen in Figure 9. The model has also been evaluated with the methane + hydrogen sulphide + carbon dioxide hydrate data reported by [START_REF] Sun | Hydrate Formation Conditions of Sour Natural Gases[END_REF]. Like methane, carbon dioxide and hydrogen sulphide, and all their mixtures are predicted to form structure I hydrate. For these systems, the thermodynamic model is in excellent agreement with these experimental data (within 0.5 K). The ratio between the water mole fraction and the mixture fraction can have a large effect if the concentration of H 2 S and CO 2 is in excess of 10 mole%.

The experimental hydrate dissociation conditions for the synthetic mixtures 2 and 3 in equilibrium with water are plotted in Figure 11. Pure CO 2 , CH 4 and synthetic mixture 2 form structure I hydrate whereas synthetic mixture 3 is predicted to form structure II because of the presence of larger hydrocarbon molecules (propane, i-butane and n-butane). It is also interesting to note that this system, depending on the water to gas ratio, should just enter the phase envelope of the system but displays at higher pressure a liquid like hydrate locus. The system is, over the full pressure range, more stable than pure CO 2 or synthetic mixture, which form structure I. At low and intermediate pressure (P<14 MPa), the system is also more stable than pure CH 4 hydrate, however at higher pressures where hydrates are in equilibrium with a denser supercritical fluid, pure CH 4 hydrates are more stable. [START_REF] Marshall | Hydrates at High Pressures: Part I. Methane-Water, Argon-Water, and Nitrogen-Water Systems[END_REF]). (), Natural gas hydrate stability zone (Chapoy and Tohidi, 2011).

Hydrate dissociation measurements in absence of a free water phase (dry systems) are difficult to measure by conventional techniques; the alternative is to measure the water content in the fluid phase in equilibrium with hydrates. Unfortunately data for acid gas systems are limited and the few data references available are not in agreement. Experimental data on water contents for CO 2 in equilibrium with hydrates from 223.15 to 263.15 K up to 10 MPa have been measured by Burgass et al. (2014). [START_REF] Chapoy | Effect of impurities on thermophysical properties and phase behaviour of a CO 2 -rich system in CCS[END_REF] have reported water content at 15 MPa from 233.15 to 288.15 K. Song and Kobayashi have reported measurements of water content in presence of CO 2 However, the reliability of these studies has been recently questioned as seen in Figure 12, some of the data reported by [START_REF] Song | The water content of CO2-rich fluids in equilibrium with liquid water and/or hydrates[END_REF] show large deviations with the developed models and the available literature data.

For multicomponent acid gas mixtures at low temperatures, only data from our laboratory are available. Experimental water content data in equilibrium with hydrates for pure CO 2 and 2 multi-component systems (synthetic mixtures 2 & 3) are plotted along with predictions of the thermodynamic model in Figure 12. As can be seen from the figure the experimental and predicted data are in good agreement with some deviation (AAD≈ 5%). As expected, less water can be dissolved in the multicomponent systems than in pure CO 2 , because the amount of water that can be dissolved in hydrocarbons is lower than liquid CO 2 at the same temperature and pressure. As seen in this figure, the water contents for the synthetic mixtures are between the water contents of pure CO 2 and pure CH 4 , with the water content in mixture 3 (lower CO 2 content) closer to pure CH 4 . The experimental and modelling results for the viscosity of the synthetic mixtures are plotted in Figures 13 and14. All experiments for the streams were conducted at pressures above saturation or in the supercritical region and then at low pressures, i.e., in the single gas phase region. The viscosity of each conducted test was calculated using the modified Pedersen model. As can be seen from the figure, the model predictions and experimental data are in good agreement. Both in the liquid and the supercritical regions, the viscosity is increasing with pressure and decreasing with temperature. In the vapour region, pressure has a weak effect on viscosity and the viscosity is increasing with temperature as expected for low pressure gaseous systems. Viscosities of a "typical" natural gas (low CO 2 concentration) were also added in Figure 14 for comparison purposes (data from [START_REF] Kashefi | Viscosity of binary and multicomponent hydrocarbon fluids at high pressure and high temperature conditions: Measurements and predictions[END_REF]. The viscosities of both the multi-component systems (synthetic mixtures 2 and 3) fall between the viscosities of pure CO 2 and of the natural gas from [START_REF] Kashefi | Viscosity of binary and multicomponent hydrocarbon fluids at high pressure and high temperature conditions: Measurements and predictions[END_REF].

Densities of the multi-component systems were measured at different pressures and temperatures in gas, liquid and supercritical regions. Both experimental and modelling results with density correction are plotted in Figure 15 for the synthetic mixture 2. By employing the density correction and the SRK-EoS the absolute average deviation is 1.7%. It is interesting to note the peculiar behaviour of density at temperatures above the critical temperature compared to the density of pure CO 2 as shown in Figure 16 in which the difference between the density of the two multicomponent mixtures and pure CO 2 is plotted at about 423.15 K. For the CO 2 mixture at the stated temperature, a maximum reduction in density from that of pure CO 2 occurs at a certain pressure. The maximum reductions for the multicomponent systems 2 and 3 are 180 and 300 kg/m 3 at a pressure of around 12 and 14 MPa, respectively (i.e., 35 % and 60% reduction).

Frost point measurements were carried out for the synthetic mixture 2. The experimental data for the system is plotted together with the predicted solid CO 2 phase boundary of the mixture and the pure CO 2 phase diagram in Figure 17. As seen in the figure, for this system the thermodynamic model using the group contribution k ij tuned on VLE data is in good agreement with the new experimental data.

Figure 16.

Predicted and experimental density difference ρ= ρ MIX -ρ CO2 , between synthetic mixture 2 (),synthetic mixture 3 ()and pure CO 2 density at 323.15 K. Lines: Predictions using the corrected SRK-EoS model (Grey line is for methane).

Conclusions

Knowledge on the phase behaviour and thermophysical properties of CO 2 -rich and acid gas systems is of great importance for CCS and developing sour gas reservoirs, as well as testing predictive models. However, limited published data sets are available for such systems. In this communication the phase behaviour and some properties of different acid gas streams have been studied, such as the phase envelope, the hydrate stability, dehydration requirement, viscosity and density of the mixture. Models have been developed to calculate and predict these properties.

Future work will concentrate on the determination/measurement and modelling of properties for other types of natural gases (different CO 2 concentrations, impact of H 2 S, etc.).

Figure 1 .

 1 Figure 1. Plot showing an example of bubble point determination from plot of change in cell pressure versus volume

Figure 2 .

 2 Figure 2. Plot showing an example of dew point determinations from equilibrium stepheating data using the isochoric method.

Figure 3 .

 3 Figure3. Dissociation point determination from equilibrium step-heating data. The equilibrium dissociation point is determined as being the intersection between the hydrate dissociation (pressure increase as a result of gas release due to temperature increase and hydrate dissociation, as well as thermal expansion) and the linear thermal expansion (no hydrate) curves.

Figure 4 .

 4 Figure 4. Schematic drawing of the viscosity -density set-up

  .

Figure 5 .

 5 Figure 5. Pxy, Phase equilibria in the carbon dioxide + water system at 298.15 K left) and 423.15 K (right). Black Lines: Model predictions. Left figure: () Experimental data from Wiebe and Gaddy (1941); () Experimental data from Gillepsie and Wilson (1982) ; () Experimental data from Nakayama et al. (1987) ; () Experimental data from King et al. (1992) ; (): Experimental data from Hou et al. (2013) ; (): Experimental data from Valtz et al. (2004). Right figure: (): Experimental data from Takenouchi and Kennedy (1964); () Experimental data from Gillepsie and Wilson (1982); (): Experimental data from Tabasinejad et al. (2011) at 422.98 K ;(): Experimental data from Hou et al. (2013).(): calculated data[START_REF] Duan | An improved model calculating CO 2 solubility in pure waterand aqueous NaCl solutions from 273 to 533 K and from 0 to 2000 bar[END_REF]

  how the viscosity can change as a function of density outside of the critical region. The excess viscosity term can be correlated as follows:

M

  are the weight average and number average molecular weights, respectively.

  of carbon dioxide, ρ 0 , is equal to 467.69 kg/m 3 . The Modified Benedict-Webb-Rubin (MBWR) equation of state has been applied for computing the reference fluid density, ρ 0 , at the desired pressure and temperature of

Figure 7 .

 7 Figure7. Experimental and predicted phase envelope of the CO 2 -rich mixture. (),Synthetic Mix 3. (), Synthetic Mix 2. Black lines: bubble lines predictions using the SRK-EoS; Dotted black lines: dew lines predictions using the SRK-EoS; Grey broken lines: predicted vapour pressure of pure CO 2 using the SRK-EoS.

Figure 12 .

 12 Figure 12. Water content of CO 2 -rich fluids. Figure on the left: Water content of CO 2 in equilibrium with hydrates at 13.79 MPa and different temperatures. (): Chapoy et al. (2011); (): data from Song and Kobayashi (1984). Figure on the right: Experimental water content data and predictions for the water content of pure CH 4 , pure CO 2 and synthetic mixures 2 &3 at 15 MPa and different temperatures.

Figure 13 .Figure 14 .Figure 15 .

 131415 Figure 13. Predicted and experimental viscosity of synthetic mixture 2. Black lines: Predictions using CSP model. Black Dotted lines: Predictions using CSP model at the bubble and dew pressures of the system. Data inside the grey box are plotted in Fig. 9. This work: (), T = 243.15 K (), T = 253.15 K (), T = 273.15 K (), T = 283.15 K (), T = 298.15 K (), T = 323.15 K (), T = 373.15 K (), T = 423.15 K Figure 14.Predicted and experimental viscosity of synthetic mixtures 2 and 3 at 323.15 K. Black and dotted lines: Predictions using the modified CSP model. Grey lines: pure CO 2 viscosity. Grey broken lines: Predictions using the original CSP model. (), synthetic mixture 2; (),synthetic mixture 3; () data from[START_REF] Kashefi | Viscosity of binary and multicomponent hydrocarbon fluids at high pressure and high temperature conditions: Measurements and predictions[END_REF] for a Natural gas (in mole% C1: 88.83; C2: 5.18; C3:1.64; iC4: 0.16; nC4: 0.27; iC5: 0.04; CO2: 2.24; N2: 1.6)

Figure 17 .

 17 Figure 17. Experimental and predicted phase envelope of the synthetic mixture 2. (), bubble/dew points. (), Frost points. Black lines: bubble lines predictions using the SRK-EoS; Dotted black lines: dew lines predictions using the SRK-EoS; Grey broken lines: phase diagram of pure CO 2 .

Table 1 .

 1 Composition, mole% each component, of the multicomponent mixtures Components

		Synthetic Mix 1	Synthetic Mix 2	Synthetic Mix 3
		(from Hajiw et al., 2014)	(from Chapoy et al., 2013)
	CO 2	-	Balance (89.83)	Balance (69.99)
	H 2 S	20.0	-	-
	Methane	80.0	-	20.02 (±0.11%)
	Ethane	-	-	6.612 (±0.034%)
	Propane	-	-	2.58 (±0.013%)
	i-Butane	-	-	0.3998 (±0.004%)
	n-Butane	-	-	0.3997(±0.004%)
	Nitrogen	-
	Total	100	100	100

  Pxy, Phase equilibria in the hydrogen sulphide + water system. Black Lines: Model predictions. Left figure: water content in the H 2 S rich phases. Right figure: H 2 S solubility in the aqueous phase. (): calculated data

			310.9 K Selleck et al. (1952)		0.1
			310.9 K Gillepsie and Wilson (1982)
							303.15 K	333.15 K
			344.3 K Selleck et al. (1952)		0.09
			344.3 K Gillepsie and Wilson (1982)
	MPa		372 K Gillepsie and Wilson (1982) 373.5 K Gillepsie and Wilson (1982) 422 K Gillepsie and Wilson (1982) 423.15 K Burgess and Germann (1969) 444.3 K Selleck et al. (1952)	fraction	0.07 0.08	363.15 K 423.15 K	393.15 K 453.15 K
	Pressure/		477.6 K Gillepsie and Wilson (1982)	x H2S / mole	0.03 0.04 0.05 0.06	513.15 K
						0.02
						0.01
						0
	1	0.2	0.3	0.4	0.5	0	5	10	15	20
								Pressure/ MPa
	y	H2O / mole fraction
	Figure 6.

Table 1 .

 1 Values of Coefficients a

	Eq. (11)	i for CO 2 in

Table 2 .

 2 Values of coefficients d ij in Eq. (13)

	i 0 1 2 3 4	a i 0.235156 -0.491266 5.211155 x 10 -2 5.347906 x 10 -2 -1.537102 x 10 -2	d ij d 11 d 21 d 64 d 81 d 82	Value 0.4071119 x 10 -2 0.7198037 x 10 -4 0.2411697 x 10 -16 0.2971072 x 10 -22 -0.1627888 x 10 -22

  Experimental and predicted hydrate stability of CO 2 , H 2 S and methane in equilibrium with liquid water. (), pure CO 2 hydrate stability zone(Chapoy et al. 2011);

		40
		35
		30
	Pressure /MPa	15 20 25
		10
		5
		0
		273.15	283.15	293.15	303.15	313.15
				Temperature / K
	Figure 8.

()

, pure H 2 S hydrate stability zone

[START_REF] Selleck | Phase behavior in the hydrogen sulfidewater system[END_REF]

); (): pure CH 4 hydrate stability zone

[START_REF] Nixdorf | Experimental determination of hydrate equilibrium conditions for pure gases, binary and ternary mixtures and natural gases[END_REF]

;(), pure CH 4 hydrate stability zone

(Marshall et al.

  Experimental and predicted hydrate stability of 80 mole% CH 4 + 20 mole% H 2 S system in equilibrium with liquid water. (), pure H 2 S hydrate stability zone[START_REF] Selleck | Phase behavior in the hydrogen sulfidewater system[END_REF]; (): 80 mole % + 20 mole% H 2 S;(), pure CH 4 hydrate stability zone[START_REF] Marshall | Hydrates at High Pressures: Part I. Methane-Water, Argon-Water, and Nitrogen-Water Systems[END_REF]. Predictions for the mixtures were performed for the reported aqueous mole fraction. Figure 11. Experimental and predicted hydrate of some selected acid gas systems in equilibrium with liquid water. (), pure CO 2 hydrate stability zone (Chapoy et al. 2011); (): pure CH 4 hydrate stability zone (Nixdorf and Oellrich, 1997) ;(), pure CH 4 hydrate stability zone

		40	12
		35	10
	Pressure /MPa	30 15 20 25	Pressure /MPa	2 4 6 8
		10		0
		5		273.15		283.15	293.15 Temperature / K	303.15	313.15
		0
		273.15	278.15	283.15	288.15	293.15	298.15
						Temperature / K
			10
				9
				8
				7
			Pressure /MPa	4 5 6				H2S
				3
				2
				1
				0
				273.15	278.15	283.15	288.15	293.15	298.15	303.15
							Temperature / K
	Figure 10. Experimental and predicted hydrate stability of the ternary CH 4 + H 2 S +CO 2
	system in equilibrium with liquid water. (), pure H 2 S hydrate stability zone (Selleck et al.
	1952); (), pure CH 4 hydrate stability zone (Marshall et al. 1964); () 75.48 mole% CH 4 +
	6.81 mole% CO 2 + 17.71 mole% H 2 S (Sun et al., 2003); ():87.65 mole% CH 4 + 7.4 mole%
	CO 2 + 4.95 mole% H 2 S (Sun et al., 2003); ():82.45 mole% CH 4 + 10.77 mole% CO 2 +
	6.78 mole% H 2 S (Sun et al., 2003) ); ():82.91 mole% CH 4 + 7.16 mole% CO 2 + 9.93
	mole% H 2 S (Sun et al., 2003) ; ():75.48 mole% CH 4 + 6.81 mole% CO 2 + 17.71 mole%
	H 2 S (Sun et al., 2003).

H2S

Figure

9

.

3.07(±0.04%) -Oxygen -5.05(±0.01%) -Argon -2.05(±0.06%) -

Acknowledgements

This research work was part of a Joint Industrial Project (JIP) conducted jointly at the Institute of Petroleum Engineering, Heriot-Watt University and the CTP laboratory of MINES ParisTech. The JIPs supported by Chevron, GALP Energia, Linde AG, OMV, Petroleum Experts, Statoil, TOTAL and National Grid Carbon Ltd, which is gratefully acknowledged. The participation of National Grid Carbon in the JIP was funded by the European Commission's European Energy Programme for Recovery. The authors would also like to thank the members of the steering committee for their fruitful comments and discussions.


N Yakdi

F Huet

K Ngo

In-situ particle sizing at millimeter scale from electrochemical noise: simulation and experiments

Keywords: Particle sizing techniques, Two-phase systems, Electrochemical noise, Electrolyte resistance

Over the last few years, particle sizing techniques in multiphase flows based on optical technologies emerged as standard tools but the main disadvantage of these techniques is their dependence on the visibility of the measurement volume and on the focal distance. Thus, it is important to promote alternative techniques for particle sizing, and, moreover, able to work in hostile environment. This paper presents a single-particle sizing technique at a millimeter scale based on the measurement of the variation of the electrolyte resistance (ER) due to the passage of an insulating sphere between two electrodes immerged in a conductive solution. A theoretical model was proposed to determine the influence of the electrode size, the interelectrode distance, the size and the position of the sphere, on the electrolyte resistance. Experimental variations of ER due to the passage of spheres and measured by using a home-made electronic device are also presented in this paper. The excellent agreement obtained between the theoretical and experimental results allows validation of both model and experimental measurements. In addition, the technique was shown to be able to perform accurate measurements of the velocity of a ball falling in a liquid.

Introduction

Many industrial processes highlight the importance of multiphase flows, in particular those involving two phases where a gaseous, liquid or solid phase is dispersed at various scales in a continuous liquid phase. Examples of such processes are encountered in a wide variety of technical applications, including inkjet printers, deposition and coating, direct injection internal-combustionengines, gas turbines, etc. The dispersed phase can lead to homogeneous or separated flows, suspensions or bubbly flows, etc. Among the two-phase systems, those with dispersed discrete particles play an important role, for example in the pharmaceutical and cosmetic industries. Hence, it is essential to seek reliable techniques to characterize the particles in mixtures because they have a direct impact on product quality [START_REF] Kleinstreuer | Two-Phase Flow: Theory and Applications[END_REF].

The various methods for particle sizing in multiphase flows have been reviewed in [START_REF] Tayali | Particle sizing techniques in multiphase flows: A review[END_REF][START_REF] Bieberle | Experimental two-phase flow measurement using ultra fast limited-angle-type electron beam X-ray computed tomography[END_REF][START_REF] Xu | Comparison of sizing small particles using different technologies[END_REF][START_REF] Snowsill | Particle Sizing[END_REF].

Recently, the techniques based on optical technologies have emerged as standard tools because of their numerous advantages such as a high spatial resolution, a moderate temporal resolution, and a large range of applicability. They are widely used for dense multiphase flows in industrial processes and in microstructures as well, for example for biological cell counting and sizing [START_REF] Gouesbet | Laser-based optical measurement techniques of discrete particles: A review[END_REF][START_REF] Honkanen | Digital imaging measurement of dense multiphase flows in industrial processes[END_REF][START_REF] Silva | Particle sizing measurements in pharmaceutical applications: comparison of in-process methods versus offline methods[END_REF][START_REF] Telleman | Cell counting and cell sizing in microstructures, in: Lab-on-a-Chip, Miniaturized Systems for (Bio)Chemical Analysis and Synthesis[END_REF][START_REF] Satinover | Single-particle optical sizing of microbubbles[END_REF][START_REF] Ma | New developments in particle characterization by laser diffraction: size and shape[END_REF]. In addition to complex and expensive equipments, the main disadvantage of these techniques is their dependence on the visibility of the measurement volume and on the focal distance. Moreover, they are difficult to adapt to hostile environments, for example at high temperature and/or high pressure.

Since particles of different size can be physically separated, non-optical techniques for particle sizing have also been developed like sieving, gravitational sedimentation, centrifuging, etc [START_REF] Allen | Particle size analysis by sieving[END_REF][START_REF] Alexander | Gravitational sedimentation of gold nanoparticles[END_REF][START_REF] Neumann | New method for density determination of nanoparticles using a CPS disc centrifuge[END_REF]. Sensors for particle sizing have been developed using electrostatic and acoustic techniques [START_REF] Bastari | Acoustic emissions for particle sizing of powders through signal processing techniques[END_REF][START_REF] Alba | Acoustic spectroscopy as a technique for the particle sizing of high concentration colloids, emulsions and suspensions[END_REF][START_REF] Tajdari | New technique to measure particle size using electrostatic sensor[END_REF]. The most popular technique based on electrical measurement is the electrical sensing zone method (Electrozone Counter or Coulter Principle), in which particles are suspended in an electrolyte and forced to flow through a small orifice. The change in electrical resistivity of the solution between two electrodes placed on either side of the orifice allows the particle size to be determined [START_REF] Hurley | Sizing particles with a Coulter counter[END_REF][START_REF] Rodriguez-Trujillo | Label-free protein detection using a microfluidic Coulter-counter device[END_REF][START_REF] Nieuwenhuis | Integrated Coulter counter based on 2-dimensional liquid aperture control[END_REF]. This method has been applied to many applications but its main disadvantage is that particles may get trapped in the orifice and cause plugging.

During the past thirty years, the electrochemical noise (EN) technique has been developed in our laboratory to study two-phase flows using electrodes immerged in a conductive electrolyte. Indeed, the presence of discrete entities (gas bubbles, oil droplets, solid particles…) near a working electrode generates fluctuations of current, potential, and also fluctuations of the electrolyte resistance (ER) between the reference electrode and the working electrode that are provoked by changes in the current distribution (screening effects) [START_REF] Gabrielli | Characterization of electrolytic bubble evolution by spectral analysis. Application to a corroding electrode[END_REF][START_REF] Gabrielli | Investigation of water electrolysis by spectral analysis, I. Influence of the current density[END_REF][START_REF] Gabrielli | Dynamic behaviour of an electrolyser with a two phase solid-liquid electrolyte. Part II: Investigation of elementary phenomena and electrode modelling[END_REF][START_REF] Bouet | Electrochemical resistance noise during composite plating[END_REF][START_REF] Bouazaze | Measurement of electrolyte resistance fluctuations generated by oil-brine mixtures in a flow-loop cell[END_REF]. The analysis of the EN signals in both time and frequency domains allows a more precise identification of the elementary events at the origin of noise, which is not possible with classical electrochemical methods, such as impedance spectroscopy or cyclic voltammetry, that give information averaged in time and over the electrode surface. This technique has been successfully applied to determine the size and departure rate of bubbles on a gas-evolving electrode [START_REF] Gabrielli | Characterization of electrolytic bubble evolution by spectral analysis. Application to a corroding electrode[END_REF][START_REF] Gabrielli | Investigation of water electrolysis by spectral analysis, I. Influence of the current density[END_REF], to determine the electrical charge exchanges between conductive beads and a current collector [START_REF] Gabrielli | Dynamic behaviour of an electrolyser with a two phase solid-liquid electrolyte. Part II: Investigation of elementary phenomena and electrode modelling[END_REF], to get a detailed view of the approach and residence of microcapsules close to the electrode during composite plating [START_REF] Bouet | Electrochemical resistance noise during composite plating[END_REF], or to characterize the composition of oil-brine mixtures in a flow-loop cell [START_REF] Bouazaze | Measurement of electrolyte resistance fluctuations generated by oil-brine mixtures in a flow-loop cell[END_REF].

The objective of this work was to assess the possibility to use the electrochemical noise as an alternative technique for particle sizing in two-phase flows. As a preliminary step, single particles were used in this first paper. Numerical simulations were performed to study the influence on the ER of an insulating sphere passing through two electrodes immerged in a conductive electrolyte and an experimental work was carried out to validate the simulation results.

Theoretical simulation

In order to determine the influence of the presence of an insulating sphere on the ER measured between two electrodes immersed in a conductive electrolyte, a theoretical model has been developed using COMSOL Multiphysics 4.1 software. As shown in Fig. 1, two disks of diameter d e represent the two electrodes positioned face-to-face. The electrolyte is contained in a cuboid of dimensions l x , D, l z , where D is the distance between the two electrodes, and d p is the diameter of the sphere that will descend along the z-axis in the experimental section. The l x and l z values are chosen large enough (100 mm) compared to the diameter of the electrodes and the interelectrode distance, so that they have no significant influence on the simulation results.

To calculate the ER (impedance at high frequency) an electrical model was chosen in this work instead of an electrochemical model to simplify the simulation. This model is valid as long as the concentration and the temperature of the electrolyte, and hence its conductivity, , are constant.

A difference of potential is applied between the two electrodes and the ER is calculated from the resulting current on the electrodes. Only the primary current and potential distributions were considered in this work since kinetic and mass transport phenomena can be neglected at the high frequency, typically 100 kHz, at which the ER is measured. The potential, , can be calculated at any point between the two electrodes using the Laplace equation:

²  = 0 (1)
The first integration of this equation allows the current density, J, to be determined:

J = -   (2)
The following boundary conditions were applied in the simulation: i) the normal component of the current, -n.J, was equal to 0 at any point of all the insulating surfaces of the sphere and around the electrodes, as well as on the boundaries of the cuboid.

ii) a difference of potential, V app , of 30 mV was applied between the two electrodes (+30 mV was applied to one electrode and 0 mV to the other one that was considered to be the grounded electrode).

The value of the electrolyte conductivity,  = 28 S/m (resistivity of 3.6  cm), was taken as that of the electrolyte used in the experimental part presented below.

As an example, Fig. 2 presents the potential and current density distributions in the presence of an insulating sphere. To improve the accuracy of the simulation, the calculations were performed using the finite element method with a more precise mesh on the surface of both electrodes. It may be noted in Fig. 2b that the current lines are strongly deflected close to the electrode edges and around the sphere.

The ER measured between the two electrodes was determined from the ratio V app /I, where I is the current flowing between the electrodes. I was calculated by integrating the current density on the surface of the non-grounded electrode. The variation of the ER, R e , due to the presence of the sphere is then:

ea ep e R R R    (3)
where R ep and R ea are the ER in the presence and in the absence of the sphere between the electrodes, respectively. The normalized ER variation, R e,norm , defined as:

100 ) ( ea ea ep norm e,     R R R R ( 4
)
depends on 6 parameters, the position (x, y, z) of the sphere centre, the sphere diameter d p , the electrode diameter d e , and the interelectrode distance D. The influence of these parameters on R e,norm is presented in the following sections in order to assess the possibility to determine the sphere diameter in single-particle sizing experiments, in which d e and D are fixed parameters.

Influence of the position of the insulating sphere

Fig. 3 shows the R e and R e values for an insulating sphere of 6 mm in diameter positioned at several positions along the vertical z axis and various x positions (y = 0) between two electrodes of diameter d e = 5 mm and an interelectrode distance D = 20 mm. Because of the symmetry of the problem, the variations of R e as a function of x or z are identical; for example, R e = 6.3664  for z = 0 and x = 5 mm, and for x = 0 and z = 5 mm. It can be noticed that R e starts increasing much before the sphere passes in front of the electrodes in the z or x direction, at a distance about three times larger than the electrode diameter, indicating that the current lines are modified on long distances by the presence of the sphere. The maximum value of R e is obtained when the centre of the sphere is positioned exactly on the axis of the two electrodes (x = z = 0), where the screening effect is the most important. Consistently, the variation of R e is symmetrical around the origin z = 0.

It must be noted that, despite a sphere diameter slightly larger than the electrode diameter, the maximum value of the ER increase due to the presence of the sphere is rather low: R e = 46.1 m at z = 0, which gives a normalized ER increase R e,norm = 0.73%.

Fig. 4 shows the ER values when the sphere is not at the same distance of both electrodes (y  0). The highest values are obtained when the sphere is moved closer to one of the electrodes (Fig. 4a) because of a stronger constriction of the current lines between the sphere and the closest electrode. For this reason, the maximum value of the ER variations is obtained when the sphere is in contact with one of the electrodes (R e = 1.16  for y = 7 mm and x = 0). It should be noticed that when the centre of the sphere is on the axis of the disk electrodes (x = z = 0), the ER values are relatively close together in the middle of the channel (R e = 53 m 15% for |y| <2 mm), as shown in Fig. 4b. This indicates that in single-particle sizing experiments, it is important to force the particle to flow at approximately the same distance of both electrodes since in that case R e,norm depends on x, z, and d p only.

Influence of the insulating sphere diameter d p

The influence of the diameter of the insulating sphere on the ER increase was investigated for spheres centred at the origin O (x = y = z = 0). The normalized ER variations, which are given in Fig. 5 for three interelectrode distances, increase with the sphere diameter, the upper limit corresponding to the case where the sphere touches both electrodes (R e,norm = 99%, 156%, and 179% for D = d p = 10, 15, and 20 mm, respectively). The analytical dependence of R e,norm on the sphere diameter is not trivial; two behaviours can be observed in Fig. 5, one for d p < D/2, the other one for d p > D/2 where the R e,norm increase with d p is stronger. A power law can be observed in Fig. 5 for small spheres and all three interelectrode distances:

b d a d R p p norm e, ) (   (5)
where a = 0.039, 0.0087, and 0.0031 for D = 10, 15, and 20 mm, respectively, and b = 3.05  0.05 in all 3 cases, according to the linear fitting performed on the (log d p , log R e,norm ) data. In the study of the ER variation induced by a spherical bubble in contact with an infinite electrode [START_REF] Sides | Primary potential and current distribution around a bubble on an electrode[END_REF] or a disk electrode [START_REF] Fransaer | Variations in the flow of current to disk electrodes caused by particles[END_REF][START_REF] Bouazaze | Electrolyte-resistance change due to an insulating sphere in contact with a disk electrode[END_REF], a power close to 3 was already obtained for the exponent of the sphere diameter. The normalized ER increment due to the presence of the sphere is given in Fig. 6b in a loglog scale. In Region 2, as soon as D is larger than 2 d p , a power law can be observed: [START_REF] Gouesbet | Laser-based optical measurement techniques of discrete particles: A review[END_REF] where the exponent d is very close to 4. For the smaller spheres (Region 1), a power law may also be observed, but at values of D much higher than 2 d p and with an exponent slightly lower than 4.

Influence of the interelectrode distance D

d D c D R    ) ( norm e,
All fitted curves shown in Fig. 6b were obtained by a linear fitting performed on the (log D, log R e,norm ) data.

Influence of the electrode diameter d e

The influence of the electrode diameter on the ER is shown in Fig. 7 for a sphere of diameter 6 mm centred at the origin O. As expected, the smaller the electrode, the higher the ER value because of a stronger constriction of the current lines close to the electrodes. It should be noted that the influence of the sphere is so small, according to the low R e,norm values, that the R e values in the absence or in the presence of the sphere would not be distinguishable in Fig. 7. In contrast, R e,norm increases with the electrode diameter, indicating that the passage of particles can be more easily detected when using electrodes larger than the particle.

Method for particle sizing

Once the electrode diameter and interelectrode distance have been chosen, the ER variation R e induced by the passage of a single spherical particle depends on 4 parameters, the position (x, y, z) of the particle centre and the particle diameter d p . According to Fig. 4b, the influence of the lateral position y on R e is not significant when the particle flows close to the middle of the channel between the electrodes, so that R e depends only on x, z and d p under this experimental condition.

For a given sphere diameter, an empirical expression of R e as a function of z can be determined for x = 0 with the help of the plot of R e (z) in a semi-logarithmic scale that indicates a linear zdependence of log(R e ) for high z values. As shown in Fig. 8a

. 0 0668 . 0 ) mm 6 , , 0 (        (z in mm) (7)
Because of the symmetry of the problem, the isovalues of R e are located on circles of equation "x 2 +z 2 = constant" when y = 0, so that the previous empirical expression can be directly extended to any position (x, z) of the sphere centre by replacing

z x z x d z x R         (x, z in mm) (8)
which is confirmed by the excellent agreement in Fig. 8a between the R e values for x  0 and the solid red lines obeying Eq. 8.

Finally, to derive an expression of R e including the dependency on d p , the power law identified in Eq. 5 for d e < D/2 was used despite the fact it was obtained for a sphere centred at the origin O.

               d z x d z x R z x (x, z, d p in mm) (9)
In single-particle sizing experiments, the diameter d p of the sphere must be determined from the experimental transient R e (z) such as those presented below when the sphere flows vertically in the channel. Still assuming that y is close to 0, the amplitude R e,max of the transient at z = 0, which can be easily determined, depends on x and d p . While a fitting procedure using Eq. 9 could give both x and d p values, an easier way consists in considering experiments in which the particle flows at a distance x close to 0. Indeed, the influence of x on R e,max is relatively low in that case (for |x| <2 mm, R e,max = 43.9 m  5% for d p = 6 mm and 5.4 m  5% for d p = 3 mm). As a consequence, the diameter d p can be evaluated with the following expression derived from Eq. 9 with x = z = 0. This expression will be used below in the experimental part. ( R e,max in , d p in mm) (10)

Experimental

All experiments were performed at room temperature in a large glass crystallizing dish of 17 cm in diameter and 10 cm in height. The two working electrodes were the cross-section of platinum rods of 5 mm in diameter embedded in epoxy resin and mounted flush in two parallel planes of a

Teflon holder located at the dish centre. The distance between the two electrodes positioned face to face was 20 mm. The electrolyte was a 0.7 M H 2 SO 4 solution prepared with deionised water. Its conductivity, which depends on the room temperature, was measured before each experiment with a conductivity meter (Meterlab CDM230) and adjusted to the value of 28 S/m by adding a few droplets of sulfuric acid or water. Nevertheless, the room temperature could change a little bit when performing long experiments, which had an effect on the amplitude R e of the experimental transients and on the mean value of the ER measured in the absence of sphere. However, this effect was not significant on the derived value of R e,norm .

Alumina spheres of various diameters, 3, 5, 6, 8, and 10 mm of excellent sphericity (tolerance 2.5 m in diameter) from Goodfellow were used as insulating spheres. To ensure that the insulating sphere could fall the most exactly possible at a given position between the two electrodes, a specific Teflon funnel with a narrow stem of diameter slightly larger than that of the sphere was used for each sphere to guide it. Moreover, the lower part of the stem was immerged in the solution to avoid any disturbance of the sphere trajectory when contacting the electrolyte.

The ER fluctuations were measured using the home-made electronic device shown in Fig. 9 [24]. The two electrodes were polarized at the same potential using a home-made potentiostat applying a 0 mV potential difference between the two electrodes, one acting as the working electrode, the other as counter electrode. No reference electrode was used. A sinusoidal voltage signal, v hf , of peak-to-peak amplitude about 50 mV and high frequency f hf = 100 kHz was applied between the two electrodes to measure the ER. Indeed, at this high frequency, the effect of the double layer is eliminated and only the ER has to be considered in the electrochemical impedance.

The sine wave current response, v hf /R e , measured across the resistor R, was passed through a capacitor to remove any residual DC offset, amplified and high-pass filtered at 20 kHz to eliminate low-frequency current fluctuations. The signal was then rectified with a diode and low-pass filtered at 10 kHz, so that the amplitude of the voltage output signal, v R e , obeyed the following relationship:

d R b v R   e 1 e (11)
The parameters b and d were determined through a calibration procedure in which the electrochemical cell was replaced by pure resistors. From Eq. 11, the mean value of the ER, R e , could be calculated from the mean value of v R e measured during the electrochemical experiments.

After DC-offset elimination, amplification (gain G 3 = -10), and analog low-pass filtering to prevent aliasing, the variations of v R e due to the passage of spheres were acquired with a real-time data acquisition system (Concurrent Computer, Maxion 9100) at different sampling frequencies. The variations of the ER, R e , could be derived from the measured variations v R e with the equation:

e 3 2 e e R v G b R R     (12)
obtained by differentiating Eq. 11 for R e variations of small amplitude. Experimental results and comparison with theory are presented in the next section for insulating spheres of various diameters falling between the two electrodes at a position x and y close to 0. with Eq. 9 by replacing z with v(t 0 -t), where v is the velocity of the sphere and t 0 the time instant at which the sphere passes at the vertical position z = 0. The agreement is so good between the experimental and theoretical curves that it is hardly possible to distinguish them for the spheres of 10 mm and 6 mm in diameter. The adjusted parameters x, v, t 0 , and d p are given in Table 1, together with characteristic parameters of the hydrodynamic flow.

Results and discussion

As expected from the experimental procedure, the x values are close to 0. On the other hand, considering that the starting point of the sphere fall was far enough of the electrodes, the sphere velocity could be assumed to be constant in front of the electrodes. This was confirmed by the perfect symmetry of the v R e transients. Moreover, according to Table 1, the sphere velocity is approximately proportional to the square root of the sphere radius. This results from the fact that the drag force of the falling sphere is quadratic with velocity at high Reynolds numbers (Re > 10 3 ) [START_REF] Owen | The effects of linear and quadratic drag on falling spheres: an undergraduate laboratory[END_REF].

Table 1 actually shows high values of the Reynolds number calculated according to the equation: [START_REF] Alexander | Gravitational sedimentation of gold nanoparticles[END_REF] where  e and  are the density and dynamic viscosity of the electrolyte, respectively ( e = 1.04 g cm -3 and  = 1.0 g m -1 s -1 ). Table 1 also shows the terminal velocity of the falling ball, v ter , calculated as follows:

   v d Re p e
) 1 ( 3 4 e d p ter     C d g v ( 14
)
where g is the acceleration of free fall (g = 9.81 m s -2 ),  is the density of the falling sphere ( = 3.95 g cm -3 ), and C d is the drag coefficient, which depends on the Reynolds number according to the relation:

4 . 0 1 6 24 d     Re Re C ( 15
)
The terminal velocities of these not-so-small spheres (compared to the interelectrode distance) have been calculated as if they were in an unbounded fluid. The wall-proximity effects in the near-wall region were not handled, so the drag coefficient was probably underestimated, which could explain why the terminal velocities had values larger than those of the measured velocities (by a factor of about 25%).

The estimated values of the sphere diameter, d p , in Table 1 are in good agreement with the actual diameter of the spheres, within an error lower than 5%, considering the fact that the position of the sphere in the x and y direction was not perfectly controlled during its fall. A non-zero x value leads to a lower R e value (Fig. 3) while, in contrast, a non-zero y value leads to a higher value.

According to the simulations, a y value of about 1 mm explains the overestimated diameter of 10.488 mm for the 10 mm-diameter sphere. The straightforward Eq. 10 gives the same values for the diameter of the spheres.

Fig. 11 shows that it would be difficult to characterize spheres of diameter much smaller than 3 mm with the present experimental set-up (D = 20 mm, d e = 5 mm) without improving the measurement of the ER variations by decreasing its background noise. Another solution consists in using closer electrodes as, for example, in microreactors. This will be presented in a future paper.

Additional experiments were performed with spheres of various sizes. In order to insure the reproducibility of the results, at least six experiments were carried out for each sphere diameter.

From the v R e transients measured due to the passage of the insulating spheres on the vertical line x = y = 0, the maximum amplitude, R e,max , of the R e transients was estimated. The standard deviation of the normalized ER increment, R e,max , norm was then calculated with the formula:

      n R R n 2 norm max, e, norm max, e, ) ( 1 1  (16)
where n is the number of experiments and norm max, e, R  the mean of the R e,max, norm values. Fig. 12 shows the theoretical and experimental R e,norm values as a function of the sphere diameter. It must noticed that the experimental errors are close (about 4-5%) for all sphere diameters, which means that the error was not due to the accuracy of the ER measurements, since otherwise it would be more important for the small spheres. This error comes from the fact that the spheres did not always fall exactly on the vertical line x = y = 0 between the two electrodes. All experimental R e,max, norm values are in close agreement with the theoretical values, as shown in Fig. 12.

Conclusions

This paper investigates the influence of an insulating sphere on the electrolyte resistance between two electrodes immerged in an electrolytic solution, to assess the method as a singleparticle sizing technique at millimeter scale. A theoretical model using the finite element method was developed to study the influence of different parameters such as the size and the position of the sphere, the diameter of the electrodes and the interelectrode distance. Experiments were performed with different insulating spheres using a home-made electronic device. The excellent agreement between the experimental and theoretical results allows validation of both the theoretical model and the experimental setup able to measure fast electrolyte resistance variations in the millisecond range. For the moment, the technique allows particle sizing for spheres of diameter larger than onetenth of the interelectrode distance. Improvement may be obtained by reducing the background noise of the ER measurement. Apart from particle sizing, this specific electrochemical noise technique is able to perform accurate measurements of the sphere velocity without the help of optical techniques, which could be interesting in some applications as for viscosity measurements with falling balls. Only results with insulating spheres have been shown in this paper but the technique also works well with metallic spheres as long as there is a significant difference between the liquid and sphere conductivities. The development of a microfluidic platform is in progress in our laboratory to extend this technique to the micrometer scale and analyze smaller particles.

Figure and Table Captions

Table 1 Parameters obtained by fitting the experimental curves in Fig. 11 with Eq. 9 and parameters of the hydrodynamic flow. with Eq. 11 and fitted parameters given in Table 1. as a function of the sphere position along the z axis. Theoretical curves (red) calculated with Eq. 9 and fitted parameters given in Table 1.

Fig. 6a shows

  Fig.6ashows the ER values as a function of the interelectrode distance for spheres of

  for D = 20 mm and d e = 5 mm, which are the values used in the experimental work, the following expression with parameters obtained by fitting the R e values presented in Fig. 3 (d p = 6 mm) gives an excellent result:

  Fig.8b shows an excellent agreement between the R e values calculated for d p = 3 mm at different x and z values and the solid red lines obeying the following expression with a value of 3.05 for the exponent of d p .

Fig. 10

 10 Fig. 10 shows three typical examples of the variations of the output voltage signal v R e due to

Fig. 1 .

 1 Fig. 1. Geometry of the model (a) and example of meshing (b). D: interelectrode distance, d p :

Fig. 2 .

 2 Fig. 2. Distribution of potential (a) and current (b) lines in the y-z plane (x = 0) in the presence of

Fig. 3 .

 3 Fig. 3. Variations of R e and R e as a function of the position of the centre of the insulating sphere

Fig. 4 .

 4 Fig. 4. Variations of R e and R e as a function of the position of the centre of the insulating sphere

Fig. 5 .

 5 Fig. 5. Normalized ER increment R e,norm due to the presence of an insulating sphere of diameter

Fig. 6 .

 6 Fig. 6. Variations of R e (a) and R e,norm (b) as a function of the interelectrode distance for

Fig. 7 .

 7 Fig. 7. Variations of R e and R e,norm as a function of the electrode diameter for an insulating

Fig. 8 .

 8 Fig. 8. Values of R e as a function of the position of the centre of the insulating sphere along the

Fig. 9 .

 9 Fig. 9. Experimental setup used for measuring the electrolyte resistance fluctuations due to the

Fig. 10 .

 10 Fig. 10. v R e variations due to the passage of single insulating spheres of diameter 10 mm, 6 mm,

Fig. 11 .

 11 Fig. 11. Experimental (black) R e variations corresponding to the v R e variations shown in Fig. 10

Fig. 12 .

 12 Fig. 12. Experimental and theoretical and values of R e,max,norm as a function of the diameter of the

Fig. 1 .

 1 Fig. 1. Geometry of the model (a) and example of meshing (b). D: interelectrode distance,

Fig. 2 .

 2 Fig. 2. Distribution of potential (a) and current (b) lines in the y-z plane (x = 0) in the presence of

Fig. 3 .

 3 Fig. 3. Variations of R e and R e as a function of the position of the centre of the insulating sphere

Fig. 4 .

 4 Fig. 4. Variations of R e and R e as a function of the position of the centre of the insulating sphere

Fig. 5 .Fig. 6 .

 56 Fig. 5. Normalized ER increment R e,norm due to the presence of an insulating sphere of diameter

Fig. 7 .

 7 Fig. 7. Variations of R e and R e,norm as a function of the electrode diameter for an insulating

Fig. 8 .Fig. 9 .

 89 Fig. 8. Values of R e as a function of the position of the centre of the insulating sphere along the

Fig. 10 .

 10 Fig. 10. v R e variations due to the passage of single insulating spheres of diameter 10 mm, 6 mm,

Fig. 11 .

 11 Fig. 11. Experimental (black) R e variations corresponding to the v R e variations shown in Fig. 10

Fig. 12 .

 12 Fig. 12. Experimental and theoretical and values of R e,max,,norm as a function of the diameter of the

Table 1

 1 Parameters obtained by fitting the experimental curves in Fig.11with Eq. 9 and parameters of the

	hydrodynamic flow.
	Sphere diameter (mm)	10	6	3
	x (mm)			0.0867	0.142	0.126
	t 0 (s)			3.9738	2.7014	27.799
	v (mm s -1 )			632.66	520.53	324.37
	d p (mm)			10.488	6.033	3.128
	v	/	sphere	radius	(mm 1/2 s -1 )	282.9	300.5	264.8
	Reynolds number, Eq. 13	6580	3250	1010
	C d , Eq. 15			0.477	0.511	0.607
	v ter (mm s -1 ), Eq. 14	876	656	425