J M López-Herrera

A M Gañán-Calvo

S Popinet

M A Herrada
email: herrada@us.es

Electrokinetic effects in the breakup of electrified jets: a Volume-Of-Fluid numerical study

Keywords: VoF, electrohydrodynamics, electrokinetics, charge conservation, two phase flows

The breakup of a charged liquid column is studied numerically using Volume-OF-Fluid (VOF) for a range of timescales where electrokinetic phenomena may become significant, i.e when the time to breakup becomes comparable or shorter than the diffusion and the electrosmotic migration times of charged species. Here we propose a conservative method to deal with the diffusion of a tracer in VOF schemes when the diffusion is limited to one of the phases. The method consists in weighing the diffusivity with the value of the volume fraction computed from the analytically reconstructed interface. In this way, the interface is made impermeable to the tracer, which is conservatively kept within one of the phases. The performance of this method is first tested by comparing simple configurations with existing analytical solutions. In the cases when the diffusion, electrosmotic motion and hydrodynamic singularities compete, the results indicate that, after breakup, charges distribute between droplets differently from models assuming homogeneous and constant electrical conductivities (i.e. no electrokinetic effects). However, such departure does not alter the main hydrodynamic balances leading to wellestablished scaling laws of breakup.

Introduction

Electrokinetic effects in liquids determine bulk charge distributions when existing ionic species in solution respond to applied electric fields, but may also contribute to the global mechanical behavior of the system when free surfaces and interfaces are present. For ordinary values of surface tension, these effects naturally take place when the length scales of the system are below the millimetric scale, and particularly at the micro-and nano-scale in general microfluidic systems.

Among these systems, electrospray is probably the most studied and exploited natural example of global hydrodynamic consequences of electrokinetic effects in the presence of free surfaces. The extensive literature on the physics and biochemical analysis applications of electrospray amounts to more than 10 5 papers and a vast, complex network of citations. However, as puzzling as it may be, this network does not necessarily reflect a complete or sufficient scientific understanding. This may be a consequence of the fast growth rate of this publication network, compared to the average rate of general scientific knowledge assimilation.

In particular, while the basics of electrokinetics were established early, a detailed analysis of relevant publications up to the present day reveals striking understanding gaps in the real sequence of electro-physical processes taking place at the smallest scales: not only in cone-jet electrospray but also in many other phenomena such as the breakup of charged capillary liquid jets. These sequences determine the macroscopic outcome of electrospray in terms of issued charges per unit time and characteristic length scales of liquid emissions (droplets or particles). Interestingly, one may also observe how basic inconsistencies and customary assumptions become fossilized in the foundations of an ample literature.

Fernández de la Mora & Loscertales [START_REF] De La Mora | The current emitted by highly conducting Taylor cones[END_REF] postulated a set of scaling laws for the electric current and characteristic scales of liquid emission in the form of a jet issued from Taylor cones. These scaling laws were derived from the assumption that electrokinetic migration, or free charge relaxation towards the surface, was halted (or "frozen") at the apex of the cone as the jet scale was reached. A more relaxed version assuming that the jet scale emerged where the relaxation times of free charges became comparable to hydrodynamic residence times led to identical results. Even earlier, Gañán-Calvo, et al. [START_REF] Gañán-Calvo | The electrodynamics of electrified conical menisci[END_REF] suggested the opposite assumption (i.e. short electrical relaxation times compared to hydrodynamic ones) to reach alternative scaling laws, subsequently revisited in [START_REF] Gañán-Calvo | Cone-jet analytical extension of Taylor's electrostatic solution and the asymptotic universal scaling laws in electrospraying[END_REF][START_REF] Gañán-Calvo | The surface charge in electrospraying: its nature and its universal scaling laws[END_REF][START_REF] Gañán-Calvo | On the general scaling theory for electrospraying[END_REF][START_REF] Gañán-Calvo | Revision of capillary cone-jet physics: Electrospray and flow focusing[END_REF] and numerically confirmed in [START_REF] Hartman | Electrohydrodynamic atomization in the cone-jet mode physical modeling of the liquid cone and jet[END_REF]. Both extreme assumptions and corresponding models recognize electrokinetic effects as the ultimate reason for the appearance of driving forces deriving from Maxwell bulk stresses in the presence of interfaces. However, their extreme nature, the striking relative proximity of their results, and the inextricably indirect way to experimentally verify their validity have hardly been of help to build true and deep scientific knowledge in the community. In particular, the early introduction in the former model [START_REF] De La Mora | The current emitted by highly conducting Taylor cones[END_REF][START_REF] De La Mora | The fluid dynamics of Taylor cones[END_REF] of a fitting function f (ε), where ε is the electric permittivity of the liquid relative to vacuum, helped many subsequent authors to fit this model to their experimental results (e.g. [START_REF] Chen | Electrospraying of conducting liquids for monodisperse aerosol generation in the 4 nm to 1.8 µm diameter range[END_REF][START_REF] Chen | Experimental investigation of scaling laws for electrospraying: Dielectric constant effect[END_REF] among hundreds of works).

The rapid evolution of computational power and the increase in efficiency and precision of numerical schemes and methods have paved the way to the widespread but bold idea of tackling scientific conundrums like the one above via numerical simulation. In reality, the physics of fluid motions at the microscopic scale exhibits many features making their study particularly appealing to numerical modeling and simulation. In general, the small scale generally characterizing fluid motions in microfluidic systems limits the relative weight of convective effects compared to diffusion. This feature, reflected in moderate to small Reynolds number values, provides the adequate traits for full numerical simulation, where a high predictive power and accuracy has already been demonstrated. In general, the Lab-On-Chip (LOC) research community has a background in biology-related issues and are less familiar with engineering aspects such as numerical simulation and its advantages. Indeed, Boy et al. [START_REF] Boy | Simulation tools for lab on a chip research: advantages, challenges, and thoughts for the future[END_REF] wrote that their focus article on available computational methods for LOC systems, could serve "... to convince the LOC community that computation is a valuable tool and should be increasingly used over the next decade ...". In particular, numerical simulations allow researchers to determine in a rapid step how a design decision can affect the performance of a particular device. This way the development cost would drastically drop by reducing the number of prototyping iterations [START_REF] Boy | Simulation tools for lab on a chip research: advantages, challenges, and thoughts for the future[END_REF].

