Hydrogen is an essential energy carrier for a successful ecological and energy transition. However, most hydrogen is produced by cracking hydrocarbons of fossil origin. Only 1% of the hydrogen currently produced comes from the electrolysis of water. Hydrogen from electrolysis is too expensive to produce. The objective of this thesis is to study electrogenerated bubbles in order to identify microfluidic aspects that could contribute to the improvement of water electrolysis. Bubbles act as an electrical insulator. By covering the electrode, they reduce the efficiency of electrolysis. Therefore, once they nucleate, they must quickly detach from the electrode. There is currently no consensus on the phenomena influencing the growth and detachment of bubbles at the microfluidic scale. Among others, there are still many uncertainties on how the wettability of the electrodes, and the Marangoni effect of thermal or solutal origin influence the nucleation, growth and detachment of bubbles. To this end, the mathematical and numerical basis for the simulation of a two-phase fluid was reviewed. In order to study such a phenomenon numerically, it is necessary to be able to simulate the surface tension variations along a liquid-gas interface, to integrate the mass transfer across the interface from the dissolved species present in the electrolyte to the gas phase, and to take into account the moving contact line. The use of the continuous surface force (CSF) model in the volume of fluid (VOF) framework is known to introduce non-physical velocities, called parasitic currents. The use of an alternative model based on the height function (HF) approach has been developed and tested. Its use limits the spurious currents and makes the VOF methodology suitable for the study of Marangoni currents at the interface of an electrogenerated bubble. A correlation for determining the growth rate of a bubble by integrating the Marangoni currents and based on the penetration theory was developed and compared to the Epstein-Plesset relation. A dimensionless study was conducted to relate the Sherwood number representing interfacial mass transfer to the Marangoni number. There are too many unknowns to draw definitive conclusions. However, the implementation of a contact line model in the numerical model could remove many uncertainties. The work done in this thesis to develop a holistic model is a first step towards Y mass fraction 6 Nomenclature Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022 Chapter 2 Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022 Chapter 2 Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022

using an inverse problem to determine the unknowns that need to be removed to optimise bubble detachment.

Résumé

L'hydrogène est un vecteur énergétique essentiel pour réussir la transition écologique et énergétique. Cependant, la majeure partie de l'hydrogène est produite par le craquage d'hydrocarbures d'origine fossile. Seul 1% de l'hydrogène actuellement produit provient de l'électrolyse de l'eau. L'hydrogène issu de l'électrolyse est trop cher à produire. L'objectif de cette thèse est d'étudier les bulles électrogénérées afin d'identifier les aspects microfluidiques qui pourraient contribuer à l'amélioration de l'électrolyse de l'eau. Les bulles agissent comme un isolant électrique. En recouvrant l'électrode, elles réduisent l'efficacité de l'électrolyse. Par conséquent, après leur nucléation, elles doivent se détacher rapidement de l'électrode. Il n'existe actuellement aucun consensus sur les phénomènes influençant la croissance et le détachement des bulles à l'échelle microfluidique. Entre autres, il existe encore de nombreuses incertitudes sur la façon dont la mouillabilité des électrodes, et l'effet Marangoni d'origine thermique ou solutale influencent la nucléation, la croissance et le détachement des bulles. Dans ce but, les bases mathématiques et numériques nécessaires à la simulation d'un fluide diphasique ont été passées en revue. Afin d'étudier numériquement un tel phénomène, il est nécessaire de pouvoir simuler les variations de tension superficielle le long d'une interface liquide-gaz, d'intégrer le transfert de masse à travers l'interface des espèces dissoutes présentes dans l'électrolyte vers la phase gazeuse, et de prendre en compte la ligne de contact mobile. L'utilisation d'un modèle CSF (continuum surface force) dans le cadre du VOF (Volume of fluid) est connue pour introduire des vitesses non physiques, appelées courants parasites. L'utilisation d'un modèle alternatif basé sur l'approche des fonctions de hauteur a été développée et testée. Son utilisation limite les courants parasites et rend la méthodologie VOF adaptée à l'étude des courants de Marangoni à l'interface d'une bulle électrogénérée. Une corrélation permettant de déterminer le taux de croissance d'une bulle en intégrant les courants de Marangoni et basée sur la théorie de la pénétration a été développée et comparée à la relation de Epstein-Plesset. Une étude sans dimension a été menée pour relier le nombre de Sherwood représentant le transfert de masse interfacial au nombre de Marangoni. Il y a trop d'inconnues pour tirer des conclusions définitives. Cependant, l'implémentation d'un modèle de ligne de contact dans le modèle numérique pourrait lever de nombreuses incertitudes. Le travail effectué dans cette thèse pour développer un modèle holistique est un premier pas vers l'utilisation d'un problème inverse pour déterminer les inconnues à lever pour optimiser le détachement des bulles.

Nomenclature

GREC λ elec electrical conductivity [ S • m -1 ] α volume fraction δ 1 diffusion layer thickness [m] ϵ void ratio η cathode anode overpotential [V] η cathode cathode overpotential [V] γ surface tension [N • m -1 ] κ curvature [m -1 ] λ thermal conductivity [W • m -1 • K -1 ] ν kinematic viscosity [m 2 • s -1 ] Φ electric potential [V] Φ B current efficiency ρ density [kg • m -3 ]
c p heat capacity [J • kg -1 • K -1 ] c bulk bulk concentration [mol • m -3 ] c sat saturation concentration [mol • m -3 ] D diffusion coefficient [m 2 • s -1 ] D conv Diffusion coefficient chaotic motion [m 2 • s -1 ] D H 2 Diffusion of H 2 in electrolyte [m 2 • s -1 ] D th thermal diffusivity [m 2 • s -1 ] F Faraday constant [C • mol -1 ] f G gas evolution efficiency H Henry's constant [N • m • mol -1 ] h electrolyte thickness of the electrolyte [m] J material flow density [mol • m 2 • s -1 ] j current density [A • m -2 ] j avg average current density [A • m -2 ] j lim limiting current density [A • m -2 ] k mass transfert coefficient [m • s -1 ] k B Boltzmann constant [J • K -1 ] k b bubble induced mass transfert coefficient [m • s -1 ] k dif f diffusion induced mass transfert coefficient [m • s -1 ] M g molar mass of gas [kg • mol -1 ] M H 2 molar mass of dihydrogenn [kg • mol -1 ] p b pressure inside the bubble [N • m -2 ] p e pressure in the electrolyte [N • m -2 ] R bubble radius [m] R c critical radius [m] R g universal gas constant [J • K -1 • mol -1 ]
V bubble volume [m 3 ] W work [N • m]
w elec energy required for the gas to form [J] Introduction 1

This thesis is the subject of a joint agreement between the University of Canterbury and the University of South Brittany. The subject is "Study of electrogenerated two-phase and microfluidic flows". It follows Damien LE BIDEAU's ADEME (Environment and Energy Management Agency) thesis which concerned the study of two-phase electrolysis for the improvement of hydrogen production. While Damien Le Bideau's thesis focused on macroscopic parameters at the scale of an electrolyzer, the aim of this thesis is to investigate the phenomena involved at the microfluidic scale in the production of electrogenerated bubbles.

The challenge of the XXI century is to make a successful ecological and energy transition from a world vision where energy resources are abundant and growth is unlimited to a world where international tensions to acquire the remaining resources are increasingly felt and where ecological constraints such as global warming are jeopardizing human development.

There are at least three main areas of research that can influence the success of the energy transition. The reduction of energy consumption in buildings and industry, the improvement of renewable energy resources, and the improvement of energy storage.

One of the main obstacles to the use of renewable energies are the spatial and temporal constraints related to their production. For example, the day-night cycles impose on photovoltaic fields a production during the day that reaches its peak at noon. This production is not in phase with consumption. Thus, the challenge here is to succeed in storing this energy during the production peaks and to restore it during the consumption peaks. As energy production is intermittent, its integration into the electricity grid cannot be done effectively. National electricity grids are not designed to receive large-scale intermittent electricity production. Production must be in line with consumption in order for the network to remain stable.

It is for this reason that the storage of the electricity produced is necessary, in order to be able to produce energy according to demand. Some of the main energy storage technologies are as follows : supracondensator (storage time around 1 second), battery (storage time around a week), "Pumped Storage Power Plants" (unlimited storage time, but storage capacity limited, depends on the capacity of the water reservoirs available in the country), and "Power to Gas" (P2G). Regarding the inter-seasonal cycles, we observe that consumption decreases in summer while photovoltaic production increases and vice versa in winter. Batteries are not adapted to store energy over such long periods. The P2G usually refers to the combined process of producing dihydrogen with the electrolysis of water from electricity from renewable energies and then methane from dihydrogen production using the Sabatier reaction [START_REF] Gruber | Integrated high-temperature electrolysis and methanation for effective power to gas conversion[END_REF]. The storage capacity is potentially unlimited, and the storage time is more than one month. The problem of self-discharge that we observe in batteries does not arise in the case of hydrogen. The properties of hydrogen are both an advantage and a disadvantage. Indeed, the amount of energy produced by the combustion of one kilogram of hydrogen is equivalent to that of about 3 kilograms of diesel, but the density of hydrogen at atmospheric pressure is extremely low compared to other fuels. It must therefore be compressed in order to store energy in an acceptable volume. A pressure of about 700 bar is required to store it in liquid form. This can lead to the transformation of hydrogen into another more easily storable gas such as methane or ammonia, but results in a loss of overall energy efficiency due to the additional transformation. The efficiency of electricity production with hydrogen is 35%, and about 25% for methane transformed from hydrogen.

As far as the alkaline electrolysis process is concerned, it can be observed that no greenhouse gases are present in the reaction products, which makes it a clean process (and the use only produces water which does not directly contribute to global warming). Hydrogen can bring flexibility to energy networks powered by intermittent renewable energies via energy storage, promote self-consumption of renewable energies by providing a local storage solution, and can also be considered for use as fuel to decarbonise the transport sector. However, the energy required as well as the cost of a production unit make hydrogen produced by water electrolysis currently expensive to be used as a large-scale storage solution. To sum up, energy efficiency is too low to be an economically competitive solution.

In 2018, the production of hydrogen amounted to 115 million tonnes worldwide, of which 70 million tonnes in pure form was used for petroleum refining and the production of ammonia for the production of fertilisers and explosives, and 45 million tonnes was used in industry without prior separation from other gases [START_REF] Iea | The future of hydrogen[END_REF]. Most of it is obtained via the cracking of hydrocarbons of fossil origin. The products of these reactions are therefore hydrogen and CO2. The production of this hydrogen contributes therefore to greenhouse gas emissions. The electrolysis of water makes it possible to obtain hydrogen and oxygen of high purity. However, its cost, which is about twice as high as that of natural gas reforming, limits its field of application. Currently only 1% of hydrogen production is supplied by water electrolysis. The production of hydrogen produces 830 million tonnes of CO2 per year or 2% of total world emissions. An example of the application of electrolysis is the Swedish "Hybrit" steel production project. By using "green hydrogen" as an ingredient and for heat production the Hybrit project produces only 25 kg of CO2 per tonne of steel instead of 1850 kg. However, to produce the equivalent of current hydrogen world production through electrolysis would require an additional 3600 TWh of electricity, which is about the European Union's electricity production or 13.5% of the world's electricity. Thus, there are many issues at stake.

The global project in which this thesis is integrated will support the design of new generation electrolyzers and fuel cells. The aim is to allow greater production capacity of dihydrogen (average current density 1 [A • m -2 ] instead of 0.5 [A • m -2 ] currently) and a cell voltage of 2[V] i.e. an energy cost of hydrogen produced around 52 [kWh.kg -1 ] of H 2 instead of 55 to 60[kWh.kg -1 ] currently. The figures quoted here are the European PAA H2020 FCH2 figures for 2017 (Game changer Water Electrolyser) and of 2018 for the large production water electrolyser . These European PAA's serve as references for the ambitions of the Brittany region.

The aim would be to give to the Brittany region an industrial and technological development path combining the arrival of new mechanical (3D printers) and electrical (pulsed processes) technologies in the field of sustainable energy based on hydrogen. In addition to the ethical aspects linked to energy sobriety and therefore innovative technological development, the economy and therefore employment, this path of progress is particularly sensitive for the region which must improve the robustness of its electrical distribution. The IRDL (Institut de Recherche Dupuy de Lôme) is in contact with off-shore wind turbine companies (Saint-Brieuc, Ailes Marines), (Groix-Belle Ile and EOLFI), hydroliennes (Naval Group, Naval Energie, SABELLA), river turbines (Guinard Energie), Photovoltaic (roof of the base of K2 submarines in Lorient with EOLFI). All these Brittany companies have a need of storage solution or conversion. Intermittent energies must be stored to ease the concerns of EDF (Électricité de France S.A., French multinationnal electric company), the French electricity supplier who claims to be "at the end of the network". At the same time Morbihan Energie acquired in 2017 a hydrogen station for power HYUNDAI cars. The concern is the price of these acquisitions: 250 k Cfor 2 kg/day of "green H 2 product".

Electrogenerated bubbles and their impact on the electrochemical process 2

The objective of this first part is to gather knowledge on how the bubbles that are produced during the electrolytic process can affect its efficiency. This thesis focuses on the microfluidic aspect and therefore only briefly addresses the influence of macroscopic aspects. The reader interested in research on larger scale aspects to improve water electrolysis can refer to the review by [START_REF] Wang | The intensification technologies to water electrolysis for hydrogen production -A review[END_REF]. There is of course a link between what happens at the bubble scale and the electrode scale. First, it is necessary to understand the basic functioning of electrolysis and then to understand the factors that influence it, by recalling the content of the work that has been done on electrolysis in the past years. This first step will lead to a detailed study of the mechanisms involved in the nucleation, growth and detachment of electrogenerated bubbles. This second step introduces the need to have a better understanding of the transport of species in the vicinity of the bubble and the interfacial mass transfer that drives the growth of the bubbles. The last section takes into account all of the micro-fluidic phenomena suspected of influencing the development of bubbles.

This first part allows us to identify the problematic of the study, and is a prerequisite to identify the needs that will allow us to carry out a numerical study of an electrogenerated bubble.

1 The electrolysis process

The history of the electrolysis process begins in 1785 when Martinus van Marum's static electricity generator was used to reduce tin and zinc from their salt. The history of water electrolysis began in 1800 when William Nicholson and Anthony Carlisle succeeded in breaking down water into hydrogen and oxygen. It was not until 1836 that Michael Faraday published his two laws of electrolysis and established the terminology associated with electrolysis (anode, cathode, electrolyte, etc.). In order to improve the efficiency of the electrolysis process, it is necessary to recall its basic principles and to know the factors that hinder its performance.

Influence of voltage and current density on process efficiency

The products of the electrolysis reaction of water are in gaseous form, which makes this process a two-phase system. The advantage of two-phase systems is that they promote mass transfer and heat transfer. Whatever the type of electrolysis studied, alkaline electrolysis of water, acid electrolysis, high temperature electrolysis, they work on the same principle.

A voltage is imposed between two electrodes located in an electrolyte (liquid or solid).

The energy supplied via electricity causes an oxidation at the anode and a reduction at the cathode. What differs between these technologies is the half-reactions and the charge carrier.

In order to understand how the efficiency of the electrolysis process is reduced by the same bubbles whose production is its purpose, it is advisable to introduce here the notions that allow to describe its functioning, taking as an example the electrolysis of water.

Energy required

By reducing the energy required to produce hydrogen, the efficiency of the process is improved. The energy required for the gas to form is governed by the following equation:

w elec = t 0 U cell Idt = t 0 U cell jS elec dt (2.1)
where w elec is the energy required for the gas to form in [J], U cell is the voltage of the electrolysis cell in [V], j the average current density in [A • m -2 ], S elec is the surface area of the electrode in [m 2 ], t is the the duration of the electrolysis [s]. Taking into account the previous equation, this amount of energy can be reduced by decreasing the voltage across the cell, or the current density.

The current density

The current density is an extensive quantity (it is a flow of charges) which can be defined at any point of the electrochemical system, i.e. at the surface of the electrodes, but also in the electrolytic medium. The current density is directly related to the amount of hydrogen produced via Faraday's law. The rate of dihydrogen production ṁH 2 in [kg • s -1 ] at the electrode can be expressed by the following equation:

ṁH 2 = j S elec 2 F M H 2 (2.2)
where F is the Faraday constant in

[C • mol -1 ], M H 2 is the molar mass of dihydrogen in [kg • mol -1 ].
Thus, it is therefore not advisable to decrease the current density, otherwise the hydrogen production will be reduced. With the current density out of the way, it remains to study the influence of the voltage.

Difference between theoretical and real cell voltage

The electrolysis of water entails a hydrogen gas release reaction on the cathode and an oxygen gas release reaction on the anode, respectively. At a temperature of 25°C and an atmospheric pressure of 1 atmosphere the reaction and standard equilibrium electrode potential E are in an alkalyne electrolyte:

• Cathode : 2H 2 O + 2e -= H 2 + 2OH -E anode = -0.83V • Anode : 2OH -= H 2 O + 1 2 O 2 + 2e -E cathode = 0.4V
By combining the two half-reactions, the total reaction is obtained:

H 2 O = H 2 + 1 2 O 2
To obtain both gaseous emissions, a theoretical electrical voltage of U theory = E anode -E cathode = 1.23V must be applied between the two electrodes. However, in practice no reaction occurs when the voltage is less than 1.6 -1.7V . The practical voltage between the electrodes is obtained from the following relationship:

U cell = U theory + |η anode | + |η cathode | + j × R ohm (2.3)
Where R ohm is the total ohmic resistance, η anode and η cathode are the two reaction overpotential at the anode and cathode, respectively.

In the absence of a bubble, the ohmic resistance of the electrolyte can be written as:

R electrolyte = h electrolyte σ (2.4)
where h electrolyte in [m] is the thickness of the electrolyte, and σ is the electrical conductivity in [S • m -1 ]. As the bubbles appear, the bubbles act as an electrical insulator and this resistance increases.

The overpotential can be defined as the difference between the theoretical thermodynamic voltage required for the half-reaction to take place and the voltage observed experimentally.

It increases with growing current density, as described by the Tafel's empirical equation [START_REF] Wang | The intensification technologies to water electrolysis for hydrogen production -A review[END_REF]:

η = a + b log(j) (2.5)
where boths a and b are Tafel constants. a depends on the properties and surface structure of electrode materials.

There are two types of overpotential: the activation overpotential, and the concentration overpotential.

Activation overpotential

The activation overpotential η act represents the activation energy of the electrochemical reaction taking place at the cathode and the anode. It is caused by the resistance against the reaction at the electrolyte-electrode interface. This overpotential increases logarithmically with the current density as shown in equation (2.5) and is dependent on the electrode material used. To reduce this overpotential, it is possible to use a suitable electrocatalyst. An electrocatalyst is a material that offers a low activation path for a given electrochemical reaction. The catalytic activity depends on the electronic configuration of the catalyst and the surface structure (nanometric, micrometric structure).

Concentration overpotential

The concentration overpotential can be caused either by a lack of reagent or too high concentration of products. In the case of a concentration overpotential due to a high concentration of product, the electrolyte close to the electrode quickly reaches the saturation of the dissolved gas and largely exceeds it, which makes the electrolyte sursaturated. This type of overpotential is easily measurable and has been studied in various publications by Vogt et al., Dukovic et al., as well as Gabrielli et al. [Vogt, 1990a;[START_REF] Dukovic | The Influence of Attached Bubbles on Potential Drop and Current Distribution at Gas-Evolving Electrodes[END_REF][START_REF] Gabrielli | Potential drops due to an attached bubble on a gas-evolving electrode[END_REF]]. The concentration overpotential is usually calculated by the following equation:

η conc = R g T n e F ln c e c sat (2.6)
where c e is the dissolved gas concentration gas closest to the electrode in [mol • m -3 ], c sat is the saturation concentration of the dissolved gas at a pressure of 1 atm in [mol • m -3 ], R g is the universal gas constant in [J • mol -1 • K -1 ], T is the temperature, and n e is the number of electrons transferred to form one gas molecule, n e = 2 for dihydrogen and n e = 4 for dioxygen. Since c e varies with the distance from the adhering bubble, c e can be assimilated to an average concentration of dissolved gas close to the electrodes and not covered by the bubbles. Assuming that the number of moles of gas produced at the electrode by the reaction per unit time and area is equal to that transported from the electrode to the bulk, the current density can be expressed as a function of the concentration at the electrode:

j n e F = k(c e -c bulk ) (2.7) c bulk is the concentration in the bulk [mol • m -3 ],
and k is the mass transfert coefficient in [m • s -1 ] describing the mass transfer of species near the electrode. The concentration c e increases until it reaches a maximum at bubble nucleation. When the current density becomes higher then the bubbles start to participate to the mass transfer. Knowing k ,j , c bulk , c e can be estimated, which allows to calculate η conc using equation (2.6). By combining the two previous equations Vogt et al. obtain [Vogt, 1990a]:

η conc = R g T n e F ln j n e F kc sat + c bulk c sat (2.8)
Assuming a constant temperature, equation (2.8) shows that the value of this overpotential is not only a function of the current density, but also of the bulk concentration and the flow of the electrolyte.

Vogt et al. studied the overpotential for current densities from 1 to 10 4 A • m -2 [Vogt, 1990a]. The overpotential increases until current density reaches a value around 1000 A • m -2 where the mass transfer is so increased by the agitation of the bubbles that the concentration of dissolved gas can no longer increase and thus the overpotential reaches a plateau. Vogt et al. define two regimes, one controlled by macroconvection (convection induced by the electrolyte flow) and the other controlled by microconvection (convection induced by the appearance of bubbles at the electrode). For currents density below 10 A • m -2 , macroconvection determines the value of the concentration overpotential because convection due to bubbles is minimal. Above 1000 A • m -2 , microconvection determines the value of the concentration overpotential.

In the case of a concentration overpotential due to a lack of reactants, the transfer of species is limited at the electrode and the current density can no longer increase and reaches a limit value j lim . In this case, the concentration overpotential can be expressed with the following equation:

η conc = R g T n e F ln 1 - j j lim
(2.9)

When the transfer of species becomes insufficient, the current density reaches a limit value.

The flux of consumption of species at the electrode becomes equal to the flux of their transport. Thus, the current density can no longer increase. This limiting current density j lim can be expressed by the following equation:

j lim = n e F kc bulk (2.10)
1.2 Bubble effect, effect of the bubbles on current density and overpotential

The energy required for the gas to form is proportional to the voltage as shown in equation (2.1). The higher the voltage to be applied, the greater the energy to be supplied to produce the same quantity of gas. For this reason, water electrolysis improvement technologies focus on reducing the overpotentials η and the ohmic drop voltage (j×R ohm ). The bubbles produced during the electrolysis process have a direct impact on η and (j × R ohm ). Understanding their influence on these quantities allows to find a way to improve the electrolysis process.

Bubbles act as an electrical insulator

Generally speaking attached-bubbles are only produced at active sites on the electrode surface where the gas molecules are supersaturated. Generally speaking, bubbles act as an electrical insulator. Bubbles attached to the surface of the electrodes disturb current distribution and isolate active sites from reaction ion during nucleation and growth preventing other bubbles from being produced. For an effective electrolysis process, the gas released must be promptly removed from the active sites in order to provide more space for the gas release reaction. Consequently, the fast elimination of these bubbles from the electrode is crucial to increase process efficiency and allow the process to operate at higher currents density, which results in higher production rates [START_REF] Wang | The intensification technologies to water electrolysis for hydrogen production -A review[END_REF][START_REF] Zeng | Recent progress in alkaline water electrolysis for hydrogen production and applications[END_REF].

In our study we focus only on the bubbles attached to the surface, so we will briefly describe the impact of bubbles dispersed in the bulk-electrolyte on the overall efficiency of the process. When water electrolysis occurs, the bubbles are not quickly released from the electrolytic system and coat the electrode area. The phenomenon is reported to lead to high reaction overpotential and large ohmic voltage drop, and is commonly referred to as the bubble effect [START_REF] Vogt | The bubble coverage of gas-evolving electrodes in stagnant electrolytes[END_REF][START_REF] Vogt | On the supersaturation of gas in the concentration boundary layer of gas evolving electrodes[END_REF][START_REF] Sides | Phenomena and effects of electrolytic gas evolution[END_REF].

Influence of the bubbles on the void ratio

A first study by Sides and Tobias experimentally evaluated the effect of attached bubbles on the conductivity of the electrolyte using insulating particles in contact with the electrode [START_REF] Sides | Resistance of a planar array of spheres: Gas bubbles on an electrode[END_REF]. The experimental results obtained by using spherical insulating particles simulating the shape of the bubbles showed the relationship between the electrical conductivity and the void ratio ϵ, which represents the gas volume fraction of the electrolyte mixture :

σ(ϵ) σ(ϵ = 0) = 1 -ϵ 1 + 0.5 * ϵ (2.11)
This relationship is valid for ϵ < 0.5, and was later taken up by Vogt, and Sides et al. [Vogt, 1983a;[START_REF] Sides | Phenomena and effects of electrolytic gas evolution[END_REF]. As the proportion of gas increases, the electrical conductivity decreases. When the void ratio is higher than 0.5 this model loses its validity and another model called constriction must be used [START_REF] Sides | Resistance of a planar array of spheres: Gas bubbles on an electrode[END_REF]. As shown by equation (2.4), a decrease in electrical conductivity increases the resistance of the electrolyte closest to the electrode, decreasing the efficiency of the process.

Main factors

A rigorous theoretical description of the electrical effects of attached bubbles taking into account the nonuniform distribution of current density and gas supersaturation was led by later by [START_REF] Dukovic | The Influence of Attached Bubbles on Potential Drop and Current Distribution at Gas-Evolving Electrodes[END_REF]. Their results come from the numerical calculation of current density distributions around truncated spherical bubbles attached to an electrode. These results show that the effect of attached bubbles on cell voltage depends on:

• the number of bubbles per unit area of the electrode;

• the contact angles of the attached bubbles;

• the rate of the electrochemical reaction;

• the conductivity of the electrolyte as shown by a previous study [START_REF] Sides | Resistance of a planar array of spheres: Gas bubbles on an electrode[END_REF].

Links between overpotential, current density and bubble coverage

The overpotential is a function of the current density as shown by equation 2.5. Understanding how the current density is distributed around the bubble helps to understand how the overpotential is being impacted.

From their simulation, Dukovic et al. find that the bubbles attached to the electrodes capture the dissolved gas produced by the electrochemical reaction, thereby decreasing the concentration overpotential [START_REF] Dukovic | The Influence of Attached Bubbles on Potential Drop and Current Distribution at Gas-Evolving Electrodes[END_REF]. This decrease in concentration overpotential predicted by calculation has been observed experimentally later by means of advanced electrochemical techniques quantifying the effects of an isolated bubble generated on a defect positioned at the edge of an electrode [START_REF] Gabrielli | Potential drops due to an attached bubble on a gas-evolving electrode[END_REF]].

Specifically, Dukovic and Tobias find that when the concentration overpotential is low and the rate of electrochemical reaction is fast, the calculation indicates that the current density is relatively lower at the anchor line of the bubble at the electrode [START_REF] Dukovic | The Influence of Attached Bubbles on Potential Drop and Current Distribution at Gas-Evolving Electrodes[END_REF].

Conversely, at the contact line between the bubble and the electrode, the capture of dissolved gas by the bubble is greater and so is the current density. The depolarisation due to the decrease in concentration overpotential can outweigh the effects of increased electrolyte resistance and activation overpotential. When the electrochemical reaction rate is slow the distribution of current lines is smoothed over the electrode and the current fluctuations at the foot of the bubble are small.

According to their analysis, in order to take into account the effect of attached bubbles on the activation overpotential and the concentration overpotential, the relationship between current density and overpotential must be calculated using a real average current density j avg . They conclude that by considering only the attached bubbles, the impact of the ohmic effect is small compared to the additional overpotential required for charge transfer with part of the electrode surface masked by the bubbles, and that a good approximation of the addition of tension due to attached bubbles can only be calculated by considering only the surface coverage Θ they make.

However, it is important to remember that these conclusions were reached under the assumption that the dissolved gas obeys a steady state of diffusion in the concentration boundary layer.

The bubble coverage Θ

As pointed out in the previous paragraphs, the covering of the electrode by the bubbles is one of the main parameters to study the influence of the overpotential on the overall efficiency.

Vogt studies this bubble coverage in more detail [START_REF] Eigeldinger | The bubble coverage of gas-evolving electrodes in a flowing electrolyte[END_REF][START_REF] Vogt | The bubble coverage of gas-evolving electrodes in stagnant electrolytes[END_REF].

Definition

The bubble coverage of an electrode is the area of the electrode masked by all attached bubbles. Θ = S bubble S elec (2.12)

where S bubble is the area of the electrode masked by an attached bubble. A commonly used definition for the fractionnal bubble coverage Θ is the fraction of the electrode area shaded by normal projection of the bubbles on the electrode surface, Figure 1 . This is a reasonable approximation of the surface actually in contact with the bubbles because the region under an adherent bubble does not participate in the electrode's reaction [START_REF] Eigeldinger | The bubble coverage of gas-evolving electrodes in a flowing electrolyte[END_REF].

Theoretical relationship between bubble coverage and overpotential

The average current density j avg can be estimated as the ratio of the intensity of the current I applied to the electrodes to their surface exposed to the electrolyte S : j avg = I S elec . By using Tafel's equation we obtain [START_REF] Wang | The intensification technologies to water electrolysis for hydrogen production -A review[END_REF]:

η = a + b log I S (2.13)
Since the fraction covered on the electrode surface is electrochemically inactive, the actual current density j cov of the electrodes is higher than the nominal current density j avg . It can be estimated as a function of bubble coverage:

j cov = I S(1 -Θ) (2.14)
A large bubble coverage reduces the effectiveness of the active area of the electrode. As a result, j cov increases and η is also higher according to Equation 2.13 .

