Multi-Criteria Shape Optimization of Flow Fields for Electrochemical Cells

We consider the shape optimization of flow fields for electrochemical cells. Our goal is to improve the cell by modifying the shape of its flow field. To do so, we introduce simulation models of the flow field with and without the porous transport layer. The latter is less detailed and used for shape optimization, whereas the former is used to validate our obtained results. We propose three objective functions based on the uniformity of the flow and residence time as well as the wall shear stress. After considering the respective optimization problems separately, we use techniques from multi-criteria optimization to treat the conflicting objective functions systematically. Our results highlight the potential of our approach for generating novel flow field designs for electrochemical cells


Introduction
To achieve the large-scale establishment of electrolysis technologies, the optimization of electrolyzers and their components is essential [1].Within these technologies, hydrogenation reactions play a central role.Electrochemical hydrogenation reaction (EChH) emerged in recent years as a more sustainable alternative to the current thermocatalytic state of the art [2][3][4][5][6][7][8][9].Instead of operating at elevated temperatures and pressures, EChH electrolyzers operate at ambient conditions using water and, ideally, renewable electricity to selectively hydrogenate unsaturated bonds.We have recently demonstrated the efficient production of the vitamin synthon 2-methyl-3-buten-2-ol in scalable zero-gap electrolyzers (ZGEs), showing how material selection plays a critical role in the observed activity [10,11].Compared to other electrolytic cell types, ZGEs use a solid polymer electrolyte to conduct ions between the two half-cell compartments and minimize ohmic resistance by mechanically compressing the two electrodes on the polymer electrolyte [12,13].Despite being one of the most promising cell architectures for EChH, there is currently little information on how to improve the efficiency and design of components of ZGEs for electrosynthetic applications [10,11,14,15].
In this work, we investigate how methods from mathematical modeling, simulation, and optimization can be applied to improve the performance of EChH reactors.We focus on the task of optimizing flow field geometries for the cathode side of the EChH reactor.The methods used to achieve this task are computational fluid dynamics (CFD), shape optimization with partial differential equations (PDEs) and multi-criteria optimization (MCO).The goal of shape optimization is to improve the quality of some system by altering its shape, in this case, the shape of the flow field.We employ techniques from shape calculus to perform shape optimization efficiently without requiring parametrizations.Such approaches have been applied in the literature to optimize aircraft [16], electromagnetic devices [17,18], polymer spin packs [19,20], and microchannel systems [21,22].Additionally, there are many recent developments in the field of shape optimization which enable the optimization of complex industrial applications.This can be seen, e.g., in [23], where space mapping methods for shape optimization are introduced, in [24], where the shape optimization software package cashocs is presented, in [25,26], where efficient solution algorithms for shape optimization are investigated, or in [27], where sophisticated mesh deformations for shape optimization are proposed.
Multi-criteria optimization in the context of shape optimization has already received some attention in the literature, e.g., in [28][29][30][31].In these publications, the shapes of the boundaries are parametrized using, e.g., B-splines, which limits the set of reachable designs.However, in this work, we consider freeform shape optimization without parametrizations of the geometry.

Methods and Models
The assembly of an EChH reactor is visualized in Figure 1.We focus our study on the cathode side of the reactor.Water is fed into the flow field geometry and some of it enters the porous transport layer (PTL, marked as cathode support in Figure 1).Then the electrochemical reaction takes place and ions are exchanged through the ion-exchange membrane.Flow field and PTL are modeled in approximation of our recently reported EChH setup with a graphite felt PTL (Sigracell GFD 2.5 EA, SGL Carbon) [10].
We first introduce a model of the flow through the combination of flow field and PTL based on CFD methods.To optimize the shapes of flow fields, we use a simplified model which only considers the flow field without the PTL.Based on this model, we perform a shape optimization with different objectives to derive new flow field shapes.Methods from multi-criteria optimization are used to investigate the combination of multiple conflicting objectives.Finally, the flow field shapes are evaluated with the detailed simulation model.