More recently, Wörner [START_REF] Wörner | Numerical modeling of multiphase flows in microfluidics and micro process engineering: a review of methods and applications[END_REF] exhaustively classified and described the foundations of the diverse numerical methods for two-phase flows and performed a complete review on the state of the art of numerical procedures to deal with challenging problems such as moving boundaries, Marangoni effects and surfactants, or heat and mass transfer across the interfaces. In spite of the strong predictive capabilities developed in this field, two-phase micro-flow problems involving the behavior of ionic species under the action of electric fields require an even deeper degree of physical insight. Contrary to what might initially appear, these problems are ubiquitous in fields handling several fluid phases at microscopic scales in predictive and consistent ways, which range from modern chemical engineering, biophysics, pharmaceutical research, to modern food processing, to name a few. Several instruments based on electrokinetic phenomena like the ζ-potential have even been developed [START_REF] Kirby | Micro-and Nanoscale Fluid Mechanics: Transport in Microfluidic Devices[END_REF]. However, those problems have traditionally been tackled following drastic electrokinetic simplifications that either assume (i) complete relaxation of all free charges at free surfaces, where bulks are neutral with homogeneous electrical conductivities (i. e., the leaky dielectric model of Melcher & Taylor [START_REF] Melcher | Electrohydrodynamics: A review of the role of interfacial shear stresses[END_REF][START_REF] Saville | Electrohydrodynamics: The taylor-melcher leaky dielectric model[END_REF], which entails having hydrodynamic times long compared to electrical relaxation), or (ii) the other extreme case where the liquids are assumed dielectric [START_REF] O'konski | The distortion of aerosol droplets by an electric field[END_REF][START_REF] Allan | Particle behaviour in shear and electric fields. i. deformation and burst of fluid drops[END_REF]. However, following the rationale of basic electrokinetics, the liquid electrical conductivity [START_REF] Saville | Electrohydrodynamics: The taylor-melcher leaky dielectric model[END_REF] can no longer be considered a homogeneously distributed value in the liquid bulk when charge relaxation is compromised by hydrodynamic motion. Therefore, the customary assumption of a constant liquid conductivity and liquid bulk electric neutrality would be inconsistent if one aims to compare numerical results based on that assumption with scaling laws such as the one in [START_REF] De La Mora | The current emitted by highly conducting Taylor cones[END_REF][START_REF] De La Mora | The fluid dynamics of Taylor cones[END_REF].

Moreover, the electrokinetic phenomena appearing in many electrohydrodynamic problems in the microfluidics field are strongly related to the presence of the electric double layers (EDL) that appear on interfaces when they are brought into contact with electrolytes [START_REF] Kirby | Micro-and Nanoscale Fluid Mechanics: Transport in Microfluidic Devices[END_REF][START_REF] Zhao | Advances in electrokinetics and their applications in micro/nano fluidics[END_REF]. The thickness of the EDL, λ D , is in the nanometric scale and can be very different from the characteristic length L o of the system being investigated. For example, in electrosmotic pumping the characteristic width of the impulsion channel is typically of the order of dozens of microns. In contrast, the EDLs present at the channel walls are nanometric. The ionic channels present in biological membranes are also of nanometric size [START_REF] Schoch | Transport phenomena in nanofluidics[END_REF][START_REF] Zheng | Second-order poisson-nernst-planck solver for ion transport[END_REF]. Consequently, very different numerical approaches have been developed depending on the particular electrokinetic system considered. In the case where the EDL can be assumed thin, λ D L o , a detailed EDL resolution can be avoided by either using the Helmholtz-Smoluchowski slip velocity or by using techniques like matched asymptotic expansions [START_REF] Squires | Induced-charge electro-osmosis[END_REF]. At the other extreme, for λ D ≥ L o , the electrokinetic systems are best described by considering all the individual atomic interactions, for example using molecular dynamics (MD) simulations [START_REF] Eijkel | Nanofluidics: what is it and what can we expect from it?[END_REF].

Between the aforementioned limits, the continuum approach in which ions are not treated as microscopic discrete entities but as continuous charged species densities [START_REF] Zheng | Second-order poisson-nernst-planck solver for ion transport[END_REF] is the proper choice.

Some of the ways used to manipulate droplets and bubbles in fluidic microsystems have their origin in electrokinetic phenomena like electrosmosis or electrophoresis [START_REF] Stone | Engineering flows in small devices: Microfluidics toward a lab-on-a-chip[END_REF]. For instance, a dielectric fluid of negligible conductivity like an oil can be continuously pumped through a channel by electrosmotic means by adding a suitable layer on an electrolyte [START_REF] Lee | Electrokinetic flow in a free surface-guided microchannel[END_REF][START_REF] Gao | Transient two-liquid electroosmotic flow with electric charges at the interface[END_REF]. Hence, investigations on the physics of electrokinetic effects on fluid-fluid interfaces such as the recent work of Pascall & Squires [START_REF] Pascall | Electrokinetics at liquid/liquid interfaces[END_REF] are of great interest. This work explains the physical reasons by which electrokinetics effects are enhanced at liquid/liquid interfaces. In this context, the work of Zholkovskij et al. [START_REF] Zholkovskij | An electrokinetic model of drop deformation in an electric field[END_REF] sheds light on the deformation of suspended droplets under an imposed axial electric field. The model of Zholkovskij et al. keeps the ionic nature of the charge and provides an analytical solution for the deformation. Interestingly, they show how the classic expressions obtained by Allan & Mason [START_REF] Allan | Particle behaviour in shear and electric fields. i. deformation and burst of fluid drops[END_REF] for pure dielectric fluids and by Taylor [START_REF] Taylor | Studies in electrohydrodynamics. i. the circulation produced in a drop by an electric field[END_REF] using the leaky-dielectric model are limit cases of the more general electrokinetic model. In the same manner, the work of Zholkovskij et al. can be used as a very complete benchmark for testing numerical models of two-phase electrokinetic problems [START_REF] Berry | A multiphase electrokinetic flow model for electrolytes with liquid/liquid interfaces[END_REF].

In the present paper, we will focus on the numerical treatment of twophase problems involving physical phenomena of electrokinetic nature using a Volume-Of-Fluid (VOF) method. In particular, our contribution here aims at the development of an accurate, efficient tool to deal with these general electrokinetic problems where the characteristic times associated to either the electrokinetics or the hydrodynamics can be comparable, with the ultimate objective to tackle the existing conundrum in electrospray physics and the electrohydrodynamic emission of extremely small charged droplets and particles from Taylor cone-jets. To this end, we use as foundation the open source code Gerris [START_REF] Popinet | The Gerris flow solver[END_REF]. Gerris, originally conceived as an incompressible Navier-Stokes equations solver [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF], also includes adaptive mesh capabilities, an accurate surface tension model [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF] and an electrohydrodynamic (EHD) solver [START_REF] López-Herrera | A charge-conservative approach for simulating electrohydrodynamic two-phase flows using volume-of-fluid[END_REF]. Thanks to its versatility, ease of use, free project character and accuracy, its uses are widespread. In particular it has been employed with remarkable success for simulating EHD problems such as electrospraying in the cone-jet mode under the hypothesis of constant conductivity (bulk ionic equilibrium) [START_REF] Herrada | Numerical simulation of electrospray in the cone-jet mode[END_REF] and tip streaming ejection from an electrified pendant drop [START_REF] Ferrera | Dynamical behavior of electrified pendant drops[END_REF].

