Air-forced flow in proton exchange membrane fuel cells: Calculation of fan-induced friction in open-cathode conduits with virtual roughness

processes Article Air-Forced Flow in Proton Exchange Membrane Fuel Cells: Calculation of Fan-Induced Friction in Open-Cathode Conduits with Virtual Roughness Dejan Brki´c 1,2,*,†and Pavel Praks 1,*,† 1IT4Innovations, VŠB—Technical University of Ostrava, 708 00 Ostrava, Czech Republic 2Faculty of Electronic Engineering, University of Niš, 18000 Niš, Serbia *Correspondence: [email protected] or [email protected] or [email protected] (D.B.); [email protected] (P.P.) †Both authors contributed equally to this article. Received: 8 April 2020; Accepted: 3 June 2020; Published: 11 June 2020   Abstract: 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%. Keywords: Colebrook equation; fuel cells; flow friction factor; open-cathode; pressure drop; symbolic regression; numerically stabile solution; roughness 1. 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. 1.1. 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]. 1 qλT(p) =−2·log10           2.51 Re(p)·1 qλT(p) + ε(p) 3.71            (1) where: λT(p)—turbulent Darcy flow friction factor for pipes (dimensionless) Processes 2020,8, 686; doi:10.3390/pr8060686 www.mdpi.com/journal/processes Processes 2020,8, 686 2 of 12 Re(p)—Reynolds number (dimensionless)—the same definition as for fuel cells ε(p)—relative roughness of inner pipe surface (dimensionless) log10—logarithmic function with base 10 p—index related to pipes In Moody’s diagram, for the turbulent regime the Reynolds number Re(p) goes from around 2320 to 10 8 while the relative roughness of inner pipe surface ε(p) 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 λL(p) does not depend of roughness ε(p) while the Reynolds number Re(p) goes from 0 to around 2320; the formula for laminar Darcy flow friction factor for pipe flow, λL(p)= 64 /Re(p) 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 Re(p) and on the relative roughness of inner pipe surface ε(p) . As can be seen from the Moody diagram [ 4 ], for the same values of the Reynolds number Re(p) the turbulent flow friction factor λT(p) is higher in the pipes with higher relative roughness ε(p) , 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 λT(p) 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]. 1.2. 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 hand-held 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 – 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 ε(FC)= 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 ε(FC)=0.03086 is from [9]. 1 qλT(FC) =−10·log10           65.6 Re(FC)·1 qλT(FC) + ε(FC) 0.1415596            (2) where: λT(FC)—turbulent Darcy flow friction factor for fuel cells (dimensionless) Re(FC)—Reynolds number (dimensionless)—the same definition as for pipes Processes 2020,8, 686 3 of 12 ε(FC)—virtual relative roughness of fuel cell (dimensionless) log10—logarithmic function with base 10 FC—index related to Fuel Cells 2. Proposed Model Proton exchange membrane fuel cells (PEMFCs) transform chemical energy from electrochemical reaction of hydrogen and oxygen to electrical energy [ 22 – 25 ]. Here we analyze fan-induced 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 h/H 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 Re(FC) 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 Re(FC)= 500. In our numerical experiments, we use Re(FC) 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, λL(p)= 64 /Re(p) is used whereas the transition from laminar to turbulent flow is around Re(p)≈ 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 Re(FC)≈ 500, as explained in [ 8 ]. Therefore, for airflow through the cathode side of the observed fuel cells, the flow friction factor λL(FC) consists of two clearly defined types of flow: 1. laminar flow λL(FC) that depends both on the Reynolds number Re(FC) and on the geometry of conduits; height h and width H of the mesh of conduits that forms a mesh of cathodic air channels, and 2. turbulent flow λT(FC) is solely the function of the Reynolds number Re(FC) for the case from the experiment of Barreras et al. [ 8 ] (in general also on virtual roughness [ 9 ], which is in this case ε(FC)=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]. 2.1. 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): 1 qλT(FC) =−10·log10           65.6 Re(FC)·1 qλT(FC) +0.218           (3) where: λT(FC)—turbulent Darcy flow friction factor for fuel cells (dimensionless) Processes 2020,8, 686 4 of 12 Re(FC)—Reynolds number (dimensionless) – the same definition as for pipes log10—logarithmic function with base 10 FC—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]: 1 qλT(FC) =−14.17 +5·          −2·log10           2.51 Re(FC)·1 qλT(FC) + ε(FC) 3.71                       (4) where: λT(FC)—turbulent Darcy flow friction factor for fuel cells (dimensionless) Re(FC)—Reynolds number (dimensionless)—the same definition as for pipes ε(FC)—virtual relative roughness of fuel cell (dimensionless) log10—logarithmic function with base 10 FC—index related to Fuel Cells For Equation (4), virtual roughness can be recalculated based on the Colebrook equation as ε(FC)= 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 ε(FC) is not a real measurable physical characteristic of the surface of the used material for conduits (on the contrary for pipes ε(p) can be measured or at least estimated accurately [32–37]). Both, Equations (3) and (4) are numerically unstable for Re(FC)< 575, which can be a critical problem knowing that turbulent zone starts for Re(FC)> 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 z=eW(z)·W(z) . However, such transformation for the Colebrook equation for pipes contains a fast-growing term ex 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 ex 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]. 1 √λT(FC) =a(FC)·W(ex(FC))−Re(FC) b(FC)· ε(FC) c(FC) x(FC)=lnRe(FC) a(FC)·b(FC)+Re(FC) a(FC)·b(FC)· ε(FC) c(FC) a(FC)=10 ln(10) b(FC)=65.6 ε(FC) c(FC)=0.218 →c(FC)=0.1415596                            (5) where: λT(FC)—turbulent Darcy flow friction factor for fuel cells (dimensionless) Re(FC)—Reynolds number (dimensionless)—the same definition as for pipes ε(FC)—virtual relative roughness of fuel cell (dimensionless) a(FC),b(FC),c(FC)—constants x(FC)—variable log10—logarithmic function with base 10 ln—natural logarithm e—exponential function W—Lambert function Processes 2020,8, 686 5 of 12 FC—index related to Fuel Cells The parameter c(FC)for fuel cell is c(FC)=0.1415596. 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): 1 √λT(FC) =10 ln(10)·hln∆(FC)+W(ex(FC))−x(FC)i 1 √λT(FC) =1 a(FC)·hB(FC)+ωx(FC)−x(FC)i x(FC)=B(FC)+0.218·∆(FC) B(FC)=ln∆(FC)=lnRe(FC)−5.652138 ∆(FC)=Re(FC)·a(FC) 65.6 =Re(FC) 284.9 1 a(FC)=10 ln(10)≈4.343                                (6) where: λT(FC)—turbulent Darcy flow friction factor for fuel cells (dimensionless) Re(FC)—Reynolds number (dimensionless)—the same definition as for pipes a(FC)—constant x(FC),∆(FC),B(FC)—variable ln—natural logarithm W—Lambert function ω—Wright function FC—index related to Fuel Cells However, symbolic regression applied on the explicit formulation, Equation (6), which involves W(ex(p))−x(p)=ωx(p)−x(p) gives very simple, but still accurate results in case of pipe flow [6,7] ([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): W(ex(FC))−x(FC)=ωx(FC)−x(FC)≈26.723 lnRe(FC) 284.9 +0.218·Re(FC) 284.9 +6.2611 −3.6795 (7) where: Re(FC)—Reynolds number (dimensionless)—the same definition as for pipes x(FC)—variable ln—natural logarithm W—Lambert function ω—Wright function FC—index related to Fuel Cells To avoid repetitive computations, parameters ∆(FC) and B(FC) are introduced, in Equation (8). Both symbolic regression analyses were performed in Eureqa, a commercial software tool, which automates the process of model building and interpretation [46,47]. 1 √λT(FC) =4.343·B(FC)+26.723 B(FC)+0.218·∆(FC)+6.2611 −3.6795 ∆(FC)=Re(FC) 284.9 B(FC)=ln∆(FC)                (8) where: λT(FC)—turbulent Darcy flow friction factor for fuel cells (dimensionless) Re(FC)—Reynolds number (dimensionless)—the same definition as for pipes ∆(FC),B(FC)—variable Processes 2020,8, 686 6 of 12 ln—natural logarithm FC—index related to Fuel Cells 2.2. Unified Model Although the expression for laminar flow through pipes is λL(p)= 64 /Re(p) , for fuel cells it is different, as given in Equation (9) [8]: λL(FC)=58.91 +50.66·e−3.4·h H Re(FC) (9) where: λL(FC)—laminar Darcy flow friction factor for fuel cells (dimensionless) h H—channel depth/channel width used only in laminar flow (dimensionless) e—exponential function FC—index related to Fuel Cells Values of h/Hare from 0.83 to 2.5. The experiment [ 8 ] shows that air flow through the cathode side of air-forced open-cathode PEMFCs are (1) laminar for the lower values of the Reynolds number, Re(FC)< 500 and (2) turbulent for the higher values, 500 < Re(FC) <4000, where the Reynolds number is in hydraulics a very well-known dimensionless parameter that is used as a criterion for foreseeing flow patterns in a fluid’s behavior (defined in the same way for air flow through pipes and here discussed air flow through fuel cells). The dimensionless Darcy’s unified flow friction factor λ(FC) , is the function of the switching function y , the laminar flow friction λL(FC) , and the turbulent flow friction λT(FC) . The unified coherent flow friction model that covers both laminar and turbulent zones is set by Equation (10) [30]: λ(FC)=y·λT(FC)+(1−y)·λL(FC)(10) where: λ(FC)—unified Darcy flow friction factor for fuel cells (dimensionless) λT(FC)—turbulent Darcy flow friction factor for fuel cells (dimensionless) λL(FC)—laminar Darcy flow friction factor for fuel cells (dimensionless) y—switching function FC—index related to Fuel Cells The novel switching function yis given in Equation (11): y=Re(FC) Re(FC)+e558−Re(FC)(11) where: Re(FC)—Reynolds number (dimensionless)—the same definition as for pipes y—switching function e—exponential function FC—index related to Fuel Cells The switching function was obtained by symbolic regression using HeuristicLab [ 47 ] and it is given in Figure 1. Processes 2020,8, 686 7 of 12 Processes 2020, 8, x FOR PEER REVIEW 7 of 12 Figure 1. Switching function given in Equation (11). 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). 𝜆(𝐹𝐶)=𝑦·𝜆𝑇(𝐹𝐶)+(1−𝑦)·𝜆𝐿(𝐹𝐶) 𝜆𝐿(𝐹𝐶)=58.91+50.66·𝑒−3.4·ℎ 𝐻 𝑅𝑒(𝐹𝐶) 1 √𝜆𝑇(𝐹𝐶)=4.343·(𝐵(𝐹𝐶)+26.723 𝑥(𝐹𝐶)+6.2611−3.6795) 𝛥(𝐹𝐶)=𝑅𝑒(𝐹𝐶) 284.9 𝐵(𝐹𝐶)=𝑙𝑛(𝛥(𝐹𝐶)) 𝑥(𝐹𝐶)=𝐵(𝐹𝐶)+0.218·𝛥(𝐹𝐶) 𝑦= 𝑅𝑒(𝐹𝐶) 𝑅𝑒(𝐹𝐶)+𝑒558−𝑅𝑒(𝐹𝐶) } (12) where: 𝜆(𝐹𝐶)—unified Darcy flow friction factor for fuel cells (dimensionless) 𝜆𝑇(𝐹𝐶)—turbulent Darcy flow friction factor for fuel cells (dimensionless) 𝜆𝐿(𝐹𝐶)—laminar Darcy flow friction factor for fuel cells (dimensionless) 𝑅𝑒(𝐹𝐶)—Reynolds number (dimensionless)—the same definition as for pipes ℎ 𝐻—channel depth/channel width used only in laminar flow (dimensionless) 𝑥(𝐹𝐶), 𝛥(𝐹𝐶), 𝐵(𝐹𝐶)—variables 𝑦—switching function 𝑒—exponential function 𝑙𝑛—natural logarithm 𝐹𝐶—index related to Fuel Cells For 216 = 65536 Sobol Quasi Monte-Carlo pairs [48], which cover 𝑅𝑒FC = 50–4100 and for ℎ/𝐻 from 0.83 to 2.45, the maximal relative error of the final calculated flow friction factor 𝜆FC 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. 3. Software Code and Measurement of Execution Speed Figure 1. Switching function given in Equation (11). The laminar flow friction λL(FC) depends on the Reynolds number Re(FC) , but also on geometry of conduits, while the turbulent flow friction λT(FC) depends only on the Reynolds number Re(FC) . In the case of fuel cells, both coefficients are empirical. In addition, the switching function y contains the exponential function, (the similar situation is for calculation of λL(FC)as already explained). To avoid numerical instability, it is recommended to use the explicit approximation which gives the following unified formula in Equation (12). λ(FC)=y·λT(FC)+(1−y)·λL(FC) λL(FC)=58.91+50.66·e−3.4·h H Re(FC) 1 √λT(FC) =4.343·B(FC)+26.723 x(FC)+6.2611 −3.6795 ∆(FC)=Re(FC) 284.9 B(FC)=ln∆(FC) x(FC)=B(FC)+0.218·∆(FC) y=Re(FC) Re(FC)+e558−Re(FC)                                          (12) where: λ(FC)—unified Darcy flow friction factor for fuel cells (dimensionless) λT(FC)—turbulent Darcy flow friction factor for fuel cells (dimensionless) λL(FC)—laminar Darcy flow friction factor for fuel cells (dimensionless) Re(FC)—Reynolds number (dimensionless)—the same definition as for pipes h H—channel depth/channel width used only in laminar flow (dimensionless) x(FC),∆(FC),B(FC)—variables y—switching function e—exponential function ln—natural logarithm FC—index related to Fuel Cells For 2 16 =65536 Sobol Quasi Monte-Carlo pairs [ 48 ], which cover ReFC =50–4100 and for h/H from 0.83 to 2.45, the maximal relative error of the final calculated flow friction factor λFC 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. 3. 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 2should be followed. Processes 2020,8, 686 8 of 12 Processes 2020, 8, x FOR PEER REVIEW 8 of 12 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. Figure 2. Algorithm for solving Equation (12). 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), as 𝐿𝑇=4.343·(𝑏+26.723./(𝑏+ 0.218·𝑎+6.2611)−3.6795). The MATLAB code also works in GNU Octave, but it can be easily translated in any programming language. The final unified fuel cell model given by Equation (12) is coded in MATLAB as: 𝑓𝑢𝑛𝑐𝑡𝑖𝑜𝑛 𝐿=𝑃𝐸𝑀𝐹𝐶𝑠(𝑅,ℎ) 𝑎=𝑙𝑜𝑔(10)/10; 𝑑=𝑅·𝑎/65.6; 𝐵=𝑙𝑜𝑔(𝑑); 𝑥=𝐵+0.218·𝑑; 𝐿𝑇=4.343·(𝐵+26.723./(𝑥+6.2611)−3.6795); 𝐿𝑇=1./𝐿𝑇.^2; 𝐿𝐿=(58.91+50.66·𝑒𝑥𝑝(−3.4·ℎ))./𝑅; 𝑦=𝑅./(𝑅+𝑒𝑥𝑝(558−𝑅)); 𝐿=𝑦.·𝐿𝑇+(1−𝑦).∗𝐿𝐿; ⌉ ⌉ ⌉ ⌉ ⌉ ⌉ ⌉ ⌉ ⌉ ⌉ in this code: Output parameters of the function: -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 216 = 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 Figure 2. Algorithm for solving Equation (12). The turbulent flow friction λT(FC)∼LT (with the intermediate step through 1/qλT(FC)∼LT ) can be expressed using “wrightOmega” function as LT =(B+wrightOmega(x)−x)/a , but it can be executed faster using approximation as given in Equation (8), as LT =4.343 ∗(b+26.723./(b+0.218 ∗a+6.2611)−3.6795) . The MATLAB code also works in GNU Octave, but it can be easily translated in any programming language. The final unified fuel cell model given by Equation (12) is coded in MATLAB as: function L =PEMFCs(R,h) a=log(10)/10; d=R∗a/65.6; B=log(d); x=B+0.218 ∗d; LT =4.343 ∗(B+26.723./(x+6.2611)−3.6795); LT =1./LT.ˆ2; LL =(58.91 +50.66 ∗exp(−3.4 ∗h))./R; y=R./(R+exp(558 −R)); L=y.∗LT +(1−y).∗LL;                                               in this code: Output parameters of the function: -flow friction factor λ(FC)is given as L Input parameters of the function: -Reynolds number Re(FC)is given as R(from the interval 50 <Re(FC)<4100) -channel depth/channel width h/His given as h(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 Re(FC) =50–4100 and for h/H 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. Processes 2020,8, 686 9 of 12 4. 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 future research, further analyses and experiments are foreseen to show if and how this value of the virtual roughness ε(FC) of fan-induced friction in a mesh of conduits for air-forced flow in the open-cathode PEMFCs can be changed and how it depends on fuel cells parameters (type of fuel cells, size, geometry of channels, etc.). The study should be extended to cover other types of fuel cell, to include water and heat management, etc. [52–60]. Author Contributions: D.B. had the initial idea for this article and he dealt with physical interpretation of the problem. P.P. conducted detailed numerical experiments (included symbolic regression) following initial inputs of D.B. D.B. wrote a draft of this article. Both authors contributed equally to this article. All authors have read and agreed to the published version of the manuscript. 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: For pipes: λT(p)turbulent Darcy flow friction factor for pipes (dimensionless) λL(p)turbulent Darcy flow friction factor for pipes (dimensionless) Re(p) Reynolds number (dimensionless)—the same definition as for fuel cells ε(p)relative roughness of inner pipe surface (dimensionless) pindex related to pipes For Fuel Cells: λ(FC)unified Darcy flow friction factor for fuel cells (dimensionless) λT(FC)turbulent Darcy flow friction factor for fuel cells (dimensionless) λL(FC)laminar Darcy flow friction factor for fuel cells (dimensionless) Re(FC)Reynolds number (dimensionless)—the same definition as for pipes ε(FC)virtual relative roughness of fuel cell (dimensionless) h H channel depth/channel width used only in laminar flow (dimensionless) x(FC),∆(FC),B(FC)variables a(FC),b(FC),c(FC)constants FC index related to Fuel Cells Functions: log10 logarithmic function with base 10 ln natural logarithm eexponential function WLambert function ωWright function yswitching function