Seeking minimum entropy production for a tree-like flow-field in a fuel cell

Abstract

Common for tree-shaped, space-filling flow-field plates in polymer electrolyte fuel cells is their ability to distribute reactants uniformly across the membrane area, thereby avoiding excess concentration polarization or entropy production at the electrodes. Such a flow field, as predicted by Murray's law for circular tubes, was recently shown experimentally to give a better polarization curve than serpentine or parallel flow fields. In this theoretical work, we document that a tree-shaped flow-field, composed of rectangular channels with T-shaped junctions, has a smaller entropy production than the one based on Murray's law. The width w₀ of the inlet channel and the width scaling parameter, a, of the tree-shaped flow-field channels were varied, and the resulting Peclet number at the channel outlets was computed. We show, using 3D hydrodynamic calculations as a reference, that pressure drops and channel flows can be accounted for within a few percents by a quasi-1D model, for most of the investigated geometries. Overall, the model gives lower energy dissipation than Murray's law. The results provide new tools and open up new possibilities for flow-field designs characterized by uniform fuel delivery in fuel cells and other catalytic systems.

---

1 Introduction

Since the start of polymer electrolyte membrane fuel cell (PEMFC) research,¹ serpentine flow-fields have been commonly used to supply reactants (oxygen and hydrogen) to the cell's catalytic layers in the membrane electrode assembly (MEA). The serpentine field is in use^2–4 also in industry, but it has become increasingly clear that other flow-fields perform better, in terms of a better polarization curve and lower operational costs.^5,6 The distance from the polarization curve to the value of the reversible potential of the cell expresses the energy dissipated as heat or the entropy production times the temperature of the fuel cell.⁷ The losses, which vary with the current density, are much larger at the electrodes and in the membrane than in the porous transport layers (PTLs).⁸ However, the losses in the PTLs will have an impact on the cell's overpotentials via the gas concentration at the catalyst. A flow-field that can deliver reactants at uniform conditions to the cell's active membrane area is therefore beneficial. An optimal field can thus be characterized by a uniform entropy production. Among numerous replacement proposals, it has been argued that bio-inspired fractal-type flow-fields would be beneficial, because they are developed for energy-efficient biofluid delivery, meaning uniform delivery of fluids to animal and plant tissues.^9,10 Recent experimental evidence confirms their superiority in this respect for PEMFCs. Fractal-type flow channels could have an open side like other conventional designs and, thus, provide a continuous non-uniform gas delivery to the membrane along the flow-field. Trogadas et al.¹¹ proposed a 3D lung-inspired design with closed channels that is perpendicular to the PEMFC membrane and is only in contact with the MEA through the openings of the smallest channels, which supplied the active area with a uniform gas mixture from above. The flow-field was better than a serpentine flow-field in terms of its polarization curve. The property of equipartition of entropy production has been connected with minimum energy dissipation in various process equipment (for an overview, see ref. 12). This altogether has led us to speculate that fractal-like (or self-similar) patterns with constant scaling parameters could be beneficial as distributors of gases in PEMFCs.¹³ Perfect equipartition can be obtained using scaling parameters obtained with the method of Gheorghiu et al.¹⁴ and Trogadas et al.,¹¹ but, in reality, there are often boundary- or material conditions that do not allow us to reach the ideal limit. Additional constraints or variables are therefore interesting. We shall study such here, comparing the outcome all the time to the performance of Murray's law for volume or area-filling conditions,¹⁵ see Appendix 1. This condition is also close to the practical condition facing a fuel cell designer. Murray¹⁵ used a volume-filling condition as a constraint for flow in tree-shaped structures, and obtained as an outcome of the optimization his famous scaling law, saying that the diameters of the branches from one generation to the next were scaling as eqn (1) shows:where d₀ is the diameter of the parent branch, and d_j are the diameters of the n daughter branches belonging to the next generation level.

Eqn (1) characterises the fluid delivery system for which the total entropy production of the flow-field pattern by viscous dissipation is minimum, given the total volume available to flow.¹⁵ An optimization problem with minimum total volume at a given entropy production, has, however, the exact same solution. Both formulations are beneficial to use in flow-field plate design. Minimum total volume at a constant channel depth means maximum gas–land contact area and, therefore, smaller total electric resistance of current collectors (smaller Ohmic losses). In this sense, it is logical to use constant volume as a constraint in the optimization problem.

Gheorghiu et al.¹⁴ studied airflows in the bronchial tree of humans and used as a constraint that the volume flow through each generation of branches was constant. They found that minimum entropy production corresponded to a constant pressure drop across each branch. In their case, minimum entropy production could be understood as equipartition of entropy production (uniformly distributed entropy production).

The geometric design principle of Murray¹⁵ was seen as a manifestation of equipartition of entropy production.¹⁴ The constraints used by Gheorghiu et al. and Murray were not the same, however. In the first case, no restriction was placed on the geometry, while in Murray's case the volume-filling geometry compatible with minimum work was of primary interest.

From studies of Nature, we know, however, that a decrease in the entropy production in one a part of a complex system could be beneficial for the overall efficiency. This was confirmed by an enormously high number of studies of fluid delivery systems in Nature (in plants, animals, sponges, fungi, etc.). They were shown to have minimum dissipated energy (entropy generation), independently of other functions.^9,10 Therefore, we shall assume that an approach to minimize the viscous dissipation in the flow-field also could be beneficial for PEMFCs, in spite of non-constant temperatures across the cell and other types of dissipation. Previous studies have already shown that the entropy production of a fuel cell can be optimized by separating the flow-field plate and MEA if the gases are delivered uniformly. In our design the fractal channels are in contact with the MEA only through the fixed number of open ends of the last generation's channels. This allows one to optimize the flow-field plate in the x–y direction and the MEA in the z-direction without any reiterations to connect both.¹³ This hypothesis needs, of course, to be checked by experiments and/or more detailed multiphysics CFD computations, but a successful starting geometry could well come from the simplest base.