Equations

For the sake of simplicity, in this work we assume that the electrolytes are fully dissociated, i.e. there are no bulk sources/sinks of ions. The conservation equations of the ionic species then take the form,

c i t + ∇ • (c i u) = ∇ • ω i k B T ∇c i -eω i z i c i E (1)
with c i the concentration (number of ions per unit volume) of the ionic ispecies, u the fluid velocity, ω i the mobility of that ion, e the elementary electron charge, z i the valence (with its sign) of that ionic specie, k B the Boltzmann constant, T the temperature and E the electric field. The above equation ( 1) expresses that the concentration of the i-species varies in time as a consequence of advection, migration under an electric field (last term on the r.h.s.) and transport by diffusion (first term on the r.h.s.) [START_REF] Saville | Electrohydrodynamics: The taylor-melcher leaky dielectric model[END_REF]. The electrical potential ,ϕ, depends on the distribution of all the charged species by means of the Poisson equation,

∇ • (εE) = ∇ • (-ε∇ϕ) = q, since q = i ez i c i (2)
where ε is the electric permittivity and q the volume charge density.

The set of equations given by Eqs. ( 1) and ( 2) leads to the Poisson-Nernst-Planck (PNP) model. This is precisely the model adopted in this work. In addition to the PNP model, the Navier-Stokes equations for the incompressible fluid motion need to be included,

∇ • u = 0, (3)
ρ ∂u ∂t + u • ∇u = -∇p + ∇ • T v + F e + σκδ s n ( 4
)
where ρ is the fluid density, σ the surface tension coefficient, κ the interface curvature, n the normal to the interface. δ s is the Dirac delta and T v is the viscous stress tensor given by,

T v = 2µD , ( 5
)
where µ is the viscosity and D the deformation tensor, D = 1 2 (∇u + ∇u T ). The bulk electric forces F e can be derived from the electrostatic Maxwell stress tensor

T e = ε EE - E 2 2 I (6)
6 by applying the divergence operator

F e = ∇ • T e = ρ e E - 1 2 E 2 ∇ε. (7)
The first term represents the electric forces exerted on the free charges present in the fluid, while the second term represents the electric forces exerted on the induced electric dipoles.

The set of equations is completed, in a VOF approach, by incorporating an additional variable, the volume fraction φ(x, t), that serves to track the interface position

∂φ ∂t + ∇ • (uφ) = 0 (8)
Thus, the entire two-phase fluid domain, formed by a fluid "o" immiscible with a fluid "e", is treated as a single fluid with properties that depend on the position (through the volume fraction φ) at each instant,

χ = χ o φ + χ e (1 -φ) , (9)
where χ stands for any of the relevant fluid properties, ρ, µ, ε, ...

Numerical scheme

The fluid domain is, for simplicity, formed by two immiscible fluids with homogeneous properties: a perfect dielectric whose properties are labeled with the subscript "o" and a fully dissociated z : z binary electrolyte solution (subscript "e"). In addition, variables and equations are made dimensionless using the fluid density ρ e , the bulk concentration c o , a characteristic length L, the surface tension σ and the permittivity ε o . The above scaling leads to the following dimensionless parameters where we will use the superscripts + and -for the cation and the anion, respectively:

• Dimensionless ion diffusivities, D +(-) = ω +(-) k B T LU
with U the capillary velocity U = (σ/ρ e L) 1/2 . Note that these dimensionless diffusivities are the inverse of the corresponding Peclet numbers P e +(-) = 1/D +(-) .

• A number γ, which measure the relative importance of the characteristic electric field to the one created by electrokinetics,

γ = E c Lez/(k B T ), with E c = (σ/ε o L) 1/2
. Dimensionless ion specific conductivities can be written as, Λ +(-) = D +(-) γ. From now on, all variables will be non-dimensional (but we will use the same symbols as in the previous section).

Interface C c cC u f Dt u f c f W F=0 L E =L E F=L E ; D E =D E F=D E F=1 Migration Flux (a) (b) L E =L E F=0 ; D E =D E F=0

Conservation of ionic species

Obviously, the conservation equation ( 1) is only applicable in regions of the domain where the free ions move, i.e. regions occupied by the liquid solvent. The concentration c i is defined as the amount (number of ions) of the ispecies per unit volume of solvent. Also, any concentration c i per unit of spatial volume can be expressed anywhere in the computational domain just by weighting it with the volume fraction of the solvent φ c i , which is more suitable for numerical methods treating two-phase flows as a single fluid with a variable property φ. To this end, Berry et al. [START_REF] Berry | A multiphase electrokinetic flow model for electrolytes with liquid/liquid interfaces[END_REF] combined (1) and (8

) to derive (c ± φ) t + ∇ • (c ± φu) = ∇ • φD ± ∇c ± ∓ ∇ • Λ ± φc ± E ( 10
)
where the interface is assumed impenetrable to ions. A physical interpretation of the factor φ in the migration terms of the r.h.s. of ( 10) is as a weighted diffusivity/conductivity ( D± = φD ± / Λ± = φΛ ± ) that nullifies the migration fluxes across boundaries out of the solvent phase (see sketch (a) of figure 1). In the present case, it is used as the control volume for the spatial integration in cell C (see figure 1, sketch (b)). The flux of species due to advection is computed with a procedure analogous to the one used for the volume of fluid phase [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF]. A similar procedure has been adopted in [START_REF] Berry | A multiphase electrokinetic flow model for electrolytes with liquid/liquid interfaces[END_REF]. As is sketched in figure 1.b, the procedure relies on computing the volume of fluid crossing the cell frontier ∂C (the dark grey area in the sketch) from the analytically reconstructed interface. Then, in order to calculate the amount of the species that leaves the cell, that volume of fluid simply has to be weighted by the value of concentration at the cell face (c ± ) f . These face values are calculated from the cell concentration and slope-limited concentration gradients.

In the spirit of using variables defined in the entire domain, the diffusivity term in [START_REF] Chen | Experimental investigation of scaling laws for electrospraying: Dielectric constant effect[END_REF]

is split in two, ∇ • (φD ± ∇c ± ) = ∇ • [D ± ∇(c ± φ)] -∇ • (D ± c ± ∇φ) , (11)
which gives the following form of the equation for the species

