# Eulerian Two-Fluid Model of Alkaline Water Electrolysis for Hydrogen Production

Damien Le Bideau, Philippe Mandin, Mohamed Benbouzid, Myeongsub Kim, Mathieu Sellier, Fabrizio Ganci, Rosalinda Inguanta

► To cite this version:

Damien Le Bideau, Philippe Mandin, Mohamed Benbouzid, Myeongsub Kim, Mathieu Sellier, et al. Eulerian Two-Fluid Model of Alkaline Water Electrolysis for Hydrogen Production. Energies, 2020, 13 (13), pp.3394. 10.3390/en13133394. [hal-05021057](https://hal-05021057)

HAL Id: hal-05021057

https://ubs.hal.science/hal-05021057v1

Submitted on 8 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Eulerian Two-Fluid Model of Alkaline Water Electrolysis for Hydrogen Production

Damien Le Bideau<fnref n="1" />, Philippe Mandin <fnref n="1" />*, Mohamed Benbouzid <fnref n="2" /> ID, Myeongsub Kim <fnref n="3" />, Mathieu Sellier ⁴ID, Fabrizio Ganci ⁵ID and Rosalinda Inguanta ⁵ID

<figure><img src="image_2.png" /></figure>

check for
updates

Received: 1 June 2020; Accepted: 30 June 2020; Published: 2 July 2020

**Abstract:** Hydrogen storage is a promising technology for storage of renewable energy resources. Despite its high energy density potential, the development of hydrogen storage has been impeded, mainly due to its significant cost. Although its cost is governed mainly by electrical energy expense, especially for hydrogen produced with alkaline water electrolysis, it is also driven by the value of the cell tension. The most common means of electrolyzer improvement is the use of an electrocatalyst, which reduces the energy required for electrochemical reaction to take place. Another efficient means of electrolyzer improvement is to use the Computational Fluid Dynamics (CFD)-assisted design that allows the comprehension of the phenomena occurring in the electrolyzer and also the improvement in the electrolyzer's efficiency. The designed two-phase hydrodynamics model of this study has been compared with the experimental results of velocity profiles measured using Laser Doppler Velocimetry (LDV) method. The simulated results were in good agreement with the experimental data in the literature. Under the good fit with experimental values, it is efficient to introduce a new physical bubble transfer phenomenon description called "bubble diffusion".

**Keywords:** hydrogen production; alkaline water electrolysis; two-phases flow; CFD; two-phase process

## 1. Introduction

A key challenge of the 21st century is to deal with climate change. The Intergovernmental Panel on Climate Change (IPCC) declares that, to stay below the 3 °C of global temperature increase from the pre-industrial era, zero net emission of greenhouse gases must be achieved at the mid-21 century. Thus, the consumption of fossil-based fuels must be reduced from the current level as much as possible. The present global energy demands can hardly be met simply by replacing conventional fossil-based energy sources (e.g., thermal or nuclear energy) by renewable ones. The integration of renewables into the national electrical grid brings the issue of energy storage because of the intermittent nature of renewable energy resources. Energy storage is achieved using several processes such as battery, Pumped Hydroelectric Energy Storage (PHES) or super capacitor. However, only a few processes allow interseasonal storage (synthetic natural gas and hydrogen (H₂)). There are several processes that do not emit greenhouse gases and produce H₂, such as thermochemical hydrogen process and electrolytic

<footnote n="1">¹ Institut de Recherche Dupuy de Lôme (UMR CNRS 6027 IRDL), University Bretagne Sud, 56100 Lorient, France; damien.le-bideau@univ-ubs.fr</footnote>

<footnote n="2">² Institut de Recherche Dupuy de Lôme (UMR CNRS 6027 IRDL), University of Brest, 29238 Brest, France; Mohamed.Benbouzid@univ-brest.fr</footnote>

<footnote n="3">³ Department of Ocean and Mechanical Engineering, Florida Atlantic University, Boca Raton, FL 33431, USA; kimm@fau.edu</footnote>

<footnote n="4">⁴ Department of Mechanical Engineering, University of Canterbury, Christchurch 8140, New Zealand; mathieu.sellier@canterbury.ac.nz</footnote>

<footnote n="5">5 Dipartimento di Ingegneria, University of Palermo, Viale delle Scienze, 90128 Palermo, Italy; fabrizio.ganci@unipa.it (F.G.); rosalinda.inguanta@unipa.it (R.I.)</footnote>

<footnote>* Correspondence: philippe.mandin@univ-ubs.fr</footnote>
H₂ process. Thermochemical processes for hydrogen production use heat to decompose water into hydrogen and oxygen [1], whereas electrolytic H₂ processes use electric energy to achieve water splitting. Electrolytic hydrogen is a promising storage technology due to its high specific potential energy and storage time, but its cost hampers its development. Another interesting advantage of hydrogen is the diversity of its use. Indeed, it can be used as a fuel for fuel cell for electricity production for buildings, as a reactant in chemical processes and as a fuel for powering vehicles. Thus, using hydrogen, especially for fueling heavy vehicles, could be one of the solutions to reduce the dependence on fossil fuels in the transport sector. There are three leading electrolysis technologies, but only two have been employed in industrial applications: they are the alkaline and Proton Exchange Membrane (PEM) electrolyzer. The latter is more efficient than the former but, due to the requirement of acidic electrolyte, rare materials must be used (such as platinum). Alkaline water electrolyzer has several advantages over PEM electrolyzer such as its robustness, cheaper costs, and system mature. One way of decreasing the cost is to increase the efficiency of the process by decreasing the cell voltage. The cell voltage is a sum of overvoltage Equation (1).

$$U_{cell}(j) = E_{rev}(P, T) + \eta(j)_{cath act} + \eta(j)_{an act} + R(T, Y_{KOH}) j + \eta(j)_{cath\backslash an conc} \quad (1)$$

$E_{rev}$ is the reversible voltage in V, $\eta_{act}$ is the activation overvoltage in V, R the sum of the electrodes, membrane and electrolyte resistance in Ω cm² and $\eta_{conc}$ the concentration overpotential in V, j the current density in A m⁻².

The $E_{rev}$ voltage represents the minimum voltage where the water electrolysis occurs. This value is around 1.23 V in terms of atmospheric pressure and 25 °C. This minimal voltage can be decreased by increasing the temperature. The majority of previous studies have focused on finding a cheap, robust and electroactive electrode material to decrease the activation overpotential, which is the second term of the Equation (1) [2–4]. There is also another research trend that suggests that another way of increasing the efficiency is to decrease the ohmic resistance through the electrolyzer design using simulations or experiments. Another interesting research trend is to find another redox couple, such as those that exist in urea electrolysis [5] or through the use of an electrocatalytically driven cyclic process involving iron oxide species that allow a low OER overvoltage [6]. The latter research trend tries to decrease the cell voltage using material science. The present approach uses the chemical engineering approach. Even if it is theoretically possible to perform monophasic alkaline water electrolysis [7], this process is performed under two-phase flow configuration in the industry. Thus, the bubble impact on the efficiency must be properly modeled or simulated. The first well-known study about this subject was performed by Tobias et al. [8] and Hine et al. [9]. Both studies focused on the bubble effect on the electrolyte resistance. More recently, the simulation has been used to study this two-phase flow phenomenon. There are two main models: Euler–Euler and Euler–Lagrange models. The latter considers the electrolyte as a continuum phase and the bubbles as a discrete phase. Mandin et al. [10] used this model in their study, and the bubble dispersion is taken into account using a horizontal source term. Hreiz et al. [11] also used this model. In order to prevent the use of the source term, the bubble injection has been shifted away from electrode surfaces. Using this technique, it has been concluded that the drag force was sufficient to describe the bubble–liquid interaction. However, the presented comparison with the experiments in their study was mainly qualitative. For Euler–Euler models, two types of models have been used: the mixture model and the two-fluid model. The former has been employed by Dahlkild [12], Wedin et al. [13] and Schillings [14]. In those three studies, the model used is called semi-empirical because an empirical correlation is chosen. In Wedin et al. [13] and Schillings et al. [14], the calculations are compared to the experiments performed by Boissonneau et al. [15]. Those models prevent the use of turbulence closing model. However, Boissonneau et al. [15] observed a bubble-induced turbulent behavior in the upper part of the narrow and small electrochemical cell. In their study, Aldas et al. [16] used a laminar two-fluid model and found that this model underestimates the void fraction distribution compared to the experiment. They concluded that local weak turbulence must be taken into account. In another
study, Mat [17] developed a two-fluid model that takes into account local turbulence. The computed results were compared with the experimental void fraction distribution results of Riegel [18] and good agreements were found. In this study, a two-fluid Euler–Euler is developed and will be compared with the experimental results of Boissonneau et al. [15] and the numerical results of Schillings et al. [14]. In this study, a new force is added to the traditional Fluent© commercial software Eulerian model. The pertinence of its use from a numerical and theoretical point of view will be discussed.

## 2. Numerical Methods

This section describes the geometry and the type of mesh used in the current study. The mathematical formulation, the boundary conditions and the numerical procedure are also introduced. Comparison of the numerical results with the experiments and discussion will be followed in the next section.

### 2.1. *Geometry and Mesh*

The geometry of the computational domain is identical to the experimental apparatus of Boissonneau et al. [15]. Figure 1 shows experimental apparatus of Boissonneau et al. [15]. It is composed of a channel with a height of 120 mm, a width of 3 mm and a thickness of 30 mm. The whole channel is submerged in the electrolyte. The electrodes are 40 mm high and are placed 40 mm away from the bottom. When the current is applied, bubbles are generated on electrode surfaces. The detached bubbles rise and trigger the electrolyte causing the bubble-induced spontaneous movement of the surrounding electrolyte. This phenomenon is called gas-lift configuration.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Geometry of the computation domain. h stands for the half-width of the electrolyte, H for the electrode height. The blue arrow symbolizes the pumping electrolyte induced by the electrogenerated bubbles.</figcaption></figure>

In numerical calculations, a grid sensitivity study must be performed in order to obtain a grid-independent solution. In this study, the grid is refined until the difference between two solutions reaches zero. On the other hand, in the 2D two-fluid model, this is not always possible. Indeed, Picardi et al. [19] suggested that there is an ideal grid resolution for the 2D two-fluid model simulation. They determined that the grid resolution must correspond to the bubble-to-cell ratio $1/\sqrt{2}$. If a finer
grid is chosen, the results become non-physical or the calculations diverge. The same results have been
reproduced by Law et al. [20]. Finally, Picardi et al. [19] stated that this problem could be resolved
when a 3D geometry was employed. However, their calculation took more time than in 2D (from
10 CPU hours to 94 days using a single processor on a quad-core 2.26 GHz Intel Nehalem processor).
Panicker et al. [21] and Vaidheeswaran et al. [22] explained this phenomenon by the fact that the
problem is ill-posed. An ill-posed problem means that the solution is elliptic and the solution of the
Nth step depends on the N+1th step. They solved the problem by adding a collision or dispersion term.
They also affirmed that although the virtual mass force has an insignificant influence on the results,
it ensures a better stability on the calculation. Panicker et al. [19] declared that other authors solved this
problem by artificially increasing the liquid viscosity or adding an interface pressure term to the model.
Both articles used a linear stability analysis and, using the results of this study, Panicker et al. [19]
added a dispersion term that uses the void fraction gradient. To ensure the hyperbolicity, a parameter
depending on the void fraction was added to the dispersion term. As a result, the resulting model gives
better results when compared with the model without the dispersion term. In this study, we choose
to add the dispersion term used in studies by Marfaing et al. [23] and Davidson [24]. The Figure 2
shows the results of the grid sensitivity with four different meshes: 100, 120, 150 and 200 µm. It can
be noted that there is a discrepancy of around 10% in the all domains with the grid size between a
200 and 150 µm. The maximum difference between the 150 µm and the other grids are located at the
velocity and void fraction peaks. Even if it seems that 100 µm is sufficient to describe the void fraction
and velocity distribution. a mesh of 60 µm has chosen in order to be sure that the results converge
with no grid-dependency.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Results of the mesh sensitivity study. The graph (a) presents the mesh dependence of the velocity. The graph (b) presents the mesh dependence of the void fraction.</figcaption></figure>

## 2.2. Mathematical Formulation

### 2.2.1. The Bubble Dispersion Problem

Lee et al. [25], Abdelouahed et al. [26] and Hreiz et al. [11] experimentally visualized and measured
the void fraction distribution in the electrolyte in a cell under forced convection [25] and not net flow
configuration cell (in this type of configuration, the gases escape at the air–liquid interface, but the
electrolyte stay confined into the inter-electrode gap) [11,26]. They all reported that the bubbles were
spread into the electrolyte. After simulating their experiments, Abdelouahed et al. [26] observed that
classical models of lift (e.g., Saffman-Mei) and drag failed to reproduce their experimental results.
They attributed this spreading to a lift force that has a negative coefficient. Although their model fits
their experimental data, the model seems to rely on numerical diffusion triggered by a coarse mesh, so its use is questionable, at least in our present configuration, due to the fact that the electrochemical cell is narrow. This is because, to obtain their results, a coarse mesh must be used, but for our present geometry, this coarse mesh does not satisfy all the mesh independence conditions (the reader can refer to Section 2.1). Hreiz et al. [11] attributed this spreading to the numerical diffusion. Hreiz et al. [11] modeled the same experimental study using a Euler–Lagrange model (gas is modeled as a discrete phase). They succeeded in quantitative reproduction of their experimental data by using the drag force only, but the mesh independence condition was not satisfied. The model of Mat [17] is a Euler–Euler model that succeeds in simulated results that fit experimental data, but the results seem odd. They measured and calculated the diphasic boundary layer and void fraction increase with increasing electrolyte flow. Those results have been invalidated experimentally by Lee et al. [25] and numerically by Schillings et al. [14]. In addition, an increase in void fraction value triggers an increase in the cell voltage and it has been measured that increasing the electrolyte flow decreases the cell voltage. Schillings et al. [14] and Wedin et al. [13] used a mixture model developed by Ishii that simulated results close to the experimental data from Boissonneau et al. [15]. In addition to the drag force and Saffman lift force, their closure term is composed of three terms of bubble–bubble interactions. Thus, in this study, we decided to introduce a bubble dispersion term. Indeed, a high concentration of bubbles increases the bubble collision and hence trigger of bubble diffusion from high concentration to low concentration. This phenomenon has been observed by Ham et al. [27]. The additional force used in this study is presented in the Equation (2). As described earlier, the additional force is inspired by the force used in Marfaing et al. [23] and Davidson [24]. However, in their model, the present term was multiplied by the drag coefficient. The void fraction gradient is used to mathematically traduce the diffusive nature of this force.

$$ \vec{F}_{BD} = - \underbrace{\varepsilon_g \rho \frac{K_g}{d_b} |U_r| \nabla \varepsilon_g}_{\text{Bubble dispersion force}} \qquad (2) $$

$\varepsilon$ is the gas or liquid fraction, the subscript k can be either g (O₂, H₂) or liq (for liquid), ρ is the density in kg m⁻³, d_b is the bubble diameter in m, U_r is the gas phase velocity minus the liquid velocity in m s⁻¹.

Using this term is expected to reproduce the turbulence-like behavior of the electrolyte flow observed by Boissonneau et al. [15].

### 2.2.2. Model

The model has been designed using the following hypotheses:

* The flow is Newtonian, viscous and incompressible;

* Due to the high heat transfer induced by the bubbles and the relatively low surface and volume inducing a low heat injection by ohmic and activation losses, the flow is considered isothermal

* At the same time, numerical simulations were carried out in order to highlight the influence of ions on the velocity and void fraction distribution. There were only very little differences when the ions distribution was taking into account. Thus, the electrolyte is considered as extremely well mixed. This hypothesis has been made also by Abdelouahed et al. [26] and Schillings et al. [14];

* Oxygen, hydrogen and electrolytes are three continuum media;

* Because the Reynolds number calculated via the data from Boissonneau et al. [15] is between 240 and 480, the flow is considered as laminar;

* The effect of the surface tension is neglected;

* Nagai et al. [28] observed a dependency of the diameter with respect with the current density and in addition to this dependency, Boissonneau et al. [15] observed that the bubble diameter increases with height. Schillings et al. [14] chose to make the diameter increase with the height and current density. In the current study and in other studies [12,13], it was decided that a constant bubble diameter is taken for a given current density. Although this hypothesis is made in a lot of studies,
it does not reflect the reality. There is a distribution of the bubble diameter, but since it is not yet possible to model this distribution, we are forced to make this assumption;

* The current density distribution does not affect the flow distribution [29,30]. Thus, the current density distribution is taken as uniform for the validation.

The accuracy of the two-dimensional (x,y) flow hypothesis needs to be discussed. Therefore, for each phase, the equation set can be written

$$ \frac{\partial \epsilon_k \rho_k}{\partial t} + \vec{\nabla} \cdot (\epsilon_k \rho_k \vec{V}_k) = S_g \quad (3) $$

$\epsilon$ is the gas or liquid fraction, the subscript k can be either g (O₂, H₂) or liq, $\rho$ is the density in kg m⁻³, V the velocity in m s⁻¹, and $S_g$ is the term source in kg m⁻³ s⁻¹.

$$ \frac{\partial}{\partial t}(\epsilon_k \rho_k \vec{V}_k) + \vec{\nabla} \cdot (\epsilon_k \rho_k \vec{V}_k \vec{V}_k) = -\epsilon_k \vec{\nabla} p + \vec{\nabla} \cdot (\epsilon_k \tau) + \epsilon_k \rho \vec{g} + \vec{F}_k \quad (4) $$

p is the pressure in Pa, $\tau$ is the stress tensor in Pa, g the gravitational acceleration in m s⁻², and $\vec{F}_k$ is the exchange term in N m⁻³.

The stress tensor is written as follows

$$ \tau = \mu_k \left[ \left( \vec{\nabla} \vec{V}_k + \vec{\nabla} \vec{V}_k^T \right) - \frac{2}{3} \vec{\nabla} \cdot \vec{V}_k I \right] \quad (5) $$

with $\mu_k$ the viscosity of the phase k in Pa s and I the unit tensor.

$$ \vec{F}_k = \vec{F}_D + \vec{F}_L + \vec{F}_{BD} \quad (6) $$

$$ \vec{F}_k = \underbrace{-\frac{3}{4}\epsilon_g \rho \frac{C_D}{d_b} |U_r| U_r}_{\text{Drag force}} - \underbrace{\epsilon_g \rho C_L |U_r| \text{rot}(\vec{V}_l)}_{\text{Lift force}} - \underbrace{\epsilon_g \rho \frac{K_g}{d_b} |U_r| \vec{\nabla} \epsilon_g}_{\text{Bubble dispersion force}} \quad (7) $$

## 2.3. Boundary Conditions

The water electrolysis performed in Boissonneau et al. [15] is neither acidic or alkaline. This electrolysis is called as aqueous in this study. The supporting electrolyte is Na₂SO₄ concentrated at 50 g L⁻¹. Therefore, the reaction occurring at the anode and the cathode are Equations (8) and (9), respectively.

$$ H_2O_{(liq)} \rightarrow 2H^+_{(aq)} + \frac{1}{2} O_2_{(g)} + 2e^- \quad (8) $$

$$ 2H_2O_{(liq)} + 2e^- \rightarrow 2OH^-_{(aq)} + H_2O_{(g)} \quad (9) $$

If the reader wants to simulate an alkaline electrolysis, he must change the thermophysical inputs from Na₂SO₄ to KOH or NaOH. The reader can refer to Le Bideau et al. [31], where the thermophysical inputs are given with their sensitivity with respect to temperature and electrolyte mass fraction.

The quantity of produced gases is directly correlated to the current density through the Faraday's law Equation (10)

$$ q_m = n \frac{j S}{F} M_g \quad (10) $$

$q_m$ is the mass flow of produced gas in kg s⁻¹, $M_g$ is the molar mass of the gas kg mol⁻¹, S is the electrode surface in m², F is the Faraday constant 96,500 C mol⁻¹, and n is the ratio of the stoichiometric number of the gas and the number of electrons exchanged during the reaction.

In the two-fluid equation, in most studies [9–11,22,25,26], the input parameters for the boundary conditions are the velocity and the void fraction. The value of the void fraction is fixed arbitrarily.
According to Alexiadis [32,33], this value does not influence the hydrodynamic characteristic of the flow. However, better results have been obtained using a source term that produces gas in the cell in the vicinity of the electrodes. This method has been used by Charton et al. [34]. This source term is written as Equation (11)

$$S_g = n \frac{j}{F \times \Delta x} M_g \qquad (11)$$

$\Delta x$ is the width of the first cell next to the electrode in m. For the bottom and top boundary condition, a pressure inlet ($P_{Tot} = 0$) and pressure outlet condition ($P = 0$) is fixed. For the other wall, a no-slip condition is fixed meaning that the velocity is set to 0 m s⁻¹. The Table 1 summarizes the boundary conditions for the present problem.

<table><caption>**Table 1.** The boundary conditions to solve the problem.</caption><thead><tr><th>Position</th><th>Boundary Conditions</th></tr></thead><tbody><tr><td>$x = 0$ $H_{elec} < y < 2$ $H_{elec}$</td><td>$\vec{V}_l = \vec{V}_{H2} = 0$</td></tr><tr><td>$x = 2h$ $0 < y < 3$ $H_{elec}$</td><td>$\vec{V}_l = \vec{V}_{O2} = 0$</td></tr><tr><td>$0 < x < 2h$ $y = 0$</td><td>$P_{Tot} = 0$</td></tr><tr><td>$0 < x < 2h$ $y = 3$ $H_{elec}$</td><td>$P = 0$</td></tr></tbody></table>

## 2.4. Numerical Procedure

Equations (3) and (4) are solved using the commercial code, Ansys Fluent. This commercial code solves equations using the finite volume method by discretizing the geometry in volume and subsequently integrating the governing equation over the volume. The governing equation is expressed as the following algebraic Equation (12)

$$a_p \varphi = \sum_{i}^{N} a_i \varphi + b \qquad (12)$$

It is considered that the convergence is met when the residuals remain stable and when the average gas and liquid velocity, as well as the average gas void fraction, reach the value of $10^{-3}$. In order to reach this convergence, for 2D simulation, the flow is initialized with a forced convection by imposing pressure at the bottom. When the convergence is reached, the obtained results are used as initial guess for the bubble-driven flow. The calculation time can vary, depending on the numerical procedure, between 20 and 30 min with an Intel Core i7-6700HQ CPU @2.60 GHz.

Table 2 gives the input data for the validation of the numerical model. It can be seen that the current density varies between 500 and 2000 A m⁻². The most interesting datapoints are those obtained at 2000 A m⁻². Indeed, industrial alkaline electrolyzers generally use a current density between 2000 and 5000 A m⁻² [35]. The input data for the electrolyte were taken from the work of Isono et al. [36]. As the bubble diameter is taken as a constant, an average value has been calculated from the correlation of Schillings et al. [14]. As aforementioned, the term K in the Equation (5) is used to fit the experimental data of Boissonneau et al. [15]. A sensitivity study has been performed to choose the parameter. The results of this study are presented in Table 3. The parameter K is always bigger for oxygen than hydrogen.
<table><caption>**Table 2.** The input values for the problem.</caption><thead><tr><th>Name</th><th>Value</th></tr></thead><tbody><tr><td><strong>Geometry inputs</strong></td><td></td></tr><tr><td>H<sub>elec</sub> (mm)</td><td>40</td></tr><tr><td>L (mm)</td><td>30</td></tr><tr><td>h (mm)</td><td>1.5</td></tr><tr><td>H<sub>Tot</sub> (mm)</td><td>120</td></tr><tr><td><strong>Physical Inputs</strong></td><td></td></tr><tr><td>ρ<sub>1</sub> (kg m<sup>-3</sup>)</td><td>1040</td></tr><tr><td>ρ<sub>O2</sub> (kg m<sup>-3</sup>)</td><td>1.3</td></tr><tr><td>ρ<sub>H2</sub> (kg m<sup>-3</sup>)</td><td>0.08</td></tr><tr><td>v<sub>1</sub> (m<sup>2</sup> s<sup>-1</sup>)</td><td>9.97 × 10<sup>-7</sup></td></tr><tr><td><strong>Two-Phase Inputs</strong></td><td></td></tr><tr><td>d<sub>b</sub> (μm)</td><td>for 500 A m<sup>-2</sup> d<sub>b</sub> = 50 μm<br>for 1000 A m<sup>-2</sup> d<sub>b</sub> = 58 μm<br>for 2000 A m<sup>-2</sup> d<sub>b</sub> = 78 μm</td></tr></tbody></table>

<table><caption>**Table 3.** Sensitivity study results for the parameter K.</caption><thead><tr><th>500 A m<sup>-2</sup></th><th>1000 A m<sup>-2</sup></th><th>2000 A m<sup>-2</sup></th></tr></thead><tbody><tr><td>K<sub>O2</sub>/d = 10.5</td><td>K<sub>O2</sub>/d = 9</td><td>K<sub>O2</sub>/d = 5</td></tr><tr><td>K<sub>H2</sub>/d = 5</td><td>K<sub>H2</sub>/d = 4</td><td>K<sub>H2</sub>/d = 2.5</td></tr></tbody></table>

In order to use the current model to other design, the Vaschy-Buckingham has been used and the sensitivity of the K parameter to dimensionless groups has been calculated in Equations (13)–(16).

$$Re_G = \frac{\rho_1 V_G H_{elec}}{\mu_L} \qquad (13)$$

$$Fr_G = \frac{g H_{elec}}{V_G^2} \qquad (14)$$

$$r^* = \frac{d}{2H} \qquad (15)$$

$$h^* = \frac{h}{H} \qquad (16)$$

Therefore, the following correlation Equation (17) was used to calculate the K parameter. However, this correlation was used for height in the order of 10 cm and for $Fr_G$ numbers higher than $10^5$.

$$\frac{K}{V_G H} = 0.197 Re_G^{0.108} + 0.5 r^{*0.124} + 0.668 h^{*-0.408} + 0.000323 Fr_G^{0.661} + 0.375 \qquad (17)$$

### 3. Results

The experimental data of Boissonneau et al. [15] are the liquid velocity at three locations: 5 mm before the electrodes (y = 35 mm), at the mid-section of the electrodes (y = 60 mm) and 5 mm before the ends of the electrodes (y = 75 mm). Figures 2–4 present the results of the current study and the experimental data. At the entrance of the channel, a Poiseuille liquid velocity distribution is observed, and this distribution is flattened at the center (also called bulk). The exchange of momentum between the gas phases and the liquid phase is easily observed next to the anode and cathode. It can be noted that the peak induced by hydrogen is bigger than the oxygen peak, because the injected volume of hydrogen is two times bigger than oxygen one. The fact is that in the bulk, a plateau is observed due to bubble induced turbulence [15]. Those figures allow for comparison of the current numerical results with experimental data. The “Abdelouahed-like” model is a model that uses the same method
as presented in Abdelouahed et al. [26]. In their study [26], the authors used the lift coefficient to fit their simulated data with their experimental values. They showed that a negative coefficient permits a good agreement between their measurements and predicted values. In the present study, a sensitivity study was performed to have the coefficient that fits the most the experimental data from Boissonneau et al. [15]. First of all, the numerical results for all the models show some discrepancies with the experimental data. Schilling’s model and Abdelouahed-like model accurately predict the velocity distribution on the cathode side, but the accuracy of the simulated data drops on the anode side. The Figure 3 shows that the “Abdelouahed-like” model describes the anode side and cathode side velocity distribution with accuracy, but the bulk velocity distribution shows larger errors. However, as shown in Figures 4 and 5, the more the current density increases, the less precise the model is.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Result of the liquid velocity in different models at *j* = 500 A m⁻² (**a**) at *y* = 60 mm and (**b**) at *y* = 75 mm. The dot are the data from Boissonneau et al. [15], the grey dashed line is the numerical results of the Schilling et al. model [14], the dotted and dashed line is a model using the method of Abdelouahed et al. [26]. Finally, the black solid line the results from the current study with 2D approximation.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Results of the liquid velocity in different models at *j* = 1000 A m⁻² (**a**) at *y* = 60 mm and (**b**) at *y* = 75 mm. The dot are the data from Boissonneau et al. [15], the grey dashed line is the numerical results of the Schilling et al. model [14], the dotted and dashed line is a model using the method of Abdelouahed et al. [26]. Finally, the black solid line the results from the current study with 2D approximation.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Result of the liquid velocity in different models at *j* = 2000 A m⁻² (a) at *y* = 60 mm and (b) at *y* = 75 mm. The dot are the data from Boissonneau et al. [15], the grey dashed line is the numerical results of the Schilling et al. model [14], the dotted and dashed line is a model using the method of Abdelouahed et al. [26]. Finally, the black solid line the results from the current study with 2D approximation.</figcaption></figure>

Figure 6 shows the void fraction evolution predicted by the current model. With increasing current density, the diphasic boundary layer of oxygen and hydrogen and the maximum oxygen and hydrogen void fraction increase. However, there are two parameters (the current density and the bubble diameter) that change between the two cases. Therefore, the calculated evolution cannot be clearly attributed to one parameter. The prediction of this evolution, depending on electrode height, is very important to predict the current density distribution. The prediction of the diphasic boundary layer thickness is also an essential output parameter because the electrolyte conductivity decreases with the void fraction.

Table 4 shows the average relative difference of the different model with the Boissonneau et al. [15] at three current density (500, 1000 and 2000 A m⁻²) at three horizontal zones and two vertical locations (*y* = 60 mm and *y* = 75 mm). The main drawback of the current model is that it does not describe the velocity at the bulk zone (between *x* = 0.5 and *x* = 2.5 mm) with high accuracy. However, for all current density and location, the current model has more accuracy than the Schillings et al. [14] model the velocity in O₂ zone (between *x* = 2.5 and *x* = 3 mm). The present model is less accurate than the model of Schillings et al. However, the model of Schillings et al. uses more semi-empiric term than the model of this study. Moreover, the mathematical formalism of the current study contains less equation than the model of Schillings et al. Thus, this model is easier to program.
<table><caption>**Table 4.** Average relative difference of the different model compared to Boissonneau et al. [15] at two vertical positions 60 and 75 mm and three horizontal zones for the three current density. H₂ is the zone between [0, 500] µm and O₂ is zone between [2.5, 3] mm.</caption><thead><tr><th>j (A m<sup>-2</sup>)</th><th>Model Comparison</th><th>Zone</th><th>60 mm</th><th>75 mm</th></tr></thead><tbody><tr><td rowspan="12">500 A m<sup>-2</sup></td><td rowspan="3">Boissonneau vs. Abdelouahed</td><td>H<sub>2</sub></td><td>6.95 x 10<sup>-2</sup></td><td>3.96 x 10<sup>-2</sup></td></tr><tr><td>O<sub>2</sub></td><td>2.48 x 10<sup>-1</sup></td><td>1.47 x 10<sup>-1</sup></td></tr><tr><td>Bulk</td><td>2.91 x 10<sup>-1</sup></td><td>2.91 x 10<sup>-1</sup></td></tr><tr><td rowspan="3">Boissonneau vs. Schillings</td><td>H<sub>2</sub></td><td>6.00 x 10<sup>-2</sup></td><td>1.60 x 10<sup>-2</sup></td></tr><tr><td>O<sub>2</sub></td><td>2.85 x 10<sup>-1</sup></td><td>2.65 x 10<sup>-1</sup></td></tr><tr><td>Bulk</td><td>1.90 x 10<sup>-1</sup></td><td>2.06 x 10<sup>-1</sup></td></tr><tr><td rowspan="3">Boissonneau vs. Current Study</td><td>H<sub>2</sub></td><td>1.76 x 10<sup>-1</sup></td><td>4.55 x 10<sup>-2</sup></td></tr><tr><td>O<sub>2</sub></td><td>2.01 x 10<sup>-1</sup></td><td>7.40 x 10<sup>-2</sup></td></tr><tr><td>Bulk</td><td>2.48 x 10<sup>-1</sup></td><td>2.36 x 10<sup>-1</sup></td></tr><tr><td rowspan="12">1000 A m<sup>-2</sup></td><td rowspan="3">Boissonneau vs. Abdelouahed</td><td>H<sub>2</sub></td><td>1.04 x 10<sup>-1</sup></td><td>6.13 x 10<sup>-2</sup></td></tr><tr><td>O<sub>2</sub></td><td>3.48 x 10<sup>-1</sup></td><td>2.69 x 10<sup>-1</sup></td></tr><tr><td>Bulk</td><td>3.83 x 10<sup>-1</sup></td><td>3.51 x 10<sup>-1</sup></td></tr><tr><td rowspan="3">Boissonneau vs. Schillings</td><td>H<sub>2</sub></td><td>1.21 x 10<sup>-1</sup></td><td>4.76 x 10<sup>-2</sup></td></tr><tr><td>O<sub>2</sub></td><td>2.99 x 10<sup>-1</sup></td><td>2.27 x 10<sup>-1</sup></td></tr><tr><td>Bulk</td><td>1.69 x 10<sup>-1</sup></td><td>1.31 x 10<sup>-1</sup></td></tr><tr><td rowspan="3">Boissonneau vs. Current Study</td><td>H<sub>2</sub></td><td>1.04 x 10<sup>-1</sup></td><td>7.95 x 10<sup>-2</sup></td></tr><tr><td>O<sub>2</sub></td><td>2.35 x 10<sup>-1</sup></td><td>6.37 x 10<sup>-2</sup></td></tr><tr><td>Bulk</td><td>2.50 x 10<sup>-1</sup></td><td>1.99 x 10<sup>-1</sup></td></tr><tr><td rowspan="12">2000 A m<sup>-2</sup></td><td rowspan="3">Boissonneau vs. Abdelouahed</td><td>H<sub>2</sub></td><td>1.43 x 10<sup>-1</sup></td><td>7.09 x 10<sup>-2</sup></td></tr><tr><td>O<sub>2</sub></td><td>2.69 x 10<sup>-1</sup></td><td>3.58 x 10<sup>-1</sup></td></tr><tr><td>Bulk</td><td>3.51 x 10<sup>-1</sup></td><td>3.70 x 10<sup>-1</sup></td></tr><tr><td rowspan="3">Boissonneau vs. Schillings</td><td>H<sub>2</sub></td><td>4.76 x 10<sup>-2</sup></td><td>4.31 x 10<sup>-2</sup></td></tr><tr><td>O<sub>2</sub></td><td>2.27 x 10<sup>-1</sup></td><td>2.90 x 10<sup>-1</sup></td></tr><tr><td>Bulk</td><td>1.31 x 10<sup>-1</sup></td><td>6.27 x 10<sup>-2</sup></td></tr><tr><td rowspan="3">Boissonneau vs. Current Study</td><td>H<sub>2</sub></td><td>7.95 x 10<sup>-2</sup></td><td>7.44 x 10<sup>-2</sup></td></tr><tr><td>O<sub>2</sub></td><td>6.37 x 10<sup>-2</sup></td><td>1.71 x 10<sup>-1</sup></td></tr><tr><td>Bulk</td><td>1.99 x 10<sup>-1</sup></td><td>2.57 x 10<sup>-1</sup></td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Void fraction evolution depending on the electrolyte width (a) and electrode height (b) for the three current density 500, 1000 and 2000 A m⁻².</figcaption></figure>
Table 5 summarizes the results of the study. It must be noticed that the diphasic boundary layer of hydrogen is thinner than the oxygen one. These results are compared with those of Schillings et al. [14]. In Schillings' study, the authors performed a dimensional study and found that the diphasic boundary layer depends on one dimensionless group, which is presented in Equations (18) and (19) (called Rayleigh-like number). The current results correspond well to the Schillings prediction. Thus, when the radius increases, the diphasic boundary layer increases, and when the equivalent injection gas velocity increases, the boundary layer decreases.

$$ \log\left(\frac{\delta}{h}\right) = -0.25 \log(Ra) + Cst \quad (18) $$

$$ Ra = \frac{\nu U_g h^5}{g r^6 H_{elec}} \quad (19) $$

The conductivity evolution depending on electrolyte gas content is well described by the Bruggeman correlation (Equation (20)). This correlation has been compared to experimental results in Hine et al. [9].

$$ \frac{\sigma(\varepsilon)}{\sigma_0} = (1 - \varepsilon)^{1.5} \quad (20) $$

$\sigma_0$ is the gas-free electrolyte conductivity in S m⁻¹.

Thus, the ohmic resistance presented in Equation (1) becomes Equation (21) [37]:

$$ R(\varepsilon) = \frac{\delta_{H2}}{\sigma_0 (1 - \varepsilon)^{1.5}} + \frac{2h - \delta_{O2} - \delta_{H2}}{\sigma_0} + \frac{\delta_{O2}}{\sigma_0 (1 - \varepsilon)^{1.5}} \quad (21) $$

In order to compare biphasic system with a gas-free system, the previous resistance can be divided by the hypothetical gas-free resistance.

$$ \frac{R(\varepsilon)}{R} = \frac{\frac{\delta_{H2}}{\sigma_0 (1-\varepsilon)^{1.5}} + \frac{2h-\delta_{O2}-\delta_{H2}}{\sigma_0} + \frac{\delta_{O2}}{\sigma_0 (1-\varepsilon)^{1.5}}}{\frac{2h}{\sigma_0}} \quad (22) $$

Table 5 shows that with an increasing current density and bubble radius, the resistance increases. This statement shows that the bubble management is an important issue.

<table><caption>**Table 5.** Maximum liquid velocity at cathode and anode sides, maximum oxygen and hydrogen void fraction and hydrogen and oxygen diphasic boundary layer for the three cases.</caption><thead><tr><th>j (A m<sup>-2</sup>)</th><th>V<sub>max liq cath</sub> (m s<sup>-1</sup>)</th><th>V<sub>max liq an</sub> (m s<sup>-1</sup>)</th><th>&epsilon;<sub>max cath</sub></th><th>&epsilon;<sub>max an</sub></th><th>&delta;<sub>cath</sub> (&mu;m)</th><th>&delta;<sub>an</sub> (&mu;m)</th><th>R(&epsilon;)/R</th></tr></thead><tbody><tr><td>500 A m<sup>-2</sup></td><td>7.8 &times; 10<sup>-2</sup></td><td>4.8 &times; 10<sup>-2</sup></td><td>0.23</td><td>0.13</td><td>515</td><td>600</td><td>1.032</td></tr><tr><td>1000 A m<sup>-2</sup></td><td>9 &times; 10<sup>-2</sup></td><td>6.8 &times; 10<sup>-2</sup></td><td>0.28</td><td>0.16</td><td>566</td><td>677</td><td>1.043</td></tr><tr><td>2000 A m<sup>-2</sup></td><td>1.25 &times; 10<sup>-1</sup></td><td>9.1 &times; 10<sup>-2</sup></td><td>0.33</td><td>0.19</td><td>690</td><td>800</td><td>1.057</td></tr></tbody></table>

## 4. Conclusions

In this study, a two-fluid multi-physics model with a new bubble transfer description has been proposed. This new description allows a good agreement with the results of Boissonneau et al.'s [15] experimental velocity profiles. It has been found out a bubble dispersion force allows a good agreement with experimental data and a better numerical convergence than the one obtained without this additional force. That observation is correlated with the studies of other studies such as Picardi et al. [19]. From the numerical point of view, the grid resolution remains a problem because even if the dispersion bubbles force improves the results, there is still a maximum grid resolution. This maximum grid resolution could lead to some inaccuracy if the fluid viscosity is too small or if the electrolyte width is too thin.
Further calculations must be performed to suppress these limitations by using a 3D geometry [19]. From the physical point of view, it can be concluded that the bubble radius and current density are two important parameters that influence the hydrodynamic. Further, calculations better properly characterize the output parameters sensitivities (velocity, void fraction etc.) to the current density and bubble diameter.

**Author Contributions:** Conceptualization, D.L.B. and P.M.; Methodology, D.L.B. and P.M.; Software, D.L.B. and P.M.; Validation, D.L.B. and P.M.; Investigation, D.L.B. and P.M.; Data curation, D.L.B. and P.M.; Writing—original draft preparation, D.L.B. and P.M.; writing—review and editing, P.M., M.B., M.K., M.S., F.G., and R.I.; Visualization, P.M. and M.B.; Supervision, P.M. and M.B.; Project administration, P.M.; Funding acquisition, P.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the national association ADEME and Bretagne Region (France).

**Acknowledgments:** We would like to deeply thank the national association ADEME for funding our research and the Bretagne region to support our research.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Castagnola, L.; Lomonaco, G.; Marotta, R. Nuclear Systems for Hydrogen Production: State of Art and Perspectives in Transport Sector. *Glob. J. Energy Technol. Res. Updates* **2014**, *1*, 4–18.

2. Schäfer, H.; Chatenet, M. Steel: The Resurrection of a Forgotten Water-Splitting Catalyst. *ACS Energy Lett.* **2018**, *3*, 574–591. [CrossRef]

3. Lavorante, M.J.; Franco, J.I. Performance of stainless steel 316L electrodes with modified surface to be use in alkaline water electrolyzers. *Int. J. Hydrogen Energy* **2016**, *41*, 9731–9737. [CrossRef]

4. Bocca, C.; Cerisola, G.; Barbucci, A.; Magnone, E. Oxygen evolution on Co2O3 and Li-doped Co2O3 coated electrodes in an alkaline solution. *Int. J. Hydrogen Energy* **1999**, *24*, 699–707. [CrossRef]

5. Boggs, B.K.; King, R.L.; Botte, G.G. Urea electrolysis: Direct hydrogen production from urine. *Chem. Commun.* **2009**, *32*, 4859–4861. [CrossRef] [PubMed]

6. Huck, M.; Ring, L.; Küpper, K.; Klare, J.P.; Daum, D.; Schäfer, H. Water Splitting Mediated by an Electrocatanytically Driven Cyclic Process Involving Iron Oxide Species; Huck, Marten, Lisa Ring, Karsten Küpper, Johann Klare, Diemo Daum, and Helmut Schäfer. *J. Mater. Chem. A* **2020**, *8*, 9896–9910. [CrossRef]

7. Babu, R.; Das, M.K. Experimental studies of natural convective mass transfer in a water-splitting system. *Int. J. Hydrogen Energy* **2019**, *44*, 14467–14480. [CrossRef]

8. Tobias, C.W. Effect of Gas Evolution on Current Distribution and Ohmic Resistance in Electrolyzers. *J. Electrochem. Soc.* **1959**, *106*, 6. [CrossRef]

9. Hine, F. Bubble Effects on the Solution IR Drop in a Vertical Electrolyzer Under Free and Forced Convection. *J. Electrochem. Soc.* **1980**, *127*, 292. [CrossRef]

10. Mandin, P.; Hamburger, J.; Bessou, S.; Picard, G. Modelling and calculation of the current density distribution evolution at vertical gas-evolving electrodes. *Electrochim. Acta* **2005**, *51*, 1140–1156. [CrossRef]

11. Hreiz, R.; Abdelouahed, L.; Fünfschilling, D.; Lapicque, F. Electrogenerated bubbles induced convection in narrow vertical cells: PIV measurements and Euler–Lagrange CFD simulation. *Chem. Eng. Sci.* **2015**, *134*, 138–152. [CrossRef]

12. Dahlkild, A.A. Modelling the two-phase flow and current distribution along a vertical gas-evolving electrode. *J. Fluid Mech.* **2001**, *428*, 249–272. [CrossRef]

13. Wedin, R.; Dahlkild, A.A. On the Transport of Small Bubbles under Developing Channel Flow in a Buoyant Gas-Evolving Electrochemical Cell. *Ind. Eng. Chem. Res.* **2001**, *40*, 5228–5233. [CrossRef]

14. Schillings, J.; Doche, O.; Deseure, J. Modeling of electrochemically generated bubbly flow under buoyancy-driven and forced convection. *Int. J. Heat Mass Transf.* **2015**, *85*, 292–299. [CrossRef]

15. Boissonneau, P.; Byrne, P. An experimental investigation of bubble-induced free convection in a small electrochemical cell. *J. Appl. Electrochem.* **2000**, *30*, 767–775. [CrossRef]

16. Aldas, K.; Pehlivanoglu, N.; Mat, M.D. Numerical and experimental investigation of two-phase flow in an electrochemical cell. *Int. J. Hydrogen Energy* **2008**, *33*, 3668–3675. [CrossRef]
17. Mat, M. A two-phase flow model for hydrogen evolution in an electrochemical cell. *Int. J. Hydrogen Energy* **2004**, *29*, 1015–1023. [CrossRef]

18. Riegel, H.; Mitrovic, J.; Stephan, K. Role of mass transfer on hydrogen evolution in aqueous media. *J. Appl. Phys.* **1998**, *28*, 10–17.

19. Picardi, R.; Zhao, L.; Battaglia, F. On the Ideal Grid Resolution for Two-Dimensional Eulerian Modeling of Gas–Liquid Flows. *J. Fluids Eng.* **2016**, *138*, 114503. [CrossRef]

20. Law, D.; Battaglia, F.; Heindel, T.J. Model validation for low and high superficial gas velocity bubble column flows. *Chem. Eng. Sci.* **2008**, *63*, 4605–4616. [CrossRef]

21. Panicker, N.; Passalacqua, A.; Fox, R.O. On the hyperbolicity of the two-fluid model for gas—Liquid bubbly flows. *Appl. Math. Model.* **2018**, *57*, 432–447. [CrossRef]

22. Vaidheeswaran, A.; de Bertodano, M.L. Stability and convergence of computational eulerian two-fluid model for a bubble plume. *Chem. Eng. Sci.* **2017**, *160*, 210–226. [CrossRef]

23. Marfaing, O.; Guingo, M.; Laviéville, J.; Bois, G.; Méchitoua, N.; Mérigoux, N.; Mimouni, S. An analytical relation for the void fraction distribution in a fully developed bubbly flow in a vertical pipe. *Chem. Eng. Sci.* **2016**, *152*, 579–585. [CrossRef]

24. Davidson, M.R. Numerical calculations of two-phase flow in a liquid bath with bottom gas injection: The central plume. *Appl. Math. Model.* **1990**, *14*, 67–76. [CrossRef]

25. Lee, J.W.; Sohn, D.K.; Ko, H.S. Study on bubble visualization of gas-evolving electrolysis in forced convective electrolyte. *Exp. Fluids* **2019**, *60*, 156. [CrossRef]

26. Abdelouahed, L.; Hreiz, R.; Poncin, S.; Valentin, G.; Lapicque, F. Hydrodynamics of gas bubbles in the gap of lantern blade electrodes without forced flow of electrolyte: Experiments and CFD modelling. *Chem. Eng. Sci.* **2014**, *111*, 255–265. [CrossRef]

27. Ham, J.; Homsy, G. Hindered settling and hydrodynamic dispersion in quiescent sedimenting suspensions. *Int. J. Multiph. Flow* **1988**, *14*, 533–546. [CrossRef]

28. Nagai, N.; Takeuchi, M.; Nakao, M. Influences of Bubbles between Electrodes onto Efficiency of Alkaline Water Electrolysis. *Pac. Cent. Therm. Fluids Eng.* **2003**, *10*, F4008.

29. Abdelouahed, L. Gestion Optimale du Gaz Électrogénéré Dans un Réacteur D’électroréduction de Minerai de fer. Ph.D. Thesis, Université de Lorraine, Nancy, France, 2013.

30. Schillings, J. Etude Numérique et Expérimentale D’écoulements Diphasiques: Application aux Écoulements à Bulles Générées par Voie Électrochimique. Ph.D. Thesis, Université Grenoble Alpes, Grenoble, France, 2018.

31. Le Bideau, D.; Mandin, P.; Sellier, M.; Benbouzid, M.; Myeongsub, K. Review of necessary thermophysical properties and their sensitivities with local temperature and electrolyte mass fraction for alkaline water electrolysis multiphysics modelling. *Int. J. Hydrogen Energy* **2019**, *44*, 4553–4569. [CrossRef]

32. Alexiadis, A.; Dudukovic, M.P.; Ramachandran, P.; Cornell, A.; Wanngård, J.; Bokkers, A. Liquid–gas flow patterns in a narrow electrochemical channel. *Chem. Eng. Sci.* **2011**, *66*, 2252–2260. [CrossRef]

33. Alexiadis, A.; Dudukovic, M.P.; Ramachandran, P.; Cornell, A.; Wanngård, J.; Bokkers, A. On the electrode boundary conditions in the simulation of two phase flow in electrochemical cells. *Int. J. Hydrogen Energy* **2011**, *36*, 8557–8559. [CrossRef]

34. Charton, S.; Rivalier, P.; Ode, D.; Morandini, J.; Caire, J.P. *Hydrogen Production by the Westinghouse Cycle: Modelling and Optimization of the Two-Phase Electrolysis Cell*; WIT Press: Bologna, Italy, 2009; pp. 11–22.

35. Grigoriev, S.A.; Fateev, V.N.; Bessarabov, D.G.; Millet, P. Current status, research trends, and challenges in water electrolysis science and technology. *Int. J. Hydrogen Energy* **2020**. [CrossRef]

36. Isono, T. Density, viscosity, and electrolytic conductivity of concentrated aqueous electrolyte solutions at several temperatures. Alkaline-earth chlorides, lanthanum chloride, sodium chloride, sodium nitrate, sodium bromide, potassium nitrate, potassium bromide, and cadmium nitrate. *J. Chem. Eng. Data* **1984**, *29*, 45–52. [CrossRef]

37. Mandin, P.; Derhoumi, Z.; Roustan, H.; Rolf, W. Bubble Over-Potential during Two-Phase Alkaline Water Electrolysis. *Electrochim. Acta* **2014**, *128*, 248–258. [CrossRef]

<figure><img src="image_9.png" /></figure>