Optimization of a Membraneless Microﬂuidic Fuel Cell with a Double-Bridge Flow Channel

: In this work, a design optimization study was conducted to improve the performance of a membraneless microﬂuidic fuel cell with a double-bridge cross-section of the ﬂow channel. Governing equations including Navier–Stokes, mass-transport, and Butler–Volmer equations were solved numerically to analyze the electrochemical phenomena and evaluate the performance of the fuel cells. Optimization was performed to maximize the peak power density using a genetic algorithm combined with a surrogate model constructed by radial basis neural network. Two sub-channel widths of the ﬂow channel were selected as design variables for the optimization. As a result, a large increase in the inner channel width and a small decrease in the outer channel width effectively increased the peak power density of the MMFC. The optimal design increased the peak power density by 57.6% compared to the reference design.


Introduction
Current combustion-based energy generation technology is very harmful to the environment and causes global problems such as ozone-layer depletion and the consistent reduction in vegetation cover resulting from climate change [1]. On the other hand, a fuel cell provides an efficient and clean mechanism as an electrochemical device that can directly convert the chemical energy of fuel into electrical energy [1,2].
Fuel cells are used to generate electricity in some applications, such as spacecraft and emergency generators, and also used as power sources for portable electronic components [3]. A research was conducted for 1-100 mW class applications of fuel cells such as distributed microsensors and wireless micro-electro mechanical systems [4], and fuel cells were introduced into small electric vehicles [5]. Fuel cells can supply electricity together with cogeneration [6]. As the use of renewable energy has been pursued, it is also being used in the field of solar heat and solar power [7]. On-chip integration of fuel cells miniaturizes the reactor and maximizes fuel storage, so they are used in applications with limited form factors, such as wireless sensor nodes [8]. A fuel cell is recommended as a high power source in portable applications due to its active system integrating single cells and stacks, and is used as a power source, such as in chargers, LEDs and small fans [9]. Due to the current increase in various diseases, implantable medical devices and wearable devices are continuously being developed, and fuel cells have been applied to overcome the problems of power generation and replacement [10]. Fuel cells are also used in military equipment because of their high power density, durability, and convenient portability [11].
Recently, as the use of portable electronics has increased, the demand for battery power has also increased, but the existing battery technology cannot keep up with this demand [12]. A fuel cell solely works to generate energy and an additional unit is coupled alongside to store the energy. Thus, in fuel cells, both units can be individually customized

MMFC Geometry
The configuration of the channel and double-bridge channel cross-section of the MMFC to be optimized in this study is shown in Figure 1. As shown in Figure 1a, fuel and oxidant enter the double-bridge channel through the inlets and proceed along the flow path. The length (L 0 ) of the microchannel is 10 mm, the half width (W) of the channel is 150 µm, and the height (H) of the channel is 50 µm. The channel half-width (W) is divided into three sub-channel widths. The widths of the outer channel, bridge channel, and inner half channel are indicated by W O , W B , and W I , respectively, as shown in Figure 1b. The three sub-channel widths are same at 50 µm, and the reference width is W r = 50 µm, which is defined to make the sub-widths dimensionless. The bridge height is h = 10 µm. For the optimization, W is fixed at 150 µm with a constraint of W B /W r ≥ 0.2. The design specifications of the reference channel are found in Table 1. Neural Network (RBNN) was used to construct the surrogate model of the objective function from the numerical results of the sample designs, and a Genetic Algorithm (GA) was utilized to find the optimal design based on the surrogate model. The power density of the MMFC was obtained from a numerical simulation solving Navier-Stokes, masstransport, and Butler-Volmer equations.

MMFC Geometry
The configuration of the channel and double-bridge channel cross-section of the MMFC to be optimized in this study is shown in Figure 1. As shown in Figure 1a, fuel and oxidant enter the double-bridge channel through the inlets and proceed along the flow path. The length (L0) of the microchannel is 10 mm, the half width (W) of the channel is 150 µm, and the height (H) of the channel is 50 µm. The channel half-width (W) is divided into three sub-channel widths. The widths of the outer channel, bridge channel, and inner half channel are indicated by WO, WB, and WI, respectively, as shown in Figure 1b. The three sub-channel widths are same at 50 µm, and the reference width is Wr = 50 µm, which is defined to make the sub-widths dimensionless. The bridge height is h = 10 µm. For the optimization, W is fixed at 150 µm with a constraint of WB/Wr ≥ 0.2. The design specifications of the reference channel are found in Table 1.

Numerical Analysis
Due to the low Reynolds number in the MMFC channel, the flow is laminar. To simplify the electrochemical analysis and the numerical model, it was assumed that the temperature and density of the working fluids are constant and the flow is steady. It was also assumed that the depletion region does not have a significant effect on the thickest mixed region at the outlet of the flow path. For the numerical investigation, the commercial computational fluid dynamics (CFD) software COMSOL Multiphysics ® (version 4.3, COMSOL, Inc., Stockholm, Sweden) [31] was used.
Due to the high energy density and convenient storage, liquid reactants have been used in numerous studies on MMFCs [14,[32][33][34]. In this study, formic acid (HCOOH) and potassium permanganate (KMnO 4 ) were used as the fuel and oxidant, respectively. Formic acid as a fuel promotes electron and proton transport within the anode compartment of the fuel cell [35]. Potassium permanganate as an oxidizing agent also improves cathode dynamics and makes it possible to generate higher power densities than air [14].
The electrochemical reaction of the anode is as follows: In addition, the primary permanganate reduction occurring at the cathode is as follows: The permanganate ion is further oxidized to insoluble MnO 2 by the following reaction: The reactants were transported downstream along the channel by convection. Fick's law of diffusion was used to decide the transport of reactants in the MMFC channel [36]. In the numerical model, it was assumed that flows are isothermal, and pressure differences are not large enough. Implementing these assumptions allows the concentration distributions of the fuel and oxidant to be described by Fick's law [37]. Fick's law predicts variation in the concentration of reactant species with time along the channel caused by the diffusion.
where u is velocity vector, D i is the diffusion coefficient, and C i is the local concentration for species i. In the electrode subdomains, a source term is added to the above equation to accommodate the consumption of species by electrochemical reactions. This results in the following equation: where the source term (S i ) is calculated by Butler-Volmer equations. The electrochemical reactions in the MMFC were modelled using the Butler-Volmer equations [36]. Additionally, the reactants consumed by the electrode were calculated using Faraday's law. The concentration of ions (i.e., H + or OH − ) was assumed to be constant in the flow channel [38]. The Butler-Volmer, Navier-Stokes, and convection-diffusion equations were combined to predict the current density for a specified cell voltage. At the channel inlet, constant concentration and cell pressure drop are specified. Additionally, at the channel outlet, constant pressure and mass flux are specified. In addition, no-slip conditions were specified on the wall of the channel.
Previous simulations of MMFCs generally introduced a single-phase flow assumption, and the present computational model also used the single-phase assumption. The formation of gas bubbles (e.g., O 2 and CO 2 ) was neglected, and it was supposed that the bubbles formatted at the anode are fully dissolved in the electrolyte [17,19,[39][40][41][42]. Wang et al. [43] and Shyu et al. [44] studied two-phase numerical models to confirm the effects of bubble formation on the performance of an MMFC. More detailed information on the numerical model can be found in ref. [30].

Design Variables
Selection of proper design variables is important in optimization. Thus, a parametric study must be performed to find the parameters affecting the objective function. When the range of each design variable is narrowed, the accuracy of the surrogate model constructed for the objective function increases for the same number of sample designs. In the previous work [30], a parametric study was conducted to find the important parameters for the design of the double-bridge MMFC among the sub-channel widths, W I , W B and W O with fixed W/W r = 3.0. The results indicated that the power density is relatively insensitive to W B among the three parameters. When W I /W r = 1.7 (W O /W r = 0.3) with fixed W B /W r = 1.0, the peak power density increased by about 18% compared to the reference design. As W I increased, the mixing layer became thinner, reducing ohmic loss. The results of the parametric study [30] with fixed W B are summarized in Table 2 and Figure 2. The overall results of the parametric study [30] indicated that W I is the most critical parameter and W O is the second important parameter in terms of sensitivity of the peak power density, and these two parameters were selected as design variables for the present optimization.   Table 3 shows the design ranges of the two design variables. In the optimization, the following two design constraints must be satisfied:   Table 3 shows the design ranges of the two design variables. In the optimization, the following two design constraints must be satisfied: The first constraint secures the minimum width of the bridge structure, which consists of the electrode and supporting structure. The second constraint indicates that the overall width of the channel remains constant for the optimization.

Optimization Methodology
The present optimization goal is to maximize the peak power density of the double-bridge MMFC. Figure 3 illustrates the present optimization procedure based on surrogate modeling.

Design of Experiment
Latin hypercube sampling (LHS) [45] was employed to select design points (or sample points) in the design space. Twenty design points for the two design variables (WI and WO) were generated using the lhsdesign function in MATLAB [46] (MathWorks, Inc., Natick, MA, USA). The LHS technique is a sampling technique that reduces time by using an m×n simulation matrix and can efficiently analyze variable space at a minimum cost. Here, m is the number of samples to be extracted and n is the number of variables. Additionally, this method is designed to represent all parts of the design space and generate random sample points. After removing sample designs that did not fit the constraints, the remaining nine sample designs were used for evaluating the objective function. Table 4 lists the numerical results of the nine sample designs. The optimization problem can be defined as: Design variable bound : Here, F(x) is the objective function, i.e., the peak power density, x(= {x i }) is a vector of n design variables, and LB i and UB i are the lower and upper limits, respectively.

Design of Experiment
Latin hypercube sampling (LHS) [45] was employed to select design points (or sample points) in the design space. Twenty design points for the two design variables (W I and W O ) were generated using the lhsdesign function in MATLAB [46] (MathWorks, Inc., Natick, MA, USA). The LHS technique is a sampling technique that reduces time by using an m×n simulation matrix and can efficiently analyze variable space at a minimum cost. Here, m is the number of samples to be extracted and n is the number of variables. Additionally, this method is designed to represent all parts of the design space and generate random sample points. After removing sample designs that did not fit the constraints, the remaining nine sample designs were used for evaluating the objective function. Table 4 lists the numerical results of the nine sample designs.  [47] was used to construct the surrogate model. RBNN is an artificial neural network that uses radial basis functions as activation functions, and has many uses, including function approximation, time series prediction, classification, and system control. The RBNN model includes the output layer of linear neurons, as well as the hidden layer of radial neurons. The hidden layer uses radial primitives to transform the input space into an intermediate space. The output of hidden layers is passed to a combiner to create targets. In MATLAB, a RBNN model was created using newrb function with the input being the obtained data of the design of the experiment.

Searching Algorithms
GA [48] was selected as the searching algorithm in the present work. GA finds the maximum of a function by mimicking biological evolutionary processes. This method randomly selects chromosomes before starting and sets the initial population size. Subsequent computations allow one chromosome to adapt in a competitive way. Chromosomes selected after the start of computation are called parental generation, and those selected after parental generation are called child generation. By genetic search operators (selection, crossover, and mutation), chromosomes that are superior to those of the previous generation are extracted. In MATLAB [46], single-objective optimization using GA can be performed using a function.

Grid Convergence Index Analysis and Validation of Numerical Results
To evaluate the grid-dependency of the numerical results, the Grid Convergence Index (GCI) analysis for the MMFC current was performed using Richardson extrapolation [49] in the previous work [30]. Table 5 summarizes the GCI analysis results. To adjust the number of nodes, the grid segmentation index was set to 1.3. The current of MMFC shows a tendency to converge gradually as the number of grid nodes increased, as can be seen in Figure 4. Extrapolated relative error (when using N 1 ) e 21 ext is 0.131%, and GCI 21 f ine is 0.163%, indicating that the numerical uncertainty is sufficiently small. Consequently, N 1 was selected as the final grid system. Additional calculations were performed based on the structure of this grid. A tetrahedral mesh was generated in the designated computational domain, and the mesh was generated using the physics-controlled mesh option built in COMSOL. A mesh created in the boundary layer was composed of fine elements. Table 5. Grid convergence index (GCI) analysis for the MMFC current [30].

Parameters Values
Number of cells  To validate the numerical results, the predicted performance curves for the singlebridge channel were compared with the experimental data measured by Montesinos et al. [17] by Oh et al. [30] using the same numerical model as in the present work, as shown in Figure 5. Compared with the experimental data, Figure 5a shows the lowest relative error of 0.65% at 0.7 V. Overall, good agreement can be seen between the numerical and experimental data. In the power density versus current density curve (Figure 5b), the numerical results somewhat differ from the previous experimental results at lower cell voltages. However, the MMFC functioned at a nominal cell voltage between 0.6 and 0.8 V. It was confirmed that it was in good agreement with the experimental results. The OCP (Open Circuit Potential) in this study is 1.4 V. However, the typical OCP is in a range of 0.6~1.2 V. The difference is due to the high oxidation and reduction potential of KMnO4. Other studies have also confirmed cases with higher OCP. Wang et al. [50] summarized MMFCs investigated with various reactants, and highlighted the impact of fuel and oxidant combination on OCP. Furthermore, it also varies with the nature of electrolyte. Various studies have highlighted that a dual-electrolyte configuration (acidic-alkaline) results in higher OCP. Liu et al. [51] employed the acid-alkaline dual-electrolyte together with KMnO4 as an oxidant. The cell OCV was as high as 1.83 V. To validate the numerical results, the predicted performance curves for the single-bridge channel were compared with the experimental data measured by Montesinos et al. [17] by Oh et al. [30] using the same numerical model as in the present work, as shown in Figure 5. Compared with the experimental data, Figure 5a shows the lowest relative error of 0.65% at 0.7 V. Overall, good agreement can be seen between the numerical and experimental data. In the power density versus current density curve (Figure 5b), the numerical results somewhat differ from the previous experimental results at lower cell voltages. However, the MMFC functioned at a nominal cell voltage between 0.6 and 0.8 V. It was confirmed that it was in good agreement with the experimental results. The OCP (Open Circuit Potential) in this study is 1.4 V. However, the typical OCP is in a range of 0.6~1.2 V. The difference is due to the high oxidation and reduction potential of KMnO 4 . Other studies have also confirmed cases with higher OCP. Wang et al. [50] summarized MMFCs investigated with various reactants, and highlighted the impact of fuel and oxidant combination on OCP. Furthermore, it also varies with the nature of electrolyte. Various studies have highlighted Energies 2022, 15, 973 9 of 13 that a dual-electrolyte configuration (acidic-alkaline) results in higher OCP. Liu et al. [51] employed the acid-alkaline dual-electrolyte together with KMnO 4 as an oxidant. The cell OCV was as high as 1.83 V. results somewhat differ from the previous experimental results at lower cell voltages. However, the MMFC functioned at a nominal cell voltage between 0.6 and 0.8 V. It was confirmed that it was in good agreement with the experimental results. The OCP (Open Circuit Potential) in this study is 1.4 V. However, the typical OCP is in a range of 0.6~1.2 V. The difference is due to the high oxidation and reduction potential of KMnO4. Other studies have also confirmed cases with higher OCP. Wang et al. [50] summarized MMFCs investigated with various reactants, and highlighted the impact of fuel and oxidant combination on OCP. Furthermore, it also varies with the nature of electrolyte. Various studies have highlighted that a dual-electrolyte configuration (acidic-alkaline) results in higher OCP. Liu et al. [51] employed the acid-alkaline dual-electrolyte together with KMnO4 as an oxidant. The cell OCV was as high as 1.83 V.

Optimization Results
From the RBNN model constructed for the objective function, GA found the optimal design of the double-bridged MMFC at WI/Wr = 1.922 (WI = 96.1 µm) and WO/Wr = 0.815 (WO = 40.7 µm) as shown in Table 6. The WB value of the optimal design (WB/Wr = 0.263) is located near the design limit, indicating the minimum thickness of the bridge structure. The optimum value of WI is more than twice as large as the optimum value of WO, which indicates that the effect of WI on the power density is far larger than that of WO. The predicted optimized peak power density is 49.85 mW/cm 2 while the CFD result is 50.24 mW/cm 2 , which yielded a relative error of only 0.77%. This small error proves that the use of the RBNN model and GA is suitable for the investigated problem. In terms of improvement from the reference design, the optimal design increased the peak power density significantly by 57.6% from that of the reference design. Figure 6 shows the peak power density curves for the reference and optimal designs of the double-bridge MMFC, which confirms the achievement of the optimization.

Optimization Results
From the RBNN model constructed for the objective function, GA found the optimal design of the double-bridged MMFC at W I /W r = 1.922 (W I = 96.1 µm) and W O /W r = 0.815 (W O = 40.7 µm) as shown in Table 6. The W B value of the optimal design (W B /W r = 0.263) is located near the design limit, indicating the minimum thickness of the bridge structure. The optimum value of W I is more than twice as large as the optimum value of W O , which indicates that the effect of W I on the power density is far larger than that of W O . The predicted optimized peak power density is 49.85 mW/cm 2 while the CFD result is 50.24 mW/cm 2 , which yielded a relative error of only 0.77%. This small error proves that the use of the RBNN model and GA is suitable for the investigated problem. In terms of improvement from the reference design, the optimal design increased the peak power density significantly by 57.6% from that of the reference design. Figure 6 shows the peak power density curves for the reference and optimal designs of the double-bridge MMFC, which confirms the achievement of the optimization.  In Figure 7, the oxidant (KMnO4) concentration contours of the reference [30] and optimized designs are plotted in the x-z plane 5 mm downstream of the channel inlet. As can be seen in this figure, it is confirmed that the optimal design has an increased inner channel half-width (WI) and a decreased outer channel width (WO) compared to the reference design, which makes the depletion region around each electrode thicker and the mixing region in the inner channel thinner. Therefore, due to the decreased thickness of the mixing region, proton transportation across the mixing layer is promoted, which reduces ohmic losses, and makes the consumption of reactants faster and more efficient. Due to the increase in the width of the inner channel, the distance between the mixing and depletion layers is also increased, and thus the possibility of crossover is reduced in the optimal design. As a sufficient oxidant reaches the activation wall of the channel with optimization, the depletion region expands, increasing mass transportation losses. However, due to the fast reaction rate and reduction of ohmic losses, protons migrated efficiently through the mixing region, resulting in higher current generation and increased peak power density. This was confirmed to be effective in overcoming mass transport losses. In Figure 7, the oxidant (KMnO 4 ) concentration contours of the reference [30] and optimized designs are plotted in the x-z plane 5 mm downstream of the channel inlet. As can be seen in this figure, it is confirmed that the optimal design has an increased inner channel half-width (W I ) and a decreased outer channel width (W O ) compared to the reference design, which makes the depletion region around each electrode thicker and the mixing region in the inner channel thinner. Therefore, due to the decreased thickness of the mixing region, proton transportation across the mixing layer is promoted, which reduces ohmic losses, and makes the consumption of reactants faster and more efficient. Due to the increase in the width of the inner channel, the distance between the mixing and depletion layers is also increased, and thus the possibility of crossover is reduced in the optimal design. As a sufficient oxidant reaches the activation wall of the channel with optimization, the depletion region expands, increasing mass transportation losses. However, due to the fast reaction rate and reduction of ohmic losses, protons migrated efficiently through the mixing region, resulting in higher current generation and increased peak power density. This was confirmed to be effective in overcoming mass transport losses.

Conclusions
The present work optimized the flow channel configuration of the double-bridge MMFC, which showed far better performance than the MMFC employing single-bridge structure for the channel cross-section. To develop clear insight into the impact of depletion regions and mixing regions on the performance of the MMFC, a single-phase numerical model was used. The surrogate-based optimization was performed using GA combined with the RBNN model to maximize the objective function, i.e., the peak power density. The inner channel width (WI) and outer channel width (WO) of the double-bridge channel were selected as design variables for optimization. The optimal design was found at WI/Wr = 1.922 and WO/Wr = 0.815. The results for the objective function predicted by the RBNN model and calculated by CFD analysis at the optimal design showed a very small error of 0.77%, confirming the validity of the employed optimization methods. Compared to the reference design, the peak power density was increased by 57.6% (50.24 mW/cm 2 ) through the optimization. The concentration contours showed that the thickness of the mixing layer in the inner channel decreased, largely reducing ohmic losses by the optimization, but the mass transportation losses increased and the depletion region expanded even though the overall performance was enhanced. Therefore, further research is needed to reduce mass transportation losses in the optimized double-bridge MMFC.

Conclusions
The present work optimized the flow channel configuration of the double-bridge MMFC, which showed far better performance than the MMFC employing single-bridge structure for the channel cross-section. To develop clear insight into the impact of depletion regions and mixing regions on the performance of the MMFC, a single-phase numerical model was used. The surrogate-based optimization was performed using GA combined with the RBNN model to maximize the objective function, i.e., the peak power density. The inner channel width (W I ) and outer channel width (W O ) of the double-bridge channel were selected as design variables for optimization. The optimal design was found at W I /W r = 1.922 and W O /W r = 0.815. The results for the objective function predicted by the RBNN model and calculated by CFD analysis at the optimal design showed a very small error of 0.77%, confirming the validity of the employed optimization methods. Compared to the reference design, the peak power density was increased by 57.6% (50.24 mW/cm 2 ) through the optimization. The concentration contours showed that the thickness of the mixing layer in the inner channel decreased, largely reducing ohmic losses by the optimization, but the mass transportation losses increased and the depletion region expanded even though the overall performance was enhanced. Therefore, further research is needed to reduce mass transportation losses in the optimized double-bridge MMFC.