η = a + b log I S(1 -Θ) = η + b log j avg 1 -Θ (2.15)
As an example, Wang reports that in a typical device, the bubble coverage generates an overpotential of 0.4 V ( current density of 300 mA.cm -2 ), and inexorably raises the energy requirements for water electrolysis [START_REF] Wang | The intensification technologies to water electrolysis for hydrogen production -A review[END_REF] . One of the usual ways to reduce bubble coverage is to have a forced flow of electrolyte. Zhang and Zeng on a vertical electrode studied the detachment of bubbles in a stagnant electrolyte and then with flow [START_REF] Zhang | Evaluating the Behavior of Electrolytic Gas Bubbles and Their Effect on the Cell Voltage in Alkaline Water Electrolysis[END_REF]. Their results showed that for a flow with a Reynolds number of 2500 the diameter of the bubbles was only slightly reduced.

Estimate

Vogt and Balzer give a relation allowing to obtain a first approximation of the bubble coverage [START_REF] Vogt | The bubble coverage of gas-evolving electrodes in stagnant electrolytes[END_REF]:

Θ = 2.310 -2 (j avg ) 0.3 (2.16)
This relationship in accordance with their observation is able to correlate the data points to a good approximation.

Parameters influencing bubble coverage

Since the bubbles grow on the electrode surface over a certain period of time, the screened area is especially time dependent [START_REF] Vogt | Superposition of microconvective and macroconvective mass transfer at gas-evolving electrodes-a theoretical attempt[END_REF]:

Θ = n S tr 0 π R 2 r t r dt (2.17)
where n is the number of bubbles in contact with the electrode, t r is the average residence time of the bubble on the electrode, and R r is the bubble radius at t = t r . Vogt relates the bubble coverage to the rate of gas produced, the average residence time of the bubbles and the average volume of the bubbles when they leave. He demonstrates the interdependence of these quantities. He establishes a semi-empirical relationship between the ratio of a nominal current density I/S to the summit value of this nominal current density (I/S) su with the bubble coverage [START_REF] Vogt | Superposition of microconvective and macroconvective mass transfer at gas-evolving electrodes-a theoretical attempt[END_REF]:

I/S -j 0 (I/S) su = 3.08Θ 1.5 (1 -Θ) 1.5 (2.18)
where j 0 is an exchange current density only effective at small values of Θ only. According to Vogt, the summit value of this nominal current density (I/S) su is mainly determined by the surface condition and wettability of the electrode. Large values of (I/S) su are characteristic of strong wettability and smooth surfaces.

The thermodynamic conditions (temperature, pressure) also play an important role in the departure of bubbles and thus on the bubble coverage. Temperature decreases the density of the gas, increases the partial pressure and decreases the solubility of the gas. Pressure decreases the size of the bubbles and therefore the bubble coverage and temperature increases it. In order to account for the influence of these thermodynamic parameters on bubble coverage, Vogt developed the following relationship [START_REF] Vogt | The rate of gas evolution of electrodes-I. An estimate of the efficiency of gas evolution from the supersaturation of electrolyte adjacent to a gas-evolving electrode[END_REF]:

Θ = πR 3 r 2 V r R r t r R 2 f G Φ B I/S n e R g T F (p -p s ) (2.19)
where V r is the bubble volume at t = t r , Φ B is the current efficiency, and f G is the gas evolution efficiency and can be described as the portion of dissolved gas produced at the electrode and participating in the growth of bubbles, with the remainder of the dissolved gas dispersed in the bulk of the electrolyte. The above equation requires knowledge of the temperature and pressure conditions near the electrode, but in addition to this the residence time and radius of the bubbles. However, there is little data on the residence time and radius of bubbles, so it is difficult to use the above equation. What this equation highlights is that in order to be able to describe bubble coverage it is necessary to have a better knowledge of the nucleation, growth and departure of the bubbles from the electrode.

The gas evolution efficiency f G f G a measure of electrolysis efficiency

The gas is produced at the electrode as a dissolved gas. This dissolved gas is transported to the core of the electrolyte by diffusion and convection. Some of it participates in the growth of bubbles. But the rest can leave the electrolysis cell without being captured by the bubbles, as shown in Fig. 2. For this gas to participate in the germination of new bubbles, its concentration must be sufficient. There are several competing mechanisms: absorption and transport by the bubbles of the gas phase, and transport by diffusion and convection within the electrolyte. To account for this phenomenon, a coefficient expressing the effective gas production is defined. The gas evolution efficiency f G is the ratio between Ṅb and •N [START_REF] Matsushima | Single bubble growth during water electrolysis under microgravity[END_REF]Vogt, 1984a;Vogt, 2011a;[START_REF] Vogt | Mechanisms of mass transfer of dissolved gas from a gas-evolving electrode and their effect on mass transfer coefficient and concentration overpotential[END_REF].

f G ≡ Y 0 d Ṅb Ṅ (2.20)
where Ṅb is the rate of dissolved gas in [mol • s -1 ] transported from the electrode to the bubbles attached to the electrode and participating in their formation by desorption from the liquid phase into the gas phase and Ṅe is total rate of dissolved substance produced at the electrode area and entering the electrolyte liquid. The distance Y from the electrode in the normal direction, is considered equal to δ 1 corresponding to the thickness of the diffusion layer when the bubble average diameter d exceed this diffusion layer thickness, and Y = d otherwise. The diffusion layer thickness defined in the vicinity of an electrode is the distance over which the concentrations are different from their value in the bulk solution.

The definition of the thickness of this diffusion layer is arbitrary because the concentration approaches asymptotically the value in the bulk solution

Bubbles population

A material balance allows us to evaluate the proportion of gas that passes from the electrode to the core of the electrolyte in the form of dissolved gas. The amount of bubbles at the electrode is necessary for the calculation of f G and the effects of the attached bubble layer on • the average bubble coverage Θ;

• the volume of gas contained in the attached bubbles per unit area ;

• the average equivalent radius of the attached bubbles ;

• the average number of bubbles per unit area;

Nucleation step is partly responsible for the number and position of bubbles attached to the electrodes and the value of the dissolved gas concentration at the electrode. The nucleation frequency f nuc , the coalescence frequency f coal , and the detachment frequency f det , allow to establish the population balance of the bubbles. This population balance allows to describe the layer of attached bubbles and the population of bubbles in the bulk. This balance is established locally on the electrode:

dN b dt = f nuc -f coal -f det (2.21)

Influencing parameters and estimate

In his calculations Vogt emphasises the interrelation between f g and Θ. He establishes several relations allowing f g to be expressed as a function of several parameters [START_REF] Vogt | Mechanisms of mass transfer of dissolved gas from a gas-evolving electrode and their effect on mass transfer coefficient and concentration overpotential[END_REF]:

• the surface of the electrode covered by the bubbles Θ;

• a Sherwood number Sh 1 simulating the effect of mass transfer to the liquid bulk by diffusion;

• Sherwood number Sh 2 simulating the effect of mass transfer to adhering bubbles;

• the concentration profile in the diffusion layer;

• the average diameter of the bubbles at the time of their departure d dep .

An estimate of the normal concentration profile at the electrode over a distance δ 1 can be obtained from f G [START_REF] Vogt | Mechanisms of mass transfer of dissolved gas from a gas-evolving electrode and their effect on mass transfer coefficient and concentration overpotential[END_REF]:

c -c bulk (c e -c bulk ) fg=0 = 1 - y δ 1 - 2 3 f G 1 - y δ 1 1.5 (2.22)
when y = 0 this relation can be simplified:

c -c bulk (c e -c bulk ) fg=0 = 1 - 2 3 f G (2.23)
which coincides with the results of a previous work [START_REF] Vogt | Mechanisms of mass transfer of dissolved gas from a gas-evolving electrode and their effect on mass transfer coefficient and concentration overpotential[END_REF]. Three models are established. The first model assumes a linear concentration profile. For big bubbles Sh 1 ≡

d dep δ 1 ≥ 1.5: f G = 1 3 Sh 2 1 ΘSh 2 + 2 3 -1
(2.24)

For small bubbles Sh 1 ≡

d dep δ 1 ≤ 1.5: f G =   1/4 Sh 1 ΘSh 2 1 -Sh1 3 + 2 3   -1
(2.25)

The second model assumes a more complex concentration profile taking into account equation (2.22). For big bubbles Sh 1 ≡

d dep δ 1 ≥ 1.5: f G = 1 3 Sh 2 1 ΘSh 2 + 0.8 -1
(2.26)

For small bubbles Sh 1 ≡

d dep δ 1 ≤ 1.5: f G = 1 -Sh1 3 1/4 Sh 1 ΘSh 2 + 2/3 -0.145Sh 1.5 1 (2.27)
In the third model the assumption of a constant Sh 2 number is replaced by a relation obtained from a sphere in an infinite solution, for big bubbles:

f G = Sh1 2 6Θ + 0.8 -1 1 + 2.036 P e Sh 1 0.8 1 - 2 3 f G 1.8
(2.28)

For small bubbles :

f G = 1 -Sh 1 1 3 + 1.53( P e Sh1 ) 0.8 (1 -2 3 f 1.8 G )1/Sh 1 [1 -(1 -2 3 Sh 1 ) 2.8 ] Sh1 8Θ + 2 3 -0.145Sh 1.5 1 (2.29)
where P e is the Péclet number which is defined to be the ratio of the flow rate of advection by the flow rate of diffusion. All three models are subject to assumptions about the diameter of the bubbles at the time of detachment, their volume, the average concentration and the concentration profile, but they give a better approximation of what was previously formulated in the work of Vogt [START_REF] Vogt | The rate of gas evolution of electrodes-I. An estimate of the efficiency of gas evolution from the supersaturation of electrolyte adjacent to a gas-evolving electrode[END_REF]. Although difficult to use as is, they highlight the need to know the bubble coverage Θ and the bubble detachment diameter at the moment of the departure d dep . Of the three models, a group of adimensional numbers influencing f G was distinguished Sh2 1 ΘSh 2 -1 . The Sh1 Sh2 ratio accounts for the competition of the transport phenomena. This ratio is itself dependent on the concentration profile. Vogt's work thus makes it possible to refine the overall view of the phenomena taking place in the boundary layer in the presence of bubbles. But it also highlights the fact that a better knowledge of the microfluidic phenomena around the bubble is necessary, hence motivating the work in this thesis.

Bubble development

To sum up, the formation of hydrogen and oxygen bubbles affects electrical efficiency in three different ways:

• the dispersed bubbles in the bulk increase the bulk-electrolyte resistivity, based on the same principle, the attached bubbles also increase the resistivity of the gas-electrolyte mixture near the electrode ;

• the attached bubbles to the electrode area raise the overpotential by insulating parts of the electrode surface and crowding the current into the remaining area;

• the attached bubbles modify the overpotential concentration by affecting the level of gas supersaturation near the electrode (capture of dissolved gases).

It is therefore advisable to speed up the detachment of the bubbles from an electrolytic system in order to improve its efficiency. To understand how to accelerate the detachment of bubbles from the electrode, it is necessary to study the physical phenomena that cause their nucleation, growth and detachment.

There are of course several ways to speed up electrolytic processes, for example by using different types of electrode materials or by circulating the electrolytic fluid. A non-exhaustive list of them is described in [START_REF] Zeng | Recent progress in alkaline water electrolysis for hydrogen production and applications[END_REF]]. However, this study focuses on the microfluidic aspects. This is why, it is necessaryto have a better knowledge of the main works on nucleation, growth, and detachment of electrogenerated bubbles. The goal being to know the phenomena related to their development and to have an overview of the models used to describe it.

Nucleation

General principle

In a pure and homogeneous liquid, bubbles can form when overheating. The liquid undergoes a phase change and the process is mainly controlled by the transfer of heat. In such a system, the concentration gradient plays no role. Another process for generating bubbles is the desorption of gas molecules from a liquid to a gas phase. This generation depends on the rate of supersaturation and is essentially controlled by the concentration gradient assuming that the heat gradient in the system is negligible [START_REF] Jones | Bubble nucleation from gas cavities -a review[END_REF]. A third type of bubble generation is cavitation and can be caused by reducing the external pressure below the vapour pressure of the pure liquid.

Another type of bubble generation is that driven by chemical processes such as electrolysis. Upon supersaturation of the reduced (or oxidised) ions at the liquid-electrode interface, the molecular species produced is transformed into gas in the form of bubbles. This phenomenon begins at specific sites on the electrode called "nucleation sites". These sites are often imperfections, scratches, and grooves on the electrode surface. After the bubble detaches, some of the gas remains trapped in these imperfections, which encourages the formation of a new bubble. As the current density increases, the sites become more and more active, allowing for an increased gas flow to the electrode. Low current density implies low supersaturation, which in turn implies a low number of active sites. In addition to the current density, the temperature and concentration gradients can be locally significant and control the efficiency of the process. The heat in the system is due to the joule effect and the presence of gas molecules is due to the chemical reaction. It should be noted as pointed out by Wang that the theoretical electrical voltage U theory is a function of temperature and may decrease by increasing the electrolytic temperature [START_REF] Wang | The intensification technologies to water electrolysis for hydrogen production -A review[END_REF]. So temperature can also have an influence on the production of gas molecules. To account for the nucleation process in the case of electrolysis, the three intensive variables temperature, concentration, and pressure are to be considered. Although in most electrochemical applications the overall system pressure can be considered equivalent to atmospheric pressure, at the microscopic scale the pressure variation is to be taken into account. As reported by Tawfik and Diez if the growth or detachment of bubbles focus attention in electrolysis publications, there is a lack of information on nucleation and the prediction of their appearance [START_REF] Tawfik | On the relation between onset of bubble nucleation and gas supersaturation concentration[END_REF].

Supersaturation

Lubetkin and Blackwell defined a saturation ratio ζ r and the supersaturation degree ζ d as [START_REF] Lubetkin | The nucleation of bubbles in supersaturated solutions[END_REF] :

ζ r = c e c b ζ d = α -1 (2.30)
where c b is the bulk concentration.

There is no direct experimental means of measuring the supersaturation degree near the electrode surface. Temperature and interfacial tension are controlled parameters that depend on the electrolysis conditions and the nature of the electrolyte. On the other hand, the rate of supersaturation is linked to numerous variables such as the wettability of the electrodes, the solubility of the gas, the transport of matter and the current density. The rate of supersaturation reached at the surface of the electrolyte is therefore a difficult parameter to predict.

However, Lubetkin distinguishes two approaches for evaluating it [START_REF] Lubetkin | The motion of electrolytic gas bubbles near electrodes[END_REF]. One is theoretical and based on nucleation theory [START_REF] Dapkus | Nucleation of electrolytically evolved hydrogen at an ideally smooth electrode[END_REF][START_REF] Lubetkin | The nucleation of bubbles in supersaturated solutions[END_REF]]. The other is experimental. From the kinetic parameters obtained experimentally, the supersaturation degree is deduced [START_REF] Westerheide | Isothermal growth of hydrogen bubbles during electrolysis[END_REF][START_REF] Shibata | The Concentration of Molecular Hydrogen on the Platinum Cathode[END_REF].

Production and transport of dissolved gas

To know the degree of supersaturation, one needs to know the amount of dissolved gas produced at the electrode and how it is transported. The number of moles of dissolved gas generated at the electrode n can be estimated with Faraday's law and is related to the electrical current I(t) flowing through the electrodes:

n = I(t)dt z F (2.31)
Where z is the valency number of ions of the substance, and F = 9.648533 10 4 [C • mol -1 ] is the Faraday constant. The dissolved gas produced is then partially transported to the nucleation site. A part of these molecules is transported into the bulk, another part contributes to the growth of other bubbles, and the rest is involved in the nucleation process. Vogt reports on the various mode of mass transfers involved [Vogt, 1990a]:

• micro-convection induced by gas bubbles forming and periodically detaching from the electrode surface;

• forced or natural macro-convection of the electrolyte superimposed on the previous micro-convection phenomenon;

• diffusion which is governed by Fick's second law.

These modes act with an intensity that depends on the operating conditions. Vogt then obtains a general expression of the concentration of the dissolved gas and makes it possible to obtain an approximation in different regions of the electrolytic device [Vogt, 1990a]. However, it does not allow to deduce a profile of the concentration gradient close to the wall.

Methodology based on Fick's second law

Tawfik and Diez deduced a profile of the concentration gradient by mesuring the nucleation time [START_REF] Tawfik | On the relation between onset of bubble nucleation and gas supersaturation concentration[END_REF]. Assuming a mass transfer without convection in a stagnant electrolyte near the electrode surface and based on Fick's second law Eq. (2.32), the authors established a dissolved dihydrogen concentration profile as a function of time and wall distance.

∂c(x, t) ∂t = D ∂ 2 c(x, t) ∂x 2 (2.32)
Where D is the diffusion coefficient. The distance over which the concentration is greater than c bulk scale with δ ∼ √ Dt. Over time, c e increases and so does the thickness of this diffusion layer. By combining the analytical solution of Eq. (2.32) with experimental measurements of the time of appearance of the first bubbles, the evolution and concentration profile in the diffusion layer can be estimated . In the case of flowing electrolyte when convection is taken into account, there is no obvious method to determine the concentration gradient in the diffusion layer.

Energy state

When the dissolved gas molecules oversaturate the solution, the environment becomes favourable to the formation of bubbles. Traditionally, it is assumed that instant nucleation occurs when this environment reaches a sufficient energetic state. Ward et al. describes this energy state [START_REF] Ward | On the Thermodynamics of Nucleation in Weak Gas-Liquid Solutions[END_REF]. Assuming that a bubble is spherical and as shown on Fig. 3, the work required for a bubble to form corresponds to a critical radius R c for which a bubble could exist in thermodynamic equilibrium, although unstable, with the surrounding liquid. Using Laplace's theorem Eq. (2.33) and the laws of thermodynamics he deduced the reversible work W Eq. (2.34) necessary to generate a bubble with a critical radius R c . The necessary work must oppose the surface energy, i.e. the pressure created by the surface tension of the future bubble.

p b -p e = 2 γ R c (2.33) W = 4π γ R 2 c 3 (2.34)
where γ is the surface tension of the liquid-gas interface and p b -p e is the pressure difference between the pressure inside the bubble and pressure in the electrolyte. This work varies according to the radius of the expected bubble, as illustrated in Fig. 3. [START_REF] Blander | Bubble nucleation in liquids[END_REF] showing the work required to form a spherical bubble according to its radius. The system seeking to achieve a minimum energy state below a radius R c the bubble cannot form, above R c the bubble can grow.

In their review Blander and Katz take up the basics of classical nucleation theory [START_REF] Blander | Bubble nucleation in liquids[END_REF]. The distribution of the number of bubbles n with a critical radius R c per unit volume containing x molecules is similar to a Boltzmann distribution and in relation to the reversible work necessary to form the bubble :

n = Ze -W k B T (2.35)
where Z is a pre-exponential factor to be determined , k B is the Boltzmann constant. The addition of a molecule to a bubble with a critical radius is sufficient to trigger the growth.

Henry's Law is a gas law that stipulates that the concentration of gas molecule dissolved in the liquid is proportional to its partial pressure k H = p c , is called the Henry's constant. It links the degree of supersaturation of a gas to its partial pressure.

p b -p e = k H (c b -c e ) = p e ( c b c e -1) = p e ζ d (2.36)
From Eq. (2.36), it comes that the higher the degree of supersaturation the higher the pressure. From Eq. (2.33) and Eq. (2.34), we conclude that the higher the pressure, the smaller the critical radius of the bubbles, and therefore from Eq. (2.35), the greater the number of bubbles per unit volume [START_REF] Jones | Bubble nucleation from gas cavities -a review[END_REF].

Rate of nucleation

From Eq. (2.34), and Eq. (2.35), the rate of formation of bubbles per unit time per unit volume can be estimated. Ward et al. give a general expression of this rate [START_REF] Ward | On the Thermodynamics of Nucleation in Weak Gas-Liquid Solutions[END_REF].

J = Z exp -4πγR 2 c 3kT (2.37)
For a heterogeneous nucleation, Lubetkin gives an expression taking into account the degree of supersaturation σ [START_REF] Lubetkin | Why Is It Much Easier To Nucleate Gas Bubbles than Theory Predicts?[END_REF].

J = Z exp -16πγ 3 Φ(θ) 3 k T (σP ′ ) 2 (2.38)
where Φ(θ) is a function dependent on the contact angle and defined according to the geometry of the surface (here a flat surface). P ′ is the pressure at which the nucleation takes place.

The pre-exponential factor Z depends on the operating conditions. According to Lubetkin,the nucleation rate J is a thousand times less sensitive to a variation in the pre-exponential factor Z than to a variation in the exponential factor. So its importance is minor compared to changes in the degree of supersaturation or surface tension. The classical nucleation theory allows to access the degree of supersaturation from the nucleation rate, the measurement of the contact angle, the pressure and the temperature.

Flaws in the classical nucleation theory

Usually nucleation is said to be homogeneous when bubbles form in the liquid and heterogeneous when they develop in contact with a surface. If the classical nucleation theory makes it possible to describe the appearance of bubbles for homogeneous nucleation, a gap between theory and experience has been reported when it comes to heterogeneous nucleation on a surface.

Dapkus observing electrogenerated hydrogen gas bubbles in aqueous sulfuric acid solutions at a mercury pool electrode found that they appeared at much lower concentrations than predicted by classical theory [START_REF] Dapkus | Nucleation of electrolytically evolved hydrogen at an ideally smooth electrode[END_REF]]. Nucleation was observed on the surface of the electrode at degrees of supersaturation at least two orders of magnitude lower than the values predicted by the theory. Other authors report this same discrepancy [START_REF] Lubetkin | Why Is It Much Easier To Nucleate Gas Bubbles than Theory Predicts?[END_REF][START_REF] Hemmingsen | Spontaneous formation of bubbles in gas-supersaturated water[END_REF][START_REF] Chen | Electrochemical Measurements of Single H2 Nanobubble Nucleation and Stability at Pt Nanoelectrodes[END_REF].

Jones et al. points out that the thermodynamic approach to describe systems that are not in equilibrium is not ideal [START_REF] Jones | Bubble nucleation from gas cavities -a review[END_REF]. At very small scales where nucleation occurs, the concept of surface tension is no longer valid, this is reported as the capillary approximation. The liquid-gas interface is not so sharp. Considering the Laplace equation, enormous pressure differences are necessary between the inside of an embryo of bubble of a few nanometres radius and the liquid. Finkelstein and Tamir have developed a reliable method to determine this pressure difference. For different gases they found pressure differences between 13 and 42MPa [START_REF] Finkelstein | Formation of gas bubbles in supersaturated solutions of gases in water[END_REF]. According to classical theory, this value is 142MPa, at least three times higher than that obtained experimentally. Moreover, the classical theory does not take into account the different types of gases.

Surfaces acting as catalyst

In his review of nucleation from gas cavities, Jones points out that solid surfaces act as catalysts [START_REF] Jones | Bubble nucleation from gas cavities -a review[END_REF]. According to the classical theory in the absence of surfaces or gaz cavities, the energy required to create a bubble must be higher than that caused by the surface tension of the liquid-gas interface that opposes it. In the presence of a solid substrate, the surface area of this interface decreases. Less interfacial free energy is needed for the bubble to grow to the critical size. Nucleation is facilitated. Lubetkin and Wilt take into account this action of the solid surface by introducing a function established according to its geometry [START_REF] Wilt | Nucleation rates and bubble stability in water-carbon dioxide solutions[END_REF][START_REF] Lubetkin | Why Is It Much Easier To Nucleate Gas Bubbles than Theory Predicts?[END_REF]]. In Eq. (2.38) the Φ(θ) function allows to reduce the nucleation rate accordingly. In addition, Jones distinguishes between bubbles nucleating from nothing and those developing from an existing gas cavity. When the bubbles come off, part of the bubbles remain attached to the solid surface. For bubbles developing on these gas cavities the required energy is even lower. Therefore, the degree of supersaturation required is lower. From these reflections Lubetkin hypothesizes that on electrodes containing asperities with adequate geometry and allowing to retain these gaseous cavities the nucleation would be greatly facilitated and would allow to improve the global electrolysis process [START_REF] Lubetkin | The motion of electrolytic gas bubbles near electrodes[END_REF] .

Lubetkin suggests that the conditions required to make a flat surface more favorable to nucleation are those that, on the contrary, make it more difficult to detach the bubbles once formed [START_REF] Lubetkin | The motion of electrolytic gas bubbles near electrodes[END_REF]. Sakuma experiments under microgravity and demonstrates that bubbles at the time of their detachment are smaller on a hydrophilic surface than on a hydrophobic surface [START_REF] Sakuma | Nucleation and growth of electrolytic gas bubbles under microgravity[END_REF]. On hydrophilic surface the contact angle approaches zero facilitating the detachment of bubbles. However, in this case the nucleation approaches the conditions of a homogeneous nucleation and makes it more difficult. On the other hand, on hydrophobic surfaces, the bubble nucleates and spreads more easily, the hydrophobic support not being favourable to the electrolyte, the gas covers this surface more easily. It becomes larger but takes longer to detach as a result. While the bubble remains attached, it successively blocks neighbouring nucleation sites as it develops.

Growth

Bubble growth in supersaturated solutions has been studied from a theoretical, experimental, and numerical point of view. Among the papers adopting the theoretical approach one can find that of Scriven and Epstein-Plesset [START_REF] Epstein | On the Stability of Gas Bubbles in Liquid-Gas Solutions[END_REF][START_REF] Scriven | On the dynamics of phase growth[END_REF]. In their work, Epstein and plesset studied the evolution of the radius of a bubble at different levels of supersaturation, assuming a spherical bubble in an infinite solution and neglecting the possible effects of convection surrounding the bubble. Scriven details this work taking into account all assumptions and limitations. Exact solutions of the equations are obtained for typical bubble growth conditions. An asymptotic expression of the evolution of the bubble radius is obtained. These early studies established the first relationships to describe bubble growth.

Mass balance at the interface

By doing a mass balance in the gas phase and assuming that the bubble is spherical by deriving, the growth rate of the bubble ṁB [kg • s -1 ] can be obtained:

ṁB = ρ g d dt 4 3 πR(t) 3 = ρ g 4πR(t) 2 dR(t) dt (2.39)
Calculating the same quantity but considering the liquid side this time and assuming transport by diffusion:

ṁB = M g 4πR 2 D dc dr R (2.40)
where r is the radial coordinate in a spherical coordinate system in which the origin coincides with the centre of the bubble.

In the diffusive regime, the concentration profile at each time and position is defined by the radially symmetric diffusion equation:

∂c(r, t) ∂t = D∇ 2 c(r, t) = D r 2 ∂ ∂r r 2 ∂c(r, t) ∂r (2.41)
Considering a single bubble in a infite domain, the boundary conditions for solving the equation (2.41) are defined as:

c(r, 0) = c 0 , r >> R (2.42a) lim r-→∞ c(r, t) = c 0 , t > 0 (2.42b) c(R, t) = c s , t > 0 (2.42c)
These boundary conditions stablish that, at the beginning of the bubble growth, the solution has an homogeneous concentration c 0 = Hp 0 which will depend on the temperature T 0 and the pressure p 0 at which the saturated solution was prepared Eq. (2.42a); very far away from the bubble, the concentration remains unaltered at all times Eq. (2.42b); and finally, the concentration at the bubble boundary remains constant at a value c s = Hp s which depends also on temperature T s (which usually coincides with T 0 ) and pressure p s at which experiments are performed Eq. (2.42c). Solving Eq. (2.41) using these boundary conditions results in:

∂c ∂r R = (c 0 -c s ) 1 R + 1 √ πDt (2.43)
This brings us to the Epstein-Plesset equation which defines the change in the growth of the bubble with respect to time:

dR dt = M D(c 0 -c s ) ρ g 1 R + 1 √ πDt = M W (c 0 -c s ) ρ l D l πt 1 + √ πD l t R (2.44)
if c s > c 0 the solution is undersaturated, the gas flow goes from the bubble to the bulk. If c s = c 0 equilibrium, if c s < c 0 the solution is supersaturated.

Radius of the bubble

In experimental approaches, usually to account for the growth of the bubble a relation of its radius as a function of time is used.

R(t) = β × t b (2.45)
where β is the growth coeficient principally dependent on the current density according to Brandon and Kelsall and b is the time coefficient [START_REF] Brandon | Growth kinetics of bubbles electrogenerated at microelectrodes[END_REF]]. These two parameters vary according to the studies and the different phases of the growth. Several authors use this model as a basis for presenting their experimental results [ • the chemical reaction which is determined by the electrolyte and surface properties;

• the gas molecule transfer process, which is determined, among other things, in a stagnant electrolyte by the diffusion of the gas molecule and the radial micro-convection effect due to the expansion of the bubble interface;

• the desorption and the absorption of the gas molecules on the bubble interface.

From these steps, Wang et al. assume that the competition between desorption and adsorption at the liquid-bubble interface is not the process limiting bubble growth but rather the combination between mass transfer and chemical reaction rate [START_REF] Wang | Investigations on bubble growth mechanism during photoelectrochemical and electrochemical conversions[END_REF]. Its importance in explaining the dynamics of growth would therefore be minor. Another phenomenon to be taken into account is the coalescence of small bubbles to large ones. According to Vogt the validity of Equation 2.45 is restricted to case where bubbles do not interfere with each other [Vogt, 1983a]. The parameter b takes different values depending on which process is limiting.

Inertial growth

During the initial tens of microseconds b = 1 the growth is inertia-controlled [START_REF] Brandon | Growth kinetics of bubbles electrogenerated at microelectrodes[END_REF]]. The growth is described by Rayleigh equation :

R(t) = 2∆p 3ρ 0.5 t (2.46)
where ρ is the electrolyte density. The bubble growth is driven by the high difference of pressure ∆p determined by the Laplace equation and due to the interfacial tension force.

Growth by diffusion

When the movement of gas molecules is slow and controlled by diffusion phenomena, Scriven showed that the radius increase with the square root of time i.e. b = 0.5 and the growth coefficient is a function of the dissolved gas concentration and the diffusion coefficient [START_REF] Scriven | On the dynamics of phase growth[END_REF].