(c ± φ) t + ∇ • (c ± φu) = ∇ • [D ± ∇(c ± φ)] term A -∇ • (D ± c ± ∇φ) term B ∓ ∇ • Λ ± φc ± E term C , (12)
where the diffusivity and conductivity are now defined in the entire fluid domain and are treated like the other fluid properties.

In Gerris, the space is discretised using an octree scheme where the unknown variables are located at the center of each cubic discretisation volume, and are interpreted as the average values of the variable in the cell. In this octree scheme, the degree of refinement of the domain is often referred by its level. Level zero correspond to a unique cell occupying the entire square domain. Each increment of the level corresponds to a generation of new 4 (8 in 3D) cells by splitting the previous (parent) once. Thus, the cell's width depends on the level, L, as h = 2 -L .

In a cell of size h, the average flux added/extracted by pure diffusion (term A) can be computed as [START_REF] Kirby | Micro-and Nanoscale Fluid Mechanics: Transport in Microfluidic Devices[END_REF] where ∇ f (c ± φ) is the normal gradient at the cell faces computed from center values at the cell of interest and its neighbors [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF]. (D ± ) f is the diffusivity value at the cell face. We compute it by averaging the diffusivity of the cells sharing the face of interest and by weighting this value by the factor φ f , (D ± ) f = D ± av φ f . The weighting factor φ f is the value of the volume fraction at the face, computed from the analytically reconstructed interface. This weighting factor represents the ratio of the face wetted by the the solvent phase (for example, in the sketch 1.b, for the upper face φ f = 0 while for the lower face φ f = 1). The weighting factor ensures that interfaces are impenetrable to ions in a way similar to solid fractions for regions occupied by solids [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF]. In this way, the concentrations will remain conservatively in the solvent phase. The correction term B of (12) has a diffusion-like structure and is computed as term A. The same holds for term C since we compute it from the electric potential ϕ, ∇ • (Λ ± φc ± E) = -∇ • (Λ ± φc ± ∇ϕ). An alternative approach to the calculation of term C could be to assimilate it into the advection term [START_REF] Berry | A multiphase electrokinetic flow model for electrolytes with liquid/liquid interfaces[END_REF].

h C ∇ • [D ± ∇(c ± φ)] = ∂C D ± ∇(c ± φ) • n = f (D ± ) f ∇ f (c ± φ),

Time integration procedure

The time discretization scheme consists in a time-splitting pressure-correction method. The time stepping integration procedure is briefly outlined below (readers can find more detailed descriptions elsewhere [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF][START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF][START_REF] Lagree | The granular column collapse as a continuum: Validity of a two-dimensional navier-stokes model with a µ(i)-rheology[END_REF]). First, anion concentrations and the volume fraction are advanced to a mid-step,

n + 1/2, φ n+ 1 2 -φ n-1 2 ∆t + ∇ • (φ n u n ) = 0 (14) c ± n+ 1 2 -c ± n-1 2 ∆t + ∇ • (c ± n u n ) = ∇ • D ± ∇c ± n+ 1 2 ∓ Λ ± c ± n E n (15)
Then the values of the fluid properties are updated,

ρ n+ 1 2 = φ n+ 1 2 + R(1 -φ n+ 1 2 ) µ n+ 1 2 = C µ [ φ n+ 1 2 + M (1 -φ n+ 1 2 )] ε n+ 1 2 = S φ n+ 1 2 + (1 -φ n+ 1 2
), [START_REF] O'konski | The distortion of aerosol droplets by an electric field[END_REF] as well as the electric potential,

∇ • (ε n+ 1 2 ∇ϕ n+ 1 2 ) = 1 2 SK 2 γ (c + n+ 1 2 -c - n+ 1 2
) .

The prediction-diffusion step is performed by solving the equation

ρ n+ 1 2 ∆t u * -∇ • (µ n+ 1 2 D * ) = ∇ • (µ n+ 1 2 D n ) + (σκδ s n) n+ 1 2 + (F e ) n+ 1 2 + ρ n+ 1 2 u n ∆t -u n+ 1 2 • ∇u n+ 1 2 , (18)
to determine the auxiliary velocity ũ * . In the above expression, the velocity advection term

u n+ 1 2 • ∇u n+ 1 2
is estimated by means of the Bell-Colella-Glaz second-order unsplit upwind scheme [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF]. The projection-correction step is then carried out by solving the Poisson equation,

∇ • ∆t ρ n+ 1 2 ∇p n+ 1 2 = ∇ • u * , (19)
and by numerically computing the divergence-free velocity field for the new instant n + 1, u n+1 = u * -∇p n+1/2 ∆t/ρ n+1/2 . The surface tension forces are computed using the Continuum-Surface-Force (CSF) approach [START_REF] Zemach | A continuum method for modeling surface tension[END_REF]. It is well known that the CSF approach can cause parasitic currents. However, it is possible to avoid them by using a balancedforce description of the surface tension and pressure gradient together with an accurate curvature estimate [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF]. The curvature is calculated using a generalised height-function technique which allows consistent and accurate estimations even at low interface resolutions.

The time integration scheme is explicit with a timestep limited by the onset of capillary, advection or diffusion instability. The most stringent limitation depends on the parameters of the problem. Note that viscosity does not appear in the list above since the viscous term is calculated implicitly.

Results and discussion

Deformation of suspended electrolytic droplets

An uncharged liquid z : z electrolyte droplet of radius a is suspended in a pure dielectric unbounded liquid atmosphere. Both fluids are immiscible with surface tension σ, and, for the sake of simplicity, have the same density and viscosity. The droplet, by the action of an imposed axial electric field, E ∞ , deforms to adopt a stationary prolate/oblate form depending on the ratio of the droplet radius to the Debye length, K = a/λ D and on the ratio of the inner to the outer electrical permittivity