Flow Field and PTL Simulation Model
The flow through a porous medium can, in general, be modeled on different scales.When the geometric structure of the porous medium is resolved, this is called a detailed or microscale model.However, the use of such a model is expensive.Therefore, the porous medium can be described as effective medium by using a permeability tensor, resulting in a so-called macroscopic model.As we consider flow in both porous and pure fluid domains simultaneously, we employ the commonly used Brinkman model [32].
The geometric setting for the detailed simulation model is shown on the left in Figure 2. The computational domain Ω is split into a pure fluid part  f and a porous part  p .In the following, we restrict our considerations to steady-state, laminar, isothermal, and incompressible flow.
Here  denotes the velocity and  the pressure.The properties  and  describe the density and dynamic viscosity of the fluid, respectively.In the porous part, the viscosity might be different and, therefore, the effective viscosity  eff is used.However, there is no consensus about the definition of the effective viscosity so that mostly the value of  is used [33].The porous medium is characterized by its permeability .The boundary of Ω is divided into the inlet Γ in , the wall boundary Γ wall , and the outlet Γ out .Due to the symmetry of the geometry, we only use half of it and denote the corresponding part of the boundary as Γ sym .Moreover, we denote by Γ int the interface between the fluid and porous part.
Here, ∂   denotes the normal derivative of a vector field , which is defined as ∂   =  , i.e., the application of the Jacobian  to , where the latter is the unit outer normal vector on Γ.
For the computation of velocity and the pressure, the software FiltEST [34] is used.

Flow Field Simulation Model for Shape Optimization
For the shape optimization of the flow field, we restrict our attention to the fluid domain Ω  only and, for the moment, neglect the porous part Ω  .Hence, we employ a no-slip boundary condition on the interface Γ int between fluid and porous part.Our simplified model for the flow in the flow field is given by ⋅  = 0 on Γ sym ,  ∂   ×  = 0 on Γ sym ,  ∂   −  = 0 on Γ out .
For the sake of faster simulation times, we solve the above model with the help of a dimension reduction technique presented in [21], which reduces the above three-dimensional problem to a two-dimensional one, significantly reducing the cost of its numerical solution while sacrificing only little accuracy.We refer the reader to [21] for more details and a thorough investigation of this dimension reduction technique.

Cost Functions for the Shape Optimization
Based on the simplified flow model (3) we optimize the shape of the flow field.To do so, we introduce three cost functions  1 ,  2 , and  3 with the following goals: •  1 aims at achieving a uniform flow distribution among the channels of the flow field, •  2 aims at achieving a uniform residence time of the flow in each channel, and •  3 aims to control the shape of the in-and outflow regions via a criterion based on the wall shear stress.

Uniform Flow Cost Function 𝑱 𝟏
Motivated by the observation that the volume flow rate varies significantly between the channels of the original domain PAR (cf.Section 3.2), we define the cost function Here, Ω channels is the part of the domain Ω f corresponding to the channels and  des is the desired velocity, which we want to achieve within the channels.The latter can be computed a-priori as Poiseuille flow profile for a single channel with the desired average volume flow rate  ̇ (cf.[21]).The cost functional aims to achieve a uniform flow distribution among the channels.

Uniform Residence Time Cost Function 𝑱 𝟐
The second cost function  2 is given by where τ  is the residence time in channel  and τ des is a desired residence time.Hence, the cost functional aims at achieving a uniform residence time, τ des , in each channel.In this paper we restrict ourselves to the case of a uniform residence time for all channels, but using different desired residence times for different channels would be possible.The residence time for channel  is computed as with   and ̅  being the length and average velocity of channel , respectively, with Here,   is the velocity component in channel direction and Ω  is the subdomain corresponding to channel .

Wall Shear Stress Cost Function 𝑱 𝟑
As third objective we consider the wall shear stress on the outer boundaries of the flow field.The goal of this criterion is to prevent the creation of unnecessarily large inflow and outflow volumes during the optimization.This is achieved by controlling the wall shear stress on the boundary part Γ wss = Γ wss in ∪ Γ wss out (cf. Figure 5).The wall shear stress can be seen as a measure of how fast the flow is near a boundary, hence, a low wall shear stress can indicate a stagnation zone.So, a criterion that ensures a sufficiently high wall shear stress can prevent stagnation zones and, thus, unnecessarily large volumes.For this reason, we introduce the constraint  ≥  thr on Γ wss , (8) where  is the wall shear stress, which is defined as Moreover,  thr is a threshold value that should be reached on Γ wss .To treat this constraint numerically, we employ a Moreau-Yosida regularization (cf.[35]) and end up with the cost functional For more details, we refer the reader, e.g., to [20], where a similar criterion is used in the context of designing polymer spin packs without stagnation zones.

