# Evolutionary Design Optimization of an Alkaline Water Electrolysis Cell for Hydrogen Production

Damien Le Bideau, Olivier Chocron, Philippe Mandin, Patrice Kiener,

Mohamed Benbouzid, Mathieu Sellier, Myeongsub Kim, Fabrizio Ganci,

Rosalinda Inguanta

► To cite this version:

Damien Le Bideau, Olivier Chocron, Philippe Mandin, Patrice Kiener, Mohamed Benbouzid, et al. Evolutionary Design Optimization of an Alkaline Water Electrolysis Cell for Hydrogen Production.

Applied Sciences, 2020, 10 (23), pp.8425. 10.3390/app10238425 . hal-05021036

HAL Id: hal-05021036

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

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
<figure><img src="image_2.png" /></figure>

Article

# Evolutionary Design Optimization of an Alkaline Water Electrolysis Cell for Hydrogen Production

Damien Le Bideau¹, Olivier Chocon², Philippe Mandin ¹,* ³, Patrice Kiener ³, Mohamed Benbouzid ⁴, Mathieu Sellier ⁵, Myeongsub Kim ⁶, Fabrizio Ganci ⁷ and Rosalinda Inguanta ⁷

¹ Institut de Recherche Dupuy de Lôme (UMR CNRS 6027 IRDL), University Bretagne Sud, 56100 Lorient, France; damien.le-bideau@univ-ubs.fr

² Institut de Recherche Dupuy de Lôme (UMR CNRS 6027 IRDL), ENI Brest, 29238 Brest CEDEX 3, France; chocon@enib.fr

³ InModelia, 5 rue Malebranche, 75005 Paris, France; patrice.kiener@inmodelia.com

⁴ Institut de Recherche Dupuy de Lôme (UMR CNRS 6027 IRDL), University of Brest, 29238 Brest, France; Mohamed.Benbouzid@univ-brest.fr

⁵ Department of Mechanical Engineering, University of Canterbury, Christchurch 8140, New Zealand; mathieu.sellier@canterbury.ac.nz

⁶ Department of Ocean and Mechanical Engineering, Florida Atlantic University, Boca Raton, FL 33431, USA; kimm@fau.edu

⁷ Dipartimento di Ingegneria, University of Palermo, Viale delle Scienze, 90128 Palermo, Italy; fabrizio.ganci@unipa.it (F.G.); rosalinda.inguanta@unipa.it (R.I.)

* Correspondence: philippe.mandin@univ-ubs.fr

Received: 30 October 2020; Accepted: 24 November 2020; Published: 26 November 2020

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

check for
updates

**Featured Application:** This article provides the strategy to find the best cell design and operating condition of a diphasic electrochemical cell using an evolutionary algorithm.

**Abstract:** Hydrogen is an excellent energy source for long-term storage and free of greenhouse gases. However, its high production cost remains an obstacle to its advancement. The two main parameters contributing to the high cost include the cost of electricity and the cost of initial financial investment. It is possible to reduce the latter by the optimization of system design and operation conditions, allowing the reduction of the cell voltage. Because the CAPEX (initial cost divided by total hydrogen production of the electrolyzer) decreases according to current density but the OPEX (operating cost depending on the cell voltage) increases depending on the current density, there exists an optimal current density. In this paper, a genetic algorithm has been developed to find the optimal evolution parameters and to determine an optimum electrolyzer design. The optimal current density has been increased by 10% and the hydrogen cost has been decreased by 1%.

**Keywords:** alkaline water electrolysis; hydrogen cost; genetic algorithm; optimization

## 1. Introduction

Mitigation of the climate change consequences is the main global challenge in the 21 century. The different IPCC (Intergovernmental Panel on Climate Change) reports stated that the climate change is mainly due to greenhouse gas emissions to the atmosphere produced by human activities. The IPCC recommends in their scenarios the decrease of energy consumption, the increase of renewable production, and sequestration of harmful CO₂ produced by the industrial activities including the production of cement and steel. The renewable energies generated from various systems such as wind turbines and photovoltaic panels do not produce greenhouse gases but suffer from time intermittency
and geographic limitations. In order to resolve this challenge, energy produced at off-peak hours must be stored. Batteries can store this energy, but their energy density is low, and they discharge the stored energy over time. Thus, if the energy needs to be stored more than a month, another storage mean must be provided. Hydrogen storage can store (renewable) energy with a high density for a long period and its use (by combustion and by reduction) only produces water. Most hydrogen is currently produced via a process that needs fossil fuels and which produces greenhouse gases. However, although there exist several other processes that can generate hydrogen without generating pollutants, water electrolysis is one of these processes and is very promising. Nevertheless, the high operation cost hampers its development. One way to decrease the cost is to choose the cheapest water electrolysis technology: Alkaline water electrolysis. There are two other technologies: PEM electrolysis and high temperature electrolysis. They both are theoretically more efficient than alkaline water electrolysis but they need noble and expensive material (PEM) or suffer from low life duration (high temperature electrolysis). Although alkaline electrolysis is the cheapest technology among all the electrolysis technology, it is still needed to decrease the hydrogen cost. Like most electrochemical processes, the cost of the hydrogen production is divided into two parts: CAPEX and OPEX. CAPEX represents the capital cost divided by the total hydrogen production. The OPEX is the production process cost. Thus, the higher the cell voltage increase, the higher the OPEX is.

The cell voltage is a sum of overvoltage as shown in Equation (1).

$$U_{cell}(j) = E_{rev}(P, T) + R(T, Y_{KOH}, \varepsilon) j + \eta_{act an}\left(\frac{j}{1 - \theta}\right) + \eta_{act cath}\left(\frac{j}{1 - \theta}\right) + R_{memb} j \quad (1)$$

$E_{rev}$ is the reversible voltage in V, $\eta_{act}$ is the activation overvoltage in V, $R$ and $R_{memb}$ respectively the electrolyte resistance and membrane resistance in $\Omega$ cm$^2$ and, $j$ the current density in A m$^{-2}$.

The $E_{rev}$ voltage represents the minimum thermodynamic cell voltage. This value is around 1.23 V at atmospheric pressure at 25 °C. Increasing temperature and decreasing the pressure decrease this minimum cell voltage. The majority of previous studies have focused on finding a cheap, robust and electroactive electrode material or electrocatalyst to decrease the activation overpotential, which are the second and third terms in Equation (1). There is also another research trend aiming at increasing efficiency, which is to reduce an ohmic resistance and surface coverage through electrolyzer design by means of simulations or experiments. Previous studies modeled the void fraction and velocity distribution [1–5] or the cell voltage [6]. There exist only few studies about optimization of this process. Villagra et al. [7] showed that there exists an optimum current density for a given design of PEM electrolyzer. Bensmann et al. [8] demonstrated that there exists an optimum of pressure for a given PEM electrolyzer.

The goal of this study is to find the optimal design and the current density for an alkaline electrolysis cell. A meta-model has been developed to simplify the physics of the void fraction and a genetic algorithm is used to find the optimal design. The optimization problem uses a complex objective function involving multiple, coupled and non-linear variables. The optimization of such a function is difficult because there are no finite time resolution algorithms. The problem is defined to be NP-complete and cannot be solved with classical methods such as the gradient descent [9]. For this reason, artificial evolution techniques (evolutionary algorithms) have been explored in the design of the electrolyzer. This paper begins with the theory, physics and economy behind the alkaline water electrolysis cell. Then, the set-up, adjustment of parameters according to the problem features and validation of the genetic algorithm are highlighted. Finally, the algorithm is used to find the optimal design of the alkaline water electrolysis.
## 2. Physics and Economy of the Alkaline Water Electrolysis Cell

### 2.1. Theory

An alkaline water electrolysis cell produces hydrogen and oxygen using electric energy via an electrochemical process. The Figure 1 presents a classical water electrolysis cell. It consists of a positive pole called the anode (where the oxidation of hydroxyl ions occurs Equation (2)), a negative pole called the cathode (where the reduction of water occurs Equation (3)) and a membrane (brown part) that prevents the mixing of hydrogen and oxygen.

$$2OH^{-}_{(aq)} \rightarrow H_{2}O_{(liq)} + \frac{1}{2}O_{2(g)} + 2e^{-} \quad (2)$$

$$2H_{2}O_{(liq)} + 2e^{-} \rightarrow 2OH^{-}_{(aq)} + H_{2}(g) \quad (3)$$

<figure><img src="image_4.png" /><figcaption>Figure 1. Design of the electrochemical cell.</figcaption></figure>

The reactions of Equations (2) and (3) occur at alkaline pH (>7). Thus, KOH or NaOH salts must be added. The classical operating conditions of alkaline electrolyzer are: Temperature around 80 °C, atmospheric pressure and 0.30 KOH. In the next paragraph, we will explain why those choices have been made.

As explained in the introduction section, the minimum voltage that must be applied to a water electrolysis cell is called reversible voltage. Hammoudi et al. [10] developed a model (Appendix A) that simulates the temperature, pressure, and mass fraction sensitivity of the reversible voltage. The reversible voltage decreases with temperature while mass fraction and the pressure increase it.

The ohmic electrolyte resistance R is depending on the electrolyte conductivity σ, the diphasic boundary thickness δ and the void fraction ε, see Equation (4) [11].

$$R(\varepsilon, \delta, T, Y) = \frac{\delta}{\sigma(\varepsilon, T, Y)} + \frac{h - \delta}{\sigma(T, Y)} \quad (4)$$
With h the thickness of the electrolyte (catholyte or anolyte) in m, σ the electrical conductivity in S m⁻¹.

Hine et al. [12] showed the electrical conductivity dependency on void fraction can be simulated using the Bruggeman’s correlation Equation (5).

$$ \frac{\sigma(\varepsilon, T, Y)}{\sigma(T, Y)} = (1 - \varepsilon)^{\frac{3}{2}} \qquad (5) $$

The Equation (6) can be simplified by assuming that the gas is dispersed throughout the electrolyte thickness.

$$ R(\varepsilon, T, Y) = \frac{h}{\sigma(\varepsilon, T, Y)} \qquad (6) $$

Thus, ohmic resistance increases with void fraction. The quantity of gas injected in the electrolyte is calculated via Faraday’s law Equation (7). According to the Equation (8), the pressure decreases the volumetric hydrogen flow and thus the void fraction. Moreover, Vogt et al. [13] and Jannsen et al. [14] showed that the pressure decreases the detachment bubble radius, thus decreasing the surface coverage. Gilliam et al. [15] and See et al. [16] measured the electrical conductivity depending on the KOH mass fraction and temperature. They showed that there exists an optimum KOH mass fraction between 0.2 and 0.4 that depends on the temperature. The different correlations for electrical conductivity, viscosity, density of KOH and NaOH have been summarized by Le Bideau et al. [17]. The temperature increases the efficiency of the process but a high temperature and an extremely alkaline (pH > 14) condition led to material deterioration and electrolyte boiling. Thus, the temperature should not exceed 100 °C.

$$ Q_{vH_2} = \frac{j S v_{mol}}{2 F} \qquad (7) $$