S = ε i /ε o , d Ca E = 9 16 × (S -1) S [2 -K 2 /(K coth K -1)] 2 -1 + K 2 S [2(S -1) + SK 2 /(K coth K -1)] 2 (20)
where d is the degree of the deformation of the droplet given by d = (aa ⊥ )/(a + a ⊥ ), and a and a ⊥ are the semi-axes parallel and normal to the external field. Ca E is the electric capillary number, Ca E = aε o E 2 ∞ /σ. The above expression is valid in the limit of small Ca E since it has been obtained by linearization. Expressions for different configurations, for example a dielectric droplet surrounded by a z : z electrolytes, can be found in Zholkovskij et al. [START_REF] Zholkovskij | An electrokinetic model of drop deformation in an electric field[END_REF]. This problem has been proposed as a benchmark by Berry et al. [START_REF] Berry | A multiphase electrokinetic flow model for electrolytes with liquid/liquid interfaces[END_REF]. In these simulations, we have set to one the Ohnesorge number, C µ , and ion diffusivities D + and D -. The fluid domain is a square box of width W = 30a. In the computations we have used the symmetry of the problem; the lower side of the square box has been set as the axis of symmetry. On the other boundaries we impose slip conditions for the velocity. We have used adaptation to refine the cells close to the interface as it deforms. The smaller cells have been used at the interface while further away the size of the cell has been increased to level L = 4 (which is equivalent to cells of width h = W 2 -L = 1.875 a). We have first checked the convergence with the grid by refining the minimum cells successively to a/h = 8.53, 17.06 and 34.12 (or equivalently, levels L = 8, 9 and 10) for K = 0.1 and K = 0.6. In these preliminary tests we have set the permittivity ratio S = 10. The results of this test are summarized in table 4.1. At level L = 8, the grid is too coarse, with relative errors larger than 10%. Increasing the resolution by a level suffices to decrease the relative error by an order of magnitude, with Table 1: Convergence of the degree of deformation, d, with the dimensionless grid size (droplet radius a divided by the minimum cell size) with S = 10, Ca E = 0.025. Parameter K is K = 0.1 and K = 6.0. The theoretical result of (20), d t for these values of K are 0.007929 and 0.0134289, respectively. The ratio of the Debye length to the minimum cell size, λ D /h is also shown. a value close to the theoretical solution given by equation ( 20) for K = 6.0. However for K = 0.1 a refinement of one level (from 8 to 9) is not so effective. In this case, the relative error decreases approximately 3.6 times, instead of a decade. In both cases, little is gained with further refinement. We have also checked if a similar agreement could be obtained for different values of the permittivity ratio S and of the dimensionless Debye parameter K. To this end, we have plotted the droplet deformation d as a function of K for electrical permittivity ratios S = 2 and S = 10 (figure 2), keeping the maximum level equal to 9. The other parameters are kept constant (as previously given). Figure 2 shows that the agreement is good for K values of the order unity, but deteriorates for low K values. The total amount of either anions and cations is very well conserved in all the simulations: in the least conservative case, the variation of the total amount of species between the first and last computational step does not exceed 10 -3 %.

