Nathalie Vanfasse

Writing Letters to Come to 'Terms' with Domestic Economy: Household Management in Dickens's Early Correspondence"

In this essay, I will show how issues related to domestic economy are spelled out in epistolary form in Dickens's correspondence. This medium plunges us in an intimate form of perception of and reaction to these economic matters. A sample of Dickens's early correspondence exemplifies how his letters, at times, dealt with delicate financial issues involving close members of his family. The following epistolary exchanges yield invaluable insights first into the economic misfortune brought about repeatedly by Dickens's father upon his family, and then into Dickens's prospective domestic economy with his future wife, Catherine Hogarth. In Dickens's letters, family business dovetails with the negotiation of publishing contracts, 1 as well as with diverse other financial transactions. These business cum private matters partake of Dickens's epistolary biography. Besides, the very form and content 2 of these letters also produce a form of economic knowledge in that they show us how Dickens comes to "terms" literally and literarily with domestic economy. What unfolds before our eyes, as we read through these letters, is household management in progress. Thus, to some extent, Dickens's letters provide a form of economic and business knowledge.

In the light of recent analyses of epistolary discourse (e.g. Banks 2013), I will consider how Dickens tackles the financial matters already mentioned through epistolary modes and forms, in other words how these monetary considerations are expressed in his letters. In this respect, epistolary theory (e.g. Matthews-Schlinzig and Socha 2017) sheds light on Dickens's writing about domestic economy. On the basis of approaches to the discourse of the absent implied by letter-writing (e.g. Arrou-Vignod 1993), and with the help of the notion of an epistolary self (e.g. Earle 1999), as well as an analysis of the intimacy of letters between an individual and his nearest and dearest (e.g. Panzera and Canonica 2017), I will show that Dickens's epistolary self may be delineated through economic narratives of loss and redemption as well as through his prospects of household management.

Recent studies on business correspondence and private correspondence,

Dickens's letters, as do inquiries into notions of time, distance, as well as into the toing and froing of letters and their conversational style (e.g. Haroche-Boyzinac 1995).

Besides, an examination of the very medium, the page outlay, the writing, paper, mobility constraints, messengers, postal networks and other physical aspects of correspondence enable us to look into material aspects of Dickens's correspondence in connexion with the economic issues at stake. I will first use these perspectives to analyse his epistolary account of a domestic tragedy of insolvency.

A domestic tragedy of insolvency in epistolary form

Dickens's early correspondence reveals much about family business. One episode in particular stretches over several letters, as Dickens struggles to deal with his father's ineradicable propensity to run into debt, and to plunge the whole family into difficultya proclivity depicted in all the major biographies of Charles Dickens. This is a paradoxical situation to be in for a budding writer, like Dickens, whose letter reveals that though he recognizes and can predict recurring features in the events that are unfolding, he is by no means in control of the global storyline. The main protagonist of this unfortunate scenario, namely his father, is both the prime mover or main actor of this plot, and its victim, in that in running into debt, he recklessly acts against his better judgment, precipitating himself and his family from financial well-being into pecuniary distress. The rest of the family, barring Dickens, are endowed with very limited agency to face this calamity. They resemble spectators at a play watching with trepidation as events unfold. However, contrary to mere theatregoers, they are also directly affected by and feature as collateral damage of these events. They partake of this episode, though they are totally powerless. Literary categories, namely the theatre, and economic theory, namely limited agency, here dovetail to yield a better understanding of the family's financial predicament.

The financial episode I am interested in begins around 20 th November 1834.

Dickens is writing to Thomas Mitton, a close friend and solicitor, 3 to inform him of his father's disappearance. This is a case where Dickens seems to have resorted to letterwriting only as a second choice. His letter clearly shows that, in this matter, he would have preferred to confer directly with Mitton rather than resort to writing. Indeed, he informs his friend that he has tried to call upon him and that, short of having managed to see him on the previous day at his office at New Inn, he is now inviting him to dinner, in the hope of being able to discuss this delicate matter with him viva voce.

The pressing nature of this invitation to dinner -in which Dickens mentions that the family is even prepared to delay the time of their supper to accommodate their guest -confirms how keen Dickens was to confer with Thomas Mitton directly. The sensitive nature of the financial issue to be dealt with accounts for Dickens's obvious reluctance to put the matter into writing. Here, epistolary form is only second best to direct oral communication, in part because of the delicate nature of the revelations, and in part because Dickens visibly deems the issue easier to deal with through direct dialogue.

Mitton is someone whom he not only trusts, but whose professional skills he hopes will help him find a way out of his present predicament. In terms of epistolarity in connection with economic issues, the question that arises here is which form of communication seems the most appropriate to express and represent the misfortunes affecting a private domestic economy.

Short of being able to communicate the matter to Mitton orally, Dickens thus falls back upon letter-writing to impart it. He informs his friend and solicitor that after an unsuccessful attempt at finding some money to pay outstanding debts, his father has vanished, and that their landlord, Robert Burr, is unwilling to wait any longer for unpaid rent:

They are all in great tribulation here; my father went out yesterday accompanied by Alfred, to endeavour to get some money, as Burr refused to wait beyond last Evening. He sent the boy home to say he had been unsuccessful, and has not made his appearance all night, or forwarded a message of any kind. (Letters 1: 44) This letter thus delineates a portrait of Dickens's father through his son's pen. John Dickens appears -though this is never said in so many words but implied by the behaviour outlined in the document -as a person who systematically shirks his responsibilities. Having failed to find the money needed to settle his accounts, he not only does not return home, but no longer gives his family any sign of life. Instead, he somewhat cowardly chooses to send one of his sons back home with the bad news of his insolvency, while he himself shows a clean pair of heels.

In the light of this letter, what primarily characterizes John Dickens is avoidance behaviour. This avoidance is qualified by Dickens as a typical stance adopted by his father in similar situations. The letter hints at a repetitive pattern, or well-known scenario that Dickens and his family have unfortunately come to be used to. As a result, Dickens confesses that he feels no immediate alarm regarding his father's disappearance, since this unfortunate event is a sadly recurring one. Repetition has led the family to anticipate a well-known sequence of events rather than experience surprise:

I own that at present his absence does not give me any very great uneasiness, knowing how apt he is to get out of the way when anything goes wrong. (Letters 1: 44)

The repeated events have created a corresponding horizon of expectation. Economic anticipations and literary reception here meet in capturing what is happening to the family. However, the situation as a whole, namely that the chief of the family has absconded and left his family at the mercy of their creditors, is by no means less stressful for all that, since it plunges the family into economic uncertainty.

Dickens himself does not lead the game here as he does in his sketches. He is at the mercy of unexpected developments that he cannot anticipate, and he finds himself reacting to circumstances rather than devising them. He knows neither whether nor when his father will show up again. Economic uncertainty here plays out as a suspenseful narrative, and the letters show how Dickens reacts to this indeterminacy:

As, however, if we do not gain some tidings of him in the course of the day I shall hardly know what to think, and as at all events it may be necessary to hold a consultation, you will greatly oblige me by complying with my request. (Letters 1: 44) All Dickens can do is play it by ear. In his letter, he depicts himself faced with the suspense and uncertainties of real life that is playing out a private drama of insolvency.

He actually uses this theatrical simile later in his letter, when he describes this financial plight in literary and theatrical terms as a "domestic tragedy": I have not yet been taken, but no doubt that will be the next act in this "domestic tragedy". (Letters 1: 45) This is a scenario in the making over which Dickens has no control, though he does anticipate and dread some of its developments. In the unfolding of the affair, Dickens foresees that he is to be a passive victim rather than an actor of his own destiny. His letters show him desperately struggling to regain control and agency over a course of events that eludes him. Dickens strives to recover financial control over his family's situation. And in order to do so, he endeavours to frame their predicament through letter-writing, in a desperate attempt to endow events with some sense of purpose or teleology akin to storytelling or plot elaboration.

Thus in the next letter to Thomas Mitton, dated 21 st November 1834, Dickens informs his friend that his father has just been arrested and incarcerated in a sponging house in Cursitor Street. The creditors this time are wine merchants. Dickens finds himself struggling between his own professional commitments and the need to find time to sort out his father's financial embarrassments. He even fears that he himself might be arrested, presumably, as indicated in the notes to the Pilgrim edition of his letters, because he backed some of his father's bills.

The tone of the letter is factual and non-accusatory, but anxiety transpires through the increasingly pressing, though polite, exhortations made by Dickens to Tom Mitton to assist him in this crisis. Dickens reiterates his invitation to have Mitton for dinner, no doubt in the hope of discussing the matter with him directly rather than in writing, and to seek advice: Dear Tom. On waking this morning, I was informed that my father (whom I have not seen, for he had gone out) had just been arrested by Shaw and Maxwell the quondam Wine People. Will you have the goodness the moment you receive this, to go over to Sloman's in Cursitor Street and see him for me, and ascertain whether anything can be done? I must work this morning, but tell him please I shall be with him about Six. -I dine at 5 if you can come up. (Letters 1: 45) Added to the anxiety related to his father's incarceration is Dickens's ignorance of the extent of his father's insolvency. His uncertainty and extreme apprehension as to the magnitude of the sums owed to diverse creditors is emphasized by the use of italics in the very body of the letter. This epistolary snippet again exemplifies the uncertainties agents can be confronted with in the real economy, and how they try to deal with this lack of economic knowledge: I fear it is an awkward business, and really I have no idea of the extent of his other embarrassments. (Letters 1: 45) What the letters have shown us so far are the intense deliberations mixed with emotions that this domestic tragedy of insolvency has given rise to. One can recognize here what John Forster in his biography called Dickens's "passionate resolve, even while he was yielding to circumstances, not to be what circumstances were conspiring to make him" (38). Crushing financial difficulties, which reveal Dickens's very limited economic agency, are fruitfully turned into deliberation followed by action through the process of letter-writing. Reflecting upon his predicament enables Dickens to outline a plan of action and its modalities. Thus letter-writing is part and parcel of Dickens's way of dealing with financial setbacks. In an effort to overcome adversity, he applies himself to repaying his father's debts. A subsequent letter dated 21 st November 1834 to Thomas Mitton shows him in the process of doing this. However the settlement described is more convoluted than one would expect. John Dickens's insolvency thus draws Dickens into fairly complicated financial schemes involving sundry intermediaries. 4 See the rest of this article in : Letter-Writer (1768), both of which profferred model letters for every occasion, he did not believe that letter-writing could be reduced to easily learned formulas" (36).

3 David Paroissien mentions letters to Thomas Mitton retained as archival materials pertaining to matters of business ( 35). Here one can see that the issue at hand is related both to Dickens's business and domestic life. 4 This kind of communication combining the standards of business letter and those of more personal epistles are of the kind analyzed by Julien Harang in his book L'épistolaire (2002).

• « Writing Letters to Come to 'Terms' with Domestic Economy : Household Management in Dickens' Early Correspondence », Nineteenth-Century Prose, vol. 46, n°1, 2019, 123-148. Notes

  1 This essay builds upon information drawn from Robert Patten's book Charles Dickens and his Publishers (2018). 2 Regarding epistolary form however, David Paroissien points out that "Although Dickens owned the 18-volume set of Elegant Extracts in Prose, Verse and Epistles (1812), and was familiar with The Complete


Sedat Tardu

ANALYSIS OF THE ELECTRIC DOUBLE LAYER EFFECT ON MICROCHANNEL FLOW STABILITY

Keywords: microfluidics, microchannels, EDL, linear stability, early transition

INTRODUCTION

