Analytic Mathematical Model for Concentration Profile in a Parallel Plate Electrolyzer with Variable Current Density

The parallel plate electrolyzer is of interest in electrochemical applications. However, the modelling of an electrolyzer often involves numerical integration or use of infinite series that is computationally expensive. Using some simplifying assumptions, a novel analytical solution to the governing equation could be obtained using a similarity transformation. To obtain the solution, a mathematical limit evaluation is required rather than a simple substitution, resulting in a novel transcendental function that describes the constants of integration. However, the constants of integration could be excellently approximated by a rational function, with the coefficients obtainable by solving equations for special cases or linear regression. The quadratic and quartic regression were found to give a better fit. It is also demonstrated that the general solution reduces to constant current density and surface concentration under special cases. Using this approach, a simple mathematical model was developed to approximate the solutions that often require cumbersome numerical methods such as the finite element method. © The Author(s) 2017. Published by ECS. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 License (CC BY, http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse of the work in any medium, provided the original work is properly cited. [DOI: 10.1149/2.0561711jes] All rights reserved.

The parallel plate electrolyzer is widely used for industrial electrochemical application due to its ease of construction and continuous operation. Because of this, the concentration profile for the boundary layer on the electrode plates of an electrolyzer is of theoretical interest in the design and optimization of a parallel plate electrolyzer. To study the boundary layer concentration profile, the parallel plate electrolyzer problem can be considered as a flat plate version of a Graetz and Leveque convection diffusion problem, where there is solute convection parallel to surface and a solute diffusion perpendicular to surface. While the plug flow (constant velocity profile) convection diffusion problem is well understood, the Levesque flow (linearly proportional velocity profile to distance from channel wall) convection diffusion problem, which is a better representation of actual flow profile, is not well studied. The standard closed form solution involves simply the uniform surface concentration: 1 On the other hand, constant current density and variable current density (or concentration) do not have closed form solutions readily available. While an excellent solution for surface reaction 2 was previously provided, the equations were cumbersome, perhaps due to the complex nature of the system defined. The model was also entirely limited to a first order reaction. For most other non-uniform current density and surface concentrations, the modelling of electrolyzer often involves cumbersome numerical simulation, such as infinite series [3][4][5][6][7][8][9][10][11] or numerical integration, [12][13][14][15][16][17][18] both of which are computationally expensive and difficult to compare with literature. The previous mentioned work discussed the possibility of a similarity solution 19,20 but did not solve it, other than 2,3,21,22 for constant concentration, however they all provided a limiting case such as the surface concentration being zero as in Ref. 22.
This electrolyzer modelling problem is encountered during glycerol carbonate electrosynthesis scale up. In this paper, it will be used as basis of the problem formulation. However, the solution derived is applicable to many reactions by changing the numerical values of the constants involved.

Problem Formulation
Consider a parallel plate electrolyzer with plane Poiseuille flow in Figure 1. The fluid consists of potassium carbonate, K 2 C O 3 , dissolved in a mixture of glycerol water, entering the parallel plate with initial carbonate concentration C 0 . The electrolyzer consists of a sealed rectangular compartment with 2 parallel plates, which are connected to power supply by electrical contact at the inlet. On the cathode, carbonate is reduced to some carbonate radical, while on the anode, glycerol is oxidized to some glycerol radical. While the actual species and mechanism of electrochemical glycerol carbonate formation is not well understood, the number of electron per cathode reactions, z 1 , can be easily determined by experiment, for any electrochemical reaction.
For Poiseuille flow, it is assumed that the fluid is Newtonian, meaning the viscosity to be constant with respect to the flow velocity. The concentration profile of the solute, which in this case is carbonate ion, under a plane Poiseuille flow, can be described by the steady state equation of continuity: Note the absence of z since it is assumed that the concentration profile is uniform in z-direction. It is also assumed the properties of the electrolyte, the diffusivity is constant. Effects like bubbling have been negligible for this pioneering work. 23 Expanding the terms: To obtain the analytic solution, multiple assumptions have been made. The first assumption is the diffusion in x-direction is negligible (plug flow assumption: Another assumption is the quadratic term to be negligible, which holds true for conventional electrolyzers with width subsequently greater than molecule size 6U Near the boundary layer, however, both the quadratic flow and the diffusion in x-direction are negligible. Such that the governing equation can be approximated as: To obtain the boundary conditions, it is assumed that electrolysis rate is insufficient to consume the concentration in bulk. This gives a Dirichlet boundary condition, Boundary Condition 1: While this boundary condition may sound counterintuitive because the bulk concentration may not be the initial concentration, the value C 0 is the first principle limit at infinite distance from the surface because a finite electrolysis rate could not deplete solute infinitely far away from electrode surface.
As such, the model would still be valid if the electrolysis rate is sufficient to affect the bulk concentration at finite width, with the bulk concentration in actual case of finite distance evaluated numerically by plugging the value of the electrolyzer width d into the concentration profile formula.
In addition, while the main aim of the reactor is to obtain high conversion rates, the low conversion assumption is usually applicable due to Faraday's law of electrolysis, which limits the amount of product to a factor of 1/96500 from the electric current through the electrode.
Next, the concentration at x = 0 is simply the inlet concentration C 0 . This gives another Dirichlet boundary condition, Boundary Condition 2: Boundary Condition 3, however, is a Neumann boundary condition obtained by a decaying surface flux at y = 0: Where K is a constant is described by [8] Note that in general, K can be a negative number for consumption flux or a positive number for generation flux. The rationale of Boundary Condition 3 is a result of obtaining the flux on cathode from charge balance: Where i(x) is the current density distribution over the electrode. While the actual current density distribution is not known normally as a hot topic of research, 13 it has been well-modeled in prior art that for constant surface concentration, the current density is inversely proportional to electrode length, such that it decreases with distance with an order of 1 3 : 24 As will be shown in discussion in later section, it must be emphasized that the order of 1 3 only applies when the surface concentration is constant. For a more general case in actual experiment other than constant surface concentration, the current may decay with different order other than factor of 1 3 . [25][26][27] To generalize, the decay order is set to a value p between 0 and 1, or p ∈ [0, 1): Where k and p are constants that can be fitted empirically. The x − p model is generally a good fit as shown in Figure 2 even though the actual current density would not go to infinity as x → 0 because mass transport cannot be infinitely fast. This x − p model is used as a compromise between simplicity and accuracy; it is simple enough to result in closed form formula for the concentration profile, yet accurate enough to represent the first principle properties of current density distribution. Because current density measurement can be obtained from electrical setup, the current density data is much easier and cheaper to obtain than multidimensional concentration data from in situ spectroscopic method. As a result, the analysis offers a way to obtain concentration profile from preliminary empirical data that can be obtained easier than actual concentration profile from experiment. While this solution may not be as accurate approximation as numerical method, it offers a way to obtain great approximation of concentration profile at substantially lower computational cost, as opposed to more complex numerical solutions that would require numerical integration point by point with some sophisticated computational facility such as MATLAB or Mathematica.
As will be shown in later section, the current density model is an interpolation between constant surface concentration and constant current density, where the two ideal extremes would not require the empirical current density distribution data. For the non-ideal cases however, this model offers a way to model the non-ideal concentration profile using convenient preliminary data, in this case simplified into simple parameter p and K . Such p can be measured easily using electrical setup without the need of any concentration profile (spectroelectrochemistry) such as UV-Vis or IR spectroscopy, both of which are much more expensive and trickier to implement than potentiostat alone.
The lower limit of p ≥ 0 arises from a constraint by first principle, the current density decreases with increasing axial distance because as axial distance increases, current need to travel longer, resulting in higher resistance. Such constraint means that the current density decays to zero or stays constant (for negligible plate resistance under operating range) as x approaches infinity because the ohmic resistance is lowest and the reactant concentration is highest at the electrolyzer inlet: [12] Because of this constraint, p must be great than or equal to 0 because the current density must either decay or stay constant with x.
The upper limit of p < 1 arises from another constraint that the average current densityī must from the model must be finite value and must agree with actual measurement from total current I and electrode length L. By the model, the average current densityī can be evaluated by the improper integral of current density from 0 to L: [13] For this improper integral to be finite value, p must be less than 1.
Both of these constraints of finite density at x → ∞ and x = 0 give rise to the choice of limited range of p, p ∈ [0, 1). Finally, plugging in the current density model into surface flux gives the Boundary Condition 3 as desired:

Derivation of Solution
To solve the governing equation with the boundary conditions, rearrangement followed by similarity transformation is used. The constants of integration are then derived by evaluation of the defined mathematical limit.
Preliminary transformation.-The following substitution is applied to reset the concentration: After applying the preliminary transformation: And the boundary conditions become, as shown in Figure 1: Similarity transformation.-The dilation similarity transformation is then applied:Ĉ = x m f (η) [21] η = y x n [22] Note that m and n are not necessarily equal.
To demonstrate the similarity transformation, the derivatives are evaluated: [25] Substituting and simplifying into the PDE gives: For the PDE to be reduced to ODE, the powers of x on both sides must be equal to cancel out. This gives: Substituting the similarity transformation into Boundary Condition 1: Substituting the similarity transformation into Boundary Condition 2: It can be seen that Boundary Condition 1 and Boundary Condition 2 merge to the same Boundary Condition A, which is a hallmark of the similarity transformation method.
Note also that x m f (∞) = 0 must hold true for any x, so f (∞) must be zero. This gives Boundary Condition A: Now, substituting the similarity transformation into Boundary Condition 3: In addition, the powers on both sides of BC must also be equal to cancel out: Such that the values of m can be obtained in terms of p: This gives Boundary Condition B: Substituting and simplifying gives the ODE: Solution to the ordinary differential equation.-The ODE is of the form: The general solution of this class of ODE using the Wolfram Computational Engine 28 is: [39] Where 1 F 1 (A; B; z) is the "Kummer's Confluent Hypergeometric Function", while C 1 and C 2 are constants of integration to be determined.
Substituting back A and B, the solution to the ordinary differential equation in this case is: However, there are some challenges to obtaining the constants of integration, because both functions approaches ∞ as η → ∞: Substituting these into the general function gives: Such that by naïve substitution, one might think that C 1 = C 2 = 0, but this would give a trivial solution! To work around this difficulty, a mathematical limit has to be taken: Interestingly, it is found that this newly defined mathematical limit exists under the range of p ∈ (0, 1). By inspection on Wolfram Alpha Computational Engine, it is found that the multiplication factor M 9 can be removed from the limit evaluation as in: This gives a convenient way to evaluate the limit on Wolfram Computational Engine: Where g( p) is a function of p. As will be shown in discussion section, while the exact values of g( p) is a transcendental function of p, it can be well approximated by a rational function. Before applying Boundary Condition B, the derivative was evaluated using Wolfram Alpha Computational Engine: [47] Applying Boundary Condition B: Given the properties of all hypergeometric functions becomes 1 as the argument z of the hypergeometric function goes to zero: ) unless CC License in place (see abstract Substituting into and simplifying gives C 2 as: Substituting back for C 1 : Substituting back into the ODE solution gives: [52] Because it has been shown that n = 1 3 , the general solution to the ODE is then: [53] To obtain a more convenient equation for the concentration profile, a conversion is performed:

Discussion
The special values of g( p) at p = 0, 1 3 , 2 3 are first presented. g( p) is then obtained for many values of p, and multiple models are developed to approximate g( p) to any required accuracy. Table I.

Values of g( p ).-The number g( p) is found to have certain values under special cases as tabulated in
Obtaining the numerical values were evaluated at different pusing Wolfram Computational Engine and plotted: It is found that g( p) shows asymptotic increase to positive infinity at values close to p = 1, as shown in Figure 3: This suggests that g( p) could be modeled involving a function having a factor of 1 1− p , which is also a function that shows asymptotic increase to positive infinity at values close to p = 1: Defining a transformed integration constant function h( p), such that: g( p) plot is now transformed to h( p) plot as follows: From Figure 4, it can be seen that h( p) can be well fitted to a quadratic and a quartic curve: Even assuming a quadratic model, the g( p) could be obtained without regression using the special values, as demonstrated: Solving this set of linear equations of 3 equations and 3 unknowns, it is obtained that: However, as shall be seen later, such analytic fit is not a perfect fit, better fits of g( p) can be obtained by quadratic and quartic regressions with the regression coefficients tabulated in Table II fit of each model, the percentage error from actual Wolfram output is first defined as: It can be seen from Figure 5 that all fits result in excellent accuracy of less than 0.5% at low p, with the error increasing drastically near p → 1 − . Such error is simply of numerical computation nature because each fit inevitably has deviation from actual value, which end up getting magnified by the multiplication factor 1 1− p that goes to infinity near p → 1 − . The accuracy and computation complexity increases together in the following order: quartic model, quadratic model and analytic fit. Note, however that the maximum error near p → 1 − for the least accurate model is still within 5%, indicating the analytic fit is still a reasonable approximation near p → 1 − .
The error for quartic fit did not rise drastically because it is a great fit that compared to the simpler models of analytic and quadratic fit (both has only fit of 2 orders). As a result, the error did not get magnified as much. Nevertheless, such increases in accuracy comes with the higher computational price of more complex quartic model.
The resulting coefficients for g( p) are tabulated as follows: Special cases of p.-Constant flux: p = 0 .-At p = 0, the hypergeometric functions reduce to: 29,30 The surface concentration is: The surface flux is: The surface concentration changes with factor of x − 1 3 as the electrode length while the surface flux changes with factor of x − 2 3 as the electrode length.
It can be seen that p = 2 3 leads to non-uniform Dirichlet and Neumann condition because both the surface concentration and the surface flux are changing with axial distance. This is an interesting case of both surface concentration and flux varying, in a way such that the hypergeometric solution reduces to Bessel functions. While this case of p = 2 3 is rare in the real world, it represents an intermediate case between the two ideal extremes between constant current density and constant surface concentration. As a result, this intermediate case is a good approximation for mediumsized electrolyzer where both the surface concentration and the current density are changing, especially without the value of p known yet from current density data.
The cases are finally tabulated in Table III for comparison.

Conclusions
A closed analytical solution for the concentration profile of a parallel plate electrolyzer has been successfully developed. The closed y 3 x )] + C 0 x )] y 3 x )] + C 0 form solution has advantages in that it is very short, and could be computed much easier than using numerical methods and infinite series. It is also much easier to compose the simulation code since the number of equations substantially reduced. While this mathematical model involves approximations, it fills the gaps of closed form solutions. It is also a more robust model that generalizes the cases of constant flux and constant surface concentration.

Future Work
The closed form solution developed is only valid for a big enough electrolyzer width such that the quadratic term becomes negligible. For cases where this is not true, as in microchannel electrolyzer, the quadratic term has to be taken account. While a dilation similarity solution of the form C = x m f (η), η = y/x n does not exist for this case, 35,36 the solution may be approximated by combining this linear term solution and a quadratic solution 37 to obtain a composite solution by regular perturbation.

My
∂C Where the perturbation variable is defined as The similarity solution is a part of general analytical method to solve partial differential equations, known as Lie Group Symmetry. 38,39 A number of sources [40][41][42][43][44] suggest that a novel solution with fewer simplification assumptions may be developed. However, it is unknown whether the symmetry solution could be represented as a known function 45 such as the closed form solution developed in this paper. A parallel plate electrolyzer for electrochemical glycerol carbonate synthesis will be constructed to validate the mathematical models with empirical data.

Acknowledgments
This work is partially supported by the Connaught Innovation Award and the Engineered Nickel Catalysts for Electrochemical Clean Energy project administered from Queen's University and supported by grant No. RGPNM 477963-2015 under the Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Frontiers Program. a 1 , a 2 , a 3 , a 4