Grid (a

Breakup of liquid, charged capillary columns

Uncharged liquid columns are unstable to perturbations because of surface tension. Capillary instabilities grow until a pinch-off occurs and droplets (often a primary and a satellite) are formed. Due to industrial implications, this phenomenon has been extensively studied since the XIX century pioneering works of Savart [START_REF] Savart | Mémoires sur la constitution des veines liquides. lances par des orifices circulaires en mince paroi[END_REF], Plateau [START_REF] Plateau | Statique experimentale et theorique des liquides soumis aux seules forces moleculaires[END_REF] and Rayleigh [START_REF] Rayleigh | On the stability of liquid jet[END_REF]. Much work has been conducted to analyze for instance, the influence of the viscosity or the surrounding ambient fluid on the size of the droplets or on the dynamics of the pinch-off process. The reader is referred to [START_REF] Eggers | Nonlinear dynamics and breakup of free-surface flows[END_REF][START_REF] Eggers | Physics of liquid jets[END_REF] for a complete overview of the state-of-the-art.

The problem is enriched if electric forces are added. A wide variety of electrical conditions emerge: the fluid could behave, for example, as a perfect conductor, a perfect dielectric, or in a milder situation, as a leaky-dielectric fluid; or some external electric field could be superimposed. Thus, a myriad of papers on these electrified ligaments can be found. Most of them are devoted to linear stability analysis [START_REF] Saville | Electrohydrodynamic stability: Fluid cylinders in longitudinal electric fields[END_REF][START_REF] Saville | Stability of electrically charged viscous cylinders[END_REF][START_REF] Mestel | Electrohydrodynamic stability of a slightly viscous jet[END_REF][START_REF] Mestel | Electrohydrodynamic stability of a highly viscous jet[END_REF][START_REF] López-Herrera | Linear stability analysis of axisymmetric perturbations in imperfectly conducting liquid jets[END_REF][START_REF] Li | Linear instability of a coflowing jet under an axial electric field[END_REF][START_REF] López-Herrera | Absolute to convective instability transition in charged liquid jets[END_REF], some to numerical non-linear analysis [START_REF] Setiawan | Nonlinear modeling of an infinite electrified jet[END_REF][START_REF] López-Herrera | One-dimensional simulation of the breakup of capillary jets of conducting liquids. application to e.h.d. spraying[END_REF][START_REF] Collins | Breakup of electrified jets[END_REF][START_REF] Wang | Breakup of a poorly conducting liquid thread subject to a radial electric field at zero reynolds number[END_REF] and a few to experimental analysis [START_REF] López-Herrera | A note on charged capillary jet breakup of conducting liquids: Experimental validation of a viscous one-dimensional model[END_REF][START_REF] Zhakin | Experimental study of the outflow of charged drops and jets[END_REF]. To our knowledge, the work of Conroy et al. [START_REF] Conroy | Breakup of an electrified viscous thread with charged surfactants[END_REF] is the only one considering electrokinetic effects in the breakup of charged threads. This work focuses on the effect of the presence at the interface of a positive insoluble surfactant interacting with ionic species. A slenderness hypothesis and the Debye-Huckel limit are used to simplify the problem.

In the present work, a net charge is induced in a slender, perfectly cylindrical column by applying a difference of voltage V between the column of diameter A (which is used as characteristic length in this problem) and a grounded concentric electrode of radius R ∞ = 15 A. The surrounding ambient is a perfect dielectric fluid with negligible dynamical effect on the column, i.e. a gas. The net charge is due to a small imbalance in the concentration of the fully dissociated z : z binary electrolyte, |c + -c -| = 0.

Since the present study is restricted to axisymmetric perturbations, we use cylindrical coordinates, (z, r; t). Initially, the concentrations of the charged species are uniformly distributed. Also, to trigger the breakup process, a sinusoidal, small perturbation of the interface is imposed,

f (z; 0) = 1 + sin(κ w z) c + (z, r; 0) = B + c -(z, r; 0) = B - ( 21
)
where f is the dimensionless interface position and and κ w are the amplitude and wavenumber of the perturbation respectively. B + and B -are the initial concentrations of the cations and the anions. As in the above problem of the deformation of the droplet, the degree of electrification is measured by the capillary electric number, Ca E , which in this case is defined as

Ca E = Aε o E 2
o /σ with E o the outer electric field at the interface of the column and ε o the permittivity of the surrounding medium.

Ca E is related to the average charge density induced in the fluid. Thus Ca E can be written in terms of the initial species concentration as

Ca E = K 4 S 2 16γ 2 (B + -B -) 2 . ( 22
)
Ca E is also called the Taylor number [START_REF] López-Herrera | Linear stability analysis of axisymmetric perturbations in imperfectly conducting liquid jets[END_REF].

To sum up, the following set of free dimensionless characteristic parameters govern the problem:

(1) for the perturbation, and κ w ;

(2) for the charged species, D + , D -, γ and K;

(3) for the electrical conditions, B + , B -and R ∞ /A; and (4) for the fluid properties, C µ , S, R and M .

Since we focus mainly on electrokinetic effects, most of the free parameters in this study are kept fixed: = 0.1, κ w = 0.6283, C µ = 0.05, R = M = 10 -2 and S = 10. Additionally, in order to have a more pronounced electrokinetic effect, we assume that the cation is of smaller size than the anion. Hence, we set a higher diffusivity for the cation than for the anion, D + = 7 and D -= 1. Also the difference of mobilities will allow to investigate the influence of the polarity in the breakup. In particular we set B + = 1.01 and B -= 0.99 to study the case of positive polarity, and we swap the values for the case of negative polarity, B + = 0.99 and B -= 1.01. The only free parameter we allow to vary is K. γ is calculated from [START_REF] Eijkel | Nanofluidics: what is it and what can we expect from it?[END_REF] with Ca E = 0.125. To fix Ca E rather than γ yields a proper comparison between the cases since the degree of electrification is the most affecting factor, after the Ohnesorge number, in the breakup process [START_REF] López-Herrera | One-dimensional simulation of the breakup of capillary jets of conducting liquids. application to e.h.d. spraying[END_REF].

In the simulations we use the axisymmetric character of the problem and the symmetry existing in the axial direction. Hence, the simulation domain occupies only half a wave length in the z direction and in the radial direction from the axis of symmetry up to the grounded electrode (see fig. 3). Consistently the following boundary conditions are used,

• At the axis of symmetry, u r (z, 0; t) = 0, v(z, 0; t) = 0 and ϕ r (z, 0; t) = 0 .

• At the left and right extremes, u(±π/(2κ w ), r; t) = 0, v z (±π/(2κ w ), r; t) = 0 and ϕ z (±π/(2κ w ), r; t) = 0 . • At the grounded electrode, u(z, R ∞ /A; t) = 0, v(z, R ∞ /A; t) = 0 and ϕ(z, R ∞ /A; t) = 0 .

with u and v the axial and radial velocity.

Simulations have been carried out using a double adaptation refinement criteria, based on the gradient of the volume fraction, φ, and on the maximum curvature of the interface κ max . The first criterion allows to guarantee that the interface is always defined with a minimum grid size of h/A 0.0097 (equivalent to 9 levels of refinement). The second criterion is introduced to get a good description of the pinching region as time proceeds. This second criterion ensures that the cell size ∆ is small enough to verify ∆κ max < 0.2. The levels are allowed to increase up to a level 14 (h/A 3.05 × 10 -4 ) before breakup. Once breakup has occurred the curvature criterion is relaxed by reducing the maximum level to 10.

In figure 4 we plot the amount of each species (expressed as the percentage of the initial one seeded in the column) that goes to the satellite droplet either for positive and negative polarity and for two pair of values of diffusivities, D+ = 7 × 10 -3 and D -= 10 -3 (green and blue lines) and D+ = 7 and D -= 1 (red and cyan lines). The first electrokinetic effect observed is that the symmetry with respect to the polarity set is broken, as expected: continuous lines corresponding to positive polarity are in all the cases above the dashed lines of negative polarity. The asymmetry is particularly intense for low diffusivities (D+ = 7 × 10 -3 and D -= 10 -3 ). Note that in this case (low diffusivities) even the amount of anion in the satellite is larger for positive polarity than for negative polarity (i.e. continuous-square line above the dashed-square one). For the high diffusivity pair (D+ = 7 and D -= 1) we can distinguish roughly three regions: region I (K < 2), region II ( 2 < K < 20) and region III (K > 20). In region I, the gap between the amount of anion and cation, which is proportional to the net charge, is greatly reduced as K is lowered. Electrokinetic effects become manifest when the polarity is switched: the charged species concentration in the satellite decays more quickly for negative polarity (cyan curve) than for positive polarity (red curve). For instance, for K = 0.5 the amount of charged species is about 2.8% for negative polarity, while it is about 3.3% for positive polarity. In region II, the amount of charged species is roughly a plateau, ranging from 3.3% for cations with negative polarity up to 3.7% for anions with positive polarity. Besides, region III is characterized by a continuous (almost linear) growth with K of the concentration of all charged species, irrespective of the polarity. For more realistic values of the diffusivity (D+ = 7 × 10 -3 and D -= 10 -3 ), region II can not be distinguished: the transition from regions I (nearly equal concentrations of anion and cation) to III (continuous growth of the amount of charge with K) takes place around K=5 without an intermediate plateau. It is worth noticing that small differences of ion concentration in the satellite (note that the y axis goes from 2.5% to 5%) yield big differences in its net charge (see fig. 7).

In figure 5 we show the dimensionless relative bulk conductivity at an instant before pinching for K = 0.5 ((a) and (b) subplot) and K = 20 ((c) and (d)). The polarity is positive and the dimensionless diffusivity pair (D + , D -) is (7.0, 1.0) in both cases. The conductivity shown has been weighted with the initial homogeneous one given by Eq. [START_REF] Gao | Transient two-liquid electroosmotic flow with electric charges at the interface[END_REF]. In both cases the relative conductivity is not homogeneous in the bulk, with the maximum and the minimum of the conductivity of the same order independently of the value of K; about 15% higher and 1.5% lower with respect to the initial weighting factor, respectively. The location of the maximum values is more interesting. For K = 20 the maximum values are located in the vicinity of the interface, where the net charge accumulates. In contrast, for K = 0.5, the maximum conductivity is located in the neck where the charges are accumulated by a relatively intense electric field. This intense electric field is a consequence of the sudden difference of electric potential that occurs through the neck region as indicated by the electric isopotential lines in figure 5(b).

In the pinch-off process of an uncharged jet, the minimum radius h min scales with time as h min ∝ (t o -t) n where t o is the breakup time [START_REF] Eggers | Physics of liquid jets[END_REF]. The exponent n depends on the relevant forces acting on the pinching process. Surface tension, viscosity and density establish a threshold value h * min such that the balance is between surface tension, viscous and inertia forces. While the instantaneous value of h min is such that h min > h * min , the balance is applied between surface tension forces and inertia, and the exponent n takes the value n = 2/3. As time approaches pinch-off and h min < h * min , viscous forces overcome inertia and the pinch-off evolves linearly (n = 1) [START_REF] Chen | Computational and experimental analysis of pinch-off and scaling[END_REF]. In figure 6 we plot the time evolution of the minimum radius for the same conditions used in figure 5. One can observe that the scalings of the last stages of the breakup are not altered by electrokinetic effects, independently of the width of the relative Debye length, measured by the dimensionless Debye parameter K. A similar behaviour is found in charged capillary jets in which electrokinetic effects are absent [START_REF] Collins | Breakup of electrified jets[END_REF].

Discussion and conclusion. The validity of the customary electrohydrodynamic, homogeneous conductivity assumption

In this section we finally analyze the validity of a general assumption made by many in the field of electrohydrodynamics and, in particular, by all investigators in the field of electrospray: the homogeneity of the electrical conductivity throughout the liquid domain. Thus, for comparison purposes we write down the equation for the charge density, q, that is obtained by combination of the equations of the cation and the anion,

q t + ∇ • (qu) = K 2 S 2γ ∇ • (D + ∇c + -D -∇c -) -∇ • D + c + + D -c - 2 K 2 SE (23)
It is customary in the electrohydrodynamic field to neglect the diffusion term since, in most cases, it is negligible compared to the electrical migration. In 24) is also plotted (dash-dot line). Thick lines for D + = 7 × 10 -3 , D -= 10 -3 , thin lines for D + = 7, D -= 1 addition, the concentrations of species tend to remain almost uniform in the bulk. It is thus sensible (and customary) to assume a constant electric conductivity. With the above characteristic simplifications [START_REF] Stone | Engineering flows in small devices: Microfluidics toward a lab-on-a-chip[END_REF] becomes

