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

#### **Damien Le Bideau 1, Philippe Mandin 1,\*, Mohamed Benbouzid 2, Myeongsub Kim 3, Mathieu Sellier 4, Fabrizio Ganci <sup>5</sup> and Rosalinda Inguanta <sup>5</sup>**


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 (H2)). There are several processes that do not emit greenhouse gases and produce H2, such as thermochemical hydrogen process and electrolytic

H2 process. Thermochemical processes for hydrogen production use heat to decompose water into hydrogen and oxygen [1], whereas electrolytic H2 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).

$$\mathcal{U}\_{\text{cell}}(j) = E\_{\text{rev}}(P, T) + \eta(j)\_{\text{cath act}} + \eta(j)\_{\text{an act}} + R(T, \,\, Y\_{\text{KOH}}) \, j + \eta(j)\_{\text{cath}|\text{an canc}} \tag{1}$$

*Erev* is the reversible voltage in *V*, η*act* is the activation overvoltage in *V*, *R* the sum of the electrodes, membrane and electrolyte resistance in Ω cm<sup>2</sup> and η*conc* the concentration overpotential in *V*, *j* the current density in A m<sup>−</sup>2.

The *Erev* 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 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.

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/ <sup>√</sup>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 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.

#### *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.

$$
\stackrel{\rightarrow}{F}\_{BD} = -\underbrace{\varepsilon\_{\mathcal{S}} \rho \stackrel{K\_{\mathcal{S}}}{d\_{\mathcal{b}}} |\mathcal{U}\_{\mathcal{l}}| \stackrel{\rightarrow}{\nabla} \varepsilon\_{\mathcal{S}}}\_{\textit{Bubble dispersion factor}} \tag{2}
$$

ε is the gas or liquid fraction, the subscript *k* can be either *g* (O2, H2) or *liq* (for liquid), ρ is the density in kg m<sup>−</sup>3, *db* is the bubble diameter in m, *Ur* is the gas phase velocity minus the liquid velocity in m s<sup>−</sup>1.

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:


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 \varepsilon\_k \rho\_k}{\partial t} + \overrightarrow{\nabla} \cdot \left( \varepsilon\_k \rho\_k \overrightarrow{V}\_k \right) = \mathcal{S}\_{\mathcal{S}} \tag{3}$$

ε is the gas or liquid fraction, the subscript *k* can be either *g* (O2, H2) or *liq*, ρ is the density in kg m<sup>−</sup>3, *V* the velocity in m s<sup>−</sup>1, and *Sg* is the term source in kg m−<sup>3</sup> s<sup>−</sup>1.

$$\frac{\partial}{\partial t} \left( \varepsilon\_k \rho\_k \stackrel{\cdot}{V}\_k \right) + \overrightarrow{\nabla} \cdot \left( \varepsilon\_k \rho\_k \stackrel{\cdot}{V}\_k \stackrel{\cdot}{V}\_k \right) = -\varepsilon\_k \overrightarrow{\nabla} p + \overrightarrow{\nabla} \cdot (\varepsilon\_k \tau) + \varepsilon\_k \rho \stackrel{\cdot}{S} + \overrightarrow{F}\_k \tag{4}$$

*<sup>p</sup>* is the pressure in Pa, <sup>τ</sup> is the stress tensor in Pa, *<sup>g</sup>* the gravitational acceleration in m s<sup>−</sup>2, and <sup>→</sup> *Fk* is the exchange term in N m<sup>−</sup>3.

The stress tensor is written as follows

$$
\pi = \mu\_k \left[ \left( \overrightarrow{\nabla V}\_k + \overrightarrow{\nabla V}\_k^T \right) - \frac{2}{3} \overrightarrow{\nabla} \cdot \overrightarrow{V}\_k I \right] \tag{5}
$$

with μ*<sup>k</sup>* the viscosity of the phase *k* in Pa s and *I* the unit tensor.

$$
\overrightarrow{F}\_k = \overrightarrow{F}\_D + \overrightarrow{F}\_L + \overrightarrow{F}\_{BD} \tag{6}
$$