R(t) = 2β(D t) 0.5 (2.47)
where β can be calculated as follows:

β = ∆c 2πρ 1 + 2πρ gaz ∆c 0.5 0.5 (2.48)
where ∆c is the supersaturation of the species under consideration near the electrode.

Wang extends his theory to the phenomenon of micro-convection due to the spreading of the bubble using the surface renewal theory to simulate the gas-liquid interfacial mass transfer [START_REF] Wang | Investigations on bubble growth mechanism during photoelectrochemical and electrochemical conversions[END_REF]. To simulate the mass transfer around the bubble they add the diffusive flow to that caused by the expansion of the bubble. The time coefficient remains unchanged, but the he finds that the effect due to micro-convection on growth may be 100 times higher than that of diffusion.

Growth limited by the reaction rate

When gas molecules move rapidly in the liquid or when the gas-liquid interface is close to the electrode, the process limiting the growth is the reaction rate. Some of the generated molecules are dispersed in the electrolyte while the other contributes to the generation of bubbles N G [mol.s -1 ]. Assuming that the bubble is spherical, we obtain:

N G = d dt 4πR 3 3V m (2.49)
where V m is the molar volume of gas. Wang deducts that b = 1/3 when the bubble growth is controlled by the surface chemical reaction [START_REF] Wang | Investigations on bubble growth mechanism during photoelectrochemical and electrochemical conversions[END_REF].

R(t) = 3V m N G 4π 1/3 t 1/3 (2.50)
It is almost impossible to study the individual behaviour of bubbles on conventional electrodes due to the high bubble coverage and void fraction at practical current densities [START_REF] Massing | Thermocapillary convection during hydrogen evolution at microelectrodes[END_REF]. This is why most experimental studies on the subject are done on microelectrodes.

These electrodes are about 100 micrometers in diameter and allow the observation of the development of a single bubble. This facilitates the study of its nucleation, growth and detachment. However, one of the main differences between a microelectrode and a conventional electrode is that the bubble completely covers its surface before detachment.

The bigger the bubble gets, the more it blocks the electrode. So there is a high current densities and high local supersaturation level at the bubble foot. Whereas for a conventional electrode part of the gas escapes into the liquid, for a microelectrode almost the entire hydrogen generated diffuses directly into the bubble at its foot. Therefore the behaviour of a bubble on a micro electrode and a conventional electrode will not be the same, and the observed growth relationships can be expected to be different.

Detachment

Break-off diameter

An appropriate analysis of the detachment is essential to determine the bubble coverage and therefore the overpotential and the efficiency of the electrolysis process. Usually a balance of forces acting on the bubble is used to determine the diameter that the bubble will have at the time of detachment and the residence time of the bubble on the electrode surface. These two parameters are usually used to account for the detachment. They are determined by the nucleation and growth steps.

Lubetkin notes that nucleation kinetics and detachment kinetics are closely related and that the classical theory of heterogeneous nucleation is only valid when the detachement is fast compared to the nucleation rate [START_REF] Lubetkin | The nucleation and detachment of bubbles[END_REF] . Vogt et al. points out that the analogy between boiling and the evolution of electrolysis bubbles is limited and that the detachement diameters of electrolysis and boiling bubbles differ greatly [START_REF] Vogt | The limits of the analogy between boiling and gas evolution at electrodes[END_REF]. In most of the relationships established by Vogt and Balzer to determine the mass transfer coefficients, a term is related to the break-off diameter of the bubble [START_REF] Vogt | The bubble coverage of gas-evolving electrodes in stagnant electrolytes[END_REF]. They use a force analysis to correlate the bubble coverage with the electrolyte flow.

Study of detachment on microelectrodes

Under normal operating conditions the bubble coverage near the electrodes is very high, which limits observation and makes it difficult to analyze bubble growth and detachment. In addition, bubbles that detach and rise near the surface create turbulence and prevent the determination of forces on an isolated bubble. This is why many electrolysis experiments [Fernández et In this case the bubble forms on a small surface. Usually before the detachment takes place, the bubble covers the entire surface of the micro or nanoelectrodes forcing the current to pass through the corner formed by the base of the bubble and the contact surface, causing an increase in current density. Its distribution is thus different compared to larger electrodes.

Force balance

Depending on the experimental conditions, stagnant electrolyte or electrolyte with flow, vertical or horizontal electrodes, micro or macro electrode, the intensity and direction of the forces involved in the balance may change. However, the principle remains the same. Their point of application is assumed to be in the center of the bubble. The bubble is supposed to be spherical. For a horizontal electrode, detachment is supposed to occur when the sum of the force projections on the vertical axis is equal to zero. The force balance usually include:

• surface tension force holding the bubble on the electrode;

• the buoyancy, the inertial, and the contact pressure force pushing the bubble away from the surface;

• the hydrodynamics forces acting in either way depending on the experimental conditions [START_REF] Hibiki | Lift force in bubbly flow systems[END_REF].

Inertial force

The inertial force is due to the growth of the bubble:

F i = 4πR 3 3 ρ g d 2 R dt 2 (2.51)
Usually neglected because the mass of the bubble is very small and the growth rate of the bubble is also very small.

Buoyancy force

The buoyancy force is the resultant of the difference between Archimedes' thrust and the force of gravity, it is given by:

F b = 4πR 3 3 (ρ g -ρ l )g (2.52)

Contact pressure force

The contact pressure force is is due to the pressure difference between the inside and outside of the bubble (see Laplace equation) and pushing on the solid surface S CP of the electrode in contact with the bubble.

F CP = S CP (p L -p G )ndS ≈ πR 2 c 2γ R (2.53)
where R c is the radius of the contact surface.

Surface tension force

The surface tension force depends on the angle θ of contact of the bubble with the surface.

F S = 2πR c γsin(θ) (2.54)
In a stagnant electrolyte on a horizontal microelectrode Chen et al. note that the surface tension and contact pressure forces are the main forces that influence the balance followed by the buoyancy force, the inertial and hydrodynamic forces being negligible [START_REF] Chen | Dynamics of single bubble departure from TiO2 nanorod-array photoelectrode[END_REF]. This means that the surface tension of the liquid-gas interface is one of the most sensitive parameters in determining the bubble departure.

A missing force ?

With only these forces taken into account, the force balance usually failed to predict the bubble departure diameter. Lubetkin reports different bubbles phenomena observed during electrolysis that cannot be explained by the usual force balance [START_REF] Lubetkin | The motion of electrolytic gas bubbles near electrodes[END_REF] . One of them is called the "bubble jump-off and return". When 2 bubbles coalesce, the new bubble comes off the surface. The surface tension force is then zero and the buoyancy force is supposed to carry the bubble into the bulk. However, the bubble immediately returns to the electrode surface. Lubetkin deduced from this the evidence of a force created by a Marangoni effect. The surface energy of a liquid depends on its composition and temperature, among other things. Thus, variations in composition or temperature along the surface of a liquid imply variations in surface tension. Areas of higher tension pull harder, resulting in tangential stress at the interface, called Marangoni stress. This stress introduces a tangential stress jump at the interface. This jump is compensated by the viscous stresses of the two fluids.The two fluids on either side of the interface are therefore carried away by viscosity, generating Marangoni flows [START_REF] Scriven | The marangoni effects[END_REF].

3 Transport of the diluted species in the vicinity of the bubble

Mass transfer Coefficient

Definition

The mass transfer coefficient k is a term that illustrates this effect of agitation in the calculations of the interfacial mass transfer or the mass transfer in the bulk . It has the dimensions of a velocity:

k = J C (2.55) with J the material flow density in [mol • m -2 • s -1 ] , and c the concentration[mol • m -3 ].

Diffusion induced transport coefficient

In the diffusion layer the mass transfer coefficient can be expressed as :

k dif f = D δ (2.56)
with D the diffusion coefficient of the species m 2 • s -1 , and δ in [m] the thickness of the diffusion layer. In the absence of gas bubble production, the material transport coefficient depends solely on the nature of the flow, the characteristics of the fluid (kinematic viscosity) and the geometry of the cell.

Expression in terms of the Sherwood number

This coefficient has been measured or calculated for many geometries and flow types. It is expressed by means of dimensionless number correlations [START_REF] Reid | Chemical engineers[END_REF][START_REF] Fahidy | Elements de genie electrochimique[END_REF].

These correlations are established between the Sherwood, Reynolds (or Grashof), Schmidt numbers and dimensional ratios characteristic of the flow.

Sh = A Re a Sc b L d c (2.57)
with A , a, b, c constants,and L/d ratio of characteristic flow dimensions. The Sherwood number for transport at the electrode is proportional to the transport coefficient:

Sh = kδ D (2.58)
The Schmidt number is defined by :

Sc = µ ρD (2.59)
The Reynolds and Schmidt numbers are used to characterise the thickness of the electrolyte layer near the electrodes in which transport takes place mainly due to diffusion. Generally, the greater the forced convection or agitation, the higher the Reynolds number, the thinner the layer and the faster the transport.

Bubble induced transport coefficient

When an electrode generates gas, the bubbles cause agitation of the liquid in the vicinity of the electrode. This bubble-induced convection strongly increases the transport coefficient [Vogt, 1983b;[START_REF] Ward | Bubble evolution in solutions with gas concentrations near the saturation value[END_REF]. The transport coefficient to describe the bubble-induced transport is denoted k b . Experimental data show that k b increases with the amount of gas produced [Vogt, 1983b]. There is therefore a relationship between the bubble-induced transport coefficient k b and the current density k b = p 1 j p 2 . A very large number of experimental relationships are available as a function of the gas produced, the electrode material, the temperature, the pressure and the electrolyte, these relationships are collected in [START_REF] Janssen | Behaviour of and mass transfer at gas-evolving electrodes[END_REF]Stephan and Vogt, 1979a]. It is shown that the transport coefficient, at a given current density, decreases as the size of the bubbles at detachment increases.

Superposition of transport coefficients

It is common practice to superimpose bubble-induced transport coefficient and diffusion-induced transport coefficient to describe electrolyte flow. In their study Bisang et al. validate this practice experimentally by showing the effects of this superposition on the current density distribution at the surface of an electrode [START_REF] Bisang | Effect of mass transfer on the current distribution in monopolar and bipolar electrochemical reactors with a gas-evolving electrode[END_REF]. The equations Eq. (2.60) and Eq. (2.61) are usually used to calculate the global transport coefficient k t from the independently calculated coefficients k dif f and k b [Vogt, 1983a].

k t = k dif f + k b (2.60) k t = k 2 dif f + k 2 b (2.61)
The choice between one or other of these equations is a choice validated a posteriori with regard to the results obtained by the model used. Therefore, to calculate k t , a model is needed to calculate the effect of the movement of the electrolyte in the vicinity of the electrode k dif f without the presence of the bubbles and a model to describe the effect of the bubbles produced on k b .

The calculation of k t allows us to evaluate the movement of the electrolyte in the area around the bubble and the electrode. It allows us to estimate how the dissolved species are transported and to evaluate the concentration overpotential from equation (2.8).

Models describing the effect of bubbles on the transport coefficient

Models estimating the average transport in the vicinity of bubbles

As shown in the diagram in Fig. 4 the amount of dissolved gas involved in the growth of the bubble will depend on three mechanisms: the production of dissolved gas at the electrode, the transport of it from the electrode to the interface, and its absorption at the interface. The transport of dissolved species is a key step, yet there is no consensus on how to model it. Different types of competing models have been developed to describe the effect of bubbles on the transport coefficient near an electrode. Each model assumes that the phenomenon it models is sufficient to describe the transport of species. Some authors do not exclude that different forms of agitation of the diffusion boundary layer can coexist and that consequently the transport coefficients calculated by the different models can add up [START_REF] Janssen | Behaviour of and mass transfer at gas-evolving electrodes[END_REF].

In the "penetration model", transport is increased by replacing the volume occupied by a bubble at the time of its detachment by the same volume of liquid from the core of the electrolyte just after detachment. The second type of model is "hydrodynamic". This model considers that the mechanism controlling the transport of material around the electrodes is the free convection brought by the detached bubbles. This model was initially developed to describe the material transport in an electrolysis cell in which gas is injected to agitate the electrolyte. Finally, the last type of model considers that the bubbles increase the transport of matter by agitating the electrolyte that surrounds them during their growth [Stephan and Vogt, 1979b;[START_REF] Vogt | Superposition of microconvective and macroconvective mass transfer at gas-evolving electrodes-a theoretical attempt[END_REF]]. This model is called "microconvective model" and predicts a transport of matter described by a correlation built from dimensionless numbers.

microconvective model

The microvonvective model assumes that the overall flow in the electrolyte can be decomposed into different microconvection phenomena. Relationships based on the Sherwood number are established theoretically and experimentally. The mass transfer coefficient is determined from this Sherwood number [Vogt, 1984b;[START_REF] Vogt | Superposition of microconvective and macroconvective mass transfer at gas-evolving electrodes-a theoretical attempt[END_REF][START_REF] Vogt | Local microprocesses at gas-evolving electrodes and their influence on mass transfer[END_REF].

A first relationship was determined by Wragg, to account for free convection on a horizontal electrode [START_REF] Wragg | Free convection mass transfer at horizontal electrodes[END_REF]. This phenomenon of microconvection taking place in the vicinity of the electrode is assumed to be caused by the difference c e which is the concentration of diluted specie adjacent to the electrode and c bulk which is the bulk concentration:

Sh = k f L D = 0.64(Gr Sc) 0.25 (2.62)
The Grashof number is :

Gr = ρβ * (c e -c bulk ) L 3 ν 2 (2.63)
where β * is:

β * = - 1 ρ ρ c T,p (2.64)
According to Vogt, the density difference in equation (2.63) is affected by the temperature difference between liquid bulk and the electrode-liquid interface in addition to the corresponding concentration differences of all subspecies [Vogt, 1993a]. The author introduces the term of single-phase free convection emphasising that the origin of this microconvection is independent of the gas phase, but not its extent, and repeats the previous equation but focuses on the temperature and concentration gradients that cause the density difference.

The dimensionless mass flux of the dissolved gas is given by:

Sh = kd D = 0.72 1 -Θ 1 -2 3 f G 0.8 j ϵα o gL 4 2F ν 3 L Sc 2 0.2 (2.65)
where the expansion coefficient α o is :

α o = 13 + 71 1 -2 3 f G × 10 -6
(2.66)

The expansion coefficient α o accounts for the difference in density caused by concentration and temperature gradients.

Another type of microconvection influence the mass transfer of the dissolved substance once the bubble begins to grow at the electrode surface is the bubble induced microconvection [Vogt, 1993a;Vogt, 1993b].

Sh = k d D = 1.89 (Re f G ) 0.5 Sc 0.487 [Θ 0.5 (1 -Θ 0 .5)] 0.5 1 -2 3 f G (2.67)
According to Matsushima et al., the above relationships were later confirmed on the growth of a single bubble under microgravity [START_REF] Matsushima | Single bubble growth during water electrolysis under microgravity[END_REF].

Mass transport within a phase depends directly on the concentration gradient of the species being transported in that phase. Mass can also be transported from one phase to another, and this process is called interphase mass transfer.

General principles

At the liquid-gas interface, the sorption phenomenon can be described as physical, whereas at the electrode it can be described as chemical. Physical absorption or non-reactive absorption is a process of mass transfer that does not involve with chemical reaction occurring in liquid phase. The rate of absorption-desorption at the interface depends on properties of gas-liquid fluid dynamics, interfacial area between phases, concentration difference, solubility of gas, temperature, pressure, and its duration time of contact.

Equilibrium at the interface

One of most important factors is the solubility of gas in liquid. It depends on temperature, pressure and the characteristics of the substance itself. There is an equilibrium between the gas and the liquid phase. This gas-liquid equilibrium can be described with Henry's Law with the assumption of an ideal liquid solution and that the perfect gas law can be applied.

c(T, P ) = p k H (T ) (2.68)
where T is the temperature , k H (T ) is the Henry's constant ( dependent on each gas-liquid couple and a decreasing function of T ) and p is the partial pressure of the gas in the liquid.

In dilute conditions, Henry's law has good capability to predict the gas-liquid equilibrium.

In any multiphase-multicomponent systems, the chemical species will always move towards equilibrium. The state of non equilibrium is what causes the mass transfer, as it is the mass transfer that drives the system towards equilibrium.

Unsteady mass balance at the interface

The basic theory of mass transfer with absorption is the two-film theory. At the interfaces of both phases, there are two films, gas film and liquid film connecting to each other. However, the interfacial mass transfer between the electrolyte and the gas bubble is an unsteady process. This non-stationary aspect is not taken into account by the film theory.

Usually to model species/mass transfer across an interface a coefficient of mass transfer k is used. The mass transfer rate through the bubble surface

ṁ in [kg • m -2 • s -1 ] is calculated as follows: ṁ = M H 2 kA I (c H 2 ,I -c H 2 ,sat ) = ρk (y H 2 ,I -y H 2 ,sat ) (2.69)
where A I is the interfacial area, c I the concentration of hydrogen around the gas bubble, c sat the saturation concentration of dissolved hydrogen, y I and c sat are the equivalent values converted in species mass fraction. In their simulation of the growth of a hydrogen bubble on the surface of an electrode, Liu et al. used this kind of relationship to model the mass transfer across the interface [START_REF] Liu | Numerical simulation of hydrogen bubble growth at an electrode surface[END_REF]. The difference (c I -c sat ) represents the driving force of the interfacial species transfer.

Estimate from Sherwood number

The species transfer coefficient k H 2 can be obtained from the Sherwood number :

Sh = 2k H2 R D H 2 (2.70)
where R is the bubble radius, and D H 2 is the diffusion coefficient of the diluted H 2 in the bulk.

The Sherwood number is used to describe the ratio of the overall species transfer to pure species transfer. Empirical or analytical relations that are used to evaluate the mass transfer through a spherical surface can be adopted to calculate the Sherwood number :

Sh = 2 + 0.6 2Rv s ρ µ 1/2 µ ρD 1/3 (2.71)
where v s is the liquid velocity around the the sphere surface [START_REF] Kashchiev | Kinetics of the initial stage of isothermal gas phase formation[END_REF], and can be evaluated as v s = dR dt ;

By combining Eq. (2.43) and Eq. (2.40) we can obtain an analytical solution of the Sherwood number [START_REF] Bird | Transport Phenomena[END_REF]:

Sh = 2 1 + R (πDt) 0.5
(2.72)

These correlations depend both on the shape of the bubble and on the operating conditions, that is why modeling mass transfer on the basis of Sh may lead to erroneous results [START_REF] Deising | Direct numerical simulation of mass transfer in bubbly flows[END_REF].

Penetration theory

Principle

The penetration theory was suggested in 1935 by Higbie who was investigating whether or not transfer resistance existed at the interface when a pure gas was absorbed into a liquid [Higbie, 1935a]. As shown in Fig. 5, Higbie assumed that each surface element of the liquid was exposed to the gas for the time required for the gas bubble to pass through it. Vortices in the fluid bring an element of the fluid to the interface where it is exposed to the second phase for a defined time interval, after which the element of fluid is mixed back into the mass of the fluid.

Thus, an element of fluid whose initial composition matches that of the main fluid away from the interface is suddenly exposed to the gas phase. An unsteady molecular diffusion process then occurs within the fluid element. Equilibrium is assumed to be reached immediately at the interface, i.e. transfer across the interface is instantaneous. The element is then displaced or remixed after a fixed time interval, i.e. each element remains in contact with the gas for the same period of time.

The existence of a velocity gradient within the liquid element is ignored and the fluid at all depths of the element is assumed to move at the same speed as the interface.

Expression of the mass transfer coefficient

Under these conditions the convection terms in the transport equation within the fluid element can be neglected and can be written as : where z is the distance from the interface in the normal direction.

∂c ∂t = D ∂c ∂z (2.73)
c(z, 0) = c 0 , z > 0, t = 0 (2.74a) c(z, t) = c i , z = 0, t > 0 (2.74b)
where c 0 is the initial concentration of the fluid element, and c i is the equilibrium concentration at the interface. Solving this differential equation in the fluid element gives a value for the mass transfer coefficient:

k(t) = D πt (2.75)
The average mass transfer coefficient during a time interval corresponding to the contact time of the fluid element with the t c interface can be obtained by integrating the previous equation:

k = 1 t c tc 0 k(t)dt = 2 D πt c (2.76)
The mass transfer coefficent is proportional to the square root of the diffusivity.

There are important differences in the implications of the theories when one must consider the impact of contaminants (surfactants) on gas transfer [START_REF] Painmanakul | Effect of surfactants on liquid-side mass transfer coefficients[END_REF].

5 Microfluidic phenomena in electrolysis

Theory -surface tension

Surface tension an expression of the bond between molecules

The molecules of condensed phases have strong attractive interactions between them. The molecules located at the interface have fewer neighbours and therefore form fewer interactions, which creates a energy deficit. In order to decrease this energy deficit, a force pulls tangentially on the interface to reduce its area. The interfaces of the condensed phases are therefore subjected to a tension, called surface tension.

We note γ the value of this tension, which is expressed in

[N • m -1 ].
It also corresponds to the energy deficit per unit area and can therefore also be expressed in

[J • m -2 ].
In a liquid, the molecules are constantly changing neighbours due to thermal agitation. The interaction energy between molecules is therefore of the order of the thermal energy k B T , where k B is the Boltzmann constant and T the temperature. The order of magnitude of the surface energy of a liquid can be estimated by relating this thermal energy to the average surface occupied by a molecule a 2 :

γ ∼ k B T a 2 ≈ 25mN • m -1
where a is the average distance between molecules, about 4 Ȧ. This calculation gives a value of 25mN • m -1 , which is the right order of magnitude for ethanol, silicone oils or alkanes. Due to their strong polarity, water molecules bind to each other by stronger interactions than van der Waals, called hydrogen bonds. The surface energy of water is higher and is about 73mN/m at 20°C.

Pressure jump

By describing the surface tension balance at the edges of a curved interface element, Young deduced that this resulted in a pressure jump on both sides of the interface [START_REF] Young | III. An essay on the cohesion of fluids[END_REF]. This pressure jump is called Laplace's pressure in honour of a contemporary scientist of Young's, Pierre-Simon de Laplace, who provided a more precise formalism:

∆p = p inside -p outside = γκ (2.77)
where κ [m -1 ] is the curvature of the interface. Depending on the sign of the curvature, the Laplace pressure jump can be positive or negative.

Marangoni stress

The curvature of a plane interface is zero, so the Laplace pressure is also zero. However, the surface tension balance is not necessarily so. The surface energy of a liquid depends on its composition and temperature, among other things. Thus, variations in composition or temperature along the surface of a liquid imply variations in surface tension. Areas of higher tension pull harder, resulting in tangential stress at the interface, called Marangoni stress. This stress introduces a tangential stress jump at the interface. This jump is compensated by the viscous stresses of the two fluids. The two fluids on either side of the interface are therefore carried away by viscosity, generating Marangoni flows [START_REF] Scriven | The marangoni effects[END_REF]. Thus, in the general case of a curved and heterogeneous interface, the surface tension balance on an elementary surface can be broken down into two terms, a normal stress, the Laplace pressure, and a tangential stress, the Marangoni stress. Due to the surface tension gradient the fluid near the interface moves from areas where the surface tension is low to areas where it is higher. In electrolysis surface tension can be found as a function of:

• the temperature, the marangoni effect is referred to as the thermocapilary effect;

• voltage, the marangoni effect is referred to as the electrocapilary effect;

• concentration of surfactants, the marangoni effect is referred to as the solutocapilary effect.

Marangoni motion near electrolytic gas bubbles

Observations

In his review, Lubetkin reports the existence of bubble-related phenomena that cannot be explained by classical theories [START_REF] Lubetkin | The motion of electrolytic gas bubbles near electrodes[END_REF]. He suspects the existence of a Marangoni effect at the interface between the electrolyte and the bubble. Based on this hypothesis, he provides an explanation for these phenomena which include :

• a phenomenon of radial adhesion, a larger bubble attracts other smaller bubbles;

• bubbles following oscillating trajectories along a vertical electrode [START_REF] Janssen | The effect of electrolytically evolved gas bubbles on the thickness of the diffusion layer[END_REF]];

• oscillations of bubbles attached to the electrode.

The Marangoni effect is a transient effect. Lubetkin assumes that these bubble shifts or oscillations are due to a force caused by the soluto-Marangoni effect. This force varies with the concentration of dissolved species. As the bubbles absorb the dissolved gases, the force varies. As a result, the bubbles leave a low concentration area after absorption and move towards the higher concentration areas, which creates an oscillatory movement if the bubble does not coalesce. According to his reasoning the soluto-Marangoni effect is more important than the themo-Marangoni effect created by a temperature gradient, but it has a shorter range than the thermo-capillary effect.

A more significant phenomenon was reported by Lubetkin [START_REF] Lubetkin | The motion of electrolytic gas bubbles near electrodes[END_REF]. When two bubbles of similar size are adjacent to each other on the horizontal surface of the electrode, and they meet and merge, the centre of the new bubble is higher than the position of either of the two original smaller bubbles. A normal impulse at the surface of the electrode is produced and may in some cases force the bubble to leave the electrode. What is interesting in this case is that the bubble subsequently returns to the electrode. This phenomenon was first observed by Westerheide and Westwater in 1964, who hypothesised that the return could be due to surface voltage gradients or electrostatic interactions [START_REF] Westerheide | Isothermal growth of hydrogen bubbles during electrolysis[END_REF]]. Lubetkhin's reasoning rejects the impact of drag and electrostatic forces, leaving only Marangoni forces to explain the phenomenon. Lubetkin hypothesises that this phenomenon was not suspected earlier because it suffers from interface masking phenomena due to surfactants.

Using Particle Tracking Velocimetry (PTV) Yang et al. were able to observe convective vortices around the bubble. An example of a result is shown in Fig. 6. From these results Yang et al. deduced that the convection movement of the electrolyte was caused by the Marangoni effect.

Origin of the surface tension gradient

However, it is still unclear whether this Marangoni effect is due to a variation in temperature, concentration of dissolved species, potential, or by pollution from surfactants. Lubetkin originally seemed to favour the soluto-Marangoni effect. Yang et al. assumed that the solutal and thermal effect were of the same order of magnitude. By comparing the experimental In their work, Meulenbroek et al. show that Marangoni convection in the vicinity of electrochemically generated bubbles is the result of thermo and solutocapillary effects at the bubble interface [START_REF] Meulenbroek | Competing Marangoni effects form a stagnant cap on the interface of a hydrogen bubble attached to a microelectrode[END_REF]. Unlike Lubetkin, they do not attribute this Marangoni effect to dissolved species but to other surfactants present and adsorbed on the liquid-gas interface. Meulenbroek et al. include the electrocapillary effect in their study. This appears to be of minor importance compared to the other effects.

The Marangoni force

Marangoni flow around a growing bubble at an electrode surface is known to delay detachment of bubbles. By examining the motion of a gas bubble at a distance z from a heated surface, and by solving the linearised Navier-Stokes equation, with appropriate boundary conditions, it was shown that the thermal Marangoni effect was linearly dependent on a temperature gradient [START_REF] Young | The motion of bubbles in a vertical temperature gradient[END_REF][START_REF] Hardy | The motion of bubbles in a vertical temperature gradient[END_REF][START_REF] Morick | Migration of Air Bubbles in Silicone Oil under the Action of Buoyancy and Thermocapillarity[END_REF][START_REF] Lubetkin | The motion of electrolytic gas bubbles near electrodes[END_REF]. The shear stress on the interface induced by the thermal Marangoni effect is:

τ M T = dγ dT dT dz (2.78)
The thermal Marangoni force F M T acting on the bubble attached to an electrode can be approximated by calculating the area integration of the shear stress over the entire spherical surface of the bubble, giving the following expression [START_REF] Lubetkin | The motion of electrolytic gas bubbles near electrodes[END_REF][START_REF] Chen | Dynamics of single bubble departure from TiO2 nanorod-array photoelectrode[END_REF]:

F M T = dγ dT ∆T πR (2.79)
where ∆T represents the difference in concentration between the surface of the electrode and the liquid bulk. By analogy it is possible to obtain a similar expression for the solutocapillary effect:

F M C = dγ dc ∆cπR (2.80)
Assuming that the two effects are independent of each other it is possible to obtain a general expression for this force:

F M = F M T + F M C = dγ dT ∆T πR + dγ dc ∆cπR (2.81)
Both dγ dT and dγ dc are negative, so the combined Marangoni force pushes the bubble towards the electrode surface.

Stagnant cap model

Meulenbroek et al. used a stagnant cap model that describes the transport of the contaminants along the interface [START_REF] Levich | Surface-Tension-Driven Phenomena[END_REF][START_REF] Griffith | The effect of surfactants on the terminal velocity of drops and bubbles[END_REF][START_REF] Sadhal | Stokes flow past bubbles and drops partially coated with thin films. Part 1. Stagnant cap of surfactant film -exact solution[END_REF][START_REF] He | The size of stagnant caps of bulk soluble surfactant on the interfaces of translating fluid droplets[END_REF][START_REF] Hosokawa | Evaluation of adsorption of surfactant at a moving interface of a single spherical drop[END_REF][START_REF] Shmyrov | Phase transitions on partially contaminated surface under the influence of thermocapillary flow[END_REF]. The competition of the Marangoni effects results in the formation of a stagnant cap at the top of the bubble. The surfactants cover the top of the bubble and prevent the Marangoni effect from taking place, while the bottom part of the bubble interface is mobile and drives a Marangoni flow. The interface is stiffened by the surfactants and prevents more minor surface tension variations such as that caused by a temperature or dissolved species gradient.

Effect of electrode size