q t + ∇ • (qu) = -∇ • (αSE) , ( 24
)
where α is the relaxation parameter [START_REF] López-Herrera | Linear stability analysis of axisymmetric perturbations in imperfectly conducting liquid jets[END_REF] that measures the relative importance of charge conduction to advection. For α 1 the perfect conductor limit is reached. On the other hand, if α 1, the "glued charge" limit is attained [START_REF] López-Herrera | Linear stability analysis of axisymmetric perturbations in imperfectly conducting liquid jets[END_REF]. Comparing [START_REF] Stone | Engineering flows in small devices: Microfluidics toward a lab-on-a-chip[END_REF] and [START_REF] Lee | Electrokinetic flow in a free surface-guided microchannel[END_REF] the following relationship between dimensionless parameters arise,

α = D + B + + D -B - 2 K 2 . ( 25
)
In figure 7 we compare the results with no assumptions on the electrical conductivity (i.e., we let each species move according to their mobility under the applied electric field) and using simplifying equations 24 and 25 with homogeneous electrical conductivity (i. e. that of the liquid in the absence of electrokinetic effects). Thus, we plot the net charge in the satellite droplet (again, as the percentage of the initial amount) either for positive and negative polarity from the general model. It is clear from figure 7 that the effect of polarity is negligible in the net charge in the droplets. It can be observed that the amount of charge in the satellite increases rapidly with K, reaching a limit value of about 20%. This is because the average conductivity is proportional to K 2 . With low values of the conductivity the main mechanism governing the movement of the charged species is convection, while the diffusion and the electrical migration are less important. Therefore, both species, the cation and the anion, are moved analogously, the concentrations are similar and the net charge is small. The presence of a stagnation point in the liquid bulk is also important since it acts like a barrier to the convection of the species. This stagnation point is a consequence of the outward pumping of the fluid occurring in the pinch-off area. On the other hand, the results of the simplified EHD model given by Eq. ( 24) for different values of α are also plotted in figure 7. K is calculated from α with Eq. [START_REF] Gao | Transient two-liquid electroosmotic flow with electric charges at the interface[END_REF]. Similar trends of lower net charge for low values of K are observed for the EHD model; however this charge is significantly overpredicted compared to that calculated with the electrokinetic model, with errors as large as 100% for low K values. These results are a warning for those using theoretical models in the field of electrospray physics in the limits of very small issued flow rates, and in particular when a first spout is issued from an electrified interface [START_REF] Collins | Electrohydrodynamic tip streaming and emission of charged drops from liquid cones[END_REF].

In conclusion, a general electrokinetic model and numerical procedure to tackle electrohydrodynamic problems has been presented with some illustrations of its validity and accuracy. In the cases when the diffusion, electrosmotic motion and hydrodynamic singularities compete, the general electrokinetic model yields results significantly different from customary models assuming homogeneous electrical conductivities with values equal to those of the liquids in the absence of electrokinetic effects. In particular, at the macroscopic level the net charge of the droplets could be overestimated specially if diffusion effects are not negligible. With the simulations (and free) tool we developed, a better insight on the distribution of charged species is gained. This is particularly valuable for the characterization of processes or techniques in which electrospray is operated with extremely small flow rates, such as in mass spectrometry. However, such differences are not sufficient to alter the already established asymptotic behavior of the minimum radius of the capillary liquid jet as the breakup time approaches.

Figure 1 :

 1 Figure 1: Sketch of an interfacial cell • The Ohnesorge number, C µ = µ e / √ σLρ e .

Figure 2 :

 2 Figure 2: Deformation d as a function of the Debye parameter K for permittivity ratios S = 2 and 10. Continuous lines correspond to the theoretical values given by (20). Symbols correspond to numerical simulations carried out with a minimum cell size equal to a/h = 17.06 (L=9). Calculations have been carried out setting C µ = D + = D -= 1 and Ca E = 0.025. The ratios of density and viscosity, R and M , are both set to one.

Figure 3 :

 3 Figure 3: Simulation domain and boundary conditions of the problem of the charged column breakup.

Figure 4 :

 4 Figure 4: Percentage of the anions and cations that go to the satellite droplet as a function of the Debye parameter, K. Rounded and square symbols denote the cation and anion species, respectively. Continuous (dashed) lines denote positive (negative) polarity results. Green and blue lines correspond to D + = 7 × 10 -3 , D -= 10 -3 , while red and cyan lines have been calculated using less realistic values D + = 7, D -= 1

Figure 5 :

 5 Figure 5: Dimensionless relative bulk conductivity, (c + D + + c -D -)/(B + D + + B -D -), for an instant before pinching. (a) and (b) corresponds to K = 0.5; (c) and (d) to K = 20. In (b) and (d) the pinch-off region is enlarged. The polarity is positive. The electric isopotential lines are also represented. Colorscale for (a) and (b) ranges from 1.1515 to 0.9886 while for (c) and (d) the range goes from 1.1277 to 0.9857. Dimensionless diffusivities are D + = 7, D -= 1 in all cases.