$$ v_{mol} = \frac{R_G T}{p} \qquad (8) $$

$Q_{vH_2}$ the volumetric hydrogen flow m³ s⁻¹, S the electrode surface in m², $v_{mol}$ the molar volume m³ mol⁻¹, F = 96,500 C mol⁻¹ the Faraday’s constant in C mol⁻¹, $R_G$ = 8.314 J kg K⁻¹ the ideal gas constant in J mol⁻¹ K⁻¹, T the temperature in K, p the pressure in Pa.

Therefore, this analysis explained the choice of mass fraction and temperature but not pressure choice. Actually, high pressure increases the manufacturing cost. However, according to Grigoriev et al. [18], the aim of pressure level for the new alkaline electrolyzer is 30 bar.

According to Equation (6), it would be tempting to decrease as much as possible the electrolyte thickness. However, Nagai et al. [19–21] observed in multiple studies that there exists an optimum of electrolyte thickness. Indeed, in an extremely narrow electrolyte the bubbles cannot leave fast enough to prevent bubble accumulation.

The activation overpotential can be modeled with different mathematical formalism depending on the current density. For the industrial current density, the Tafel law is used, see Equation (9).

$$ \eta_{act}(j) = a + b \log\left(\frac{j}{1 - \theta}\right) \qquad (9) $$

With a and b the Tafel parameter in respectively V and V dec⁻¹ and θ the surface coverage ratio.

The parameter a and b depend on the electrode material and electrocatalyst. The surface coverage depends on a lot of parameters (temperature, pressure, mass fraction, current density, presence or absence of electrolyte flow). However, there is no numerical model that can simulate the surface coverage ratio depending on all these parameters. Vogt et al. [22] developed a model that roughly estimates the surface coverage depending on current density Equation (10).

$$ \theta = 0.023 j^{0.3} \qquad (10) $$
Using the data from Jannsen et al. [14], the correlation that gives the bubble radius according to the pressure is given by Equation (11).

$$r(p) = 0.876 r_0^{-0.471 p} \qquad (11)$$

$r_0$ the mean bubble radius at atmospheric pressure in m.

The pressure, current density and temperature sensitivities of cell voltage have been given. One of the key parameters is the void fraction. This parameter depends on several parameters. The sensitivities will be explained in the next section.

## 2.2. Physics of the Void Fraction

### 2.2.1. Mathematical and Numerical Tools

#### Finite Volume Model (CFD)

The void fraction can be calculated using numerical simulations. A model has been developed and validated by Le Bideau et al. [1] in a previous study. The following hypothesis have been chosen:

* The flow is Newtonian, viscous and incompressible;

* the flow is considered isothermal;

* ions distributions are neglected;

* the flow is considered laminar;

* bubble diameter is constant for a given operating condition; and

* The current density is taken as constant.

For more explanations about the hypotheses, see Le Bideau et al. [1]. The model is a two-fluid Eulerian model Equations (12)–(16).

$$\frac{\partial \varepsilon_k \rho_k}{\partial t} + \vec{\nabla} \cdot (\varepsilon_k \rho_k \vec{V}_k) = S_k \qquad (12)$$

$\varepsilon$ 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_k$ is the term source in kg m⁻³ s⁻¹

$$\frac{\partial}{\partial t} (\varepsilon_k \rho_k \vec{V}_k) + \vec{\nabla} \cdot (\varepsilon_k \rho_k \vec{V}_k \vec{V}_k) = - \varepsilon_k \vec{\nabla} p + \vec{\nabla} \cdot (\varepsilon_k \bar{\tau}) + \varepsilon_k \rho \vec{g} + \vec{F}_k \qquad (13)$$

p is the pressure in Pa, $\bar{\tau}$ is the stress tensor in Pa, $\vec{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:

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

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} \qquad (15)$$

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

Equation (15) is the exchange term. It is the sum of the usual drag force and lift force. The study of Le Bideau et al. [1] showed that it was mandatory to take into account the bubble diffusion force
$\vec{F}_{BD}$) to reproduce the experimental data. The coefficient $K_g$ has been used to fit the experimental data of Boissonneau et al. [23].

Using the Vaschy-Buckingham theorem, four dimensionless numbers have been determined Equations (17)–(20).

$$ Re_{V_G} = \frac{\rho_l V_G H_{elec}}{\mu_L} \qquad (17) $$

$$ Fr_{V_G} = \frac{g H_{elec}}{V_G^2} \qquad (18) $$

$$ r^* = \frac{r}{H} \qquad (19) $$

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

## Artificial Neural Network (ANN)

These four dimensionless numbers (Equations (21)–(24)) have been used to correlate the void fraction using an artificial neural network (ANN). Appendix B shows the design of experiments used to train the ANN. The numerical search field is summarized in Table 1. This numerical search field allows the finding an optimum design for an average current density $j$ in $[10^3–10^4]$ A m$^{-2}$, an electrode height $H_{elec}$ in [0.01–0.1] m, a temperature $T$ in [293–350] K, a pressure $p$ in [1–30] bar, a mass fraction $Y$ in [0.2–0.4], an electrolyte thickness $h$ in $[1.5 \times 10^{-4}–1.5 \times 10^{-3}]$ m.

<table><caption>**Table 1.** Numerical search field for every dimensionless parameters.</caption><thead><tr><th>Dimensionless Parameter</th><th>x<sub>min</sub></th><th>x<sub>max</sub></th><th>x<sub>max</sub>+x<sub>min</sub>/2</th><th>x<sub>max</sub>-x<sub>min</sub>/2</th></tr></thead><tbody><tr><td>Re<sub>V_G</sub></td><td>1 × 10<sup>-2</sup></td><td>300</td><td>150.005</td><td>149.995</td></tr><tr><td>Fr<sub>V_G</sub></td><td>5.85 × 10<sup>4</sup></td><td>5.30 × 10<sup>10</sup></td><td>2.65 × 10<sup>10</sup></td><td>2.65 × 10<sup>10</sup></td></tr><tr><td>r<sup>*</sup></td><td>1.14 × 10<sup>-5</sup></td><td>1.50 × 10<sup>-3</sup></td><td>7.75 × 10<sup>-2</sup></td><td>7.25 × 10<sup>-2</sup></td></tr><tr><td>h<sup>*</sup></td><td>5 × 10<sup>-3</sup></td><td>1.50 × 10<sup>-1</sup></td><td>7.557 × 10<sup>-4</sup></td><td>7.442 × 10<sup>-4</sup></td></tr></tbody></table>

The dimensionless numbers have been reduced using Equation (21).

$$ \bar{X} = \frac{x - \left( \frac{x_{max} + x_{min}}{2} \right)}{\frac{x_{max} - x_{min}}{2}} \qquad (21) $$

With $\bar{X}$ the reduced dimensionless number, $x$ the dimensionless number

The weights of the ANN are obtained using the software NeurOne©. Because the evolutions of the void fraction according to the four dimensionless parameters are non-linear, the learning is done by a Levenberg-Marquardt (L-M) algorithm [24] with parameters initialized according to a normal law N(0, 0.1) over about 70 iterations. The L-M algorithm consists in reducing a Sum Square Error (SSE) defined by Equation (22).

$$ SSE\left(\left(\begin{array}{c} \overline{Re}_{V_G} \\ \overline{Fr}_{V_G} \\ \overline{h}^* \\ \overline{r}^* \end{array}\right), \vec{w}_{N,M}\right) = \frac{1}{2} \sum_{i=0}^{i=O} e_i^2 = \frac{1}{2} \sum_{i=0}^{i=O} (\varepsilon_{fluent,i} - \varepsilon_{ANN,i})^2 \qquad (22) $$

$\vec{w}_{N,M}$ the weight vector with the subscript $N$ the number of layer (from 0 to $N$) and $M$ the total number of weights per layer and $O$ the number of outputs.
The SSE is reduced by a combination of two algorithm: the steepest descent algorithm and the Gauss–Newton algorithm [24]. Thus, the weight vector is updated following the Equation (23).

$$J_k = \begin{bmatrix} \frac{\partial e}{\partial w_{0,0}} & \frac{\partial e}{\partial w_{1,0}} & \cdots & \frac{\partial e}{\partial w_{N,0}} \\ \frac{\partial e}{\partial w_{0,1}} & \frac{\partial e}{\partial w_{1,1}} & \cdots & \frac{\partial e}{\partial w_{N,1}} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial e}{\partial w_{0,M}} & \frac{\partial e}{\partial w_{1,M}} & \cdots & \frac{\partial e}{\partial w_{N,M}} \end{bmatrix} \qquad (23)$$

$$\vec{w}_{N,Mk+1} = \vec{w}_{N,Mk} - (J_k^T J_k + \lambda I)^{-1} J_k (\vec{\epsilon}_{fluent,k} - \vec{\epsilon}_{ANN,k}) \quad (24)$$

With $J_k$ the Jacobian matrix, the $k$ subscript for the number of update and $\lambda$ the combination coefficient.

If $\lambda$ is close to 0 then Gauss–Newton algorithm is used whereas when $\lambda$ is very large the L–M algorithm switches to the steepest descent algorithm.

The Figure 2 gives the architecture of the artificial neural network used to estimate the void fraction, which correspond to the Equations (25)–(29).

$$\text{for } i \text{ in } [0;3] \qquad H_{1,i} = w_{0,i} \times H_{0,i} + c_{0,i} \qquad (25)$$

$$\text{for } i \text{ in } [0;2] \qquad H_{2,i} = \text{atan}\left(\left(\sum_{k=0}^{3} w_{1,(k \times 3)} + i \times H_{0,k}\right) + c_{1,i}\right) \qquad (26)$$

$$H_{3,0} = \left( \sum_{i=0}^{2} w_{2,i} H_{2,i} \right) + c_{2,0} \qquad (27)$$

$$H_{4,0} = w_{3,0} H_{3,0} + c_{3,0} \qquad (28)$$

$$H_{5,0} = \exp(w_{4,0} H_{4,0}) \qquad (29)$$

<figure><img src="image_5.png" /><figcaption>Figure 2. Architecture of the artificial neural network used to estimate the void fraction d.</figcaption></figure>

