Air ‐ Forced Flow in Proton Exchange Membrane Fuel Cells: Calculation of Fan ‐ Induced Friction in Open ‐ Cathode Conduits with Virtual Roughness

: Measurements of pressure drop during experiments with fan ‐ induced air flow in the open ‐ cathode proton exchange membrane fuel cells (PEMFCs) show that flow friction in its open ‐ cathode side follows logarithmic law similar to Colebrook’s model for flow through pipes. The stable symbolic regression model for both laminar and turbulent flow presented in this article correlates air flow and pressure drop as a function of the variable flow friction factor which further depends on the Reynolds number and the virtual roughness. To follow the measured data, virtual inner roughness related to the mesh of conduits of fuel cell used in the mentioned experiment is 0.03086, whereas for pipes, real physical roughness of their inner pipe surface goes practically from 0 to 0.05. Numerical experiments indicate that the novel approximation of the Wright ‐ω function reduced the computational time from half of a minute to fragments of a second. The relative error of the estimated friction flow factor is less than 0.5%.


Introduction
Flow friction is a complicated physical phenomenon and it is not constant but depends on flow rate and pressure drop. Because of its complexity, equations which describe the flow friction are mostly empirical [1]. The most used empirical formulation for turbulent pipe flow is given by Colebrook's equation [2]. Here we will evaluate flow friction factor caused by air flow in the cathodic side of the observed PEMFCsand we will provide accurate and consistent solution.

Colebrook Equation for Pipe Flow Friction
The standard Colebrook's friction equation for turbulent pipe flow [2]; Equation (1), follows logarithmic law and is based on an experiment performed by Colebrook and White in the 1930s [3], while its graphical interpretation was given by Moody in 1944 [4].
where: -turbulent Darcy flow friction factor for pipes (dimensionless) -Reynolds number (dimensionless)-the same definition as for fuel cells -relative roughness of inner pipe surface (dimensionless) -logarithmic function with base 10 -index related to pipes In Moody's diagram, for the turbulent regime the Reynolds number goes from around 2320 to 10 8 while the relative roughness of inner pipe surface from 0 for smooth surfaces to about 0.05 for very rough surfaces (on the other hand, flow friction for laminar regime for pipe flow does not depend of roughness while the Reynolds number goes from 0 to around 2320; the formula for laminar Darcy flow friction factor for pipe flow, 64 ⁄ is not empirical but theoretically founded).
Colebrook and White experimented with airflow through pipes of different roughness of inner surface [3]. They used a set of pipes from which one left with a smooth inner surface, while others were covered with sand of different size of grains. For each pipe, one uniform grain size with glue as adhesive material was used. Thus, the pipes in the experiment were gradually smooth to the fully rough. The experiment revealed that the turbulent flow friction depends on the Reynolds number and on the relative roughness of inner pipe surface . As can be seen from the Moody diagram [4], for the same values of the Reynolds number the turbulent flow friction factor is higher in the pipes with higher relative roughness , where subsequently for the same flow, the corresponding pressure drop is higher, too.
The Colebrook equation is implicitly given in respect to the unknown turbulent flow friction factor and it can be rearranged in explicit form only in terms of Lambert W-function [5] or its cognate Wright ω-function, and even then further evaluation can be only approximated. Very accurate approximate formulas for the Colebrook equation for pipe flow based on the Wright ωfunction are available [6,7].

Modified Colebrook Equation for Flow Friction
The Colebrook equation is empirical [3] and therefore possible modifications based on different conditions of flow can be done. We will evaluate here a modification for fuel cells [8,9]. Further, for example, US Bureau of Mines published a report in 1956 [10] that introduced a modified form of the Colebrook equation for gas flow where coefficient 2.51 was replaced with 2.825. Also, very recent modifications of a variety of empirical equations for pipe flow are available [11]. Here we analyze one modification for air flow friction in the open-cathode side of the observed type of PEMFCs [8], and subsequently, we give a stable and computationally efficient explicit solution which is valid in this case. Analogous analyses can be done for air flow for cooling of electronic products from handheld devices to supercomputers [12,13]. Here we discuss the only influence of hydraulic effects of flow through the open-cathode conduits of fuel cells while available literature [14][15][16][17][18][19][20][21] should be consulted for thermodynamic aspects.
The Colebrook equation can be rearranged following data obtained from the experiment by Barreras et al. with fuel cells [8]. Based on this data and following analogy with pipe flow, it is estimated that virtual roughness is fixed by value 0.03086. Therefore, the Colebrook equation can be rewritten for the observed fuel cell as Equation (2). We will further analyze this equation to provide its stable numerical solution. The value of virtual roughness 0.03086 is from [9].