The particularity of the previous studies is that they are conducted on microelectrodes. Hossain et al. performed several numerical simulations of the thermocapillary effect at the level of electrogenerated bubbles on electrodes of varying size [START_REF] Hosokawa | Evaluation of adsorption of surfactant at a moving interface of a single spherical drop[END_REF]. Depending on the size of the electrode, the hot spot described earlier by Massing et al. moves. As the size of the electrode changes, the bubble covers the electrode to a greater or lesser extent, which changes the current density profile. The Joule effect is thus modified, which moves the hot spot. In the case of microelectrodes, the hot spot is close to the electrode. As the surface area of the electrode increases the hot spot moves upwards. This reveals a double vortex structure of the thermocapillary flow, which has not been taken into account before because the lower vortex is small at the microelectrode. When several bubbles develop simultaneously, the electric current must flow through the inter-bubble space. This causes an increase in the maximum current density near the inter-bubble equatorial region. The temperature hotspot is located near the equator and thus an almost symmetrical double vortex structure is generated near the bubble interface.

Influence of surface tension on electrogenerated bubbles

Surface tension as a parameter influencing nucleation

As a general rule, surface tension is an essential parameter to understand the formation of electrogenerated bubbles. Surface tension is a function of temperature , concentration of different surfactants, pressure and voltage. These empirically determined relationships are assumed to be linear and independent of each other [START_REF] Lubetkin | The motion of electrolytic gas bubbles near electrodes[END_REF][START_REF] Massoudi | Effect of pressure on the surface tension of water. Adsorption of low molecular weight gases on water at 25.deg[END_REF][START_REF] Weissenborn | Surface Tension of Aqueous Solutions of Electrolytes: Relationship with Ion Hydration, Oxygen Solubility, and Bubble Coalescence[END_REF][START_REF] Zhang | Evaluating the Behavior of Electrolytic Gas Bubbles and Their Effect on the Cell Voltage in Alkaline Water Electrolysis[END_REF]. Starting from the fact that Lubetkin pointed out that the increase in the concentration of different gases dissolved in water lowers the surface tension (including dihydrogen and oxygen), the nucleation phenomenon can be interpreted as follows:

• The concentration of dissolved gases increasing, the surface tension of the liquid gas interface decreases, which reduces the work required Eq. (2.34) and allows nucleation;

• After the nucleation event the concentration of dissolved gases decreases rapidly and the surface tension increases.

This brief reasoning need to be validated but shows the dependence of the evolution of surface tension on nucleation.

Surface tension as a parameter influencing growth

In the same register as shown in Fig. 7, the evolution of the surface tension influences the growth of the bubble. Surface tension is a function of temperature and concentration, among other things. As these two quantities vary at the interface, a variation of the surface tension is created. This variation has two plausible consequences. Firstly, it will generate a Marangoni effect which will modify the temperature and concentration gradients. Secondly it will modify the contact angle which will change the microconvection phenomenon generated by the spreading of the bubble. These two microconvection phenomena will affect the distribution of dissolved species around the bubble, thus modifying its growth. In the studies conducted on the nucleation, the growth, and the detachment of electrogenerated bubbles, only the surface tension of the gas-liquid interface is taken into account.

The resultant forces due to the surface tension of the 3 interfaces are estimated from a measurement of the contact angle. The problems of contact angles hysteresis are rarely addressed [START_REF] Brussieux | Controlled electrochemical gas bubble release from electrodes entirely and partially covered with hydrophobic materials[END_REF], which prevents a distinction from being made between the wettability of the material used and the roughness of the surface. With regard to the impact of surface tension on bubble formation, there is a need to specify its use in models on the formation of electrogenerated bubbles. As described above the dynamic contact angle cannot be predicted, and therefore must be considered as an input parameter of the model. The contact angle depends on the wettability of the electrode and on the surface tension of the hydrogen-electrolyte interface. The measurement of the contact angle requires particular attention and great rigour in the setting up of the experimental device. Any form of pollution that could affect the measurements must be avoided. The surface of the electrode must be clean and free of all forms of asperity that could affect the contact line [De Gennes and Brochard-Wyart, 2015]. The contact angle could be defined as the angle between the tangent to the interface at the point of triple contact with the solid surface. However, this measurement depends on the zoom with which the contact line is observed.

Yang et al. defined the contact angle according to geometrical parameters of the bubble, as it is shown on Fig. 8 . They measure the contact angle for one configuration in which the platinum electrode is embedded in a hydrophilic surface (glass) and another in which it is embedded in a hydrophobic surface (epoxy). The results are shown in the Fig. 9. As shown by the study conducted by Sakuma , the wettability of the surface will influence the effectiveness of the microconvection effect. A bubble on a hydrophobic surface will spread more easily, thus modifying the dynamic of growth. For oxygen electrogenerated bubbles Sakuma et al. observed initial sizes in the range of 10 -30 µm depending on whether the electrode surface is hydrophobic or hydrophilic. During these first moments the degree of supersaturation of the liquid as close as possible to the bubble varies greatly. The time coefficient does not vary from one electrode to another (hydrophilic or hydrophobic), so this means that in Eq. (2.45) β depends on the wettability of the surface. Consequently, the contact angle is another parameter to consider when considering bubble growth. 6 How to improve electrodes?

In the theoretical overview of the nucleation, growth and detachment of a bubble described in the previous section, it should be noted that because of the connection between these 3 steps, they must be studied together in order to be able to predict the detachment of the bubble. The conditions required to make an electrode surface more favorable to nucleation are those that, on the contrary, make it more difficult to detach the bubbles once formed.

Wilt gives several relationships to determine the nucleation rate for different surface geometries Fig. 11 [START_REF] Wilt | Nucleation rates and bubble stability in water-carbon dioxide solutions[END_REF]. These relationships have been established within the framework of the classical theory of heterogeneous nucleation. The discrepancy between the experience and the results provided by classical theory discourages its use as it stands. However, it makes it possible to understand that modifying the geometry makes it possible to modify the contact surface of the bubble with the electrode, to modify the ratio of the radius of curvature of the bubble with respect to its volume and thus the pressure difference and thus can make it easier to nucleate and detach it. By modifying the contact surface, we modify the contact pressure force F CP and the surface tension force F S , which are the main forces involved in the force balance [START_REF] Chen | Dynamics of single bubble departure from TiO2 nanorod-array photoelectrode[END_REF].

Using these two levers to facilitate detachment (wettability, and design of the electrode), Lubetkin imagines an electrode design that allows rapid nucleation and detachment of dihydrogen bubbles from the electrode surface, see Fig. 12. However, he suggests that another type of design could be more efficient and that this sophistication is not necessary

to achieve what he calls the phenomenon of "rapid-fire emission". This phenomenon has been observed by [START_REF] Weissenborn | Surface Tension of Aqueous Solutions of Electrolytes: Relationship with Ion Hydration, Oxygen Solubility, and Bubble Coalescence[END_REF]. He assumes that the rapid departure of a first bubble could cause sufficient disruption and reduce the concentration gradient, thus minimizing the marangoni effect.

Most current nucleation models do not take into account the contribution of the Marangoni effect or the impact of surface wettability and often focus on only one of its stages, so they are not able to predict the bubble detachment. A holistic approach is required. Most of the studies on the subject are experimental and require the use of high-speed cameras and microelectrodes to be able to observe the vortex flow of the fluid around the bubble which is characteristic of the marangoni effect. On the other hand, there are no numerical studies that model the development of an electrogenerated bubble and include the Marangoni effect on the liquid-gas interface of the bubble.

Integration into an inverse resolution problem

From the previous sections several questions emerge:

• How to understand the growth of a bubble with the Marangoni effect? Hypotheses about microconvection currents have been formulated, but there is no direct way to measure them. Based on this observation, other methods should be used to access information on bubble development that is not available through experiments. The presence of a Marangoni effect invalidates the hypothesis of transport of dissolved species by diffusion from the electrode to the bubble. • What is the origin of this Marangoni effect? Through experiments it is possible to measure: the diameter of the bubble, or the velocity of the fluid around the bubble. [START_REF] Lubetkin | The motion of electrolytic gas bubbles near electrodes[END_REF] : "A hypothetical composite conical or crack-like nucleation site. The region A is hydrophobic, possibly as a result of capillary condensation of a volatile organic material. This deposition mechanism would restrict its presence to the narrowest region of the pore. Nucleation is easy in this environment, with its high contact angle. As the nascent bubble grows past the hydrophobic region, and into the hydrophilic region B, a small remnant may be detached, remaining in the hydrophobic region. This residual gas would promote the instantaneous growth of the next bubble. As it grows, it will be shot out of the conical site, its zero contact angle ensuring the absence of friction. The rapidity of its departure might cause sufficient stirring in the vicinity C, to disrupt the Nernst layer, thus eliminating the concentration Marangoni effect, which would otherwise tend to hold the bubble near the electrode until it was larger. Alternatively, the rapid growth of the first bubble during its rise, and the immediate following of others might cause depletion of the Nernst layer in a narrow bubble chimney, again destroying the concentration gradient." However, it is not possible to access the concentration gradient of dissolved species, temperature, or to identify the presence of surfactant. From this point of view it is impossible to distinguish between a Marangoni effect of thermal solutal or surfactant origin.

• The Marangoni effect delays the detachment of the bubbles, but to what extent? While the presence of a Marangoni effect has been proven, it is currently impossible to measure the impact on bubble growth and bubble detachment. The diameter of the bubble at the time of detachment is an essential parameter for predicting the efficiency of the electrolytic process. However, the uncertainties on the phenomena influencing the evolution of this diameter cannot be resolved with the current state of theory and the available experimental techniques.

One of the solutions proposed in this thesis to answer the questions mentioned above is to pose an inverse problem. The objective here is to present quickly the inverse methodologies and how they can be used in the work related to this thesis. Inverse problems are problems whose formulation is incorrect. The calculated "solutions" are often quite different from the real solution. Nevertheless, the need to solve inverse problems is becoming more and more necessary.

Usually, in order to analyse the behaviour of a system as well as possible, it is necessary to build models to represent reality. These are determined from equations derived from physical laws and allow the behaviour of a system to be predicted under the effect of a known stimulus. A characteristic of these models is their causality: subsequent conditions depend on previous ones. When the input data and parameters are assumed to be known, solving the modelling equations is used to predict the output of the model.

To validate a model it is usual to compare the experimental results y mes with the modelling results y mod . If both fit, the assumptions are considered validated.

Inverse problems are the opposite of these direct problems. They are non-causal problems. They describe the situation in which one tries to determine the causes of a phenomenon from experimental observations of its effects. In other words, the objective may be to identify one or more of the elements that define the model.

The principle of an inverse methodology is to test several sets of parameters x in a model A to find those for which the results obtained by the model A(x) best fit the experimental results y mes : A(x) ≈ y mes (2.82)

Without going into the details of an inverse methodology it is necessary to discuss the measurable quantities in the case of an electrogenerated bubble on a micro-electrode and the unknown or supposedly known parameters. As described above, the experiments of Yang et al. allowed the observation of the temperature and velocity field around the bubble [START_REF] Yang | Marangoni convection at electrogenerated hydrogen bubbles[END_REF]. The temperature field is one of the quantities that can be used to determine the solutocapillary effect. The velocity field in a stagnant electrolyte is the result of the Marangoni effect. The other measurable quantity is the evolution of the bubble radius. These three quantities are the only ones that can be measured. The local concentration and the current density are determined by the model. One of the unknown parameters of the model is the presence of surfactants. This parameter is a priori very sensitive as it influences both the interfacial mass transfer and the surface tension gradient by inhibiting a large part of the bubble interface. As described above there is no certainty about the value of ∂γ ∂T , ∂γ ∂Φ , or ∂γ ∂c S the nature of the surfactant being unknown. The uncertainty in these parameters makes it impossible to determine the surface tension gradient correctly.

Usually in the framework of an inverse methodology it is advisable to carry out a sensitivity study on each of the parameters taken independently, as a small error on one parameter can strongly modify the output values. The diagram in Fig. 13 provides a non-exhaustive inventory of the sources of error that can occur in the process. In the spirit, beyond being an optimisation problem, the problem of an inverse methodology is to estimate the sources of error as well as possible, and their influence on the result. The aim of this section on inverse methods is to illustrate that one of the objectives of this thesis has been to build some of the necessary components of the inverse problem.

8 Overview of the chapter and objectives of the thesis

It has been shown in the previous sections that the major lever for improving the efficiency of the electrolysis process is to limit the coverage of bubbles on the electrode surface. For this it is necessary to accelerate the departure of the bubbles from the electrode. Consequently, a literature review on bubble nucleation, growth, and detachment was conducted in order to evaluate the phenomena involved at the microfluidic scale in the production of electrogenerated bubbles. It was found that the hypothesis of a Marangoni effect should be taken into account in order to understand the development of bubbles and thus be able to facilitate their detachment. However, several questions remain. It is impossible to determine the origin of the Marangoni effect. Is it solutal, thermal or due to surfactants? Furthermore, how can we understand the growth of a bubble with the Marangoni effect? In the course of the literary research carried out, it became clear that the experimental means are not sufficient to resolve these uncertainties. This is why the research track proposed in this thesis is to use inverse methods.

The work presented in this thesis is part of the inverse problem approach. Therefore, the objectives of this thesis are the following:

• To provide a mathematical or numerical model Numerical modeling 3

To model an electrogenerated two-phase fluid at the microfluidic scale, it is necessary to take into account the marangoni effect and mass transfer at the interface and to include the problem of the contact lines with the electrode. For this, multiphase microfluidic modeling ingredients are needed, in particular to model phases, mass transfer and the marangoni effect at interfaces. A review on the subject was done by [START_REF] Wörner | Numerical modeling of multiphase flows in microfluidics and micro process engineering: a review of methods and applications[END_REF]. There are several types of methods for modelling a multiphase fluid. They all have their advantages and disadvantages. Some are more suitable for microfluidic flows.

The first objective of this chapter is to describe these methods and to evaluate their suitability for the problem at hand. Then, the mathematical basis for the development of the numerical model is given. Finally, the numerical model created and the reasons for the choices made during its design will be presented.

Toward a direct numerical simulation of the phenomenum

In order to explain the reasons that led to the choice of the specific numerical method, it is appropriate first to give a brief description of all the methods for modelling a multiphase system and then to state the necessary requirements of the model.

Overview of the numerical methodology

To model an electrogenerated two-phase fluid at the microfluidic scale, it is necessary to take into account the marangoni effect and mass transfer at the interface and to include the problem of the contact lines with the electrode. For this, multiphase microfluidic modeling ingredients are needed, in particular to model phases, mass transfer and the marangoni effect at interfaces. A review on the subject was done by [START_REF] Wörner | Numerical modeling of multiphase flows in microfluidics and micro process engineering: a review of methods and applications[END_REF]. There are several types of methods for modelling a multiphase fluid. They all have their advantages and disadvantages. Some are more suitable for microfluidic flows. The aim here is to describe these methods and to assess their suitability for the problem at hand.

Modelling an interface

Reproducing and following an interface with complex shape and dynamics that can develop large deformations, singularities and topological changes is a numerical challenge. The interface to be described (gas-liquid) has material properties, i.e. it changes the flow behaviour, as is the case with surface tension. There is a zone of discontinuity for pressure and properties between the two phases. If the interface has material properties, there may be diffusion of these properties along the interface, or transfer between phases. Several difficulties must be overcome when dealing with the numerical simulation of interfaces:

• the domain of each phase is unknown;

• the deformation of the interface is governed by a discontinuous pressure jump condition;

• the very shape of the interface has a direct influence on the action of surface tension;

• in the case of the phenomena of adsorption-desorption of dissolved species on the interface or coalescence, we are at the limit of the validity of macroscopic description.

Interface characteristics

Tryggvason et al. point out that in order to simulate gas-liquid multiphase flows it is necessary to make several hypotheses [START_REF] Tryggvason | Direct Numerical Simulations of Gas-Liquid Multiphase Flows[END_REF]]:

• the length scale of the problem is much larger than the mean free path of the molecules, so the continuum hypothesis is valid and the flow equations take the classical form;

• the interface separating two or more fluids, which in fact has a finite thickness of the order of a nanometer and constitutes a transition zone for the properties of the fluid, is assumed to be sharp with a negligible thickness;

• the intermolecular forces, which determine the dynamics of the interface, are modelled on a continuum scale as a force localised at the interface which is proportional to a surface tension coefficient, which represents the rate of change of the excess free surface energy produced by a unitary increase in the surface area of the interface resulting from the deformation.

Explicit or implicit representation of the interface

An overview of the methods for simulating a multiphase fluid is shown in Fig. 14 [START_REF] Wörner | Numerical modeling of multiphase flows in microfluidics and micro process engineering: a review of methods and applications[END_REF]. There are two main classes of interface processing methods: Lagrangian methods and Eulerian methods (volume tracking). The choice of an efficient and robust method to take into account the interface depends on the physical problem to be studied, as each method has its strengths and weaknesses. The difference between these two classes of methods lies in the representation of the interface: explicit or implicit.

In Lagrangian methods the interface is explicitly tracked. Indeed, in these methods, the mesh adapts over time so as to merge with the interface. While Eulerian methods use a fixed mesh, with an interface that is not explicitly tracked but reconstructed from a phase indicator function or color function. Lagrangian methods maintain the interface as a discontinuity, and explicitly track its evolution. No modeling is required to define the interface or its effects on the flow. Moreover, boundary conditions can be applied exactly to the interface. However, these methods require a remeshing at each time step. It should be noted that in case of strong distortions of the interface the mesh may be strongly altered and not uniformly distributed, which may degrade the accuracy of the resolution methods. The main disadvantage of these methods lies in their difficulty to take into account topological changes and in particular the ruptures or coalescence.

Eulerian methods, also called fixed grid methods, front capturing methods or volume tracking methods, require modelling where additional equations are needed to obtain information on the location of phases and discontinuities. Indeed, unlike the Lagrangian methods, the interface is not explicitly tracked in the Eulerian methods. To locate the different phases and impose the interfacial conditions, a phase indicator function or color function is introduced. This phase indicator function is defined over the whole computational domain and allows to locate the different phases. At each time step, the interface can be located and reconstructed from this phase indicator function. The phase-indicator function allows boundary conditions to appear in the flow equations, but these boundary conditions are altered. Indeed in these methods the interface is diffuse and of non-zero thickness. Therefore the information on the interface is smooth, leading to a spreading and dispersion of the information. Eulerian methods have the advantage of not needing remeshing procedures. In addition, these methods automatically take into account topology changes such as coalescence and fragmentation. However, due to the smoothing of the interface, the physical phenomena related to the interface are not described precisely, as the interface is not explicitly represented.

Eulerian methods

Eulerian methods include: Volume Of FLuid (VOF), Level-set (LS). These are the most common methods used in CFD software (FLUENT, COMSOL, OpenFOAM). VOF methods naturally ensure the conservation of volume and mass in incompressible flows and, with some improvements, in compressible flows. However, the description of the interface is not precise, which makes it difficult to evaluate the curvature of the interface and impose boundary conditions. The level-set method, like the VOF methods, automatically takes into account topological changes. It describes the interface implicitly using a signed distance function which gives a more precise definition of the interface than in standard VOF methods.

But the signed distance function has to be reset frequently alterating mass conservation.

Finally the Lagrangian methods are very precise with a thickness free interface and boundary conditions that are imposed exactly on the interface. However, changes in topology and highly deformed interfaces are not easily accessible by this type of methods because the remeshing procedure which allows to preserve the adequate mesh size can become very complicated in this case.

Requirements of the simulation

Towards a holistic approach

Numerical models to simulate the Marangoni effect along the wall of an electrogenerated bubble have already been presented [START_REF] Massing | Thermocapillary convection during hydrogen evolution at microelectrodes[END_REF]. But these consider the interface as fixed, and assume that the bubble does not grow by considering only a small time interval during its development. Liu et al. have modeled the growth and detachment of an electrogenerated bubble [START_REF] Liu | Numerical simulation of hydrogen bubble growth at an electrode surface[END_REF]. But they point out that their numerical results do not match the experimental results, partly because the marangoni effect is not modelled.

The real added value of a numerical model is to combine bubble growth with the marangoni effect.

Dealing with numerical errors

The model must be reliable, and must, among other things, limit problems related to mass conservation and spurious currents. Spurious currents are proportional to capillary effects, and in this case, in a simulation where the capillary effects are at the origin of the main forces governing the evolution of the system, they compete directly with the currents caused by the Marangoni effect, thus distorting the results. It is therefore important, if it is not possible to eliminate them completely, to at least know their importance and to reduce their effects. The spurious currents depend, among other things, on surface tension, viscosity, and discretization techniques used. They increase proportionally with decreasing capillary number Ca = µv 0 γ [START_REF] Ren | Heterogeneous multiscale method for the modeling of complex fluids and micro-fluidics[END_REF][START_REF] Herrmann | A balanced force refined level set grid method for two-phase flows on unstructured flow solver grids[END_REF][START_REF] Dupont | Numerical simulation of static and sliding drop with contact angle hysteresis[END_REF] . In the simulation we want to run, the capillary effects are predominant. The problem is that most of the CFD software available uses models close to the Brackbill model [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF]].

Choice of a numerical method

Choosing an appropriate numerical method turned out to be complex because of the difficulties in modelling an interface at the microfluidic scale. There are a multitude of choices and several criteria to consider. If it is possible to express the qualities required for a simulation in terms of equations, the choice of a numerical method calls for another type of knowledge. What is needed is a numerical method capable of "reasonably" simulating the problems we are interested in. It is not simply a question of modeling an interface but of being able to simulate the marangoni effect with mass transfers while managing the problems related to contact angles. The growth and detachment of a bubble from the wall of an electrode involves large deformations and a change of topology, which excludes moving mesh methods. The final choice was a VOF method. The advantage of this method is that it limits errors due to volume conservation. Even if the curvature calculations are not very accurate, there are solutions to improve it and thus limit the spurious currents. Combining the methods described by Guo et al. with those of Ivana Seric et al. provides a good method of mass conservation (VOF) with a means of accurately calculating the curvature (Height function), while being able to model the growth of a bubble with the Marangoni effect [START_REF] Guo | Implementation of a height function method to alleviate spurious currents in CFD modelling of annular flow in microchannels[END_REF][START_REF] Seric | Direct numerical simulation of variable surface tension flows using a Volume-of-Fluid method[END_REF]. This is a first step towards a holistic model.

Mathematical model

In order to understand the choices of the numerical methodology used, it is necessary to recall the foundations of the mathematical model. In a mono-fluid system the Navier-Stokes equations and boundary conditions make it possible to describe the whole flow, however in a two-phase liquid-gas system, it is necessary to reapply the principle of conservation of mass and momentum to be able to describe what happens at the interface.

Mathematical operators

Before the conservation principle can be applied to the two-phase system, a jump operator and a surface gradient operator must first be defined. The Reynolds transport theorem must also be developed for systems with mass transfer at the interface.

Jump operator

The jump operator, as shown in Fig. 15, describes the passage of a quantity Q through an interface between two distinct volumes. Assuming that Q has a limit on each side of the interface I, we define for any point x belonging to I, the jump relation for Q by :

Q = lim h →0 + Q l (x + h • n) -Q g (x + h • n) (3.1)
where n is the vector normal to the interface. The direction of the vector n is chosen here arbitrarily from gas to liquid.

Surface gradient operator

The surface gradient operator ∇ s Q is defined as the projection of the gradient ∇Q onto the surface :

∇ s Q = ∇Q -n(n • ∇Q) (3.2)
The diagram in the fig. 16 gives a geometric representation of this operator. The gradient ∇Q is projected onto the tangent at the interface.

74

Chapter 3

Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022

Two-phase particle

The Fig. 17 shows a fluid particle consisting of a liquid phase and a gas phase (e.g. liquid l, gas g). The two phases are separated by an interface I. In order to show the boundary conditions at the interface, the reynolds transport theorem is first applied to this bi-fluid volume. The total volume is divided into two sub-volumes:

V(t) = V g (t) + V l (t) (3.3)
Green-Ostrogradsky bi-fluid theorem

For a fluid particle the Green-Ostrogradsky theorem allows to link the flux of a vector field E(x, t) to the divergence of this field.

∂V(t) E(x, t) • n ∂V dS = V ∇ • E(x, t) dV (3.4)
The surface integral applied on the boundary ∂V of the scalar product of E with n ∂V can be decomposed. The Green-Ostrogradsky theorem can then be applied to the two boundaries of the gas and liquid volumes. This decomposition into two fluids reveals a jump relationship at the interface :

∂V(t) E(x, t) • n ∂V dS = ∂V g E • n ∂V g dS - I E g • (-n)dS + ∂V l E • n ∂V l dS - I E l • n I dS = V g+l ∇ • E dV + I E • n dS (3.5)
Two-phase reynolds transport theorem with mass transfer at the interface For a volume V with a single-phase fluid, the time derivative of the volume integral of a scalar field Q over V(t), dependent on the time t, is equal to the sum of the volume integral of the time partial derivative of Q and the surface integral over the integration volume boundary ∂V(t) , which is time-dependent, of the product of Q with the normal displacement velocity of the boundary v ∂V • n ∂V : Then using Eq. (3.6) on the two separate volumes we obtain:

d dt V(t) Q(x, t)dV = V(t) ∂Q(x, t) ∂t dV + ∂V(t) Q(x, t) v ∂V • n ∂V dS (3.6)
d dt V(t) Q(x, t)dV = d dt V g (t) Q(x, t)dV + d dt V l (t) Q(x, t)dV (3.7) d dt V(t) Q(x, t)dV = V g (t)
∂Q(x, t) ∂t dV

+ V l (t) ∂Q(x, t) ∂t dV + ∂V g (t) Q(x, t) v ∂V g • n ∂V g dS + ∂V l (t) Q(x, t) v ∂V l • n ∂V l dS (3.8)
A surface integral on the interface I between the two phases appears knowing that :

∂V(t) = ∂V l (t) -I(t) + ∂V g (t) -I(t) (3.9)
The vectors n ∂V l and n ∂V g are directed outwards from the volumes V l and V g , and are normal to both surfaces ∂V l and ∂V g . At the interface level they are opposite. The previous equation becomes:

d dt V(t) Q(x, t)dV = V g+l (t)
∂Q(x, t) ∂t dV

+ ∂V(t) Q(x, t) v ∂V • n ∂V dS + I(t) Q(x, t) g v ∂V g • n ∂V g dS + I(t) Q(x, t) l v ∂V l • n ∂V l dS (3.10)
The speed of movement of an interface is generally identical to the speed of movement of adjacent continuous media. However, when there is mass transfer or phase change at its location, the interface no longer moves at the same speed as the surrounding continuous media, so v ̸ = v I ,where v is the speed of the continuous media, and v I is the speed of the interface. Locally at the closest to the interface the following equalities apply:

v ∂V g = v ∂V l = v -v I (3.11)
Elsewhere the velocity of the volume boundary is the same as that of the continuous medium.

Thus we have v = v ∂V l , v = v ∂V g , v = v ∂V .
The vector n is defined as being normal to the interface and being directed from the liquid phase to the gaseous phase:

n ∂V l = -n ∂V g = n (3.12)
From the Eq. (3.10) we get:

d dt V(t) Q(x, t)dV = V g+l (t) ∂Q(x, t) ∂t dV + ∂V(t) Q(x, t) v • n ∂V dS + I(t) Q(x, t) l (v -v I ) • n I dS - I(t) Q(x, t) g (v -v I ) • n I dS (3.13)
The Green-Ostrogradsky's theorem applied to the surface integral yields, i.e. the second right hand term of Eq. (3.13):

∂V(t) Q(x, t) v • n ∂V dS = V(t) ∇ • (Q(x, t) v)dV (3.14)
Then the jump relationship for Q(x, t) (vv I ) is formed with the surface integrals applied to the interface:

d dt V(t) Q(x, t)dV = V(t) ∂Q(x, t) ∂t + ∇ • (Q(x, t) v) dV + I(t) Q(x, t) (v -v I ) • n dS (3.15)

Conservation principle

Mass conservation

In the same way as in the single-phase case, the mass of a small volume element is preserved, it is considered that there is no mass source term in the volume or on the surfaces:

d dt [m(t, x)] = d dt V(t) ρ(x, t)dV = 0 (3.16)
The equation (3.15) then implies that:

d dt V(t) ρ(x, t)dV = V(t) ∂ρ ∂t + ∇ • (ρ v) dV + I(t) ρ (v -v I ) • n dS = 0 (3.17)
The localization principle states that a null integral for any volume V (t) implies that its integrand is null. This establishes the principle of local conservation of the mass. For the volume V(t) we obtain :

∂ρ ∂t + ∇ • (ρ v) = 0 (3.18)
And at the interface I(t) for mass conservation the jump relation is :

ρ (v -v I ) • n = 0 (3.19)
This last relation reflects the equality of the mass flows on each side of the interface. By this means an input parameter of the model is introduced, the mass flow rate of mass transfer at the interface ṁ in [kg.m -2 .s

-1 ]. ṁ = ρ g (v g -v I ) • n = ρ l (v l -v I ) • n (3.20)

Momentum conservation

The temporal variation in the momentum mv of a fluid particle results from the forces acting on this particle. As in the single phase, the forces affecting the particle are divided into surface forces F s and volumic forces F v . However, a new type of force is added, the linear forces F l , which are exerted exclusively on the edge of the interface ∂I. By defining the surface tension γ, we obtain for the three forces:

F l = ∂I γ n ∂I dL (3.21) F s = ∂V T • n ∂V dS (3.22) F v = V ρ f v dV (3.23)
With T the stress tensor, and f v the density of the volume forces. Before establishing the law of conservation, the linear force is rewritten as a surface integral. By applying Stokes' theorem, which is the surface version of Green-Ostrogradsky's theorem, we transform the linear integral into a surface integral :

F l = I (n × ∇) × (γ n) dS (3.24)
Then , the operators surface gradient ∇ s , surface divergence ∇ s • and the definition of the curvature of the interface κ such as κ = ∇ s • (-n) allow to write the linear force as :