The parameters of Equations (25)–(29) are given in the Table 2.
<table><caption>**Table 2.** Parameters for Equations (25)–(29).</caption><thead><tr><th></th><th></th><th>i = 0</th><th>i = 1</th><th>i = 2</th><th>i = 3</th></tr></thead><tbody><tr><td rowspan="2">H<sub>1,i</sub></td><td>w<sub>0,i</sub></td><td>1.3965</td><td>1.3532</td><td>1.4242</td><td>1.4926</td></tr><tr><td>c<sub>0,i</sub></td><td>0.61711</td><td>0.6445</td><td>0.5920</td><td>0.1461</td></tr><tr><td rowspan="5">H<sub>2,i</sub></td><td>w<sub>1,i</sub></td><td>41.63</td><td>-0.848</td><td>-0.444</td><td>-</td></tr><tr><td>w<sub>1,(i+3)</sub></td><td>-2.6544</td><td>-0.837</td><td>114</td><td>-</td></tr><tr><td>w<sub>1,(i+6)</sub></td><td>-1.5778</td><td>-9.13</td><td>-0.015</td><td>-</td></tr><tr><td>w<sub>1,(i+9)</sub></td><td>0.2068</td><td>-0.154</td><td>-0.347</td><td>-</td></tr><tr><td>c<sub>1,i</sub></td><td>29.78</td><td>-13.57</td><td>81.04</td><td>-</td></tr><tr><td rowspan="4">H<sub>3,i</sub></td><td>w<sub>2,0</sub></td><td>-0.4</td><td>-</td><td>-</td><td>-</td></tr><tr><td>w<sub>2,1</sub></td><td>-0.959</td><td>-</td><td>-</td><td>-</td></tr><tr><td>w<sub>2,2</sub></td><td>12.54</td><td>-</td><td>-</td><td>-</td></tr><tr><td>c<sub>2,0</sub></td><td>19.84</td><td>-</td><td>-</td><td>-</td></tr><tr><td rowspan="2">H<sub>4,i</sub></td><td>w<sub>3,0</sub></td><td>0.99</td><td>-</td><td>-</td><td>-</td></tr><tr><td>c<sub>3,0</sub></td><td>-2.03</td><td>-</td><td>-</td><td>-</td></tr><tr><td>H<sub>5,i</sub></td><td>w<sub>4,0</sub></td><td>2.30</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

Figure 3 shows that the model approximates more or less 40% the numerical calculations.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Results of the training of the artificial neural network (ANN), the results of void fraction of the ANN $\epsilon_{ANN}$ are compared to the void fraction calculated to $\epsilon_{fluent}$. The blue points represent the learning samples and the red points the validation samples. The red line is $\epsilon_{fluent} + 40\%$ and the blue line $\epsilon_{fluent} - 40\%$.</figcaption></figure>

### 2.2.2. Sensitivity Analysis of Void Fraction and Ohmic Resistance

In order to observe the effect of the different parameters on the electrolysis performance, the artificial neural network has been used to predict the void fraction $\epsilon$ and the electrolyte ohmic resistance $R$. Therefore, we performed a fractional factorial design of experiment (DOE). The DOE used is a “one factor at a time” (OFAT), referred to the minimal parameter-values curve (red). In all other curves, a single parameter among (T, Y, j, H, and h) changes from min to max value to assess the effect of each one on the void fraction and electrolyte ohmic resistance. Table 3 gives the minimum and maximum values of different parameters in this DOE. The bubble radius $r$ is set at 25 µm. Because the ANN
cannot evaluate the void fraction for $r^*$ above $1.5 \times 10^{-3}$, the electrode height $H_{elec}$ has been set to $5 \times 10^{-2}$.

<table><caption>**Table 3.** Minimum and maximum of the explored parameters for the sensitivity analysis.</caption><thead><tr><th>Variables</th><th>Min</th><th>Max</th></tr></thead><tbody><tr><td>$T$ (K)</td><td>293</td><td>363</td></tr><tr><td>$Y$ (-)</td><td>0.2</td><td>0.3</td></tr><tr><td>$h$ (m)</td><td>$4 \times 10^{-4}$</td><td>$10^{-3}$</td></tr><tr><td>$H_{elec}$ (m)</td><td>$5 \times 10^{-2}$</td><td>$10^{-1}$</td></tr><tr><td>$j$ (A m<sup>-2</sup>)</td><td>$10^3$</td><td>$10^4$</td></tr><tr><td>$p$ (bar)</td><td>1</td><td>30</td></tr></tbody></table>

The results of the predictions are represented in Figures 4–7. In the Figure 4a,b, it can be observed that the pressure decreases for all the cases the void fraction $\varepsilon$ at least until 5 bar and the ohmic resistance $R$. However, the artificial neural network predicts that the ohmic resistance reaches a plateau around 10 bar for most of the cases. In the case where $H_{max}$ = 10 cm, the minimum value is obtained near 5 bar but the void fraction rises again after this limit. This case, mass fraction and temperature seem to have little effect on the void fraction but have more important effects on the resistance.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Void fraction $\varepsilon_{ANN}$ (a) and ohmic resistance $R$ (b) depending on the pressure in [1–30] bar.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Void fraction $\varepsilon_{ANN}$ (a) and ohmic resistance $R$ (b) depending on the electrolyte thickness $h$ in [0.3;1.5] mm.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Void fraction $ε_{ANN}$ (a) and ohmic resistance R (b) depending current density in [1;10] kA m⁻².</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Void fraction $ε_{ANN}$ (a) and ohmic resistance R (b) depending on electrode height $H_{elec}$.</figcaption></figure>

In Figure 5a, it can be clearly seen that the electrolyte thickness has a great influence on the void fraction. By decreasing the electrolyte thickness *h* by few hundred micrometers, the void fraction increases abruptly from 0.1 to 0.6. Therefore, the existence of an optimum is explained. In the Figure 5b, the optimum of electrolyte thickness *h*<sub>opt</sub> can be observed. This optimum mainly depends on the electrode heights *H*<sub>elec</sub> and pressure *p*. The pressure reduces the radius of the bubble *r* and thus diminishes the optimal thickness of the electrolyte. Increasing the electrode height means increasing the injected gas flow so it increases the void fraction and the electrolyte ohmic resistance.

In Figure 6a,b, it can be observed that the current density increases the void fraction until 1500 A m⁻² and reached a plateau. The value of the plateau depends on the geometry of the cell especially the ratio electrolyte thickness/electrode height. However, the kinematic viscosity *ν* seems to play a role also.

In Figure 7a,b, the electrode height sensitivity impact on the ohmic resistance and void fraction can be observed. As expected, the electrode height increases the void fraction and thus the ohmic resistance. However, for *p*<sub>max</sub> (30 bar), it seems that electrode height does not impact the void fraction values. This is a surprising result that must be taken with caution.

In conclusion, the neural network gives roughly the sensitivity to all parameters even if for some case the results are surprising. Indeed, increasing the pressure must strictly decreases the void fraction and the void fraction should always increase according to electrode height. However, even if the sensitivities are strictly respected, the neural network approach save a lot of computational time. Indeed, to train the ANN, 61 CFD evaluation has been needed (corresponding to approximately 20 h of calculation on an Intel Core i7-6700HQ CPU @2.60 GHz). To make the Figures 4–7, at least 480
evaluations must be performed corresponding to approximately 30 s using the ANN and 160 CPU hours using the CFD code. Moreover, the genetic algorithm requires at the bare minimum 500 evaluations so the optimization would take 167 CPU hours if the CFD code is used against 30 CPU seconds using the ANN.

## 2.3. Economy

The hydrogen cost has two contribution: The CAPEX that represents the cost of the investments and the OPEX that represents the cost during electrolyzer operations (Equation (30)). CAPEX is estimated by using Faraday’s law, see Equation (31). The investment cost is divided by the total amount of hydrogen produced during the electrolyzer lifetime. The OPEX depends on the cell voltage and electricity cost.

$$OPEX = \frac{U_{cell} \cdot F}{3600} \cdot EC \qquad (30)$$

$$CAPEX = \left( \frac{2 F \cdot IC_S}{t \cdot M_{H_2}} \right) \cdot \frac{1}{j} \qquad (31)$$

$$F_c = CAPEX + OPEX \qquad (32)$$

With OPEX and CAPEX in € kg⁻¹, Ucell the cell voltage in V, F Faraday’s constant in A s mol⁻¹, EC the energy cost in € kWh⁻¹, ICs the initial capital cost per unit area in € m⁻², t the electrolyzer lifespan in s, M_H₂ is the molar mass in kg mol⁻¹.

The life of the electrolyzer, the initial capital cost per unit area and the cost of electricity cannot be calculated using a numerical model, but estimates do exist or can be calculated. Grigoriev et al. [18] have indicated that the life of the electrolyzer is approximately 90,000 h. The cost of electricity varies greatly depending on the country or the means of power generation used to produce electricity. In metropolitan France, the cost of electricity is about 0.17 € kWh⁻¹ and in Germany about 0.30 € kWh⁻¹. The effect of these two costs will be studied in the optimization section.

To estimate the ICs, the cost per kWh given by Grigoriev et al. [18] can be used. They stated that the cost per kW is between 500 and 1400 today and will be between 400 and 850 in 2030 and will reach the interval of 200–700 € kW⁻¹. According to Grigoriev et al. [18], the electrolyzer power density is between 4 and 10 kW m⁻². Thus, ICs will be explored between 800 and 14,000 € m⁻². The Figure 8 presents the OPEX, CAPEX, and total hydrogen cost for an electrolyzer depending on the current density for a naïve design (hH₂ = hO₂ = 1.5 mm, Helec = 10 cm, p = 1 bar, T = 353 K, Y = 0.30, ban = 0.15 V dec⁻¹, bcath = 0.15 V dec⁻¹, IC = 14,000 € m⁻², and EC = 0.17 € kWh⁻¹) optimum current density is 2923 A m⁻² and the minimum cost is 12.48 € kg⁻¹.

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Investment cost (CAPEX), operating cost (OPEX), and Cost function for an electrochemical cell with hH₂ = hO₂ = 1.5 mm, Helec = 10 cm, p = 1 bar, T = 353 K, Y = 0.30, ban = 0.15 V dec⁻¹, bcath = 0.15 V dec⁻¹, r = 25 µm, ICs = 14,000 € m⁻², and EC = 0.17 € kWh⁻¹.</figcaption></figure>
## 3. Set up of the Genetic Algorithm

### 3.1. Optimization Problem

We also define the cost function $F_c(CAPEX, OPEX)$ representing the addition of investment cost (CAPEX) and operating cost (OPEX) in € kg⁻¹.

$$F_c(CAPEX, OPEX) = CAPEX + OPEX \qquad (33)$$

The set of design parameters used in the rest of the article will be referred as DP.

$$DP = \{h_{H_2}, h_{O_2}, H_{elec}, j, p, T, Y\} \qquad (34)$$

This type of optimization problem is classified as a NP-complete (non-polynomial solving time) problem that has not yet found an effective deterministic solving algorithm.

$$ \text{Given } F_c(DP), \text{ find } DP^* \in (\tilde{E} \subseteq E) \text{ such as :} \\ F_c(DP^*) \le F_c(\tilde{DP}), \forall \tilde{DP} \in \tilde{E} \\ \text{with } g_i(\tilde{DP}) = 0 \text{ et } h_j(\tilde{DP}) \ge 0 \qquad (35) $$

With $DP^*$ the set of design parameters that gives the minimal cost value (i.e., the best design), $\tilde{E}$ the set of feasible solutions wrt $g_i$ and $h_j$ of constraint and $E$ the set of general solutions.

Taking into account the dimensionality, the non-linearity and the couplings in this problem, classical methods such as gradient descent are not likely to produce a successful optimization. Therefore, we propose to rely on stochastic optimization such as evolutionary algorithms. This method allows to find good solution when the problem is very difficult, and the search space is huge.

