Topology optimization for uniform flow distribution in electrolysis cells

In this paper, we consider the topology optimization for a bipolar plate of a hydrogen electrolysis cell. We present a model for the bipolar plate using the Stokes equation with an additional drag term, which models the influence of fluid and solid regions. Furthermore, we derive a criterion for a uniform flow distribution in the bipolar plate. To obtain shapes that are well‐manufacturable, we introduce a novel smoothing technique for the fluid velocity. Finally, we present some numerical results and investigate the influence of the smoothing on the obtained shapes.


Introduction
Hydrogen plays an important role for achieving a climate neutral industry as well as climate neutral means of transportation.For that the hydrogen needs to be produced, for example, by means of hydrogen electrolysis cells.Typically, those cells split water into oxygen and the desired hydrogen by using (green) electrical energy.We focus on a special kind of electrolysis cell here, the so called proton exchange membrane (PEM) electrolysis cell [1], which works as follows: First, water enters the cell on the anode side and is distributed over the cell by the bipolar plate.The water wanders through the so-called porous transport layer (PTL) to the membrane where the electrochemical reaction takes place.The oxygen is then discharged from the cell again by the bipolar plate.Whereas positive hydrogen ions can pass through the membrane and form hydrogen at the cathode side.To guarantee the efficiency of such a PEM electrolysis cell, the flow through the anode side bipolar plate is essential.The water needs to be distributed uniformly over the whole cell so that the entire membrane area is reached.For more details on PEM electrolysis cells, we refer the reader, e.g., to [1].
Here, we only investigate the anode side of the bipolar plate.For simplicity, we neglect the electrochemical reaction and phase changes and only consider the distribution of a homogeneous fluid.To obtain a uniform flow through the bipolar plate, we use techniques from topology optimization based on topological sensitivity analysis [2].
Topology optimization considers the optimization of a domain by changing its topological features by either adding or removing material.It was first introduced in the context of solid mechanics, but has since been applied to a wide variety of applications, for example, in elasticity [3,4] and fluid mechanics [5,6].We use the concept of the topological derivative, which was first introduced in [3] and later mathematically justified in [2].The topological derivative measures the sensitivity of a shape functional under infinitesimal topological changes.There are different approaches for solving topology optimization problems numerically and here we utilize a level-set approach which was introduced in [4].
This paper is structured as follows.In Section 2, we introduce a model for the anode side bipolar plate using the Borvall-Petersson model from [5].Additionally, we derive a criterion for a uniform flow distribution in the bipolar plate.To obtain geometries suitable for manufacturing, we consider a smoothing of the fluid velocity which is incorporated in the cost functional for the optimization.In Section 3, we summarize some basic concepts of topology optimization including the topological derivative as well as a gradient-based algorithm from [4].Further, we formally present the topological derivative of the cost functional for our problem.Finally, in Chapter 4, we numerically solve the topology optimization problem and present our numerical results.Particularly, we focus on the influence of the smoothing on the optimized geometries.Our results show that topology optimization can be used to create novel bipolar plate designs.

Model problem for the anode side bipolar plate
Let us begin with presenting our mathematical model for the anode side of the bipolar plate as well as the corresponding topology optimization problem for a uniform flow distribution.

The Borvall-Petersson model
First, let the hold-all domain D ⊂ R d be an open and bounded set.In this paper, we focus on the case d = 2, but an extension to d = 3 is also possible.Additionally, let Ω ⊂ D be measurable and open.Then, we identify Ω as the fluid region and D \ Ω consequently as the solid part of D. To model the fluid flow, we use the Borvall Petersson model, which was first introduced in [5] for topology optimization of fluids.We assume that the boundary Γ = ∂D of the hold-all domain D is divided into three parts: The inlet Γ in , where the water enters the plate with velocity u in , the boundary Γ wall , where we have a no-slip condition and Γ out , where a natural outflow condition is applied.The model reads where u and p denote the fluid velocity and pressure, respectively.Here, n denotes the outer unit normal vector on Γ.To distinguish between the solid and fluid part of the domain the inverse permeability α is used, which is given by where α L and α U are positive constants.In particular, α is chosen small inside the fluid region and large in the solid part.For a detailed derivation of the model we refer to [5].

A criterion for uniform flow distribution
To increase the efficiency of the electrolysis cell, we want to assure that the water gets distributed uniformly all over the cell.Without knowledge of the shape of Ω, it is not obvious how to characterize a flow as uniform mathematically.Therefore, we take a different approach and introduce a threshold velocity magnitude u t > 0. The goal is now to achieve a flow whose magnitude reaches at least the threshold.This criterion ensures that each part of the bipolar plate receives a sufficiently large flow and no dead spots occur.Additionally, as our model also yields, albeit small, flow velocities in the solid part, we have to restrict our criterion to the fluid part Ω.The optimization goal then reads where ∥•∥ denotes the Euclidean norm on R d .There are two aspects to consider here.First, the target velocity constraint might not be fulfillable close to the boundary ∂Ω of the fluid domain due to small velocities in the solid part.Additionally, as we do not have any restrictions on the solid part, large solid areas might appear in the design structures.In fact, we observed this behavior in numerical tests.Usually, we want to avoid solid structures above a certain size as the flow inside the porous transport layer is hindered by such objects and this leads to a degradation of the cell efficiency.