Figure 6 :

 6 Figure 6: Evolution of the minimum radius as the breakup time approaches for K = 20 and K = 0.5. The figure also shows the scalings h min ∝ (t o -t) 2/3 (blue line) and h min ∝ (t o -t) (magenta line) for uncharged jets.

Figure 7 :

 7 Figure 7: Percentage of the net charge going to the satellite droplet as a function of the Debye parameter, K. The net charge predicted by the EHD model given by Eq. (24) is also plotted (dash-dot line). Thick lines for D + = 7 × 10 -3 , D -= 10 -3 , thin lines for D + = 7, D -= 1

Acknowledgements

Partial support from the former Ministry of Science and Education and Junta de Andalucía (Spain) through Grants Nos. DPI2010-21103 and P08-TEP-04128, respectively, is gratefully acknowledged. We also thank the reviewers for their valuable help in improving the present work.


Christine Nazaret

Sophie Colombié

Bertrand Beauvoit

Martine Dieuaide-Noubhani

Jean-Pierre Mazat

Yves Gibon

How to optimize

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

The modelling of metabolic pathways is developed here for application to sugar metabolism of tomatoes during fruit expansion. In the now classic calculation of network fluxes, we add concentrations. Under assumptions of steady state and of "minimal effort" as proposed in [START_REF] Holzhütter | The principle of flux minimization and its application to estimate stationary fluxes in metabolic networks[END_REF], fluxes and concentrations, respecting constraints of experimental measurements and of thermodynamics, are searched. In order to obtain a well-posed problem, a preliminary mathematical study is in progress.

Backgrounds and Basics

The development time of tomato fruit is about 50 days. This period is subdivided in 50 daily phases numbered by l. At the scale of a phase (day), we suppose that the flux vector is at steady state and follows some fixed constraints. But at the scale of the development, these constraints will change. So we formulate a general problem at a phase scale l and solve it for l = 1 • • • 50.

Notations 1. X = (X 1 , • • • , X m ) the vector of internal metabolites, 2. V = (v 1 , • • • , v n ) the flux
vector composed by the rate v j of each reaction j in the network.

3. S = (s ij ) i=1•••m,j=1••
•n the stoichiometry matrix (only refers to the internal metabolites) where s ij is the stoichiometric coefficient of metabolite X i in reaction j.

S = (s

ij ) i=1•••m,j=1•••n i a submatrix
extracted from the stoichiometry matrix containing only the stoichiometric coefficients of internal reactions (reactions between internal metabolites).

Previous approach and results: problem n o 1

Constraints Steady-state (for all internal metabolites) assumption

dX/dt = SV = 0
Bounds on the fluxes (the values v l,j,min and v l,j,max are constants)

v l,j,min ≤ v j ≤ v l,j,max for j = 1 • • • n Proposition. min V ∈A l
||V || 2 admits a unique solution where A l is the non empty following feasible set

A l = {V ∈ R n , SV = 0, v l,j,min ≤ v j ≤ v l,j,max for j = 1 • • • n}

Thermodynamics Basics

Let us recall that the Gibbs free energy change for a reaction r i

n 1 X 1 + n 2 X 2 ←→ n 3 X 3 + n 4 X 4
is defined by

∆G r i = ∆G 0 r i + RT ln X n 3 3 X n 4 4 X n 1 1 X n 2 2
where ∆G 0

r i = n 3 ∆G 0 f X 3 + n 4 ∆G 0 f X 4 -n 1 ∆G 0 f X 1 -n 2 ∆G 0 f X 2
is a constant, R the gas constant, and T the absolute temperature in Kelvin.

Generalisation: With m metabolites et n i internal reactions, let us set

∆G r (X) =   ∆G r 1 (X) ... ∆G r n i (X)   Y (X) =    ∆G 0 f X 1 + RT ln(X 1 ) ... ∆G 0 f X m + RT ln(X m )    Proposition. ∆G r (X) = t SY (X)
Simple Example: problem n o 1

S =        1 0 0 0 -1 -1 -1 0 0 0 0 1 0 0 1 0 0 -1 -1 0 0 0 0 0 0 1 0 1 0 -1 0 0 0 0 0 0 1 0 0 -1 0 0 0 -1 0 0 0 0 0 1 0 0 -1 0 0 0 0 0 1 1       
if v l,j,min = 0 ∀j = 1, 2, 8 and v l,j,min = v l,j,max = 1 for j = 1, 2 then the unique solution is V ≈ (1, 1, 1.31, 0.69, 0, 0.31, 0.69, 0.38, 0.61, 0.69)

S =        -1 -1 -1 0 0 0 1 0 0 -1 -1 0 0 1 0 1 0 -1 0 0 1 0 0 -1 0 0 0 0 0 1 0 0 0 0 1 1       

New Approach

In this new modelling, in addition to the fluxes, we introduce metabolite concentrations (as in [START_REF] Hoppe | Including metabolite concentrations into flux balance analysis: thermodynamic realizability as a constraint on flux distributions in metabolic networks[END_REF]), which are also supposed to be at steady state. In order to ensure thermodynamic reliability, we add the Gibbs free energy change of each reaction. Still assuming some "minimal effort", we simultaneously minimize ||V || 2 and ||∆G r (X)|| 2 . We add an other set of constraints B l for the metabolites to take into account some experimental measurements.

Mathematical Formulation

min (V,X)∈K l ||V || 2 , ||∆G r (X)|| 2 with K l = {(V, X) ∈ A l × B l , ∆G r i (X) v i < 0 for i = 1 • • • n i } and B l = {X ∈ R m + , M l X = T l }.
This problem is a multi-objective optimization with uncoupled functions but with constraints (respect the second law of thermodynamics ∆G r i (X) v i < 0) that may be conflicting. So a rigorous answer to this problem comes from understanding the separated problems.

Study of problem n o 2

Proposition. min X∈R m + ||∆G r (X)|| 2 admits a space of solutions of dimension equals to dim(ker( t S)).

The problem without constraints is solved...

Simple example: problem

n o 2 dim(ker( t S)) = 1, so X = (X 1 , • • • , X 6 ) with X i = exp ( α -∆G 0 f X i

RT

), α ∈ R are the only solutions.

Conclusion

Why a mathematical study ? With only numerical computations, the set of optimal solutions may not be exhaustively described. To find all the solutions (if they exist) of the multiobjective problem, understand each separated problem is a pre-requisite. Moreover this formal outcome will fit to any other network.