A numerical study on effects of current density distribution, turbulence, and magnetohydrodynamics (MHD) on electrolytic gas flow with application to alkaline water electrolysis (AWE)

Chemical Engineering Research and Design 208 (2024) 731–739 Available online 20 July 2024 0263-8762/© 2024 The Author(s). Published by Elsevier Ltd on behalf of Institution of Chemical Engineers. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). A numerical study on effects of current density distribution, turbulence, and magnetohydrodynamics (MHD) on electrolytic gas flow with application to alkaline water electrolysis (AWE) E. Karimi-Sibaki a,b , A. Vakhrushev a , M. Wu b , J. Bohacek c , A. Kharicha a,b,* a Christian-Doppler Laboratory for Metallurgical Applications of Magnetohydrodynamics, Montanuniversitaet of Leoben, Franz-Josef-Str. 18, Leoben A-8700, Austria b Chair of Simulation and Modeling of Metallurgical Processes, Montanuniversitaet of Leoben, Franz-Josef-Str. 18, Leoben A-8700, Austria c Heat Transfer and Fluid Flow Laboratory, Faculty of Mechanical Engineering, Brno University of Technology, Technicka 2896/2, Brno 616 69, Czech Republic ARTICLE INFO Keywords: Alkaline Water Electrolysis (AWE) Electrolytic gas flow Three-phase Eulerian model Magnetohydrodynamics (MHD) Numerical simulation Multiphase velocity field ABSTRACT A three-phase Eulerian model is proposed to investigate the induced flow due to the generation of gas bubbles between two parallel plates without forced convection with application to alkaline water electrolysis (AWE). Earlier models, assuming a laminar regime, accurately predicted the multiphase flow near electrodes but struggled to calculate bulk liquid electrolyte flow away from them. Herein, we study the influences of electric current density distribution, turbulence effects, and the interaction between flow and the magnetic field known as magnetohydrodynamics (MHD). Based on our modeling results, the traditional method using an averaged uniform current density along electrodes (e.g. here 2000 A m −2 ) is feasible, as incorporating calculated nonuniform current distribution minimally affects the multiphase velocity field. The Lorentz force, originating from flow interaction with the (self-induced) magnetic field, is negligible compared to forces like drag or bubble dispersion. Consequently, MHD effects only become relevant upon introducing an external magnetic field. Including turbulence in the model, being minor in magnitude but non-negligible, significantly improves the predicted velocity profile. Modeling results are validated against an experiment. 1. Introduction Generation and transport of gas bubbles, such as oxygen or hydrogen through electrolysis, plays a pivotal role in various chemical and metallurgical processes such as chloro-alkali (Karlsson and Cornell, 2016), sodium chlorate (Endr˝ odi et al., 2017), electroslag remelting (Kharicha et al., 2018), aluminium electrolysis (Cooksey et al., 2008; Sun et al., 2023), and water electrolysis (Angulo et al., 2020). Water electrolysis, in particular, has prompted extensive interest as a net-zero method for producing hydrogen, which is a clean and versatile energy carrier crucial for transitioning towards a low-carbon economy. Prominent water electrolysis technologies are alkaline water electrolysis (AWE), proton exchange membrane water electrolysis (PEM), and solid oxide electrolysis cell (SOEC) (Chi and Yu, 2018). Among these processes, AWE stands out as a mature technology for large-scale green hydrogen production (Hu et al., 2022). In AWE, water undergoes electrolysis with electricity to produce hydrogen and oxygen gases in the presence of an alkaline electrolyte like Potassium Hydroxide (KOH), facilitating the hydrogen evolution reaction (HER) at the cathode and the oxygen evolution reaction (OER) at the anode (see Fig. 1). The complex multiphase nature of gas-liquid flows in AWE presents significant challenges for understanding and enhancing the process. Accordingly, extensive experiments and modeling studies have been conducted over the past decades (Hreiz et al., 2015a; Daoudi and Bounahmidi, 2024). Computational modeling, as a powerful tool to unravel the underlying dynamics of the gas-liquid two-phase flow, provides insights that are often difficult to obtain through sole experiments. Modeling approaches are divided into Euler-Lagrange, Euler-Euler (Mixture), and Euler-Euler (Two-Fluid) (Hreiz et al., 2015a). The Euler-Lagrange model treats the liquid electrolyte as a continuous phase and gas bubbles as discrete entities (Nierhaus et al., 2009; Mandin et al., 2005; Hreiz et al., 2015b). In contrast, the Euler-Euler models, including both Mixture and Two-Fluid approaches, focus on understanding the interplay between these phases (Rivera et al., 2021). * Corresponding author at: Christian-Doppler Laboratory for Metallurgical Applications of Magnetohydrodynamics, Montanuniversitaet of Leoben, Franz-Josef-Str. 18, Leoben A-8700, Austria. E-mail address: [email protected] (A. Kharicha). Contents lists available at ScienceDirect Chemical Engineering Research and Design journal homepage: www.elsevier.com/locate/cherd https://doi.org/10.1016/j.cherd.2024.07.042 Received 25 April 2024; Received in revised form 20 June 2024; Accepted 15 July 2024 Chemical Engineering Research and Design 208 (2024) 731–739 732 The Mixture model uses the relative (slip) velocity between the liquid and gas phases to balance pressure and viscous drag forces, simplifying the complexity of two-phase interactions by treating them as a combined single-phase mixture (Hreiz et al., 2015b). Dahlkild (2001), Wedin and Dahlkild (2001) suggested hydrodynamic terms considering self, shear-induced, and Saffman’s force to modify the slip velocity term. Schillings et al. (2015) proposed a comprehensive Mixture model that exploits a thermal buoyancy-driven flow analogy to describe different bubble dispersion mechanisms. These mechanisms are suggested to be strongly dependent on flow conditions, imposed current density, and physical parameters such as bubble diameter and electrolyte viscosity. Essentially, the Mixture model facilitates computational efficiency by treating the phases as a combined single-phase mixture. However, the Mixture model may struggle to accurately capture phase distribution, especially in scenarios with significant phase segregation. Thus, the Two-Fluid model has become the most popular Eulerian volume-averaged approach to study the dynamics of gas-liquid two-- phase flow (Hreiz et al., 2015b). Several studies have utilized the Two-Fluid model (Daoudi and Bounahmidi, 2024; Hreiz et al., 2015b). A shared pressure field is considered in the Two-Fluid approach, which employs two continuity and two momentum equations to calculate volume fractions and velocity fields for each phase (Zarghami et al., 2020). Notable studies utilizing the Two-Fluid model are discussed in the following. Aldas (2004), Mat (2004), Aldas et al. (2008) incorporated a term into the continuity equations of two-phase flows to describe bubble dispersion, considering scenarios with the assumption of turbulent or laminar flows. Nomenclature AMagnetic vector potential, V sm−1 BMagnetic field, T CdDrag coefficient, - C1 ε Coefficient for turbulent dissipation, - C2 ε Coefficient for turbulent dissipation, - C μ Coefficient for turbulent viscosity, - dbBubble diameter, m FFaraday constant, A smol−1 FdVolumetric drag force, Nm−3 Fbd Volumetric bubble dispersion force, Nm−3 FLVolumetric Lorentz force, Nm−3 Gk,mSource turbulent kinetic energy, kgm−1s−3 gGravity constant, ms−2 H2Hydrogen, - IIdentity matrix, - ITurbulent intensity, % jElectric current density, Am−2 j0Apparent exchange current density, Am−2 kiTurbulent kinetic energy of i-th phase, m2s−3 kTurbulent kinetic energy of mixture, m2s−2 KgCoefficient for dispersion force, m2s−1 lLiquid electrolyte, - MH2Molecular weight Hydrogen, gmol−1 MO2Molecular weight Oxygen, gmol−1 nH2Hydrogen stoichiometric to electron ratio, - nO2Oxygen stoichiometric to electron ratio, - O2Oxygen, - pPressure, Pa RUniversal gas constant, JK−1mol−1 SiVolumetric mass source i-th phase, kgm−3s−1 SH2Volumetric mass source Hydrogen, kgm−3s−1 SO2Volumetric mass source Oxygen, kgm−3s−1 tTime, s TTemperature, K uiVelocity of i-th phase, ms−1 ugVelocity of gas phase, ms−1 ulVelocity of liquid phase, ms−1 umVelocity of mixture, ms−1 ∇TuiTranspose i-th phase velocity gradient, s−1 ∇TumTranspose mixture velocity gradient, s−1 Velec Applied potential at electrode, V WsWidth of adjacent computational cell, m X,YCoordinates, m α iVolume fraction of i-th phase, - α gVolume fraction of gas, - α H2Volume fraction of hydrogen, - α O2Volume fraction of oxygen, - ε Turbulent dissipation rate of mixture, m2s−3 μ iViscosity of i-th phase, kgs−1m−1 μ t iTurbulent viscosity of i-th phase, kgs−1m−1 μ t mTurbulent viscosity of mixture, kgs−1m−1 μ 0Magnetic permeability, Hm−1 ρ iDensity of i-th phase, kgm−3 ρ lDensity of liquid, kgm−3 ρ mDensity of mixture, kgm−3 σ Electrical conductivity, Sm−1 σ 0Electrical conductivity of liquid, Sm−1 σ kTurbulent kinetic energy Prandtl number, - σε Turbulent dissipation rate Prandtl number, - τ iShear stress tensor of i-th phase, Pa τ Re iTurbulent shear stress tensor of i-th phase, Pa τ Re mTurbulent shear stress tensor of mixture, Pa φElectric potential, V φcElectric potential adjacent to electrode, V Fig. 1. Schematic representation of the computational domain, including dimensions and boundaries. Illustratively, the Y direction is scaled down for clarity. E. Karimi-Sibaki et al. Chemical Engineering Research and Design 208 (2024) 731–739 733 Charton et al. (2009) developed a 3D steady-state model including the effect of turbulence, considering that the gas-liquid flow undergoes forced convection by pumping the electrolyte between electrodes. Alexiadis et al. (2011a, 2011b, 2012a, 2012b) developed 2D and 3D transient models assuming laminar flow, which consider the free surface of the liquid positioned at the end of electrodes and allow for the recirculation of gas bubbles between the electrodes. Also, assuming laminar flow and allowing for the recirculation of gas bubbles between electrodes, Alexiadis et al. (2012b), Abdelouahed et al. (2014a), Abdelouahed et al. (2014b) developed 3D models to study the evolution of oxygen gas bubbles. El-Askary et al. (2015) investigated the geometrical aspects and conditions of forced convection of the electrolyte, including turbulence. They reported that decreasing the main flow velocity, increasing the current density, and reducing the gap distance between the cathode and the anode improves the hydrogen production process. Detailed analysis of drag force and turbulent dispersion force considering various turbulent models for gas-liquid multiphase flow was reported by Zarghami et al. (2020). Also, including the turbulent model considering forced convection of the electrolyte, Rodríguez and Amores (2020) studied the impacts of gap width, temperature, and electrical conductivity of the solution. Le Bideau et al. (2020) employed the Two-Fluid model under laminar flow conditions, proposing a bubble dispersion force to study flow induced by gas bubble generation. The flow advection in such scenarios has highlighted the potential for designing cost-effective, energy-efficient electrolyzers without membranes for gas separation (Pang et al., 2020). Despite these advancements, existing models often fail to capture flow characteristics in the bulk liquid electrolyte away from the electrodes. The present study aims to address this issue by integrating considerations for current density distribution, turbulence effects, and magnetohydrodynamics (MHD). By systematically examining these factors and their collective influence, this study seeks to provide valuable insights into the mechanisms governing flow induced by electrogenerated gas bubbles in electrolytic cells, thereby facilitating the design and optimization of electrochemical systems for sustainable energy production. 2. Modeling All symbols are listed in “Nomenclature”. 2.1. Governing equations We consider the Eulerian (Two-Fluid) multiphase model to simulate the flow of electrolyte and electrogenerated oxygen and hydrogen. The flow is incompressible, Newtonian, and isothermal. The liquid electrolyte and involved gases are considered as interpenetrating continua. Herein, iis the subscript that refers to liquid (l) or gas (g), including hydrogen (H2) or oxygen (O2) phases.The total sum of the volume fraction of all three phases must remain equal to one. α l+ α H2+ α O2=1.(1) The continuity equation for each phase is solved. ∂α i ρ i ∂ t+ ∇ • ( α i ρ iui) = Si.(2) The source term in Eq. (2) describes the generation of hydrogen near the cathode, Eq. (3), and the generation of oxygen near the anode, Eq. (4), as follows: SH2=nH2MH2‖j‖ FWs .(3) SO2=nO2MO2‖j‖ FWs .(4) Following Charton et al. (2009) and Le Bideau et al. (2020), the source terms are released adjacent to electrodes, considering Wsas the width of the computational cell next to the electrodes. Momentum equation, Eq. (5), for each phase is solved to predict the velocity fields as follows: ∂ ( α i ρ iui) ∂ t+ ∇ • ( α i ρ iui⊗ui) = − α i∇p+ ∇ •( α i[ τ i+ τ Re i])+ α i ρ ig+Fd+Fbd +FL. (5) Reynolds-Averaged Navier-Stokes (RANS) formulation is used, in which flow properties are averaged over time, to calculate mean flow quantities such as velocity and pressure. The drag force is described using Eq. (6), considering Schiller and Naumann correlation (Schiller and Naumann, 1935), Fd= − 3Cd ρ l α g 4dbug−ul(ug−ul).(6) Electrogenerated gas bubbles create turbulent-like fluctuations within the multiphase gas-liquid flow due to the combined actions of generated eddies and interphase drag (Zarghami et al., 2020). While the impacts of these fluctuations are minimal in the vertical direction compared to the mean drag and buoyancy forces, they have a significant influence in the spanwise direction, resulting in the dispersion of bubbles (Zarghami et al., 2020). Thus, the bubble dispersion force, Eq. (7), must be included as follows (Le Bideau et al., 2020): Fbd = − Kg ρ l α g dbug−ul∇ α g.(7) The proportionality of the bubble dispersion force with the gradient of the gas fraction indicates the tendency of gas bubbles to transport from regions of high (gas) concentration to regions of low concentration (Zarghami et al., 2020). The interaction of electric current density and the self-induced magnetic field creates the Lorentz force given by, FL=j×B.(8) Momentum equation, Eq. (5), for each phase also involves stress tensor expressed as follows: τ i= μ i[∇ui+ ∇Tui−2 3∇ • (ui)I].(9) Generally, the flow induced by electrogenerated gas bubbles is considered to remain globally laminar, assuming that including the bubble dispersion force, as described in Eq. (7), is sufficient to characterize turbulent-like structures in the flow (Hreiz et al., 2015b; Le Bideau et al., 2020). The validity of this assumption is scrutinized in the present study by considering one case study in which the turbulence effects on the mean flow are included in the model. For that purpose, the Reynolds stress tensor is included as follows: τ Re i= μ t i[∇ui+ ∇Tui−2 3∇ • (ui)I]−2 3 ρ ikiI.(10) Herein, the mixture turbulence model is employed, which is recommended for situations when phases separate, such as the presence of a plume in the liquid phase (ANSYS Fluent 14.5 User’s Guide, 2012). Thus, the mixture turbulent Reynolds stress is expressed as, τ Re m= μ t m[∇um+ ∇Tum−2 3∇ • (um)I]−2 3 ρ mkI.(11) Considering the significant difference in density between the liquid and gas phases, the mixture turbulence approach indeed characterizes the turbulence in the liquid phase (Zarghami et al., 2020). The required parameters to describe the mixture, including density, ρ m=∑3 i=1 α i ρ i, velocity, um=∑3 i=1 α i ρ iui ∑3 i=1 α i ρ i , and turbulent viscosity, μ t m= ρ mC μ k2 ε , are E. Karimi-Sibaki et al. Chemical Engineering Research and Design 208 (2024) 731–739 734 calculated. The latter is predicted by employing the two transport equations (kε ) for turbulent kinetic energy (k), Eq. (12), and turbulent dissipation rate ( ε ), Eq. (13), for the mixture. ∂ ( ρ mk) ∂ t+ ∇ • ( ρ mkum) = ∇ • ( μ t m σ k ∇k)+Gk,m− ρ m ε .(12) ∂ ( ρ m ε ) ∂ t+ ∇ • ( ρ m ε um) = ∇ • ( μ t m σ ε ∇ ε )+ ε k(C1 ε Gk,m−C2 ε ρ m ε ).(13) The model constants in Eqs. (12)-(13). are commonly recognized and provided (ANSYS Fluent 14.5 User’s Guide, 2012). It must be emphasized that no extra term to describe the turbulent interaction between phases within the mixture is required, as turbulent dispersion is already accounted for in the bubble dispersion force through Eq. (7). Boundary conditions for the flow, as extensively described by Le Bideau et al. (2020), involve no-slip conditions, and pressure conditions at entry and outlet. The total pressure is set to zero at the entry, whereas the static pressure is zero at the outlet (Le Bideau et al., 2020). Special care must be taken to define the boundary conditions for turbulent kinetic energy and dissipation rate. The standard wall function approach (ANSYS Fluent 14.5 User’s Guide, 2012) is employed for all boundaries except at the entry point, where the liquid electrolyte is drawn into the gap. We assume that no turbulence is generated at the entry, where the magnitudes of both turbulent kinetic energy and turbulent dissipation rate are set to zero. In literature (Hreiz et al., 2015b), it is common to assume a uniform current density distribution and, consequently, a uniform generation of gases along electrodes, as described by Eqs. (3)-(4). The validity of this assumption is also examined in the present study. For that purpose, the conservation of electric current density (j= − σ ∇φ) is calculated using Eq. (14), while considering the Laplace equation for electric potential given by, ∇ • ( − σ ∇φ) = 0.(14) The electrical conductivity, which varies spatially, reflecting the heterogeneous nature of the system, is dependent on the phase fractions following Bruggeman’s relationship (Cooksey et al., 2008; Le Bideau et al., 2020): σ = σ 0(1− α H2)1.5(1− α O2)1.5.(15) Boundary conditions for electric potential are assigned. All walls except electrodes are assumed to be insulating, thereby setting the flux of electric potential to zero. Following Dahlkild (2001) to account for generated gas bubbles, the electric current density at electrodes, including anode (+) and cathode (-), is as follows: j= ±j0(1− α g)e[−F(Velec−φc) 2RT ](16) Of note, the apparent current density j0in Eq. (16), differs from the exchange current density (Newman and Balsara, 2019; Karimi-Sibaki et al., 2018; Kharicha et al., 2021), which describes the kinetics of redox reactions. Herein, the applied potential (Velec) and the apparent current density (j0)are arbitrary values (Dahlkild, 2001); they must be chosen to ensure that the total average current density along the electrodes is equivalent to the operational current density of the electrolyzer (e.g. here, 2000 A m −2 ). The electric current density flowing through the system induces a magnetic field, which is calculated through Ampere’s law considering vector potential formulation (Song and Ida, 1991; Karimi-Sibaki et al., 2021; Preis et al., 1991) as follows: ∇ × [1 μ 0 ∇ × A]=j.(17) To ensure a unique solution for the magnetic field (B= ∇ × A), the Coulomb gauge (∇ • A=0) is employed (Preis et al., 1991). Eventually, the calculated current density and magnetic fields are used to compute the Lorentz force through Eq. (8). Considering that the magnetic field is predominantly in the polar direction, the flux of the vector potential is set to zero at all walls except for the entry and outlet, where the value of the vector potential is set to zero. 2.2. Other settings All transport phenomena, including flow, electric potential, and magnetic fields, are calculated using the well-established finite volume method (FVM) to discretize the governing equations (Versteeg and Malalasekera, 2007). The spatial discretization scheme is third-order MUSCL (Versteeg and Malalasekera, 2007); the temporal discretization scheme is first-order implicit (Versteeg and Malalasekera, 2007). Several user-defined functions (UDF), e.g. for bubble dispersion force, electrogenerated gas source, and calculation of electromagnetic field, are implemented in the commercial CFD software, ANSYS FLUENT v.17.1, to accurately model the multiphase system. A very fine mesh involving ca. 10 5 quadrilateral square-shaped mesh elements with the size of 60 μ m is generated. The width of the computational domain is 3 mm, and the length is 120 mm. The computational domain, dimensions, and boundaries are schematically shown in Fig. 1, where, for illustrative purposes, the Y direction is scaled down. The computational domain is configured according to the experiment conducted by Boissonneau and Byrne (Boissonneau and Byrne, 2000), who measured the velocity profiles in the inter-electrode gap marked by lines A-A and B-B in Fig. 1 using laser Doppler velocimetry (LDV). Of note, the flow of electrolyte is not due to forced convection (Boissonneau and Byrne, 2000). The generation and rising of bubbles near the electrodes cause spontaneous movement of the surrounding electrolyte, known as the gas-lift phenomenon (Boissonneau and Byrne, 2000). Consequently, the liquid electrolyte is drawn into the gap between electrodes through “Entry”to ensure mass conservation, as shown in Fig. 1. We validated our modeling results against their experimental measurements, as well as analyses presented by Le Bideau et al. (2020) and Schillings et al. (2015). The parameters used in the calculations, extracted from Refs (Schillings et al., 2015; Le Bideau et al., 2020; Boissonneau and Byrne, 2000), are summarized in Table 1. Table 1 Parameters used in our calculations. Parameter Unit Value μ 0Hm−14 π x10 −7 RJK−1mol−18.314546 gms−29.81 FA smol−196485 TK 298 σ 0Sm−15 ρ lkgm−31040 ρ H2kgm−30.082 ρ O2kgm−31.29 μ lkgs−1m−10.00103 μ H2kgs−1m−18.41 ×10 −6 μ O2kgs−1m−11.92 ×10 −5 db μ m 78 j0(Anode) Am−23770 j0(Cathode) Am−24900 Velec (Anode) V 1.5 Velec(Cathode) V 0 Kg/db(H2)ms−12.5 Kg/db(O2)ms−15 nH20.5 nO20.25 MH2gmol−12 MO2gmol−132 E. Karimi-Sibaki et al. Chemical Engineering Research and Design 208 (2024) 731–739 735 3. Results and discussions We structured our research into four distinct case studies to systematically understand the influences of various factors on the overall results. A summary of the different case studies is provided in Table 2. We started with a uniform current density (Case I) to establish a baseline under idealized conditions. Then, we introduced calculated current density (Case II) to observe the impact of more realistic variations. Next, we included magnetohydrodynamics (MHD) (Case III) to account for the Lorentz force arising from the interaction between current density and the magnetic field. Finally, we integrated turbulence (Case IV) to explore the flow characteristics in the bulk of electrolyte. This sequential approach ensures a clear progression from simplified to more comprehensive models, highlighting the contribution of each factor and their collective influence on the results. 3.1. Case study I Case I investigates a uniform distribution of electric current density along the electrodes. Only drag and bubble dispersion forces are considered, and the flow regime is assumed to be laminar, as described in Table 2. Given that the operating current density of industrial alkaline water electrolysis (AWE) typically falls within the range of 2000–4000 A m −2 (Hu et al., 2022), we conduct our investigations with an electric current density of 2000 A m −2 . This choice allows us to validate our findings against results reported in Refs (Schillings et al., 2015; Le Bideau et al., 2020; Boissonneau and Byrne, 2000)., making it particularly relevant for industrial applications. Case I is indeed identical to the analysis presented by Le Bideau et al. (2020). As shown in Fig. 2, the calculated velocity profiles in the inter-electrode gap at two different heights (see Fig. 1) are compared. The proposed models adequately forecast the flow behavior near the electrodes. However, they are deficient in capturing flow characteristics within the bulk of the liquid electrolyte, away from the electrodes. 3.2. Case study II Case II investigates the scenario in which the drag and bubble dispersion forces are considered, and the flow regime is assumed to be laminar, as described in Table 2. However, Case II involves a nonuniform distribution of electric current density along the electrodes. When considering Eq. (16) to calculate electric current density, the mass sources represented by Eqs. (3)-(4), corresponding to electrogenerated gas bubbles, exhibit a non-uniform distribution along the height of the electrodes. Of note, while electric current density becomes non-uniform, the total average current density along the electrodes must remain equivalent to the operational current density of the electrolyzer (2000 A m −2 ). Field structures associated with electric current density are shown in Fig. 3. Utilizing Bruggeman’s relationship, Eq. (15), we observe a nearly 40 % reduction in the electrical conductivity of the liquid-gas mixture near the electrodes, as illustrated in Fig. 3(a). Within the bulk of the electrolyte, particularly in the inter-electrode region, the electric current density maintains a relatively uniform distribution, as demonstrated in Fig. 3(b). Additionally, a minor amount of electric current flows both above and below the electrodes. Iso-lines of electric potential in black are shown, which are always perpendicular to the streamlines of electric current density. Fig. 3(b) further presents the distributions of electric current density along the height of the electrodes. Notably, the peak current density is situated upstream near the edge of the anode/ cathode surface. Subsequently, the current density gradually diminishes along the electrode due to the generation and upward transport of bubbles. The minimum current density is observed at the edges of the electrodes near the downstream walls. In Fig. 4, we compare the calculated velocity profiles in the interelectrode gap at two different heights (see Fig. 1), taking into account the effect of current density distribution. The uniform current represents Case I, whereas the non-uniform current indicates Case II. Although incorporating non-uniform current distribution in the model slightly modifies velocity profiles near the electrodes, it has minimal influence on the velocity in the bulk of the electrolyte. This implies that the assumption of uniform current density along electrodes, which significantly reduces the computational complexity, is practically useful and valid. Our findings are consistent with the analysis provided by Schillings et al. (2015), who suggested that the non-uniformity of current density below 40 % of gas fraction has little effect on the overall results. 3.3. Case study III Assuming a laminar flow regime, Case III accounts for Lorentz forces arising from the interaction between electric current density and the self-induced magnetic field, Eq. (8), as described in Table 2. Our analysis indicates that incorporating the Lorentz force into the model does not affect the overall results. Fig. 5 illustrates the field structures associated with electromagnetic aspects. The magnetic field generated in the polar direction is very weak in the order of the earth’s magnetic field (ca. 0.05 mT). The sign denotes the direction: a positive value indicates that the magnetic field points outward, while a negative value indicates that the magnetic field points inward in the plane, as shown in Fig. 5(a). The distribution of the Lorentz force is illustrated in Fig. 5(b), where the force becomes negligible at the mid-plane location of the electrodes (60 mm), corresponding to the zero magnetic field region. Below the mid-plane, the force acts upward, while above the mid-plane, the force acts downward. For comparison purposes, Fig. 5(c) illustrates the distribution of the bubble dispersion force, as described by Eq. (7). This force exhibits significant strength near the electrodes and downstream walls due to its dependency on the gradient of gas fractions. Notably, the Lorentz force is approximately three orders of magnitude smaller than the bubble dispersion force. Thus, it is safe to assume that the Lorentz force arising from the selfinduced magnetic field is negligible for industrial AWE, which typically operates within the range of 2000–4000 A m −2 as the operating current density. However, it is essential to note that dismissing MHD effects in AWE applications may not be entirely accurate. Magnetized electrodes or an external magnetic field can be deliberately introduced in AWE setups, enhancing the Lorentz force near the electrodes, improving bubble desorption and reducing electrode coverage, as discussed in several studies (Baczyzmalski et al., 2015; Lin et al., 2012;Angulo et al., 2020;Matsushima et al., 2013). Therefore, the self-induced magnetic field/Lorentz force may be negligible under typical operating conditions, yet the intentional use of external magnetic fields can have significant practical implications for the efficiency and effectiveness of AWE systems. 3.4. Case study IV Case IV incorporates the turbulent regime for the gas electrolytic flow in the present study, as described in Table 2. The liquid velocity vector overlaid on the contour of velocity magnitude is shown in Fig. 6 (a). Near the surface of electrodes, gases form a plume, and they are transported upward due to buoyancy, leading to an acceleration of downstream flow. As anticipated, there is an increase in gas fraction (corresponding to a decrease in liquid fraction) along the height of the electrodes, as illustrated in Fig. 6(b). The turbulent viscosity ratio, Fig. 6 (c), and turbulent intensity, Fig. 6(d), are also demonstrated. While the Table 2 Summary of case studies. Case study Electric current distribution Forces Flow regime IUniform Fd,Fbd Laminar II Calculated Fd,Fbd Laminar III Calculated Fd,Fbd,FLLaminar IV Calculated Fd,Fbd,FLTurbulent E. Karimi-Sibaki et al. Chemical Engineering Research and Design 208 (2024) 731–739 736 values predicted for turbulent viscosity ratio and turbulent intensity may be relatively small in magnitude in the inter-electrode region, they are regarded as non-negligible due to their potential impact on the overall behavior of the system. Fig. 7 illustrates the comparison of calculated velocity profiles in the inter-electrode gap at two different heights, considering the effect of turbulence. Incorporating turbulence effects, albeit relatively minor, into the model noticeably enhances the accuracy of the predicted velocity profile. As previously discussed, extensive research has been conducted on the flow induced by electrogenerated gas bubbles between two parallel plates (Hreiz et al., 2015b; Schillings et al., 2015; Le Bideau et al., 2020). Fig. 2. Comparison of predicted velocity profiles for Case I, based on our modeling results and literature data, at positions (a) A-A (Y=60 mm) and (b) B-B (Y=75 mm). Fig. 3. (a) Ratio of electrical conductivity of mixture to the liquid electrolyte, (b) electric current density field, iso-lines of labeled normalized electric potential in black, and plots of calculated current density along the cathode and anode. E. Karimi-Sibaki et al. Chemical Engineering Research and Design 208 (2024) 731–739 737 While the proposed models, accounting for the laminar regime, successfully predict flow behavior near the electrodes, they exhibit limitations in capturing characteristics of the bulk liquid electrolyte away from the electrodes. This study highlights the importance of incorporating the turbulent regime for modeling flow induced by electrogenerated gas bubbles in scenarios without forced convection. 4. Summary Electrolytic gas generation observed, for instance, in alkaline water electrolysis (AWE) to produce hydrogen has been extensively studied over past decades. Understanding the flow induced by electrogenerated gas bubbles between two parallel plates without forced convection is of Fig. 4. Comparison of predicted velocity profiles for Case II, based on our modeling results and literature data, at positions (a) A-A (Y=60 mm) and (b) B-B (Y=75 mm). Fig. 5. (a) Magnetic field, (b) Lorentz force, (c) bubble dispersion force. Equisized arrows indicate directions. E. Karimi-Sibaki et al. Chemical Engineering Research and Design 208 (2024) 731–739 738 great importance when designing cost-effective and energy-efficient (membrane-less) electrolyzers. The process is a complex multiphase flow influenced by numerous factors. Prior models assuming a laminar regime accurately predicted flow near electrodes but failed to predict bulk liquid electrolyte flow away from them. In this study, by employing a three-phase Eulerian model, we examined the influences of electric current density distribution, turbulence effects, and the interaction between flow and the magnetic field known as magnetohydrodynamics (MHD) on the behavior of the process. Through a systematic examination of these factors and their collective influence on the predicted velocity field, the following conclusions are drawn. Incorporating nonuniform current distribution into the model slightly alters electrodeadjacent velocity profiles but minimally affects bulk electrolyte velocity. Hence, assuming an average (uniform) current density along Fig. 6. (a) Velocity field, (b) volume fraction of liquid phase, (c) turbulent viscosity ratio, (d) turbulent intensity. Equisized arrows indicate directions. Fig. 7. Comparison of predicted velocity profiles for Case IV, based on our modeling resultsand literature data, at positions (a) A-A (Y=60 mm) and (b) B-B (Y=75 mm). E. Karimi-Sibaki et al. Chemical Engineering Research and Design 208 (2024) 731–739 739 electrodes is viable. The Lorentz force, arising from flow interaction with the (self-induced) magnetic field, is negligible compared to other forces like bubble dispersion. Therefore, MHD effects can be disregarded unless an external magnetic field is introduced. Incorporating turbulence, being relatively minor in magnitude but substantial between electrodes, significantly improves the predicted velocity profile. Modeling results are validated against experimental measurements. Declaration of Competing Interest We have no conflicts of interest to disclose. Acknowledgements The authors acknowledge financial support from the Austrian Federal Ministry of Economy, Family and Youth and the National Foundation for Research, Technology and Development within the framework of the Christian-Doppler Laboratory for Metallurgical Applications of Magnetohydrodynamics. References Abdelouahed, L., Hreiz, R., Poncin, S., Valentin, G., Lapicque, F., 2014b. Hydrodynamics of gas bubbles in the gap of lantern blade electrodes without forced flow of electrolyte: Experiments and CFD modelling. Chem. Eng. Sci. 111, 255–265. https:// doi.org/10.1016/j.ces.2014.01.028. Abdelouahed, L., Valentin, G., Poncin, S., Lapicque, F., 2014a. Current density distribution and gas volume fraction in the gap of lantern blade electrodes. Chem. Eng. Res. Des. 92, 559–570. https://doi.org/10.1016/j.cherd.2013.10.003. Aldas, K., 2004. Application of a two-phase flow model for hydrogen evolution in an electrochemical cell. Appl. Math. Comput. 154, 507–519. https://doi.org/10.1016/ S0096-3003(03)00731-8. Aldas, K., Pehlivanoglu, N., Mat, M., 2008. Numerical and experimental investigation of two-phase flow in an electrochemical cell. Int. J. Hydrog. Energy 33, 3668–3675. https://doi.org/10.1016/j.ijhydene.2008.04.047. Alexiadis, A., Dudukovic, M.P., Ramachandran, P., Cornell, A., Wanngård, J., Bokkers, A., 2011a. On the electrode boundary conditions in the simulation of two phase flow in electrochemical cells. Int. J. Hydrog. Energy 36, 8557–8559. https:// doi.org/10.1016/j.ijhydene.2011.04.149. Alexiadis, A., Dudukovic, M.P., Ramachandran, P., Cornell, A., Wanngård, J., Bokkers, A., 2011b. Liquid–gas flow patterns in a narrow electrochemical channel. Chem. Eng. Sci. 66, 2252–2260. https://doi.org/10.1016/j.ces.2011.02.046. Alexiadis, A., Dudukovic, M.P., Ramachandran, P., Cornell, A., 2012a. On the stability of the flow in multi-channel electrochemical systems. J. Appl. Electrochem 42, 679–687. https://doi.org/10.1007/s10800-012-0426-0. Alexiadis, A., Dudukovic, M.P., Ramachandran, P., Cornell, A., Wanngård, J., Bokkers, A., 2012b. Transition to pseudo-turbulence in a narrow gas-evolving channel. Theor. Comput. Fluid Dyn. 26, 551–564. https://doi.org/10.1007/s00162011-0248-4. Angulo, A., Van Der Linde, P., Gardeniers, H., Modestino, M., Fern´ andez Rivas, D., 2020. Influence of bubbles on the energy conversion efficiency of electrochemical reactors. Joule 4, 555–579. https://doi.org/10.1016/j.joule.2020.01.005. ANSYS Fluent 14.5 User’s Guide, ANSYS Fluent 14.5 User’s Guide, (2012). Baczyzmalski, D., Weier, T., K¨ ahler, C.J., Cierpka, C., 2015. Near-wall measurements of the bubbleand Lorentz-force-driven convection at gas-evolving electrodes. Exp. Fluids 56, 162. https://doi.org/10.1007/s00348-015-2029-0. Boissonneau, P., Byrne, P., 2000. An experimental investigation of bubble-induced free convection in a small electrochemical cell. J. Appl. Electrochem. 30, 767–775. https://doi.org/10.1023/A:1004034807331. S. Charton, P. Rivalier, D. Ode, J. Morandini, J.P. Caire, Hydrogen production by the Westinghouse cycle: modelling and optimization of the two-phase electrolysis cell, in: Bologna, Italy, 2009: pp. 11–22. https://doi.org/10.2495/ECOR090021. Chi, J., Yu, H., 2018. Water electrolysis based on renewable energy for hydrogen production. Chin. J. Catal. 39, 390–394. https://doi.org/10.1016/S1872-2067(17) 62949-8. Cooksey, M.A., Taylor, M.P., Chen, J.J.J., 2008. Resistance due to gas bubbles in aluminum reduction cells. JOM 60, 51–57. https://doi.org/10.1007/s11837-0080019-x. Dahlkild, A.A., 2001. Modelling the two-phase flow and current distribution along a vertical gas-evolving electrode. J. Fluid Mech. 428, 249–272. https://doi.org/ 10.1017/S0022112000002639. Daoudi, C., Bounahmidi, T., 2024. Overview of alkaline water electrolysis modeling. Int. J. Hydrog. Energy 49, 646–667. https://doi.org/10.1016/j.ijhydene.2023.08.345. El-Askary, W.A., Sakr, I.M., Ibrahim, K.A., Balabel, A., 2015. Hydrodynamics characteristics of hydrogen evolution process through electrolysis: Numerical and experimental studies. Energy 90, 722–737. https://doi.org/10.1016/j. energy.2015.07.108. Endr˝ odi, B., Simic, N., Wildlock, M., Cornell, A., 2017. A review of chromium(VI) use in chlorate electrolysis: Functions, challenges and suggested alternatives. Electrochim. Acta 234, 108–122. https://doi.org/10.1016/j.electacta.2017.02.150. Hreiz, R., Abdelouahed, L., Fünfschilling, D., Lapicque, F., 2015a. Electrogenerated bubbles induced convection in narrow vertical cells: a review. Chem. Eng. Res. Des. 100, 268–281. https://doi.org/10.1016/j.cherd.2015.05.035. Hreiz, R., Abdelouahed, L., Fünfschilling, D., Lapicque, F., 2015b. Electrogenerated bubbles induced convection in narrow vertical cells: PIV measurements and Euler–Lagrange CFD simulation. Chem. Eng. Sci. 134, 138–152. https://doi.org/ 10.1016/j.ces.2015.04.041. Hu, S., Guo, B., Ding, S., Yang, F., Dang, J., Liu, B., Gu, J., Ma, J., Ouyang, M., 2022. A comprehensive review of alkaline water electrolysis mathematical modeling. Appl. Energy 327, 120099. https://doi.org/10.1016/j.apenergy.2022.120099. Karimi-Sibaki, E., Kharicha, A., Abdi, M., Vakhrushev, A., Wu, M., Ludwig, A., Bohacek, J., 2021. A Numerical Study on the Influence of an Axial Magnetic Field (AMF) on Vacuum Arc Remelting (VAR) Process. Met. Mater. Trans. B 52, 3354–3362. https://doi.org/10.1007/s11663-021-02264-w. Karimi-Sibaki, E., Kharicha, A., Wu, M., Ludwig, A., Bohacek, J., 2018. Contribution of an electro-vortex flow to primary, secondary, and tertiary electric current distribution in an electrolyte. J. Electrochem. Soc. 165, E604–E615. https://doi.org/ 10.1149/2.1201811jes. Karlsson, R.K.B., Cornell, A., 2016. Selectivity between Oxygen and Chlorine Evolution in the Chlor-Alkali and Chlorate Processes. Chem. Rev. 116, 2982–3028. https://doi. org/10.1021/acs.chemrev.5b00389. Kharicha, A., Karimi-Sibaki, E., Vakhrushev, A., Wu, M., Ludwig, A., Bohacek, J., 2021. Hydrodynamically enhanced electrochemical mass transfer on the surface of an electrically conductive droplet. Heat. Mass Transf. 57, 1697–1705. https://doi.org/ 10.1007/s00231-021-03071-4. Kharicha, A., Karimi-Sibaki, E., Wu, M., Ludwig, A., Bohacek, J., 2018. Review on modeling and simulation of electroslag remelting. Steel Res. Int. 89, 1700100 https://doi.org/10.1002/srin.201700100. Le Bideau, D., Mandin, P., Benbouzid, M., Kim, M., Sellier, M., Ganci, F., Inguanta, R., 2020. Eulerian two-fluid model of alkaline water electrolysis for hydrogen production. Energies 13, 3394. https://doi.org/10.3390/en13133394. Lin, M.-Y., Hourng, L.-W., Kuo, C.-W., 2012. The effect of magnetic force on hydrogen production efficiency in water electrolysis. Int. J. Hydrog. Energy 37, 1311–1320. https://doi.org/10.1016/j.ijhydene.2011.10.024. Mandin, P., Hamburger, J., Bessou, S., Picard, G., 2005. Modelling and calculation of the current density distribution evolution at vertical gas-evolving electrodes. Electrochim. Acta 51, 1140–1156. https://doi.org/10.1016/j.electacta.2005.06.007. Mat, M., 2004. A two-phase flow model for hydrogen evolution in an electrochemical cell. Int. J. Hydrog. Energy 29, 1015–1023. https://doi.org/10.1016/j. ijhydene.2003.11.007. Matsushima, H., Iida, T., Fukunaka, Y., 2013. Gas bubble evolution on transparent electrode during water electrolysis in a magnetic field. Electrochim. Acta 100, 261–264. https://doi.org/10.1016/j.electacta.2012.05.082. J.S. Newman, N.P. Balsara, Electrochemical systems, Fourth edition, Wiley, Hoboken, NJ, 2019. Nierhaus, T., Van Parys, H., Dehaeck, S., Van Beeck, J., Deconinck, H., Deconinck, J., Hubin, A., 2009. Simulation of the Two-Phase Flow Hydrodynamics in an IRDE Reactor. J. Electrochem. Soc. 156, P139. https://doi.org/10.1149/1.3155423. Pang, X., Davis, J.T., Harvey, A.D., Esposito, D.V., 2020. Framework for evaluating the performance limits of membraneless electrolyzers. Energy Environ. Sci. 13, 3663–3678. https://doi.org/10.1039/D0EE02268C. Preis, K., Bardi, I., Biro, O., Magele, C., Renhart, W., Richter, K.R., Vrisk, G., 1991. Numerical analysis of 3D magnetostatic fields. IEEE Trans. Magn. 27, 3798–3803. https://doi.org/10.1109/20.104929. Rivera, F.F., P´ erez, T., Casta˜ neda, L.F., Nava, J.L., 2021. Mathematical modeling and simulation of electrochemical reactors: A critical review. Chem. Eng. Sci. 239, 116622 https://doi.org/10.1016/j.ces.2021.116622. Rodríguez, J., Amores, E., 2020. CFD modeling and experimental validation of an alkaline water electrolysis cell for hydrogen production. Processes 8, 1634. https:// doi.org/10.3390/pr8121634. Schiller, L., Naumann, A., 1935. A drag coefficient correlation. Z. Des. Ver. Dtsch. Ing. 77, 318–320. Schillings, J., Doche, O., Deseure, J., 2015. Modeling of electrochemically generated bubbly flow under buoyancy-driven and forced convection. Int. J. Heat. Mass Transf. 85, 292–299. https://doi.org/10.1016/j.ijheatmasstransfer.2015.01.121. Song, H., Ida, N., 1991. An eddy current constraint formulation for 3D electromagnetic field calculations. IEEE Trans. Magn. 27, 4012–4015. https://doi.org/10.1109/ 20.104981. Sun, M., Li, B., Liu, Z., 2023. Experimental and numerical analysis on electrolytic bubble group dynamics in the aluminum electrolysis process using the slotted anode. Powder Technol. 428, 118854 https://doi.org/10.1016/j.powtec.2023.118854. H.K. Versteeg, W. Malalasekera, An introduction to computational fluid dynamics: the finite volume method, 2nd ed, Pearson Education Ltd, Harlow, England; New York, 2007. Wedin, R., Dahlkild, A.A., 2001. On the transport of small bubbles under developing channel flow in a buoyant gas-evolving electrochemical cell. Ind. Eng. Chem. Res. 40, 5228–5233. https://doi.org/10.1021/ie001073u. Zarghami, A., Deen, N.G., Vreman, A.W., 2020. CFD modeling of multiphase flow in an alkaline water electrolyzer. Chem. Eng. Sci. 227, 115926 https://doi.org/10.1016/j. ces.2020.115926. E. Karimi-Sibaki et al.