Electrochemical transport modelling and open-source simulation of pore-scale solid-liquid systems

The modelling of electrokinetic flows is a critical aspect spanning many industrial applications and research fields. This has introduced great demand in flexible numerical solvers to describe these flows. The underlying phenomena is microscopic, non-linear, and often involving multiple domains. Therefore often model assumptions and several numerical approximations are introduced to simplify the solution. In this work we present a multi-domain multi-species electrokinetic flow model including complex interface and bulk reactions. After a dimensional analysis and an overview of some limiting regimes, we present a set of general-purpose finite-volume solvers, based on OpenFOAM(R), capable of describing an arbitrary number of electrochemical species over multiple interacting (solid or fluid) domains. We provide a verification of the computational approach for several cases involving electrokinetic flows, reactions between species, and complex geometries. We first present three one-dimensional verification test-cases for single- and multi-domain cases and then show the capability of the solver to tackle two- and three-dimensional electrically driven flows and ionic transport in random porous structures. The purpose of this work is to lay the foundation of a general-purpose open-source flexible modelling tool for problems in electrochemistry and electrokinetics at different scales.

The involvement of electrical forces leads to a number of interesting electrokinetic phenomena.One such phenomena is electroosmosis [1], where an applied electric field induces fluid movement due to formation of electric double layers, see Fig. 1.Another interesting set of phenomena is known as induced-charge electrokinetics.Whilst similar to traditional electrokinetic phenomena, the difference comes from the double layer being induced by an applied electric field [36].The main difficulty of modelling electrokinetic flows stems from the microscopic scale of phenomena.For applications with large scale domains, solving at the microscopic scale can be computationally taxing.For systems with complex geometries, such as porous media, this difficulty is furthered.
To describe electrokinetic flows we require modelling the ion concentrations, electric potential and fluid velocity, with the Stokes-Nernst-Planck-Poisson (SPNP) model.In this work, for the velocity field we will, in fact, consider Stokes' flow as, at the micro-and meso-scale, the viscous forces far outweigh the inertial due to small velocities and length scale [35].The Stokes equation is coupled with the Poisson's equation for the electric potential, relating the electric fields composition to the variation in ion charge density.Finally these are both coupled with the Nernst-Planck equation to describe the ion transport.In using Nernst-Planck, we neglect any ion-ion interaction by assuming the ionic solution is sufficiently dilute.The Nernst-Planck ionic flux was first formulated in its steady form for a one-dimensional cylinder by Nernst [26].It was later extended by Planck [30] to a transient setting, furthermore introducing the continuity and Poisson's equations [22].In doing this the two paved the way in helping develop SPNP, providing a simple yet accurate description of electrokinetic flow for dilute solutions.
A common similarity in many applications containing electrokinetic flows is the involvement of multiple regions, such as fluid electrolyte and solid electrodes in batteries.These regions exchange ions with each other and in some instances contain chemical reactions exchanging mass between species.For batteries these reactions are a necessary process for operational use.However in examples such as chloride corrosion in reinforced concrete, it is a detriment [25].As such, modelling these reactions has just as much importance as the flow they reside in.
Whilst used extensively in a wide range of physical settings by the research community, SPNP does come with its own caveats.For one, SPNP neglects any ion-ion interactions that may occur by assuming a dilute solution.This may not hold true for solutions with many ionic species [32].Also, as SPNP is a continuum model, any steric effects are ignored.As such, many efforts have been made extending SPNP to be include other physical processes.To cover steric effects, free energy functionals using density functional theory (DFT) accounting for long-range Coulomb correlation and hard sphere (HS) interactions of ions [14,21,37] are formulated.Extensions to make Nernst-Planck more thermodynamically consistent under non-equilibrium thermodynamics [18,8] have also been proposed.To model non-ideal solutions, [34] proposes an added term to the Nernst-Planck ionic flux considering varying chemical activities solved by the Debye-Hückel model.
Whilst these extensions do further the physical realism of the original SPNP model, this often exacerbates other challenges of SPNP.One of the foremost being the non-linear coupling between fields.For example, in [34] the Debye-Hückel model equates the chemical activity of a species to the solutions ionic often used electro-neutrality approximation.For §4 we outline what is required to capture reactions at a multi-region interface given different restriction on ion movement.In §5 the implementation of our solvers for single and multi-region is discussed, as well as the iterative algorithm performed when introducing our reactive conditions.To verify accuracy of our solvers and reactive conditions, we provide necessary numerical examples in §6, with concluding remarks given in §7.

Stokes-Poisson-Nernst-Planck model
Here we discuss the equations that make up Stokes-Poisson-Nernst-Planck, modelling the advective, diffusive and electrostatic forces of an ionic solution.As many real-world applications of electrochemistry involve interacting solids and fluids we consider a multi-domain scenario of a whole domain Ω split, without loss of generality, into two sub-domains Ω f , a fluid, and Ω s , a solid, such that Ω = Ω f ∪ Ω s , and with Γ being the solid-fluid interface and ∂Ω = ∂Ω f ∪ ∂Ω s the external boundaries.More in general, in some applications and in our computational framework, we have allowed for an arbitrary number of solid and fluid regions, separated by different interfaces.We consider N ionic species, with concentrations and valencies c i and z i respectively, and i = 1, ..., N .To describe the ionic transport we must define equations for the electric field E, ion concentrations c i within Ω and fluid velocity u within Ω f .