Numerical Solution of the Shape Optimization Problems
First, we consider the single criteria shape optimization problem subject to (3), for   being any of the previously defined cost functions  1 ,  2 or  3 .
To solve the optimization problem (11) numerically, we use our shape optimization software cashocs [24], which is a software package for the automated solution of PDE-constrained shape optimization problems.We impose the following geometrical constraints on the problem: The outer boundary of the geometry, including inlet, outlet, and symmetry boundary, is not allowed to change.Additionally, the channels of the geometry are only allowed to change in length.As cashocs is based on the finite element software FEniCS [36], we use the finite element method to discretize the state and adjoint systems with LBB-stable Taylor-Hood elements, i.e., piecewise quadratic and linear Lagrange elements for velocity and pressure, respectively.For the solution of the optimization problem, we use a limited memory BFGS method from [25] which is implemented in cashocs.

Multi-Criteria Shape Optimization
Finally, we can formulate the PDE-constrained multi-criteria shape optimization problem as follows: subject to (3), where we consider the above multi-criteria optimization problem (12) in the sense of Pareto-efficiency (cf.[37]).
To approximate the Pareto front, the calculation of single Pareto points is performed by scalarizing the cost functions to arrive at a single objective function to be minimized.To do so, the weighted-sum scalarization is used [37]: subject to (3), where λ  are non-negative weights.To ensure that the problem is well-scaled, each cost functional is divided by its respective value for the initial domain Ω f,0 with corresponding velocity  0 , where the original flow field shape PAR is used as Ω f,0 .An effective approximation of the Pareto front is obtained by applying the so-called sandwiching algorithm [38,39].Sandwiching relies on an iterative approach that provides, at each step, a weight vector  = [λ 1 , λ 2 , λ 3 ] to be used in (13).New Pareto points are added as long as the desired approximation quality [39] is not achieved.For the sake of brevity, we do not give the details of the sandwiching technique here but refer the reader to [38,39].
The decision maker can then interactively explore the Pareto set by navigating with a graphical user interface [38].Graphical sliders are used for this purpose, each corresponding to one cost function.By moving one slider, the other sliders are updated in real time, i.e., information on the trade-offs between the best compromises is directly visualized.We refer the reader to, e.g., [38][39][40] for more details.

Results and Discussion
We begin our study by investigating the original flow field shape PAR using our detailed simulation model.We then use shape optimization and multi-criteria shape optimization to derive optimized flow field shapes which are compared based on criteria derived from the detailed simulation model.

Simulation Setup
For now, we only consider the flow of water through the flow field and PTL and neglect the electrochemical reaction.Note that, in reality, the electrochemical reaction produces gas as a side product, which mixes with the water and thus changes the flow properties when advancing through the cell.It is up for future work to potentially include this effect into the simulation model.

Simulation Results for Original Flow Field PAR
Simulation results obtained from the detailed simulation model ( 1) are shown in Figure 3.The left plot shows streamlines in the flow field and the right plot shows streamlines in the PTL.Both plots are colored by the velocity magnitude.Note that in the plot on the right, there is an area with few streamlines.This area is still covered by flow, but only reached by a few of the streamlines used for the visualization.
These simulations show an uneven flow distribution: while the flow rate in the middle two channels is high, the remaining channels receive significantly less flow.This also translates to the flow in the PTL, which is also much higher in the middle.Based on this observation, we focus on the improvement of flow uniformity in the following.In addition, we also investigate the influence of other criteria, such as residence time and wall shear stress.

Shape Optimization of Flow Field
To improve the flow properties of the flow field, we solve the single criteria shape optimization problem (11) for each of the previously defined cost functions: we use the uniform flow cost function  1 but also test the influence of the two other cost functions  2 and  3 .Solving the shape optimization problem gives rise to the following shapes: • PAR J1: Uniform flow per channel, minimizing the cost function  1 , • PAR J2: Uniform residence time per channel, minimizing the cost function  2 , and • PAR J3: Sufficiently high wall shear stress on Γ wss , minimizing the cost function  3 .
The resulting flow field shapes are shown in Figure 4. 1) The first column in Figure 5 shows the volume flow rate through each of the 18 channels of the flow field, corresponding to cost function  1 .The volume flow rate is obtained from the velocity integrated over the cross section of each channel.
2) The second column in Figure 5 shows the residence time within each of the 18 channels of the flow field, corresponding to cost function  2 .3) The third column in Figure 5 shows the wall shear stress on the upper wall Γ wss in of the flow field, which contributes to cost function  3 .4) The fourth column in Figure 5 shows the wall shear stress on the lower wall Γ wss out of the flow field, which contributes to cost function  3 .