The "macro effects" should not be confused with real "micro effects" when explaining the deviations from the classical Navier-Stokes approach observed in microchannel flows. The roughness, for instance, is not a micro effect, since on the same scale it produces similar behavior in macro flows. Microscale effects can be neglected in conventional situations, but become increasingly important as the characteristic length of the flow decreases. Molecular, interfacial, and surface effects enter into this category.

Contrary to gas flows, the molecular effects in micro liquid flows are difficult to model. The slip length and the temperature jump can be related to the Knudsen number through the momentum and temperature accommodation coefficients and first-or secondorder models in the rarefied regime of gas flows. Discontinuous boundary conditions in liquid flows and significant molecular effects hold only in extremely small devices having extremely high speeds [START_REF] Gad-El-Hak | The Fluid Mechanics of Microdevices[END_REF]. Tardu [START_REF] Tardu | Transferts Thermiques Dans Les Microcanaux[END_REF] used the universal relation of slip length as a function of the shear rate given by Thompson and Troian [3] to estimate the critical width of microchannels to have slip liquid flows. He concluded that the hydraulic diameters of such devices have to be smaller than 0.52 μm. Even then, the slip velocity is only 0.04% of the centerline velocity at the transitional Reynolds number. The interfacial effects are presumably most responsible for the deviations observed from the classical macro theory. One of the micro effects that may play an important role is the electric double layer (EDL). Most solid surfaces have an electric surface potential. The most common mechanism for the charging of surface layers in microfluidics is the depronotonation of surface groups including silica, glass, acrylic, and polyester [START_REF] Sharp | Liquid Flows in Microchannels[END_REF]. The electrostatic charges present on the solid surface attract the counterions to establish an electrical field. The ions are immobile in the compact layer next to the wall and less than 1 nm thick. They can, however, be transported by an imposed flow field in the diffuse EDL layer. Consequently, the counterion concentration near the wall is larger than in the bulk of the fluid. The imposed pressure gradient accumulates the mobile ions downstream and sets up an electrical field whose potential is defined as the streaming or electrokinetic potential. The streaming potential and the net charge density induces a streamwise external force. In the steady state, the streaming current is in equilibrium with the conduction current in the opposite direction. This allows the determination of the streaming potential and the velocity profiles under the EDL effect. In macroscale flows, these effects are negligible because the thickness of the EDL is minor compared to the height of the channel. However, in micro flows the EDL plays a rather significant role.

NOMENCLATURE

Well-controlled recent experiments have clearly confirmed that the EDL can explain the behavior of the Poiseuille number in laminar regime, providing that the liquid contains a very small amount of ions. Kulinsky et al. [START_REF] Kulinsky | Electroviscous Effects in Microchannels[END_REF], for example, reported an increase of 70% of the friction factor under the EDL effect in planar channels of 4 μm heights.

Early transition to turbulence in microchannel flows has been reported by several investigations in the past. Wang and Peng [START_REF] Wang | Experimental Investigation on Liquid Forced Convection Heat Transfer through Microchannels[END_REF] indicated that the transition occurred at Re < 300 in microchannels 0.2-0.8 mm wide and 0.7 mm deep, which is significantly smaller than 750 corresponding to macroscale flows. (The Reynolds number throughout this article is based on the centerline velocity and half height of the channel.) A transitional number of about 500 has also been reported in the experiments conducted in 1-mm-wide trapezoidal microchannels with channel depths ranging between 79 and 325 μm [START_REF] Rahman | Experimental Measurements of Fluid Flow and Heat Transfer in Microchannel Cooling Passages in a Chip Substrate[END_REF][START_REF] Rahman | Design, Fabrication, and Testing Microchannel Heat Sinks for Aircraft Avionics Cooling[END_REF][START_REF] Gui | Enhanced Heat Transfer in the Entrance Region of Microchannels[END_REF]. The micro tube experiments of Mala and Li [START_REF] Mala | Flow Characteristics of Water in Microtubes[END_REF] indicate that there is an early transition from laminar to turbulent flows for Re > 112-338. Pointing at a plausible EDL effect, they have clearly stated that "There is an early transition from laminar to turbulent flow, and the range of Re t values varies somewhat, depending on the diameter and the material of the wall." These effects may also be attributed to the roughness that may both influence the transition in the entry region, or directly affect the flow behavior through an implied roughness viscosity [START_REF] Qu | Pressure-Driven Water Flows in Trapezoidal Silicon Microchannels[END_REF]. It may, however, also be presumed that direct interfacial effects play some role in early transition.

The main goal of this study is to directly investigate the EDL effect on the linear stability of planar channel flow and estimate indirectly the resulting transitional Reynolds numbers. It is asked here whether electrokinetic effects on liquid flows may cause the early transition or not and whether the linear stability EDL flow characteristics may explain the small transitional Reynolds numbers reported in some recent microchannel experiments.

EQUATIONS GOVERNING 2-D CHANNEL FLOW UNDER THE EDL EFFECT

The ion distribution in a symmetric plane channel flow is given by the equilibrium Boltzmann equation:

n ± = n 0 exp ± zeψ k b T (1)
where n 0 is the bulk ionic concentration, ψ is the electrical potential, z is the valence of positive or negative ions, e is the electron charge, k b is the Boltzmann constant, and T is the absolute temperature. The effects of finite ion size, the gradients of the dielectric strength and viscosity are neglected and the fluid properties (viscosity and permittivity) are constant. The net charge density in a unit volume of fluid is:

ρ = (n + -n -)ze = -2n 0 sinh zeψ k b T (2)
The Poisson equation relating the electrical potential and the net charge density per unit volume reads as

∇ 2 ψ = - ρ ε 0 ε (3)
with ε and ε 0 standing, respectively, for the dielectric constant of the medium and the permittivity of vacuum. Equation (3) combined with Eq. ( 2) may be written in nondimensional form as

∇ 2 ψ = -κ 2 ρ(y) = κ 2 sinh(ψ) (4)
The scaling reference potential and charge density are, respectively, E S = k b T ze and ρ S = n 0 ze. Equation (4) makes the nondimensional Debye parameter κ = ka appear with:

k = (2n 0 z 2 e 2 /εε 0 k b T ) 1/2 (5)
whose inverse is the EDL thickness, and a is the channel half height. Equation ( 4) may be linearized when the electrical potential is small compared with the thermal energy of ions, i.e., when |zeψ| |k b T |, or ψ 1 to give:

∇ 2 ψ = κ 2 ψ (6)
Equation ( 6) holds remarkably well for ψ < 2. One has:

d 2 ψ dy 2 = κ 2 ψ (7)
for a two-dimensional channel, where y = Y/a stands for the nondimensional wall normal distance.

The resulting charge density is ρ(y) = -ψ = -1

κ 2 d 2 ψ
dy 2 . The boundary conditions for Eq. ( 7) are:

• At the wall ψ = ψ 0 . The electrical potential at the wall can be approximated by the zeta potential ψ 0 = ζ . The general solution of the resulting EDL-velocity profile is not affected by this approximation. However, it is generally difficult to measure ψ 0 . It is easy to determine experimentally the zeta potential at the boundary between the compact and diffuse layers by determining the streaming potential and the pressure drop [START_REF] Lykelma | On Surface Conduction and Its Role in Electrokinetics[END_REF]. • At the centerline ψ = 0, providing that the EDL diffuse layers developing on both sides of the channel do not overlap. The solution of Eq. ( 7) is:

ψ = ζ sinh(κ) | sinh(κy)| (8)
The charge density ρ(y) and the streaming electrical field strength caused by the ions mobilized under the imposed pressure gradient acts as an external force in the streamwise momentum equation of a fully developed 2-D channel flow:

0 = μ d 2 U dY 2 - dP dX + E x ρ(Y ) (9)
The electrical field strength is related to the streaming potential by E x = E s l , where l is the length of the channel. The nondimensional form of Eq. ( 9) is:

0 = d 2 u dY 2 + 2 - 2 E s κ 2 d 2 ψ dy 2 (10)
where = ζ n 0 zea 2 lμU s and the scaling velocity is the centerline velocity of the Poiseuille component, i.e., U s = -a 2 dp/dx 2μ . Integrating Eq. ( 10) with the appropriate boundary conditions gives the velocity profile:

u = 1 -y 2 -4 E s ζ κ 2 1 - sinh κy sinh κ (11)
The last equation is the nondimensional form of the velocity profile already obtained by Mala et al. [START_REF] Mala | Heat Transfer and Fluid Flow in Microchannels[END_REF]. The streaming current due to the flux of the net density charge is given by the following cross integral by definition:

I s = Uρ(Y )dY dZ ( 12
)
The streaming current produces the conduction current in the reverse direction:

I c = E x λ 0 S = E s l λ 0 S ( 13
)
where λ 0 is the electrical conductivity of the fluid, which is proportional to the ionic concentration. In the steady state there is no net current I s + I c = 0. This allows the closure of the equations and the determination of the streaming potential. The final version of the nondimensional velocity profile is:

u = 1 -y 2 -4 I 1 -I 2 κ 2 sinh κ ζ 2 G + 4 I 3 - I 4 sinh κ 1 - sinh κy sinh κ (14)
The parameter G is given by G = (n 0 zea) 2 λ 0 μ . The quantities I in Eq. ( 1) are:

I 1 = I 3 = cosh κ -1 κ I 2 = 1 κ + 2 κ 3 cosh κ - 2 κ 2 sinh κ - 2 κ 3 I 4 = sinh κ cosh κ 2κ - 1 2 (15) 1
1 There is a small printing error in Mala et al. [START_REF] Mala | Heat Transfer and Fluid Flow in Microchannels[END_REF], p. 3082

It is seen that the velocity profile can be decomposed into a macro Poiseuille component plus an EDL effect component. A closer inspection of Eq. ( 14) allows expressing the latter in another closed form as:

u EDL = -4 I 1 -I 2 λ 0 μ sinh κ ζ 2 + 4 I 3 - I 4 sinh κ 1 - sinh κy sinh κ (16)
The parameter λ 0 μ is the nondimensional electrical conductivity times the dynamic viscosity of the fluid. It can be seen that the scaling parameter for λ 0 μ is the square of the scaling charge density times the scaling length, i.e., (ρ S /k) 2 leading to:

λ 0 μ = λ 0 μ (n 0 εε 0 k b T /2) ( 17
)
There are three parameters governing the EDL effect, namely

u EDL = u EDL (κ, ζ , λ 0 μ) (18)
Due to the grouping of the parameters this relation is equivalent to:

u EDL = u EDL (κ, Gζ 2 ) ≡ u EDL κ, λ 0 μ ζ 2 (19)
as clearly suggested by Eqs. ( 14) and ( 16). Gradients of the dielectric strength, viscosity, and conductivity should be incorporated in more realistic models. Yet these effects are presumably negligible in significantly dilute solutions, such as ultra-deionized, ultra-filtered water and pure organic liquids.

The main EDL effects may be summarized as follows:

• An increase of the friction constant and apparent viscosity.

• A decrease of the Nusselt number.

These effects are persistent, yet not significantly important, at least for large values of κ.

The apparent viscosity increases, for example, by a factor of nearly 3 at κ = 2, but the EDL effect on the wall shear stress disappears quickly when κ ≥ 10 [START_REF] Mala | Heat Transfer and Fluid Flow in Microchannels[END_REF]. The Nusselt number decreases by nearly 40% at κ = 5, and less than 5% at κ = 50 [START_REF] Tardu | Transferts Thermiques Dans Les Microcanaux[END_REF].

INFLEXIONAL INVISCID INSTABILITY

One important aspect of the EDL velocity profile is the existence of an inflexional, or, in other words, maximum vorticity point. One has indeed:

∂ 2 u ∂y 2 = -2 + 4 κ 2 (I 1 -I 2 ) sinh κ λ 0 μ sinh κ ζ 2 + 4 I 3 - I 4 sinh κ sinh κy ( 20
)
which can also be written as:

u = ∂ 2 u ∂y 2 = -2 - Rκ 2 sinh κ sinh κy (21)
where R ≤ 0 is the ratio of the EDL to the Poiseuille velocity components at the centerline of the channel, i.e., R = u EDL u Pois at y = 0. The maximum vorticity takes place at the inflexion point:

y I + ± 1 κ arc sinh -2 sinh κ Rκ 2 ≈ ± 1 + 1 κ ln - 2 Rκ 2 (22)
The right-hand side of the last relationship is a good approximation to y I for κ ≥ 10.

Clearly, the EDL flow is inviscidly unstable. The inviscid instability is of Fjortoft type because after the inflexion point near the wall, u > 0 and the velocity is smaller than the velocity at the inflexion point u I , resulting in u (uu I ) < 0 in this zone. The present investigation is a parametrical one based on nondimensional quantities. Yet we will introduce a numerical case to fix the ideas. Thus, we choose a "reference EDL flow" with κ = 41, G = 12,720, and ζ = 2.1254. This case corresponds to the flow of an infinitely diluted KCl aqueous solution (n 0 = 3.764 × 10 19 m -3 ) through a microchannel of height 100 μm subject to a zeta potential of 50 mV. The Debye length is 1.2 μm, and the bulk conductivity is λ 0 = 8 × 10 -9 S/m. The parameters G and ζ and are identical to those of Mala et al. [START_REF] Mala | Flow Characteristics of Water in Microtubes[END_REF] wherein detailed electro-viscous effects can be found. More severe cases, leading to large EDL effects have been treated in the literature [START_REF] Yang | Modeling Forced Liquid Convection in Rectangular Microchannels with Electrokinetic Effects[END_REF].

Figure 1 shows the distribution of the ratio of the EDL to Poiseuille velocity components at the centerline, together with the location of the inflexion point, versus κ. It appears that R increases when κ decreases, showing that the decelerating EDL effect becomes important with either a small separation distance between the plates or a large EDL thickness. At the same time the inflexion point moves away from the wall toward a region wherein viscous effects are less significant. Thus, it may be expected that the microflow becomes increasingly susceptible to inviscid instability with decreasing κ. These arguments are logical, but a definite picture of the EDL linear stability mechanism can only be given by a detailed Orr-Sommerfeld analysis.

STABILITY ANALYSIS

The linear hydrodynamic stability under the EDL effect is studied through classical methods. The Orr-Sommerfeld equation is solved by a Galerkin-like procedure [START_REF] Von Kerczek | The Instability of Oscillatory Poiseuille Flow[END_REF]. The normal mode solutions of the disturbance equation are:

(u, v, p) = R( v, v, p) exp(iαx) ( 23
)
where R is the real part, α is the dimensionless wave number of the disturbance, and x is the streamwise coordinate. Introducing the stream function:

(x, y, t) = φ(y, t) exp(iαx) (24)
the Orr-Sommerfeld equation takes the form:

∂ ∂t Lφ = 1 Re L 2 φ -iα uLφ - d 2 u dy 2 φ ( 25
)
with the boundary conditions φ = ∂φ ∂y = 0 at y = ±1. The operator L is L ≡ ∂ 2 ∂y 2α 2 . The stream function is expanded in a Chebyshev polynomial series:

φ(y, t) = N n=1 a n (t)T 2n-2 (y) ( 26
)
where T m (y) = cos(m cos -1 y) denotes the Chebyshev polynomials of the first kind. We took N = 256 through this study. Equation ( 21) takes the form:

Q da dt = (P -iαJ) • a (27)
and the matrices are determined by using the τ -method described by Orszag [START_REF] Orszag | Accurate Solution of Orr-Somerfeld Stability Equation[END_REF]. The last equation can be written as:

db dt = D • b (28)
by introducing B the matrix that diagonalizes

Q -1 • (P -iαJ), b = B -1
• a, and the diagonal matrix

D = λ i δ ij = B -1 • Q -1 • (P -iαJ) • B.
The eigenvalues are denoted by λ i and the flow is unstable when. The disturbances with symmetric streamfunctions are considered only. The method had similar results when compared to Grosch and Salwen [START_REF] Grosch | The Stability of Steady and Time-Dependent Plane Poiseuille Flow[END_REF], who used different sets of expansion functions.

RESULTS

Effect of the Debye-Hückel Parameter

The EDL flow depends on the groups of the parameters κ and Gζ 2 (Eq. ( 19)).

When the effect of the nondimensional Debye-Hückel length is first analyzed, the parameter Gζ 2 is kept constant for a given κ. The sensitivity of the hydrodynamic stability characteristics to the zeta potential and the parameter G will be analyzed in the next section. Recall that Gζ 2 = 57,460 and κ = 41 in the reference case introduced previously.

For subsequent computations, the zeta potential, the ionic concentration, and the physical characteristics of the fluid (the permittivity, conductivity and viscosity) are kept constant and the half channel height a is varied alone. Therefore, κ and G decrease proportionally with a and a 2 , respectively. Thus, Gζ 2 is 8740 and 2185 at, for example, κ = 16 and κ = 8, respectively.

The neutral curves deduced from the hydrodynamic stability analysis are summarized in Figure 2. The friction factor increases by only 10% at κ = 41, but the critical Reynolds number decreases by a factor nearly equal to 2: the critical wave and Reynolds numbers of the micro flow are, respectively, α c = 1.10 and Re c = 3190, to be compared with α c = 1.02 and Re c = 5772 of the conventional Poiseuille flow. The destabilization under EDL disappears quickly as κ increases. Thus, at κ = 164 (a minichannel of 400 μm height) the critical Reynolds number is Re c = 5730 (triangle in Fig. 2), which only slightly differs from Re c = 5772. The micro EDL effects for the micro/minichannels of height larger than typically 100 μm are not significant. Note, however, that, in a purely The distribution of Re c and α c as a function of κ is shown in Figure 3. The critical Reynolds number changes linearly versus log(κ) in a large range of the Debye-Hückel parameter before reaching a plateau region near κ = 4 (Re c = 279, α c = 1.9). The critical wave number, on the other hand, increases linearly with log(κ) when κ < 41. The analysis is deliberately stopped at κ = 4, since there may be some adverse effects due to the overlap of the EDL diffuse layers developing on both sides of the channel for κ < 4.

There are two main features of the neutral curves shown in Figure 2. First, the band of unstable wave numbers of the EDL micro flow is significantly larger compared with the Poiseuille macro flow due to the inviscid inflexional instability (unstable for Re → ∞ for a given α). The smaller the parameter κ, the broader the neutral curve is. Second, the upper branch of the neutral curves diverges slowly as the Reynolds number increases, while the lower branch of the micro and macro flows have similar behaviors. This can be seen with a careful look at Figure 2 of the EDL, the upper and lower branches diverge so that the neutral curves are open. This means that the inviscid instability overcomes the shear layer instability at Re → ∞ as expected. When κ is low, the inflexion point leaves the viscosity-dominated near wall region, so that the inviscid instability mechanism can become more operative. The socalled divergence is slow. The wave number begins to increase at Re = 200 × 10 3 for κ = 41, and earlier when the EDL layer is thicker, i.e., at approximately Re = 10 × 10 3 for κ = 16, and κ = 8.

Spectrum of Eigenvalues

Figure 5 shows the spectrum of eigenvalues as a function of Re for a single fixed wave number α = 1 and κ = 16.4. Only the first three eigenvalues corresponding to the complex wave speed c = -λ/iα are shown for the sake of clarity. The modes are ordered in such a way that ascending values correspond to descending values of Im{c}. Remember that the flow is stable (unstable) for R{λ} < 0 (> 0), or equivalently when Im{c} < 0 (> 0). The Re was varied from 300 to 20 × 10 3 . The curves corresponding to Poiseuille and EDL flows are denoted respectively by Pi and Ei, where i is the mode number.

Figure 5b shows that only the first mode has a different Im{c} under the EDL effect. The imaginary parts of higher modes collapse entirely with those of the macroscale flow. Note that the imaginary part of the first mode in micro flow denoted by Im{c E-1 } is larger than that of the macroscale Poiseuille flow Im{c P -1 }, showing the strong destabilizing effect of the EDL.

The EDL modes are slower than the Poiseuille ones, in the sense that R{c E-i } is smaller than R{c P -i } for all modes (Fig. 5a). At moderate Reynolds number, they may be classified as slow modes since the real part of the complex wave speed is smaller than the average velocity of the channel, i.e., R{c E-i } < 2 3 for roughly Re < 5 × 10 3 . Modes 2 and 3 of macroscale flow are clearly fast modes, since their phase speed is close to the centerline velocity.

The eigenvalue spectrum is further analyzed for a more destabilizing case at κ = 8. Figure 6 shows the resulting distribution of the complex wave speed real and imaginary parts. It is seen in Figure 6a that all EDL modes are now slow modes in the range of Reynolds numbers investigated here. The differences between R{c E-i } and R{c P -i } are larger compared to the case κ = 16. The parameter κ does not appreciably affect the imaginary part of the spectrum as can be seen in Figure 6b, except a slight increase with decreasing κ in the first mode.

Stability Sensitivity to the Zeta Potential and the Parameter G

We have previously shown that the EDL velocity component depends on the group of parameter Gζ 2 for a given Debye-Hückel parameter κ (Eq. ( 14)). The parameter G can be expressed as a function of the nondimensional conductivity and viscosity product G = κ 2 λ 0 μ , showing that G varies like the inverse of λ 0 μ = λ 0 μ (n 0 εε 0 k b T /2) for fixed κ. Thus, G increases with decreasing conductivity and/or viscosity.

To each couple (κ, Gζ 2 ) corresponds a different set of the critical wave and Reynolds numbers (α c , Re c ). It is difficult to illustrate the interdependency of Re c via Gζ 2 , since α c changes continuously as Gζ 2 varies. Therefore, we fixed the wave number at its macroscale flow critical value α c = 1.02 to give a general tendency and extract clear information. [START_REF] Qu | Pressure-Driven Water Flows in Trapezoidal Silicon Microchannels[END_REF] Trapezoidal silicon microchannels 71-23.5 Gao et al. [START_REF] Gao | Scale Effects on Hydrodynamics and Heat Transfer in Two-Dimensional Mini and Microchannels[END_REF] 2-D microchannels 100-53 with κ (Table 1) and it is approximately two decades large. The lower and upper limits of the Gζ 2 range increase with κ. It is clear from both Figure 7 and Table 1 that Re c is more sensitive to Gζ 2 when the Debye-Hückel parameter is small enough, i.e., for κ ≤ 40. The first and main condition to expect significant EDL effects on hydrodynamic linear stability is to deal with liquids of low ionic concentration. The second condition requires large enough Gζ 2 value. This is fulfilled by low conductivity, which is generally associated with wall/liquids interactions that lead to high zeta potential ζ . The second condition is less critical than the first. To show this, consider a median value of κ = 20.

Figure 7 shows that Re c decreases by a factor of 6 at Gζ 2 = 10 3 , corresponding, for example, to G = 250 and ζ = 2; that is, a liquid that is 50 times more conductive than the reference case.

DISCUSSION AND CONCLUSION

The electro-viscous effects modify considerably the linear stability characteristics of microchannel flows. Two main conditions are required for a significant alteration of the critical Reynolds number:

1. The nondimensional Debye-Hückel parameter should be smaller than κ ≤ 40.