### 3.2. Genetic Optimization Algorithm (GA)

Historically, genetic algorithms (GAs) are the first and the simplest ones to implement [25]. GAs are particularly well adapted for parametric optimization as is our design problem (see Equation (36)).

First of all, the design parameters (DP) must be put under genetic material (chromosome). The fitness function is based on the chromosome and cost function. To begin the genetic algorithm, an initial population of individuals (i.e., design candidates) is randomly generated and their associated fitnesses are calculated. Next, individuals are selected for pairing according to their fitness to produce offspring. During the pairing, chromosome may undergo crossover or mutation. Eventually, we look if the best solution is satisfying if not, we generate a new population using the same process (evaluation, selection, cross-over, mutation). There are different types of genetic algorithms, but they can all be summarized in one general flowchart (Figure 9).
<figure><img src="image_12.png" /><figcaption>Figure 9. Flowchart of the genetic algorithm.</figcaption></figure>

### 3.3. Coding

The method of genetic algorithms takes its vocabulary from genetics. A design solution is called a genotype (composed of a single chromosome in the Simple Genetic Algorithm) and the parameters that make up this case are called genes. In our case, a chromosome containing all the genes would take the following form Equation (36).

$$ genotype = chromosome = \underbrace{[\bar{h}_{H2}, \bar{h}_{O2}, \bar{p}, \bar{H}_{elec}, \bar{j}, \overline{T,Y}]}_{\text{genes ($n_p$ parameters)}} \quad (36) $$

All genes are normalized using Equation (37) and therefore their value is between 0 and 1.

$$ \bar{x}_{norm} = \frac{x - x_{min}}{x_{max} - x_{min}} \quad (37) $$
The number of bits in a gene depends on the desired precision (or quantification) for a parameter. The input parameters for bits determination are summarized in Table 4. To determine the number of bits, Equations (38) and (39) are used.

<table><caption>**Table 4.** Design parameter their intervals and their resolution.</caption><thead><tr><th>Design Parameter</th><th>x<sub>min</sub></th><th>x<sub>max</sub></th><th>Q<sub>wtd</sub></th></tr></thead><tbody><tr><td>h<sub>H2</sub> (m)</td><td>1.5 × 10<sup>-4</sup></td><td>1.5 × 10<sup>-3</sup></td><td>1 × 10<sup>-4</sup></td></tr><tr><td>h<sub>H2</sub> (m)</td><td>1.5 × 10<sup>-4</sup></td><td>1.5 × 10<sup>-3</sup></td><td>1 × 10<sup>-4</sup></td></tr><tr><td>p (bar)</td><td>1</td><td>30</td><td>1</td></tr><tr><td>H<sub>elec</sub> (m)</td><td>5 × 10<sup>-2</sup></td><td>10<sup>-1</sup></td><td>5 × 10<sup>-3</sup></td></tr><tr><td>T (K)</td><td>298</td><td>358</td><td>10</td></tr><tr><td>j (A m<sup>-2</sup>)</td><td>10<sup>3</sup></td><td>10<sup>4</sup></td><td>10<sup>2</sup></td></tr><tr><td>Y (-)</td><td>0.2</td><td>0.4</td><td>2.5 × 10<sup>-2</sup></td></tr></tbody></table>

We need a number of bits for parameter DP<sub>i</sub>:

$$ \frac{x_{max,i} - x_{min,i}}{2^{n_{b,i}} - 1} = Q_{obti} \qquad (38) $$

$$ n_{b,i} \in \mathbb{N}^{+} \qquad n_{b,i} \geq \left( \ln \left( \frac{x_{max, i} - x_{min, i} + Q_{wtd,i}}{Q_{wtd,i}} \right) \frac{1}{\ln(2)} \right) \qquad (39) $$

Finally, the total chromosome number of bits is determined via the Equation:

$$ n_{bt} = \sum_{i}^{np} n_{b,i} \qquad (40) $$

Table 5 gives the number of bits, the wanted and obtained quantification for each design parameters. The total number of bits is 33.

<table><caption>**Table 5.** Design parameter and their quantification.</caption><thead><tr><th>Design Parameter</th><th>Q<sub>wtd</sub></th><th>n<sub>b</sub></th><th>Q<sub>obt</sub></th><th>Example</th><th>Coded Example</th><th>Decoded Example</th></tr></thead><tbody><tr><td>h<sub>H2</sub> (m)</td><td>10<sup>-4</sup></td><td>4</td><td>9 × 10<sup>-5</sup></td><td>3 × 10<sup>-4</sup></td><td>[0 0 0 1 1]</td><td>2.765 × 10<sup>-4</sup></td></tr><tr><td>h<sub>O2</sub> (m)</td><td>10<sup>-4</sup></td><td>4</td><td>9 × 10<sup>-5</sup></td><td>3 × 10<sup>-4</sup></td><td>[0 0 0 1 1]</td><td>2.765 × 10<sup>-4</sup></td></tr><tr><td>p (bar)</td><td>1</td><td>7</td><td>7.79 × 10<sup>-1</sup></td><td>3</td><td>[0 0 0 0 1 0]</td><td>2.54</td></tr><tr><td>H<sub>elec</sub> (m)</td><td>5 × 10<sup>-3</sup></td><td>4</td><td>3.33 × 10<sup>-3</sup></td><td>6 × 10<sup>-2</sup></td><td>[0 0 1 1]</td><td>5.94 × 10<sup>-2</sup></td></tr><tr><td>T (K)</td><td>10</td><td>3</td><td>7.85</td><td>3.08 × 10<sup>2</sup></td><td>[0 0 1]</td><td>3.04 × 10<sup>2</sup></td></tr><tr><td>j (A m<sup>-2</sup>)</td><td>10<sup>2</sup></td><td>7</td><td>7.08 × 10<sup>1</sup></td><td>2 × 10<sup>3</sup></td><td>[0 0 0 1 1 0]</td><td>1.98 × 10<sup>3</sup></td></tr><tr><td>Y (-)</td><td>2.5 × 10<sup>-2</sup></td><td>4</td><td>1.33 × 10<sup>-2</sup></td><td>0.3</td><td>[0 1 1 1]</td><td>2.87 × 10<sup>-1</sup></td></tr></tbody></table>

## 3.4. Evaluation

In optimization, we need a normalized function to be maximized, so we will use an inverse function. In order to have a solution between 0 and 1, the inverse function is scaled using the theoretical minimum of the CAPEX and the OPEX. The theoretical minimum of OPEX is obtained by neglecting the diphasic effect. The calculated minimum and maximum of CAPEX and OPEX are given in the Table 6.

$$ U_{cellmin} = E_{rev}(T_{max}, p_{min}) + \left( \frac{h_{H2min}}{\sigma_{max}} + \frac{h_{O2min}}{\sigma_{max}} \right) j_{min} + b_{cath} \log(j_{min}) + b_{an} \log(j_{min}) + R j_{min} \\ OPEX_{min} = \frac{U_{cellmin} F}{3600} EC \qquad (41) $$

$$ U_{cellmax} = E_{rev}(T_{min}, p_{max}) + \left( \frac{h_{H2max}}{\sigma_{min}} + \frac{h_{O2max}}{\sigma_{min}} \right) j_{max} + b_{cath} \log(j_{max}) + b_{an} \log(j_{max}) + R j_{max} \\ OPEX_{max} = \frac{U_{cellmax} F}{3600} EC \qquad (42) $$

$$ CAPEX_{min} = \left( \frac{2 F \cdot IC_S}{t \cdot M_{H_2}} \right) \cdot \frac{1}{j_{max}} \qquad (43) $$
$$ CAPEX_{max} = \left( \frac{2 F \cdot IC_S}{t \cdot M_{H_2}} \right) \cdot \frac{1}{j_{min}} \qquad (44) $$

<table><caption>**Table 6.** CAPEX and OPEX minimum and maximum and their associated value of fitness.</caption><thead><tr><th rowspan="4"></th><th rowspan="4"></th><th rowspan="4"></th><th colspan="4">EC (€ kWh<sup>-1</sup>)</th></tr><tr><th colspan="2">0.17</th><th colspan="2">0.3</th></tr><tr><th>OPEX<sub>min</sub></th><th>OPEX<sub>max</sub></th><th>OPEX<sub>min</sub></th><th>OPEX<sub>max</sub></th></tr><tr><th>9.48</th><th>13.46</th><th>16.74</th><th>23.76</th></tr></thead><tbody><tr><td rowspan="2">800</td><td>CAPEX<sub>min</sub></td><td>0.02</td><td>F&#773; = 1<br/>f = 1</td><td>F&#773; = 0.70<br/>f = 0.04</td><td>F&#773; = 1<br/>f = 1</td><td>F&#773; = 0.70<br/>f = 0.02</td></tr><tr><td>CAPEX<sub>max</sub></td><td>0.24</td><td>F&#773; = 0.98<br/>f = 0.92</td><td>F&#773; = 0.69<br/>f = 0</td><td>F&#773; = 0.99<br/>f = 0.96</td><td>F&#773; = 0.69<br/>f = 0</td></tr><tr><td rowspan="2">3000</td><td>CAPEX<sub>min</sub></td><td>0.09</td><td>F&#773; = 1<br/>f = 1</td><td>F&#773; = 0.71<br/>f = 0.12</td><td>F&#773; = 1<br/>f = 1</td><td>F&#773; = 0.71<br/>f = 0.07</td></tr><tr><td>CAPEX<sub>max</sub></td><td>0.89</td><td>F&#773; = 0.92<br/>f = 0.77</td><td>F&#773; = 0.67<br/>f = 0</td><td>F&#773; = 0.95<br/>f = 0.86</td><td>F&#773; = 0.68<br/>f = 0</td></tr><tr><td rowspan="2">IC<sub>s</sub> (€ m<sup>-2</sup>)</td><td>CAPEX<sub>min</sub></td><td>0.17</td><td>F&#773; = 1<br/>f = 1</td><td>F&#773; = 0.71<br/>f = 0.20</td><td>F&#773; = 1<br/>f = 1</td><td>F&#773; = 0.71<br/>f = 0.13</td></tr><tr><td>CAPEX<sub>max</sub></td><td>1.78</td><td>F&#773; = 0.85<br/>f = 0.61</td><td>F&#773; = 0.63<br/>f = 0</td><td>F&#773; = 0.91<br/>f = 0.74</td><td>F&#773; = 0.66<br/>f = 0</td></tr><tr><td rowspan="2">6000</td><td>CAPEX<sub>min</sub></td><td>0.27</td><td>F&#773; = 1<br/>f = 1</td><td>F&#773; = 0.71<br/>f = 0.27</td><td>F&#773; = 1<br/>f = 1</td><td>F&#773; = 0.71<br/>f = 0.18</td></tr><tr><td>CAPEX<sub>max</sub></td><td>2.68</td><td>F&#773; = 0.80<br/>f = 0.50</td><td>F&#773; = 0.60<br/>f = 0</td><td>F&#773; = 0.88<br/>f = 0.65</td><td>F&#773; = 0.64<br/>f = 0</td></tr><tr><td rowspan="2">9000</td><td>CAPEX<sub>min</sub></td><td>0.42</td><td>F&#773; = 1<br/>f = 1</td><td>F&#773; = 0.71<br/>f = 0.35</td><td>F&#773; = 1<br/>f = 1</td><td>F&#773; = 0.71<br/>f = 0.25</td></tr><tr><td>CAPEX<sub>max</sub></td><td>4.17</td><td>F&#773; = 0.72<br/>f = 0.37</td><td>F&#773; = 0.56<br/>f = 0</td><td>F&#773; = 0.82<br/>f = 0.53</td><td>F&#773; = 0.61<br/>f = 0</td></tr><tr><td>14,000</td><td>CAPEX<sub>min</sub></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>CAPEX<sub>max</sub></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