There are presently few quantitative reports on bio-inspired designs,^11,13,16 and the present work attempts to change this situation, contributing to analytic and numerical studies of biomimetic designs of the flow-fields in PEMFCs. The problem is challenging, as the transport of oxygen can be severely hindered by water formation in the catalytic layers and water clogging in the PTL. In addition, the heat production at the cathode can shift the ratio of water vapor to liquid.⁸

But, in this context, a tree-like structure for supply of gases may offer advantages. Trogadas et al.¹¹ used a 3D lung-inspired tree perpendicular to the PEMFC membrane, supplying the active area with the same gas mixture from above. We shall use a “flat” tree, or a quasi-3D-tree that is parallel to the membrane, see Fig. 1, as in the original proposal of Kjelstrup et al.¹³ The inlet is set at the symmetry line of the area close to its perimeter (Branch 0, 1 in Fig. 1), while the outlets are (only) at the ends of all the branches of the last generation (here 4 generations are shown). Both types of trees will, by construction, deliver gas at (nearly) uniform conditions to the membrane active area. The “flat” tree may also be embedded into the end-plates of the cell, and possibly help us avoid the PTL all together. There is a possibility to construct parallel outlet channels for water, with a similar structure operating in the reverse direction.

Fig. 1 Top: Schematic picture of the tree-shaped flow-field in a layer parallel to the membrane. The picture shows the numbering of branches i and generation levels j from the inlet 0 to the outlet. When the gas outlet is restricted to the ends of the branched tree, gas is delivered at similar conditions. Bottom: Self similar fractal-like pattern with k-parameter 2 (left) and 3 (right).

These potential options have motivated this first study of a flat tree-like flow-field. We are concerned with a gas inlet distributor where hydrodynamic flows take place in a plane, with a certain channel depth and width. The Peclet number will be used to assess how far we are from the diffusional regime. The model applies so far to the anode side of the PEMFC, where hydrogen and water vapour flow in the same direction, and where there is negligible water condensation. In order to use it for the cathode, we need to introduce two-phase flow of water and gas with the possibility of water clogging the pores in an outlet system with the same shape. This problem will be postponed, however. The paper is organized as follows. We describe first (Section 2) the quasi 3-dimensional (3D) tree-shaped flow-field pattern, with a finite number of branches, to set a practical stage for the investigation. Afterwards we define the equations which describe the system in use. The 3D system describing the flow in the fractal flow-field is solved numerically, while a simplified 1D representation of the same is solved analytically. The flow-fields are studied under flow and boundary conditions typical for PEM fuel cells. The impact of geometry is pointed out. In our study, we find that there are other geometries than those given by Murray that have a smaller energy dissipation. We proceed to explain the benefits/disadvantages of various case studies, followed by a discussion of the results. We believe that our systematic use of the entropy production as a tool to observe the performance of a flow-field is new in the context of fuel cells.

2 System

The flow-field of a fuel cell should preferably deliver gas at uniform conditions all over the 2D-catalytic layer. This motivated our choice of a tree-like structure where each branch always splits into 2 sub-branches. The cross-section of the branch is rectangular or square. Such patterns, which are the focus here, can be characterized by geometric series. The system can be designed to fill the space available, be it a certain volume or, in the quasi 3D-case, an area.

The uniform supply of a fluid to the catalyst is the first essential requirement of an optimal PEMFC.¹³ A working quasi-3D flow-field needs inlet- and outlet-patterns. Only one pattern leads to dead-end flow channels. However, we shall now look at the inlet pattern only for the sake of simplification, where only the branches in the last generation are connected to the gas diffusion layer. The situation shares similarities with an interdigitated flow-field design, where the flow needs to go through the gas diffusion layer before the gases and water are able to leave the fuel cell.

The outlet pattern can be machined onto a new plate, making it a 2-layer flow-field plate. Machining can be done with conventional milling methods since we propose to use a 2D channel structure with constant channel depth. Constant depth allows an easy milling process and reduces pressure loss effects from a step-wise decrease of the channel depth at T-junctions. Any 3D printing will, therefore, not be needed. The branches of the last generation of the inlet pattern will be in contact with the openings drilled in the second plate. One advantage of this type of flow-field plate is the modularity, meaning that the outlet pattern can be easily changed independently of the inlet pattern and vice versa. The first generation of both inlet and outlet patterns will need to be connected to the fuel inlet and gas outlet of the fuel cell, and this can be done by a straight channel.

Take the symbol l₀ for the length of the first branch which is not yet split (level 0). When N_j is the sum of the number of branches at level j, the ratio N_j/N_j+1 is therefore 1/2. The parameter j_f is defined as the generation level number of the last generation. In the first step we are looking at a self-similar system, meaning that the length, l_j,i, and width, w_j,i, of each branch i at generation level j are scaled down by the same factor, k and a, respectively. Widths and lengths have the unit m.

The numbering of branches is illustrated in Fig. 1. We show in Appendix 1 how the length scaling parameter can be derived for volume-filling or area-filling conditions. The branch length will scale as:

with k = 3 for volume-filling constraints, and k = 2 for area-filling constraints. The width scales according to a geometric series, cf.eqn (3):

w_j,i = a^jw₀ (3)

where a is the dimensionless diameter scaling parameter and w₀ is the width in the first branch of generation level 0. Another and quite common way to scale the diameter or width in biological flow systems is to use Murray's law¹⁵ (eqn (4))where d₀ is the diameter of the first branch (generation level 0) in m and d_j,i is the diameter of the branch at generation level j in m. The flow rate Q_j,i of each branch can be calculated with the following equation:where Q₀ is the flow rate in the first branch of generation level 0 in m³ s⁻¹.

For the case of a rectangular-shaped (quasi 3D-) flow channel, we chose to replace d₀ and d_j,i with w₀ and w_j,i. From eqn (4), we calculated a value for a, which is 0.5^1/3 ≈ 0.79. This value will be used for all calculations where we refer to Murray's law. The k-value has an impact on how the pattern fills out the free space in the plane. We shall study values where 2 < k < 3 in connection with Murray's law for width-scaling cases. The endpoint values have clear geometric interpretations. The shape of other fractal-like patterns will also be investigated, changing the k-parameter from 2 up to 3.

Fig. 1 (bottom two figures) shows the space-filling ability of such a tree as a function of the k-parameter (the width was not scaled in this figure). The number of generations in the tree was set to 9 for visual purposes, which means that j_f = 8. Values below k = 2 created patterns that did not fill the space. Values of k > 2.2 create cross-overs. This means that the pattern will turn out to be problematic in a manufacturing process that uses 3D techniques like 3D-printing. The analysis (Fig. 1) of the length scaling parameter k demonstrated the biggest disadvantage of this way to scale the pattern: this always results in a rectangular area. The PEMFC area to be investigated was however a square. Therefore, we set the branch length at a predefined value, constant for all cases, when we asked for the space filling pattern. The only parameter to be scaled was then the width of the rectangular channel. Table 1 shows the length of the different generation levels which were used in the calculations (j_f set to 4), chosen in such a way that we obtain an area-filling pattern for a 25 cm² active fuel cell membrane area.

Table 1 Lengths chosen to give a pattern that fills a 25 cm² active fuel cell membrane area

Generation level j	Length [mm]
0	24
1	12
2	12
3	6
4	6

---

3 Theory

The “flat” tree flow-field will be modelled in several ways. In the most precise manner, we consider all three dimensions of each channel, and compute the flow-field in the full network of branches. The entropy production due to viscous dissipation at isothermal conditions in the 3D-case is:Here is the entropy production of branch i at generation level j of the flow-field pattern in J s⁻¹ K⁻¹, Π is the viscous stress tensor in Pa, V_i,j is the volume of branch i at generation level j in m³, and v is the barycentric velocity in m s⁻¹. The total entropy production is the sum of the entropy production of all branches:For the 1D-tree, the expression for the entropy production simplifies. For one branch, the driving force is minus the pressure drop across the branch divided by the temperature. This gives:Here Q_j,i is the volume flow rate in m³ s⁻¹ in branch i at level j, Δp_j,i the pressure drop across branch i at level j from the outlet to inlet in Pa and T the temperature of the system in K. The entropy production of one generation level is calculated withwhere Δp_j,i is the pressure drop in branch i at level j in Pa and N_j is the number of branches at generation level j, which is 2^j. The total entropy production (TEP) of the flow-field pattern is then simply the sum of the values obtained from eqn (9) over all generation levels in the fractal-like pattern.Here j_f is the maximum number of generation levels in the tree-shaped pattern. To compare different geometries, the specific entropy production in J m⁻¹ s⁻¹ K⁻¹ is introduced:The total specific entropy production (TSEP) is then again the sum of the specific entropy production of each branch in the flow system:We introduce the assumption of Poiseuille flow for cylindrical flow channels^11,14 and obtain:where μ is the dynamic viscosity of the fluid in Pa s. We shall compare two expressions for the diameter in eqn (13). In the first case (Case D_h) we used the common hydraulic diameter (d_j,i = D_j,i):¹⁷In the second case (Case Equivalent A) we used the diameter of a cylindrical cross-section, which gives the same area as the rectangular or squared channel, leading to:In addition to the use of Hagen–Poiseuille flow for the cylindrical flow channel, combined with the first two cases, we also used the analytical solution of Hagen–Poiseuille flow for a rectangular channel¹⁷ (Case Rectangular):The total pressure drop along the tree branches, Δp, is the sum of the pressure drops across the generations j. See Fig. 1 for the branch notation.

4 Methods

4.1 Case studies

Table 2 shows a summary of the conducted studies on the chosen tree-shaped pattern. To start, we studied the variation in the entropy production with flow rate, and the impact of the pressure drop calculation method (Study 1 in Table 2), followed by an analysis of Murray's law and the impact of it on the entropy production (Study 2). We studied the impact of the geometry and scaling properties on the TEP (eqn (9)) and TSEP (eqn (12)) of the flow-field pattern to answer the question: is Murray's law the most efficient way to scale the pattern in terms of entropy production?

Table 2 List of conducted case studies

Study	Investigation	Variables	Method
1	Q 0 and Δp calculation method dependency on TEP and TSEP	Q 0, Δp calculation method, w0	1D
2	Murray's law and entropy production	a, w0	1D
3	1D Δp	a, Δp calculation method	1D
4	Flow rate distribution	a, w0	3D
5	3D TEP and TSEP	a, w0, Q0	3D
6	a, comparison of 1D and 3D Δp	a, Δp calculation method, 3D Δp	3D, 1D
7	Peclet number	L, w0, a, Q0	1D