Stokes' flow
Consider Ω f , with velocity profile u(x, t) governing the advective dynamics of the ions.Assume a negligible Reynolds numbers Re defined by the fluid density ρ f , characteristic velocity U , characteristic length scale L and dynamic viscosity µ such that viscous forces within Ω f outweigh the inertial.This common assumption in ionic transport [35,15,7] leads to linear Stokes flow.Furthermore we assume the fluid to be incompressible, i.e., where p = p(x, t), ρ el = ρ el (x, t) = F N i=1 z i c i and E are, respectively, the static fluid pressure, electric charge density and electric field respectively, and F is Faraday's constant.Compared to the standard Stokes equation, we have the presence of the body force term ρ el E, describing the Coulomb forces acted on the fluid by the ions [15,31].We neglect any magnetic contribution by assuming our ions move slowly such that E is irrotational, i.e. ∇ × E = 0.This body force term may be set to zero if the fluid is unaffected by E or the fluid is approximated as electrically-neutral, ρ el = 0.This will be further discussed in §3.The first coupling term, between the variables u, E and c i appears here, showing one of the significant difficulties of describing electrokinetic flows.The coupling terms (particularly if non-linear) often add significant numerical difficulties.First of all, they make the velocity field time-dependent, although.neglecting the time derivative, we assume instantaneous relaxation to equilibrium.For segregated approaches, the disparity of relaxation time scales between Stokes, Poisson and Nernst-Planck, leading to mixed parabolic-elliptic systems, can pose severe instability problems or slow convergence of the coupled system.

Poisson equation
To model the electric field E, neglecting magnetic forcses, We may then write E = −∇φ and focus on the electric potential φ = φ(x, t).Assume for each sub-domain their respective electric permittivity ε is spatially constant.From Maxwell's equations we obtain Poisson's electrostatic equation denoting variations in φ by changes in the charge density ρ el , where ε = ε s in the solid and ε = ε f in the fluid.Again we have direct coupling between our variables, here between c i and φ, although this time the former appears linearly in the source term.Like for Stokes flow, this equation depends on time only through the source/coupling terms, in particular the time-dependent ionic concentration c i (x, t).

Ionic transport
Assume the ionic fluid in Ω f is sufficiently dilute to ignore ion-ion interactions and diffusion is isotropic.Under these assumptions we may use the Nernst-Planck flux [15,30,22,10] as where denote R, T , D i,f and D i,s are, respectively, the ideal gas constant, absolute temperature and diffusion coefficient of species i in the fluid and in the solid.Taking (4) in conjunction with the continuity equation for mass conservation we arrive at the set of equations modelling transport of c i , As mentioned, we assume a dilute solution to ignore ion-ion interactions.This may not hold true in some cases.Alternatively the Stefan-Maxwell equations [32,15,16] may be used in lieu of Nernst-Planck.In short, Stefan-Maxwell balances the driving forces exerted on a species with the frictional forces between species.This introduces cross diffusivities D ij describing the drag between species i and j due to these frictional forces [15].The added complexity from such explicit coupling between species is difficult to describe, as the exact description of D ij is hard to determine [32].Even in situations where the ionic fluid is bordering on dilute, it is common for Nernst-Planck to still be used [33] due to its simplicity.
To close the system of equations listed here, we must include conditions along the boundaries of Ω f and Ω s .We briefly discuss some common choices below and their physical representation.We also require conditions at the interface Γ to describe interaction between sub-domains.This is discussed further below for c i and φ when considering reactions and (non-)conductive interfaces.