Proposed Model
Proton exchange membrane fuel cells (PEMFCs) transform chemical energy from electrochemical reaction of hydrogen and oxygen to electrical energy [22][23][24][25]. Here we analyze faninduced air-forced flow, based on data from the experiment with pressure drop in the cathode side of air-forced open-cathode PEMFCs by Barreras et al. [8]. Their experiment with fuel cells can be compared with the experiment performed by Colebrook and White with air flow through pipes [3]. Barreras et al. [8] use three different cathode configurations with aspect ratios ℎ/ from 0.83 to 2.5 to form a mesh of cathodic channels to supply the fuel cell with enough air for cooling and with oxygen for a chemical reaction (roughness is real physical characteristic of pipe surface [26], but not of cathodic channels of fuel cells in terms of hydraulic performances).
Value of the Reynolds number during the experiment was from 45 to 4000, while as a difference from pipes, during flow through the observed fuel cell, the transition from laminar to turbulent flow occurred around 500. In our numerical experiments, we use between 50 and 4100.
As already mentioned, the original Colebrook equation is valid for turbulent flow of air, water, or natural gas through pipes. On the other hand, for laminar flow, 64 ⁄ is used whereas the transition from laminar to turbulent flow is around 2320. This transitional border at the Moody's plot [4] is sharp where the equivalent sharp transition from laminar to turbulent flow for the observed fuel cell starts at around 500, as explained in [8]. Therefore, for airflow through the cathode side of the observed fuel cells, the flow friction factor consists of two clearly defined types of flow: 1. laminar flow that depends both on the Reynolds number and on the geometry of conduits; height ℎ and width of the mesh of conduits that forms a mesh of cathodic air channels, and 2. turbulent flow is solely the function of the Reynolds number for the case from the experiment of Barreras et al. [8] (in general also on virtual roughness [9], which is in this case 0.03086).
For solving implicitly given equations, instead of iterative procedures [27,28], appropriate explicit approximations which are accurate but also computationally efficient can be used (review of appropriate explicit approximations for pipe flow friction is available in [29]). A computationally efficient and stabile unified equation for the observed type of fuel cells which is valid both for laminar and turbulent regime will be given in this article [30].

Turbulent Flow
In case of turbulent airflow during experiments with open-cathode PEMFCs, measurements show that pressure drop during turbulent flow at its cathode side follows logarithmic law, which form is comparable to the Colebrook's flow friction equation for pipe flow, but with different numerical values [8]. The flow friction related to air flow is given by Equation (3): where: -turbulent Darcy flow friction factor for fuel cells (dimensionless) -Reynolds number (dimensionless) -the same definition as for pipes -logarithmic function with base 10 -index related to Fuel Cells During turbulent flow, numerical values for the flow friction factor in pipe and fuel cells are different and that difference can go up to 60%. To make a direct connection between Equation (1) for pipe flow and Equation (3) for the observed fuel cell, Equation (4) can be used [25]: where: -turbulent Darcy flow friction factor for fuel cells (dimensionless) -Reynolds number (dimensionless)-the same definition as for pipes -virtual relative roughness of fuel cell (dimensionless) -logarithmic function with base 10 -index related to Fuel Cells For Equation (4), virtual roughness can be recalculated based on the Colebrook equation as 0.03086 for the observed fuel cell in the experiment [8]. This fuel cell was tested with three different cathode configurations [8]. As noted in [31], this roughness is not a real measurable physical characteristic of the surface of the used material for conduits (on the contrary for pipes can be measured or at least estimated accurately [32][33][34][35][36][37]).
Both, Equations (3) and (4) are numerically unstable for 575, which can be a critical problem knowing that turbulent zone starts for 500. However, the novel solution proposed in this article is numerically stable. Generally, implicitly given equations can be transformed in explicit form through the Lambert W-function [38,39]. The Lambert W-function [5] is defined as the multivalued function W that satisfies . However, such transformation for the Colebrook equation for pipes contains a fast-growing term and because of that, overflow error in computers is possible [40,41]. Happily, results with fuel cells show that the solution is not in the zone where is too big to be stored in registers of computers. The model for fuel cells is given in Equation (5), while the related model for pipe flow friction model can be seen in [42]. After procedures from [6,7,43], the following form for fuel cells expressed through the Lambert W-function and its cognate Wright ω-function is given in Equation (6) [44,45] confirm these results), but unfortunately these analytical formulas, which are optimized for pipes, cannot be directly applied on the fuel cell equation. Fortunately, symbolic regression gives also very promising results for fuel cells as given in Equation (7)