PAR J1: Flow Field with Uniform Flow
The flow field shape PAR J1 is obtained by solving the single criteria shape optimization problem (11) with cost function  1 (uniform flow).The desired flow rate per channel is chosen to be  ̇des = 6.9 × 10 −9 m 3 s =  ̇ 18 , with 18 being the number of channels.The analysis in Figure 5 shows that it is indeed possible to achieve nearly uniform volume flow rates through all channels.This can also be seen from Figure 4: The velocity in the channels is very uniform, which is not the case for the original flow field shape PAR.The inlet and outlet distributor volumes, before and after the channels, respectively, have increased when comparing PAR J1 with the original flow field PAR.The flow pattern in the inlet distributor is different than in the outlet distributor volume due to inertia effects.Especially in the inlet distributor volume of PAR J1 there is a large blue region (see Figure 4) of stagnating flow, which is not desirable in a flow field shape.This stagnation zone can also be seen in the plot of wall shear stress in Figure 5, where the wall shear stress of PAR J1 is significantly lower compared to the original shape PAR.To prevent such a stagnation zone, we later use a combination of cost functions  1 and  3 .Before doing so, we first continue to study the influence of each cost function on the geometry, individually.

PAR J2: Flow Field with Uniform Residence Time
The shape PAR J2 is obtained by solving the shape optimization problem (11) with cost function  2 (uniform residence time).The desired residence time in the channels is τ des = 3.4 s. Figure 5 shows a very uniform residence time in all channels.However, the inlet distributor volume is greatly enlarged with a large stagnation area which can also be seen from the low wall shear stress on Γ wss in (see Figure 5).The flow rate through the outer channels is even lower than for the PAR flow field.The reason is that the outer channels are shorter than the inner channels and a uniform residence time is achieved by reducing the flow rate in the shorter channels.So, if there are channels with different lengths, the criteria  1 and  2 are contradicting.

PAR J3: Flow Field with Sufficiently High Wall Shear Stress
The shape PAR J3 is obtained by solving the shape optimization problem (11) with cost function  3 .The goal of this cost function is to achieve a sufficiently high wall shear stress on Γ wss , which should be above the chosen threshold value of σ thr = 0.025 Pa.The plots in Figure 5 show that the wall shear stress for the PAR J3 shape is indeed mostly above the threshold value.However, the differences in flow rate and residence time have greatly increased.In Figure 4 we observe that particularly the inlet distributor volume decreases slightly to achieve the threshold value for the wall shear stress, confirming our assumptions in Section 2.3.3.Therefore, the wall shear stress criterion can be combined with another criterion to prevent the creation of stagnation zones, which we do in the following.

Multi-Criteria Shape Optimization
The Pareto front of ( 12) is approximated by solving multiple, i.e., 16, scalarized shape optimization problems (13) with different weight vectors to achieve an approximation quality of 0.01 [39] .The  approximation of the Pareto front is visualized in Figure 6, where the optimized flow field shapes are shown together with the corresponding point on the Pareto front.
We have encountered some issues when computing the Pareto front: For some combinations of weights, the shape optimization algorithm did not fully converge.This happened, e.g., when transitioning from cost function  1 to  2 , and, hence, the point corresponding to PAR J2 is rather isolated.Usually, the reason for the failed convergence of the algorithm is a decreasing mesh quality.The impact of mesh quality is further discussed, e.g., in [27,42,43].In case of a non-converging scalarized shapeoptimization problem (13), the corresponding weight vector provided by the sandwiching algorithm is slightly modified by hand, and then the optimization problem is solved again with the modified weights.

PAR Pareto: Selected Pareto Optimal Compromise
To find an optimal compromise between the considered cost functions, we explore the computed Pareto front and investigate the obtained Pareto optimal shapes.We select the point with objective value  = [0.14, 0.42, 0.73] from the Pareto front and re-solve the scalarized shape optimization problem (13) to obtain the corresponding shape which we denote by PAR Pareto.We choose this Pareto solution, i.e., reactor shape, as it seems to show a good uniform flow distribution, due to the cost function  1 , but with smaller distributor volumes at the in-and outlet compared to PAR J1, due to the wall shear stress functional  3 .