Inversing the cost function, having the lowest hydrogen cost ($F_c = F_{cmin}$) corresponding to the highest fitness ($\bar{F}_{max} = 1$) gives:

$$ \bar{F} = \frac{F_{cmin}}{F_c} = \frac{(CAPEX_{min} + OPEX_{min})}{OPEX + CAPEX} \qquad (45) $$

Unfortunately, $\bar{F}(F_{cmax})$, corresponding to the highest cost (then the lowest fitness $F_{min}$) is not 0, as it should be:

$$ \bar{F}_{min} = \frac{F_{cmin}}{F_{cmax}} = \frac{(CAPEX_{min} + OPEX_{min})}{(OPEX_{max} + CAPEX_{max})} \qquad (46) $$

Thus, it is interesting to make a linear transformation to take advantage of the whole range 0 to 1. With:

$$ \begin{aligned} A_0 &= \frac{1}{1 - \frac{CAPEX_{min} + OPEX_{min}}{CAPEX_{max} + OPEX_{max}}} \\ B_0 &= -A_0 \frac{CAPEX_{min} + OPEX_{min}}{CAPEX_{max} + OPEX_{max}} \end{aligned} \qquad (47) $$

The value of these coefficients depends on the initial assumptions. Finally, the Equations (48) and (49) gives the form of the fitness function.

$$ f(CAPEX, OPEX) = A_0 \bar{F} + B_0 \qquad (48) $$

$$ f(CAPEX, OPEX) = A_0 \frac{CAPEX_{min} + OPEX_{min}}{CAPEX + OPEX} + B_0 \qquad (49) $$
### 3.5. Operators

The Flowchart shown in the Figure 9 shows there are several operators: Initial random population generation, selection, crossing-over, mutation, and population replacement.

#### 3.5.1. Initial Random Population

The first generation is created from scratch, randomly drawing bits values for all parameters of each individuals in the population. The population includes $k_{pop}$ candidate solutions (individuals). The following generations will undergo successively the genetic operators all along the evolution that include $k_{gen}$ generations.

#### 3.5.2. Roulette Wheel Selection

To simulate the natural selection, the candidates from the population are selected using the roulette wheel selection [25]. Each candidate is given a portion of an imaginary wheel proportionally to its fitness (so the bigger fitness, the bigger portion). The wheel is spun $k_{pop}$ times ($k_{pop}$ being the total number of candidate solutions in the population). The selected candidates are paired to be mated. Each couple will be the parents producing offspring for the next generation.

#### 3.5.3. Single Point Crossing-Over

After the selection, a second operator that mimics the biological cell meiosis is used: The crossing-over. This operator is used to combine the parents' genetic material in the chromosome (genotype). A crossing-over is carried out with a probability occurrence of $p_c$. A crossing-over point (locus) is chosen randomly in the chromosome. Before this point, the chromosome is unchanged and beyond, the chromosome sequence is swapped between the parents. This produces two offspring that will be new candidates for the next generation. The artificial evolution if not very sensitive to the crossing-over probability [26]. It is why $p_c$ is generally set to 50%.

#### 3.5.4. Bitwise Mutation

The mutation operator is used to maintain genetic diversity in the population over the generations. It can be compared to the real biological mutation. In a living organism, the mutation occurs due to external factors (radioactivity, environmental turmoil) or just due to DNA recopying error. In the bitwise mutation, a mutation probability is drawn for each bit of the chromosome, the bits are flipped from 1 to 0 or from 0 to 1. The probability of mutation $p_m$ must be set very low to avoid too much genotype instability in the population. A good solution would be disrupted easily and its chance to survive mutation will be null. If $p_m$ reaches 50%, the genetic algorithm is equivalent to a random search, which must be prevented.

#### 3.5.5. Replacement Generator (Generational)

The children obtained from parents through selection, crossing-over and mutation replace totally the previous generation. A specific place in reserved for the best individual (in all previous generations). This requires scarifying the first newborn children or having an odd population size $k_{pop}$.

#### 3.5.6. Presentation of the Evolution Parameter

The evolution parameters have been introduced briefly in the previous subsection. In total, there are 4 of them:

* $k_{gen}$ the number of generations

* $k_{pop}$ the number of populations per generation

* $p_m$ probability of mutation

* $p_c$ probability of crossing-over is definitively set to 50%
The values of these parameters can greatly impact the results of the optimization simulation, they must be chosen wisely. The choice of $k_{gen}$ and $k_{pop}$ is constrained first by the time resource (CPU time) allocated to an optimization and then by the respect of infinitesimality constraint. Indeed, the advantage of the genetic algorithm is that it evaluates only a tiny (infinitesimal) part of the possible solution. To comply with this constraint:

$$ \frac{k_{gen}k_{pop}}{card(E)} \ll 1 \qquad (50) $$

$$ card(E) = 2^{n_{bt}} = 8.59 \times 10^{9} \qquad (51) $$

To validate the algorithm, we use the ad hoc parameters $k_{gen} = 50$ and $k_{pop} = 10$, we define the two dimensionless parameters:

$$ R_{GP} = \frac{k_{gen}}{k_{pop}} = 5 \qquad (52) $$

$$ P_{GP} = k_{gen}k_{pop} = 500 \ll Card(E). \qquad (53) $$

According to Chocron [26], $p_m$ must be chosen in such a way that, in average, one bit mutates per generation while having an average proportion of mutant individuals per generation equal to 50% (Equations (54) and (55)). By this setting, the aim is that the mutation should have the same impact as the crossing-over.

$$ E(nbmp) = k_{gen}E(nbmi) = k_{gen} \sum_{i=1}^{n_p} n_{bi} = \frac{k_{gen}}{2} \qquad (54) $$

$$ p_m = \frac{1}{2 \sum_{i=1}^{n_p} n_{bi}} = \frac{1}{66} = 0.015 = 1.55\% \qquad (55) $$

With $E(nbmp)$ and $E(nbmi)$ the expected value that a mutation affects the entire population and a single individual.

### 3.6. Validation of the Algorithm

Hypothesis:

* $IC_s = 800 \text{ € m}^{-2}$

* $EC = 17c\text{€ kWh}^{-1}$

* $b_{an} = b_{cath} = 0.15 \text{ V dec}^{-1}$

* $R_{memb} = 3.23 \times 10^{-5} \text{ Ohm m}^2$ [27]

The calculation has been performed using an Intel Core i7-6700HQ CPU @2.60 GHz. The optimization simulation time varies around 5 s. To validate the algorithm, the physics of the model presented in Section 2 will be simplified. The diphasic effect will not be taken into account. In doing so, we can analyze the equations to extract the sensitivities and find the optimal design if the process was not two-phase. The next step is to see if the algorithm finds a better candidate. By analyzing Equations (30) and (31), we can deduce that the current density must be around 1000. As temperature increases the electrical conductivity and decreases the reversible voltage, the algorithm must find the highest possible temperature. Since two-phase phenomena are not taken into account, the pressure proposed by the algorithm must be as low as possible (because it increases the reversible voltage). The same applies to the electrolyte thickness. Indeed, the latter decreases the ohmic resistance.

Figure 10 shows the evolution of the fitnesses (a) and the cost function (b) according to generations. In (a) we observe that $f_{best}$ and $f_{avg}$ increase greatly in the first 10 generations, we call this phase the exploration. During this phase, the GA searches good locations in the DP space. It finds very quickly good regions in DP, which explains the quick fitness growing, thus the cost dwindling.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Averaged over 100 evolutions (a) fitness and (b) cost function evolution according to the generation.</figcaption></figure>

Between generations 10 and 50 both fitnesses continue to increase but moderately, this in the exploitation phase. During this phase the algorithms attracts most of the population towards the best solution, which explain the reduction of $f_{best}$ and $f_{avg}$ difference. Usually, there exist a last phase where the $f_{avg}$ stagnates below the $f_{best}$ that continues to slightly increase until reaching a plateau around the end of evolution. This phase is called the convergence one and pursues the search around the best solutions to find slightest increases for $f_{best}$ (eventually ending up with the global optimum), maintaining $f_{avg}$ just below $f_{best}$ using the mutation as permanent genetic stirrer. Unfortunately, we cannot see this last phase that means the evolution is not completed.

In Figures 11 and 12 are presented the design parameters evolution path from random (first population) to overall best (last population) values. In (a), both *h* parameters ($H_2$ and $O_2$) converge to the smallest available value. This can be explained by the Ohm's law in the monophasic model.

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Design parameters evolution for the best DP set solution (in full DP space and for a single evolution) of the best candidate depending on the generation (a) for the anolyte and catholyte thickness $h_{H2}, h_{O2}$ (b) for electrolyte temperature *T*.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Design parameters evolution for the best DP set (in full DP space and for a single evolution) of the best candidate depending on the generation (a) current density *p*, (b) electrode height *H*<sub>elec</sub>, (c) pressure *p* and (d) KOH mass fraction *Y*.</figcaption></figure>

In Figures 11 and 12, we note that almost all DP values converge during the evolution towards their best value, but the parameter *H*<sub>elec</sub> (b). Actually, in this simplified model (monophasic), *H*<sub>elec</sub> has no influence on the cost function (gaz effects neglected).

This first tentative of artificial evolution of alkaline electrolyzer is promising, but we know that it
could be improved in two ways:

* The simplified model does not take into account the two-phase phenomena, which can greatly influence the hydrogen costs (*F*<sub>c</sub>).

* The evolution parameters used in this first try are not well adapted to the problem (missing convergence phase to achieve evolution).

In the next section we answer these two issues by using the more elaborate two-phase model and
by performing a sensitivity analysis of the genetic optimization w.r.t. the evolution parameters.

**4. Sensitivity Study of the Results to the Evolution Parameters**

In this section, we analyze the evolution algorithm on the more realistic diphasic problem to
tune the evolution parameters. The goal of this study is to adjust the evolution parameters (*k*<sub>gen</sub>, *k*<sub>pop</sub>,
$p_m$, $p_c$) in order to have the best genetic optimization performances. The hypothesis are the same than previously but because we take into account the diphasic effect, the bubble radius must be set. Thus, the following hypothesis are used.

* $IC_s = 800 \text{ € m}^{-2}$

* $EC = 17 \text{ € kWh}^{-1}$

* $b_{an} = b_{cath} = 0.15 \text{ V dec}^{-1}$

* $R_{memb} = 3.23 \times 10^{-5} \text{ Ohm m}^2$

* $r = 25 \times 10^{-6} \text{ m}$

### 4.1. Impact of $k_{gen}$ and $k_{pop}$

The first study should make it possible to determine the value of $k_{gen}$ and $k_{pop}$, the numerical experiments carried out are gathered in Table 7. $R_{GP}$ represents the ratio between $k_{gen}$ over $k_{pop}$. It allows to set (inversely) the level of “parallelism” of the genetic algorithm, to answer our high dimension and strongly coupled optimization problem. $P_{GP}$ is the product of $k_{gen}$ by $k_{pop}$ (i.e., the number of evaluations then of simulations), the higher this parameter is, the longer the evolution will take. To reduce the optimization computational costs, it is necessary to reach a compromise between fitness results and calculation time.

<table><caption>**Table 7.** Design of experiments for the two variables $P_{GP}$ and $R_{GP}$. The optimal number of generation $k_{gen}$ and population size $k_{pop}$ are in bold.</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="6">P<sub>GP</sub></th></tr><tr><th>500</th><th>1000</th><th>1500</th><th>2000</th><th>2500</th><th>3000</th></tr></thead><tbody><tr><td rowspan="2">1</td><td>k<sub>pop</sub></td><td>22</td><td>32</td><td>38</td><td>44</td><td>50</td><td>54</td></tr><tr><td>k<sub>gen</sub></td><td>23</td><td>31</td><td>40</td><td>45</td><td>50</td><td>56</td></tr><tr><td rowspan="2">3</td><td>k<sub>pop</sub></td><td>13</td><td>19</td><td>23</td><td>26</td><td>29</td><td>32</td></tr><tr><td>k<sub>gen</sub></td><td>38</td><td>53</td><td>66</td><td>76</td><td>85</td><td>93</td></tr><tr><td rowspan="2">5</td><td>k<sub>pop</sub></td><td>10</td><td>14</td><td>17</td><td>20</td><td>22</td><td>24</td></tr><tr><td>k<sub>gen</sub></td><td>50</td><td>70</td><td>87</td><td>100</td><td>112</td><td>122</td></tr><tr><td rowspan="2">7</td><td>k<sub>pop</sub></td><td>8</td><td>12</td><td>15</td><td>17</td><td>19</td><td>21</td></tr><tr><td>k<sub>gen</sub></td><td>59</td><td>84</td><td>102</td><td>118</td><td>132</td><td>145</td></tr><tr><td rowspan="2">10</td><td>k<sub>pop</sub></td><td>7</td><td>10</td><td>12</td><td>14</td><td>16</td><td>17</td></tr><tr><td>k<sub>gen</sub></td><td>71</td><td>100</td><td>123</td><td>141</td><td>158</td><td>173</td></tr></tbody></table>

The results that really interest us, the final best fitness, are summarized in Figure 13. It can be seen in Figure 13a that for any $R_{GP}$, it is necessary to choose a $P_{GP}$ of at least 1500 but for $R_{GP}$ of 1. However, in Figure 13b, we really observe a stagnation of the fitness value from $P_{GP} = 2000$ and $R_{GP} = 7$. If $P_{GP} = 3000$, the necessary $R_{GP}$ is around 3, and this $R_{GP}$ increases when $P_{GP}$ decreases. We can therefore conclude that the optimum setting is $P_{GP} = 2000$ is necessary for a $R_{GP} = 7$. If a larger $P_{GP}$ is available (if more CPU time is available), then we should decrease $R_{GP}$ to 3 for $P_{GP} = 3000$. This amounts to increasing the number of individuals in the population, i.e., taking advantage of resources to increase parallelism, which can be enhanced materially by using multiprocessors. So we choose $k_{gen} = 118$ and $k_{pop} = 17$.
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Result of the sensitivity study of the fitness to the $R_{GP}$ (a) and $P_{GP}$ (b).</figcaption></figure>

## 4.2. Impact of $p_m$

The mutation probability $p_m$ is important for the initial population draw, the exploration phase and the genetic stirring all along the evolution. If $p_m$ is too low, no new genetic material is found, and if too high it is disruptive and prevent keeping good solutions. Here we study its impact upon the fitnesses.

Initially, $p_m$ = was set at 1.55% as recommended by Chocron [26]. In order to refute or affirm this choice, a sensitivity study based on $p_m$ was conducted. In order to observe only the operator effect of the mutation, elitism was deactivated (i.e., the fact that the best individual is always re-injected into the population). The results of the evolutions are shown in Figure 14. The maximum fitness increases from 0 to about 0.5%, then reaches a plateau of up to 2% and then gradually decreases. For $p_m$ = [0.5–2]%, there is a peak at $p_m$ = 0.7%. However, the difference between the fitness value at $p_m$ = 0.7% and 1.5% is very small. The choice of $p_m$ recommended by Chocron is therefore relevant as far as we focus on the best fitness.

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Result of the sensitivity study of the fitness to the mutation rate. In (a) the whole research space has been plotted and (b) only the mutation rate between 0 and 2% has been plotted.</figcaption></figure>
In Figure 15, the evolution parameters are the ones just defined. We see in this 100-averaged run that the exploration phase stops at gen = 10, followed by the exploitation phase ending for gen = 90. This time, we observe a convergence phase (gen = 90 to 117) where $f_{avg}$ stagnates just below $f_{best}$ that increases slightly. The convergence gap between $f_{best}$ and $f_{avg}$ is only 3% of the $f_{best}$ value, as expected in the canonical GA (<5%). In conclusion, we can consider these evolution parameters to be validated.

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Averaged over 100 evolutions (a) Fitness and (b) cost function evolution according to the generation with $k_{gen}$ = 118, $k_{pop}$ = 17, $p_m$ = 0.5% ($IC_s$ = 800 € m⁻², $EC$ = 0.17 € kWh⁻¹, $b_{an}$ = $b_{cath}$ = 0.15 V dec⁻¹).</figcaption></figure>

Eventually, the evolution parameter are summarized in Table 8.

<table><caption>**Table 8.** Optimum evolution parameters.</caption><thead><tr><th></th><th>Value</th></tr></thead><tbody><tr><td>$k_{pop}$</td><td>17</td></tr><tr><td>$k_{gen}$</td><td>118</td></tr><tr><td>$p_c$</td><td>50%</td></tr><tr><td>$p_m$</td><td>0.5%</td></tr><tr><td>$R_{GP}$</td><td>7</td></tr><tr><td>$P_{GP}$</td><td>200</td></tr></tbody></table>

These evolution parameters allow to get the optimum DP set summarized in Table 9. The main results are the difference in catholyte and anolyte thickness and electrode heights. Due to the presence of gases, the minimum of electrolyte thickness cannot be chosen. The electrode height is set to the minimum value available. There is twice as much gas injected in the catholyte than in anolyte. Thus, the catholyte thickness $h_{H2}$ is bigger than the anolyte thickness $h_{H2}$. The lower the electrode height, the lesser gases injected resulting in a lower void fraction. The pressure stays at 1 bar. It can be explained by the fact that the surface coverage sensitivity regarding the pressure has not been taken into account. Thus, at relatively low current density, the pressure effect increases the cell voltage.
<table><caption>**Table 9.** Optimum DP with their associated hydrogen cost for $IC_s = 800 \ \text{€ m}^{-2}$, EC = 17 € kWh$^{-1}$, $r = 25 \ \mu\text{m}$ and $R_{\text{memb}} = 3.23 \times 10^{-5} \ \text{Ohm m}^2$.</caption><thead><tr><th></th><th>Value</th></tr></thead><tbody><tr><td>$h_{H2}$(m)</td><td>$3 \times 10^{-4}$</td></tr><tr><td>$h_{O2}$(m)</td><td>$2 \times 10^{-4}$</td></tr><tr><td>$H_{\text{elec}}$(m)</td><td>$5 \times 10^{-2}$</td></tr><tr><td>$j$ (A m$^{-2}$)</td><td>1000</td></tr><tr><td>$Y$ (-)</td><td>0.2</td></tr><tr><td>$T$ (K)</td><td>350</td></tr><tr><td>$p$ (bar)</td><td>1</td></tr><tr><td>CAPEX (€ kg$^{-1}$)</td><td>0.22</td></tr><tr><td>OPEX (€ kg$^{-1}$)</td><td>9.85</td></tr><tr><td>$F_c$ (€ kg$^{-1}$)</td><td>10.07</td></tr></tbody></table>

### 4.3. Optimization of the Naïve Solution

In the Section 2.3, we defined a naïve solution. In this section, we will compare this naïve solution with a GA optimized solution with the same hypothesis. The Figure 16a presents the average and best fitness according to the generations. The evolution is correctly achieved around $g = 97$. In the Figure 16b, it can be observed that the genetic optimization allows to decrease the random hydrogen cost from $12.70 \ \text{€ kg}^{-1}$ to $12.30 \ \text{€ kg}^{-1}$.

<figure><img src="image_19.png" /><figcaption>**Figure 16.** Averaged over 100 evolutions (a) Fitness and (b) cost function evolution according to the generation with $k_{\text{gen}} = 118$, $k_{\text{pop}} = 17$, $p_m = 0.5\%$ ($IC_s = 14,000 \ \text{€ m}^{-2}$, EC = 0.17 € kWh$^{-1}$, $b_{\text{an}} = b_{\text{cath}} = 0.15 \ \text{V dec}^{-1}$).</figcaption></figure>

The Table 10 summarized the DP of the two solutions (the naïve and optimized solutions) and their corresponding hydrogen cost. The optimization allowed to find the optimal $H_{\text{elec}}$, $h_{O2}$, $h_{H2}$ and $Y$. The electrode height is decreased at 5 cm and the anolyte and catholyte thickness are decreased at 400 µm. The KOH mass fraction is decreased to 0.23 (the loss of electrical conductivity is balanced by the decrease of void fraction allowed by a less viscous electrolyte). This optimization allows to decrease the price of 0.18 € kg$^{-1}$ (1% of decrease) and an increase the current density of 290 A m$^{-2}$ (10% of increase).
<table><caption>**Table 10.** Naïve and Optimum DP with their associated hydrogen cost for EC = 17 € kWh⁻¹, r = 25 µm and Rₘₐ₈₅ = 3.23 × 10⁻⁵ Ohm m².</caption><thead><tr><th></th><th>Naïve</th><th>GA1</th><th>GA2</th></tr></thead><tbody><tr><td>IC<sub>s</sub> (€ m<sup>-2</sup>)</td><td>1.4 × 10<sup>4</sup></td><td>1.4 × 10<sup>4</sup></td><td>8 × 10<sup>2</sup></td></tr><tr><td>h<sub>H2</sub> (m)</td><td>1.5 × 10<sup>-3</sup></td><td>4 × 10<sup>-4</sup></td><td>3 × 10<sup>-4</sup></td></tr><tr><td>h<sub>O2</sub> (m)</td><td>1.5 × 10<sup>-3</sup></td><td>4 × 10<sup>-4</sup></td><td>2 × 10<sup>-4</sup></td></tr><tr><td>H<sub>elec</sub> (m)</td><td>10<sup>-1</sup></td><td>5 × 10<sup>-2</sup></td><td>5 × 10<sup>-2</sup></td></tr><tr><td>j (A m<sup>-2</sup>)</td><td>2923</td><td>3214</td><td>1000</td></tr><tr><td>Y (-)</td><td>0.3</td><td>0.23</td><td>0.2</td></tr><tr><td>T (K)</td><td>350</td><td>350</td><td>350</td></tr><tr><td>p (bar)</td><td>1</td><td>1</td><td>1</td></tr><tr><td>CAPEX (€ kg<sup>-1</sup>)</td><td>1.42</td><td>1.30</td><td>0.22</td></tr><tr><td>OPEX (€ kg<sup>-1</sup>)</td><td>11.06</td><td>11</td><td>9.85</td></tr><tr><td>F<sub>c</sub> (€ kg<sup>-1</sup>)</td><td>12.48</td><td>12.30</td><td>10.07</td></tr></tbody></table>