There is no effect on stability for liquids with high ionic concentration and/or minichannels. 2. Large enough zeta potential and liquids with low conductivity/viscosity leading typically to Gζ 2 ≥ 5000 for 1 k ≈ 1 μm are necessary. This condition is less critical. With a Debye length 1 k = 2 μm, which is just two times larger than the reference case, the critical Reynolds number decreases by nearly 80% in a 32 μm width channel at Gζ 2 as small as Gζ 2 = 100. The key question: Can this effect be experimentally verified and validate early transition reported in some recent investigations?

To give preliminary answers to these questions and discuss related phenomena, we first have to note that the hydrodynamic linear stability and transition mechanisms are different in plane Poiseuille flows. What we are interested in from an applications point of view is indeed the transitional Reynolds number Re t related but not strictly equal to Re c .

According to Orszag and Patera [START_REF] Orszag | Secondary Instability of Wall Bounded Shear Flows[END_REF], the instability process of wall bounded macroscale shear flows involves three steps:

1. Primary (linear) instability of the basic shear flow 2. Nonlinear saturation of the primary instability and formation of a secondary flow 3. Secondary instability, which is the linear stability of the secondary flow The first step gives Re c = 5772 for macroscale plane Poiseuille flows, but transitional Reynolds numbers as low as Re t = 400 have been reported in experiments. This is because the second step results in a Landau-type mechanism leading to subcritical instability or the metastability. The instability occurs with finite amplitude when all infinitesimal disturbances are stable at a significantly lower critical value compared to Re c [START_REF] Drazin | Hydrodynamic Stability[END_REF]. The quasi-steady nonlinear analysis conducted by Orszag and Patera [START_REF] Orszag | Secondary Instability of Wall Bounded Shear Flows[END_REF] suggests a critical Reynolds number for the secondary instability of R ∼ = 400. Three-dimensional disturbances decay for R < 1000; however, 1000 cannot be taken as a precise number because the three-dimensional instability requires that a threshold two-dimensional amplitude be achieved, which depends on the perturbation level. The experiments indicate that the ratio of the transitional Reynolds number to the critical Reynolds number of the primary instability is approximately r = Re t Re c ≈ 750 5772 ≈ 1 8 [START_REF] Schlichting | Boundary Layer Theory[END_REF] but this is not an exact universal value either. Some arguments suggest that the EDL can reinforce the subcritical nature of macro Poiseuille flow through the nonlinear saturation of the primary instability (step 2). Stuart [START_REF] Stuart | On the Non-Linear Mechanics of Wave Disturbances in Stable and Unstable Parallel Flows. Part 1. The Basic Behaviour in Plane Poiseuille Flow[END_REF] has shown that there are two main mechanisms generating the subcritical equilibrium state. The first mechanism prevents the subcritical instability through the distortion of the mean motion. This effect is proportional to Re c α 2 c , which is significantly small in EDL flow. The second mechanism reinforces the metastability by the wall normal distortion of the fundamental and it is inversely proportional to Re c . EDL directly enhances the subcritical state through the second mechanism and indirectly by retarding the first mechanism. One of the simple explanations of the secondary instability (step 3) is the creation of local inflexional streamwise profiles with strong shear susceptible to inviscid instability. The basic flow in EDL is already inflexional. Interactions of the first and second steps may lead to significantly smaller transitional Reynolds numbers in the presence of EDL. The reasonable conclusion of these arguments is that the EDL flow is subcritical and that presumably r = Re t Re c ≤ 1 8 . It is somewhat difficult to reach high Reynolds numbers in microchannels because high pressure can cause channels to break. We show in Figure 8 the estimated transi- 2.

tional Reynolds as a function of κ, together with some microchannel and microtube data reported in the literature. Recall that Hagen-Poiseuille flow is linearly stable, and some aspects of the transition in pipe flows can be explained only by secondary instabilities. We supposed that r = Re t Re c = 1 8 as in macroscale flow. We need a Debye length to express data in terms of Reκ distribution. We took 1 k = 1 μm to be consistent. The open circles in Figure 8 correspond to the reference case. A second case, shown by bold circles, is also introduced. It consists of a liquid that is 11.5 times less conductive than the reference case leading to Gζ 2 = 5000 at κ = 41. It is clearly seen in Figure 8 that the experimental Reynolds numbers are larger than Re t , concluding that experimental verification of the present study is plausible.

Early transition in microchannel flows has been reported in several experimental studies as we already indicated in the introduction. There are two data points related to transition in Figure 8. The first is from Mala and Li [START_REF] Mala | Flow Characteristics of Water in Microtubes[END_REF], who indicated Re t ≈ 300 in microtubes with hydraulic diameters D h = 130 μm. The second is from trapezoidal silicon microchannel experiments of Qu et al. [START_REF] Qu | Pressure-Driven Water Flows in Trapezoidal Silicon Microchannels[END_REF]Fig. 8b,p. 362] and is shown by a bold square in Figure 8. The results obtained by these authors at D h = 142 μm suggest a transitional number of 400. There is a satisfactory agreement between the estimated transitional numbers and the experiments in Figure 8, but that may be coincidental. Controlled experiments in a way similar to those reported by Ren et al. [START_REF] Ren | Electro-Viscous Effects on Liquid Flow in Microchannels[END_REF] have to be conducted. That may be achieved, for instance, by keeping the same channel with the same roughness distribution and by changing the ionic concentration of the liquid. Somewhat unfortunately, the maximum Debye length in Ren et al.'s [START_REF] Ren | Electro-Viscous Effects on Liquid Flow in Microchannels[END_REF] experiments is four times smaller than the reference case taken here. Ren et al. [START_REF] Ren | Interfacial Electrokinetic Effects on Liquid Flow in Microchannels[END_REF] and Sze et al. [START_REF] Sze | Zeta-Potential Measurement Using Smoluchowski Equation and the Slope of the Current-Time Relationship in Electroosmotic Flow[END_REF] have measured the zeta potential (ζ = 245 mV), Debye length ( 1 k = 0.305 μm), and the bulk conductivity (λ 0 = 1.053 × 10 -4 S/m) in microchannels of height varying from 14.1 μm to 40.5 μm. Note that the conductivity is 1.3×10 3 times larger than the reference case. The present estimation is only qualitatively valid in this case because the zeta potential is too high for the Poisson equation linearization. Furthermore, the important contribution of the corners to the EDL field [START_REF] Yang | Electrokinetic Effects on Pressure-Driven Liquid Flows in Rectangular Microchannels[END_REF][START_REF] Yang | Analysis of Electrokinetic Effects on the Liquid Flow in Rectangular Microchannels[END_REF] and to the local flow instability has to be considered to model Ren et al.'s [START_REF] Ren | Electro-Viscous Effects on Liquid Flow in Microchannels[END_REF] experiments. With these restrictions in mind, one has κ = 23.6 and Gζ 2 = 500 for the 14.1 μm height channel reported by these authors. Up to 20% higher flow resistance is found and the transitional Reynolds number decreases by some 63%, with Re t ≈ 270 in this particular situation. There is a significant effect on the transition, but the resulting Re t is presumably too high to be experimentally reached in such a small device. In Wang and Peng [START_REF] Wang | Experimental Investigation on Liquid Forced Convection Heat Transfer through Microchannels[END_REF] the value of κ is presumably larger than 120 and the low transitional Reynolds numbers of approximately 300 they report cannot be explained by the EDL effect. Strong effect of roughness in the inlet region is suspected in these experiments.

Early transition in microchannels can have interesting applications for enhancement of mixing together with heat and mass transfer processes. A simple estimation shows that the ratio of the Nusselt numbers in EDL and macroscale flow is roughly inversely proportional to the ratio of the transitional Reynolds numbers. A decrease of Re t by a factor f therefore implies an enhancement of the heat transfer by the same f in some giving working conditions.

Fully developed turbulent EDL microchannel flow reveals several questions. A transitional Reynolds number of Re t = 817.5 (κ = 16, G = 1937) would result in a half height channel in wall units (i.e., scaled by the shear velocity and viscosity) of a + ≈ 30. Recalling that the logarithmic equilibrium layer begins at 30 wall units in macroscale conventional turbulent channel flows, the structure of "micro turbulence" in such situations is open to question. Our team tries to clarify these points by using direct numerical simulations.

Figure 1 .

 1 Figure 1. Distribution of R, the ratio of the EDL to Poiseuille velocity components (open circles and ordinate at left) and of y I , the inflexion point position (bold circles and ordinate at right) vs. the nondimensional Debye-Hückel parameter κ. The zeta potential is ζ = 2.1254. The parameters κ and G change as the half width a of the channel is varied (see the text). One has G = 12,720 at κ = 41 (reference case).

Figure 2 .

 2 Figure 2. Neutral curves of the EDL flow compared with the Poiseuille flow. The open circles correspond to Poiseuille flow with κ = ∞ and/or ζ = 0. The results of Grosh and Salwen [17] are shown by small bold squares. Bold circles correspond to κ = 41, G = 12,720, and ζ = 2.1254. The rest of the results are obtained by changing the microchannel height and keeping constant the rest of the parameters. The triangle is obtained for κ = 164 and G = 203,520. The neutral curves for κ = 8, G = 484 (bold squares), κ = 16, G = 1937 (open squares), and κ = 41, G = 12,720 (bold circles) are compared with the macroscale flow (open circles).

Figure 3 .

 3 Figure 3. Distribution of the critical wave (open symbols, ordinate at left) and Reynolds numbers (closed symbols, ordinate at right) deduced from Figure 2 vs. κ.

  for κ = 16 and κ = 8. The increase of the critical wave number corresponding to the upper branch of the neutral curves is seen more clearly in Figure4. The neutral curve in the macro Poiseuille flow is closed for large Reynolds numbers, whereby the upper branch joins the lower one when Re is large enough. This is simply because the Poiseuille flow is inviscidly stable. Under the presence

Figure 4 .

 4 Figure 4. Increase of the wave number along the upper branch of the neutral curves for three different values of κ deduced from Figure 2.

Figure 5 .

 5 Figure 5. Real (a) and imaginary parts (b) of the first three eigenvalues vs. Re, α = 1, and κ = 16. The eigenvalue i in EDL and Poiseuille flow is denoted by Ei and Pi, respectively.

Figure 6 .

 6 Figure 6. Real (a) and imaginary parts (b) of the first three eigenvalues vs. Re, α = 1, and κ = 8. The eigenvalue i in EDL and Poiseuille flow is denoted by Ei and Pi, respectively.

Figure 7 Figure 7 .

 77 Figure 7. Distribution of the critical Reynolds number vs. Gζ 2 for α = 1.02.

Figure 8 .

 8 Figure 8. Estimated values of the transitional numbers under the EDL effect compared with some experiments and the range of Reynolds numbers that can be reached in experiments in microchannels. Experimental configurations are resumed in Table2.

Table 1 .

 1 Sensitivity of the critical Reynolds number to Gζ

		2 variations
		Range of Gζ	2 within
		Re c decreases sharply
	κ	Lower value-upper value
	41	10 3 -10 5
	20	10 2 -10 4
	16.4	10 2 -10 4
	8.2	10-10 3

Table 2 .

 2 Experimental conditions related to the results presented in Figure8

	Authors	Configuration	Hydraulic radius (μm)
	Mala and Li [10]	Microtubes	102.5-25
	Qu et al.


Interfacial Electrokinetic Effect on the Microchannel Flow Linear Stability