Unified Model
Although the expression for laminar flow through pipes is 64 ⁄ , for fuel cells it is different, as given in Equation (9)  -switching function -index related to Fuel Cells The novel switching function is given in Equation (11): (11) where: -Reynolds number (dimensionless)-the same definition as for pipes -switching function -exponential function -index related to Fuel Cells The switching function was obtained by symbolic regression using HeuristicLab [47] and it is given in Figure 1. The laminar flow friction depends on the Reynolds number , but also on geometry of conduits, while the turbulent flow friction depends only on the Reynolds number . In the case of fuel cells, both coefficients are empirical. In addition, the switching function contains the exponential function, (the similar situation is for calculation of as already explained). To avoid numerical instability, it is recommended to use the explicit approximation which gives the following unified formula in Equation (12) [48], which cover = 50-4100 and for ℎ/ from 0.83 to 2.45, the maximal relative error of the final calculated flow friction factor using Equation (12) is 0.46% compared with the original Equation (2). The accuracy and speed of execution are tested through the code given in the next section.

Software Code and Measurement of Execution Speed
The unified equation for laminar and turbulent fan-induced air flow through open-cathode side of the observed PEMFCs is given in Equation (12), where the algorithm from Figure 2 should be followed. The turbulent flow friction ~ (with the intermediate step through 1 ⁄~ ) can be expressed using "wrightOmega" function as ℎ / , but it can be executed faster using approximation as given in Equation (8) -flow friction factor is given as Input parameters of the function: -Reynolds number is given as (from the interval 50 < < 4100) -channel depth/channel width ℎ/ is given as ℎ (from 0.83 to 2.5) In MATLAB, symbol log() denotes the natural logarithm. Implicitly given equations can be solved using iterative procedures [49,50], but also using appropriate accurate but also computationally efficient explicit approximations especially developed for the observed purpose. Our numerical results show that the computationally efficient approximation does not contain the time-consuming evaluation of the Wright ω-function [6,7], but simple polynomials found by symbolic regression [51], which can be easily evaluated on computers. It is because simple functions such as +, −, * and / are executed directly in the CPU and hence they are very fast [43].
Using 2 16 = 65536 Sobol Quasi Monte-Carlo pairs [48], which cover = 50-4100 and for ℎ/ from 0.83 to 2.45, the evaluation of the MATLAB built-in "wrightOmega" function required 30.9 s, while our novel approximation required only 0.0022 s (our approximation is around fourteen thousand times faster). Consequently, numerical tests show that the novel approximation presented here is very suitable for modeling of fan-induced flow friction in a mesh with virtual roughness for air-forced flow in the open-cathode PEMFCs, as it is very fast and still very accurate.

Conclusions
This paper gives a novel numerically stable explicit solution for flow friction during airflow in cathode side of open-cathode PEMFCs. Symbolic regression is successfully used for obtaining a cheap but still accurate approximation of the Wright-ω function for fuel cells-based explicit friction model and also for approximations of the switching function, which is needed for the unified formulation of fuel cell friction model. Numerical experiments indicate that the novel approximation of the Wright-ω function reduced the computational time from half of a minute to fragments of a second. The relative error of the estimated friction flow factor is less than 0.5%.
In Funding: The work of D.B. has been supported by the Ministry of Education, Youth and Sports of the Czech Republic through the National Programme of Sustainability (NPS II) project "IT4Innovations excellence in science-LQ1602", while the work of P.P. has been partially supported by the Technology Agency of the Czech Republic under the project "Energy System for Grids" TK02030039. Additionally, the work of D.B. has been partially supported by the Ministry of Education, Science, and Technological Development of the Republic of Serbia.

Conflicts of Interest:
The authors declare no conflict of interest.

Notations
The following symbols are used in this article: