A Numerical Simulation of Evolution Processes and Entropy Generation for Optimal Architecture of an Electrochemical Reaction-Diffusion System: Comparison of Two Optimization Strategies

In this study, a mathematical optimization process is employed to seek the optimal architecture for a porous electrochemical reactor. Such an engineered structure obtained by TO is shown to be capable of enhancing the system performance substantially. This enhancement is obtained by improving the utilization of the constituents materials that is achieved thanks to a heterogeneous layout. This research work introduced a general framework for optimization of constituents placement in an electrochemical reactor. Hence, similar procedure might be used for real applications, such as catalyst layer of a FC. By revisiting the definition of performance enhancement in an electrochemical system, the present work compares two different strategies in optimizing electrochemical systems.

While previous optimization studies (e.g. ^13,14,22 ) predominantly concentrated on the comparative analysis of designs pre- and post-optimization, this research takes a deeper dive into the physical aspects driving design improvements. To achieve this, the entire course of optimization is thoroughly explored, rather than just looking at the starting and ending designs. Furthermore, this study introduces two distinct electrode optimization strategies and examines their physicochemical differences. To the best of the authors' knowledge, such distinctions have not been observed in prior research. As previously highlighted, the full-scale commercialization of electrochemical energy devices hinges on advancing their performance and cost efficiency. In this context, the application of topological optimization to electrode structures may yield highly efficient designs. Furthermore, a more comprehensive understanding of the optimization process through NET analysis could pave the way for uncovering the upper limits of performance achievable with current materials. This insight could also shed light on areas ripe for further improvement through research and development efforts.

The governing equations presented here are based on a volume-averaged approach. The generic differntial form of any balance equation (e.g. mass, energy, and charge) on a control volume is expressed as:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}={\rm{\nabla }}\cdot {\bf{f}}+S\end{eqnarray} \tag{ 2 }$$

where ρ is any quantity of interest, such as mass, energy, and charge. Moreover, f and S represent the flux and source terms of the relevant quantity. In the present study, it is assumed that the distribution of the reactant concentration throughout the porous reactor is controlled by the Fick's second law of diffusion ^44–46 and electrochemical reaction rate. Therefore, conservation of species in this ERD system reads as: ⁴⁴

$$\begin{eqnarray}&&\displaystyle \frac{\partial C}{\partial t}={\rm{\nabla }}\cdot ({D}^{\mathrm{eff}}{\rm{\nabla }}C)+{R}_{{\rm{C}}}\end{eqnarray} \tag{ 3 }$$

in which C is reactant species concentration (mol m⁻³), D^eff is effective diffusivity (m² s⁻¹), and R_C is the reaction source term (mol m⁻³ s⁻¹). This source term is defined as ${R}_{{\rm{C}}}=-\tfrac{{a}^{\mathrm{eff}}{j}_{\mathrm{src}}}{F}$ for an oxidation reaction and ${R}_{{\rm{C}}}=\tfrac{{a}^{\mathrm{eff}}{j}_{\mathrm{src}}}{F}$ for a reduction reaction. Also, j_src is current density source (A m⁻²), a^eff is effective active specific surface area (m² m⁻³), and F is Faraday constant (96485 C mol⁻¹). Under a steady-state condition, the term on the left-hand side of Eq. 3 is zero. Therefore, in the calculations of the present study, the species accumulation term is disregarded. However, since the unsteady form of this equation is used for derivation of entropy production rate in the subsequent section, it is shown in this form. The conservation of electric charges in electrolyte and solid phases follows Ohm's law and is expressed as follows for an electroneutral system: ⁴⁷

$$\begin{eqnarray}&&\displaystyle \frac{\partial {q}^{+}}{\partial t}={\rm{\nabla }}\cdot ({\sigma }_{{\rm{l}}}^{\mathrm{eff}}{\rm{\nabla }}{\phi }_{{\rm{l}}})+{a}^{\mathrm{eff}}{j}_{\mathrm{src}}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {q}^{-}}{\partial t}={\rm{\nabla }}\cdot ({\sigma }_{{\rm{s}}}^{\mathrm{eff}}{\rm{\nabla }}{\phi }_{{\rm{s}}})-{a}^{\mathrm{eff}}{j}_{\mathrm{src}}\end{eqnarray} \tag{ 5 }$$

where q⁺ and q⁻ are positive and negative charge densities (C m⁻³), ${\sigma }_{{\rm{l}}}^{\mathrm{eff}}$ and ${\sigma }_{{\rm{s}}}^{\mathrm{eff}}$ are effective electrolyte and solid phases charge conductivity (S m⁻¹), and ϕ_l and ϕ_s are electrolyte and solid phases potentials (V), respectively. It should be emphasized that the left-hand side of Eqs. 4 and 5 represent accumulation rate of electric charge in the system. In electrochemical reactors (e.g. electrode of batteries or FCs), no charge is accumulated and therefore, this term is usually not shown in the literature. Moreover, since this study assumes a steady-state condition, these terms are zero (similar to what exists in the literature). However, as aforementioned, the unsteady form of governing equations is used for derivation of local entropy generation rate. The current density source, j_src, might be expressed by well-known Butler-Volmer equation as follows. ⁴⁷

$$\begin{eqnarray}&&{j}_{\mathrm{src}}={i}_{0}\left(\displaystyle \frac{C}{{C}^{\mathrm{ref}}}\right)\left\{\exp \left[\displaystyle \frac{{\alpha }_{{\rm{a}}}F\eta }{{RT}}\right]-\exp \left[\displaystyle \frac{-{\alpha }_{{\rm{c}}}F\eta }{{RT}}\right]\right\}\end{eqnarray} \tag{ 6 }$$

In Eq. 6, i₀, C^ref, α_a, α_c, η, R, and T are exchange current density (A m⁻²), reference concentration (mol m⁻³), anodic charge transfer coefficient, cathodic charge transfer coefficient, activation overpotential (V), gas constant (8.3145 J mol⁻¹ K⁻¹), and temperature (K). The correction term $\left(\tfrac{C}{{C}^{\mathrm{ref}}}\right)$ is considered to include the effect of mass transport limitations. In this study, one of the terms on the right-hand-side of Eq. 6 is neglected by assuming a largely-polarized condition. The activation overpotential is measured based on the local electrolyte and solid phases potentials as follows: ⁴⁷

$$\begin{eqnarray}&&\eta ={\phi }_{{\rm{s}}}-{\phi }_{{\rm{l}}}-{E}_{\mathrm{eq}}\end{eqnarray} \tag{ 7 }$$

where E_eq is equilibrium potential (V). As indicated in Eq. 6, in such an ERD system, current density and overpotential are correlated through two rate mechanisms, including: (1) charge transfer rate and (2) mass transport rate. The former mechanism dictates how fast the reactant species is consumed (or likewise, the speed at which the products are produced). On the other hand, the mass transport rate determines the pace of reactant supply and products discharge from the reaction site. The overall process kinetics depends on the slowest mechanism. In electrochemical applications, sluggish mass transport is a serious challenge when the device is working at a high current density and causes a huge drop in output power. Therefore, not only sufficient supply of reactant species but also adequate discharge of products is vital to improve the performance. While thanks to previous experimental studies, the bulk magnitude of mass diffusivity and electric charges conductivity are known for many substances, estimating the effective values of these transport properties in porous media is a long-running dispute among researchers. ^48–50 The Bruggeman equation is an enduring correlation relating the effective properties to the material volume fraction through a power-law. Despite numerous arguments around its validity, ^51,52 the simplicity of the Bruggeman equation led to frequent utilization and recognition of this model in various disciplines, including electrochemistry applications. ^53–57 Therefore, in this study, the effective properties are computed based on the Bruggeman model with a different exponent from that of the original equation as follows:

$$\begin{eqnarray}\begin{array}{rcl}{D}^{\mathrm{eff}} & = & {D}^{0}{\varepsilon }_{{\rm{v}}}^{{\beta }_{1}};\quad {\sigma }_{{\rm{l}}}^{\mathrm{eff}}={\sigma }_{{\rm{l}}}^{0}{\varepsilon }_{{\rm{l}}}^{{\beta }_{2}};\\ {\sigma }_{{\rm{s}}}^{\mathrm{eff}} & = & {\sigma }_{{\rm{s}}}^{0}{\varepsilon }_{{\rm{s}}}^{{\beta }_{3}};\quad {a}^{\mathrm{eff}}={a}^{0}{\varepsilon }_{{\rm{s}}}^{{\beta }_{4}}\end{array}\end{eqnarray} \tag{ 8 }$$

in which D⁰, ${\sigma }_{{\rm{l}}}^{0}$, ${\sigma }_{{\rm{s}}}^{0}$, and a⁰ are bulk mass diffusivity (m² s⁻¹), bulk electrolyte ionic conductivity (S m⁻¹), bulk solid phase electric conductivity (S m⁻¹), and bulk active specific surface area (m² m⁻³), respectively. These bulk values are corrected according to the volume fractions of void (ε_v), electrolyte phase (ε_l), and solid phase (ε_s) together with the relevant given penalty exponents, β₁ to β₄.

The system of partial differential equations (PDEs) given by Eqs. 3 to 5 might be solved in the calculation domain, Ω, if sufficient B.C.s are provided on the boundary of the domain, ∂Ω. While the proposed model could work with any set of sufficient B.C.s, two of them are our interest in this study. The chosen B.C.s, according to Fig. 1, are given by:

$$\begin{eqnarray}\begin{array}{rcl}C & = & {C}^{\mathrm{bnd}}\quad \mathrm{on}\,{{\rm{\Gamma }}}_{1}\quad \mathrm{and}\quad -{D}^{\mathrm{eff}}{\rm{\nabla }}C\cdot \hat{{\bf{n}}}=0\quad \mathrm{on}\,\partial {\rm{\Omega }}\setminus {{\rm{\Gamma }}}_{1}\\ {\phi }_{{\rm{l}}} & = & {\phi }_{{\rm{l}}}^{\mathrm{bnd}}\quad \mathrm{on}\,{{\rm{\Gamma }}}_{3}\quad \mathrm{and}\quad -{\sigma }_{{\rm{l}}}^{\mathrm{eff}}{\rm{\nabla }}{\phi }_{{\rm{l}}}\cdot \hat{{\bf{n}}}=0\quad \mathrm{on}\,\partial {\rm{\Omega }}\setminus {{\rm{\Gamma }}}_{3}\\ {\phi }_{{\rm{s}}} & = & {\phi }_{{\rm{s}}}^{\mathrm{bnd}}\quad \mathrm{on}\,{{\rm{\Gamma }}}_{1}\quad \mathrm{and}\quad -{\sigma }_{{\rm{s}}}^{\mathrm{eff}}{\rm{\nabla }}{\phi }_{{\rm{s}}}\cdot \hat{{\bf{n}}}=0\quad \mathrm{on}\,\partial {\rm{\Omega }}\setminus {{\rm{\Gamma }}}_{1}\end{array}\end{eqnarray} \tag{ 9 }$$

and

$$\begin{eqnarray}&&\begin{array}{l}C={C}^{\mathrm{bnd}}\quad \mathrm{on}\,{{\rm{\Gamma }}}_{1}\quad \mathrm{and}\quad -{D}^{\mathrm{eff}}{\rm{\nabla }}C\cdot \hat{{\bf{n}}}=0\quad \mathrm{on}\,\partial {\rm{\Omega }}\setminus {{\rm{\Gamma }}}_{1}\\ \quad -\,{\sigma }_{{\rm{l}}}^{\mathrm{eff}}{\rm{\nabla }}{\phi }_{{\rm{l}}}\cdot \hat{{\bf{n}}}={j}^{\mathrm{bnd}}\quad \mathrm{on}\,{{\rm{\Gamma }}}_{3}\quad \mathrm{and}\quad -{\sigma }_{{\rm{l}}}^{\mathrm{eff}}{\rm{\nabla }}{\phi }_{{\rm{l}}}\cdot \hat{{\bf{n}}}=0\quad \mathrm{on}\,\partial {\rm{\Omega }}\setminus {{\rm{\Gamma }}}_{3}\\ \quad {\phi }_{{\rm{s}}}={\phi }_{{\rm{s}}}^{\mathrm{bnd}}\quad \mathrm{on}\,{{\rm{\Gamma }}}_{1}\quad \mathrm{and}\quad -{\sigma }_{{\rm{s}}}^{\mathrm{eff}}{\rm{\nabla }}{\phi }_{{\rm{s}}}\cdot \hat{{\bf{n}}}=0\quad \mathrm{on}\,\partial {\rm{\Omega }}\setminus {{\rm{\Gamma }}}_{1}\end{array}\end{eqnarray} \tag{ 10 }$$

where C^bnd, ${\phi }_{{\rm{l}}}^{\mathrm{bnd}}$, ${\phi }_{{\rm{s}}}^{\mathrm{bnd}}$, and j^bnd are constant concentration (mol m⁻³), electrolyte potential (V), solid potential (V), and ionic current density (A m⁻²) on boundary, respectively. Also, $\hat{{\bf{n}}}$ denotes the outward unit normal vector. In the first set specified by Eq. 9, a combination of Dirichlet and homogeneous Neumann B.C.s are introduced for all variables. The values of all variables (concentration, electrolyte potential, and solid potential) are given on one of the system boundaries. This means that the system total overpotential (i.e., the difference of solid and electrolyte phases potentials on the boundaries that the current density is flowing) is given. Based on these settings, the total electric current that passes the system boundaries is computed by solving the system of PDEs. This current represents the total electrochemical reaction rate that is taking place within the ERD system. In Eq. 10, the B.C.s for reactant concentration and solid phase potential are kept the same as in previous conditions. In contrast, the ionic current density passing Γ₃ is prescribed. Accordingly, the total electrochemical reaction rate, that is happening the entire porous reactor, is known. As a result, in such kind of setup, the system overpotential would be determined as a solution of governing equations.

The reason behind these choices for B.C.s are explained as follows. As stated previously, one of the aims of this study is to elucidate the best material distribution (structure) in the porous reactor that enhances the system performance. Hence, it is essential to clearly interpret the meaning of "performance". In electrochemical systems, performance enhancement is achieved either by: (1) maximization of total reaction rate (or likewise, the current that is flowing through the system boundaries) at a given overpotential level or, (2) minimization of system overpotential (which represents the energy loss) for a given total reaction rate. Although the combination of both these strategies might be a third alternative for performance improvement, it is out of the scope of this study. Based on the given definitions of "performance" in an ERD system, the optimization is formulated as both maximization and minimization problems. The B.C.s given by Eqs. 9 and 10 are chosen in a way to fulfill the aforementioned optimization goals. Using the first set of B.C.s, the solution of governing equations (i.e., total current density flowing the system boundaries), might be used as the objective function for the maximization problem. On the other hand, by employment of B.C.s given in Eq. 10, the system overpotential could be used as the objective function for minimization formulation. The details of maximization and minimization problems are comprehensively discussed in Optimization algorithm section. To distinguish the two formulations described above, hereinafter, the governing equations (Eqs. 3 to 5) together with B.C.s of Eq. 9 are called "overpotential-based ERD system" and the combination of the governing equations and B.C.s of Eq. 10 are called "current-based ERD system". The names are assigned based on the electrolyte phase condition on boundary Γ₃ that is given as the model input. Thus, "overpotential-based ERD system" is associated with the maximization optimization and "current-based ERD system" is associated with the minimization optimization. It is also noteworthy that by a proper choice of values for B.C.s given in Eqs. 9 and 10, the electrochemical reaction might be switched between oxidation and reduction. For instance, a positive value for j^bnd in B.C.s of Eq. 10, provides an oxidation reaction. A negative value will result in a reduction reaction.

A system is at an equilibrium state if all macroscopic flows, such as heat or mass fluxes, are zero within its boundaries. Equilibrium thermodynamics, as part of contemporary thermodynamics, successfully handles the systems at an equilibrium state as well as those involving low equilibrium, relatively slow, quasiequilibrium processes. Non-equilibrium thermodynamics (NET) extends conventional equilibrium thermodynamics to characterize the systems which are not in an equilibrium state. ^29,30 The concept has been widely used in the literature to analyze the entropy generation in various engineering applications and thermofluid systems. ^58–64 Since the ERD system in this study is in a non-equilibrium condition from a global standpoint, the concept of NET is recruited to develop a model for evaluation of local and global entropy generation rate. A lower entropy production is only favorable if the output of the system (in ERD system, the total reaction rate) is kept at a certain amount. If a system is working under a finite time/size thermodynamics, doing more (i.e. a higher electrochemical reaction rate in the case of ERD system) will inevitably result in a higher amount of entropy production. Therefore, maximization of current in the "overpotential-based ERD system" is expected to increase the total amount of entropy generation rate. However, at a same level of reaction rate (i.e. happening in "current-based ERD system"), an optimized design solution will minimize the system losses and consequently reduce the total entropy generation of the system. The difference between the two optimization approaches, including maximization and minimization formulations, might be recognized in terms of entropy generation. This contrast is highlighted in Results and discussion section where the results of both strategies are discussed. Development of entropy generation model starts with the extended form of Gibbs equation, which is a fundamental thermodynamics equation and reads as follows: ^29,65,66

$$\begin{eqnarray}&&{dU}={TdS}-{pdV}+\mu {dN}+{\phi }_{{\rm{l}}}{{dQ}}^{+}+{\phi }_{{\rm{s}}}{{dQ}}^{-}\end{eqnarray} \tag{ 11 }$$

where U, S, V, N, Q⁺ and Q⁻ are internal energy (J), entropy (J K⁻¹), volume (m³), the number of moles (mol), positive electric charge (C), and negative electric charge (C) respectively. Also, T, p, and μ are temperature (K), pressure (Pa), and chemical potential (J mol⁻¹), respectively. It is assumed that the reactant substance is an ideal gas, and the system is working under a constant temperature (isothermal) condition. Since the internal energy of an ideal gas is a function of temperature and the volume of the porous reactor is not changing, dU and pdV are zero. By dividing Eq. 11 by the volume, it is simplified into:

$$\begin{eqnarray}&&{Tds}=-\mu {dC}-{\phi }_{{\rm{l}}}{{dq}}^{+}-{\phi }_{{\rm{s}}}{{dq}}^{-}\end{eqnarray} \tag{ 12 }$$

where s = S/V is the entropy density (J k⁻¹ m⁻³). Given the fact that the system does not have any convective flow, the time derivative of the entropy is given by:

$$\begin{eqnarray}&&\displaystyle \frac{\partial s}{\partial t}=-\displaystyle \frac{1}{T}\left(\mu \displaystyle \frac{\partial C}{\partial t}+{\phi }_{{\rm{l}}}\displaystyle \frac{\partial {q}^{+}}{\partial t}+{\phi }_{{\rm{s}}}\displaystyle \frac{\partial {q}^{-}}{\partial t}\right)\end{eqnarray} \tag{ 13 }$$

By substituting the partial derivatives on the right-hand-side of Eq. 13 with Eqs. 3 to 5, the rate of entropy accumulation is as follows:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial s}{\partial t} & = & -{\rm{\nabla }}\cdot \left[\displaystyle \frac{\mu {D}^{\mathrm{eff}}{\rm{\nabla }}C+{\phi }_{{\rm{l}}}{\sigma }_{{\rm{l}}}^{\mathrm{eff}}{\rm{\nabla }}{\phi }_{{\rm{l}}}+{\phi }_{{\rm{s}}}{\sigma }_{{\rm{s}}}^{\mathrm{eff}}{\rm{\nabla }}{\phi }_{{\rm{s}}}}{T}\right]\\ & & +\displaystyle \frac{1}{T}\left(\Space{0ex}{3.2ex}{0ex}{D}^{\mathrm{eff}}{\rm{\nabla }}C\cdot {\rm{\nabla }}\mu \right.\\ & & +{\sigma }_{{\rm{l}}}^{\mathrm{eff}}{\rm{\nabla }}{\phi }_{{\rm{l}}}\cdot {\rm{\nabla }}{\phi }_{{\rm{l}}}+{\sigma }_{{\rm{s}}}^{\mathrm{eff}}{\rm{\nabla }}{\phi }_{{\rm{s}}}\cdot {\rm{\nabla }}{\phi }_{{\rm{s}}}\\ & & \pm \left.\displaystyle \frac{\mu {a}^{\mathrm{eff}}{j}_{\mathrm{src}}}{F}-{\phi }_{{\rm{l}}}{a}^{\mathrm{eff}}{j}_{\mathrm{src}}+{\phi }_{{\rm{s}}}{a}^{\mathrm{eff}}{j}_{\mathrm{src}}\right)\end{array}\end{eqnarray} \tag{ 14 }$$

in which plus and minus signs of term $\pm \tfrac{\mu {a}^{\mathrm{eff}}{j}_{\mathrm{src}}}{F}$ are related to oxidation and reduction reactions, respectively. By comparing the general form of entropy balance $\left(\tfrac{\partial s}{\partial t}=-{\rm{\nabla }}\cdot {{\bf{J}}}_{{\rm{s}}}+\sigma \right)$ ^29,30 with Eq. 14, the entropy flux (W K⁻¹ m⁻²) and entropy generation rate per unit volume (W K⁻¹ m⁻³) are respectively given by:

$$\begin{eqnarray}&&{{\bf{J}}}_{{\rm{s}}}=\displaystyle \frac{\mu {D}^{\mathrm{eff}}{\rm{\nabla }}C+{\phi }_{{\rm{l}}}{\sigma }_{{\rm{l}}}^{\mathrm{eff}}{\rm{\nabla }}{\phi }_{{\rm{l}}}+{\phi }_{{\rm{s}}}{\sigma }_{{\rm{s}}}^{\mathrm{eff}}{\rm{\nabla }}{\phi }_{{\rm{s}}}}{T}\end{eqnarray} \tag{ 15 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}\sigma & = & \displaystyle \frac{1}{T}\Space{0ex}{3ex}{0ex}({D}^{\mathrm{eff}}{\rm{\nabla }}C\cdot {\rm{\nabla }}\mu +{\sigma }_{{\rm{l}}}^{\mathrm{eff}}{\rm{\nabla }}{\phi }_{{\rm{l}}}\cdot {\rm{\nabla }}{\phi }_{{\rm{l}}}+{\sigma }_{{\rm{s}}}^{\mathrm{eff}}{\rm{\nabla }}{\phi }_{{\rm{s}}}\cdot {\rm{\nabla }}{\phi }_{{\rm{s}}}\\ & & \pm \displaystyle \frac{\mu {a}^{\mathrm{eff}}{j}_{\mathrm{src}}}{F}-{\phi }_{{\rm{l}}}{a}^{\mathrm{eff}}{j}_{\mathrm{src}}+{\phi }_{{\rm{s}}}{a}^{\mathrm{eff}}{j}_{\mathrm{src}}\Space{0ex}{1ex}{0ex})\end{array}\end{eqnarray} \tag{ 16 }$$

The first three terms on right-hand-side of entropy generation rate equation, σ, are related to the entropy produced by the transport processes. The other three terms correspond to the source terms in the governing equations (Eqs. 3 to 5). To quantify the individual contribution of each process to the entropy generation rate, the terms are split into four contributions as follows:

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{\mathrm{mt}} & = & \displaystyle \frac{1}{T}{D}^{\mathrm{eff}}{\rm{\nabla }}C\cdot {\rm{\nabla }}\mu \\ {\sigma }_{\mathrm{it}} & = & \displaystyle \frac{1}{T}{\sigma }_{{\rm{l}}}^{\mathrm{eff}}{\rm{\nabla }}{\phi }_{{\rm{l}}}\cdot {\rm{\nabla }}{\phi }_{{\rm{l}}}\\ {\sigma }_{\mathrm{et}} & = & \displaystyle \frac{1}{T}{\sigma }_{{\rm{s}}}^{\mathrm{eff}}{\rm{\nabla }}{\phi }_{{\rm{s}}}\cdot {\rm{\nabla }}{\phi }_{{\rm{s}}}\\ {\sigma }_{\mathrm{rx}} & = & \displaystyle \frac{1}{T}\left(\pm \displaystyle \frac{\mu {a}^{\mathrm{eff}}{j}_{\mathrm{src}}}{F}-{\phi }_{{\rm{l}}}{a}^{\mathrm{eff}}{j}_{\mathrm{src}}+{\phi }_{{\rm{s}}}{a}^{\mathrm{eff}}{j}_{\mathrm{src}}\right)\end{array}\end{eqnarray} \tag{ 17 }$$

where σ_mt, σ_it, σ_et, and σ_rx correspond to the entropy production by mass diffusion, ion transport, electron transport, and reaction (W K⁻¹ m⁻³), respectively. Moreover, to evaluate the total amount of entropy production throughout the entire system (W K⁻¹ m⁻¹), the local entropy generation rate, σ, should be integrated as:

$$\begin{eqnarray}&&P={\int }_{{\rm{\Omega }}}\sigma d{\rm{\Omega }}\end{eqnarray} \tag{ 18 }$$

Likewise, the total entropy production is also broken into four parts, including P_mt = ∫_Ω σ_mt dΩ, P_it = ∫_Ω σ_it dΩ, P_et = ∫_Ω σ_et dΩ, and P_rx = ∫_Ω σ_rx dΩ, each of which is associated with the integral of corresponding local term over the entire domain. The chemical potential of an ideal gas at a constant temperature is computed by: ²⁹

$$\begin{eqnarray}&&\mu ={\mu }^{\mathrm{ref}}+{RT}\mathrm{ln}\left(\displaystyle \frac{C}{{C}^{\mathrm{ref}}}\right)\end{eqnarray} \tag{ 19 }$$

in which μ^ref and R are chemical potential at the reference point (J mol⁻¹) and universal gas constant (8.314 J K⁻¹ mol⁻¹), respectively. For a given electrochemical reaction, the reference chemical potential correlated to the equilibrium potential by the following equation. ⁶⁷

$$\begin{eqnarray}&&{E}_{\mathrm{eq}}=\displaystyle \frac{{\mu }^{\mathrm{ref}}}{F}\end{eqnarray} \tag{ 20 }$$

To generalize the applicability of this study, the governing equations are converted into a dimensionless form. Non-dimensionalization might reduce the number of model parameters, expand the applicability of the obtained results, increase the stability and accuracy of the simulated results, and make the comparison between various scenarios easier. ^68,69 The system of governing equations, as well as the B.C.s, are non-dimensionalized by substitution of proper dimensionless variables to remove the units. The dimensionless form of Eqs. 3 to 5 at a steady-state condition reads as:

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{* }\cdot \left({\theta }_{1}{\varepsilon }_{{\rm{v}}}^{{\beta }_{1}}{{\rm{\nabla }}}^{* }{C}^{* }\right)\mp {\varepsilon }_{{\rm{s}}}^{{\beta }_{4}}{j}_{\mathrm{src}}^{* }\,=\,0\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{* }\cdot \left(\displaystyle \frac{{\theta }_{2}}{1+{\theta }_{3}}{\varepsilon }_{{\rm{l}}}^{{\beta }_{2}}{{\rm{\nabla }}}^{* }{\phi }_{{\rm{l}}}^{* }\right)+{\varepsilon }_{{\rm{s}}}^{{\beta }_{4}}{j}_{\mathrm{src}}^{* }\,=\,0\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{* }\cdot \left(\displaystyle \frac{{\theta }_{2}{\theta }_{3}}{1+{\theta }_{3}}{\varepsilon }_{{\rm{s}}}^{{\beta }_{3}}{{\rm{\nabla }}}^{* }{\phi }_{{\rm{s}}}^{* }\right)-{\varepsilon }_{{\rm{s}}}^{{\beta }_{4}}{j}_{\mathrm{src}}^{* }\,=\,0\end{eqnarray} \tag{ 23 }$$

in which C^*, ${\phi }_{{\rm{l}}}^{* }$, and ${\phi }_{{\rm{s}}}^{* }$ are dimensionless concentration, electrolyte potential, and solid potential, respectively. Moreover, ∇^* is differentiation with respect to dimensionless coordinates x^* and y^*. The dimensionless dependent and independent variables are given by:

$$\begin{eqnarray}\begin{array}{rcl}{x}^{* } & = & \displaystyle \frac{x}{L};\quad {y}^{* }=\displaystyle \frac{y}{L};\quad {C}^{* }=\displaystyle \frac{C}{{C}_{0}};\quad {\mu }^{* }=\displaystyle \frac{\mu }{{RT}}\\ {\phi }_{{\rm{l}}}^{* } & = & \displaystyle \frac{{\phi }_{{\rm{l}}}}{{\phi }_{0}};\quad {\phi }_{{\rm{s}}}^{* }=\displaystyle \frac{{\phi }_{{\rm{s}}}}{{\phi }_{0}};\quad {a}^{* }=\displaystyle \frac{{a}^{\mathrm{eff}}}{{a}_{0}};\quad {j}_{\mathrm{src}}^{* }\,=\,\displaystyle \frac{{j}_{\mathrm{src}}}{{i}_{0}}\end{array}\end{eqnarray} \tag{ 24 }$$

In addition, the dimensionless groups are expressed as:

$$\begin{eqnarray}\begin{array}{rcl}{\theta }_{1} & = & \displaystyle \frac{{D}^{0}}{{L}^{2}{a}^{0}{i}_{0}/({C}_{0}F)}=\displaystyle \frac{{D}^{0}}{{D}^{\mathrm{exchange}}}\\ {\theta }_{2} & = & \displaystyle \frac{({\sigma }_{{\rm{l}}}^{0}+{\sigma }_{{\rm{s}}}^{0})}{{L}^{2}{a}^{0}{i}_{0}/{\phi }_{0}}=\displaystyle \frac{({\sigma }_{{\rm{l}}}^{0}+{\sigma }_{{\rm{s}}}^{0})}{{\sigma }_{\mathrm{ohm}}^{\mathrm{exchange}}}\\ {\theta }_{3} & = & \displaystyle \frac{{\sigma }_{{\rm{s}}}^{0}}{{\sigma }_{{\rm{l}}}^{0}}\end{array}\end{eqnarray} \tag{ 25 }$$

where L, C₀ and ϕ₀ are the characteristic length (m), characteristic concentration (mol m⁻³), and characteristic potential (V) of the system, respectively. Based on these definitions, θ₁ is the dimensionless conductivity. This dimensionless number could be interpreted as the ratio of species intrinsic diffusivity and system exchange diffusivity (D^exchange = L² a⁰ i₀/(C₀ F)). From another perspective, θ₁ is attributed to inverse of Damköhler number (or inverse square Thiele modulus) in chemical engineering literature, ^70–73 which indicates the ratio of the bulk diffusion rate and reaction rate. Moreover, θ₂ represents the ratio of intrinsic conductivity and system exchange ohmic conductivity (${\sigma }_{\mathrm{ohm}}^{\mathrm{exchange}}={L}^{2}{a}^{0}{i}_{0}/{\phi }_{0}$). In electrochemistry literature, dimensionless parameter θ₂ is known as Wagner number, ^47,67 which shows the proportion of kinetic resistance to ohmic resistance. Finally, θ₃ represents the ratio of solid and electrolyte phases conductivity. In many real electrochemical applications, such as polymer electrolyte membrane fuel cell (PEMFC), this ratio is relatively large (θ₃ ≫ 1). While there might be different alternatives for the characteristic parameters of the system, a reasonable choice is given:

$$\begin{eqnarray}&&{C}_{0}={C}^{\mathrm{bnd}};\quad {\phi }_{0}=\displaystyle \frac{{RT}}{F}\end{eqnarray} \tag{ 26 }$$

Based on Eq. 26, the characteristic concentration is assumed to be the constant concentration on the inlet boundary. This way, the concentration B.C. would be easily scaled to unity. Also, characteristic potential is chosen to be the thermal voltage of the system that is correlated to the temperature at which the system is working. The characteristic length is considered to be equal to the length of the domain, as shown in Fig. 1, for the sake of simplicity. Using the dimensionless parameters introduced before, the B.C.s of "overpotential-based ERD system" and "current-based ERD system" are non-dimensionalized as follows:

$$\begin{eqnarray}\begin{array}{rcl}{C}^{* } & = & 1\quad \mathrm{on}\,{{\rm{\Gamma }}}_{1}\quad \mathrm{and}\quad -{D}^{* }{{\rm{\nabla }}}^{* }{C}^{* }\cdot \hat{{\bf{n}}}=0\quad \mathrm{on}\,\partial {\rm{\Omega }}\setminus {{\rm{\Gamma }}}_{1}\\ {\phi }_{{\rm{l}}}^{* } & = & {\phi }_{{\rm{l}}}^{\mathrm{bnd},* }\quad \mathrm{on}\,{{\rm{\Gamma }}}_{3}\quad \mathrm{and}\quad -{\sigma }_{{\rm{l}}}^{* }{{\rm{\nabla }}}^{* }{\phi }_{{\rm{l}}}^{* }\cdot \hat{{\bf{n}}}=0\quad \mathrm{on}\,\partial {\rm{\Omega }}\setminus {{\rm{\Gamma }}}_{3}\\ {\phi }_{{\rm{s}}}^{* } & = & {\phi }_{{\rm{s}}}^{\mathrm{bnd},* }\quad \mathrm{on}\,{{\rm{\Gamma }}}_{1}\quad \mathrm{and}\quad -{\sigma }_{{\rm{s}}}^{* }{{\rm{\nabla }}}^{* }{\phi }_{{\rm{s}}}^{* }\cdot \hat{{\bf{n}}}=0\quad \mathrm{on}\,\partial {\rm{\Omega }}\setminus {{\rm{\Gamma }}}_{1}\end{array}\end{eqnarray} \tag{ 27 }$$

and

$$\begin{eqnarray}&&\begin{array}{l}{C}^{* }=1\quad \mathrm{on}\,{{\rm{\Gamma }}}_{1}\quad \mathrm{and}\quad -{D}^{* }{{\rm{\nabla }}}^{* }{C}^{* }\cdot \hat{{\bf{n}}}=0\quad \mathrm{on}\,\partial {\rm{\Omega }}\setminus {{\rm{\Gamma }}}_{1}\\ \quad -{\sigma }_{{\rm{l}}}^{* }{{\rm{\nabla }}}^{* }{\phi }_{{\rm{l}}}^{* }\cdot \hat{{\bf{n}}}={j}^{\mathrm{bnd},* }\quad \mathrm{on}\,{{\rm{\Gamma }}}_{3}\quad \mathrm{and}\quad -{\sigma }_{{\rm{l}}}^{* }{{\rm{\nabla }}}^{* }{\phi }_{{\rm{l}}}^{* }\cdot \hat{{\bf{n}}}=0\quad \mathrm{on}\,\partial {\rm{\Omega }}\setminus {{\rm{\Gamma }}}_{3}\\ \quad {\phi }_{{\rm{s}}}^{* }={\phi }_{{\rm{s}}}^{\mathrm{bnd},* }\quad \mathrm{on}\,{{\rm{\Gamma }}}_{1}\quad \mathrm{and}\quad -{\sigma }_{{\rm{s}}}^{* }{{\rm{\nabla }}}^{* }{\phi }_{{\rm{s}}}^{* }\cdot \hat{{\bf{n}}}=0\quad \mathrm{on}\,\partial {\rm{\Omega }}\setminus {{\rm{\Gamma }}}_{1}\end{array}\end{eqnarray} \tag{ 28 }$$

in which, the boundary values are converted to a dimensionless form according to:

$$\begin{eqnarray}&&{\phi }_{{\rm{l}}}^{\mathrm{bnd},* }=\displaystyle \frac{{\phi }_{{\rm{l}}}^{\mathrm{bnd}}}{{\phi }_{0}};\quad {\phi }_{{\rm{s}}}^{\mathrm{bnd},* }=\displaystyle \frac{{\phi }_{{\rm{s}}}^{\mathrm{bnd}}}{{\phi }_{0}};\quad {j}^{\mathrm{bnd},* }=\displaystyle \frac{{j}^{\mathrm{bnd}}}{{i}_{0}}\end{eqnarray} \tag{ 29 }$$

Furthermore, the dimensionless entropy generation rate, σ^* might be determined through diving the left-hand-side of Eq. 16 by the exchange entropy generation rate (W K⁻¹m⁻³) and is formulated as:

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }^{* } & = & \displaystyle \frac{\sigma }{{a}^{0}{i}_{0}{\phi }_{0}/T}=\displaystyle \frac{\sigma }{{\sigma }^{\mathrm{exchange}}}={\theta }_{1}{\varepsilon }_{{\rm{v}}}^{{\beta }_{1}}{{\rm{\nabla }}}^{* }{C}^{* }\cdot {{\rm{\nabla }}}^{* }{\mu }^{* }+\displaystyle \frac{{\theta }_{2}}{1+{\theta }_{3}}{\varepsilon }_{{\rm{l}}}^{{\beta }_{2}}{{\rm{\nabla }}}^{* }{\phi }_{{\rm{l}}}^{* }\cdot {{\rm{\nabla }}}^{* }{\phi }_{{\rm{l}}}^{* }\\ & & +\displaystyle \frac{{\theta }_{2}{\theta }_{3}}{1+{\theta }_{3}}{\varepsilon }_{{\rm{s}}}^{{\beta }_{3}}{{\rm{\nabla }}}^{* }{\phi }_{{\rm{s}}}^{* }\cdot {{\rm{\nabla }}}^{* }{\phi }_{{\rm{s}}}^{* }\pm {\varepsilon }_{{\rm{s}}}^{{\beta }_{4}}{\mu }^{* }{j}_{\mathrm{src}}^{* }-{\varepsilon }_{{\rm{s}}}^{{\beta }_{4}}{\phi }_{{\rm{l}}}^{* }{j}_{\mathrm{src}}^{* }+{\varepsilon }_{{\rm{s}}}^{{\beta }_{4}}{\phi }_{{\rm{s}}}^{* }{j}_{\mathrm{src}}^{* }\end{array}\end{eqnarray} \tag{ 30 }$$

in which the exchange entropy generation rate is defined as σ^exchange = a⁰ i₀ ϕ₀/T. The individual contributions, σ_mt, σ_it, σ_et, and σ_rx, could be non-dimensionalized to ${\sigma }_{\mathrm{mt}}^{* }$, ${\sigma }_{\mathrm{it}}^{* }$, ${\sigma }_{\mathrm{et}}^{* }$, and ${\sigma }_{\mathrm{rx}}^{* }$ after dividing by the exchange entropy production. Finally, the total dimensionless entropy generation is computed as:

$$\begin{eqnarray}&&{P}^{* }={\int }_{{{\rm{\Omega }}}^{* }}{\sigma }^{* }d{{\rm{\Omega }}}^{* }\end{eqnarray} \tag{ 31 }$$

where Ω^* represents the non-dimensional domain. The related dimensionless breakdowns, ${P}_{\mathrm{mt}}^{* }$, ${P}_{\mathrm{it}}^{* }$, ${P}_{\mathrm{et}}^{* }$, and ${P}_{\mathrm{rx}}^{* }$, are evaluated by a similar integration.

Since the current density that passes the system boundaries is given in a "current-based ERD system", the total rate of electrochemical reaction is prescribed in this system setting as an input. In this formulation, the goal is to minimize the system total overpotential for a certain reaction rate. Total overpotential indicates the excess required potential to drive a Faradaic reaction at a specific current density. Hence, it is favorable to reduce the total overpotential in an electrochemical system for a given current density. Minimization of total overpotential will lead to maximization of output power for a power generating system and will reduce the input power for a power consuming system. Based on these explanations, the optimization problem is defined as:

$$\begin{eqnarray}&&\begin{array}{l}\mathop{\min }\limits_{{\varepsilon }_{{\rm{l}}},{\varepsilon }_{{\rm{s}}}}\quad {F}_{\mathrm{obj},\ 1}=\left|{\eta }^{\mathrm{tot},* }\right|=\left|\left({\phi }_{{\rm{s}}}^{* }{| }_{{{\rm{\Gamma }}}_{1}}-{\phi }_{{\rm{l}}}^{* }{| }_{{{\rm{\Gamma }}}_{3}}-{E}_{\mathrm{eq}}^{* }\right)\right|\\ {\rm{s}}.{\rm{t}}.\quad \mathrm{Eqs.}\ 21\ \mathrm{to}\ 23\\ \quad {\varepsilon }_{{\rm{l}}}^{\min }\leqslant {\varepsilon }_{{\rm{l}}}({{\bf{x}}}^{* })\leqslant {\varepsilon }_{{\rm{l}}}^{\max }\quad \forall {{\bf{x}}}^{* }\in {{\rm{\Omega }}}^{* }\\ \quad {\varepsilon }_{{\rm{s}}}^{\min }\leqslant {\varepsilon }_{{\rm{s}}}({{\bf{x}}}^{* })\leqslant {\varepsilon }_{{\rm{s}}}^{\max }\quad \forall {{\bf{x}}}^{* }\in {{\rm{\Omega }}}^{* }\\ \quad {\varepsilon }_{{\rm{l}}}({{\bf{x}}}^{* })+{\varepsilon }_{{\rm{s}}}({{\bf{x}}}^{* })+{\varepsilon }_{{\rm{v}}}({{\bf{x}}}^{* })=1\quad \forall {{\bf{x}}}^{* }\in {{\rm{\Omega }}}^{* }\end{array}\end{eqnarray} \tag{ 32 }$$

in which F_{obj, 1} is the objective function and ${E}_{\mathrm{eq}}^{* }=\tfrac{{E}_{\mathrm{eq}}}{{\phi }_{0}}$ is non-dimensional equilibrium potential. Knowing that the total overpotential for a reduction reaction is negative, the absolute value of overpotential is used in the objective function definition. It is noteworthy that in an ideal condition, the reactor potential (i.e., ${\phi }_{{\rm{s}}}^{* }{| }_{{{\rm{\Gamma }}}_{1}}-{\phi }_{{\rm{l}}}^{* }{| }_{{{\rm{\Gamma }}}_{3}}$) should be as close as possible to the thermodynamic equilibrium potential (${E}_{\mathrm{eq}}^{* }$). Volume fractions of electrolyte and solid phases are considered as design variables in this optimization problem. These volume fractions are the system parameters that are required for solving the system of governing equations. Since the summation of all volume fractions should be equal to unity, by knowing two of them, the third one might be calculated automatically. Hence, the porosity (volume fraction of voids) is calculated based on the values of two design variables. Moreover, the optimization problem is constrained by the governing equations of the system. As an additional constraint, it is assumed that the local volume fractions of the electrolyte and solid phases could only alter within a predefined range. These ranges are given by ${\varepsilon }_{{\rm{l}}}^{\min }$, ${\varepsilon }_{{\rm{l}}}^{\max }$, ${\varepsilon }_{{\rm{s}}}^{\min }$, and ${\varepsilon }_{{\rm{s}}}^{\max }$ that are minimum and maximum allowed volume fractions of electrolyte and solid phases, respectively.

In an "overpotential-based ERD system", the total overpotential (${\eta }^{\mathrm{tot}}={\phi }_{{\rm{s}}}{| }_{{{\rm{\Gamma }}}_{1}}-{\phi }_{{\rm{l}}}{| }_{{{\rm{\Gamma }}}_{3}}$) is given as the problem input through the B.C.s and the total electrochemical reaction rate might be calculated by solving the governing equations. Total electrochemical reaction rate is equivalent to the total current that passes the system boundary. Since the system is electroneutral, the total amount of negative and positive electric charges that leaves (or come into) the system is equal. Hence, in this system setting, the aim is to maximize the total electrochemical reaction rate. Increasing the electrochemical reaction rate (or system current) at a given overpotential is equivalent to maximizing the output power of the system (for a power generating system). Improving the active specific surface area and discharge of products from the reaction sites will boost the reaction rate. This could be obtained by increasing the volume fractions of electrolyte and solid phases. However, to keep the reaction rate high, it is necessary to assure sufficient reactant delivery to the reaction sites through pores. By increasing the volume fractions of electrolyte and solid phases, the volume fraction of voids will consequently be reduced. This means that the mass transport resistance will be increased and the delivery of reactant spices will be more complicated. Therefore, there should be a trade-off between reactant consumption and delivery rates. This trade-off leads to an optimization problem that might be addressed by controlling the local composition of the porous reactor. This optimization problem is mathematically indicated in a non-dimensional form as:

$$\begin{eqnarray}&&\begin{array}{l}\mathop{\max }\limits_{{\varepsilon }_{{\rm{l}}},{\varepsilon }_{{\rm{s}}}}\quad {F}_{\mathrm{obj},\ 2}={\displaystyle \int }_{{{\rm{\Omega }}}^{* }}\left|{\varepsilon }_{{\rm{s}}}^{{\beta }_{4}}{j}_{\mathrm{src}}^{* }\right|d{{\rm{\Omega }}}^{* }\\ {\rm{s}}.{\rm{t}}.\quad \mathrm{Eqs.}\ 19\ \mathrm{to}\ 21\\ \quad {\varepsilon }_{{\rm{l}}}^{\min }\leqslant {\varepsilon }_{{\rm{l}}}({{\bf{x}}}^{* })\leqslant {\varepsilon }_{{\rm{l}}}^{\max }\quad \forall {{\bf{x}}}^{* }\in {{\rm{\Omega }}}^{* }\\ \quad {\varepsilon }_{{\rm{s}}}^{\min }\leqslant {\varepsilon }_{{\rm{s}}}({{\bf{x}}}^{* })\leqslant {\varepsilon }_{{\rm{s}}}^{\max }\quad \forall {{\bf{x}}}^{* }\in {{\rm{\Omega }}}^{* }\\ \quad {\varepsilon }_{{\rm{l}}}({{\bf{x}}}^{* })+{\varepsilon }_{{\rm{s}}}({{\bf{x}}}^{* })+{\varepsilon }_{{\rm{v}}}({{\bf{x}}}^{* })=1\quad \forall {{\bf{x}}}^{* }\in {{\rm{\Omega }}}^{* }\end{array}\end{eqnarray} \tag{ 33 }$$

where F_{obj, 2} is the objective function. Given the fact that ${j}_{\mathrm{src}}^{* }$ is negative in case of a reduction reaction, the absolute value of reaction rate is used in this equation to generalize the applicability of the objective function for any kind of electrochemical reaction. As represented in Eq. 33, the goal is to maximize the overall dimensionless rate of electrochemical reaction over the entire domain. The term ${a}^{* }{j}_{\mathrm{src}}^{* }$ is the local dimensionless rate of reaction and by integrating it over the non-dimensional domain, the overall value might be computed. Evidently, this objective function could be expressed in a dimensional form with the corresponding dimensional parameters. However, since this study focuses on a non-dimensional analysis of the ERD system, the objective function is also expressed in a dimensionless form.

In this research work, COMSOL Multiphysics ^® (version 5.6) is used to solve the governing equations and to implement the optimization procedure. The calculation domain, system parameters, governing equations, and B.C.s are specified in the software. In the first step of the optimization process, the design variables are initialized and Eqs. 21 to 23 are solved using a finite element method. By solving these equations, the distribution of concentration, electrolyte potential, and solid potential are obtained. Next, the objective function is evaluated according to F_{obj, 1} or F_{obj, 2}. The choice of objective function depends on the system of interest. Afterward, the sensitivity (or gradient) of objective function is computed with respect to the design variables using an adjoint state method. ⁸⁰ Since the adjoint method is independent of the number of the design variables, it is computationally much cheaper than other approaches, such as the forward method. Therefore, it is very useful for optimization problems which contain numerous decision variables. Next, the design solutions are regularized using a Helmholtz filter ⁸¹ and hyperbolic tangent project. ⁸² These regularization helps to prevent checkerboard pattern, ¹⁶ which is a common challenge in TO. This problem has been addressed in the literature comprehensively and the readers are referred to the papers in this field, such as Refs. 81, 83, 84, for further information. In the next step, the design variables are updated locally in throughout the entire calculation domain. The updating procedure is conducted using a globally convergent method of moving asymptotes (GCMMA) algorithm. ⁸⁵ GCMMA is an efficient algorithm for handling problems with abundant design variables. As a post-processing step, the local and total entropy production rate are also calculated. This step is separate from the optimization process and the entropy generation is assessed in each iteration of optimization to track the changes. The explained process is repeated until the convergence criterion (in this case, maximum number of iteration) is reached. The optimization procedure is summarized as follows:

• Step 1:   
  The calculation domain, governing equations, and B.C.s are defined and the design variables are initialized.
• Step 2:   
  The ERD system is solved based on Eqs. 21 to 23 and the objective function is evaluated using Eq. 32 or 33.
• Step 3:   
  The gradient of objective function with respect to the design variables is evaluated using the adjoint method.
• Step 4:   
  The design solutions are regularized using a Helmholtz filter and a hyperbolic tangent projection to obtain smoothed solutions.
• Step 5:   
  The design variables are updated using the GCMMA algorithm.
• Step 6:   
  The local and global entropy generation rates are evaluated using Eqs. 30 and 31.
• Step 7:   
  If the convergence criteria are met, the iteration is stopped. Otherwise, steps 2 to 6 are repeated.

As explained before, in this problem setting, the objective is to reduce the total overpotential of the system. The total overpotential refers to the amount of energy necessary for an electrochemical reaction to proceed at a given conversion rate (current density). Figures 2 and 3 depict the optimized distribution of constituents under various system settings. The optimal spatial distribution of electrolyte, solid, and void phases are illustrated for two different dimensionless current densities (j^bnd,*). In each figure, the optimized results for four different combinations of dimensionless numbers (θ₁ to θ₃) are shown to reflect the impact of each number on the final optimization outcome. Analyzing and contrasting the results from each system configuration elucidates the impact of diverse transport, rate, and characteristic parameters on the ultimate optimal structure. Under all conditions, the optimization process led to complex tree-root-like shapes. Comparing the final objective function values of each scenario across Figs. 2 and 3 shows that in general, a higher working current density corresponds to a higher total overpotential. For instance, as j^bnd,* is changed from 30 to 5 at the first θs settings (Figs. 2a vs 3a), the non-dimensional total overpotential is altered from 4.64 to 0.589. This is due to the intrinsic relationship between current density and overpotential in an electrochemical system. The higher the working current density of the system, the higher total overpotential. The initial layout used for all the cases reported in these two Figures are the same (uniform ε_l = 0.25 and ε_s = 0.25). However, it is evident that the solid volume fraction in final layouts of Fig. 3 (j^bnd,* = 5) is higher than those of Fig. 2 (j^bnd,* = 30). At a low current density, the activation overpotential is the major contributor to the total overpotential. Therefore, it is favorable to increase the average volume fraction of the solid phase. On the contrary, at a high current density(j^bnd,* = 30), where sluggish mass transport is dominant, larger diffusion channels are preferred to assure a sufficient reactant delivery. As it can be understood from Eq. 25, the system mass transport is controlled by dimensionless number θ₁. Therefore, changing the value of this parameter should affect the optimized distribution of voids. As it can be seen in Fig. 2, decreasing the value of this parameter influences the optimized porosity distribution in the ERD system. However, this impact is more significant at a high current density because of a serious mass transport limitation under this operational condition. The same changes in θ₁ at a low current density cause a slight change in the final void distribution. This is due to the fact that at j^bnd,* = 5, the concentration overpotential is small and a variation of θ₁ could be compensated by a tiny change in the optimized layout (see Figs. 3a and 3b).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In addition, dimensionless number θ₂ dictates the ohmic conductivity of the system. Hence, diminishing this parameter is associated with new layouts (Figs. 2c and 3c) which provide a better pathway for conduction of electric charges. Finally, decreasing parameter θ₃ boosts the ionic conductivity of the system in comparison to the electric conductivity. In many real applications, the electrochemical devices suffer from a high ionic resistance, but not electric one. Therefore, decreasing θ₃ (by keeping θ₂ constant) improves the balance between electron and cation transport in the ERD system. Consequently, a lower total overpotential could be achieved after optimization with a smaller θ₃. For j^bnd,* = 5, this means a less requirement for the electrolyte phase (see Fig. 3d). The optimized layouts shown in Figs. 2 and 3 reduce the overpotential of the ERD system between 39.6% and 64.2% compared to the initial uniform configuration (before optimization). The improvements owe to the heterogeneous structures obtained from TO. These tree-root-like structures facilitate diffusion of reactant species and electric charges while augmenting the reaction rate throughout the system. For instance, in the case of mass diffusion, all optimized structures involve some primary diffusion channels together with secondary channels. The primary channels are extended throughout the system in x-direction, which facilitates the delivery of the reactant material to the regions far from the inlet boundary (∂Ω₁). Evidently, the system of interest in this study is symmetrical with respect to the y-axis (see Fig. 1). Therefore, it is necessary to direct the reactant in x-direction, especially at a high j^bnd,*, which mass transport limitation is more severe. Additionally, the secondary channels connected to the primary ones help a better dissemination of the reactant over the entire porous reactor. On a similar note, the electrolyte and solid phases are also arranged in a way that give a greater access to electric charge currents. Since the left boundary is isolated to the transport of ions, the electrolyte phase is more concentrated on the right half of the reactor. It is notable that the solid phase not only transfers the produced electrons to the outside but also provides the necessary surface area for the reaction. Therefore, both these roles are reflected in the optimal layout. The heterogeneous distribution of materials within the ERD system strikes a balance between multiple transport and rate processes, consequently leading to enhanced overall performance.