The electrostatic double layer (EDL) effect on the linear hydrodynamic stability of micro-channel flows is investigated. It is shown that the EDL destabilizes the Poiseuille flow considerably . The critical Reynolds number decr eases by a factor five when the nondimensional Debye-Huckel parameter is ar ound ten. Thus, the transition may be quite rapid for microchannels of a couple of microns heights in particular when the liquid contains a very small number of ions. The EDL effect disappears quickly for у150 corr esponding typically to channels of heights 400 m or larger . These r esults may explain why significantly low critical Reynolds numbers have been encountered in some experiments dealing with microchannel flows.

Introduction

Both the rapid development of MEMS, and the need for high rates of heat removal of micro electronic devices require a clear understanding of micro flows in channels of heights ranging from some 100 to a couple of m. It is not clear whether the classical macro-approaches are still valid or not in such small devices. Several experimental studies published so far are quite inconsistent, and no clear answer to these questions has been given yet. ͑Gad-el Hak ͓1͔, Tardu ͓2͔͒.

In liquid flows, the molecular effects inducing slip boundary conditions become important only for microchannels of hydraulic diameters smaller than 1 m ͓1,2͔. Thus, it is rather difficult to explain the large departures from conventional transport theory observed in some microchannel flows of larger heights ͓3,4͔. One of the micro-effects that may play an important role is the interfacial effect caused by the electric double layer ͑EDL͒ at the solid/ liquid interface. The electrostatic charges present on the solid surface attract the counterions to establish an electrical field. In the compact layer next to the wall and less than 1 nm thick the ions are immobile. In the diffuse EDL layer however, the ions are less affected by the electrical field and can move. The counterion concentration near the wall is larger than in the bulk of the fluid. That results in a net charge density in a unit volume resulting from the concentration difference between cations and anions, according to the Boltzmann equation. The electrostatic potential at any point near the surface, provided that it is small compared to the energy of ions, may be obtained by a linear approximation of the Poisson-Boltzmann equation. Its value at the wall can be related to the Zeta potential between the compact layer and diffuse layer, when the EDLs near the opposite walls do not overlap. The Zeta potential is a property of the solid-liquid pair and can be determined experimentally ͓5͔. The imposed pressure gradient accumulates the mobile ions downstream and sets up an electrical field whose potential is called the streaming or electrokinetic potential. The streaming potential and the net charge density induces a streamwise external force. In the steady state, the streaming current due to the transport of charges is in equilibrium with the conduction current in the opposite direction. That allows the determination of the streaming potential and of the velocity profiles under the EDL effect. In macro scale flows, these effects are negligible, as well as the thickness of the EDL is very small compared to the height of the channel.

The velocity profile under the electrokinetic EDL effect obtained by Mala et al. ͓6͔ in the Debye-Hu ¨ckel approximation, can be put in non-dimensional form as ͓2͔;

uϭ1Ϫy 2 Ϫ4 I 1 ϪI 2 2 sinh ¯2G ϩ4 ͩ I 3 Ϫ I 4 sinh ͪ ͭ 1Ϫ ͯ sinh y sinh ͯͮ (1)
where the scaling velocity is the centerline velocity of the Poiseuille component, i.e. Ϫ a 2 dp/dx/2 and the scaling length is the half channel height a. There are several parameters in this equation, for instance Gϭ (n 0 z e a) 2 / 0 , with n 0 standing for the ionic number concentration, z for the valence of positive or negative ions, e for the electron charge 0 the electric conductivity of the fluid, and for its dynamic viscosity. One of the most important quantities involving in ͑1͒ is the non dimensional Debye-Hu ¨ckle parameter ϭa kϭa(2 n 0 z 2 e 2 / 0 k b T) 1/2 with and 0 being respectively the dielectric constant of the medium and the permittivity of vacuum, k b the Boltzmann constant and T the absolute temperature. The characteristic EDL thickness is 1/k. The non-dimensional Zeta potential reads for ¯ϭ z e/k b T . The quantities I in ͑1͒ are given by:

I 1 ϭI 3 ϭ cosh Ϫ1 , I 2 ϭ ͩ 1 ϩ 2 3 ͪ cosh Ϫ 2 2 sinh Ϫ 2 3 , I 4 ϭ sinh cosh 2 Ϫ 1 2 (2)
Fig. 1a compares the velocity profile of an ''EDL flow'' with the conventional Poiseuille flow, when ϭ41, Gϭ12720 and ϭ2.1254. This case corresponds to the flow of an infinitely diluted aqueous KCl solution (n 0 ϭ3.764ϫ10 19 m Ϫ3 , 0 ϭ7.89 ϫ10 Ϫ9 S/m) through a microchannel of height 100 m subjected to a Zeta potential of 50 mV. The parameters G ¯and are exactly those used by Mala et al. ͓6͔ who analyzed the electrokinetic effects of the double layer. The Debye length is 1.2 m which is close to the value corresponding to the de-ionized ultrafiltered water experiments of Ren et al. ͓5͔ The first impression one has from Fig. 1a is the close similarity between the EDL and Poiseuille profiles with a decrease of the centreline velocity typical to the EDL flows. The important difference however is the presence of an inflexional point at yϷ 1/ arcsin h ͕Ϫ 2/r 2 sinh()͖ in the EDL profile where r is the ratio of the 1

Sedat Tardu

Laboratoire des Ecoulements Géophysiques et Industriels, LEGI B.P. 53 X, 38041 Grenoble, Cédex France e-mail: Sedat.Tardu@hmg.inpg.fr EDL and Poiseuille flows centerline velocities. This makes the flow inviscidly unstable, according to the Fjortoft's criteria.

The linear hydrodynamic stability under the EDL effect is studied through classical methods. The Orr-Sommerfeld equation is solved by a Galerkin-like procedure ͓7͔. The stream function is expanded in a Chebyshev polynomial series with up to 256 terms. The -method described by Orszag ͓8͔ is used. The disturbances with symmetric streamfunctions are considered only. The eigenvalues of the matrices were numerically computed by the QR algorithm. The method gave very close results to Grosch and Salwen ͓9͔ who used different sets of expansion functions.

Figure 2 shows the neutral curve corresponding to the microchannel-flow with the parameters given below, together with a ''macro'' Poiseuille flow ͑i.e. ¯ϭ0 or →ϱ). It is clearly seen that the critical Reynolds number decreases by a factor nearly equal to 2 under the effect of EDL: the critical wave and Reynolds numbers of the microflow are respectively ␣ c ϭ1.10 and Re c ϭ3190 to be compared with ␣ c ϭ1.02 and Re c ϭ5772 of the conventional Poiseuille flow. Due to its inviscid inflexional instability ͑unstable for Re→ϱ for a given ␣͒, the band of unstable wave numbers of the EDL-microflow is significantly larger compared with the Poiseuille-macroflow. The destabilizing EDL effect dis-appears quickly when the height of the channel is increased by a factor of 4 ͑400 m͒ ͑triangle in Fig. 2 corresponding to ϭ163). This goes in the same line as previous experimental results showing the lack of micro-effects for the microchannels of height larger than typically 100 m. For smaller values of , in return, the effect of the interfacial effects caused by EDL on the transition may be much more severe. For instance, the critical Reynolds number decreases up to Re c ϭ1042 for a channel with a 40 m separation distance subject to the same conditions ͑the square in Fig. 2 with ϭ16), and to a value as small as Re c ϭ496 when ϭ8. This effect, combined with the influence of the roughness that can be critical in microchannel flows, may be the explanation of significantly small transitional Reynolds numbers reported in some studies, and the decrease of Re c with the decrease of channel size ͓3,4͔.

The second parameter that plays an important role in the stability mechanism under the EDL effect is G ¯2 ͑Eq. 1͒. It is clear that EDL will play a negligible role as G ¯2→0, that is, when the Zeta potential becomes very small. Yet, the inflexional instability effect is so strong that the critical Reynolds number is affected even when the electroviscous effects are negligibly small. Figure 3 shows the critical Reynolds number variation versus G ¯2 for ϭ16 and ␣ϭ1.02. It is seen that there is a decrease of 30% in Re c already at G ¯2ϭ100. The EDL and macro-Poiseuille velocity profiles coincide almost perfectly at this low value of G ¯2 ͑not shown here͒. The second striking feature of Fig. 3 is the rapid decrease of the critical Reynolds number as G ¯2 increases, before reaching a plateau region at G ¯2Ϸ10 4 . The macroscale Poiseuille flow is metastable; i.e., the corresponding stability is subcritical. Nonlinear analysis shows that the instability may occur with finite amplitude when all infinitesimal disturbances are stable ͓10͔. There is a significantly lower critical value of the Reynolds number Re G compared to Re c , above which the flow is unstable and below which there is no bifurcation. The exact theoretical determination of Re G is still a matter of research. Experiments show that Re G /Re c Ϸ1000/5772 Ϸ1/6. The transitional Reynolds number depends on the shape and shape factors of the channels.

The Poiseuille flow is monotonically stable only for ReϽ100 to be compared with 5772. The lower limit of the transitional Reynolds number Re t of channel flows is about 400 according to experiments and the growth rate of three-dimensional disturbances of linear secondary instability ͓11͔ ͑1100 based on the hydraulic diameter and channel averaged velocity͒. Thus, the ratio Re t /Re c Ͻ 1 15 in macroscale Poiseuille flows. The main non-linear stability mechanism in macro-Poiseuille flow is due to a secondary instability with the development of mean spanwise profiles inflexion with strong shear ͓12,13͔. In the EDL flow there is already a streamwise inflexion point in the base flow. Thus, the Reynolds numbers Re G and Re t should be much lower under the EDL effect. A nonlinear analysis is necessary to check out this point but the first stage is then the determination of the marginal curve which is done here. The nonlinear saturation of the primary stability and formation of a secondary flow, together with the secondary instability processes have to be analyzed in EDL flow similarly to the Poiseuille macro-flow ͓11͔. Some arguments on the reinforcing effect of the EDL on the subcritical nature of the macro Poiseuille flow may however already be given. The square of the amplitude of a finite disturbance is given by:

d͉A 1 2 ͉ dt ϭ2␣c i ͉A 1 ͉ 2 ϩ͑k 1 ϩk 2 ϩk 3 ͉͒A 1 ͉ 4 (3)
according to Stuart ͓14͔. The flow reaches a subcritical equilibrium state when k 1 ϩk 2 ϩk 3 Ͼ0. The coefficient k 1 represents the distortion of the mean motion: it is related to the eigenfunctions of the linear stability problem, and it is negative. The coefficient k 2 is linked to the generation of the harmonic of the fundamental and is also likely negative. The wall normal distortion of the fundamental (k 3 ) must ''be positive and outweigh the combined negative effect of k 1 and k 2 to reach a subcritical state'' ͓14͔. Now, k 1 is proportional to Re c ␣ c 2 ͑Eq. 6.3 of ͓14͔͒. It has therefore a significantly smaller negative contribution to k 1 ϩk 2 ϩk 3 under the EDL effect. Furthermore, part of the terms involving in the coefficient k 3 is inversely proportional to Re c according to the Eq. 6.5 of ͓14͔ and the EDL presumably reinforces the positive character of k 3 in the subcritical state.

The experimental verification of the linear stability results is difficult because of the subcritical character of the stability mechanism, and the difficult control of the level of turbulence at the inlet. That also rises the question of the EDL effect on the stability in the developing region. It is also somewhat difficult to reach such high Reynolds numbers in microchannels, even though one may consider the possibility of Debye lengths larger than 1 m, thus higher channels compared with the numerical example given before. Well-controlled experiments are certainly more difficult in microchannels. Yet, it should be possible to detect the EDL effects on the transitional Reynolds number by classical Repressure gradient curves. Table 1 shows the expected transitional Reynolds numbers under the EDL effect as a function of . The ratio Re t /Re c ϭ 1/15 is taken same as in macro channels, although, it is more than likely that Re t /Re c is presumably smaller under the EDL effect according to the arguments given above. Table 1 also shows the Reynolds numbers achieved in some experiments. It is seen that the expected transitional Reynolds number is within the range of experimental possibilities. For larger channel of typical height 100 m or so, the critical Reynolds number should also be measured, at least in principle. The experi-ments are limited to small Reynolds numbers for microchannels of height smaller than 20 m, because very high pressure is then necessary and the channels will then break.