F l = I ∇ s γ + γκn I dS (3.25)
This expression gives rise to the Marangoni term ∇ s γ which is expressed in the case of spatial variations in surface tension, and the term γκn I which is the Laplace pressure responsible for the pressure difference between the two phases. Now that the linear force is expressed as a surface integral, Newton's second law is applied to the fluid particle :

d dt V(t) ρ(x, t) v(x, t)dV = ∂I γ n ∂I dL + ∂V T n ∂V dS + V ρ f v dV (3.26)
Then, the Eq. (3.15) is applied to the left-hand term, the linear force term is rewritten, and the Eq. (3.5) is applied to the surface force term:

d dt V(t) ρ(x, t)v(x, t)dV = V(t) ∂ρv ∂t + ∇ • (ρ v ⊗ v) -∇ • (T) -ρf v dV + I(t) ( ρv ⊗ (v -v I ) • n -T • n -γ κn -∇ s γ) dS (3.27)
Locally in each phase we have:

∂ρv ∂t + ∇ • (ρ v ⊗ v) = ∇ • (T) + ρf v (3.28)
And at the interface :

ρv ⊗ (v -v I ) • n -T • n -γ κn -∇ s γ = 0 (3.29)
The stress tensor can be decomposed into a pressure component and a viscous stress tensor:

T = -pI + 2µD (3.30)
where µ is the dynamic viscosity, and D = ∇v + (∇v) T /2. As a reminder getting mass conservation at the interface ρ (vv I ) • n = 0 and so ρv ⊗ (vv I ) • n can be expressed as ṁ v . At the interface by displaying the mass flow rate and developing the stress tensor we finally obtain:

ṁ v + pI • n -2µD • n = γκn + ∇ s γ (3.31)
Projections along a normal axis and an axis tangential to the interface results in:

ṁ v • n + p -2µD • n • n = γκ (3.32) ṁ v • t -2µD • n • t = ∇ s γ • t (3.33)
As shown in Eq. (3.32) and Eq. (3.33) the flow at the interface is influenced by surface tension, the Marangoni effect, and mass transfer. The term ṁ v refers to the consequences of mass transfer across the interface on the tangential and normal stress balance. This means that it is possible to have slippage at the interface caused by the mass transfer.

3 From continuous to discrete model

Having laid the foundations of the mathematical model, it is necessary to describe how the transition from the mathematical model to the numerical model is made.

FVM -volume of fluid

Averaging of continous quantities on a finite volume

In order to move from a mathematical model to a discrete model, it is necessary to transform the partial differential equations representing the conservation laws into discrete algebraic equations. In the finite volume method the continuous quantities (or functions of space and time in mathematical terms) are averaged over each volume represented by the cells of the mesh.

β = 1 V cell V cell β i dV (3.34)
The dynamics and continuity equations resulting from the conservation principle are adapted to the single field approach using a weighted average of the quantities followed by these equations. In the process of calculating the weighted average, the operation carried out for the calculation of the weighted average on the divergence and gradient operators is not simple and straightforward. The integrals and derivatives can not be interchanged in volumes that contain the interface I. The spatial volume averaging theorem has to be used [START_REF] Whitaker | Theory and applications of transport in porous media: The method of volume averaging[END_REF]:

∇β i = ∇β i + 1 V cell A I,cell β i ndA (3.35) ∇ • β i = ∇ • β i + 1 V cell A I,cell β i • ndA (3.36)
where A I is the part of the surface of the interface in the cell.

Volume fraction

Instead of considering two fluids (gas-liquids), the the Volume Of Fluid (VOF) method assume a single fluid, and solve a single conservation equation for mass and momentum. A continuous indicator function 1 i is used to distinguish the phase i from the others. It takes the value 1 when phase i is present at a given point in the system and 0 otherwise. In the VOF methodology the integral on the cell of the function defines the volume fraction α i :

α i = 1 V cell V cell 1 i dV (3.37)
In the cells that contain the phase i α i = 1, and 0 otherwise. A notable consequence of this kind of multiphase modelling is that the interface I is diffuse and spans several cells of the mesh.

Transport of the volume fraction

A transport equation is added to the Navier Stokes equations to track the volume fraction of phase i, which makes it possible to distinguish the phases and identify their interface :

∂α(t) ∂t + ∇ • (α(t) v) = 0 (3.38)
In cells crossed by the interface I, we get 0 < α i < 1. This equation is common in CFD and is routinely solved. This allows precise control of the volume of each phase in each cell of the mesh at a lower computational cost. Therefore, the VOF methodology is known to have good volume conservation properties. Several immiscible fluids are considered as one effective fluid throughout the domain. This is referred to in the literature as a one-fluid formulation. A notable consequence of this modelling is that the interface I is diffuse and spreads over several cells of the mesh. This diffusivity of the interface is a handicap of the VOF methodology. Surface tension is a local force whose points of application are directly on the interface and which varies along the interface. In this context, it may be considered counter-intuitive and inacurate to choose a multiphase model that diffuses the interface to account for the impact of surface tension on it.

Continuum surface force formulation

There are several formulations for modelling surface tension. In a literature review Popinet describes the recent development of the Eulerian surface tension formulations [START_REF] Abu-Al-Saud | A conservative and well-balanced surface tension model[END_REF]. One of the most common formulation of surface tension for VOF was first introduced by Brackbill et al., and is called "Continuum Surface Force" (CSF) [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF]].

The formulation adopted in this thesis is based on it. It is therefore necessary to describe it first to understand its shortcomings and to understand the modifications made to overcome them.

Surface tension expressed as a volumetric force

In this formulation, the surface tension force is modified into a volumetric force and introduced into the momentum equation. This surface tension source term can be calculated from the values of the volume fraction and the surface tension coefficient. There are different models. One of the most commonly used is the continuum surface force model CSF. By using a Dirac delta function, the surface integral from Eq. (3.25) taking into account the surface tension becomes a volume integral :

F l = I ∇ s γ + γκndS = V (∇ s γ + γκn) δ S dV (3.39)
This volumetric formulation makes it possible to integrate surface tension as a source term in the momentum equation Eq. (3.28):

f γ = (∇ s γ + γκn) δ S (3.40) ∂ρv ∂t + ∇ • (ρ v ⊗ v) = ∇ • (T) + ρf v + f γ (3.41)

Interface density

The operator δ S represents the interface density present in the cell. V ∩ I is the proportion of the interface present in the volume V:

V ∩ I = (V∩I) dS = V |∇α|dV (3.42)
In three spatial dimensions this volume integral gives the interface area, in two dimensions the corresponding area integral gives the interface length. |∇α| being constant within each cell of the mesh, V cell the previous integral is transformed into :

V cell |∇α|dV = V cell |∇α| = V cell ∩ I = I cell (3.43)
I cell represents the part of the interface present in the cell. In each cell of index (i, j), |∇α (i,j) | = I cell V cell represents the interface density of the cell and so:

δ S = |∇α| (3.44)

Normal and curvature

The geometrical properties of the interface, the normal vector and the curvature can be calculated from the gradient of α.

n = ∇α |∇α| (3.45) κ = ∇ • n = ∇ • ∇α |∇α| (3.46)
The volume surface tensionforce f γ can be expressed using the CSF formulation (neglecting the Marangoni term), that is, using |∇α| :

f γ = γκ n δ S = γ κ ∇α |∇α| |∇α| = γκ∇α (3.47)

CSF Marangoni effect

Very often authors seeking to model the effects of surface tension assume that it is constant. And consequently, they do not express the marangoni effect, and remove the expression ∇ s γ from the Eq. (3.40). Surface tension is a function dependent on temperature and concentration, amongs other things. A value of the surface tension can thus be calculated in each cell of the mesh. In other words, a value of the surface tension coefficient can be calculated in each cell of the mesh even if the interface does not pass through it. One can thus calculate a gradient of surface tension. From the definition of surface tension we get:

∇ s γ = ∇γ -n(n • ∇γ) = ∇γ - ∇α |∇α| ∇α |∇α| • ∇γ (3.48)
The CSF is one of several methods to model surface tension. However, whatever the chosen method, several difficulties have to be overcome :

• appearance of spurious currents;

• failure to preserve mass;

• contact line limitations.

Numerical errors

In order to properly assess the errors that could be generated by the model, it is necessary to state in this section the common errors encountered in multiphase fluid modelling: parasitic currents, and conservation of volume. This description helps to confirm the choices made in the design of the model.

Spurious currents

Observation

Spurious currents are currents observed in a numerical simulation that have reached a steady state of equilibrium when no energy is injected into the system. These currents result from errors in the discretization of the surface tension and have serious consequences on the results of the calculations. The first consequence is that they make certain calculations impossible [START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF][START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF] • the existence of mechanical equilibrium when the fluids are at rest;

• the normal and curvature are estimated accurately;

• the consistency of the discretization operators. Starting from Eq. (3.32), and assuming v = 0 on either side of the interface, we get the Young-Laplace equation -p = γκ . This corresponds to the case of a bubble or drop in zero gravity in a stagnant fluid. By applying the hypotheses related to this systems, the Navier-Stokes equations are simplified: Newtonian, considered incompressible, no mass transfer at the interface, no gravity, constant surface tension, the flow is isothermal. The volumetric forces term cancels out in the momentum equation. Since there is no mass transfer at the interface, the velocity of the two fluids on either side of the interface is equal to that at the interface. The two phases are in equilibrium, which theoretically means that the velocity is zero on both sides of the interface. This means that if the fluid is moving near the interface it can only be due to a numerical error. As shown in fig. 18, spurious currents appear when the CSF formulation is used and they can be significant.

Well-balanced relation

It is assumed that there is a discretized interface geometry such that a zero velocity field is the solution to the Navier-Stokes equations. The mechanical equilibrium of the discrete system at zero velocity is characterized by the momentum equation where all velocity-dependent terms are removed .

∂ρv ∂t + ∇ • (ρ v ⊗ v) + ∇p -∇ • (2µD) = f γ (3.49) -→ ∇p = f γ (3.50)
This leaves a balance between the discrete pressure gradient, and the surface tension term.

In the absence of gravity and for a constant surface tension , we obtain [START_REF] Abu-Al-Saud | A conservative and well-balanced surface tension model[END_REF] :

-∇p + γκn I δ I = -∇p + γκ∇α = 0 (3.51)
Assuming a constant curvature we obtain :

-∇ * (p -γκα) = 0 (3.52)
where ∇ * is a numerical approximation of the gradient. The exact discrete numerical solution, which guarantees exact balance between surface tension and pressure in the case of constant curvature and surface tension, is then simply :

p -γκα = constant (3.53)
In the case of a numerical simulation of a static bubble in a liquid in a zero gravity environment and giving the curvature (which is constant) as an input parameter of the simulation rather than trying to calculate it, this relationship should be logically verified. If this is not the case, an imbalance is created by the discretization system.

Popinet specifies that for the equilibrium condition to be met, the pressure gradient should be estimated using the same discrete operator as that used to estimate the gradient of the indicator function used in the volumetric surface tension force calculation [START_REF] Abu-Al-Saud | A conservative and well-balanced surface tension model[END_REF].

He points out that Brackbill et al. in the original CSF article uses two different operators to calculate the pressure gradient whose values are taken from the centre of the faces and the gradient of the volume fraction whose values are taken from the centre of the cells [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF]. They calculate the surface tension force at the center of the faces by averaging the values taken at the center of the cell to perform the force balance. The values used to calculate the discretized gradient of pressure and volume fraction must be taken at the same location in the mesh in order to get a well-balanced relation. If this is not the case, the discrete operators used to calculate the pressure gradient and volume fraction are not the same and an imbalance is created. However, the use of this discrete operator does not give a sufficient approximation of the gradient of the volume fraction to be able to estimate the interface normal and the curvature accurately. Popinet suggests calculating the curvature using another method than the one appearing in the CSF method, i.e. without using the gradient of the volume fraction [START_REF] Abu-Al-Saud | A conservative and well-balanced surface tension model[END_REF].

Abadie et al. point out that any discrete vector field is not the gradient of a scalar field [START_REF] Abadie | On the combined effects of surface tension force calculation and interface advection on spurious currents within Volume of Fluid and Level Set frameworks[END_REF]. Its curl operator must be null. Looking at Eq. (3.50), the curl of the second member must be equal to 0 (i.e. ∇ × (f γ + ρg) = 0 ) in order to be exactly equal to the discrete gradient pressure which is a scalar field. This mathematical condition is a prerequisite for equilibrium. However, it is not fulfilled if we combine the approaches mentioned in the CSF model to calculate f γ . However, the gravity term ρg is not the gradient of a scalar, even in the discrete sense. Only the sum f γ + ρg is written as a gradient at equilibrium. These two terms must therefore be discretized together in order to hope to achieve a balance of the discrete system. Some authors have addressed this issue and implemented well-balanced algorithms in the VOF, Level-set, front-tracking framework [Francois et

Time step condition

For explicit schemes to ensure that spurious currents do not develop over time, a stability condition on the time step introduced by Brackbill et al. must be applied [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF]:

∆t < ρ avg (∆ x ) 3 2 πγ (3.54)
where ∆ x is the grid spacing, γ the surface tension, ρ avg is the average density of the both phases. The physical reason given by Brackbill et al. is that the time step must be small enough to resolve the fastest capillary waves. The value of the time step is limited by the size of the mesh. As shown by Eq. (3.54) there is a proportional relationship between the time step and the grid spacing, ∆t ∝ (∆ x ) 3/2 . See [START_REF] Abu-Al-Saud | A conservative and well-balanced surface tension model[END_REF] for a detailled discussion on the subject.

Volume conservation

Transporting materials by an incompressible flow results in the conservation of volume and mass. This property is a very important stake for the VOF method. By nature, the VOF method has good volume conservation properties, but some steps of the algorithm are approximate and the conservation is inaccurate in most existing algorithms [START_REF] Aulisa | A geometrical area-preserving Volume-of-Fluid advection method[END_REF] .

A natural definition of mass conservation is a method which conserves the total area at each time step so that:

ij A α n ij = ij A α n+1 ij
The advection of the interface requires the calculation of cell boundary fluxes. The use of volume of fluid data and fluxes should lead directly to exact mass conservation but it is in fact not so.

Moreover, interface advection algorithms may produce some systematic errors, such as volume fractions that do not satisfy :

0 ≤ α n ij ≤ 1
The above-mentioned inconsistencies are difficult to correct: it is not obvious where the excess or missing mass should be disposed of, or retrieved. Code writers then routinely redistribute it in the surrounding cells with some diffusion algorithm or reset the volume fraction to 1 or 0 thus destroying mass conservation.

Various attemps have been made in order to assure boundedness and conservativeness of the phase fraction [START_REF] Scardovelli | DIRECT NUMERICAL SIMULATION OF FREE-SURFACE AND INTERFACIAL FLOW[END_REF][START_REF] Cummins | Estimating curvature from volume fractions[END_REF][START_REF] Afkhami | Height functions for applying contact angles to 2D VOF simulations[END_REF].

The error on the volume balance has important consequences on the simulation of bubble growth on walls. Even if the error is smaller, the growth time and the diameter at detachment can be poorly predicted. Therefore, the accuracy of the phase volume balance is an indispensable ingredient of the numerical method.

By nature, the VOF method has good volume conservation properties, but some steps in the algorithm are approximate and conservation not satisfied in most existing algorithms .

The first challenge to overcome in the context of multiphase fluids is the modelling of the pressure jump at the interface. After having described in the previous section the methods usually used to manage this discontinuity at the interface and the spurious currents they can generate, it is appropriate to present here the method that has been chosen to carry out the simulations of chapter 4 and 5.

The objective of this section is to calculate γκn by an alternative and more efficient method than the CSF methodology. For this purpose a first method based on height functions has been described. The use of this methodology allows to determine the curvature, the normal and the tangent at the interface. However, the use of the height function methodology is not straightforward. A first necessary step is to define local coordinate systems to retrieve the data needed to calculate the height functions, as described in subsection 5.2. Furthermore, the height function method loses accuracy when the slope of the interface is too steep compared to the local coordinate system. The methodology for choosing between the vertical or horizontal local coordinate system is described in subsection 5.3. Despite a good choice of the coordinate system the use of a second method based on a polynomial fit described in subsection 5.4 is necessary to obtain an accurate calculation of the curvature. The choice of the transition from one method to the other is not straightforward and requires tests which will be described in chapter 4. Finally, the last subsection describes how the term γκn is integrated into the Navier-Stokes equations.

Height definition

Definition of the height

The height function method method makes the calculation of normal and curvature more accurate, thus reducing parasitic currents [START_REF] Poo | A computational method for determining curvatures[END_REF] As shown in Fig. 19 the integral is represented by the air under the curve. Typically on a mesh crossed by an interface this value is given exactly by the volume fraction. The volume fraction gives the part of the cell area occupied by a phase. A stencil of several cells around the cell through which the interface passes and for which the curvature is to be calculated is used as the basis for a coordinate system. In each cell of width ∆ x and height ∆ y , we have the value of the volume fraction occupied by one of the phases. Summing the volume fractions of a column of the stencil and multiplying it by the width of the cells makes it possible to carry out a calculation equivalent to the calculation of an integral of a function whose graph is represented by the interface in a local coordinate system represented by the stencil :

f (x 0 ) = H(x 0 ) = x 0 + ∆x 2 x 0 -∆x 2 f (x)/∆ x dx (3.55) H(x 0 ) = H(i) = j=+∞ j=-∞ α ij • ∆ y (3.56)

Curvature calculation

From this quantity we can obtain an approximation of the first and second derivative of the function in the ith-column by using the value of the heigh in left (i -1) and the right (i + 1) column of the stencil.

H ′ (i) = H(i + 1) -H(i -1) ∆ x (3.57) H ′′ (i) = H(i + 1) -2H(i) + H(i -1) ∆ 2 x (3.58)
An estimation of the tangential, the normal vector, and the curvature can be obtained from there :

t = 1 1 + H ′ (i) 2   1 H ′ (i)   (3.59) n = 1 1 + H ′ (i) 2   H ′ (i) 1   (3.60)

92

Chapter 3

Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022

The curvature is calculated as the negative of the divergence of the normal vector :

κ = -∇n = - H ′′ (i) (1 + H ′ (i) 2 ) 3/2 (3.61)
The key to the success of this method is having access to sufficiently accurate discrete values of the height function [START_REF] Francois | Interface curvature via volume fractions, heights, and mean values on nonuniform rectangular grids[END_REF]. So the smaller the cell width, the better the approximation.

An issue arises when the slope of the interface tends to infinity. In mathematical terms the function is no longer defined. In a more general way according to the Popinet analysis, the error on the calculation of the curvature increases with the magnitude of H(i) ′ . For each cell the coordinate system is oriented so that H(i) ′ ≤ 1. Usually a vertical stencil and a horizontal stencil are used. The stencil used is the one in which the interface is represented by the function graph with the lowest slope.

Definition of a local coordinate system

The use of the height function methodology is not straightforward. A first necessary step is to define local coordinate systems to retrieve the data needed to calculate the height functions.

The goal here is to capture the interface to be able to assimilate it to a mathematical function. By using stencils the data necessary to calculate the curvature of the interface in each cell can be retrieved. However, it is necessary to adapt the size and direction of the stencil, and in some cases modify the data to fit the mathematical formalism of the analytical calculation.

To be able to use Eq. (3.58), it is necessary to retrieve the values of the volume fractions of the cells allowing to calculate 3 contiguous heights. This is why the stencil is usually made up of 3 cells in width. For the choice of the number of cells to be considered in height, the choices differ according to the authors. Guo et al. chose to start on a 3 × 7 stencil base.

Popinet adapts the height of the number of cells of each column according to the need. The essential thing is to be able to get the integral of the interface in each column correctly. The volume fraction of the lower cell must be equal to 1 and the volume fraction of the upper cell must be zero, as shown on Fig. 20.

One of the problems to take into account when switching from an interface to the graph of a function is that for one x-coordinate the interface can have several y-coordinates. Two cases oriented horizontally with respect to the mesh, i.e. a stencil for which the height is calculated along the x-axis of the mesh.

Offset procedure

The interface extends over several cells. Each cell for which the interface density is not zero must be taken into account. In other words, as long as ∇γ ̸ = 0, a source term is calculated for the cell. Since the interface is diffuse in the VOF method, the volume fraction gradient represents the portion of the interface present in the cell. In cells where the volume fraction is equal to 0 or 1 and its gradient is not zero, the calculation of the height does not fall within the cell concerned. For example, in the diagram in Fig. 19 the gradient of the volume fraction of the cells in the top row is different from zero, so the source term needs to be calculated. However, the volume fraction is zero. When the height is calculated the value is outside the cell concerned. The calculation of the curvature is then false. To overcome this defect it is necessary to recover the calculation of the curvature of another nearby cell. An offset procedure originally thought up by Magnini can be used to recover this value in one of the cells of the stencil column [Magnini, 2016b]. However this target cell sometimes can be 4

Algorithm 1: Calculation of the curvature of the cells with α = 0 or 1 Result: return κ (i,j) num (i,j) = 0; den (i,j) = 0; if ∇ α (i,j) > 0 then if α (i,j) = 0 or 1 then for k1 = -1 to 1 do for k2 = -1 to 1 do if (k1 ̸ = 0 and k2 ̸ = 0) and ∇α (i+k1,j+k2) > 0 and (α (i+k1,j+k2) ̸ = 0 or 1) then num (i,j) = num (i,j) + κ (i+k1,j+k2) × ∇α (i+k1,j+k2) ; den (i,j) = den (i,j) + ∇α (i+k1,j+k2) ; end end end end end

κ (i,j) = num (i,j)
den (i,j) ; cells away from the original cell (usually when the interface is diagonal to the mesh). Results can be misleading when testing the code on a stagnant bubble in a zero-g environment. The curvature on a circle is constant. So when the algorithm retrieves data from a remote cell, it looks better. But this does not work for an arbitrary interface. This is why another procedure averaging the curvature weighted by the gradient of the volume fraction of the cells adjacent to the original cell has been used, as shown in Alg. 1.

Orientation of the local coordinate system

As mentioned above, the HF method loses accuracy when the interface slope is too steep in relation to the local coordinate system. The approximation of the first and second derivatives becomes weaker. This is why, it is necessary to use two types of local coordinate systems, one horizontal and the other vertical, and whose axes are parallel to those of the mesh, and then to switch from one to the other as needed. According to the mesh coordinate system, the value of the components of the volume fraction gradient allows to estimate a first orientation of the interface. This orientation estimation allows the selection of the most favourable stencil for the calculation of the heights. In other words, the one in which the interface will be oriented most parallel to the abscissa axis of the local coordinate system. However, when the interface is diagonal to the mesh axes, another method must be used.

As shown in Fig. 23 one solution is to used a mesh decoupled method and computes heights within columns not aligned with the computational mesh but rather aligned with the interface normal vector [START_REF] Ito | A high-precision calculation method for interface normal and curvature on an unstructured grid[END_REF][START_REF] Owkes | A mesh-decoupled height function method for computing interface curvature[END_REF]]. The problem is that these methods need a first computation of the interface normal, and thus of the volume fraction gradient. Another point to emphasize is that these methods cut the cells of the mesh to compute the integrals. Data from several columns of the mesh are used to reconstruct the integral under the curve. Thus the equality relationship between the volume fraction and the integral demonstrated in Fig. 19 under the interface is lost. The heights are approximated by relying on the volume fraction gradient. The approximation of the volume fraction gradient is the first source of error in the calculation of the curvature in more traditional methods such as CSF. The advantage of the HF method is to be able to calculate the curvature without using this gradient. The disadvantage is that its usage is restricted to rectangular meshes.

Selection of the heights

In the previous step, heights were calculated for the horizontal and vertical stencil. These heights will define points in a local coordinate system to be determined in which the interface will be approximated by curvature of a second degree polynomial. It is not clear which method is preferable between using height functions directly to determine the curvature and using a polynomial fit. The answer to this question will be discussed in a later chapter, as it requires a special numerical study. It should be noted that the polynomial fit technique is more complex and computationally intensive than the one based on height functions. In a way, the heights are used to reconstruct the interface. The diagram in Fig. 24 shows a representation of the interface. In the two vertical and horizontal stencils a total of six heights were calculated, each stencil having 3 columns. However, four heights were retained, i.e. four points representing the interface. One of the tools to help select these points is the normal to the interface n. This was calculated previously using the gradient of the volume fraction. Stencil columns in which the interface is too sharp are excluded. An orthonormal system whose origin O is the median of the two points closest to the centre of the cell and whose ordinate is directed along the normal to the interface (O, i ′ , n). Then the coordinates of each point P are computed in the new sytem x m , y m .

Fit of a parabola

The parameters (a 0 , a 1 , a 2 ) of the equation of a parabola are estimated from the minimisation of an objective function:

f pol (a i , x m ) = a 0 x 2 m + a 1 x m + a 2 (3.62) F objective = m [y m -f pol (a i , x m )] 2 (3.63)
The curvature is estimated from these estimated parameters:

κ = 2 a 0 (1 + a 2 1 ) 3/2 (3.64)
As pointed out by Popinet, while the least-square minimisation is not particularly complex or computationally expensive, the difficulty is to select the good points [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF]. In a more general way, the main difficulty of an HF algorithm is for the programmer to adapt to the many particular cases, which can make the code with many conditional statement to implement.

Identification of parameters

Two types of methodologies were used to identify the parameters. One based on the direct minimisation of Eq. (3.63), and the other based on the Levenberg-Marquardt algorithm. The second method is more flexible and can be adapted to both linear and non-linear problems. The coding of this second methodology was done firstly with the intention of validating the two methodologies by comparing their results, and secondly with the objective of making a tool for minimising quantities for testing purposes.

The objective function reaches a minimum when F ′ objective = 0. This corresponds to solving the following system of equations:

m y = a 2 m + a 1 m x + a 0 m x 2 (3.65) m xy = a 2 m x + a 1 m x 2 + a 0 m x 3 (3.66) m x 2 y = a 2 m x 2 + a 1 m x 3 + a 0 m x 4 (3.67)
Where m is the number of coordinates considered. The methodology of Levenberg-Marquardt is shown in Alg. 2.

Algorithm 2: Determination of parameters -Levenberg-Marquardt d LM is the direction of descent of the algorithm J is the Jacobian matrix Result: return crit quad n iter = 0; while n iter < cste do crit quad = 0

for i = 0; i < m; i + 1 do dev[i] = y m [i] -f poly (x m , a 0 , a 1 , a 2 ) crit quad = dev 2 [i] + crit quad end if n iter == 0 then λ = cste; λ rec = λ crit quad,rec = crit quad a 2,rec = a 2 , a 1,rec = a 1 , a 0,rec = a 0 else if crit quad < crit quad,rec then λ = 0.1λ; λ rec = λ crit quad,rec = crit quad a 2,rec = a 2 , a 1,rec = a 1 , a 0,rec = a 0 else λ = 10λ a 2 = a 2,rec , a 1 = a 1,rec , a 0 = a 0,rec end end d LM = (J t J + λΩ) -1 (J t [y -f poly (x)]) for i = 0; i < m; i + 1 do a[i] = a r ec[i] + d LM [i] crit quad = dev 2 [i] + crit quad end n iter = n iter + 1 end 100 Chapter 3
Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022

Insertion in the Navier-Stokes equations

Averaging of physical properties at the interface

The physical properties of each fluids are calculated as weighted averages based on the distribution of a phase volume fraction, thus being equal to the properties of each fluid in their corresponding occcupied regions and varying only across the interface. Single-field quantity are defined all over the computational domain. In the case of a system with two fluids gas and liquid, in each cell we get :

β = α g β g + α l β l (3.68) β = (1 -α l )β g + α l β l (3.69)
To simplify the reasoning a single variable α is considered such that :

α = α l = 1 -α g (3.70)
Thus density is considered a function of volume fraction::

ρ(α) = (ρ l -ρ g ) α + ρ g (3.71)
However, in order to correctly capture the viscosity term at the interface µ ∇v + (∇v) T , an average of the dynamic viscosity µ performed as for the density can be approximative. At the interface, gross and abrupt changes in viscosity take place. An accurate evaluation of µ is crucial to reproduce the correct free surface. The use of an arithmetic mean actually causes an artificial acceleration of the fluid in the less dense phase, resulting in speeds that are too high due to a non-physical viscous term. According to Kothe et al., the relative interface/cell face orientation must be taken into account [START_REF] Kothe | Perspective on Eulerian finite volume methods for incompressible interfacial flows[END_REF].

µ = η [(1 -α)µ g + αµ l ] + (1 -η) 1 -α µ g + α µ l -1 η = |n • n F | (3.72)
where n F is a unit vector normal to the face cell. The formulations of the averaged quantities vary from one model to the other to ensure code stability and reduce discretization errors.

Density shift procedure

In the original CSF method the final momentum source is defined as shown in the equation below

f γ ′ = ρ 1 2 (ρ g + ρ l ) f γ (3.73)
Surface tension should act uniformly, regardless of the instantaneous liquid fraction and density in the cell. It may be argued that the force should act on the mass in a shifted cell, which is ∆x/2 thick on both sides of the interface. In this cell, the density is constant and equal to the mean density ρ = 1 2 (ρ g + ρ l ). By dividing by ρ rather than ρ (i,j) , spurious currents can indeed be reduced.

Guo et al. used modified a density shifting procedure to improve the stability of the momentum sources [START_REF] Guo | Implementation of a height function method to alleviate spurious currents in CFD modelling of annular flow in microchannels[END_REF]:

f (γ,i,s) ′ = N (i,s) ρf (γ,i,s) (3.74) N (i,s) = V cell f (γ,i,s) V cell ρf (γ,i,s) (3.75)
where i = x, y denotes the coordinates, s = ps, ns denotes the positive and negative sources. This shifting procedure causes the surface tension force to be applied more to the phase with the higher density and therefore improves the numerical stability. ρ is the bulk density and N (i,s) is is a normalisation factor used to conserve the source. The shifting for positive and negative sources needs to be kept separate otherwise a negative value of the normalisation factor may occur and therefore cause the signs of the local source values to be reversed.