The Figure 17 shows this result. It can be seen that the difference of price increase with an increasing current density until reaching a difference of 0.72 € kg⁻¹ at j = 10⁴ A m⁻².

<figure><img src="image_20.png" /><figcaption>**Figure 17.** CAPEX, OPEX, and total cost for the naïve solution (dashed line) and optimum solution (solid line) depending on the current density (ICs = 14,000 € m⁻², EC = 0.17 € kWh⁻¹, bₐ₅ = bₐ₆ = 0.15 V dec⁻¹).</figcaption></figure>

Figure 18a,b explains why there exist a difference of price. We can see that the ohmic overpotential is strongly decreased. In Figure 18b, it can be observed that the major cost of OPEX is mainly due to the activation overpotential and reversible voltage.
<figure><img src="image_21.png" /><figcaption>**Figure 18.** (a) Cell voltage of the naïve solution (dashed line) and optimum solution (solid line) (b) proportion of the different voltage of the cell voltage depending on the current density.</figcaption></figure>

The other overpotential are the same because they only are depending on the hypothesis.

**5. Conclusions**

A CFD model has been used to create an artificial neural network (ANN). This ANN allows the estimation of the void fraction very quickly from 166 CPU hours to around 30 s for 500 evaluations. However, further experimental studies must be performed to increase the accuracy of the CFD model. Indeed, the CFD model is based on only 3 experimental cases. Thus, the bubble diffusivity coefficient is not well correlated to the dimensionless parameters. With this new model, it would be possible to predict with accuracy the diphasic boundary layer thickness. The genetic algorithm parameters have been determined by making sensitivity analyses. To ensure a good optimization, the number of evolutions *P*<sub>GP</sub> = 2000 and a *R*<sub>GP</sub> = 7 have been chosen. For the mutation rate *p*<sub>m</sub>, the only based on the best fitness mutation rate (proposed by Chocron [26]), has been replaced by a more suitable *p*<sub>m</sub> = 0.5% (using also the average fitness). For the electrolyzer design, the algorithm preconizes the lowest height possible and find an optimum anolyte and catholyte thickness between 200 µm and 400 µm, depending on the optimum current density. It has been noticed that increasing the surface cost *IC*<sub>s</sub> increases the optimal current density. The optimum pressure has been determined at 1 bar but the surface coverage dependency according to the pressure (or the other parameters) has not been taken into account. Moreover, Bensmann et al. [8] showed that the optimal pressure depends on also of the compression, cooling, purity needed. Two further studies should be performed. The first one must focus on the determination of the surface coverage w.r.t dimensionless parameter (Reynolds number *Re*<sub>VG</sub>, Froude number *Fr*<sub>VG</sub>, dimensionless electrolyte thickness and radius *h*<sup>*</sup>, *r*<sup>*</sup>). Another one should address the influence of electrolyte and hydrogen temperature, pressure(s) of compression stage on the global hydrogen production system (cooling, compression, purification of hydrogen).

**Author Contributions:** Methodology, D.L.B., O.C. and P.M.; software, O.C. and P.K.; validation, D.L.B.; formal analysis, D.L.B., O.C. and P.K.; investigation, D.L.B.; data curation, D.L.B., O.C. and P.K.; writing—original draft preparation, D.L.B.; writing—review and editing, P.M., M.B., M.S., M.K., F.G. and R.I.; visualization, D.L.B. and O.C.; supervision, P.M.; 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.

## Appendix A

The model of Hammoudi et al. [10] used to calculate the reversible voltage in this written as Equations (A1)–(A4).

$$ E_{rev}(T, P) = E_{rev}(T, 1\text{bar}) + \frac{RT}{nF}\ln\left((P - Pw)^{1.5} \frac{Pw^*}{Pw}\right) + (P - Pw)\left(21.661 \times 10^{-6} - \frac{5.471 \times 10^{-3}}{T}\right) + (P - Pw)^2\left(-\frac{6.289 \times 10^{-6}}{T} + \frac{0.135 \times 10^{-3}}{T^{1.5}} + \frac{2.547 \times 10^{-3}}{T^2} - \frac{0.4825}{T^3}\right) \quad (A1) $$

$$ E_{rev}(T, 1\text{bar}) = 1.50342 - 9.956 \times 10^{-4} \times T + 2.5 \times 10^{-7} \times T^2 \quad (A2) $$

$$ Pw = T^{-3.498} \exp\left(37.93 - \frac{6426.32}{T}\right) \exp\left(0.016214 - 0.13802 \times C + 0.1933 \times C^{\frac{1}{2}}\right) \quad (A3) $$

$$ Pw^* = T^{-3.4159} \exp\left(37.043 - \frac{6275.7}{T}\right) \quad (A4) $$

## Appendix B

The following numerical experiments has been performed to train the ANN.