The last column of Table 1 shows transitional Reynolds numbers reported by some authors indicating early transition in microchannel liquid flows. All of these experiments have been conducted with DIUF-W but the channel wall material may be different from one experiment to the other. In lack of experimental details it is hard to give definitive conclusions. Yet, there is a satisfactory agreement between the estimated theoretical Re t and the measurements, in particular with those of Mala and Li ͓15͔. There is only one value reported by them on which we are sure, but unfortunately other details are not available. They, however, clearly indicated that ''there is an early transition from laminar to turbulent flow at Re t у56Ϫ168'' and that ''the range of Re t values varies somewhat, depending on the diameter and the material of the wall'' pointing at a plausible EDL effect. Note that these experiments have been conducted in micro tubes, and that the Poiseuille-Hagen flow is linearly stable. That does not matter for the qualitative comparisons made here, since the nonlinear character of the stability is similar in planar and cylindrical channels.

To conclude, the EDL destabilizes the linear modes of the Poiseuille channel flow and early transition in microchannels is plausible. This effect can be experimentally checked, provided that the liquid contains a very small amount of ions, the Zeta potential is high enough and the channel height is sufficiently small. In practice that would require the use of DIUF-W or organic liquids, and channel heights larger than 20 m ͑but smaller than typically 100 m͒ for feasibility problems ͑microchannel failure at high pres-sures͒. There is no effect on stability for liquids with high ionic concentration and/or mini channels. The nonlinear stability analysis of the EDL flow is necessary, although physical considerations indicate a much more rapid transition compared to macro flows. Direct Numerical Simulations can also be helpful to this end. Controlled experiments in a way similar to those reported by Ren et al. ͓5͔ have to be conducted, by keeping the same channel with the same roughness distribution, and changing the ionic concentration of the liquid. The fact that the channels cross-section shape, in particular the corners have important contribution to the EDL field have also to be considered ͑Yang et al., ͓16,17͔͒. These results may imply important industrial applications, in particular in the control and manipulation of the heat transfer process in

Fig. 1

 1 Fig. 1 EDL and Poiseuille velocity profiles "a…, and the inflexional instability of Fjortoft type "b…. The parameters of the EDL flow are given in the text. The broken line in "b… corresponds to Poiseuille flow with Äؕ and the circles to the Ä41 EDL flow. The inflexion point is shown by I. The inviscid instability is of Fjortoft type because after the inflexion point near the wall, uЉ Ì0 and the velocity is smaller than the velocity at the inflexion point u I , resulting in uЉ"uÀu I …Ë0 in this zone.

Fig. 2

 2 Fig. 2 Neutral curves of the EDL flow compared with the Poiseuille flow. The open circles correspond to Poiseuille flow with Äؕ. The results of Grosh and Salwen †9 ‡ are shown by small bold squares. Bold circles correspond to Ä41, G Ä12720 and ¯Ä2.1254. The rest of the results are obtained by changing the microchannel height and keeping constant the rest of the parameters.

Table 1

 1 Estimated values of the transitional numbers under the EDL effect compared with some experiments and the range of Re numbers that can be reached in microchannels. 1: Mala and Li †15 ‡, microtubes, the transitional Re number of 150 "second line, last column… correspond to D h Õ2Ä65 m. 2: Gao et al. †18 ‡, 2-D microchannels, 3: Qu et al. †19 ‡, trapezoidal silicon microchannels. Profiles at D h Õ2-71 m "their Fig. 8… suggest a transitional number of 200. 4: Wang and Peng †4 ‡ microchannels. This is a particularly low transitional Reynolds number. The exact value of D h Õ2 is 150 m. The value of is larger than 120 for other microchannels these authors investigated. Although their experiments have been conducted with distilled water "so that roughly 1ÕkÄ1m… the low transitional Reynolds numbers they report for higher channels cannot be explained by the EDL effect. Strong effect of roughness in the inlet region is suspected in these experiments.

microchannels. The complete parametrical study including the EDL effects on amplification and decay rates, and the eigenvalue spectrum is under investigation. The nonlinear stability analysis of EDL flow is also currently investigated.


Éphane St

Brull

AN ELLIPSOIDAL STATISTICAL MODEL FOR GAS MIXTURES

Keywords: Kinetic theory, gas mixtures, BGK models, moments systems subject classifications 35Q20, 35Q35

In this paper, we propose a construction of a new BGK model generalizing the Ellipsoidal Statistical Model ([2], [19]) to the context of gas mixtures. The derivation of the model is based on the introduction of relaxation coefficients associated to some moments and the resolution of a minimization problem as in ([7], [8], [10]

). We obtain in this work, an ESBGK model for gas mixtures satisfying the fundamental properties of the Boltzmann collision operator (conservation laws, H theorem, equilibrium states, . . . ) and that is able to give a range of Prandtl numbers including the indifferentiability situation.

Introduction

The complexity of the nonlinear Boltzmann operator suggests to introduce simpler kinetic models. Hence, the BGK model ( [START_REF] Gross | Model for collision processes in gases: small-amplitude oscillations of charged two-component systems[END_REF]) which is a well known simplified kinetic relaxation model has been introduced. Its features are to replace the complicated integral of collisions by a relaxation model while keeping some important physical and mathematical properties of the interaction term (conservation laws, H theorem, equilibrium states, . . . ). The interest of these models is that they are easier to handle numerically and less costly at a computational point of view. However the extension of the BGK model to multi-component gases meets fundamental difficulties. For example, the hydrodynamic limit is much more complicated since phenomena such as diffusion or thermal diffusion must be considered.

Many BGK models have been proposed up to now, but they are still unsatisfactory and have been constructed on the same manner. That is they must reproduce the rates of exchanges of impulsion and energy of the Boltzmann operator for Maxwell molecules between the different species. For example one polynomial model has been proposed in [START_REF] Garzo | A kinetic model for a multicomponent gas[END_REF]. But its relaxation function being a polynom, the nonnegativity of the distribution function is not satisfied. However this model satisfies the indifferentiability principle stated in [START_REF] Garzo | A kinetic model for a multicomponent gas[END_REF] inherited from the Boltzmann operator. That is when all the molecules have the same masses and their cross sections are equal, then the system of equations reduces to a single one by adding the distribution functions. A BGK model enjoying good mathematical properties has been derived in [START_REF] Andries | A consistent BGK-type model for gas mixtures[END_REF]. The main idea is to introduce only one BGK operator per species whose macrosopic parameters reproduce the interaction between each species with the others. But this model leads to uncorrect transport coefficients at the hydrodynamic limit ( [START_REF] Kosuge | Model Boltzmann equation for gas mixtures: Construction and numerical comparison[END_REF]). In [START_REF] Kosuge | Model Boltzmann equation for gas mixtures: Construction and numerical comparison[END_REF], the authors constructed a BGK model which coincides with the Grad moments of the linearized Boltzmann operator. However the relaxation function of the model being polynomial, the nonnegativity of the distribution function is not sure.

Hence the aim of this paper is to construct a BGK model that is able to recover correct transport coefficients and enjoying good mathematical properties. The question is important because in general the authors compute the hydrodynamic limit from their own model and eventually compare with the correct fluid model ( [START_REF] Hamel | Two-fluid hydrodynamic equations for a neutral, disparate-mass, binary mixture[END_REF]). In the situation of single and monospecies a BGK model leading to the correct Prandtl number has been proposed in ([19]). This model is called Ellipsoidal Statistical Model (ESBGK). But the proof of the H theorem has been shown later on ( [START_REF] Andries | Entropy condition for the ES BGK model of Boltzmann equation for mono and polyatomic gases[END_REF]). Recently, a new approach of such a model has been proposed in ( [START_REF] Brull | A new approach of the Ellipsoidal Statistical Model[END_REF]) based on the an entropy minimization principle ( [START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF], [START_REF] Schneider | Entropic approximation in kinetic theory[END_REF], [START_REF] Junk | Maximum entropy for reduced moment problems[END_REF]). The philosophy of the derivation consists in defining different relaxation rates associated to some moments of the distribution function. These rates intuitively describe the way in which these moments vanish for the trend to equilibrium. For example, defining a relaxation rate on the viscous stress tensor allows to reproduce viscosity and heat flux. Generally, the moment that is dissipated is considered with respect to the polynomial function that is involved in the computation of the corresponding transport coefficient that we aim to fit. In ( [START_REF] Brull | On the Ellipsoidal Statistical Model for polyatomic gases[END_REF]), this approach has been generalized to the polyatomic context to propose a new construction of the Ellipsoidal Statistical Model ( [START_REF] Andries | Entropy condition for the ES BGK model of Boltzmann equation for mono and polyatomic gases[END_REF]). In that case, another constraint of equalization of translational and internal temperature allows to fit simultaneously shear viscosity, volume viscosity and heat flux. In the context of gas mixtures, things are more complicated. Indeed, there are more transport coefficients, that is the Fick, Soret, Dufour, Fourier and Newton coefficients (see 4.6). In ( [START_REF] Brull | Derivation of BGK models for mixtures[END_REF]), this method has been considered for gas mixtures by deriving a BGK model that is able to recover the correct Fick law. In this situation, several relaxation parameters are involved and are shown to be directly related to the non zero eigenvalues of the diffusion part of the Onsager matrix. Moreover even if Soret, Dufour and Fourier coefficients are not correct, the structure of the fluxes is recovered. In ( [START_REF] Brull | Derivation of a BGK model for reacting gas mixtures[END_REF]), this model has been generalized the slow chemistry framework where the mechanical process dominates the chemical one. The BGK model is then defined as a splitting between the mechanical model of ( [START_REF] Brull | Derivation of BGK models for mixtures[END_REF]) and the chemical one constructed in ( [START_REF] Groppi | A kinetic relaxation approach to fast reactive mixtures: shock wave structure[END_REF]). Next a Navier-Stokes system is constructed for the weak reaction regime from a Chapman-Enskog procedure and the Fick matrix is shown to be recovered by the model.

The present paper is organized as follows. In section 2, we clarify some notations and recall some backgrounds about the Boltzmann operator for gas mixtures. The third section is devoted to the construction and the fundamental properties of the ESBGK model for gas mixtures. Firstly, we introduce the set of constraints. Next, we solve a minimization problem under moments constraints and we prove that the resulting model satisfies a H theorem and the indifferentiability principle. In section 4, we perform a Chapman-Enskog expansion up to Navier-Stokes level. At this step, the two relaxation coefficients are computed in order to get an exact viscosity and an exact heat flux. Finally section 5 deals with a conclusion and perspectives to this work.

2. Notations and some backgrounds about kinetic equations for gas mixtures 2.1. Boltzmann equation for gas mixtures In this paper, we consider a gas mixtures with p components. This mixture can be describded at the kinetic level by its distribution function f = (f 1 , . . . , f p ), where f i (t, x, v) represents the distribution function of the species i. Here, t corresponds to the time variable, x to the space variable and v to the velocity variable. The distribution function for the mixture evolves according to the Boltzmann equation for a p component gas