---

All three cases described in Section 3 were used and compared to each other (Study 3). It was of interest to see how the pressure varied along the branches for different scaling parameters. The 1D-calculations of the pressure drop in the tree-shaped pattern (with the equations given in Section 3) do not include entrance length effects or flow retardation at the T-junctions, thus promoting the 3D calculations (Studies 4–6).

Due to branching, the flow inside the channels can become asymmetric (see ref. 18 and 19), and this may cause an undesirable nonuniform flow distribution over the PEMFC. Therefore we computed the outlet flow rates from the last generation and compared them to each other, in order to see if we were able to reproduce the results of Ramos et al.,¹⁹ and to prove that tree-shaped patterns deliver the fuels uniformly (Study 4).

To answer how well the simple 1D method can capture a more advanced result, we performed 3D simulations. We examined the pressure drop prediction of the 1D-model and evaluated its impact on the entropy production (Study 5 and 6). Finally (Study 7) we computed the Peclet number (cf.eqn (18)) at the outlet of the last branch. The Peclet number defines the ratio between convective and diffusive flow. A Peclet number lower than one means that diffusion dominates the flow. We have:

where L is the characteristic length in m, v is the flow velocity in m s⁻¹ and D is the diffusion coefficient of gas relative to the wall in m s⁻². A Peclet number smaller than 1 in the last branch helps provide a uniform distribution of the fluids to the catalyst.²⁰

4.2 Computational methods

4.2.1 1D-calculations. The 1D calculations were done using eqn (6)–(17) in Section 3. The equations were solved in MATLAB R2019a for various w₀ and scaling parameters a. We used pure oxygen as a flow medium. The temperature was set at 353 K, a common temperature in fuel cell experiments. The system was assumed to be isothermal working at a constant flow rate condition.¹¹ The flow rates, which can be seen in Table 3, were calculated using Faraday's law. The fuel cell active area was 25 cm², the density of oxygen was 1.09 kg m⁻³, the stoichiometric coefficient was 3, and the molar mass was 36 g mol⁻¹. The flow pattern outlet has atmospheric pressure and the inlet pressure is then calculated with the pressure drop along the branches. The viscosity of oxygen at the given operating conditions was calculated with Sutherland's formula and was 2.1 × 10⁻⁵ Pa s. To investigate the impact of the channel width, the channel depth was kept constant at 1 mm.

Table 3 Flow rates at the inlet of the branch in generation 0 used in the computations

Current density [A m^−2]	Flow rate [m^3 s^−1]
10 000	5.7104 × 10^−6
5000	2.8552 × 10^−6
1000	5.7104 × 10^−7
100	5.7104 × 10^−8

---