<table><caption>**Table A1.** Training set.</caption><thead><tr><th>N°</th><th><span style="text-decoration: overline;">Re</span><sub>V<sub>G</sub></sub></th><th><span style="text-decoration: overline;">Fr</span><sub>V<sub>G</sub></sub></th><th><span style="text-decoration: overline;">h</span><sup>*</sup></th><th><span style="text-decoration: overline;">r</span><sup>*</sup></th><th>ε<sub>fluent</sub></th></tr></thead><tbody><tr><td>1</td><td>-0.375</td><td>1</td><td>0.625</td><td>-0.25</td><td>2.07 × 10<sup>-4</sup></td></tr><tr><td>2</td><td>-0.875</td><td>-0.5</td><td>0.75</td><td>0.125</td><td>3.52 × 10<sup>-4</sup></td></tr><tr><td>3</td><td>-0.75</td><td>-0.125</td><td>-0.875</td><td>-0.5</td><td>1.95 × 10<sup>-3</sup></td></tr><tr><td>4</td><td>-0.625</td><td>0.25</td><td>-0.375</td><td>1</td><td>3.868 × 10<sup>-4</sup></td></tr><tr><td>5</td><td>0.5</td><td>0.875</td><td>-0.125</td><td>-0.75</td><td>3.972 × 10<sup>-4</sup></td></tr><tr><td>6</td><td>1</td><td>-0.375</td><td>-0.25</td><td>0.625</td><td>4.543 × 10<sup>-4</sup></td></tr><tr><td>7</td><td>0.25</td><td>-0.625</td><td>1</td><td>-0.375</td><td>3.884 × 10<sup>-4</sup></td></tr><tr><td>8</td><td>0.125</td><td>0.75</td><td>0.5</td><td>0.875</td><td>1.737 × 10<sup>-4</sup></td></tr><tr><td>9</td><td>0</td><td>0</td><td>0</td><td>0</td><td>3.624 × 10<sup>-4</sup></td></tr><tr><td>10</td><td>0.375</td><td>-1</td><td>-0.625</td><td>0.25</td><td>1.477 × 10<sup>-1</sup></td></tr><tr><td>11</td><td>0.875</td><td>0.5</td><td>-0.75</td><td>-0.125</td><td>8.782 × 10<sup>-4</sup></td></tr><tr><td>12</td><td>0.75</td><td>0.125</td><td>0.875</td><td>0.5</td><td>1.897 × 10<sup>-4</sup></td></tr><tr><td>13</td><td>0.625</td><td>-0.25</td><td>0.375</td><td>-1</td><td>5.931 × 10<sup>-4</sup></td></tr><tr><td>14</td><td>-0.5</td><td>-0.875</td><td>0.125</td><td>0.75</td><td>7.161 × 10<sup>-4</sup></td></tr><tr><td>15</td><td>-1</td><td>0.375</td><td>0.25</td><td>-0.625</td><td>5.067 × 10<sup>-3</sup></td></tr><tr><td>16</td><td>-0.25</td><td>0.625</td><td>-1</td><td>0.375</td><td>3.289 × 10<sup>-3</sup></td></tr><tr><td>17</td><td>-0.125</td><td>-0.75</td><td>-0.5</td><td>-0.875</td><td>1.515 × 10<sup>-3</sup></td></tr><tr><td>18</td><td>1</td><td>-1</td><td>-1</td><td>1</td><td>1</td></tr><tr><td>19</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1.219 × 10<sup>-4</sup></td></tr><tr><td>20</td><td>-1</td><td>-1</td><td>1</td><td>1</td><td>1.856 × 10<sup>-1</sup></td></tr><tr><td>21</td><td>1</td><td>-1</td><td>1</td><td>-1</td><td>2.096 × 10<sup>-2</sup></td></tr><tr><td>22</td><td>-1</td><td>1</td><td>-1</td><td>1</td><td>1.063 × 10<sup>-1</sup></td></tr><tr><td>23</td><td>-1</td><td>1</td><td>1</td><td>-1</td><td>2.935 × 10<sup>-3</sup></td></tr><tr><td>24</td><td>1</td><td>1</td><td>-1</td><td>-1</td><td>3.986 × 10<sup>-3</sup></td></tr><tr><td>25</td><td>-1</td><td>-1</td><td>-1</td><td>-1</td><td>1</td></tr><tr><td>26</td><td>-0.9325</td><td>-0.98</td><td>-0.95</td><td>-1</td><td>4.058 × 10<sup>-2</sup></td></tr><tr><td>27</td><td>-0.9325</td><td>-1</td><td>-0.875</td><td>1</td><td>6.778 × 10<sup>-1</sup></td></tr><tr><td>28</td><td>-0.99</td><td>-0.99</td><td>-0.95</td><td>1</td><td>9.828 × 10<sup>-2</sup></td></tr><tr><td>29</td><td>-0.875</td><td>-0.99</td><td>-0.875</td><td>-1</td><td>2.301 × 10<sup>-2</sup></td></tr><tr><td>30</td><td>-0.99</td><td>-0.98</td><td>-0.875</td><td>0</td><td>4.274 × 10<sup>-2</sup></td></tr><tr><td>31</td><td>-0.875</td><td>-1</td><td>-0.95</td><td>0</td><td>8.424 × 10<sup>-1</sup></td></tr></tbody></table>
<table><caption>Table A1. *Cont.*</caption><thead><tr><th>N°</th><th>Re<sub>V<sub>G</sub></sub></th><th>Fr<sub>V<sub>G</sub></sub></th><th>h<sup>*</sup></th><th>r<sup>*</sup></th><th>ε<sub>fluent</sub></th></tr></thead><tbody><tr><td>32</td><td>-0.946</td><td>-0.990</td><td>-0.728</td><td>-0.427</td><td>1.460 × 10<sup>-2</sup></td></tr><tr><td>33</td><td>-0.983</td><td>-0.997</td><td>-0.704</td><td>-0.200</td><td>4.440 × 10<sup>-2</sup></td></tr><tr><td>34</td><td>-0.973</td><td>-0.996</td><td>-1.023</td><td>-0.578</td><td>4.539 × 10<sup>-1</sup></td></tr><tr><td>35</td><td>-0.964</td><td>-0.994</td><td>-0.925</td><td>0.328</td><td>6.174 × 10<sup>-2</sup></td></tr><tr><td>36</td><td>-0.881</td><td>-0.991</td><td>-0.876</td><td>-0.729</td><td>2.390 × 10<sup>-2</sup></td></tr><tr><td>37</td><td>-0.844</td><td>-0.997</td><td>-0.900</td><td>0.101</td><td>3.623 × 10<sup>-2</sup></td></tr><tr><td>38</td><td>-0.900</td><td>-0.998</td><td>-0.655</td><td>-0.503</td><td>1.993 × 10<sup>-2</sup></td></tr><tr><td>39</td><td>-0.909</td><td>-0.992</td><td>-0.753</td><td>0.252</td><td>1.350 × 10<sup>-2</sup></td></tr><tr><td>40</td><td>-0.918</td><td>-0.995</td><td>-0.851</td><td>-0.276</td><td>2.892 × 10<sup>-2</sup></td></tr><tr><td>41</td><td>-0.890</td><td>-0.999</td><td>-0.974</td><td>-0.125</td><td>8.122 × 10<sup>-1</sup></td></tr><tr><td>42</td><td>-0.853</td><td>-0.993</td><td>-0.999</td><td>-0.352</td><td>1.044 × 10<sup>-1</sup></td></tr><tr><td>43</td><td>-0.863</td><td>-0.994</td><td>-0.679</td><td>0.026</td><td>1.155 × 10<sup>-2</sup></td></tr><tr><td>44</td><td>-0.872</td><td>-0.996</td><td>-0.778</td><td>-0.881</td><td>2.310 × 10<sup>-2</sup></td></tr><tr><td>45</td><td>-0.955</td><td>-0.999</td><td>-0.827</td><td>0.177</td><td>9.900 × 10<sup>-2</sup></td></tr><tr><td>46</td><td>-0.992</td><td>-0.993</td><td>-0.802</td><td>-0.654</td><td>6.528 × 10<sup>-2</sup></td></tr><tr><td>47</td><td>-0.937</td><td>-0.992</td><td>-1.048</td><td>-0.050</td><td>6.000 × 10<sup>-1</sup></td></tr><tr><td>48</td><td>-0.927</td><td>-0.998</td><td>-0.95</td><td>-0.805</td><td>1.431 × 10<sup>-1</sup></td></tr><tr><td>49</td><td>-0.8</td><td>-0.87</td><td>-0.85</td><td>-0.875</td><td>5.296 × 10<sup>-3</sup></td></tr><tr><td>50</td><td>-0.8</td><td>-0.97</td><td>-0.55</td><td>-0.125</td><td>3.516 × 10<sup>-3</sup></td></tr><tr><td>51</td><td>-0.9</td><td>-0.92</td><td>-0.85</td><td>-0.125</td><td>6.844 × 10<sup>-3</sup></td></tr><tr><td>52</td><td>-0.7</td><td>-0.92</td><td>-0.55</td><td>-0.875</td><td>2.807 × 10<sup>-3</sup></td></tr><tr><td>53</td><td>-0.9</td><td>-0.97</td><td>-0.7</td><td>-0.875</td><td>7.510 × 10<sup>-3</sup></td></tr><tr><td>54</td><td>-0.7</td><td>-0.87</td><td>-0.7</td><td>-0.125</td><td>2.482 × 10<sup>-3</sup></td></tr><tr><td>55</td><td>-0.9</td><td>-0.87</td><td>-0.55</td><td>-0.5</td><td>2.279 × 10<sup>-3</sup></td></tr><tr><td>56</td><td>-0.7</td><td>-0.97</td><td>-0.85</td><td>-0.5</td><td>1.068 × 10<sup>-2</sup></td></tr><tr><td>57</td><td>-0.8</td><td>-0.92</td><td>-0.7</td><td>-0.5</td><td>3.584 × 10<sup>-3</sup></td></tr><tr><td>58</td><td>-0.8</td><td>-0.87</td><td>-0.85</td><td>-1</td><td>1.740 × 10<sup>-2</sup></td></tr><tr><td>59</td><td>-0.8</td><td>-0.97</td><td>-0.55</td><td>-0.9</td><td>9.525 × 10<sup>-3</sup></td></tr><tr><td>60</td><td>-0.7</td><td>-0.92</td><td>-0.85</td><td>-0.9</td><td>1.082 × 10<sup>-2</sup></td></tr><tr><td>61</td><td>-0.9</td><td>-0.92</td><td>-0.55</td><td>-1</td><td>1.636 × 10<sup>-2</sup></td></tr><tr><td>62</td><td>-0.7</td><td>-0.97</td><td>-0.7</td><td>-1</td><td>3296 × 10<sup>-2</sup></td></tr><tr><td>63</td><td>-0.9</td><td>-0.87</td><td>-0.7</td><td>-0.9</td><td>3899 × 10<sup>-3</sup></td></tr><tr><td>64</td><td>-0.7</td><td>-0.87</td><td>-0.55</td><td>-0.95</td><td>7417 × 10<sup>-3</sup></td></tr><tr><td>65</td><td>-0.9</td><td>-0.97</td><td>-0.85</td><td>-0.95</td><td>1508 × 10<sup>-2</sup></td></tr></tbody></table>

## References

1. Le Bideau, D.; Mandin, P.; Benbouzid, M.; Kim, M.; Sellier, M.; Ganci, F.; Inguanta, R. Eulerian Two-Fluid Model of Alkaline Water Electrolysis for Hydrogen Production. *Energies* **2020**, *13*, 3394. [CrossRef]

2. 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]

3. Mat, M. A two-phase flow model for hydrogen evolution in an electrochemical cell. *Int. J. Hydrog. Energy* **2004**, *29*, 1015–1023. [CrossRef]

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

5. Wedin, R.; Dahlkild, 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]

6. Ulleberg, O. Modeling of advanced alkaline electrolyzers: A system simulation approach. *Int. J. Hydrog. Energy* **2003**, *28*, 21–33. [CrossRef]

7. Villagra, A.; Millet, P. An analysis of PEM water electrolysis cells operating at elevated current densities. *Int. J. Hydrog. Energy* **2019**, *44*, 9708–9717. [CrossRef]

8. Bensmann, B.; Hanke-Rauschenbach, R.; Müller-Syring, G.; Henel, M.; Sundmacher, K. Optimal configuration and pressure levels of electrolyzer plants in context of power-to-gas applications. *Appl. Energy* **2016**, *167*, 107–124. [CrossRef]

9. Chocron, O. Evolutionary design of modular robotic arms. *Robotica* **2008**, *26*, 323–330. [CrossRef]
10. Hammoudi, M.; Henao, C.; Agbossou, K.; Dubé, Y.; Doumbia, M. New multi-physics approach for modelling and design of alkaline electrolyzers. *Int. J. Hydrog. Energy* **2012**, *37*, 13895–13913. [CrossRef]

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

12. Hine, F.; Murakami, K. Bubble Effects on the Solution IR Drop in a Vertical Electrolyzer under Free and Forced Convection. *J. Electrochem. Soc.* **1980**, *127*, 292–297. [CrossRef]

13. Vogt, H. The Quantities Affecting the Bubble Coverage of Gas-Evolving Electrodes. *Electrochim. Acta* **2017**, *235*, 495–499. [CrossRef]

14. Jannsen, L.J.J.; Sillen, C.W.M.P.; Van Stralen, S.J.D. Bubble behaviour during oxygen and hydrogen evolution at transparent electrodes in KOH solution. *Electrochim. Acta* **1984**, *29*, 633–642. [CrossRef]

15. Gilliam, R.; Graydon, J.; Kirk, D.; Thorpe, S. A review of specific conductivities of potassium hydroxide solutions for various concentrations and temperatures. *Int. J. Hydrog. Energy* **2007**, *32*, 359–364. [CrossRef]

16. See, D.M.; White, R.E. Temperature and Concentration Dependence of the Specific Conductivity of Concentrated Solutions of Potassium Hydroxide. *J. Chem. Eng. Data* **1997**, *42*, 1266–1268. [CrossRef]

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

18. Grigoriev, S.; Fateev, V.; Bessarabov, D.; Millet, P. Current status, research trends, and challenges in water electrolysis science and technology. *Int. J. Hydrog. Energy* **2020**, *45*, 26036–26058. [CrossRef]

19. Nagai, N. Existence of optimum space between electrodes on hydrogen production by water electrolysis. *Int. J. Hydrog. Energy* **2003**, *28*, 35–41. [CrossRef]

20. Nagai, N.; Ito, T.; Nishijiri, N. Visualization and Numerical Simulation of Two-Phase Flow between Electrodes on Alkaline Water Electrolysis. *J. Flow Vis. Image Process.* **2010**, *17*, 99–112. [CrossRef]

21. Nagai, N.; Takeuchi, M.; Furuta, T. Effects of Bubbles Between Electrodes on Alkaline Water Electrolysis Efficiency Under Forced Convection of Electrolyte. In Proceedings of the 16th World Hydrogen Energy Conference, Lyon, France, 13–16 June 2006.

22. Vogt, H.; Balzer, R. The bubble coverage of gas-evolving electrodes in stagnant electrolytes. *Electrochim. Acta* **2005**, *50*, 2073–2079. [CrossRef]

23. 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]

24. Yu, H.; Wilamowski, B.M. Levenberg-Marquardt Training. In *The Industrial Electronics Handbook*; CRC Press: Boca Raton, FL, USA, 2011.

25. Goldberg, D.E. *Genetic Algorithms in Search, Optimization and MachineLearning*; Addison-Wesley Professional: Boston, MA, USA, 1989.

26. Chocron, O. Conception Evolutionnaire de Systèmes Robotiques. Ph.D. Thesis, Université Paris, Paris, France, 2000.

27. Seetharaman, S.; Davidson, D.J.; Vasudeva, S.; Sozhan, G. Polyvinyl Alcohol Based Membrane as Separator for Alkaline Water Electrolyzer. *Separ. Sci. Tech.* **2011**, *46*, 1563–1570. [CrossRef]

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

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