The problem with Eq. (3.73) is that it may not strictly conserve the source in the computational domain although it can improve the numerical stability compared to the direct implementation of f γ . However, by using the normalisation factor, the density shifting procedure shown by Eq. (3.74) achieves accurate source conservation.

Marangoni model

Having described in the previous section how the pressure jump at the interface was calculated, it is appropriate to describe how ∂γ ∂s t, i.e. the Marangoni stress, can be calculated.

The methodology presented here is based on height functions. The calculations use the data collected by the part of the algorithm presented in section 5.2.

Deficiency of the surface gradient operator

It is underlined by Seric et al. that using the surface gradient operator as it is defined in Eq. (3.48) can result in inaccuracies when implemented in the VOF method [START_REF] Seric | Direct numerical simulation of variable surface tension flows using a Volume-of-Fluid method[END_REF].

This definition of the surface gradient can result in inaccuracies when implemented in the VOF method for general variable surface tension for two reasons. First, the discontinuities of the material properties across the interface represented by Eq. (3.1) can result in Q having a large jump across the interface: for example, in the case of surface tension dependence on the temperature where the fluids on each side of the interface have large difference in the conductivity. The second reason is that, in general, surface tension can depend on the concentration: for example, in the case of the mixing of two liquids with different surface tension, or in the case of surface tension dependent on the surfactant concentration.

The Marangoni effect is caused by a tangential stress located on the interface due to a variation in surface tension. Surface tension is a concept that only makes physical sense at the interface and its variation only makes sense along that interface. An estimation of the surface tension gradient is necessary to be able to calculate the surface gradient operator. The surface tension is a function of, among other things, the temperature and the concentration of chemical species. However, these quantities generally vary abruptly when they are considered on either side of the interface.

Therefore in VOF and similar methods the gradient of surface tension is a vector whose direction is close to the normal of the interface, which makes no physical sense. The surface gradient operator subtracts its interface normal component to recover its tangential component. This approximation constructed from data from few cells (depending on how the gradient operator is calculated) is very sensitive to errors.

Calculation Marangoni term with HF

Seric et al. propose to implement the variation of surface tension using a method inspired by height functions [START_REF] Seric | Direct numerical simulation of variable surface tension flows using a Volume-of-Fluid method[END_REF]. The algorithm for implementing ∇ s γ(x) in the VOF method starts with the approximation of the interfacial values of the surface tension in each cell containing an interface segment. More precisely, the idea of constructing the columns of cells inspired by the computation of interfacial curvature and normals using height functions is used.

In this method the derivative along the interface of the surface tension is calculated directly:

f st = ∂γ ∂s δ s t (3.76)
where s is the arc length. The derivative ∂γ ∂s , δ s , and t are evaluated using the cell-centers values. Firstly the surface tension values are defined at the interface, then the derivatives of γ are computed along the interface.

In this method the tangential component at the interface of the surface tension gradient is directly obtained which avoids the projection along the tangent to the interface thus avoiding the calculation errors generated by the surface gradient operator. Its other advantage is that the diffuse data due to the VOF method are averaged in a direction close to the normal to the interface, the columns of the stencil being oriented in the direction most perpendicular to the interface.

Calculation of the average surface tension coefficient in each cell

The first step is to determine the surface tension evaluated from the concentration at the center of all interfacial cells γ(pos, c h 2 ), with the volume fraction α pos . Where pos is the the coordinate of the targetted cell center. The surface tension in each column, denoted by γ x (pos x , c h 2 ), is defined so that it has only one value in each column, regardless of how many interfacial cells are contained in that column. The superscripts x, y represent the column direction.

For columns with only one interfacial cell , the surface tension of the interfacial cells is not averaged among the neighbouring cells. The value is directly taken from the surface tension of the cell.

If there is more than one interfacial cell in the column, then γ x (pos x , c h 2 ) is approximated by the volume weighted average of the suface tension values belonging to the same column.

γx i (pos x , c h 2 ) = α i,j γ i,j + α i,j+1 γ i,j+1 + • • • j α i,j
(3.77)

γy j (pos y , c h 2 ) = α i,j γ i,j + α i+1,j γ i+1,j + • • • i α i,j (3.78)
In this implementation γx i (pos x , c h 2 ) is first defined for all interfacials cells. For certain coordinate , it is posible to define γ for more than just one interfacial cell.

The direction of the column as to be determined along whit the computation of the surface forces. The choice of the direction is based on the interface orientation: x or y is chosen to be the same as the largest component of the normal vector to the interface. The same choice is made for computing curvature and the interface normal using height functions.

Derivative of surface tension

The next step is to evaluate the derivative along the interface ∂γ ∂s .

The derivative of the surface tension along the interface, ∂γ ∂s , is approximated by the derivative of the interfacial value, γ in the column which is formed in the direction x or y.

The derivative is computed as :

∂γ ∂s i,j = γ(i + 1) -γ(i -1) ds (3.79)
γ(i) is the weighted average of the values taken by the surface tension for the cells of column i. Therefore the same value of surface tension will be used in the calculation of the derivative for the cells of the same column of the stencil.

γ(i) = j=+∞ j=-∞ ∇α ij × γ i,j j=+∞ j=-∞ ∇α ij (3.80)
In each interfacial cell, the derivative is computed along the interface using central difference, i.e. the finite difference of the ~γ in the two neighboring columns.

∂γ x ∂s i,j = γj+1 -γj-1 ds (3.81)
γj is the interfacial value of the surface tension in the column j constructed in the x direction.

Derivative of the arc length

Then the derivative of the arc length is calculated from the derivative of the height function:

ds = 2∆ x 1 + H ′ (i) (3.82)
The arc length ds is computed from height function in the same direction as ∂γx ∂s . For the example given previous equation , the arc length is :

ds = 2∆ y 1 + H ′ (j) (3.83)
where H ′ (j) is the derivative of the height function and ∆ y is the cell size. The next part of the surface gradient implementation is the choice of the tangent vector, t, which is computed so that it satisfies t • n = 0. The direction of t depends on the direction used for computing ∂γ ∂s . points in the direction of the positive component orthogonal to the x, y direction. As an example, t points in the positive x direction if column is built in the y direction.

An intermediate value as to be considerate when choosing the tangent vector.

M G x = ∂γ ∂s sgn(t x ) (3.84) M G y = ∂γ ∂s sgn(t y ) (3.85)
The Fig. 25 show an examples ∂γ ∂s , M G x and M G y respectively, are computed in all interfacial cells, where a positive uniform gradient of the surface tension is imposed in the y direction. On the left figure ∂γ ∂s changes sign in the first and third quadrant at the angles, defined from the positive x axis, of π/4 and 5π/4, respectively. At these points the direction of the columns used in gradient computation changes. Hence, the two neighboring cells have opposite sign of ∂γ ∂s . However, once the correct sign of the tangent vector components is included and each component is considered separately, as in the two previous equations , this inconsistency in the sign is corrected.

The numerical expression of the component of the surface force are then :

F st,x = M G x |t x |δ s (3.86) F st,y = M G y |t y |δ s (3.87)
where δ s is the dirac function defined previously , which represent from the numerical point of view the part of the cell present in the interface The two values M G x and M G y have to be evaluated in all the cells where δ ̸ = 0 . The same approach as the one used to define the curvature of the cells neighboring the interface is used. The values in the cells neighboring the interfacial cells are defined by averaging the values in the direct neighbors that already have the curvature value defined. This procedure is repeated twice, insuring that the curvature values for the corner neighbors to the interfacial cells are defined as well.

We use an identical approach for defining the x and y components of M G in the cells around the interface which are subsequently used in the two previous equation.

Numerical model mass transfer across the interface

Firstly, this section presents how interfacial mass transfer is modelled and why it is necessary to perform an accurate calculation of the interfacial area to ensure mass conservation on both sides of the interface. Then in a second step the numerical method used to calculate this area is presented.

Conservation of dissolved species

Capturing the fluid behavior of a multiphase and multicomponent system composed of a gas phase and a liquid phase with the gas diluted in the liquid requires the modeling of its multiphase and multicomponent characteristics. From a numerical point of view, the multiphase behavior can be captured by solving the VOF equations and the multicomponent aspect by solving the species transport equations. However, certain considerations must be taken into account in order to develop a numerical and holistic model simulating simultaneously multiphase and multicomponent physics.

Firstly the dissolved species must be tracked only within the phase it exists in, i.e. the solution to the species transport equation for dissolved concentration must only have non-zero values within the liquid phase and not the gas phase. Otherwise, having a non-zero solution in the gas phase would imply the existence of dissolved gas species within the gas phase itself, which is non-physical.

Secondly the mass must be preserved for the species studied during mass transfer. For example, to capture the absorption of dissolved species in the gas phase, a source/sink term S α,sp-trsf t should be added to the VOF equation, and a source/sink term S Y,sp-trsf t should be added to the species transport equation:

∂α(t) ∂t + ∇ • (α(t) v) = S α,sp-trsf t (3.88) ∂ρY i ∂t + ∇ • ρ (vY i -D i ∇Y i ) = S Y,sp-trsf t (3.89)
where Y i is the mass fraction of species i, D i is the diffusion coefficient of species i.

The mass transfer rate must be determined by the interface jump conditions so that there is a chemical equilibrium between the species.

In the bulk of each phase, these source terms must be equal to 0, and must have an influence only at the interface. These source terms are volumetric and from a numerical point of view will be applied to the centre of each cell of the mesh. However, interfacial mass transfer is a phenomenon that occurs on a surface and is localised at the interface. It is therefore necessary to ensure from a numerical point of view that the species are well conserved.

In a multiphase, multi-component system the chemical species will move to bring the system to equilibrium. The state of non-equilibrium is at the origin of mass transfer. It is the mass transfer that pushes the system towards equilibrium. In many numerical simulations the rate required for the system to reach equilibrium is used as the basis for formulating the rate of mass transfer across the interface. The objective of the mass transfer is to calculate the coefficient ṁ of the mass transfer flow rate.

Interfacial area calculation