Smoothing of the velocity
To avoid large solid structures and diminish the effect of small flow velocities in the solid parts, our goal is to extend the target velocity goal (2) to the hold-all domain D. We consider a smoothed flow velocity which is defined on the entire hold-all domain D. For that we use the smoothing characteristics of the heat equation.We consider the following equation Here, we have the velocity solution u of the Stokes Darcy equation (1) as the initial condition.To reduce the numerical effort, we discretize (3) with a single implicit Euler step with step length ∆t > 0 and arrive at By modifying the time step length ∆t we can control the influence of the smoothing, where the smoothing effect gets larger for increasing step lengths.With the smoothed velocity u s , we can now extend the criterion (2) to the entire hold-all domain by considering Finally, we apply the Moreau-Yosida regularization [7] to (5) to arrive at our objective function Remark 2.1 Using the smoothed velocity has several advantages for the topology optimization.First, when we consider a small solid inclusion, we expect that it does not severely effect the flow in the underlying porous transport layer.In the smoothed velocity u s with a suitable time step ∆t, the effect of the small solid inclusion will not be resolved anymore, so that the smoothed velocity u s is not influenced anymore by the small velocity inside the solid region and the target velocity goal (5) will be fulfilled.On the other hand, if we consider a large solid inclusion, it will still be resolved by the smoothed velocity u s and the criterion (5) will not be satisfied.Therefore, large solid structures can be avoided with this approach.

Topology optimization problem
We summarize (1), ( 4) and (6) to state our topology optimization problem such that (1) as well as (4) are fulfilled in the weak sense, Here, we introduce an additional volume constraint for the fluid volume with V L and V U being the lower and the upper border, respectively.The goal is now to perform a topology optimization for problem (7) using the topological derivative.

Topology optimization
We seek to apply a gradient-based topology optimization approach for (7).For that we introduce the topological derivative and state the corresponding topological derivative for the objective (6).Furthermore, we briefly present the gradient-based topology optimization algorithm we use for our numerical experiments.

Topological derivative
We consider an analogous setting to before, where D is an open subset of R d .We define A = {Ω ⊂ D s.t.Ω open} as the set of all open subsets of D. To present the topological derivative, we consider a general shape functional S, given by Additionally, we assume to have a fixed shape Ω and a fixed point z ∈ D \ ∂Ω.The topological derivative then measures the sensitivity of the shape functional S with respect to infinitesimal topological changes of the shape Ω around the point z.We characterize the perturbed shape by where ω z,ϵ = z + ϵω with ω ⊂ R d and 0 ∈ ω represents the shape of the perturbation.Furthermore, we assume to have a positive function l with lim ϵ→0 l(ϵ) = 0.The topological derivative, see for example [2], is then defined as if the limit exists.Typical choices are w = B 1 (0) with l(ϵ) = |ω z,ϵ |, where B 1 (0) is the open unit ball with center 0.

Topological derivative for our model problem
We formally present the topological derivative of our model problem (7), which can be derived, for example, with an averaged adjoint approach [8].The topological derivative for (6) is given by for all z ∈ D \ ∂Ω.Here, u is the weak solution of (1) and v is the adjoint velocity which solves and v s is the adjoint of the smoothed velocity which solves the equation The actual computation of the topological derivative is beyond the scope of this paper and a topic of future research.For more details regarding topological sensitivity analysis, we refer the reader to [2].

A gradient-based solution algorithm for topology optimization
We present an algorithm using the topological derivative as a search direction in order to solve topology optimization problems.
The method was introduced by Amstutz and Andrä [4].For the sake of brevity, we only give a short presentation here and refer the reader to [8] for an overview over established topology optimization algorithms as well as novel quasi-Newton methods for topology optimization.First, we introduce a continuous level-set function ψ : D → R to characterize the fluid as well as the solid part of the hold-all domain D by Additionally, we define the generalized topological derivative The algorithm then updates the level-set function consequently by a linear combination of itself and the generalized topological derivative on the L 2 -sphere.The update formula reads Here, θ n is the L 2 -angle between the level-set function and the generalized topological derivative .
Additionally, κ n is chosen by a line search approach to guarantee that the objective function value decreases each iteration.This procedure is carried out until a local optimality condition is reached.Such an optimality condition has been derived e.g. in [9].Numerically, the algorithm is stopped, when the angle θ n becomes smaller than a certain numerical tolerance ϵ θ .This then implies a local optimality condition, again the reader is referred to [4] or [9].We want to conclude by stating that all ψ n fulfill ∥ψ n ∥ L 2 (D) = 1 if the condition holds for the initial iteration.