4.2.2 3D-calculations. The 3D calculations were done using OpenFoam 4.1. The simpleFoam solver was used to solve the Navier–Stokes equation for isothermal, incompressible, single-phase and steady-state flow. Meshes of the tree-shaped patterns were created in Ansys Workbench with a fully hexahedral mesh. Models for five different widths (1, 1.5, 2, 2.5 and 5 mm) and 3 different width scaling parameters (0.79, which was scaling according to Murray's law, 0.9 and 1) were created and computed for the same 4 flow rates (Table 3) which were used in the 1D calculations. The same viscosity and temperature as in the 1D studies were used. To reduce the computational time, the tree-shaped pattern was split in half, with the symmetry plane being at the axisymmetric line of generation 0. Afterwards the pressure drop was evaluated and the specific entropy production computed in ParaView-5.6.0. The entropy production was computed from eqn (6).

5 Results and discussion

5.1 TEP dependence on flow rate and total pressure difference

The results for the TEP and TSEP of the flow-field pattern (Study 1) are shown in Fig. 2 as a function of channel width, for flow rates corresponding to a current density of 10⁴ and 100 A m⁻², and the three ways to compute the pressure drop. The value of a was kept constant (0.79). The pressure drop was calculated using the hydraulic diameter (full line), the analytic solution of Hagen–Poiseuille flow for a rectangular channel (broken line) and the same for an equivalent cross-sectional area (dotted line). We saw an expected impact of the flow rate. An increase or decrease in the rate led to an increase or decrease respectively in the TEP and TSEP of the flow-field pattern. The variation between the flows was large, four orders of magnitude for the values chosen, since the entropy production is scaled with the flow to the power of 2. Fig. 2 shows furthermore the change of the TEP and TSEP with different ways to compute the pressure drop. Generally speaking, the use of the hydraulic diameter gave the highest entropy production (pressure drop values), whereas the method which uses a diameter calculated from the equivalent cross-sectional area gave the lowest values. The difference between the choice of areas was not large, however, considering the variation in the gas flow rates.

Fig. 2 TEP as a function of channel width for a = 0.79, at a flow rate equivalent to current densities of 10 000 A m⁻² (upper region in the graph) and 100 A m⁻² (lower region in the graph). The pressure drop was calculated with the three different cases, namely Case D_h (full line), Case Equivalent A (dotted line) and Case Rectangular (broken line).

The width of the channels (Fig. 2 shows the channel width of generation 0 in mm) had a big impact on the TEP (and TSEP) of the flow-field pattern, especially at the initial increase of w₀. The higher w₀ was, the lower was the decrease in entropy production. This could be explained by a much lower pressure drop at higher channel widths (at constant channel depth). The results will be compared to full hydrodynamic 3D simulations below.

5.2 Murray's law and entropy production

The impact of the width scaling parameter a on the TEP and TSEP of the flow-field pattern (Study 2) is shown in Fig. 3. The figure shows the TEP and TSEP as a function of the width of generation 0 for different values of a, at a flow rate which corresponds to 10 000 A m⁻², where the pressure drop was calculated with the hydraulic diameter D_h. We saw that the value of a obtained from Murray's law (a = 0.79) did not give the lowest total entropy production (total specific entropy production) of the flow-field pattern, independent of the way we computed the pressure drop and the flow rate. An increase in a led to a decrease of the entropy production values. This was caused by the faster decline in the pressure drop with a higher a, see Fig. 4 (Study 3). If a > 0.79 (Fig. 4 bottom left) the pressure along the branches became a convex function, whereas if a < 0.79 (Fig. 4 top left) the pressure dropped in a concave way. Using Murray's law led (Fig. 4 top right) to a close-to-linear pressure variation. The same effects could also be found for a variation in w₀. Gheorghiu et al.¹⁴ described that if the pattern is scaled with Murray's law, a Lagrange multiplier optimization applied on the entropy production equation will show that the minimum entropy production will be reached when the pressure drop is constant along all branches. This was no longer the case here because we were setting the lengths manually to obtain a square area-filling pattern, contrary to the work of Gheorghiu et al. where the lengths were scaled with the assumption of a volume filling pattern (k = 3).

Fig. 3 TEP and TSEP as a function of channel width for different values of a. The flow rate was equivalent to a current density of 10 000 A m⁻². The pressure drop was calculated with Case D_h.

Fig. 4 Pressure along the branches of different generation levels at flow rates equivalent to current densities of 10 000 A m⁻² (black lines), 5000 A m⁻² (red lines) and 1000 A m⁻² (green lines) for a = 0.7 (top left), a = 0.79 (top right) and a = 1.1 (bottom left) and w₀ = 1 mm, where the pressure drop was calculated using Case D_h (full line), Case Equivalent A (dotted line) and Case Rectangular (broken line).

5.3 Flow rate distribution

As a consistency check, the overall outlet flow rate (summed flow rates at each branch outlet in the fourth generation) was compared to the set inlet flow rate (Study 4). The deviation was maximum 0.45% and minimum 0.19%. The differences were caused by numerical errors and we assumed that the results are therefore acceptable. The results for the analysis of the outlet flow rates can be found in Table 4. Furthermore, Fig. 5 shows the asymmetric effects caused by the branching in the tree-shaped pattern. As it can be seen, the difference was in most cases negligible, except for the cases with w₀ = 5 mm, a = 0.9 and a = 1. For the case of a = 0.9, differences up to 2.5% could be seen, but only at flow rates equivalent to 10⁴ A m⁻². The explanation could be the asymmetry in the system: the branches did not have enough time to fully develop the hydrodynamic flow, leading to an asymmetric flow. For the case of a = 1, variations up to 8% could be seen even at lower flow rates, making this geometry less suitable for use in fuel cells due to non-uniform fuel distribution. Also Fan et al.¹⁸ documented non-uniform flow rates, and proposed that a good flow distribution could be achieved at certain flow rates. In their study, however, only one geometry was analysed. Ramos et al.¹⁹ also worked on this type of investigation with similar results.

Table 4 Maximum differences in % between the highest flow rate at the outlet of the branches in the fourth generation and the actual flow rate at those outlets. The numerical errors in the numbers were within 0.5%

	10 000 A m^−2	5000 A m^−2	1000 A m^−2	100 A m^−2
w 0 = 1 mm, a = 0.79	0.02	0.03	0.04	0.08
w 0 = 1 mm, a = 0.9	0.11	0.20	0.27	0.53
w 0 = 1 mm, a = 1	0.23	0.36	0.62	0.84
w 0 = 1.5 mm, a = 0.79	0.03	0.01	0.01	0.01
w 0 = 1.5 mm, a = 0.9	0.16	0.09	0.14	0.36
w 0 = 1.5 mm, a = 1	0.37	0.14	0.17	0.46
w 0 = 2 mm, a = 0.79	0.24	0.17	0.31	0.59
w 0 = 2 mm, a = 0.9	0.29	0.03	0.04	0.09
w 0 = 2 mm, a = 1	0.77	0.10	0.04	0.09
w 0 = 2.5 mm, a = 0.79	0.20	0.02	0.03	0.06
w 0 = 2.5 mm, a = 0.9	0.61	0.06	0.07	0.14
w 0 = 2.5 mm, a = 1	0.75	0.11	0.18	0.34
w 0 = 5 mm, a = 0.79	0.28	0.14	0.40	0.53
w 0 = 5 mm, a = 0.9	2.54	0.41	0.70	0.90
w 0 = 5 mm, a = 1	3.78	3.82	7.77	8.53

---

Fig. 5 Flow velocity distribution in the geometry with w₀ = 1 mm and a = 0.79, at a flow rate equivalent to 10 000 A m⁻².

5.4 TEP and TSEP from 3D calculations

The results of the TEP and TSEP computations are shown in Fig. 6 (Study 5). The same conclusion could be drawn from this as we did in the 1D calculations. An increase in a and w₀ led to a lower TEP and TSEP of the flow-field pattern. Scaling the width according to Murray's law did not give the minimum entropy production. A quantitative comparison between the 3D and 1D results showed that the entropy production values, both the TEP and TSEP, were on the same order of magnitude. The maximum difference between the 1D and 3D calculations was dependent on the pressure drop calculation method. For the most suitable one, we had a maximum deviation of around 10%, whereas the worst one gave differences up to 61%. This emphasized the importance of the selection of the pressure drop calculation method. Also, the same trends could be observed.

Fig. 6 TEP and TSEP as a function of channel width for all simulated geometries for flow rates equivalent to current densities of 10 000 and 1000 A m⁻² for a = 0.79 (solid lines), a = 0.90 (dashed lines) and a = 1.00 (dotted lines). The two sets of three lines at the top are the TEP and TSEP respectively for 10 000 A m⁻², whereas the lower two sets of three lines are the TSEP and TEP respectively for 100 A m⁻².

5.5 1D- and 3D pressure drop calculations compared

The pressure, plotted along the center of the flow channels, for the 1 mm and 2.5 mm geometries at different width scaling parameters and a flow rate equivalent to 10⁴ A m⁻² is presented in Fig. 7. Here the 3D results are compared to the ones from the 1D calculations (Study 6). Again there was a lower entropy production at higher a values than the one given by Murray's law, which could be explained by looking at the pressure drop (Fig. 7). The pressure drop was nearly linear with a = 0.79, whereas there was a non-linear pressure drop with higher scaling parameters, leading to a higher gradient. In this sense, we reproduced the results of Gheorghiu et al.¹⁴

Fig. 7 Pressure plotted along the centerline of the channels for different a parameters, at a flow rate equivalent to a current density of 10 000 A m⁻², the pressure drop calculated with the three different cases and for the geometry with w₀ = 1 mm (top) and w₀ = 2.5 mm (bottom).

The pressure peaks appear only in the 3D simulations, which have a continuously connected flow channel. The peaks are due to the branching. In the bifurcations, flow is hampered and that increases their hydraulic resistance and, therefore, produces higher pressure drops to maintain the constant flow rate. 1D calculations based upon analytic equations do not produce such peaks. Fig. 7 gives a comparison of this branching effect on the pressure drop. The pressure peaks may be computed from the 1D model if the known approximate hydraulic formulas for the inlet flows, and flows in curved or T-shaped geometries, are used. Since CFD computations are quite fast and reliable, we have not reproduced the distributions with pressure peaks in the 1D models with approximate formulas. However, except for the first 2 branches, there was only a negligible pressure drop added to the system from this branching. Even though Ramos et al.¹⁹ used a slightly different geometry than here, the pressure drop values computed in the 3D simulations were of the same order of magnitude. The use of the hydraulic diameter overestimated the pressure drop in all simulated geometries. For the w₀ = 1 mm geometry, the pressure drop in generation 0 and 1 could be well estimated with the analytic solution of Hagen–Poiseuille for rectangular channels. However, the results started to deviate at higher generation levels. For the other geometries (e.g.Fig. 7 bottom), the most accurate way to estimate the pressure drop at higher generation levels was the method of equivalent cross-sectional area. However, the interesting part is that even though the pressure drop of the branches itself deviated from the method using the analytic solution, the overall pressure drop of the complete pattern (pressure at x = 0) was quite accurate up to a width of generation 0 of 2 mm. It seemed that the more rectangular a channel got, the more difficult it was to estimate the pressure drop in the channels with 1D calculations.

5.6 Peclet number

As discussed in Section 4.1, the Peclet number (eqn (18)) could be used as some form of design criterion for the channel geometries (Study 7). Trogadas et al.¹¹ approached this problem by looking at the number of generation levels needed to achieve a Peclet number lower than one, using the thickness of the PTL as the characteristic length. In our case, we have a different approach. We were looking at two different characteristic lengths. One was the hydraulic diameter of a rectangular channel, and the other one was the width of the last branch. By doing so we were able to find the Peclet number directly at the outlet of the flow pattern, independent of the thickness of the PTL. It was then possible to conclude on the status of the flow at the end of the tree-shaped pattern. The Peclet number at the last branch (here generation 4) was calculated for both mentioned characteristic lengths, a varying width between 1 and 5 mm, a width scaling parameter a of 0.79, 0.9 and 1, and a current density between 10 A m⁻² and 10⁴ A m⁻². We chose a diffusivity coefficient of 3.5 × 10⁻⁵ m s⁻² for oxygen in water vapour according to Trogadas et al. The results can be seen in Fig. 8.

The graphs illustrate the following points: both characteristic lengths give nearly the same results. At the last branch, the surface plot showing the Peclet number showed non-linear behavior, but neither the width of the channel nor the width scaling parameter had any significant influence on the Peclet number. The biggest impact came from the flow rate. For the case where the Peclet number was calculated with the width of the last branch as the characteristic length, the Peclet number was influenced by neither the width scaling parameter nor the width itself. This was due to the fact that the widths in eqn (18) cancelled each other out (if L = w_4,i). It only depended on the depth of the channel, which was constant in our case.

This means that it will be difficult to use the Peclet number as a design criterion for the dimension/shape of the flow pattern. The huge influence of the flow rate within a typical regime of flow variation prevented this. It was still possible to use the Peclet number of a given channel width and width-scaling parameter to find a current density range where the Peclet number was smaller than 1. This will then enable us to obtain a diffusion dominated flow at the last branches. An increased number of generations will make the condition more likely.

6 Discussion

We know already from the experiments of Trogadas et al.¹¹ that a tree-shaped field gives a better PEMFC performance than serpentine fields. The present flow-field gave also a uniform distribution of reactants, but at an entropy production that was lower than the one in the field predicted by Murray. The entropy production was reduced from branch to branch. The work pointed at further improvements, in terms of geometric choices, and possibilities to keep the Peclet number at a reasonably low value.

We have seen that much of the 3D-behaviour could be captured in a simpler quasi 1D-model. This has the interesting aspect that a model of the total cell obtains a good quasi 1D-representation. A simplified, but realistic, fuel cell model can, therefore, be obtained by combining the present 1D tree calculation with 1D-calculations, see e.g.ref. 8 and 21. This allows for simple tests of, say, the impact of boundary conditions. The value of the flow rate is more important for the concrete result than any approximation used to compute the hydraulic diameter. All this is possible due to the uniform flow distribution and because the inlet pattern is only in contact with the MEA over the 2^j_f branches of the last generation. The pattern used by Trogadas et al.¹¹ was similar to ours, as it was also inspired by a natural design (the lung). While their design was a 3D-design perpendicular to the MEA, ours had the flow-field tree in-plane with respect to the MEA. Our design can be machined with conventional milling techniques, where the outlet branches are located on a different layer. This may facilitate production and lead to decreased production costs for the flow-field plates.

A quasi-3D pattern allows for easy adjustment of the area filling property, through an increase of the generation levels. The pressure drop, as shown in Fig. 4, will not increase significantly when the generation level increases, due to the ever-decreasing gradient when a > 0.79. This creates room for the adjustment of generation levels and area filling properties. The pattern of Trogadas et al.¹¹ does not have this flexibility.

The computed pressure drop across our flow-field could be compared to measured results for the conventional serpentine or parallel patterns. In order to do so, we used the results from 3D simulations of Su et al.²² With similar operating conditions to ours, we found that our tree-shaped field had a lower pressure drop than the serpentine field. Depending on the chosen width scaling parameter a, the difference could be as big as one order of magnitude. It has to be noted, however, that the pressure drop across the tree-shaped pattern will slightly increase if a more space-filling pattern is used. On the other hand, the pressure drop was bigger than the one across the simulated parallel pattern (up to 1 order of magnitude). To see the full impact on the performance of our flow-field on the fuel cell, experiments need to be conducted, and compared to industry-standard flow-field patterns. Another point of interest in connection with the tree-shaped pattern is the use of different branching rules. For example, different branching angles or numbers of branches per generation level (e.g. 4 instead of 2, to create some form of “H” branching) can be introduced. Effects on the uniformity of fuel and entropy production need to be investigated and compared with the presented flow-field. The fractal FFP has been designed for uniform distribution of reactant gases over the PEM, with a constant density of outlet channels. This implies that the total entropy production in a fuel cell can be lowered by optimizing the flow-field separately from the MEA. This means that to find the minimum entropy production of the fuel cell, we have to minimize the total entropy production both in the MEA and the FFP. In this case, the flow-field optimization would take place in the x–y plane, whereas the MEA optimization would be over the z-direction due to the uniform fuel distribution. This separation of problems was performed already in the original article of Kjelstrup et al.,¹³ which also Cho et al.²³ draw inspiration from. The optimal catalytic layer did not depend (much) on the flow-field geometry, at least when the water accumulation had no impact on the hydraulic resistivity. According to Cho et al.²³ the multilayer 3D FFP may lead to high water accumulation. Thick 3D FFPs are not beneficial for FCs because of their high mass and cost. The present fractal flow-field can be milled into a flow-field plate of standard thickness and hopefully give an essentially lighter and cheaper alternative.

In a complete fuel cell system, there are many other different sources of entropy production. These sources and values for different current densities can be found in the paper of Sauermoser et al.²¹ which describes the 1D modelling of a PEMFC with non-equilibrium thermodynamics. The results show the main contributions to the entropy production, namely the electrode overpotentials, ohmic resistances and losses due to concentration-overpotentials. A general aspect of the above study should finally be pointed out. The flow-field evaluated here is by no means restricted to PEMFCs. Also, other fuel cells may benefit from a similar structure. Supply systems like this are common in nature, and may thus have more applications also in technology, say in catalysis.

7 Conclusion and perspectives

We have documented a tree-shaped flow-field based on Murray's law's constraints, but with different geometric variables, for the supply of gases to a PEMFC membrane. The gases were delivered at the membrane at uniform conditions, but with smaller entropy production than Murray's tree. Entropy production was due to viscous dissipation only. A quasi-1D model was used and found to represent a 3D model within a maximum deviation of 10% for the most suitable pressure drop calculation method and most flow magnitudes relevant to PEMFCs. This applied when the channel cross-section was rectangular, not circular. Flow channels which had a high width to depth ratio produce a small asymmetric behaviour, which explained the deviations.

We have thus established a basis for a 1D analysis of a whole single cell PEMFC, to be expanded to non-isothermal conditions in the future. 3D computations are time-consuming and can only be conducted for a restricted number of geometries. The 1D-results give an opportunity for fast quantitative estimations on different fractal-type geometries.

The main difference between this design and that of Trogadas et al. is that Trogadas et al. used a 3D structure, whereas we propose a two-layer flow-field plate with 2D patterns of constant channel depth. We hope that the new design may help reduce the weight and costs of the flow-field plate. Experimental proof beyond that offered by Trogadas et al.¹¹ remains to be obtained. Additional sources of entropy production from thermal, electrical and chemical sources exist in the MEA.²¹ Studies indicate that the present design allows us to optimize these sources independently of the flow-field optimization.¹³

Conflicts of interest

There are no conflicts of interest to declare.

Appendix 1: filling properties of self-similar networks

West et al.²⁴ required that the spherical volume around one branch was equal to the sum of the volumes of all spheres around the sub-branches. This is expressed in eqn (19). The conditions follow from the need to accommodate the entire volume flow in the network.²⁴ Each sphere has as the diameter the length of the corresponding branch.Here V_{spherical,i,j} is the spherical volume around one branch i at level j and V_spherical is the volume around the branch at level 0, which is constant.Rearrangement leads to eqn (21):Due to the premise used here that each branch splits into two branches, the ratio between N_j and N_j+1 equals 1/2. After reformulating eqn (21), the length scales as follows:Eqn (22) is similar to eqn (2), apart from the length scaling parameter k of 3. It is easy to show that the relation changes when the constraint changes. For constant occupation in terms of areas we obtainAny value for k can be used between constraints given by 3 or 2-dimensions, imply a particular space-filling structure. When k = 2, the flow-field is optimized to deliver the media to a (quasi-) two-dimensional area. For k = 3 it is optimised for a 3-dimensional one.

Fig. 8 Peclet number plotted over different flow rates equivalent to a certain current density in A m⁻² and channel width of generation 0 for different width scaling parameter a. Top: Calculated with the hydraulic diameter D_h as characteristic length L. Bottom: Calculated with w₀ as characteristic length L.

Acknowledgements

We gratefully acknowledge the support from NTNU in Trondheim and the Research Council of Norway through its Centre of Excellence funding scheme with Project No. 262644 (PoreLab). The simulations were performed on resources provided by UNINETT Sigma2 – the National Infrastructure for High Performance Computing and Data Storage in Norway.

Notes and references

1. C. Spiegel, Designing and building fuel cells, McGraw-Hill, New York, 2007.
2. P. J. Hamilton and B. G. Pollet, Fuel Cells, 2010, 10, 489–509.
3. H. Kahraman and M. F. Orhan, Energy Convers. Manage., 2017, 133, 363–384.
4. H. W. Wu, Appl. Energy, 2016, 165, 81–106.
5. H. Liu, P. Li, D. Juarez-Robles, K. Wang and A. Hernandez-Guerrero, Front. Energy Res., 2014, 2, 2.
6. T. Monsaf, B. M. Hocine, S. Youcef and M. Abdallah, Int. J. Hydrogen Energy, 2017, 42, 1237–1251.
7. O. Burheim, P. J. Vie, S. Moller-Holst, J. Pharoah and S. Kjelstrup, Electrochim. Acta, 2010, 55, 935–942.
8. S. Kjelstrup and D. Bedeaux, Non-equilibrium Thermodynamics Of Heterogeneous Systems, World Scientific, 2008.
9. M. LaBarbera, Science, 1990, 249, 992–1000.
10. N. Kizilova, Computational Science and Its Applications – ICCSA 2004, Berlin, Heidelberg, 2004, pp. 476–485.
11. P. Trogadas, J. I. S. Cho, T. P. Neville, J. Marquis, B. Wu, D. J. L. Brett and M.-O. Coppens, Energy Environ. Sci., 2018, 11, 136–143.
12. S. Kjelstrup, D. Bedeaux, E. Johannessen and J. Gross, Non-equilibrium thermodynamics for engineers, World Scientific, 2nd edn, 2017.
13. S. Kjelstrup, M. O. Coppens, J. G. Pharoah and P. Pfeifer, Energy Fuels, 2010, 24, 5097–5108.
14. S. Gheorghiu, S. Kjelstrup, P. Pfeifer and M.-O. Coppens, Fractals in Biology and Medicine, Basel, 2005, pp. 31–42.
15. C. D. Murray, Proc. Natl. Acad. Sci. U. S. A., 1926, 12, 299–304.
16. A. Arvay, J. French, J. C. Wang, X. H. Peng and A. M. Kannan, Nature inspired flow field designs for proton exchange membrane fuel cell, 2013.
17. F. M. White, Fluid Mechanics, McGraw-Hill, 7th edn, 2011.
18. Z. Fan, X. Zhou, L. Luo and W. Yuan, Exp. Therm. Fluid Sci., 2008, 33, 77–83.
19. B. Ramos-Alvarado, A. Hernandez-Guerrero, F. Elizalde-Blancas and M. W. Ellis, Int. J. Hydrogen Energy, 2011, 36, 12965–12976.
20. C. Hou, S. Gheorghiu, M.-O. Coppens, V. H. Huxley and P. Pfeifer, Fractals in Biology and Medicine, Basel, 2005, pp. 17–30.
21. M. Sauermoser, G. Fossati, N. Kizilova and S. Kjelstrup, Electrochem. Soc. Trans., 2019, 92, 279–292.
22. A. Su, Y. C. Chiu and F. B. Weng, Int. J. Energy Res., 2005, 29, 409–425.
23. J. Cho, T. Neville, P. Trogadas, Q. Meyer, Y. Wu, R. Ziesche, P. Boillat, M. Cochet, V. Manzi-Orezzoli, P. Shearing, D. Brett and M.-O. Coppens, Energy, 2019, 170, 14–21.
24. B. J. West, G. B. Brown and J. H. Enquist, Science, 1997, 276, 122–126.

---