For the mass transfer rate to be calculated accurately A I must be evaluated correctly [Soh et al., 2016;[START_REF] Schlottke | Direct numerical simulation of evaporating droplets[END_REF]. The gradient of α gives a biased representation of the interface. At the local level, |∇α| can have non-zero values in cells of the mesh where in the continuous model the interface is not present. These cells are adjacent to the interface

Algorithm 3: Calculation of A I in 2D Result: return A I n iter = 0 A I = 0 step x = |x 2 -x 1 | n L x a = 0, x b = x a + step x while n iter < n L do L = 0 y a = f pol (x a ) y b = f pol (x b ) L = (x b -x a ) 2 + (y b -y a ) 2 A I = A I + L x a = x a + step x x b = x b + step x n iter = n iter + 1 end
cells, |∇α| is computed using the α values of the neighboring cells, including the interface cells. Thus |∇α| may not be zero for cells where α = 0 and α = 1. These cells can generate an artificial mass transfer and after the transfer phase the calculation can result in a value of α negative or greater than unity.

With the developments made previously on the height functions it is possible to calculate A I without using the gradient of the volume fraction. For a 2D simulation, by determining the coordinates of the points where the interface crosses the axes of the mesh and knowing f pol the length of the interface can be determined by using the alg. 3 allowing estimation of the length of the curve of a function. The prerequisite is to know the intersection points of the parabola with the cell boundaries named here arbitrarily P a (x a , y a ) and P b (x b , y b ). Meier et al. used a similar method [START_REF] Meier | A novel technique for including surface tension in PLIC-VOF methods[END_REF]. The calculation of A I can be extended to a 3D simulation [Soh et al., 2016]. The use of such a computational method is necessary to be able to integrate a mass transfer model across the interface [START_REF] Soh | A CFD model for the coupling of multiphase, multicomponent and mass transfer physics for micro-scale simulations[END_REF].

Integration of the code in Fluent

This section presents how the numerical methods described above have been integrated into the Ansys Fluent code.

User Defined Function

The commercial software Ansys Fluent 2020 R2 was used to perform the numerical simulations. The model was coded and implemented using User-Defined Functions (UDFs).

The coding was done in C language. The code was inserted within a DEFINE ADJUST macro. This macro is executed by Fluent just before solving the mass, momentum, volume fraction and species transport equations. The main function of the macro code is to calculate the source terms that can be added to the conservation and transport equations. Fluent allows its users to calculate the surface tension terms with different types of methodology through a graphical interface and in particular with the CSF method. For comparison and testing purposes this module was used, but disabled in the general case so as not to interfere with the calculations done by the DEFINE ADJUST macro.

The code has been adapted to parallel calculation, in order to reduce the calculation time. One of the main difficulties of the coding was to adapt the custom-made calculations to Fluent's calculation system. Fluent offers a turnkey CFD solution, which means that the user can customise the software to his needs but is not expected to replace an entire computing system. In order to identify each cell Fluent uses its own referencing system. The calculations performed on the height functions require an exact knowledge of the position of the coordinates of each cell. This is why one of the first functions of the DEFINE ADJUST is to establish its own coordinate and referencing system different from that of fluent, and to establish the correspondence between the two systems if necessary. This explains the use of a particular mesh. The DEFINE ADJUST code can only work with this type of mesh. The mesh consists of square and orthonormal cells as shown on fig. 26. Different mesh sizes were used for the simulation: 40 × 40 , 80 × 80, 120 × 120, 160 × 160. Increasing the mesh size increases the computation time. To make the calculation of surface tension forces more accurate at the interface it is interesting to refine the mesh size locally. Fluent provides macros for local refinement of the mesh size. However, this has the effect of changing the cell referencing system, which makes it difficult to transfer the information to the DEFINE ADJUST macro referencing system. From a coding point of view, the simplest solution is to globally refine the mesh to obtain the desired level of accuracy of the surface tension force calculation.

Discretisation

The gradients of scalars are calculated as cell centroid values from the centroid values of faces surrounding the cell. The Green-Gauss node-based method is used for this calculation. The PRESTO scheme is used for pressure interpolation. The QUICK scheme is used for the discretisation of the momentum and the energy equations. The Piecewise-Linear Interface Calculation (PLIC) scheme is used for the discretisation of the volume fraction equation. When simulations are made using the CSF method, the default node based smoothing of the volume fraction field prior to calculation of the curvature was enabled. No smoothing of the calculated curvatures was performed. A first order implicit scheme is used for the temporal discretisation of the transient terms. Finally, for the pressure-velocity coupling, the SIMPLE algorithm is used.

Presentation of the code used in the simulations

The numerical methods presented above are part of the puzzle pieces that were used to code the DEFINE ADJUST macro used in the simulations presented in this thesis. The diagram in Fig. 27 describes how the DEFINE ADJUST macro is integrated into the Fluent algorithm. The source term f γ Eq. (3.90) is calculated at each iteration before the Navier-stokes equations are solved. Suitability of the VOF approach to model electrogenerated bubble with Marangoni micro-convection flow 4

In order to disentangle the effects of spurious currents from Marangoni currents, the overall approach presented in this chapter is to evaluate spurious currents when there is no simulated Marangoni effect. First, the error in the curvature calculation is evaluated. Then tests on bubbles in stagnant fluids are performed. These allow the evaluation of spurious currents as theoretically no current should be generated in the case of a stagnant bubble. These tests are performed for different mesh resolutions, the errors generated being dependent on the mesh resolution. The objective of this chapter is to evaluate the errors generated by the algorithm presented in the previous chapter, to compare it with the CSF methodology, and to evaluate its suitability to model electrogenerated bubble with Marangoni micro-convection flow.

1 Curvature calculation errors

Choice of curvature calculation method

It has been shown by Cummins et al. that starting from an exact volume fraction value, calculations with the standard method of height functions estimate the curvature asymptotically with second order accuracy [START_REF] Cummins | Estimating curvature from volume fractions[END_REF]. But, as previously mentioned, the height function method loses its effectiveness when the interface approaches the diagonal of the mesh axes, in which case a polynomial fit method is more appropriate.

However, from a theoretical point of view there is no obvious way to decide whether to switch from one of these methods to the other. An orientation angle θ switch of the interface must be determined to know when to switch from one method to another. Numerical tests are needed. A first step is to evaluate the performance of each method. Then by comparing the errors of each method the most suitable angle θ switch can be chosen.

Study parameter and evaluation criteria

The first parameter to take into account is the size of the mesh, several tests have been performed for different spatial resolutions. The curvature calculation tests were performed on circular interfaces. The ratio 1 κ ∆ is used to compare the mesh size to the curvature. Another criterion to consider is the theoretical position of the interface within the cell. Depending on this position the value of the volume fraction is modified. This can have an influence on the curvature results provided.

To prevent the curvature calculation from being interfered with by errors in the numerical calculation of the volume fraction, the curvatures have to be evaluated from analytically calculated exact volume fractions, so that only the curvature calculation method is evaluated. In this study the interface is made up of arcs. The exact volume fraction can be determined from the integral of the function representing these arcs of circles. The circles have as their centre the origin of the reference frame for determining the coordinates of the mesh and have the function: x 2 + y 2 = R 2 . In this trivial case, the exact curvature is the inverse of the radius of the circle. The evaluated curvature is compared with the exact curvature:

∆κ error = 1 κ exact |κ -κ exact | (4.1)
The quantity ∆κ error is used to determine the accuracy of the numerical algorithm for calculating the curvature.

Calculation of the exact volume fraction

As shown in Fig. 29, each cell of the mesh is used as a basis for determining a local coordinate system whose centre on the diagram has the coordinates (x 0 , y 0 ).This gives the circle equation in the local system: In order to calculate the integral of the function, the coordinates are expressed in an explicit form:

(x + x 0 ) 2 + (y + y 0 ) 2 = R 2 (4.2)
y = R 2 -(x + x 0 ) 2 -y 0 (4.3) x = R 2 -(y + y 0 ) 2 -x 0 (4.4)
To simplify the following calculations only the cases where x > 0 and y > 0 are considered. The calculations and reasoning remain equivalent to the nearest sign in the other cases. As shown in Figure 1, the arc of the circle intersects the y-axis at m y and the x-axis at m x in the local coordinate system:

m x = R 2 -y 2 0 -x 0 (4.5) m y = R 2 -x 2 0 -y 0 (4.6)
In Fig. 29 the surface A c , allowing to determine the volume fraction, is the intersection of the surface of the cell of the mesh with the surface represented by the integral of the curve and satisfying the equation:

A c = A 0 -A 1 -A 2 .
The integral under the curve A 0 is cut into three areas, in order to determine the area representing the volume fraction. By determining A 0 , and the two surfaces A 1 and A 2 , the surface A c is obtained:

A 0 = mx 0 y dx (4.7) A 1 =      0 if |m x | < |∆ x | mx ∆x ydx |m x | > |∆ x | (4.8) A 2 =        0 if |m y | < |∆ y | my ∆y xdy |m y | > |∆ y | (4.9)
In the case where |m x | < |∆ x | or |m y | < |∆ y |, the surfaces A 1 and A 2 do not exist. This is why it is necessary to include a condition in Eq. (4.8) and Eq. (4.9). From these areas, A c can be determined, which makes it possible to find the volume fraction by taking up the definition in Eq. (3.37):

α = A0 -A 1 -A 2 ∆ x ∆ y = A c ∆ x ∆ y (4.10)

Determination of θ switch

Beyond a certain angle of inclination θ the height function method becomes less efficient. The polynomial fit method is more expensive to compute. A system of equations must be solved for each cell of the mesh. There is the computational cost and precision to consider. An angle to define when to switch from using one methodology to another to calculate curvature must be determined. A suitable angle value that assures that the results will not diverge has to be found. The choice of this angle is more a question of safety than of optimization.

Curvatures are evaluated for interfaces with a tangent inclined at an angle θ of 0 to 45°with respect to the horizontal axis of the mesh. The usual height function method as mentioned before presents good results for interfaces whose inclination is close to the axes of the mesh. However, from a certain inclination, the results obtained diverge from the real value of the curvature, as shown in Fig. 30. As the inclination of the interface increases, the error increases. The choice of the alternative method of polynomial fit is to be considered. In view of the results obtained, the transition from one method to the other is to be considered for an interface tilt around θ switch = 22.5°. The main difficulty of the polynomial fit method is to choose the interface points to use. When θ is less than θ switch , the height function method is used. The transition to the polynomial fit method is made for θ greater than θ switch . With regard to the results presented in Fig. 31 for resolutions where the radius of the circular interface considered has a length equivalent to 5 cell widths of the mesh, it appears that the use of the two techniques used is inconsistent. The points approximating the position of the interface are not close enough to be able to correctly estimate the polynomial parameters.

On the other hand, when the spatial resolution becomes finer and for θ greater than θ switch the fitting method gives better results.

These tests on the curvature allowed us to establish suitable value for θ switch . The center of the circular interface was moved to test the robustness of the methodology. In general this test has no influence on the results obtained. However, in cases where the part of the interface present in the cell is too small, the calculation error on the curvature diverges. The choice of the weighted average calculation of the curvatures on the adjacent cells was considered instead. This choice allows much better results to be obtained. The electrogenerated bubbles having almost circular interfaces so the choice of an average seems coherent. Generally speaking, for resolutions for which the interface radius is equivalent to 15 times the width of the mesh, calculation errors of less than 0.3% are obtained. This preliminary test allows us to be confident about the curvature calculation methodology.

Static bubble test case

Pressure jump at the interface

The analytical solution for the simulation of a stationary bubble in a zero velocity field, and the analytical curvature can be easily obtained from the bubble radius. A circular interface with surface tension should remain at rest, with the pressure jump at the interface exactly balancing the surface tension force (Laplace's law). The velocity field being zero, Eq. (3.28) reduces to:

-∇ • (pI) = 0 (4.11)
In this test case within each fluid the pressure is constant. The jump relation Eq. (3.31) which is applicable only at the interface reduces to the mathematically exact formulation:

pI • n I + γκn I = 0 (4.12)
This brings us back to the relation of Laplace. In each phase the pressure is constant and a pressure jump occurs at the interface. As shown in Fig. 32 the pressure jump created by the CSF model at the interface is less direct, which deviates from the real conditions, while the height function methodology gives a better approximation.

Time to reach equilibrium

In practice, depending on the method used to discretize the pressure gradient and the surface tension force, parasitic currents appear. The exact numerical balance is difficult to obtain [START_REF] Abu-Al-Saud | A conservative and well-balanced surface tension model[END_REF]. The numerical imbalance created is at the origin of these currents.

Similar to what was done by Popinet, it is appropriate to first test the model by imposing the exact curvature in the entire domain for the calculation of f γ [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF]. This tests the adequacy of the model by excluding the curvature calculation, and thus verifies that the balance calculation between the pressure term and the surface tension term is indeed achieved. The time required for the momentum to diffuse over a distance L is proportional to t ν , where

t ν ∝ L 2 ν (4.13)
and ν is the kinematic viscosity of the liquid. As noticed by Popinet, the time scale needed to reach the numerical equilibrium solution is comparable with the time scale of viscous dissipation t ν , as expected from physical considerations [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF]. In practice, this means that test cases designed to evaluate the accuracy of a given surface tension model (for a stagnant bubble equivalent problem) must ensure that simulations are run for time scales comparable with t ν . In our case t ν is close to 2 ms. The other quantity to consider is the velocity associated with the capillary wave u γ .

u γ ∝ γ ρ L (4.14)
It can be interpreted as the scale of the velocities associated with a capillary wave of wavelength comparable with L. As shown in Fig. 33, the average velocity obtained decays for a time equivalent to the viscous dissipation time. The velocity and time have been scaled using u γ and t ν . Thus, the numerical calculation verifies the theoretical equilibrium and the spurious currents observed in the following can be attributed to errors in the curvature calculation.

Errors caused by spurious currents

In this second test the curvature is calculated by the model. In order to evaluate the impact that spurious currents could have on a simulation with Marangoni effect, the maximum speed of the spurious currents u max,spurious currents obtained are scaled using an average speed u average, M arangoni . Several simulations were performed for different mesh resolutions for each of the two tested models: CSF and and HF the algorithm based on height functions presented in section 3.9. The simulations were carried out for a time equivalent to the experimentally observed growth time of a bubble. As shown in Fig. 34 and in Fig. 35, for the CSF method the spurious currents increase when ∆x decreases . This is consistent with the analysis presented by Harvie et al. [START_REF] Harvie | An analysis of parasitic current generation in volume of fluid simulations[END_REF]. The CSF method is therefore clearly not suitable for the simulation that is the objective of this study. This validates the use of a more efficient interface representation method. The error generated using HF decreases with the reduction of the grid spacing. The results obtained in this section show that the method used is balanced and allows estimation of the curvature sufficiently accurate to obtain a solution close to the exact equilibrium (for the velocity). The numerical equilibrium obtained is very close to the theoretical equilibrium. Even for coarse resolutions the error generated on the final simulation is less than 1 %.

Interfacial area error calculation

In order to evaluate the accuracy of the calculation of the interfacial area of each cell of the mesh, several simulations were conducted. The calculated value is compared to the exact value of the interface A I,exact . The test case of a static bubble as described below was used. The interface being circular the exact value of the interfacial area in each cell can be calculated analytically. As the position of the interface in each cell can influence the calculated numerical value, the values of all cells through which the interface passes were averaged, as described by the following equation:

E(A I ) = N |A I,exact -A I | N (4.16)
where N is the number of cells for which the calculation was performed. The results are shown in the graph in Figure 36. The calculations were performed as a function of the parameter n L , which determines the accuracy of Alg.3, and as a function of the ratio R/∆ x , which determines the accuracy of the mesh. The graph in Fig. 37 shows the maximum error found, as described by the following equation:

E max (A I ) = max|A I,exact -A I | (4.17)
The error decreases with a finer mesh, and by increasing the value of the parameter n L .

Surface gradient error calculation

Next, the efficiency of the surface tension gradient calculation should be tested, as shown in Eq. (3.79). The static bubble is subjected to different temperature gradients over a given distance as shown in Fig. 38. The objective here is to expose the interface of the bubble to variations in surface tension similar to what it might encounter as it grows, the bubble is exposed to surface tension variations ranging from 0.1N • m -2 to 50N • m -2 . As the interface is circular the exact value of the surface tension gradient can be calculated.

For each cell crossed by the interface the length of the interface is known and every two cells the temperature difference can be calculated. For each simulation the surface tension gradient is averaged along the interface in order to compensate for uncertainties concerning the influence of the position of the interface within the cell. The relative error found is calculated according to the following equation: where N is the number of cells used for the calculation. The maximum errors found are also recorded, and calculated using the equation:

E(∇ s γ) = N |∇ s γ exact -∇ s γ|/|∇ s γ exact | N (4.18)
E max (∇ s γ) = max|∇ s γ exact -∇ s γ| |∇ s γ exact | (4.19)
The errors found are counted and plotted in Fig. 39. The finer the mesh, the smaller the error.

As shown in the graph in Fig. 39, the calculated errors are not very sensitive to the value of the surface tension gradient, but they decrease rapidly as the mesh is refined. The maximum errors are mainly due to cells where the interface share is small compared to the total cell volume.

Translating bubble test case

While the case of a stagnant bubble allows us to test the equilibrium of the model by referring to an exact solution of the velocity field, it does not allow us to evaluate the combined accuracy of the interface advection and surface tension model. As proposed by Popinet, the horizontal translation of a bubble carried by a uniform flow field is a more robust and realistic test [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF]. In the case of electrogenerated bubbles, when the bubble grows, the interface translates at a vertical speed of a few millimeters per second. This is a preliminary test before testing mass transfer models across the interface. A uniform horizontal velocity u 0 is imposed in the whole domain with periodic boundary conditions on lateral sides and symmetry boundary conditions on the top and bottom. As already reported by Popinet, the absolute error on the velocity does not depend on u 0 and is weakly dependent on the Laplace number La = ρLγ µ . It is thus the transport scheme of the interface that is directly tested and therefore the impact of the spatial resolution. In our study the Laplace number varies between 5, 000 and 20, 000. A new time scale is introduced, to account for the time needed for the bubble to cross a length L, and is given by:

t u 0 ∝ L u 0 (4.20)
The velocity has been scaled with u M and the time with t u 0 . The Laplace number was fixed at 12, 000. The results presented here as an example reflect a general trend in the evolution of the spurious velocity over time as observed by Abadie et al. [START_REF] Abadie | On the combined effects of surface tension force calculation and interface advection on spurious currents within Volume of Fluid and Level Set frameworks[END_REF]. Popinet notes that these oscillations are proportional to u 0 /∆x. The advection errors of the models continuously disturb the calculations related to the interface [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF].

As shown in Fig. 40, the model has similar behavior to the previous studies [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF][START_REF] Abadie | On the combined effects of surface tension force calculation and interface advection on spurious currents within Volume of Fluid and Level Set frameworks[END_REF]. The drastic drop in performance can be noticed. Even if the height function method allows accurate curvature calculations, the flaw of the method comes essentially from the advection method. The previous simulation was repeated for different mesh resolutions. The dimensionless quantity used in the abscissa is R/∆x. As for the static case the method converges when the mesh is refined. We obtain an error of 2.5% for R/∆x = 70.

Conclusion of the chapter

As shown in the example in Fig. 41, the use of the height function method greatly reduces the error rate due to spurious currents. The objective of this simulation is to visually represent the impact that spurious currents could have on a solutal Marangoni type effect on an electrogenerated bubble attached to an electrode. The interface is fixed, the contact line model is not implemented and so the bubble is not in contact with any edge. To replicate what might exist in an electrode-generated bubble, a gradient of dissolved species is initiated along the interface, where the Marangoni effect takes place as shown in Fig. 42.

The gradients of scalars are calculated as cell centroid values from the centroid values of faces surrounding the cell. The Green-Gauss node-based method is used for this calculation.

The PRESTO scheme is used for pressure interpolation. The QUICK scheme is used for the discretisation of the momentum and the energy equations. The Piecewise-Linear Interface Calculation (PLIC) scheme is used for the discretisation of the volume fraction equation. When simulations are made using the CSF method, the default node based smoothing of the volume fraction field prior to calculation of the curvature was enabled. No smoothing of the calculated curvatures was performed. A first order implicit scheme is used for the temporal discretisation of the transient terms. Finally, for the pressure-velocity coupling, the SIMPLE algorithm is used. Globally scaled residuals are used and the residual targets for all the equations are set to 1 × 10 -6 . The mesh consists of square and orthonormal cells.

The mesh with R/∆ x = 20 ratio is shown in Fig. 42. The properties of the electrolyte and the gas considered (H 2 ) are shown in the

ρ l [kg • m -3 ] ρ g [kg • m -3 ] µ l [kg • m -1 • s -1 ] µ g [kg • m -1 • s -1 ] 1000 0.0899 1.2 × 10 -3 8.79 × 10 -6
tension varies with the concentration of dissolved gases [START_REF] Lubetkin | The motion of electrolytic gas bubbles near electrodes[END_REF]. To be consistent with the experiments performed on microelectrodes, the initial surface tension was set at [START_REF] Glas | Measurements of the growth of electrolytic bubbles[END_REF][START_REF] Liu | Numerical simulation of hydrogen bubble growth at an electrode surface[END_REF].

0.075 [N • m -1 ] [Glas
The work of [START_REF] Massoudi | Effect of pressure on the surface tension of water. Adsorption of low molecular weight gases on water at 25.deg[END_REF] has been able to establish relationships between the variation in partial pressure of gas and surface tension . At low pressures, the concentration of dissolved hydrogen and the partial hydrogen pressure can be related through Henry's law:

c h 2 = p K H (4.21)
where K H the constant of proportionality is dependent on the temperature and pressure, but in our case can be estimated as [START_REF] Wiebe | The Solubility of Hydrogen in Water at 0, 50, 75 and 100°from 25 to 1000 Atmospheres[END_REF][START_REF] Sander | Compilation of Henry's law constants (version 4.0) for water as solvent[END_REF]:

K H = 7.8 × 10 -6 [mol • m -3 • Pa -1 ] (4.22)
As a result, the variation of surface tension as a function of concentration can be obtained:

∂γ ∂c H 2 = 1 K H ∂γ ∂p = -3.2 × 10 -5 [N • m 3 • m -1 • mol -1 ] (4.23)
A concentration gradient was initially applied along the bubble interface as shown in Fig. 42. From this concentration gradient results a surface tension gradient along the bubble interface which generates a Marangoni current as shown in Fig. 41 .

This artificial situation is only intended to illustrate visually the impact of spurious currents on a simulation, and the improvement that can be made by modelling the surface tension with height functions as described in the previous chapters. In the image on the left with the CSF method, spurious currents appear along the interface, which disrupt the Marangoni currents. This is not the case on the left image with the height function method.

In order to disentangle the effects of spurious currents from Marangoni currents, the overall approach presented in this chapter has been to evaluate spurious currents when there is no simulated Marangoni effect as described earlier in the case of the static bubble. This last simulation of this chapter tends to illustrate the error that spurious currents could have generated on the calculation of Marangoni currents if a generic methodology such as the CSF methodology had been used. Overall, the surface tension variation part does not create spurious currents if the curvature calculation is correct in the first place. 1 The Marangoni effect as an alternative to diffusive transport 1.1 Alternative to the concept of diffusion in the Nernst's layer

The Marangoni effect is a transient phenomenon, but in the case of electrogenerated bubbles it is sustained by the Joule effect and the production of gas at the electrode. The observation of vortices by Yang et al. is most likely due to the Marangoni effect as hypothesised by Lubetkin [START_REF] Lubetkin | The motion of electrolytic gas bubbles near electrodes[END_REF][START_REF] Yang | Marangoni convection at electrogenerated hydrogen bubbles[END_REF]. This calls into question the hypothesis of transport by pure diffusion of the species produced at the electrode towards the bubble interface.

The microscopic chaotic motions that take place in the viscous sublayer adjacent to the surface of an electrode have long been regarded as diffusive motion. Moreover the bubble growth laws assuming a pure diffusive motion give a relatively good approximation to what is observed. The work of Amatore et al. sheds light on this interpretation [START_REF] Amatore | The real meaning of Nernst's steady diffusion layer concept under non-forced hydrodynamic conditions. A simple model based on Levich's seminal view of convection[END_REF]. In particular, they introduce a diffusion coefficient D conv that depends on the species flux and the convective movement over a given distance:

D conv = v < ∆ > (5.1)
where < ∆ > is the length over which the dissolved species travel during a time interval, and v is the average velocity of chaotic motion that develops over this length. The difference with a classical diffusion coefficient obeying the Nernst relationships is that this diffusion coefficient is spacially dependent. On the basis of this analysis and without the possibility of observing the currents in the vicinity of the bubble, it is understandable that chaotic movements that were in fact generated by a Marangoni effect could have been interpreted in the past on the basis of a diffusive phenomenon. Yang et al. observed vortex flows around the bubble as shown in Fig. 6 and the only possible interpretation is the presence of stress on the interface caused by a variation in the surface tension [START_REF] Yang | Marangoni convection at electrogenerated hydrogen bubbles[END_REF].

The presence of a Marangoni effect around electrogenerated bubbles is still a recent discovery.

While the vortex currents resulting from the Marangoni effect can be visualised, it is not yet clear how they can affect the growth and detachment of a bubble on a conventional electrode. There is also no consensus on their origins.

Existence of interfacial gradients

The Marangoni effect is due to a variation of the surface tension. This variation as a function of temperature and concentration of different surfactants or pollutants, is essential to understand the formation of electrogenerated bubbles [START_REF] Lubetkin | The motion of electrolytic gas bubbles near electrodes[END_REF][START_REF] Massoudi | Effect of pressure on the surface tension of water. Adsorption of low molecular weight gases on water at 25.deg[END_REF][START_REF] Weissenborn | Surface Tension of Aqueous Solutions of Electrolytes: Relationship with Ion Hydration, Oxygen Solubility, and Bubble Coalescence[END_REF][START_REF] Zhang | Evaluating the Behavior of Electrolytic Gas Bubbles and Their Effect on the Cell Voltage in Alkaline Water Electrolysis[END_REF].

The interfacial gradients of these physical properties vary the surface tension along the interface. Therefore, due to the unbalanced forces at the interface, fluid elements experience a net shear stress there and move to regions of the interface where the interfacial tension is higher. The assumptions made about the strength of the Marangoni effects depend directly on the knowledge of the gradients of the quantities involved. The important physical prerequisite is the presence of a sufficiently long-lasting gradient to allow observation of the resulting motion. If the presence of these gradients at the interface is proven, their magnitudes remain a model hypothesis.

The convection generated by these ends very quickly, as it promotes rapid equilibration of the temperature or concentration distribution along the interface. To sustain convection for a long time, a mechanism must be provided to maintain the surface tension gradient.

In the case of the thermocapillary effect, the Joule effect maintains the phenomenon. A distinction must then be made between the solutal capillary effect caused by pollutants or surfactants and that caused by dissolved gas species. In the first case, the surfactants present in the electrolyte only modify the value of the surface tension of the interface when they are adsorbed. There are no sources of surfactants that could maintain the phenomenon over the long term. In the second case, dissolved gas species are generated at the electrode and then absorbed by the bubble as it grows. There is therefore a driving force here that could sustain the Marangoni effect over the long term.

To understand how temperature or concentration gradients occur it is worth mentioning the couplings that take place, as they highlight the complexity of the competing phenomena and their influence on the Marangoni effect. The couplings between the model equations are shown in Fig. 43. The implications of Fig. 43 are explained in section 3. Before commenting on this diagram and as a pre-requisite, the heat equation and the transport equation for dissolved species must be stated here.

2 Heat, transport, and surfactant models

Heat equation

As the surface tension is a function of the temperature, in order to study this variation it is necessary to include in the model the heat equation Eq. (5.2) as well as the Laplace equation Eq. ( 5.3) to calculate the Joule effect Eq. (5.5), and the source term S T of the heat equation: where j is the curent density in A • m -2 . As shown in Fig. 43, at the interface, as the temperature varies, the value of the surface tension will be modified. This generates the Marangoni effect. The velocity at the interface is modified which changes the temperature profile around the bubble.

ρC p ∂T ∂t + v • ∇T = ∇ • (λ∇T ) + S T (5.2) ∇ 2 Φ = 0 (5.3) j = λ elec ∇Φ (5.
ṁ v + pI • n I -2µD • n I = γκn I + ∇ s γ(T, c) (5.6)
The Marangoni effect will therefore depend on the Joule effect. The conductivity of the electrolyte and the size and shape of the electrode will therefore have a direct

Transport of dissolved species

In order to evaluate the interfacial mass transfer it is necessary to know the concentration near the interface, to do this the transport equation of dissolved species must be introduced:

∂c ∂t + v • ∇c = ∇ • (D∇c) + S electrode + S interf ace (5.7)
Two source terms are to be taken into account in this equation: the production of dissolved species at the electrode S electrode and their absorption at the interface S abs . The source term at the electrode can be calculated by Faraday equation

ṁelectrode = |j|M F v (5.8)
where ṁelectrode [kg • m -2 • s -1 ] is the production rate of dissolved species at the electrode.

The second source term will depend on ṁ which is found in the jump relation established from the conservation of mass at the interface:

ṁ = ρ g (v g -v I ) • n I = ρ l (v l -v I ) • n I (5.9)
This transfer rate will depend on the assumptions made about interfacial mass transfer and the concentration of dissolved species closest to the interface. As the bubble grows, the interface is supersaturated on average. It is the difference between the interface concentration and the saturation concentration of the dissolved gas that drives this growth. This changes the transport equation, which in turn changes the value of the surface tension. As with the thermo-capillary effect, this variation in surface tension creates a Marangoni effect which modulates the flow and thus the transport of dissolved species.

Presence of contaminants

The presence of surfactants/pollutants in the electrolyte is a very likely hypothesis. Their impact on the surface tension cannot be neglected. Levich et al. show that the adsorption of an insoluble surfactant stiffens the interface [START_REF] Levich | Surface-Tension-Driven Phenomena[END_REF]. Such a rigid interface is unable to transfer tangential stress, which of course blocks the effects of the surface tension gradient at the interface. The surfactants move along the interface and clump together to form this stiffened area. A first modelling approach is to consider two zones on the interface, one stiffened by the surfactants and a mobile one where other sources of surface tension variation are not inhibited by the surfactant influence. The conclusions are limited to the surfactant distribution for a fixed bubble size whose diameter corresponds to a stage in the bubble cycle. The proportion of the surface covered by surfactants is unknown and can only be estimated a posteriori by comparing the modelling results with the experimental results. As shown in Fig. 35, in the case of a hydrogen bubble growing on the surface of a microelectrode, Meulenbroek et al. assume that the surface stiffened by the surfactants covers the top of the bubble and forms a cap [START_REF] Meulenbroek | Competing Marangoni effects form a stagnant cap on the interface of a hydrogen bubble attached to a microelectrode[END_REF]. They define a stagnation angle θ S to give the position of the surfactant stagnation point on the interface, below which the interface is free of surfactant. Using the contact angle θ c , the area for which surface tension changes are not inhibited can be calculated.

Couplings

Model uncertainties

The Marangoni effect is a phenomenon generated at the interface that has an influence on what happens in the bulk. Including it in the analysis of the development of electrogenerated bubbles profoundly changes the way in which bubble behaviour should be interpreted. In view of what has been written above, the diagram in Fig. 43 contains all the elements for a holistic direct numerical simulation. Four equations must be solved in the electrolyte, and two jump relations at the interface. To this must be added a model to take into account the presence of surfactants on the interface. In the case of electrogenerated bubbles the surface tension of the bubble can be affected by either the temperature, the concentration of any surfactants and the concentration of dissolved gases. For the sake of completeness, the electrocapillary effect should be added. The electrocapillary effect is a phenomenon that occurs at the gas-liquid interface due to the presence of a surface charge. The electrocapillary effect can cause currents [START_REF] Johnson | Electrocapillary Flows[END_REF]. However, compared to the thermocapillary and solutocapillary effects, the electrocapillary effect is less well understood. It should also be added here that no indication of the presence of surfactants within the electrolyte can be measured. It remains an unknown which can only be deduced a posteriori from the experimental results [START_REF] Sellier | Unraveling surfactant transport on a thin liquid film[END_REF]. This leaves temperature and dissolved gas concentration for which the uncertainties in the equations for monitoring them are smaller. The quantities involved in these equations are known. However, as previously mentioned, solutions for monitoring the Marangoni effect on a free interface are not common. Numerical studies of the Marangoni effect on an electrogenerated bubble assume that the interface is fixed, i.e. that the growth of the bubble is negligible with respect to the time interval of study considered [START_REF] Massing | Thermocapillary convection during hydrogen evolution at microelectrodes[END_REF][START_REF] Hosokawa | Evaluation of adsorption of surfactant at a moving interface of a single spherical drop[END_REF][START_REF] Meulenbroek | Competing Marangoni effects form a stagnant cap on the interface of a hydrogen bubble attached to a microelectrode[END_REF]. These studies study the bubble at one stage of its growth and aim to reproduce the velocity and temperature field observed experimentally. In this context, it is legitimate to question the initial conditions used.

The argument developed by the above authors is that the temperature and velocity field evolve rapidly enough to consider that the steady solution of a numerical simulation assuming a fixed interface can account for the development of an electrogenerated bubble from a microelectrode in a real situation. The assumption made here is that the parameter that can make the numerical system evolve is the growth of the bubble. By deriving the relationship of bubble diameter versus time for reaction rate limited bubble growth Eq. ( 2.50), it can be shown that for long times at the end of the bubble's development the growth rate of the bubble becomes less and less important:

d R dt = β 1 3t 2/3
(5.11)

Using a value of β = 360 found experimentally by Massing et al. and from a bubble radius of R = 560 µm, a characteristic bubble growth time in the experiment can be expressed:

t (growth,r=560 µm) = R dR/dt ≈ 6.7s
(5.12)

By comparing this time with the recirculation time within the vortex of 0.12s, it can be concluded that the flow field will grow rapidly relative to the bubble growth rate. In other words, it is possible to neglect the bubble growth rate by using a fixed bubble size. However, this assumption is only valid at the end of the bubble's growth, in the first few moments, for example, for a bubble radius of 50µm, we obtain a characteristic bubble growth time of about 1s. The growth rate of the bubble is thus no longer negligible.

Mass transfer and surface tension variation

Concerning the transport of dissolved species it is appropriate to examine the corresponding Peclet number:

P e = u R D H 2 ≈ 758 (5.13)
Microconvection movements clearly dominate diffusive phenomena in the case of an electrogenerated bubble on a microelectrode. One of the objectives of a model based on the diagram in Fig. 43 is to be able to determine the growth rate of the bubble. This growth rate depends on the interfacial mass transfer. This depends on the transport of dissolved species in the electrolyte which is a result of microconvection currents. To perform a direct simulation it is necessary to be able to associate the interfacial mass transfer to the microconvective currents in the vicinity of the bubble, i.e. to the Marangoni effect. As previously discussed, from a numerical point of view, it is necessary to be able to combine mass transfer and surface tension variation at the interface and more particularly at the contact line between the solid, liquid and gaseous phase where the problem of a moving contact line must be taken into account. At the level of the contact line, deformations are important and it is probable that the hypothesis of a spherical interface does not hold anymore. This deformation could influence the velocity profile and the pressure distribution. In addition, it has been observed that a mat of microbubbles is formed under the detaching large bubble, which could also significantly influence the detachment and the force balance [START_REF] Bashkatov | Oscillating Hydrogen Bubbles at Pt Microelectrodes[END_REF]. Solving these problems, in the form of simulations with a growing and deformable bubble interface, is a crucial step in understanding bubble detachment.

Competing Marangoni effect

The other emerging issue is the competition between the different Marangoni effects. Lubetkin assumed that the magnitude of the solutal Marangoni force is greater than thermal marangoni [START_REF] Lubetkin | The motion of electrolytic gas bubbles near electrodes[END_REF]. A first approach to assess the competition between the different effects is to estimate and compare the different surface tension gradients ∇ s γ along the interface according to their source. The surface tension gradient of the thermocapillary effect can be approximated by:

(∇ s γ) T ≈ ∂γ ∂T ∂T ∂z (5.14)
As shown in Tab. 5.2, the value of ∂γ ∂T varies according to the references considered [START_REF] Young | The motion of bubbles in a vertical temperature gradient[END_REF][START_REF] Hardy | The motion of bubbles in a vertical temperature gradient[END_REF][START_REF] Prigogine | Statistical mechanics of surface tension and adsorption[END_REF][START_REF] Morick | Migration of Air Bubbles in Silicone Oil under the Action of Buoyancy and Thermocapillarity[END_REF][START_REF] Vazquez | Surface Tension of Alcohol Water + Water from 20 to 50 .degree.C[END_REF]. The low value and the typical value were reported as low and typical by Lubetkin in his review [START_REF] Lubetkin | The motion of electrolytic gas bubbles near electrodes[END_REF]. The high value comes from the work of Vazquez et al. [START_REF] Vazquez | Surface Tension of Alcohol Water + Water from 20 to 50 .degree.C[END_REF]. This value was cited by Yang et al By carrying out a numerical simulation for a conventional electrode the authors obtain ∆T ≈ 1.5K for Θ = 0.87 and ∆T ≈ 0.5K for Θ = 0.31. The greater the bubble coverage, the greater the Joule effect. The temperature difference between a microelectrode and an electrode decreases by a factor of 10. By taking into account the length at the interface for which this temperature difference can be observed, we can obtain an approximation of the temperature gradient. The surface tension gradient can then be calculated:

(∇ s γ) T [N • m -2 ] Low High Microelectrode 6 22 Regular electrode 0.1 0.4
The surface tension gradient of the solutocapillary effect of dissolved gas can be approximated by:

(∇ s γ) c ≈
∂γ ∂c H 2 = 1 k H ∂γ ∂p = -3.2 × 10 -5 N m 3 m -1 mol -1
(5.17)

As shown in Tab.5.4, the value of ∂c ∂z varies according to the references considered [START_REF] Westerheide | Isothermal growth of hydrogen bubbles during electrolysis[END_REF][START_REF] Shibata | The Concentration of Molecular Hydrogen on the Platinum Cathode[END_REF][START_REF] Glas | Measurements of the growth of electrolytic bubbles[END_REF][START_REF] Sides | A Close View of Gas Evolution from the Back Side of a Transparent Electrode[END_REF] The surface tension gradient can then be calculated:

(∇ s γ) c [N • m -2 ]
Low typical High 0.1 14 500

In view of the values obtained for the surface gradients, it is difficult to say which capillary effect may dominate. By examining the ratio (∇ s γ) c /(∇ s γ) T it can vary from 0.5 10 -3 to 5 10 3 . The dominant effect will be determined essentially from the experimental conditions.

Bubble growth rate with Marangoni Flow

There is currently no model to simulate the growth of a bubble in conjunction with the Marangoni effect. It is therefore not possible to establish a relationship of radius growth as a function of time R(t). On the other hand, there are data on the velocity field around the bubble created by the Marangoni effect for a given bubble size.

Growth with the diffusion hypothesis

To close the bubble growth model it is necessary to express a value for interfacial mass transfer. Indirectly it is necessary to be able to estimate the interfacial mass transfer coefficient k [m • s -1 ]. The equation Eq. (2.43) modelling a bubble with no contact with the electrode allow the variation of the concentration as a function of radius to be obtained. Under this assumption, and by combining this equation with Eq. (2.40), it is possible to obtain an expression for interfacial mass transfer ṁdiff,hyp by dividing by the area of a sphere:

ṁdiff,hyp = M D ∂c ∂r R = M (c 0 -c s ) D R 1 + R √ πDt (5.18)
Then the mass transfer coefficient k can be estimated by:

k = D R + D πt (5.19)
This mass transfer coefficient is used to express the mass transfer rate ṁ appearing in the interfacial mass transfer equation, thus closing the model:

ṁ = M k(c i -c s ) (5.20)
This relationship was used by Lui et al. to simulate the growth of a hydrogen bubble in a numerical simulation on the surface of an electrode [START_REF] Liu | Numerical simulation of hydrogen bubble growth at an electrode surface[END_REF]. The authors do not include in this study the Marangoni effect. However the results concerning the growth of a bubble can be considered to be in poor agreement with the measurements made by Glas and Westwater [START_REF] Glas | Measurements of the growth of electrolytic bubbles[END_REF]. The mass transfer coefficient here can be considered as a global variable determined at the bubble scale. The assumptions made here are strong and vary greatly from reality. The bubble is considered as a sphere. The analytical solution of the equation starts from the assumption that the bubble is in the middle of an infinite electrolyte of homogenous concentration c 0 in which the transport of dissolved species takes place only by diffusion. In contrast to this theoretical value c 0 , the interfacial concentration c i is not homogeneous along the interface.

Growth with the penetration theory

The penetration theory was suggested by Higbie who was investigating whether or not a resistance to transfer existed at the interface when a pure gas was absorbed in a liquid [Higbie, 1935b]. In the penetration theory, the mass transfer coefficient of the fluid element with the residence time is expressed by the following equation:

k = 2 D πt c (5.21)
where t c is the contact time of a fluid element containing the dissolved species with the bubble interface.

On the electrolyte side we can evaluate the growth rate of the bubble ṁB [kg • s -1 ] from k:

ṁB = M k(c i -c s )S interf ace (5.22)
Where S interf ace is the surface on which the mass transfer is performed. In other words, it is the surface not covered by the surfactants. As shown in Fig. 45 the assumption of a stagnant surfactant cap prevents interfacial mass transfer.

By combining this relationship with what is happening on the gas side Eq. (5.23):

ṁB = ρ g 4πR(t) 2 dR(t) dt [kg • s -1 ] (5.23)
we can determine the growth rate of the bubble:

dR(t) dt = M ρ g 2 D πt c (c i -c s ) S interf ace 4πR(t) 2 (5.24)
This relationship is close to that of the Eq. (2.44), and helps to validate the reasoning that the use of the penetration theory is appropriate.

The objective here is to establish that Eq. ( 5.24) is suitable to study of bubble growth. With the experimental protocol used by Yang et al. it is possible to measure the fluid velocity as close to the interface as possible [START_REF] Yang | Marangoni convection at electrogenerated hydrogen bubbles[END_REF]. Similarly, this value can be deduced from a numerical simulation such as that of [START_REF] Meulenbroek | Competing Marangoni effects form a stagnant cap on the interface of a hydrogen bubble attached to a microelectrode[END_REF].

Using the data from this study an average velocity along the interface can be calculated

v i,avg = 12[m • s -1 ].
The transfer area S interf ace can be calculated from the bubble radius, the contact angle θ c = 4.2°and the stagnation angle θ S = 57°estimated by the same authors. The time t c is the contact time of the fluid with the interface. This contact time can be deduced from the average velocity along the interface. From these data we can estimate an average contact time along the interface of t c ≈ 0.0414s. Tab. 5.5 shows the other data used. The main unknown here is the average concentration at the interface c i,avg . The Tab. 5.5: Parameters common to both equations, Eq. (2.44) and Eq. (5.24)

Quantity Unit Estimate R m 5.6 • 10 -4 D H 2 m 2 • s -1 7.38 • 10 -9 ρ H 2 kg • m -3 0.09 M H 2 kg • mol -2 0.002
growth rate of the bubble measured by Massing et al. is dR(t) dt ≈ 84 [µm s -1 ]. Using this value it is possible to find c i,avg = 37.5 [mol • m -3 ] from Eq. (5.24). Which is a plausible value [START_REF] Westerheide | Isothermal growth of hydrogen bubbles during electrolysis[END_REF][START_REF] Shibata | The Concentration of Molecular Hydrogen on the Platinum Cathode[END_REF][START_REF] Vogt | On the supersaturation of gas in the concentration boundary layer of gas evolving electrodes[END_REF]. However, many parameters are approximated. Eq. (5.24) needs to be tested with more experimental results to be validated. This relationship is thus a new tool to better understand what is happening around the bubble.

6 Non-dimensional study

Simplification of the model

There are few studies that relate the Marangoni effect to interfacial mass transfer. The aim of this dimensionless study is to evaluate the interfacial mass transfer by means of a Sherwood number from an estimated or measured Marangoni effect as in the Yang et al. experiment [START_REF] Yang | Marangoni convection at electrogenerated hydrogen bubbles[END_REF]. In a previous section we were able to estimate the errors of the numerical model, which will allow us to obtain a reliable evaluation of the Sherwood number. The main assumption used here is the use of the penetration theory. The other assumption made here is that the quantity influencing the value of the surface tension at the interface can be compared to a temperature, i.e. the coupling between the absorption of dissolved species at the interface and the resulting change in surface tension is not taken into account. The surface tension gradient is generated from the temperature profile deduced from the heat equation. The variable of interest here is the surface tension along the interface. This type of study is thus in line with the hypotheses made by Massing et al. Meulenbroek et al. and Hossain et al. which • L interf ace the length of the interface over which the interfacial mass transfer and the Marangoni effect take place;

• κ the curvature of this interface;

• ∇ s T the temperature gradient along the interface. From this temperature surface gradient can be deduced the surface tension gradient, the surface tension in this model being a function of temperature only.

Dimensionless equations

A surface tension is a force per unit length, so the resulting stress must scale as ∆γ/L , while the viscous stress scales as µu/L, for u the speed of the Marangoni flow. Equating the two we have a flow speed u = ∆γ/µ. As Ma is a type of Péclet number, it is a velocity times a length, divided by a diffusion constant D. The Marangoni number can be defined as :

M a = v char L char D ≈ |∇ s γL 2 char | µD T ≈ |∆γL char | µD T (5.25)
M a = advective transport rate due to surface tension gradient diffusive transport rate of the quantity of interest (5.26)

In order to make the Marangoni number appear, the Navier-Stokes and transport equations and the corresponding jump relations are made dimensionless. As illustrated in Fig. 47 and Fig. 48, the aim here is to expose the perimeter of the bubble to different gradient values to generate a Marangoni effect. Due to blocking by surfactants, only a portion of the interface will participate in the surface tension gradient, and interfacial mass transfer. The objective of this study is to create a tool to relate a surface tension gradient through a Marangoni number to the interfacial mass transfer through a Sherwood number, the common characteristic length of these two dimensional numbers is the surfactant free interface length L interf ace . Locally in an experiment it is possible to estimate this length. As mentioned earlier the Marangoni effect is intermittent, so the transfer of interfacial mass is only retained when the velocity of the Marangoni effect is highest. The graph in Fig. 47 details the correspondence between the observation that can be made experimentally and the numerical study performed. The interface is assumed to be exposed to a temperature gradient over a length of interface.

Mass transfer is only simulated along the length L interf ace . It is assumed that Surfactants block interfacial mass transfer, covering the bubble interface, which prevents interfacial mass transfer of dissolved species. A condition has been added so that interfacial mass transfer is blocked outside the length L interf ace . In Fig. 47 the diagram on the left shows a bubble attached to an electrode. However, the numerical model presented in Chapter 3 does not include a contact line module. Therefore, the choice was made to perform a simulation on a bubble in a stagnant liquid. The idea here is to simulate a surface tension variation along the interface to obtain local information on the interfacial mass transfer. As suggested by Hossain et al. the surface tension variation along the interface is not necessarily uniform, nor linear [START_REF] Hosokawa | Evaluation of adsorption of surfactant at a moving interface of a single spherical drop[END_REF]. The interface can be divided into several sections where the surface tension stresses are different. The Marangoni currents can be of opposite direction depending on the portion of the interface considered as suggested in Fig. 44. The objective of this dimensionless study is to develop a tool to determine the local interfacial mass transfer as a function of the interface section considered.

Fig. 48 shows the result of one of the simulations carried out to produce the graphs in Fig. 49 and Fig. 50. In each cell of the mesh the fluid will be exposed to the interface for a time The graphs in Fig. 49 and Fig. 50 show the results of many simulations. As an example, the Tab.5.6 shows the data of the point indicated by the label data in Fig. 49 and Fig. 50. The Sherwood numbers and the mass transfer coefficient are calculated as with the methodology indicated just above. The Marangoni number is calculated from the input parameters as indicated by Eq. (5.25). Marangoni currents were generated along a spherical interface over a distance L interf ace . These graphs allow to obtain local information on the interfacial mass transfer, and can be used accordingly for experimental observations. In accordance with the observations made by Golovin, the Sherwood number increases as a function of the Marangoni number [START_REF] Golovin | Mass transfer under interfacial turbulence: kinetic regulaties[END_REF]. The graph in Fig. 49 shows the Sherwood number as a function of the Marangoni number. In the framework of a simulation made on a fixed liquid-bubble interface this allows to estimate the Sherwood number describing the interfacial mass transfer. The graph in Fig. 50 shows the mass transfer coefficient as a function of the Marangoni effect. The three curves allow to distinguish different ratios of bubble length to bubble radius. The overall objective here was to provide a tool to better understand the mass transfer due to the Marangoni effect.

Conclusion of the chapter and discussion of model uncertainties

The models described above do not confirm the real situation of interfacial mass transfer at the bubble interface. There are too many unknowns to form definite conclusion. However, the implementation of a contact line model in the numerical model could remove many uncertainties. At this stage it is appropriate to list the main quantities for which these uncertainties exist:

• Concentration gradients around the bubble, i.e. the local value of the concentration of dissolved species along the interface;

• the value of the surface tension gradient along the interface;

• local value of the interfacial mass transfer coefficient;

• presence and impact of surfactants;

• contact angle.

A direct numerical simulation of a holistic model such as the one shown in Fig. 43 would provide a microscopic approach to the problem. One of the major objectives of this thesis has been to develop a model to address these uncertainties. In most models, the interfacial mass transfer coefficient remains an input parameter of the model. The advantage of the penetration theory is to integrate the mass transfer coefficient as a local value calculated directly by the model, and not as an average quantity calculated a priori. The validity of the penetration theory in the case of electrogenerated bubbles would remove the uncertainty regarding the mass transfer coefficient. The coupling between dissolved species transport, interfacial mass transfer and the momentum jump relation would remove the uncertainty in the concentration gradient around the bubble. Removing this uncertainty would provide a value for the soluto Marangoni effect, which could distinguish its importance from other effects, and thus provide a value for the surface tension gradient along the interface. With the development of a holistic model, however, the value of the contact angle and the impact of surfactants would remain unknown. There is currently no agreement on the theory to predict the evolution of the contact angle. The impact of surfactants can only be assessed after an analysis of the results of a model with regard to the experimental results.

Conclusion 6

1 Perspective: integration of a moving contact line model

As previously stated, in order to simulate a holistic model like the one shown in Fig. 43, a moving contact line model must be included. As such it is consistent to introduce the problem of simulating these moving contact lines, and how this can be integrated into the numerical model presented in this thesis through height functions. The contact line is a singularity of an even higher order than the interfaces: it is a linear singularity located at the edge of a surface singularity. Not only are the physical quantities not continuous, which poses numerical problems, but even a simple physical model is not always available to describe them [START_REF] Legendre | Comparison between numerical models for the simulation of moving contact lines[END_REF].

Experience shows that the behaviour of contact lines depends on a large number of factors, few of which are controlled in practice (wall roughness, surface condition and chemical contamination, composition of the fluid and possible contamination by surfactants, etc). These lines of contact play a crucial role in the phenomenon we want to study. If we want to study the influence of a hydrophobic surface compared to a hydrophilic surface, it is essential that the model used can account for their differences. The Navier-Stokes equations include some aspects of intermolecular forces. Viscosity, pressure or surface tension are all quantities that result directly from the interactions between the molecules of the fluid. The application of these interactions to the contact lines leads very naturally to the surface energy of the Young-Dupré relationship and the link with the equilibrium contact angle is obvious.

γ lg cos θ E = γ sg -γ sl (6.1)
where θ E is the stactic contact angle and can be defined as the angle that the interface makes with the wall when the contact line is stationary. Unfortunately, , this relationship is only relevant for a static contact line without mass transfer. In the presence of mass transfer, not only is the contact angle modified, but the balances and exchanges are modified by the presence of the wall.

It is easy to find theoretical studies on the dynamics of contact lines and contact angle. These studies are based on classical continuous medium mechanics and their two essential ingredients are the Navier-Stokes equation and the Laplace relationship. There are also other approaches to model contact lines, where the microscopic aspects are more emphasized [START_REF] Hadjiconstantinou | Hybrid Atomistic-Continuum Formulations and the Moving Contact-Line Problem[END_REF][START_REF] Gouin | The wetting problem of fluids on solid surfaces: Dynamics of lines and contact angle hysteresis[END_REF][START_REF] Gouin | Energy of Interaction between Solid Surfaces and Liquids[END_REF][START_REF] Pomeau | Recent progress in the moving contact line problem: a review[END_REF]. Experiments on the subject often use perfect surfaces (usually glass, sometimes coated to modify its wettability) and equally ideal fluids. In [START_REF] Gennes | Wetting: statics and dynamics[END_REF], de Gennes draws a particularly complete panorama of the different aspects of wetting by incorporating a very large number of physical ingredients. Studies of contact angle hysteresis are rarer.

The common feature of these models is the calculation of the liquid-gas interface profile from a simplified solution of the Stokes flow near the contact line. This flow gives rise to a non-integrable singularity [START_REF] Hocking | A moving fluid interface. Part 2. The removal of the force singularity by a slip flow[END_REF] which is solved by introducing a cut-off scale and a physical mechanism at the molecular scale to close the system of equations. When a contact line moves at a velocity U cl along a wall and a non-slip condition is imposed, a stress is generated:

τ singularity ≈ µ U cl ∆ (6.2)
where ∆ is the grid spacing and µ is the fluid viscosity. When ∆ tends towards zero this constraint τ singularity diverges. Refining the mesh makes the calculations diverge. To deal with the singularity several authors introduce the navier slip condition in their model. In fact, many models can be interpreted as variants of the Dussan model [START_REF] Dussan | The moving contact line: the slip boundary condition[END_REF], where microscopic phenomena are summarized by a microscopic contact angle and a slip length.

In such model the tangential component to the wall of the velocity is estimated using the following relation:

U w = λ N ∂u ∂n w (6.3)
where U w is the fluid velocity at the wall, n w is the normal to the wall, and λ N is the slip length, which is usually estimated to be of nanometer scale. Convergence of the grid is then achieved by solving the complete hydrodynamic problem within the moving contact line region.

However, λ N values are unrealistic in most simulations, they are too large. The reason for this is the limited refinement of the grid. So in practice λ N becomes an adjustable parameter for the simulation. The calculations converge but the slip condition becomes unphysical.

This boundary condition for the velocity field is used in conjunction with a dynamic contact angle model. The first step is to determine the dynamic contact angle. One of the solution is to assume the dynamic contact angle to be constant and equal to the static contact angle.

Then in the VOF framework, the idea is to impose on the contact line a normal to the interface depending on the contact angle.

n I = sinθn w,∥ + cosθn w,⊥ (6.4)

where n w,∥ and n w,⊥ are the components of the normal vector, parallel and normal to the wall. This value is used into the surface tension model. Then it is imposed as a boundary condition .

Most models vary depending on how the contact angle is determined and the description of the slip condition. However the solution in this type of model depends on the size of the mesh [START_REF] Afkhami | A mesh-dependent model for applying dynamic contact angles to VOF simulations[END_REF]. Afkhami et al. based on the theory of [START_REF] Cox | The dynamics of the spreading of liquids on a solid surface. Part 1. Viscous flow[END_REF] proposed the following expression for modeling the contact angle [START_REF] Afkhami | A mesh-dependent model for applying dynamic contact angles to VOF simulations[END_REF][START_REF] Sheng | Immiscible-fluid displacement: Contact-line dynamics and the velocity-dependent capillary pressure[END_REF] : cos(θ num ) = cos(θ app ) + 5.63Ca ln( K ∆/2 ) (6.5)

where θ num is an angle defined according to the mesh. They point out that there is a linear dependance of cos(θ num ) -cos(θ app ) on Ca ln( K ∆/2 ) when applying both no-slip and Navierslip boundary conditions. Ca is the capillary number, Ca = µU cl /γ (U cl is the contact line velocity, µ the fluid viscosity ). K is a constant which can be determined by fitting numerical data to data obtained experimentally. The big advantage of this method is that it eliminates the stress singularity at the contact line; the solutions converge with mesh refinement. A shown on Fig. 51, in order to model the numerical contact angle, Afkhami et al. used a methodology based on height functions [START_REF] Afkhami | Height functions for applying contact angles to 2D VOF simulations[END_REF]. This model implemented in the VOF framework has been used in a more recent publication [START_REF] Afkhami | Transition in a numerical model of contact line dynamics and forced dewetting[END_REF].

Conclusion of the work

This work was carried out with the aim of gathering the knowledge necessary to model bubbles generated from electrolysis cells.

The work started with a thorough literature review of the behaviour of bubbles in electrolysis cells and the effect of these bubbles on the efficiency of the process. The outgassing produced in electrolysis cells has been studied extensively by several generations of researchers. This study revealed that knowledge of outgassing characteristics, such as size at detachment, growth rate or bubble movement, is fragmentary. The aim of this literature review is twofold, we have sought to describe the knowledge needed to model an electrolysis cell but also to highlight aspects that could improve the electrochemical process. The covering of the electrode by the bubbles is an essential parameter in the calculation of the efficiency of the electrolysis process. Therefore, without knowledge of the residence time of the bubbles on the electrode, their growth rate, and their diameter at the moment of detachment, it is not possible to model the electrolysis process correctly. It turned out that a better understanding of the transport of dissolved species, the interfacial mass transfer, and the Marangoni effect contributing to delay the detachment of bubbles from the electrode, would allow to elucidate the behaviour of bubbles in certain experimental cases. The desire to produce a holistic model of an electrogenerated bubble following this bibliographical study stems from the observation of the gaps concerning these three points. Without the implementation of a model capable of simulating a Marangoni effect on a free surface, an essential simulation tool is missing. In this manuscript, the mathematical model to simulate two-phase flows with a Marangoni effect and interfacial mass transfer has been presented. Then, the numerical methods needed to apply this model have been outlined. An algorithm based on VOF methodology and height functions was detailed. Finally, it was shown that this tool successfully reduces parasitic currents on a static and moving isolated bubble and is suitable for the study of the Marangoni effect around an electrogenerated bubble. This tool is proving to be a first step towards a holistic model. The couplings necessary to perform a direct simulation of the phenomenon were presented. In order to achieve this objective, the implementation of a moving contact line module is missing. However, a parametric study to account for the interfacial mass transfer as a function of the Marangoni effect was carried out using the numerical model. In the absence of a direct holistic simulation, an equation giving the growth rate of a bubble, based on the penetration theory, allowing the Marangoni effect to be taken into account, as well as the covering of the bubble by surfactants, was established. This relationship is close to the Eptein-Plesset equation, and gives comparable results. It was pointed out that the amount of experimental data is still insufficient to validate a numerical model of an electrogenerated bubble including the Marangoni effect. The work carried out in this manuscript has laid the theoretical and numerical foundations for an improved understanding of dissolved gas transport, interfacial transfer and the Marangoni effect around an electrogenerated bubble.

As mentioned earlier, the work presented in this thesis is part of an approach based on inverse problems. In this perspective, among all the sources of errors mentioned in the diagram in Fig. 13 , this thesis has focused on reducing three of them:

• Errors in hypothesis and models: They have been extensively discussed in this thesis. Different models have been presented to account for the transfer of species from the electrode to the bubble interface, interfacial mass transfer, surface tension variation at the interface, bubble growth and detachment. The model shown in the diagram in Fig. 43 presents the elements necessary to correctly describe an electrogenerated bubble.

• Errors due to assumed known parameters: it was noted that depending on the references considered ∂γ ∂T could vary by a factor 10. From this point of view, it is likely that the thermocapillary effect is badly evaluated. This is why evaluating the surface tension variation on the interface excluding other sources of effect may lead to a wrong interpretation.

• Errors in direct or numerical calculations: in order to limit the direct calculation errors, several numerical methods have been tested (height functions), the aim being to estimate these errors and their impact on the simulation result. The aim has been to limit these errors due to numerical calculation. As a reminder, in the case of a moving bubble, the error due to parasitic currents was limited to 2.5% for a ratio R/∆x = 70.

In conclusion, this thesis :

• presented and discussed the assumptions and structural choices of a holistic model for simulating electrogenerated bubbles;

• provided a numerical model for simulating the Marangoni currents around an electrogenerated bubble;

• evaluated and reduced the errors (spurious currents) of this model using innovative methodologies;

• established a new relation allowing to evaluate the growth rate of a bubble as a function of the intensity of the Marangoni effect around this bubble;

• provided a tool allowing to evaluate the interfacial mass transfer through a mass transfer coefficient as a function of a Marangoni number.

Fig. 1 :

 1 Fig. 1: Definition of bubble coverage, the grey surface represents the normal projection of the bubble on the electrode surface.

Fig. 2 :

 2 Fig. 2: Extract from [Vogt, 2011b], gas evolution efficiency

Fig. 3 :

 3 Fig. 3: Diagram inspired by[START_REF] Blander | Bubble nucleation in liquids[END_REF] showing the work required to form a spherical bubble according to its radius. The system seeking to achieve a minimum energy state below a radius R c the bubble cannot form, above R c the bubble can grow.

Fig. 4 :

 4 Fig. 4: Local microprocesses at gas-evolving electrodes and their influence on mass transfer

Fig. 5 :

 5 Fig. 5: Penetration Theory

Fig. 6 :

 6 Fig. 6: Extract from [Yang et al., 2018], particle trajectories and corresponding velocities aroundthe growing bubble at a potential of -8V .

Fig. 7 :

 7 Fig. 7: Influence of surface tension variation on mass transfer

Fig. 8 :Fig. 9 :

 89 Fig. 8: Extract from [Yang et al., 2015], bubble geometry

Fig. 10 :

 10 Fig. 10: Extract from [Sakuma et al., 2014], bubble growth model for the absorption of the dissolved oxygen gas on (A) hydrophobic and (B) hydrophilic electrode surfaces

Fig. 11 :

 11 Fig. 11: Figure and comments extracted from[START_REF] Wilt | Nucleation rates and bubble stability in water-carbon dioxide solutions[END_REF] -"Definition of geometrical parameters for heterogeneous nucleation in the cases of conical and spherical cavities and projections."

Fig. 12 :

 12 Fig. 12: Figure and comments extracted from[START_REF] Lubetkin | The motion of electrolytic gas bubbles near electrodes[END_REF] : "A hypothetical composite conical or crack-like nucleation site. The region A is hydrophobic, possibly as a result of capillary condensation of a volatile organic material. This deposition mechanism would restrict its presence to the narrowest region of the pore. Nucleation is easy in this environment, with its high contact angle. As the nascent bubble grows past the hydrophobic region, and into the hydrophilic region B, a small remnant may be detached, remaining in the hydrophobic region. This residual gas would promote the instantaneous growth of the next bubble. As it grows, it will be shot out of the conical site, its zero contact angle ensuring the absence of friction. The rapidity of its departure might cause sufficient stirring in the vicinity C, to disrupt the Nernst layer, thus eliminating the concentration Marangoni effect, which would otherwise tend to hold the bubble near the electrode until it was larger. Alternatively, the rapid growth of the first bubble during its rise, and the immediate following of others might cause depletion of the Nernst layer in a narrow bubble chimney, again destroying the concentration gradient."

Fig. 13 :

 13 Fig. 13: Diagram showing the errors that can influence the inverse process.

Fig. 14 :

 14 Fig.14: Illustration of the different continuum methods for describing the evolution of deforming interfaces extract from[START_REF] Wörner | Numerical modeling of multiphase flows in microfluidics and micro process engineering: a review of methods and applications[END_REF]

Fig. 15 :

 15 Fig. 15: The jump operator describes the passage of a quantity Q through an interface between two distinct volumes.

Fig. 16 :

 16 Fig. 16: Diagram of the surface gradient operator ∇ s Q. It is defined as the projection of the gradient ∇Q onto the surface.

Fig. 17 :

 17 Fig. 17: Diagram of a fluid particle. The two gaseous V g and liquid V l volumes are separated by an interface I.

Fig. 18 :

 18 Fig. 18: Example of a 2D simulation with Ansys Fluent of a bubble within a stagnant liquid in a zero-gravity system. The classical CSF formulation was used to calculate the surface tension forces. Spurious currents appear at the interface.

Fig. 19 :

 19 Fig. 19: Methodology of height functions. The sum of the volume fractions present in a column is equal to the average value of a function f on an interval ∆x of which the curve represents the interface.

Fig. 20 :Fig. 21 :Fig. 22 :

 202122 Fig. 20: Stencil adapted according to the interface

Fig. 23 :

 23 Fig. 23: Extracted from [Owkes and Desjardins, 2015] , example of mesh-decoupled columns and heights used to compute the interface curvature.

Fig. 24 :

 24 Fig. 24: Diagram showing the interface. In the columns and rows of the mesh the height functions are estimated.The local coordinate is identified by the two vectors n and

Fig. 25 :

 25 Fig. 25: a : example gradient of ∂γ ∂s , when switching from one to another direction the values of ∂γ ∂s

Fig. 26 :

 26 Fig. 26: 40 × 40 structured mesh used for the simulation

Fig. 27 :Fig. 28 :

 2728 Fig. 27: Algorithm used by Fluent to solve the Navier-Stokes equations. Integration of the DEFINE ADJUST macro to calculate the source term f γ

Fig. 29 :

 29 Fig. 29: Calculation of the exact volume fraction

Fig. 30 :

 30 Fig. 30: Relative curvature errors with the height function methodology along the circular interface as a function of the interface inclination angle with respect to the horizontal.

Fig. 31 :

 31 Fig. 31: Relative curvature errors along the circular interface as a function of the interface inclination angle with respect to the horizontal for two interface resolutions.

Fig. 32 :

 32 Fig. 32: Comparison between analytically calculated Young-Laplace pressure , and numerically evaluated pressures around the bubble for the continuous surface force (CSF) model and the height function (HF) model. The points represent the average pressure at the center of the cells and the x-axis represents the distance from the center of the bubble.

Fig. 33 :Fig. 34 :

 3334 Fig. 33: Evolution of the maximum intensity of the spurious currents observed around the bubble.With the use of an exact curvature for the simulation, the equilibrium is reached for a time t ν .

Fig. 35 :

 35 Fig. 35: Spurious current generated using two mesh densities. On the left R/∆x = 10 on the right R/∆x = 20

Fig. 36 :

 36 Fig. 36: Relative interfacial area errors averaged along a circular interface as a function of n L and the mesh refinement.

Fig. 37 :

 37 Fig. 37: Maximal relative interfacial area errors along a circular interface as a function of n L and the mesh refinement.

Fig. 38 :

 38 Fig. 38: Surface gradient calculation. The interface of the bubble is exposed to a temperature difference along the interface. The value calculated with Eq. (3.79) is compared with the exact value.

Fig. 39 :Fig. 40 :

 3940 Fig. 39: Relative error of the surface tension gradient as a function of mesh resolution. The maximum E max and average E avg error are calculated for two values of the surface tension gradient, 0.1N • m -2 and 50N • m -2 .

Fig. 41 :

 41 Fig. 41: Numerical modelling made with Ansys Fluent of a bubble exposed to a surface tension gradient. The image on the left shows the velocity vectors of the current generated by the Marangoni effect using the method described in this thesis. The Marangoni currents in the image on the right are obtained using the CSF methodology.

Fig. 42 :

 42 Fig. 42: Gradient of a dissolved species initiated along the interface and the origin of the soluto capillary effect shown in Fig.41. The mesh is represented by the grey grid with R/∆ x = 20.

Fig. 43 :

 43 Fig. 43: Diagram representing the different couplings between the different relationships involved in the production of electrogenerated bubbles.

Fig. 44 :

 44 Fig. 44: Diagram extract from Hossain et al. showing the flows around the bubble as a function of electrode size [Hossain et al., 2020].The red lines to the left of each bubble represent the electrical current lines. The narrowing of these electrical current lines makes it possible to determine the position of a hot spot which is at the origin of a temperature gradient and consequently of a surface tension gradient. The black lines to the right of the bubbles represent the lines of fluid movement generated by the surface tension gradient.

  influence on the fluid currents around the bubble. It is worth recalling here the study by Hossain et al. on the influence of electrode size on the thermo-Marangoni effect [Hossain et al., 2020].Fig.44summarises a key point of this study. The temperature gradient along the interface depends on a hot spot created by the Joule effect. The more the electric current lines tighten around the bubble, the greater the Joule effect. The profile of the electric current lines will depend on the size of the electrodes. For micro electrodes the hot spot is at the foot of the bubble. For larger electrodes the point where the current lines narrow is near its equator. Around this hot spot two vortices form. The hot spot reduces the surface tension. The fluid along the interface moves from areas of low surface tension to areas of high surface tension. The results obtained by Hossain et al. exclude other sources of the Marangoni effect[START_REF] Hosokawa | Evaluation of adsorption of surfactant at a moving interface of a single spherical drop[END_REF]. Although these results do not capture the complexity of the real phenomenon, they clarify the influence of the Joule effect on the Marangoni effect, and allow to put into perspective the experimental results obtained from a microelectrode compared to those obtained from a larger electrode.

Fig. 45 :

 45 Fig. 45: Diagram of an electrogenerated bubble on a microelectrode, coated with surfactants. The angle θ S is used to define the free surface of surfactants. Surfactants bind to the top of the bubble and prevent interfacial mass transfer.

  . and then adopted by other authors[Yang et al., 2018; Massing et al., 2019; Hossain et al., 2020; Meulenbroek et al., 2021]. The temperature gradient will depend on the size of the Tab. 5.2: Comparaison of some key references -1 • m -1 5.5 × 10 -5 6.5 × 10 -5 1.6 × 10 -4 electrode (micro or conventional electrode), the bubble coverage, and the current density. The observations of Yang et al. allow to estimate ∆T ≈ 10[K] for a microelectrode [Yang et al., 2018]. In their work Hossain et al. obtained similar results for a microelectrode.

  mol • m -4 4 × 10 3 4 × 10 5 2 × 10 7

  neglect the soluto Marangoni effect due to absorbed dissolved species [Hossain et al., 2020; Massing et al., 2019; Meulenbroek et al., 2021]. This assumption greatly simplifies the model. The assumptions outlined above are shown in the graph in Fig. 46. The model will take into account 3 input parameters :

Fig. 46 :

 46 Fig. 46: Simplified model of the dimensionless study. The model depends on 3 input parameters, the interface length L interf ace over which the interfacial mass transfer and the Marangoni effect takes place, the curvature of the interface κ, and the surface temperature gradient∇ s T .

  is thus found in the transport equation and in the momentum jump relation: m ṽ + pI • n Iof the mass transfer coefficient and the Sherwood number as a function of the Marangoni number

Fig. 47 :

 47 Fig. 47: Diagram explaining the correspondence between the experimental observation and the numerical modelling carried out. In each case the interface is exposed to the same assumed surface tension gradient.

Fig. 48 :Fig. 49 :

 4849 Fig. 48: Numerical modelling made with Ansys Fluent of a bubble exposed to a gradient of an extensive quantity (∆T = 0.1K). The first image on the left shows the initial gradient applied. The second image shows the evolution of this gradient due to the Marangoni effect. The image on the right shows the velocity vectors of the current generated by the Marangoni effect.

Fig. 50 :

 50 Fig. 50: Average mass transfer coefficient over a length L as a function of bubble radius R and Marangoni number M a.

Tab. 5 . 6 :

 56 Data from the point identified on the graphs in Fig.49and Fig.50Ratio of interface length to radius of curvature 0.7D T [m 2 • s -1 ] Thermal diffusivity 1.39 • 10 -7 µ [kg • m -1 • s -1 ] Dynamic viscosity 1.02 • 10 -3 ∆γ ≈ ∂γ ∂T ∆T [N • m -1 ]Surface tension change 2.84 • 10 -3

Fig. 51 :

 51 Fig. 51: Figure extracted from [Afkhami and Bussmann, 2008], using phantom cells below the electrode position it is possible to calculate the curvature at the contact line from height functions.

Glas and West- water, 1964; Fernández et al., 2014; Sakuma et al., 2014; Yang et al., 2015; Liu et al., 2016; Wang et al., 2016; Van Der Linde et al., 2018]. Processes driving bubble growth

	Important aspects of the growth are not yet fully understood. According to Wang in the case
	of electrolysis, growth is controlled by 3 processes [Wang

et al., 2016] :

al., 2012; Luo and White, 2013; Chen et al., 2014; Fernández et al., 2014; Sakuma et al., 2014; Liu et al., 2015; Yang et al., 2015; Baczyzmalski et al., 2016; Karnbach et al., 2016; Liu et al., 2016; Massing et al

., 2019] conducted on bubble use microelectrodes.

; Popinet, 2009; Abadie et al., 2015; Abu-Al-Saud et al

., 2018]. When discretizing and choosing the surface tension model, one must make sure of :

; Cummins et al., 2005; Popinet, 2009; Guo et al., 2015; Magnini, 2016a; Owkes and Desjardins, 2015].

  It can be integrated into surface tension models such as the Brackbill's model. An interface can always be described locally as a graph of a function. The principle of the height function method is to use a local coordinate system to be able to find the curvature of the interface. By calculating the integral of this function and dividing it by the interval over which it has been integrated we obtain the mean value or height of this function over the interval.

  Table 4.1. As noticed by Lubetkin the surface

	Tab. 4.1: Physical parameters used

  Yang et al. found that the thermal and solutal effect had an amplitude of the same order of magnitude[START_REF] Yang | Marangoni convection at electrogenerated hydrogen bubbles[END_REF]. Massing with a microelectrode in a 1 M H 2 SO 4 solution found that the solutal effect was negligible[START_REF] Massing | Thermocapillary convection during hydrogen evolution at microelectrodes[END_REF]. The study by Meulenbroek et al. also excludes the soluto-Marangoni effect caused by dissolved gases but introduces that caused by surfactants not being absorbed by the bubble but on the contrary remaining fixed on the interface and partially blocking the interfacial mass transfer[START_REF] Meulenbroek | Competing Marangoni effects form a stagnant cap on the interface of a hydrogen bubble attached to a microelectrode[END_REF]. This uncertainty cannot be resolved if the effects are not studied simultaneously.

  As shown in Tab. 5.3, Massoudi established how the surface tension of water varies as a function of the partial pressure of different dissolved gases[START_REF] Massoudi | Effect of pressure on the surface tension of water. Adsorption of low molecular weight gases on water at 25.deg[END_REF]. At Linear approximation of the surface tension variation for different gases, γ = γ 0 + bp + cp 2 + dp 3 with γ in mN/m, p in bar, and γ 0 = 71.98 mN/mWhere k H the constant of proportionality is dependent on the temperature and pressure, but in our case can be estimated as k

	Tab. 5.3: Gases	b	c	d
	H 2	-0.025	-	-
	O 2	-0.0779 +0.000104	-
	CO 2	-0.7789 +0.00543 -0.000042
	low pressures, the concentration of dissolved hydrogen and the partial hydrogen pressure
	can be related through Henry's law:
		c h 2 = p k H	(5.16)
		∂γ ∂c	∂c ∂z	(5.15)