Numerical results
We turn back to optimization problem (7) and describe our numerical setting.First, the hold all domain is given by D = (0, 1) × (0, 1).We define a quadratic inflow profile u in given by u in = − 400 9 (y − 0.35)(y − 0.65) 0 for x = 0 and 0.35 ≤ y ≤ 0.65 on the inflow boundary Γ in .The situation is displayed in Figure 1.The inverse permeability α for the fluid as well as the solid part as is defined as Furthermore, the stopping criterion for the algorithm is given by ϵ θ = 0.035 which is equivalent to an angle of 2 degrees.The target velocity is chosen as u t = 0.1 and the lower and upper values for the volume constraint are given by V L = 0.5 and V U = 0.7, respectively.

Numerical implementation
We give a short overview of the software used for the numerical implementation.First, all underlying PDEs are solved using the finite element package FEniCS [10].We use a uniform triangular mesh consisting of 10.000 elements.The Stokes Darcy equation (1) is discretized using LBB-stable Taylor-Hood finite elements which consist of quadratic Lagrange elements for the velocity and linear Lagrange elements for the pressure.Additionally, for the smoothing equation (4), continuous quadratic Lagrange elements are used.Furthermore, we use the software package cashocs [11] for adjoint computations.
Cashocs is an open source software package with a python interface which is based on FEniCS.The software can be used to solve optimization problems constrained by partial differential equations in the context of shape optimization and optimal control in an automated fashion.Furthermore, cashocs implements a discretization of the continuous adjoint approach which allows for adjoints to be computed automatically.
As we also introduced a volume constraint for (7), we give an explanation on how the constraint is handled numerically.After each update of the level-set function we perform a projection of the level-set function onto the set of the admissible shapes.First, we assume to have a newly computed update of the level-set function ψ n+1 , which was generated with (12), and the corresponding shape Ω n+1 .Additionally, we assume that the volume constraint is not fulfilled, otherwise nothing needs to be done.Without loss of generality we assume that |Ω n+1 | > V U .To perform the projection, we move the level-set function upwards by a strictly positive constant until the upper border V U of the volume constraint is reached.This means we search for a c > 0, to move the level-set function ψn+1 = ψ n+1 + c, such that the corresponding shape Ωn+1 fulfills | Ωn+1 | = V U .The case |Ω n+1 | < V L is handled analogously by taking a strictly negative constant to alter the level-set function.Numerically, we compute the constant c by a bisection approach with numerical accuracy ϵ c = 10 −4 .As we do not need to solve any additional PDEs, the computational cost for the projection procedure is reasonable.

Numerical results
We start by choosing ∆t = 10 −3 as the step length for the smoothing equation (4).The resulting optimized geometry and flow velocity are displayed in Figure 2. The algorithm took 137 iterations to reach the target accuracy.The objective function value is given by 8.6•10 −8 .From the obtained velocity fields and the previous investigations, we observe that the optimization algorithm indeed produced a geometry with a very uniform flow distribution.To do so, six fluid channels have been formed and the flow is distributed rather uniformly between them.In more detail, the target velocity goal (5) for the smoothed velocity u s is reached on 93.84% of the whole domain.Furthermore, this extends to a 91.79% fulfillment of (2) for u on the fluid area.Additionally, all obstacle sizes in fluid flow direction are about the same order of magnitude.
We want to investigate the obstacle sizes for different values of the smoothing parameter ∆t.For that we choose two additional values for the step length given by ∆t 1 = 5 • 10 −3 as well as ∆t 2 = 5 • 10 −4 .The results can be seen in Figure 3 for ∆t 1 and in Figure 4 for ∆t 2 .It is apparent that the solid structures for the higher time step length get larger which can be explained by the higher influence of the smoothing.The larger the smoothing parameter is chosen, the larger a solid inclusion has to be in order to be resolved in the smoothed velocity u s .In particular, we now observe that only 3 major fluid channels can be found in the optimized geometry.Here, the target velocity goal (5) is fulfilled on 98.26% of the domain D. For the smaller time step length ∆t 2 on the contrary, the solid structures get smaller, which is again a consequence of the smaller influence of the smoothing and we observe that the geometry branches out into several channels which are then merged back together towards the outlet.We achieve a 89.46% realization of the target velocity goal for the smoothed velocity u s .

Conclusion
Summarizing, we presented a model for reaching a uniform flow distribution in the anode side bipolar plate for hydrogen electrolysis cells.We showed in numerical tests that a uniform flow distribution is indeed reached in the optimized shapes.Additionally, we investigated the influence of the time step length ∆t on the size of the solid structures, giving us a tool to control the object sizes in the resulting shapes.

Fig. 3 :Fig. 4 :
Fig. 3: Result of the topology optimization problem (7) with ∆t1 = 5 • 10 −3 .In Figure a the final shape is displayed, in Figure b the corresponding norm of the velocity field and in Figure c the norm of the smoothed velocity field.