∀i ∈ [1, p] , ∂ t f i + v • ∇ x f i = k=p k=1 Q ki (f k , f i ) := Q i (f, f) , (2.1)
where

Q ki (f k , f i ) = R 3 ×S 2 (f k (w * ki ) f i (v * ki ) -f k (w) f i (v)) σ ik (ω ω ω.V, V ) V dwdω ω ω.
Here Q ki is the Boltzmann collision operator between molecules of species i and k and σ ik = σ ki is the differential cross section which depend on the interaction potential between species i and k. Finally V = w-v is the relative velocity. The post collisional velocities are given by

v * ki = v -2 m k m i + m k ((v -w) • ω ω ω) ω ω ω, w * ki = w + 2 m i m i + m k ((v -w) • ω ω ω) ω ω ω.
Those equations satisfy the conservation of momentum and energy at a microscopic level

m i v + m k w = m i v * ki + m k w * ki , m i v 2 + m k w 2 = m i v * ki 2 + m k v * ki 2 .

Macroscopic quantities

We denote with n i , ρ i , u i , E i , T i and E i the macroscopic quantities representing respectively the number density, density, average velocity, energy per unit volume, energy per particle and finally temperature of a given species i. They are defined by the following relations

n i = R 3 f i dv, ρ i = m i n i , n i u i = R 3 v f i dv, E i = 1 2 ρ i (u i ) 2 + n i E i , E i = 3 2 k B T i = m i 2n i R 3 (v -u i ) 2 f i dv,
where m i represents the molecular mass of the species i and k B is the Boltzmann constant. In the same way, we define the hydrodynamic quantities for the mixture

n = p k=1 n k , ρ = p k=1 ρ k , ρu = p k=1 ρ k u k , nE + ρ 2 u 2 = E = p k=1 E k , E = 3 2 k B T.
n i , u et T being given, we define the Maxwellian distribution functions as

∀i ∈ [1, p], M i = n i (2πk B T /m i ) 3 2 exp(- m i (v -u) 2 2k B T ).
The entropy associated to the distribution function f of the mixture is defined by

H(f ) = p i=1 R 3 (f i ln(f i ) -f i ) dv. (2.2)

Other considerations.

Using the above notations we note as L 2 (M) the set of measurable functions Ψ = (ψ 1 , • • • , ψ p ) such that:

Ψ 2 := i=p i=1 R 3 ψ 2 i M i < +∞.
This space is equipped using its natural dot product:

Ψ, Φ = i=p i=1 R 3 ψ i φ i M i dv.
In the sequel we often use the dot product notation. This notation is valid both for list of p scalar functions and for p tensorial functions. For instance if Ψ, Φ are two lists of p (symmetrical) tensorial functions ψ i , φ i , i ∈ [1, p], then the dot product notation ψ, φ should be understood as

ψ, φ := i=p i=1 R 3 ψ i ⊗ φ i M i dv,
where ⊗ denotes the usual tensorial product. Note also the following convention: assume that V is a list of p vectors and T a list of p square matrixes. Then a notation like αV, αT (where α is a scalar) means that the scalar α is distributed on each line of the list of vectors or matrix. Besides if β β β is a vector, then a notation like β β β • V means that we distribute the dot product by β β β on each line of the vector V. The same way, if γ γ γ is a square matrix, then a notation like γ γ γ : T means that we distribute the (total) dot product (between matrix) by γ γ γ on each line of the line of T. Finally, if s is a tensor and S a list of p tensors, then a notation like s ⊗ S means that we distribute the tensor product by s on the left on each line of the list S. Definition 2.1. The set K belonging to L 2 (M) of the collision invariants of the Boltzmann operator (2.1) is spanned by the following list of functions:

     1 0 . . . 0      , • • • ,      0 0 . . . 1      ,      m 1 v x m 2 v x . . . m p v x      ,      m 1 v y m 2 v y . . . m p v y      ,      m 1 v z m 2 v z . . . m p v z      ,      m 1 v 2 m 2 v 2 . . . m p v 2      . (2.3)
This space is of dimension p + 4. We denote by φ φ φ l , l ∈ [1, p + 4] the list of these p + 4 elements. Contrarily to the case of monatomic gas there exists a space C of moments of degree 1 in velocity which is not conserved. This space is defined as follows Definition 2.2. Let C i be the vector whose i th component is vu and others are 0. Denote by P K the orthogonal projection on K and I the identity operator. Then we define C as the space generated by the vectors (

I -P K ) (C i ) , i ∈ [1, p].
In ( [START_REF] Brull | Derivation of BGK models for mixtures[END_REF]), we have shown the following lemma which exhibits a basis of C. Lemma 2.3. The family (I

-P K ) (C i ) , i ∈ [1, p] is composed of p -1 independant
"vectors" and as a consequence the dimension of C is 3 (p -1). For the Chapman-Enskog procedure, introduce the Sonine polynomials A and B whose i th component write

A i = m i v ⊗ v - 1 3 v 2 Id , B i = (v -u) 1 2 m i (v -u) 2 - 5 2 k B T .
(2.4)

3. Construction and properties of the model 3.1. Fundamental properties of kinetic operators As in ( [START_REF] Andries | A consistent BGK-type model for gas mixtures[END_REF], [START_REF] Brull | Derivation of BGK models for mixtures[END_REF]), the BGK model R (f) constructed in this paper consists in one operator per species as

R (f) = (R 1 (f) , . . . , R p (f)). (3.1)
We present now the fundamental properties coming from the Boltzmann operator for gas mixtures that our BGK model has to satisfy.

1. Collisional invariants ∀f, f i ≥ 0, ∀φ φ φ, i=p i=1 R 3 R i (f) φ i dv = 0 ⇔ φ φ φ ∈ K. (3.2)
2. H-theorem: for any list of nonnegative functions f = (f 1 , . . . , f p ) there holds

i=p i=1 R 3 R i (f) ln (f i ) dv ≤ 0. (3.3)
3. Equilibrium states: the equality holds in the above equation if and only if f is at thermodynamic equilibrium i.e there exists macroscopic values n 1 , ..., n p , u, T such that

∀i ∈ [1, p] , f i = M i .
In such a case we denote f = M. Moreover M is the only set of functions such that

R i (f) = 0. (3.4)
3.2. Space of constraints f being given, G is researched as the solution of an entropy minimisation problem under moments constraints. The space of constraints that is used in this paper is given in the following definition. Definition 3.1. The space of constraints K(f ) associated to f , is defined by the set of functions g = (g 1 , . . . , g p ) s.t. g ≥ 0 a.e. satisfying ∀l ∈ {1, p + 4},

p i=1 R 3 φ l i (g i -f i ) dv = 0 (3.5)
to obtain the conservation properties and the relaxation constraint

p i=1 R 3 λ(g i -f i ) A i (v -u) dv = -λ 1 p i=1 R 3 f i A i (v -u) dv, (3.6)
where λ > 0 and λ 1 > 0 are free parameters that are devoted to fit the transport coefficients. Like in ( [START_REF] Brull | A new approach of the Ellipsoidal Statistical Model[END_REF], [START_REF] Brull | On the Ellipsoidal Statistical Model for polyatomic gases[END_REF], [START_REF] Brull | Derivation of BGK models for mixtures[END_REF]), they are fixed during the Chapman-Enskog procedure.

The constraint (3.6) can be rewritten as

1 ρ p i=1 R 3 m i g i (v -u) ⊗ (v -u) dv = (1 - λ 1 λ )Θ + λ 1 λ T, (3.7)
where Θ represents the pressure tensor defined by

Θ = p i=1 R 3 m i (v -u) ⊗ (v -u) f i dv. (3.8)
In order to compare our model with the one derived in ([2], [START_REF] Brull | A new approach of the Ellipsoidal Statistical Model[END_REF]), we use the change of variable: ν = 1 -λ1 λ . So the constraint (3.7) on g can be rewritten

1 ρ p i=1 R 3 m i g i (v -u) ⊗ (v -u) dv = νΘ + (1 -ν)T = τ. (3.9)
Remark 1. In ( [START_REF] Brull | Derivation of BGK models for mixtures[END_REF]), the relaxation constraints (3.6) are written on a basis of C instead of the Sonine polynomials (2.4). This basis is obtained from the diagonalization of the diffusion part of the Onsager matrix.

G i = ρ i 2π det(τ ) exp - 1 2 v -u, τ -1 (v -u) . (3.11)
Conversely, if the minimization problem (3.10) has a unique solution for any nonnegative function