To further analyze the optimization process, the scenario with θ₁ = 5, θ₂ = 50, and θ₃ = 10 is assessed in details. The convergence history of the optimization process is plotted in Fig. 4a along with the system layout at two intermediate iterations (iteration number 20 and 70). As expected, during the optimization, the objective function monotonically decreased. A sharp decline is observed at the beginning, which is due to the higher sensitivity (or gradient) of the objective function in these iterations. As the optimization proceeds, the change rate slows down. The objective function (total overpotential) reaches a constant value of 4.64 after around 125 iterations, where the changes in objective function drops under 0.01% afterwards. This assures a confident convergence. Figure 4b demonstrates the current density-overpotential (I-V) relationship of the system before and after optimization. Indeed, an I-V curve shows the performance of the system for a wide range of current densities. To demonstrate how the results of optimization at various current densities are different, the I-V curves are plotted for two structures obtained from optimization at j^bnd,* equal to 30 (point A on the plot) and 5 (point A on the plot). From the first optimization, point A is shifted to ${\rm{A}}^{\prime}$ and the structure shown in Fig. 2a is obtained. Then, this optimal layout is used to plot the I-V curve over a range of 5 ≤ j^bnd,* ≤ 30. A similar procedure has been conducted for optimization at point B. In both cases, the optimal I-V curve is vertically shifted downward (vertical optimization). However, the final I-V characteristics of the two optimal structures are different. Evidently, at medium and high current densities (j^bnd,* ≥ 15), I-V curve of ${\rm{A}}^{\prime}$ is superior. However, in low current density regions (j^bnd,* ≤ 5), the I-V characteristics of ${\rm{B}}^{\prime}$ surpasses that of the other optimization. The reason for this difference is when the optimization is performed at a high current density, the optimal layout comes with large diffusion channels (see Fig. 2a). At j^bnd,* ≥ 15, the concentration depletion is the prevailing mechanism. Hence, larger diffusion channels (higher average porosity) is beneficial to make a compromise between the competing processes. On the other hand, at j^bnd,* ≤ 10, sufficient reactant is supplied to the reaction site. Hence, a higher volume fraction of the solid phase is needed to increase the effective electrochemical surface area and consequently decrease the system total overpotential. What stands out in this figure is that when optimizing an electrochemical system at a specific working point, the optimal result might not be equally good over the entire I-V curve. Therefore, the optimization point should be adjusted according to the application and purposes. The average effective transport coefficients of the system before and after optimization are compared in Table II. The transport properties after optimization are divided by the relevant exchange properties (D^exchange and ${\sigma }_{\mathrm{ohm}}^{\mathrm{exchange}}$) to obtain the dimensionless properties. As it can be seen in this table, all transport properties are considerably enhanced after optimization. The only exception is average diffusivity after optimization at j^bnd,* = 5. Since the overall reaction rate is low at j^bnd,* = 5, the reactant delivery is sufficient even at a relatively lower effective diffusivity. Therefore, TO favors to reduce the porosity (effective diffusivity) and instead increases the volume fraction of solid phase to augment the electrochemical surface area.

Zoom In Zoom Out Reset image size

Table II. Average values of dimensionless transport properties and volume fractions before and after optimization for θ₁ = 5, θ₂ = 50, and θ₃ = 10; the values in parentheses shows the changes with respect to the before optimization case.

Optimization	D*, ave	${\sigma }_{{\rm{l}}}^{* ,\ \mathrm{ave}}$	${\sigma }_{{\rm{s}}}^{* ,\ \mathrm{ave}}$	${\varepsilon }_{{\rm{l}}}^{\mathrm{ave}}$	${\varepsilon }_{{\rm{s}}}^{\mathrm{ave}}$	${\varepsilon }_{{\rm{v}}}^{\mathrm{ave}}$
Before	1.25 (−)	0.284 (−)	2.84 (−)	0.25	0.25	0.5
After - jbnd,* = 30	1.62 (+29.6%)	0.679 (+139%)	5.05 (+77.8%)	0.264	0.313	0.423
After - jbnd,* = 5	0.766 (−38.7%)	0.537 (+89.1%)	9.78 (+244%)	0.274	0.446	0.280

To clarify the benefits of a heterogeneous structure, the optimal layout (obtained from optimization of θ₁ = 5, θ₂ = 50, and θ₃ = 10 at j^bnd,* = 30) is compared to two uniform distributions. As aforementioned, the optimization process started with a homogeneous material distribution in which volume fractions of both electrolyte and solid phases were 0.25. In this case, the non-dimensional total overpotential is 9.18. After optimization, this value is declined to 4.64, which shows a 49.5% reduction in system losses. Such a substantial improvement proves the advantage of the topologically-optimized layout over the uniform distribution. In the final optimum layout, the average volume fraction of various phases are ${\varepsilon }_{{\rm{l}}}^{\mathrm{ave}}=0.264$, ${\varepsilon }_{{\rm{s}}}^{\mathrm{ave}}=0.313$, and ${\varepsilon }_{{\rm{v}}}^{\mathrm{ave}}=0.423$ as shown in Table II. Considering these average values with a uniform distribution gives a total overpotential of $\left|{\eta }^{\mathrm{tot},* }\right|=9.83$, which is considerably higher than 4.64. In fact, it is even higher than that of before optimization. This reveals that the placement of constituents is as important as their overall average amounts. In other words, TO adjusts the overall volume fraction of constituents as well as their distribution in a manner to reduce the total overpotential.

A closer inspection of system evolution is required to obtain a deeper grasp about how TO lead to a better design solution. To do so, some of the system parameters are traced during optimization. Since the system is symmetrical, the parameters of interest are projected on x-axis to reduce their distribution to 1D for the sake of simplicity. The projected dimensionless concentration (C^*) and current source ($\left|{a}^{* }{j}_{\mathrm{src}}^{* }\right|$) along with entropy production evolution are set out in Fig. 5. What stands out in Fig. 5a is during the optimization process, the reactant concentration is elevated in the system. The concentration at x^* = 1 is increased by two orders of magnitude after the optimization. This is attributed to a better porosity distribution that assist delivery of reactant to the regions far from the system inlet. The emerged diffusion channels in x-direction prevent significant concentration overpotential. This is especially essential for high operational current densities where mass transport limitation is dominant. What is striking about these results is that the increase in concentration occurs in a condition where the average porosity in the optimal structure (0.423) is lower than the initial one (0.5). Checking the general trend of projected current source, that is depicted in Fig. 5b, shows that the optimization algorithm favors design solutions with a more uniform current source distribution. This means that to lessen the total overpotential, it is necessary to scatter the electrochemical reaction all over the ERD system. From a global standpoint, spreading the reaction throughout the design domain assists to reduce the entropy that is generated in the system. Figure 5a illustrates the changing progress of dimensionless total entropy generation, P^*, in the course of optimization along with local distribution at some selected iteration numbers. As it can be confirmed by both local and global analysis, the total entropy generation rate is decreasing as the optimization proceeds. This observation is in full agreement with the principle of Entropy Generation Minimization (EGM). ⁴³ Hence, for a minimization problem set-up, the total overpotential could be substituted with the total entropy generation as an objective function. It will be later discussed that this replacement does not hold for a maximization formulation. Since the entropy generation rate is reduced after optimization, it can be concluded that the second law efficiency of the system is increased. To understand how the entropy production reduced by TO, the distribution of projected entropy is plotted in Fig. 5b for some iterations. It can be seen that the entropy curve at iteration zero is at a higher level than other iterations. After some early changes in the system layout, the entropy curve shifted downward considerably. This huge drop in entropy generation happened thanks to an improvement in transport processes. In the subsequent steps, the optimization algorithm attempted to further reduce entropy generation rate by providing a more uniform distribution of entropy production throughout the system (see the curves for iterations 35, 50, and 300 in Fig. 5b). A more uniform distribution of entropy production might (but not necessarily) be helpful in reducing the total entropy production and consequently increasing the system efficiency. This principle is known as Equipartition of Entropy Production (EoEP) ^86,87 in the literature. It has been applied to many heat exchanger and membrane systems. ^88–95