Simulation-Based Benchmarking of Flow Field Shapes
The optimized flow fields investigated in this paper are obtained with the simplified model (3) which does not consider the flow in the PTL.For this reason, we now simulate these flow fields with the detailed simulation model (1) to investigate the influence of the PTL.Results for the comparison of the flow rates per channel are shown in Figure 8. There, some differences in flow rates between the two models are visible.For the flow field PAR J1 the distribution with the more realistic model is not as uniform as for the simplified model, but still good.The reason for the difference in the outer channels is that part of the flow passes through the PTL and can shortcut the way from inlet to outlet.The flow in the middle channels might be higher in some cases since the no-slip condition at the interface is not used here.Better results can be obtained if the model used for shape optimization does also account for the flow in the PTL, which is up for future work.Alternatively, a space mapping approach as demonstrated in [23] can be used.
Finally, Figure 7 shows a comparison of the flow in the PTL between the original flow field PAR and the flow fields PAR J1 and PAR Pareto which have been optimized to make the flow more uniform.The flow velocity is indicated by color on a logarithmic scale.Based on the results shown there, we can conclude that a uniform flow distribution in the flow field also leads to a more uniform flow through the PTL.

Conclusions
We have shown that our shape optimization approach works well and is a flexible tool for designing novel flow field shapes.As optimality criteria we considered the uniformity of the flow and residence time as well as the wall shear stress.We have seen that methods from multi-criteria optimization are well-suited to select compromises from multiple objective criteria.With the Pareto front at hand, one can easily explore and compare Pareto optimal shapes before making an informed decision.The next step is to experimentally test how the flow fields optimized here affect cell efficiency.For the future, we plan to implement criteria that are more directly related to cell efficiency.For this, the flow and possibly reactions within the PTL must be considered in the optimization problem.
2021-0023 and "KataSign" Ko-2021-0016).D. S. is grateful for financial support by the BMBF (Federal ministry of education and research) within the NanoMatFutur Project "H2Organic" (No. 03XP0421)).K. P. gratefully acknowledges a PhD fellowship granted by the Fonds of the Chemical Industry.

𝐷𝑣
[-] Jacobian of a function

Figure 2 .
Figure 2. Geometry for the detailed (left) and simplified (right) flow field simulation model.This is the original flow field shape with parallel channels which is denoted by PAR throughout this manuscript.
The fluid is modeled using a constant viscosity μ = 3.547 ⋅ 10 −4 Pa ⋅ s and a constant density ρ = 971.79  3 , which corresponds to water at 80°C.The inlet velocity is chosen so that we have an inlet flow rate of  ̇in = 7.5  .The permeability of the PTL is estimated as  = 10 −11  2 from[41].

Figure 5
Figure 5 compares the individual quantities which give rise to the cost functions:

Figure 3 .
Figure 3. Simulation of the original flow field PAR using the detailed simulation model.Streamline visualization of velocity magnitude for flow field (left) and PTL (right).

Figure 4 .
Figure 4. Overview of the flow field shapes considered in this paper.The flow velocity is indicated by color using a logarithmic scale.PAR -original shape, PAR J1 -uniform flow, PAR J2 -uniform residence time, PAR J3 -wall shear stress, PAR Pareto -Pareto optimal compromise.

Figure 5 .
Figure 5.Comparison of the flow quantities (velocity, residence time, and wall shear stress) corresponding the optimization criteria  1 ,  2 , and  3 for the five flow field shapes considered.

Figure 7 .
Figure 7.Comparison of the flow in the PTL for three different designs.

Figure 8 .
Figure 8.Comparison of the flow rate per channel between simplified model (denoted 2D) and detailed model (denoted 3D).
s] fluid velocity   [m/s] -component of  ̅ [m/s] average velocity  ̇ [m³/s] volumetric flow rate   [-] set of weights for the scalarization of the multi-criteria optimization problem interface between fluid and porous medium out belonging to the outlet p belonging to the porous part of the domain sym belonging to the symmetry boundary thr threshold wall belonging to the wall boundary wss belonging to the part of the boundary where the wall shear stress is considered