f ∈ L 2 1 , then ν ∈ [-1 2 , 1[. Proof.
To prove the existence of the distribution function (3.11), we have to solve a minimization problem under moments constraints. As in ( [START_REF] Brull | Derivation of BGK models for mixtures[END_REF]), we follow the methodology of Junk ([20]). Hence we consider the set

Λ = {(α l ) l , (λ r ) r / exp p+3 l=1 α l φ l i + λ 1 v 2 1 + λ 2 v 2 2 + λ 3 v 2 3 + λ 4 v 1 v 2 + λ 5 v 1 v 2 + λ 6 v 2 v 3 < +∞}. (3.12)
If this set is open, then the minimization problem has a unique solution as soon as K(f ) = ∅. This last condition is fullfilled because G defined by (3.11) belongs to K(f ).

Hence it remains to show that (3.12) is open. Proving this assertion is equivalent to show that the signature of the quadratic form

q(v) = λ 1 v 2 1 + λ 2 v 2 2 + λ 3 v 2 3 + λ 4 v 1 v 2 + λ 5 v 1 v 2 + λ 6 v 2 v 3 (3.13)
is (0, 3). Introducing the matrix associated to the quadratic form q, M = (m ij ) ij , the Gauss reduction of q (3.13) writes

q(v) = m 11 (v 1 + m 12 m 11 v 2 + m 13 m 11 v 3 ) 2 + α(v 2 + 1 α (m 23 - m 13 m 11 )v 3 ) 2 + γv 2 3 , m 11 < 0, α < 0, γ < 0, (3.14)
where

α = m 22 - m 2 12 m 11 , γ = m 33 - m 2 13 m 11 - 1 α (m 23 - m 12 m 13 m 11 ) 2 .
Hence for ε > 0 small enough, we can show that the signatures of the quadratic form associated to the matrix whose terms are m ij + ε and m ij -ε is (0, 3). This means that the set (3.12) is open. The proof of the second part of the theorem is analogous to the one given in the monospecies setting. So we refer to ( [START_REF] Brull | A new approach of the Ellipsoidal Statistical Model[END_REF]) for more details.

H Theorem

In this section, we show that the present BGK model is entropic for -1 2 ≤ ν < 1. Theorem 3.3. For -1 2 ≤ ν < 1 and G i defined by (3.11), we have the following entropy dissipation property

D(f ) = p i=1 (G i -f i ) ln(f i ) dv ≤ 0. Moreover, D(f ) ≤ 0 for -1 2 ≤ ν < 1 with equality iff f i = M i , i = {1, .
. . , p}. Proof. The proof is similar to the one given in ( [START_REF] Brull | A new approach of the Ellipsoidal Statistical Model[END_REF]) but adapted to a gas mixture framework. For the sake of clarity, we give the proof in details. By using the convexity of x → x ln(x) -x, we get that

p i=1 R 3 (G i,ν -f i ) ln f i dv ≤ H(G ν ) -H(f ),
where G i,ν is the function defined in (3.11) and ν represents its dependance w.r.t

ν ∈ [-1 2 , 1[. Moreover, G 1-ε being the result of a minimization problem under constraints of f , H(G 1-ε ) ≤ H(f ). So p i=1 R 3 (G i,ν -f i ) ln f i dv ≤ H(G ν ) -H(G 1-ε ). But a direct computation of H(G ν ) gives H(G ν ) = p i=1 ρ i ln ρ i √ 2π detT - 5 2 ρ.
Moreover the set of constraints that is used for the definition of G ν are linear, the application ν → G i,ν is convex. So to obtain the entropy dissipation property, it is enough to show that

H(G ν ) ≤ H(G 1-ε ). H(G ν ) -H(G 1-ε ) = 1 2 ρ (ln(det((1 -ε) Θ + εT Id)) -ln(det(νΘ + (1 -ν)T Id))) .
But this term has the same sign as

p(ν) = det((1 -ε) Θ + εT Id)) -det(νΘ + (1 -ν)T Id).
p(ν) can be rewritten as

p(ν) = 3 i=1 (θ i + ε(T -θ i )) + 3 i=1 (ν(T -θ i ) -T ) .
Hence by setting θi = θ i /T , we get

p(ν) = T 3 3 i=1 ( θi + ε(1 -θi )) + 3 i=1 (ν(1 -θi ) -1) . Then p(- 1 2 ) = 3 i=1 ( θi + ε(1 -θi )) + 3 i=1 (- 1 2 (1 -θi ) -1) .
Therefore after some computations p(-

1 2 ) = 3T 3 8 3( θ1 θ2 θ3 -θ1 θ2 -θ1 θ3 -θ2 θ3 ) + εl( θ1 , θ2 , θ3 ) with l( θ1 , θ2 , θ3 ) = (1 -θ1 ) θ2 θ3 + θ1 (1 -θ2 ) θ3 + θ1 θ2 (1 -θ3 ) + ε((1 -θ1 )(1 -θ2 ) θ3 + (1 -θ1 ) θ2 (1 -θ3 ) + θ1 (1 -θ2 )(1 -θ3 )) + ε 2 (1 -θ1 )(1 -θ2 )(1 -θ3 ).
By using Lagrange multipliers, we get that the maximum of p(-1 2 ) is obtained for θ1 = θ2 = θ3 = 1 and yields p(- 1 2 ) = 0. 3.4. Indifferentiability condition In this section, we prove that our model satisfies the indifferentiability property stated in [START_REF] Garzo | A kinetic model for a multicomponent gas[END_REF] and recalled in introduction. Proposition 3.4. If all the molecular masses are equal [START_REF] Brull | A new approach of the Ellipsoidal Statistical Model[END_REF]). Hence the model satisfies the indifferentiability principle. The proof is trivial and performed by adding term by term each component of kinetic equations.

(∀ i ∈ {1; p}, m i = m), then the distribution function f = p i=1 f i is solution of the Ellipsoidal Statistical Model for monospecies ([2],

Chapman-Enskog expansion

4.1. The Navier-Stokes system for gas mixtures The Navier-Stokes system for a mixture of p components reads

∀i ∈ [1, p] , ∂ t n i + ∇ • (n i u + J i ) = 0, (4.1)
∂ t (ρ u) + ∇ • (P + ρ u ⊗ u + J u ) = 0, (4.2)
∂ t E + ∇ • (Eu + P [u] + J u [u] + J q ) = 0, (4.3)
where J i , J u J q are respectively the mass, momentum and heat fluxes. These fluxes can be expressed either from thermodynamics of irreversible process (TIP) ( [START_REF] De Groot | Mazur Nonequilibrium Thermodynamics[END_REF]) or from the phenomenological point of view. Of course these two approaches are strictly equivalent according to the definition of the chemical potential (4.5). As it has been explained in ( [START_REF] Brull | Derivation of BGK models for mixtures[END_REF]), the Chapman-Enskog procedure leads naturally to fluxes defined from TIP assumption whereas coefficients appearing in phenomenological approach are obtained from experiments. We firstly describe the transport coefficients obtained from the TIP point of view. In this case, mass, impulsion and energy fluxes write

J i = j=p j=1 L ij ∇ -µj T + L iq ∇ 1 T , J q = j=p j=1 L qj ∇ -µj T + L qq ∇ 1 T , J u = L uu D (u) , (4.4)
where

D(u) = 1 2 ∇u + (∇u) T - 1 3 (∇ • u)I
represents the Reynolds tensor and

- µ i T = k B ln n i - 3 2 ln 2πk B T m i (4.5)
is the chemical potential.

By using the phenomenological approach, mass and impulsion fluxes write w.r.t densities and temperature gradients as

J i = j=p j=1 D ij ∇n j + D iT ∇T, J q = j=p j=1 D qj ∇n j -D qq ∇T, (4.6)
where D ij and D iT , D qj and D qq respectively denote the Fick, Soret, Duffour and Fourier coefficients.

Computation of the relaxation coefficients λ and λ 1

We determine now the coefficients λ and λ 1 in such a way that the Chapman-Enskog expansion up to order one leads to the true viscosity and the true heat flux. That is the ones that have been computed from the Boltzmann equation or that have been obtain experimentally.

Hence we consider the system of kinetic equations

∂ t f i + v • ∇f i = λ ε (G i -f i ), i ∈ {1, . . . , p}, (4.7)
where ε is proportional to the Knudsen number. For each component i, the distribution function f i is expanded according to

f i = M i + εf (1) i
+ . . . . For more details on Chapman-Enskog expansion for gas mixtures, we refer to ( [START_REF] Chapman | The mathematical theory of non uniform gases[END_REF], [START_REF] Brull | Derivation of BGK models for mixtures[END_REF]). According to ( [START_REF] Brull | Derivation of BGK models for mixtures[END_REF]), the fluxes write

J i = f (1) , C i , J q = f (1) , B , J u = f (1) , A . (4.8)
From (4.7), G which depends on f is expanded around the Maxwellian distribution M according to G = M + εG (1) + . . . .

Using that G i (M) = M i , the expression of G

(1) i is given by

G (1) i = DG i (M).f (1) i = lim τ →0 G i (M + τ f 1 ) -M i τ . (4.9)
So the order 0 in (4.7) leads to 1) f (1) ). (4.10)

∂ t M + v • ∇M = λ(G (
But, from ([10]), we get the relation

(∂ t + v.∇)M = p j=1 k -1 B (I -P K )(C j ).∇(- µ j k B T ) + 1 k B T A : D(u) + 1 k B B. ∇T T 2 M, (4.11)
the Sonine polynomial B has been defined in (2.4). Therefore f (1) has the following expression

f (1) = G (1) - p j=1 k -1 B (I -P K )(C j ).∇(- µ j k B T ) + 1 k B T A : D(u) + 1 k B B. ∇T T 2 M.
(4.12)

Relation (3.6) considered at order 1 for g = G leads to

p i=1 R 3 λ(G (1)
i -f

(1)

i ) A i (v -u) dv = -λ 1 p i=1 R 3 f (1) i A i (v -u) dv.
From orthogonality reasons, it comes that

p i=1 R 3 f (1) i A i (v -u) dv = - 1 k B T λ 1 p i=1 R 3 A i (v -u) : A i (v -u)M i dv D(u) = - nk B T λ 1 D(u).
So we define λ 1 as λ 1 = nk B T µ . By setting

λ = λ 1 1 -ν = ρT µ(1 -ν) , (4.13)
we get the relaxation operator (3.1)

R i (f ) = ρT µ(1 -ν) (G i -f i ), i ∈ {1, . . . , p}.
In order to determine the heat flux, we use (4.8). Multiply (4.10) by B(vu) and integrate. From relation (4.9) defining G (1) , we can show that

R 3 G (1) B(v -u) dv = 0.
By using the orthogonality of B w.r.t. A and (I -P K )(C j ), it holds that

p i=1 R 3 B i (v -u)f (1) i dv = - 1 3k B λ p i=1 R 3 B i (v -u).B i (v -u) M i dv ∇T T 2 .
A direct computation gives

p i=1 R 3 B i (v -u).B i (v -u) M i dv = 15 2 k 3 B T 3 p i=1 n i m i .
Hence, from (4.4) the heat flux L qq has the following expression

L qq = 5 2λ k 2 B T 3 p i=1 n i m i . ( 4

.14)

In order to compare the model constructed in the present paper, with the monospecies model derived in ( [START_REF] Brull | A new approach of the Ellipsoidal Statistical Model[END_REF]), we set κ = Lqq T 2 . Hence, Prandtl number writes

P r = 5 2 Rµ κ = R λ λ 1 n k B p i=1 n i m i .
Therefore by using the definition (4.13) of the relaxation coefficient λ and R = k B /m, we get

P r = R 1 -ν n k B p i=1 n i m i = 1 1 -ν 1 p i=1 n i n m m i , ( 4
]), P r = 1/(1 -ν). Hence P r = 2/3 is obtained for ν = -1/2 and this totally determines λ.

Next the other transport coefficients are computed. Their values are not correct, but the shape of the fluxes is recovered. In property 1, we collect the values of all the transport coefficients. Property 1. The transport coefficients are given by

L ij = 3T n i k B m i δ ij , L iq = L qi = 0, L qq = 5 2λ k 2 B T 3 p i=1 n i m i , L uu = -µ.
Proof. According to Chapmann-Enskog expansion, mass flux writes J i = f (1) , C i , where f (1) is computed from (4.12). But, we have A, C i = 0, B, C i = 0. Then L iq = L qi = 0 and

L ij = - 1 k B R 3 C i .C j M i dv = - 1 k B R 3 (I -P K ) (C i ). (I -P K ) (C j ) M i dv.
So the expression of L ij follows.

Conclusion and perspectives

In this paper, we have introduced an Ellipsoidal Statistical Model for gas mixtures following the moments relaxation method introduced in ([7], [START_REF] Brull | On the Ellipsoidal Statistical Model for polyatomic gases[END_REF], [START_REF] Brull | Derivation of BGK models for mixtures[END_REF]). The resulting model fullfills two important features. Firstly, it enjoys the fundamental properties of the Boltzmann operator (H theorem, nonnegativity of the distribution function, conservation laws, indifferentiability condition, . . . ) which is only satisfied by few models ( [START_REF] Andries | A consistent BGK-type model for gas mixtures[END_REF]). Moreover, its construction is based on the right hydrodynamic coefficients and the model is able to reach a certain rage of Prandl numbers. This range includes the correct Prandtl number in the situation of indifferentiability condition.

For concrete application, this methodology should be generalized to the situation of polyatomic gases in order to fit Prandtl number. Hence the following step of this work should be the derivation of an Ellipsoidal Statistical Model for polyatomic gas mixtures. For example the atmosphere can be considered as a binary mixture of Oxygen and Nitrogen. In that case, the distribution function depends on an additional continous variable, the internal energy ( [START_REF] Bourgat | Microreversible collisions for polyatomic gases and Boltzmann's theorem[END_REF]), collecting rotational, vibrational, translational energy and belonging to R + . In ( [START_REF] Brull | On the Ellipsoidal Statistical Model for polyatomic gases[END_REF]), we have derived from the methodology described in this paper an ESBGK model for a polyatomic single gas. So by taking advantage of this work, we plan to generalize this construction to a gas mixture framework. Before this, the Chapman-Enskog expansion for polyatomic gas mixtures has to be clarified. This work is postponed to a future paper.

Theorem 3 . 2 .

 32 For any nonnegative function f = (f 1 , . . . , f p ) and ν ∈ [-1 2 , 1[, the tensor τ defined by (3.9) is symetric defined positif and the minimization entropy problemG = min g∈K(f ) H(g)(3.10)has a unique solution G = (G 1 , . . . , G p ) where G i writes

( 4 . 4 )

 44 and (4.6) are shown to be equivalent by plugging the expression (4.5) of -µi T into (4.4).