Zoom In Zoom Out Reset image size

These behaviors can be observed in the individual entropy contribution by each process. Figure 6 illustrates the breakdown of entropy production by various processes, including mass diffusion, ion transfer by electrolyte phase, electron transfer by solid phase, and electrochemical reaction. The projected distribution of entropy production by each mechanism is shown in this figure. From these breakdowns, one can understand that the entropy curves of transport processes (Figs. 6a to 6c) lay quickly down hinged at one boundary. Especially, the entropy production by ion transport (Fig. 6b) shows a significant drop. This is because of the enhanced pathways formed by TO for transport of reactant, ion, and electron. Meanwhile, the flattening of these entropy curves causes a sharp decrease in total entropy production, as shown in Fig. 5a. The behavior of entropy contribution by the electrochemical reaction is, however, different from that of transport phenomena. According to the curves of Fig. 6d, the entropy production distribution becomes more equipartitioned during the optimization. This is mainly because of a more uniform distribution of the current source that is discussed before. To achieve this flatness, the entropy production in part of the system is reduced, while it is increased in the other part. Therefore, the total amount of entropy production by the reaction does not change very much. This could be confirmed from the column chart demonstrated in Fig. 6e. This chart compares the total entropy production by each process. It can be seen that the contribution by electrochemical reaction stays almost constant during the optimization process. The slight reduction is related to the more equipartitioned distribution of entropy generation, which prevents excess dissipation caused by non-uniformity of thermodynamics driving force in the system. ^86,87 Since the total reaction rate is given by the system current density, most entropy production by the reaction is the inevitable entropy generation that cannot be prevented in a finite time/size context. ^18,77 Based on the data of Fig. 6e, the electron transport has the smallest contribution to the total entropy compared to other transport phenomena. This is because of the high intrinsic electronic conductivity of the system.

Zoom In Zoom Out Reset image size

In this section, the optimization results of the overpotential-based ERD system are briefly discussed. In contrast to the previous strategy, the goal of this optimization is to increase the current density that is passing the system boundary at a given overpotential level. The optimized system layout for an ERD system with θ₁ = 5, θ₂ = 50, and θ₃ = 10 are demonstrated in Fig. 7 for two different overpotentials (${\phi }_{{\rm{l}}}^{\mathrm{bnd},* }=1$ and 9). The structure obtained for ${\phi }_{{\rm{l}}}^{\mathrm{bnd},* }=9$ (Fig. 7a) increased the system current density by 116.7%. A closer look at this optimal layout reveals that it is pretty similar to Fig. 2b. While the exact distributions are different in the two figures, the overall shapes look alike. However, the system settings (θ₁) are different. The reason for this similarity is as follows. The value of θ₁ in the optimized system of Fig. 2b has a small value, meaning that the system is diffusion-controlled (diffusion-limited) by nature. Therefore, the optimization algorithm attempts to compensate for the weak reactant delivery by forming a distributed high porosity region in half of the reactor (far from the inlet boundary). On the other hand, the system in Fig. 7a is not diffusion-limited by nature. However, increasing the current density necessitates a faster reactant delivery to the reaction sites. As stated before, the output power of a power-generating electrochemical system could be increased in two ways, including reduction of total overpotential and increment of current density. The latter comes at the cost of a higher overall electrochemical reaction rate. To increase the reaction rate, it is required to supply more reactant. This explains why the optimized layouts in Fig. 2b and Fig. 7a are similar despite their different settings (different θ₁). The results of optimization at a lower overpotential (${\phi }_{{\rm{l}}}^{\mathrm{bnd},* }=1$) are presented in Fig. 7a. At a lower overpotential level, the activation overpotential prevails. As a result, a higher volume fraction of the solid phase is beneficial to increase the current density.

Zoom In Zoom Out Reset image size

From the convergence diagram shown in Fig. 8a, it is understandable that the objective function reaches a constant maximum value after around 115 iterations. The insets in Fig. 8a illustrate how the distributions of decision variables are changing over the optimization course before fully converging. The I-V curves obtained from the optimized systems for the two cases (${\phi }_{{\rm{l}}}^{\mathrm{bnd},* }=1$ and 9) are plotted in Fig. 8b. Opposed to a minimization problem, this optimization approach shifted the initial I-V curve horizontally to the right (horizontal optimization). Performing an optimization at a high overpotential point (${\phi }_{{\rm{l}}}^{\mathrm{bnd},* }=9$, point C in Fig. 8b) leads to a significant improvement in terms of I-V characteristics for medium and high overpotentials (${\phi }_{{\rm{l}}}^{\mathrm{bnd},* }\geqslant 5$). However, this optimization is mainly focused on the improvement of reactant deficit as a result of increased reaction rate. Since the volume fraction of solid phase is relatively low in the optimized structure, the I-V performance of this structure is comparatively weak at low overpotentials (${\phi }_{{\rm{l}}}^{\mathrm{bnd},* }\leqslant 3$). In addition, optimizing the system at a lower overpotential (${\phi }_{{\rm{l}}}^{\mathrm{bnd},* }=1$) increases the current density by 70.5%. However, it is not as beneficial when the system is working at a high overpotential level.

Zoom In Zoom Out Reset image size

To obtain some insight about maximization approach and identify the differences with the previous formulation, the evolution of projected concentration and convergence of entropy production are investigated for optimization at ${\phi }_{{\rm{l}}}^{\mathrm{bnd},* }=9$. As it can be seen in Fig. 8c, the concentration changes non-monotonically during the optimization process. In the primary iterations, when there is an increase in the reaction rate as a result of an increase in ionic conductivity, the reactant concentration suddenly drops. To compensate for the concentration depletion and raise the reaction rate, TO forms some diffusion channels, which helped the concentration to be almost retrieved to the value before optimization. Increasing the concentration through facilitated mass transport helps the system to work at a higher reaction rate without experiencing any reactant shortage. The behavior of dimensionless entropy generation differs from that of the minimization approach. As shown in Fig. 8d, the total entropy production is increasing during the optimization. Before the optimization, the local distribution of entropy generation was relatively uniform and its quantity was also low. However, the final optimized solution comes with a higher entropy generation rate, especially in areas with a high volume fraction of the solid phase. This behavior might be assigned to the fact that the electrochemical reaction occurs in those areas. In a finite time/size thermodynamics context, any system that is "doing" something is doomed to produce some amount of entropy. Therefore, when the "doing" of a system is increased (in this case the reaction rate is increased), it is expected that total entropy production to rise. As discussed in Refs. 18, 77, this increase is associated with the "inevitable entropy production". Another striking point in Fig. 8d is that despite the objective function shows a monotonic convergence behavior (see Fig. 8a), the convergence curve of total entropy production shows some noises. This means that the maximization formulation may produce some design solutions that are desirable in terms of objective function but are not completely favorable from the perspective of the second law of thermodynamics. Hence, maximization of current density is not equivalent to the maximization entropy production and the objective function cannot be substituted. To rephrase it, a lower entropy production corresponds to a lower dissipation; nonetheless, a higher entropy production does not imply a higher reaction rate. This is the major difference between the two optimization strategies. A vertical optimization attempts to reduce the losses in the system; however, a horizontal optimization aims to increase the conversion capacity of the system. Therefore, the results obtained from the latter approach might not have the highest possible second law efficiency. In summary, the findings of this study proves that introducing a heterogeneous material distribution within an electrode can effectively enhance material utilization, thereby improving overall performance. Thisenhancement is achieved through a balance among competing transport and rate processes. Unlike parametric optimization, which primarily controls design variables globally, TO fine-tunes material distribution locally, thus surpassing conventional optimization methods. Moreover, the results of this study show that a proper choice of optimization approach depends on the application of interest. In principle, both strategies might be used for any application by sufficient considerations. Yet, a vertical optimization may be recommended in power generating systems (e.g. FCs) where there is a need to reduce the energy losses. Moreover, a horizontal optimization may be employed for power consuming systems (e.g. electrolyzers) in which the conversion rate is a major challenge. It is also important to choose the optimization point on I-V curve in accordance with the practical working conditions. Although the optimized design enhances the performance at the chosen optimization point, it might not be beneficial for other working conditions. In the case of an ERD system, the I-V characteristic does not improve equally over a wide operating range.