Boundary conditions
Below we list some relevant conditions for different physical situations often seen in electrokinetic problems: • Solid walls: Either a conductive or insulating wall.The no-penetration (u • n = 0, and zero normal stresses), where n is the outward unit normal to the boundary, or no-slip condition (u = 0) may be used.Ions may not pass through the wall so their fluxes are zero (j i • n = 0).If conductive we impose a fixed normal current density I w along the wall (i If insulating (non-conductive) we may simply set I w = 0.
• Permeable membrane and reactive boundaries: For a membrane we may fix the flux of each species (j i • n = j m ), non-zero for those capable of passing through, zero for who cannot.We will discuss the case of reactive boundaries later in §4.2.
• Inlet/Outlet: Often used to model an in or outflow of fluid or ions.For a momentum driven inlet, pressure or flow velocity may be fixed (u = u in , p = p in/out ).For ions we may apply fixed fluxes (j i •n = j in/out ), resulting in a fixed electric current (i•n = I in/out ).More complicated is to determine suitable boundary conditions for the potential.Typically either a Dirichlet φ = φ in/out or a Neumann condition can be imposed.The latter imposes the total electric current to be either equal to zero or a fixed value.
• Periodicity: when dealing with large (quasi-)periodic (or homogeneous) structures, it is often impossible to solve for the entire domain of interest.In these cases, smaller representative unit cells can be solved with (quasi-)periodic external boundary conditions [5] are imposed.In these cases, an additional driving forces (as a bulk source term or a modification of the periodic BC) need to be added.

Dimensional analysis
When dealing with systems of coupled transport equations it is useful to perform a dimensional analysis to better understand the relationships between the different transport phenomena and identify the limiting regimes and possible approximations.We denote dimensionless variables with a hat symbol, e.g x, and reference values with a bar unless stated otherwise.We use the reference values L, U , φ and ci for the length scale, velocity, electric potential and concentrations respectively.For pressure we take p = µU L as this the appropriate form when under Stokes flow.For time, we take the diffusive timescale t = L 2 D i .With these choices and defining C = 1 2 N i=1 z 2 i ci as the reference total ionic strength, we obtain the dimensionless system of equations Substituting these dimensionless variables into (1), (3) and (5) the dimensionless system of equations is The underlying dimensionless numbers may be found by dividing all other terms by the reference values of one term.For ( 7)- (9) we divide by z 2 i ci is the reference ionic strength needed to resolve the issue of not being able to factor out the reference concentrations ci : where we have defined four dimensionless numbers N , Pe, L and P, defined in table 1 Note that whilst (7) had four terms, we only arrive at two numbers due to our choice of reference time.The same can be said with (12) and reference pressure.1.: List of dimensionless numbers and their forms In table 1 we list the forms of these numbers in terms of the reference values where we denote e, k B and N A to be the elementary charge, Boltzmann's constant and Avogadro's number respectively.N represents the ratio of electrostatic over diffusive forces, with N 1 indicating diffusion is dominant.We also obtain the Péclet number Pe, i.e., the ratio of advective over diffusive phenomena and L = λ D L , with the dimensionless form of the Debye length λ D defined as and approximating the distance at which a charge's electrostatic effect persists.Typically L 1 stating how the Debye length is much smaller than the reference length L. Finally we have P denoting the ratio of viscous and electric forces upon our fluid.λ D is also the approximate width of a common electrokinetic phenomena known as the electric double layer (EDL), see Fig. 1, that forms on boundaries.The EDL consists first of ions adsorbed at the boundary, known as the Stern layer, and another of free ions moved by electrical attraction and diffusive motion by the Stern layer, deemed the diffuse layer.
For sufficiently thin EDLs there is a common model reduction known as the electro-neutrality assumption, often employed in development of macroscopic models [33,35,17].This reduction can in fact be determined through dimensional analysis and will be briefly discussed next.

Asymptotics and electro-neutrality
The electro-neutrality approximation states that for a sufficiently dilute electrolyte all charges of ionic species within the solution roughly cancel each other out, leaving the solution electrically neutral.Electroneutrality for a solution containing N species is defined as and is often used as model reduction when modelling ionic flows.Where ( 14) becomes invalid however is in the thin charged double layers, or EDL, mentioned above and often seen in real-world settings.
To get an understanding of where electro-neutrality comes from and how it relates to EDL formation, we perform an asymptotic analysis.Without any lack of generality, we limit here to a one-dimensional electrolyte and a binary electrolyte.In the numerical framework described above an arbitrary number of ionic species in three-dimensions can be considered.Let x ∈ Ξ where Ξ = [0, 1] is a binary electrolyte solution with ionic concentrations c 1 and c 2 of opposing valencies.For now we omit any boundary conditions, only requiring those chosen result in a boundary layer formation near x = 0. To arrive at (14) we start by considering the asymptotics of the outer (bulk) layer away from x = 0.For simplicity we only consider the leading order terms..

Outer (bulk) layer:
To arrive at an asymptotic leading order solution in the bulk layer of Ξ we take the following expansions of û, φ, p, ĉ1 and ĉ2 .These denote the fluid velocity, electric potential, pressure and ion concentrations respectively, all dimensionless.Expansions are taken in powers of = L 2 , the squared dimensionless Debye screening length since 1, or, λ D L: Substituting these expansions into the dimensionless transport equations ( 10)-( 12) for one-dimension, alongside the incompressibility condition, we arrive at By considering only the leading order terms i.e., terms of O( 1), these equations reduce to, Note how ( 24) is the electro-neutrality approximation mentioned before, for a binary solution, and a direct consequence of the asymptotics.This implies said electro-neutrality is only accurate up to leading order.What's more the electric body force of Stokes vanishes as a consequence of (24).To close the system an equation for the leading order potential φ (0) is required.This is possible by multiplying (23) by their respective valencies z i , summing over i and utilizing (24), resulting in which represents a steady equation for φ (0) .This can be interpreted as a modified Ohm's law with new conductivity accounting for all contributions to the total current: advective, diffusive and electrical, from left to right.Now that we have constructed the asymptotic solution, up to O(1), for the outer (bulk) layer of Ξ, we move on to determine the inner (boundary) layer asymptotics.As mentioned at the start, we state said inner layer formation is near x = 0 of Ξ.Much like the outer layer, we make no case of boundary conditions, only that a boundary layer near x = 0 forms as a result of them.
Inner solution: To construct a solution for the inner layer near x = 0 of Ξ we define the variable y to span the inner layer and be a 'fast' variable counterpart to x, changing more rapidly, By the definition of x ∈ Ξ y therefore has domain [0, ∞), as y = 0 for x = 0 and y → ∞ for x → 1, since 1.Such change of variable gives derivatives via chain rule as Substitution of the above derivatives into the original transport equations ( 10)-( 12) results in, where all variables u, p, φ and c i are in terms of y e.g., c i = c i (y, t).Just like the outer layer, we substitute the asymptotic expansions ( 15)-( 18) into the above equations to obtain Considering only the leading order terms of the four equations above: O( −1 ), O(1), O( −1 ) and O(1) terms respectively, we arrive at the following leading order set of equations for the inner layer: Like the outer layer equations, to close the system we need another equation, this time for p (0) .Considering the O( (1/2) ) terms of (36) we can retrieve such an equation as As mentioned earlier, when considering situations involving multiple sub-domains we must also have appropriate interface conditions.In many real applications of electrokinetic flows there are chemical reactions at interfaces between sub-domains, exchanging mass across the ionic species involved.In these cases, neglecting the EDL might lead to significant errors.In the next section we will consider such reactions occurring on our interface Γ, formulating appropriate conditions to capture them whilst retaining mass conservation.

Multi-domain formulation and reactions
Here we formulate conditions to model heterogeneous reactions which are crucial for many electrochemical and electrokinetic problems.We write a general reaction rate that ensures total mass conservation and apply it to form reactive interface conditions, where we consider the scenarios of species that exist in the whole domain (unrestricted ) or only in a specific sub-domain (restricted ).We then discuss conditions on φ when the interface is conductive or non-conductive.

Reaction model
Consider a general elementary reaction transferring mass between reactants i = J, ..., K and products i = K + 1, ..., M with exchanged molar masses B i and valencies z i .We denote ν i > 0 to be the stoichiometric coefficients determining the number of moles of species i is lost or gained and n the number of released (n < 0) or absorbed (n > 0) electrons, alongside the electron mass e − .The mass balance of the reaction reads: To formulate (43) as conditions on Γ we first determine the rate r i at which each species in our reaction is exchanging mass.More complex reactions involving several intermediate reactions can be decomposed into elementary steps, i.e a set of elementary reactions.Typically the overall reaction rate is then given by the rate of the slowest elementary reaction [6].One option to determine the rate is the law of mass action, also known as the rate law [2,6].Given an ideal solution and reaction involving chemical species [j] with stoichiometric coefficients ν j , at dynamic equilibrium where the net reaction rate (density) r with units mol m 2 s , is given by: Here k f , k r are the rate constants for the forward and reverse reaction and can be empirically modelled by the Arrhenius equation [31,6].Note that r is always positive, so we cannot simply take r i = r as we must allow r i < 0 for reactant species.Instead, we define the coefficient α i as and, multiplying by r and ν i , we obtain: For species not involved in the reaction we set r i = 0. Assuming a closed reactive system, i.e no trace ion species, we can represent the mass conservation as the balance across reaction rates r i : We weight by the molar masses m i to convert from moles (a non-conserved quantity between species) to grams.We use the rate law as the example here as it is valid for many reactions [6].It is also the more general form for the commonly used Butler-Volmer equation for faradaic reactions, which uses the energy dependence of k f and k r to give explicit electric potential dependency.It is important to notice that the units of measure of the reaction rate are [mol/s] for bulk reaction and [m mol/s] for surface reactions.Therefore, for the case of linear reactions, the reaction constants k i can be either [1/s] or [m/s].

Interface conditions
The reaction models above can be applied in the bulk or on a interface.Here we apply them on Γ for φ and c i given conductive or non-conductive interface and surface reactions respectively.We employ general reaction rates r i to balance the ionic fluxes j i through Γ with mass exchanged by the reaction.For φ we use conservation of charge to find conditions on the current passing through Γ when acting conductively or non-conductively.

Interface conditions for ions concentration
Here we outline reactive conditions along Γ to model (43).We use a set of general, possibly non-linear, reaction rates r i for i = J, . . ., M of ion species involved in (43).We make no assumptions on the form of r i , other than making sure (49) holds true.For species not involved we simply take flux continuity.the results is jump conditions in ionic flux between Ω s and Ω f , equal to their respective reaction rates r i : Here we denote j i • n Γ f to be the normal flux evaluated at Γ from Ω f 's side.Alongside this we assume continuity of concentrations to provide our second condition: So, for all species involved in the reaction (43) we have the difference in ionic flux from Ω s and Ω f to be the rate r i of the respective species.In some scenarios one or more of the ion species may be restricted to reside in a single sub-domain of Ω.This can however be easily modelled by simply modifying the conditions of those restricted species, setting both j i and c i to be zero in the inaccessible domains.

Interface conditions for the potential
We consider conditions for the electric potential φ given two situations.Suppose our interface Γ is acting as a conductor such that we see a electric current flowing through.Assume we are under the most general case where all ion species are unrestricted.Conservation of charge and Ohm's law then implies the change in current through Γ is proportional to the sum of changes in ionic fluxes: Here we denote σ to be the electric conductivity which is discontinuous across the interface.Alongside this we assume continuity of the potential: We see therefore that the conditions on the current are closely linked to the previous reactive conditions.As if no reaction occurs at Γ, we simply have continuity of current due to ionic flux continuity.If a reaction is present, then (50) tells us our jump in current is proportional to the sum of reaction rates r i alongside their respective charge numbers z i .
If Γ is non-conductive instead we fix separately the current on each side of Γ and no longer impose continuity:

Numerical implementation
As we have mentioned our goal is to construct simple, effective numerical solvers to the Stoke-Poisson-Nernst-Planck (SPNP) model at the pore scale.To do so, we employ the computational fluid dynamics (CFD) package OpenFOAM ® due to its open-source nature, large active community and robust handling of complex geometries.OpenFOAM ® is a finite volume library for general unstructured mesh.Coupled equations are solved iteratively in a segregated/splitted approach, solving sequentially each discretised equation.This removes the explicit coupling between equations, as well as the need of linearising multilinear terms (terms linear in each variable but where multiple variables appear) at the expense of having internal iterations.These would be unavoidable also if we adopted a monolithic approach due to the non-linearities of the model.Two separate solvers have been implemented.The first, pnpFoam, models electrokinetic flow of a single ionic fluid, modelled by SPNP.The second, pnpMultiFoam, is more generalised, modelling a general set of ionic fluids and solids following SPNP and diffusion respectively.Furthermore we developed the numerical counterpart, named mappedChemicalKinetics, to the reactive conditions in §4.2.1 for the case of a binary reaction.In this section we will describe in detail the structure of the solvers and the corresponding boundary conditions.

Single-and multi-domain solvers
OpenFOAM ® discretise each equation using a sparse matrix M with entries relating to the cell centres of the mesh.For example, for the fluid momentum the discrete form reads: where [u] is the vector of unknown velocities at cell centres and M the coefficient matrix scaling the effect of neighbouring cell velocities.Note the right hand terms are left as source terms and are computed explicitly and the gradients operators can be discretised with different schemes (details about the schemes will be presented in the results section).To solve for the pressure-velocity coupling we employ the PIMPLE algorithm with an extra term due to the electric body force.This works by decomposing M into its diagonal, A, and off-diagonal, H, parts M[u] = Au − H.This leads to the velocity correction equation: Interpolating u to the cell faces and taking the dot product with cell face area vectors S f leads to the flux U correction equation where subscript f denotes values at cell faces.Discretising the incompressibility condition ∇ • u = 0 gives ∇ • U = 0 which when applied to the flux correction equation forms the pressure correction equation: This is solved iteratively until convergence.In an external loop, momentum and pressure equations are coupled with the concentration and potential equations.When the mesh is highly skewed, or for complex discretisation schemes with implicit-explicit terms, additional iterations can be added for each single equation.The pseudo-code algorithm for our single ionic fluid solver pnpFoam is presented in Algorithm 1.

Algorithm 1: pnpFoam algorithm
Construct fields: u, p, φ, c i ; while t < end time do while !converged do Initialise u n+1 , p n+1 to previous values: u n+1 = u n , p n+1 = p n ; Solve momentum equation for u n+1 : −µ∇ 2 u n+1 + ρ el ∇φ = −∇p n+1 ; while Pressure correction do Solve pressure correction equation: ∇ The same method is employed within the algorithm of pnpMultiFoam, our solver for ionic transport over a general set of ionic fluids and solids, where each region is solved separately, and an additional outer loop is added to ensure the coupling between regions.The algorithm is presented in Algorithm 2. Both solvers here are presented for the case of non-electroneutral solution.In case the electroneutrality is assumed, the algorithm is slightly modified to solve a modified potential equation eq. ( 27) (with ionic conductivity instead of permittivity) and with the last species calculated to ensure electro-neutrality.

Boundary and interface conditions
To make use of the object orientation of OpenFOAM ® , all conditions are first reformulate as effective Robin conditions.We consider here, as an example, the inhomogeneous Robin BC for a variable c with coefficients D , K and F along a boundary Γ with normal n: In OpenFOAM ® the boundary values c| Γ and ∇ n c| Γ are approximated using the value of the cell centres with faces along Γ: Here c c and c f are the values of c at the cell centre and face respectively.The α's are the interpolation weights, and B is the inverse distance between the cell centre and the boundary.We can then rearrange (60) using ( 61)-( 62) to find the α values that allow us to approximate (60) using the cell centres c c : This forms the basis of the Robin BC implemented.Reformulating all other conditions into a form like (60) lets us reuse the same equations in (63) to approximate all other conditions.The interface conditions are implemented as a derived class named mappedChemicalKinetics.If we consider here the limiting case of two reacting species c s and c f , each restricted to their respective subdomains Ω s and Ω f respectively.Along the interface Γ connecting the two sub-domains is the reaction: The general reactive conditions eq. (50) in this case reduces to the following condition along Γ: where denote here Γ s to be evaluation at the Ω s side of Γ and n the unit normal of Γ facing into Ω f .Note the reaction rate r is kept general and not necessarily linear.Equation ( 65) is solved iteratively using the Newton-Raphson method, linearised about the previous solution c N = (c N s c N f ) .Here the previous solution could be the solution at the previous time step (if an explicit time stepping is chosen) or at the previous internal iteration (for fully implicit time stepping).This allows us to rewrite (65) into two decoupled effective Robin conditions for c N +1 s and c N +1 f that we solve separately.The coupling between the two sides of the interface conditions (and therefore the two domains) is achieved through the internal iterations but, for stiff reactions, additional sub-looping to update both boundary values whenever one of the two domains is solved for.In algorithm 3 we detail the mappedChemicalKinetics pseudo-code to solve the effective Robin conditions.As both resulting effective Robin conditions are solved in the same manner we only write here solving of the condition for c f .
More details about the linearisation can be found in appendix A. Alongside the reactive conditions we have also implemented a number of simpler conditions, such as the continuity of total fluxes, continuity of value or continuity of derivatives, often used within applications.Just as with the non-linear reactive condition, we rewrite all conditions into effective Robin conditions.

Numerical examples
Here we present four numerical examples of electrokinetic flows.To verify accuracy of results we compare with the spectral Matlab ® toolbox Chebfun [9] with machine-precision accuracy.The first case verifies Call mappedChemicalKinetics BC on the solid side; end Increment iterator: N → N + 1 end accuracy of the flow description given a single ion species using a pressure-driven infinite ion channel similar to [4].Next we verify accuracy when considering multiple ion species.Afterwards we verify the implemented reactive interface conditions counterpart to §4.2.1.The final case displays the capabilities of our solver(s), simulating ionic transport within a randomised solid-fluid porous medium.
To show spatial convergence of OpenFOAM ® results we use the following normalised L 2 error point norm.We use here the dummy variables v and v ref to denote OpenFOAM ® and Chebfun results respectively as an example:

SPNP in an infinite channel
To verify pnpFoam and pnpMultiFoam for a single ion species take Ω, of length L, with two boundaries Γ t and Γ b denoting the top and bottom channel walls respectively.Take boundaries Γ in and Γ out denoting the inlet and outlet of fluid in Ω.To mathematically describe this channel as infinite in length we take periodic conditions on Γ in and Γ out . to induce transport across the channel length.To allow electro-migration of c between the two channel walls, fix φ = 0.1 along Γ t and φ = 0 on Γ b .For c apply no-flux along Γ t and Γ b , denoting n to be the outward unit normal of either Γ t or Γ b .We consider the channel at steady state and, by assuming L is sufficiently large, take all derivatives in x to be zero, i.e. ∂a ∂x = 0 for a ∈ {u, φ, c, p}.Just like Poiseuille flow we assume a uni-directional velocity, such that u = u 1 0 .These assumptions in turn produce a one-dimensional reduce system of equations in y.
The SPNP system with a single ionic species in an infinite channel is considered with an added pressuredriven driving force J: The complete system is solved in OpenFOAM ® in a periodic channel, and compared with the the equivalent one-dimensional model.By applying the assumptions of uni-directional flow and zero derivatives in x we obtain a one-dimensional reduced set of equations in y: where we excluded the incompressibility condition as this is trivially satisfied.
To compare results between Chebfun and our single ionic fluid solver pnpFoam we the following dimensionless numbers, geometrical parameters and transport properties, according to table 2.
When running pnpFoam we solve it in steady state and take a large number of PIMPLE corrector iterations to insure convergence of the coupled system.The mesh used in OpenFOAM ® is up to N = 420 cells (100 in x, 320 in y).Comparison of results and the L 2 error norm convergence are shown in fig. 3. We observe an accumulation of c along Γ b , causing a large pressure gradient to form.sThe velocity profile follows a parabolic arc, as expected, whilst we see a non-linear profile for φ.Overall we find good agreement, with linear order spatial convergence O(N −1 ).

Multi-component ionic fluid
To verify pnpFoam and pnpMultiFoam for a multi-component fluid we take two species c 1 and c 2 with opposite valencies z 1 = −z 2 = 1 over the domain Ω ∈ [0, L] where L = 1 × 10 −6 .We set a zero-flux BC, i.e., j(c i ) = 0 at x = 0, L for both species and fix φ = 0.1 at x = L and φ = 0 at x = 0. We fix here a constant velocity u = 0.001 for simplicity and we by-pass the solution of the Stokes system (frozenFlow flag).Concentrations and potential are initially set to c i = 1 × 10 −3 and φ = 0.05(1 − cos πx L ).We consider here the dimensionless number values in table 3, indicating electrostatic forces dominate.From L we find the Debye length as approximately λ D = 308 nm.Since we are fixing u and neglecting Stokes we have P = 0: Table 3.: Dimensionless numbers and parameters for the multi-component ionic fluid testcase.
When running pnpFoam we use a time step ∆t = 10 −7 s, end time t = 1 × 10 −5 s, third-order implicit time scheme (backward keyword) and a mesh of up to N = 1000 cells.Results and L 2 norm convergence plots are depicted in Fig. 4 where we find good agreement of results and linear spatial convergence of all fields.
We see the ions are transported to the outer walls due to the high electric potential gradient.Most of the ions then accumulate within λ D from the walls forming two overlapping EDLs.Non-linear behaviour between c i and φ is observed through the slight shift in φ's profile due to the clustering of ions at the walls.

Reactive interface
Here we verify the accuracy of mappedChemicalKinetics, the numerical counterpart to the conditions found in §4.2.1.Consider the domain , each containing the species c f and c s respectively.Both species are restricted to their respective sub-domain, i.e. c s = 0 in Ω f and vice versa.Consider the following elementary reaction at the x = 0 interface: As c s and c f are restricted to Ω s and Ω f respectively the reactive conditions of (50) and (51) become where we have used the linear rate law to model the reaction rate.Note c f and c s in (83) and (84) are evaluated on Ω f and Ω s 's side of x = 0 respectively.We set both species as uncharged, i.e. z i = 0, to remove electrostatic effects and u = 1 in Ω f to ignore Stokes.We initially set c s = 0 in Ω s and c f = e −200(x+0.5) 2 in Ω f so that a clear increase in c s can be seen.
Looking at the dimensionless numbers for the problem in table 4 we can notice here that advection dominates.We also include two more numbers from the reactive interface conditions at x = 0.These are Da II,f and Da II,r , the Damköhler numbers given as the ratio between the diffusion rate and reaction rates, forward and reverse, respectively: Symbol N Pe L P Da II,f Da II,r Value 0 1 0 0 10 100

Table 4.: Table of dimensionless values
When running pnpMultiFoam to show qualitative comparison we use a time step ∆t = 10 −3 s, secondorder implicit time scheme backwards and mesh discretisation of N = 1000 cells (500 in Ω f , 500 in Ω s ).To show convergence of results and compute the L 2 error point norm we move to the steadyState time scheme to compare steady states.Results of the transient case are shown on the left in Fig. 6 where we find good agreement.First an initial diffusion and advection of c f is seen towards the interface, once c f reaches x = 0 it reacts to form c s in Ω s , where afterwards c s diffuses through Ω s .After t = 15s a steady state is reached implying chemical equilibrium of the reaction.As for the L 2 error norm of the steady case, seen in Fig. 6, we find a second-order convergence of mappedChemicalKinetics to Chebfun.This is because the interface conditions are linear and therefore exactly approximated by the linear approximation.

Random porous REV
To demonstrate the capabilities of our solvers, we consider a randomly generated porous solid-fluid domain Ω = Ω s ∪ Ω f , see Fig. 7.The random generation is done by a truncated Gaussian random field [12,24], sampled at the mesh points and categorised into two bins, denoting solid and fluid cells, by a threshold.Raising or lowering this threshold alters the porosity of the domain.
To have Ω be a microscopic representative elementary volume (REV) of a much larger macroscopic porous medium, we apply periodic conditions along the outer boundaries To generate movement, we fix a jump in φ between Γ in,s and Γ out,s , where the subscript s denotes the section of Γ in neighbouring Ω s .This is to mimic a fixed applied potential difference across the macroscopic medium, such as an applied voltage across a battery cell.We consider two ion species c 1 and c 2 with opposite valencies z 1 = −z 2 = −1.The height and length of the region is set as H = 1 × 10 −4 m.In Ω f we start with a uniform concentration of c 2 , the same is done for Ω s with c 1 .Fluid is initially taken at rest.Along Γ sf we set no flux for c 1 and flux continuity for c 2 .For φ we assume continuity of the electric displacement D = εE along Γ sf .In table 5 we list the dimensionless number values of the case: Symbol N Pe L P Value 0.4 9.65 × 10 −4 0.95 1 Table 5.: Dimensionless numbers for the random porous media testcase.
The system is initialised with constant c 1 and c 2 in Ω f and Ω s respectively.All fields are periodic on the outer boundaries Γ ext = Γ in ∪ Γ out ∪ Γ t ∪ Γ b , apart from φ where we introduce a quasi-periodic condition with a fixed jump between Γ in,s and Γ out,s , denoting the parts of Γ in and Γ out neighbouring Ω s .We apply no flux and flux continuity for c 1 and c s along Γ sf respectively.For φ we assume continuity of the electric displacement D = εE along Γ sf : To show consistency in the results we run four simulations in total.This is done for two randomly generated meshes, each discretised over two levels of refinement for N = 1 × 10 4 and N = 4 × 10 4 cells.We will use letters A and B to differentiate between the random generations and subscripts 1 and 2 for the levels of refinement.E.g, A 2 denotes the first random mesh, refined using 4 × 10 4 cells.All runs use a time step ∆t = 5 × 10 −4 and second-order implicit time scheme backwards.To observe the advective and electrostatic effects on our ions we define the following velocities in such that the continuity equation (5) in Ω f may be written as Figure 8 shows the results of the first random realisation with a finer mesh (N = 1 × 10 4 ), where we observe a gradient in φ formed between regions of Ω s connected at Γ in and Γ out .This is primarily due to our jump condition for φ applied along Γ in,s and Γ out,s .Whilst movement of our ion species have an effect on this gradient, for the most part it remains steady.The applied electrical driving force causes the uniform concentration c 1 in Ω f to accumulate around the regions of Ω s connected at Γ out , where a higher electric potential is present.Conversely, around Ω s connected at Γ in we find a much lower, but not zero, concentration of c 1 .A similar observation is found with c 2 , once diffused out of Ω s , accumulating where there is a lower potential.As with c 1 , at the region of Ω s connected at Γ out there is a much lower, but again not zero, concentration of c 2 .We observe therefore the formation of EDL's, see §3, around these connected solids.Given more time to evolve these EDL's would become further apparent, with more of c 2 diffusing out of Ω s .As for v 1 and v 2 we find their magnitudes ||v 1 || and ||v 2 || are highest within these EDL's located at the connected solids due to higher gradients in φ.
We can compare these results with coarse mesh (N = 1 × 10 4 ) results as well as with other random generation seed.Results of first realisation with a coarse mesh, and of a coarse and fine simulation of a second random realisation are reported respectively in fig. 9

Conclusions
In this paper we present the development of open-source solvers in OpenFOAM ® capable of simulating microscopic electrokinetic flows of dilute ionic fluids.The underlying system, known as Stokes-Poisson-Nernst-Planck, is thoroughly reviewed where we discuss the assumptions on the fluid properties and    limitations these bring.We later apply dimensional analysis to characterise the effects of dominating physical forces.This analysis is later used to give better understanding to a common model reduction, known as the electro-neutrality assumption, as a result of such analysis.Many real-world applications of these flows involve some form of reactions between involved ionic species at the interface between fluid and solid.As such, we formulate a fully general reaction model capable of describing said heterogeneous reactions as a balance of fluxes across all reacting species.
To verify flow descriptions obtained from our solvers are accurate, we compare against highly accurate solutions obtained, for simplified cases, with the Matlab toolbox Chebfun for a number of cases.With each case we find good agreement, showing spatial grid convergence of the results.We later use our solvers across a randomly generated solid-fluid porous cell, similar to a representative elementary volume (REV) used in homogenisation theory, to show the solvers ability in handling complex micro-structures.
In the future we plan on constructing full physically realistic cases, such as modelling the electrokinetic flow through a porous battery half-cell.Furthermore we later plan to apply uncertainty quantification in conjunction with these randomly generated porous cells, quantifying the effects on the flow for varying geometrical properties.
whereby setting c = c N +1 is the next iterative solution.We denote here ∂J ∂c to be the Jacobian matrix of J. To determine the Jacobian we define the directional derivative for some function G(c) along the direction δc = c − c N as: Using this definition of the directional derivative we can then determine the Jacobian within (97).Since the method is similar for all components of J, we will focus on the first element J 1 .Taking the directional derivative of J 1 with respect to c s along δc s at c = c N gives: Taking the following Taylor expansions about ε = 0 of r and the concentration gradient then substituting into (99) we obtain: The same can be done to find the directional derivative of J 1 with respect to c f , giving: Repeating the same process for J 2 of (96) and collecting all elements together and evaluating at c N gives ( where r and its derivatives are evaluated at c N .Re-arranging the first element of (105), using the first order approximation ∇ n c s = B s (c s − c s,c ) where c s,c is the value at the cell centre of cells sharing a face with Γ, and substituting into the second element then gives: (u

A.2. General non-linear unrestricted reaction rate
Here we will formulate a set of iterative Robin conditions required to solve the reactive conditions for unrestricted species and non-linear reaction.Consider a general non-linear reaction with i = J, ..., K reactants and i = K + 1, ..., M products: We model (114) using the following conditions with general, possibly non-linear, reaction rates r i : We begin by writing j i in full and use n = −n to give • n , a more acceptable form for OpenFOAM ® .Since formulating the necessary Robin conditions is the same for all species i, we omit the subscript from here onward unless necessary.Instead we introduce two new subscripts s and f denoting evaluation on the solid or fluid side of Γ: We next use the first order approximation ∇ n c f = B f (c f − c f,c ), where the coefficient B f is determined by OpenFOAM ® and c f,c denotes the cell centre value of cells sharing a face with Γ. Alongside the continuity condition (116) this allows us to rearrange (117) into a form close to Robin: Since r i may be non-linear, we linearise about an initial guess c N = (c N J , c N J+1 , ..., c N i , ..., c N M ) of concentrations and define δc j = c j − c N j : r i (c j=J,...,M ) = r i (c N ) + δc J ∂r i ∂c J c N + δc J+1 ∂r i ∂c J+1 c N + ... + δc i ∂r i ∂c i c N + ... + δc M ∂r i ∂c M c N + O(δc 2 ) (119) Substitute (119) into (118) and take c s = c N +1 s to be the next iterative solution.To decouple the iterative conditions, all other concentrations we set as c j = c N j for i = j.The result is the following iterative Robin condition for c N +1 s , the concentration of species i on the solid side of Γ: The same process can be done to determine the iterative Robin condition for c f = c N +1 f , the concentration of species i on the fluid side of Γ: ∂c i c N (127)

Figure 1 .
Figure 1.: A graphical representation of the electric double layer along a positively charged surface, comprised of a thin layer of adsorbed negative ions (Stern layer) and loosely connected positive and negative ions (Diffuse layer).

Algorithm 3 :
mappedChemicalKinetics Set tolerance: tol; Set max Newton iterations: N max ; Reset iterator: N = 0; while ||c N f − c N +1 f || > tol ∩ N ≤ N max do Store previous iteration: c N f = c N +1 f ; Extract the face and internal values of the neighbouring field (solid side); Update D * f , K * f and F * f : see (108)-(110) ; Evaluate iterative effective Robin BC: −D * f ∇ n c N +1 f = −K * f c N +1 f − F * f ; Apply the boundary condition to obtain the new value: c N +1 f ; if Stiff then

Figure 2 . 2 0
Figure 2.: Graphical representation of the infinite ion channel with channel width H.

Figure 4 .
Figure 4.: Left: Results of c 1 , left curve, and c s , right curve.Blue dots denote Chebfun and red line for pnpFoam at t = 1 × 10 −5 s with initial conditions given in black., Middle: Results of φ for Chebfun (blue dots) and pnpFoam (red line) at t = 10 −5 s, with initial condition given in black., Right: L 2 point error norm convergence plot.

Figure 5 .
Figure 5.: Graphical representation of domain Ω containing a solid and fluid with zero flux outer boundaries and reactive flux conditions at the shared interface.

Figure 6 .
Figure 6.: Left: Molar concentration results of c f in fluid (dashed for pnpMultiFoam, markers for Chebfun) and c s in solid (solid for pnpMultiFoam, markers for Chebfun).Blue for initial concentrations, black at intermediate times, red for final time (i.e steady state), Right: L 2 norm convergence plot between Chebfun and mappedChemicalKinetics/pnpMultiFoam

Figure 7 .
Figure 7.: Visual representation of the random porous domain containing solid Ω s (grey) and fluid Ω f (white) with shared interface Γ sf .
Figure8shows the results of the first random realisation with a finer mesh (N = 1 × 10 4 ), where we observe a gradient in φ formed between regions of Ω s connected at Γ in and Γ out .This is primarily due to our jump condition for φ applied along Γ in,s and Γ out,s .Whilst movement of our ion species have an effect on this gradient, for the most part it remains steady.The applied electrical driving force causes the uniform concentration c 1 in Ω f to accumulate around the regions of Ω s connected at Γ out , where a higher electric potential is present.Conversely, around Ω s connected at Γ in we find a much lower, but not zero, concentration of c 1 .A similar observation is found with c 2 , once diffused out of Ω s , accumulating where there is a lower potential.As with c 1 , at the region of Ω s connected at Γ out there is a much lower, but again not zero, concentration of c 2 .We observe therefore the formation of EDL's, see §3, around these connected solids.Given more time to evolve these EDL's would become further apparent, with more of c 2 diffusing out of Ω s .As for v 1 and v 2 we find their magnitudes ||v 1 || and ||v 2 || are highest within these EDL's located at the connected solids due to higher gradients in φ.We can compare these results with coarse mesh (N = 1 × 10 4 ) results as well as with other random generation seed.Results of first realisation with a coarse mesh, and of a coarse and fine simulation of a second random realisation are reported respectively in fig.9, fig. 10, and fig.11.
s ∇ n δc s − δc s ∂r ∂cs − δc f ∂r ∂c f δc s ∂r ∂cs + (u f • n )δc f − z f D f F RT δc f ∇ n φ f + δc f r are evaluated at c N .We may then substitute (104) into (97) to give  −D s ∇ n c s − r − δc f ∂r ∂c f − δc s ∂r ∂cs

P=
(c f ) = ∂r ∂c s   −D s B s c s,c + r + δc f ∂r ∂c f − c N s ∂r ∂cs D s B s + ∂r ∂cs   (107) Setting c = c N +1 to indicate the next iterative solution results in an iterative Robin condition of the form −D* f ∇ n c N +1 f = −K * f c N +1 f − F * f , with effective coefficients dependent on the previous iteration c N : f F RT ∇ n φ f + (u f • n ) s B s c s,c + r − cand its derivatives are evaluated at the previous iteration c N .Similar steps can also be applied to the accompanying iterative Robin condition for c s of the form−D * s ∇ n c N +1 s • n ) − D f B f − z f D f F RT ∇ n φ f + ∂r ∂c f f B f c f,c − r + c N f ∂r ∂c f + c N s ∂r ∂cs

F
f = D i B s c s,c − r i (c N ) + c N f ∂r i