H = 7.8 × 10 -6 mol m -3 P a -1

[START_REF] Wiebe | The Solubility of Hydrogen in Water at 0, 50, 75 and 100°from 25 to 1000 Atmospheres[END_REF][START_REF] Sander | Compilation of Henry's law constants (version 4.0) for water as solvent[END_REF]

. As a result we obtain :

al., 2019; Meulenbroek et al., 2021].

  Based on the velocity field observed by Massing et al. and numerically recovered by Meulenbroek et al. it is possible to establish a relationship that allows the growth ratedR(t) dt to be recovered[Massing et

Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022 Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022

Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022

NomenclatureStudy of electrogenerated two-phase and microfluidic flows Florent Struyven 2022

Chapter 1Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022

The electrolysis process 13 Study of electrogenerated two-phase and microfluidic flows Florent Struyven

The electrolysis process

Chapter 2Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022

Bubble development

Transport of the diluted species in the vicinity of the bubble

Interfacial mass transfer

Microfluidic phenomena in electrolysis

Microfluidic phenomena in electrolysis 57 Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022

How to improve electrodes?

Integration into an inverse resolution problem

Overview of the chapter and objectives of the thesis

Chapter 3Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022

Toward a direct numerical simulation of the phenomenum

Mathematical model

From continuous to discrete model

Numerical errors

Surface tension model 91 Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022

Surface tension model

Chapter 3

Marangoni model 103 Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022

Marangoni model

Numerical model mass transfer across the interface

Integration of the code in Fluent

Presentation of the code used in the simulations

Chapter 4Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022

Curvature calculation errors

Chapter 4Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022

Static bubble test case

Static bubble test case 123 Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022

Chapter 4Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022

Interfacial area error calculation

Surface gradient error calculation

Chapter 4Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022

Translating bubble test case 129 Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022

Translating bubble test case 131 Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022

Conclusion of the chapter

Conclusion of the chapter 135 Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022

Chapter 5Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022

Heat, transport, and surfactant models

Couplings

Competing Marangoni effect

Bubble growth rate with Marangoni Flow

Non-dimensional study

Conclusion of the chapter and discussion of model uncertainties

Chapter 6Study of electrogenerated two-phase and microfluidic flows Florent Struyven 2022

Perspective: integration of a moving contact line model

Conclusion of the work

BibliographyStudy of electrogenerated two-phase and microfluidic flows Florent Struyven 2022

Fixed interface hypothesis

A first consideration is to examine the time scales considered. The time corresponding to the bubble cycle varies according to the studies from 2 to 5s [START_REF] Sakuma | Nucleation and growth of electrolytic gas bubbles under microgravity[END_REF][START_REF] Yang | Dynamics of Single Hydrogen Bubbles at a Platinum Microelectrode[END_REF][START_REF] Guo | Implementation of a height function method to alleviate spurious currents in CFD modelling of annular flow in microchannels[END_REF][START_REF] Yang | Marangoni convection at electrogenerated hydrogen bubbles[END_REF][START_REF] Massing | Thermocapillary convection during hydrogen evolution at microelectrodes[END_REF]. This time must be compared with the evolution of the temperature and velocity profile near the interface, which evolve quite rapidly in less than a second. Data on the thermal properties of the electrolyte when growing an electrogenerated bubble on a microelectrode according to the experiments by Massing et al. are given in Tab. 5.1 [Massing et al., 2019]. The value of the radius R = 560 µm chosen by these authors for their numerical simulations and corresponding to an advanced stage of the bubble's development serves as a basis for the following reasoning.

Symbol

Descritpion Value

Thermal conductivity 0.58

Taking the radius of the bubble as the characteristic length, and starting from the first instant when the Joule effect starts to modify the temperature of the electrolyte, theoretically close to the ambient temperature, the heat will diffuse within the electrolyte for a time

The heat is thus diffused slowly with respect to the bubble cycle. However, heat is also advected by the thermocapillary movement of the electrolyte at the interface. The ratio of convective to diffusive heat transfer is expressed by the Péclet number. Based on the results of Yang et al. a mean velocity u = 10 mm s -1 is a good compromise for estimating the convective motions within the electrolyte. The Péclet number is expressed as: (5.10) where D th = k ρ cp is the thermal diffusivity. Convective heat transport is therefore much more important than diffusive transport. Another time scale to consider is the recirculation time of the electrolyte within the vortex near the bubble. A characteristic size of this vortex of 200 µm seems to be a good approximation in view of the observations made by [START_REF] Yang | Marangoni convection at electrogenerated hydrogen bubbles[END_REF]. By assimilating this vortex to a circle, we obtain a a characteristic recirculation time within the vortex of 0.12 s.