$$\overrightarrow{F}\_{k} = \underbrace{-\frac{3}{4}\varepsilon\_{\mathcal{S}}\rho\frac{\mathbf{C}\_{D}}{d\_{b}}}\_{\text{Drag force}}\underbrace{\mathrm{l}\mathrm{L}\_{l}\mathrm{l}\mathrm{l}\mathrm{L}\_{r}}\_{\text{Lift force}} - \underbrace{\varepsilon\_{\mathcal{S}}\rho\mathrm{C}\_{L}\mathrm{l}\mathrm{l}\mathrm{l}\mathrm{l}\mathrm{r}\mathrm{f}\mathrm{r}\mathrm{s}\mathrm{f}\mathrm{d}\mathrm{l}\mathrm{s}}\_{\text{Lift force}} - \underbrace{\varepsilon\_{\mathcal{S}}\rho\frac{\mathrm{K}\_{\mathcal{S}}}{d\_{b}}\mathrm{l}\mathrm{L}\_{r}\mathrm{l}\overrightarrow{\nabla}\mathrm{c}\_{\mathcal{S}}}\_{\text{Bubble dispersion force}}\tag{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 Na2SO4 concentrated at 50 g L<sup>−</sup>1. Therefore, the reaction occurring at the anode and the cathode are Equations (8) and (9), respectively.

$$2H\_2O\_{(liq)} \to 2H^+\_{(aq)} + \frac{1}{2} \text{ O}\_{2(g)} + 2e^- \tag{8}$$

$$2H\_2O\_{(liq)} + 2e^- \to 2OH^-\_{(aq)} + H\_{2(g)}\tag{9}$$

If the reader wants to simulate an alkaline electrolysis, he must change the thermophysical inputs from Na2SO4 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 = \mathbf{n} \frac{\mathbf{j} \cdot \mathbf{S}}{F} \mathbf{M}\_{\mathbf{S}} \tag{10}$$

*qm* is the mass flow of produced gas in kg s−1, *Mg* is the molar mass of the gas kg mol−1, *S* is the electrode surface in m2, *F* is the Faraday constant 96,500 C mol<sup>−</sup>1, 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\_{\mathcal{S}} = \mathbf{n} \frac{j}{F \times \Delta x} M\_{\mathcal{S}} \tag{11}$$

Δ*x* is the width of the first cell next to the electrode in m. For the bottom and top boundary condition, a pressure inlet (*PTot* = 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−1. The Table 1 summarizes the boundary conditions for the present problem.

**Table 1.** The boundary conditions to solve the problem.


#### *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 \wp = \sum\_{i}^{N} a\_i \wp + b$$

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<sup>−</sup>2. The most interesting datapoints are those obtained at 2000 A m<sup>−</sup>2. Indeed, industrial alkaline electrolyzers generally use a current density between 2000 and 5000 A m−<sup>2</sup> [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 2.** The input values for the problem.



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).

$$\text{Re}\_G = \frac{\rho\_l \, V\_G \, \, H\_{\text{clec}}}{\mu\_L} \tag{13}$$

$$Fr\_G = \frac{\text{g } H\_{\text{elac}}}{V\_G^2} \tag{14}$$

$$r^\* = \frac{d}{2\ H} \tag{15}$$

$$h^\* = \frac{h}{H} \tag{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 *FrG* numbers higher than 105.

$$\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 \tag{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 3.** Result of the liquid velocity in different models at *j* = 500 A m−<sup>2</sup> (**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.

**Figure 4.** Results of the liquid velocity in different models at *j* = 1000 A m−<sup>2</sup> (**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.

**Figure 5.** Result of the liquid velocity in different models at *j* = 2000 A m−<sup>2</sup> (**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.

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 different of the different model with the Boissonneau et al. [15] at three current density (500, 1000 and 2000 A m<sup>−</sup>2) 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 O2 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 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. H2 is the zone between [0, 500] μm and O2 is zone between [2.5, 3] mm.

**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<sup>−</sup>2.

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 the 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) + \text{Cst}\tag{18}
$$

$$Ra = \frac{\nu}{\mathcal{g}} \frac{\mathcal{U}\_{\mathcal{S}} h^5}{r^6 \, H\_{\text{clcc}}} \tag{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} = \left(1 - \varepsilon\right)^{1.5} \tag{20}$$

σ<sup>0</sup> is the gas-free electrolyte conductivity in S m<sup>−</sup>1.

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

$$R(\varepsilon) = \frac{\delta\_{H2}}{\sigma\_0 \left(1 - \varepsilon\right)^{1.5}} + \frac{2h - \delta\_{O2} - \delta\_{H2}}{\sigma\_0} + \frac{\delta\_{O2}}{\sigma\_0 \left(1 - \varepsilon\right)^{1.5}} \tag{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}}{\varepsilon\_0 \left(1 - \varepsilon\right)^{1.5}} + \frac{2\hbar - \delta\_{\rm CO} - \delta \mu\_2}{\varepsilon\_0} + \frac{\delta \varepsilon\_{\rm CO}}{\varepsilon\_0 \left(1 - \varepsilon\right)^{1.5}}}{\frac{2\hbar}{\varepsilon\_0}}\tag{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 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.


#### **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**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
