Computational framework for resolving boundary layers in electrochemical systems using weak imposition of Dirichlet boundary conditions

https://doi.org/10.1016/j.finel.2022.103749Get rights and content

Highlights

  • Finite element approach to solve equations that model charged species transport.

  • Dirichlet-to-Neumann transformation for weak implementation of species concentration and potential boundary condition.

  • Accurate boundary flux calculation even with poor mesh resolution.

  • Dirichlet-to-Neumann transformation provides reliable simulation results for thin boundary layer.

  • Illustrate this approach on canonical 3D problems in electrochemical separations.

Abstract

We present a finite element based computational framework to model electrochemical systems. The electrochemical system is represented by the coupled Poisson–Nernst–Planck (PNP) and Navier–Stokes (NS) equations. The key quantity of interest in such simulations is the current (flux) at the system boundaries. Accurately computing the current flux is challenging due to the small critical dimension of the boundary layers (small Debye layer) that require fine mesh resolution at the boundaries. We present a numerical framework which resolves this challenge by utilizing a weak imposition of Dirichlet boundary conditions for Poisson–Nernst–Planck equations. In this numerical framework we utilize a block iterative strategy to solve NS and PNP equations. This allows us to efficiently and easily implement the weak imposition of Dirichlet boundary conditions. The results from our numerical framework shows excellent agreement when compared to strong imposition of boundary conditions (strong imposition requires a much finer mesh). Furthermore, we show that the weak imposition of the boundary conditions allows us to resolve the fluxes in the boundary layers with much coarser meshes compared to strong imposition. We also show that the method converges as we refine the mesh near the boundaries at a much faster rate compared to strong imposition of the boundary layer. We present multiple test cases with varying boundary layer thickness to illustrate the utility of the numerical framework. We illustrate the approach on canonical 3D problems that otherwise would have been computationally intractable to solve accurately. Lastly, we simulate electrokinetic instabilities near a perm-selective membrane with weakly imposed boundary conditions on the membrane. This approach substantially reduces the computational cost of modeling thin boundary layers in electrochemical systems.

Introduction

An understanding of charged species transport is critical to the development of electrochemical and electrokinetic systems relevant to a wide range of disciplines (engineering, chemistry, physics) and applications (sensing, energy, water purification). Systems that employ non-linear electrokinetics, in which the electric field is varied spatially (and also temporally in some cases), are especially difficult to model. For example, charged species can be electrokinetically focused along a steep electric field gradient formed near an ion-selective membrane or a bipolar electrode (BPE) [1], [2]. In both cases, the electric field gradient results from the local depletion of charge carriers at one end of the membrane (by selective charge transport) or BPE (by Faradaic reactions) [3]. The formation of an ion depletion zone (IDZ) and ion enrichment zone (IEZ) at opposite sides of the membrane or BPE is called ion concentration polarization (ICP). A few prominent applications include water purification and desalting [4], [5], biomedical engineering [5], and enrichment and detection of trace analytes [6], [7]. In all of these applications, the stability of the IDZ drastically limits the volumetric throughput of these devices. Therefore, the ability to simulate species transport in these systems is critical to their advancement.

Species transport in electrochemical systems, such as ICP, is a complex multi-physics problem driven by diffusion, electromigration, and convection [8]. Experimental approaches employed to characterize this multi-physics problem are generally limited to the measurement of electrical current or to the visualization of fluorescent tracer molecules. As a result, it is difficult to fully understand the mechanism of ICP using these methods alone. Therefore, there have been several numerical studies of transport in such systems to compliment experimental results. For example, Zangle et al. [9], [10] derived 1-dimensional (1D) governing equations for ICP in a system comprising a micro–nano–micro junction and calculated shock wave-like IDZ and IEZ propagation along the microchannel segments, originating at the nanochannel. Using this approach, they found that Dukhin number (surface conductivity over bulk fluid conductivity) and the electrophoretic mobility of charged species are the parameters that dictate the rate and extent of the propagation. Numerical results obtained by Mani and collaborators showed that chaotic fluid motion originates from the locally high electric field [11] or alternating current (AC) [12] even in the low Reynolds number regime. ICP is made further complex when, in addition to convection, diffusion, and migration, chemical reactions are involved. To address such a case, Kler et al. included reaction terms to simulate electrophoresis accompanied by acid and base reactions [13]. Similarly, Tallarek and coworkers included acid/base and Faradaic reaction terms in the simulation of ICP at BPEs [14].

Although these studies exemplify successful simulation of charged species transport in non-linear electrokinetics, it is still challenging to obtain reliable results with a reasonable computational cost. The primary reason for this difficulty is the multi-scale nature of the problem [15]. The smallest scale feature that impacts the physics in electrochemical and electrokinetic systems is the electrical double layer ( 10nm) or EDL, which comprises electrical potential and ion concentration gradients in the boundary layer present at a liquid–solid interface. In contrast, species transport relevant to most applications of such systems extend over length scales of 10μm to 1000μm. Importantly, unresolved boundary layers can result in unfavorable oscillations extending outside of the boundary layers into the bulk domain. Refining the mesh near the boundary is a reasonable approach to address challenges from multi-scale characteristics [11], [16] and provides reliable results in the entire domain. However, the computational cost incurred by the increased mesh density in the boundary layer can be prohibitive for very small Debye lengths. Considering that most applications are interested in what happens in the ‘bulk’ of the fluid domain (and its impact on current flux), not in the vicinity of the boundary, resolving the mesh near this layer is not computationally economical. Jia et al. [17], [18], using a commercial code, simplified complex boundary physics with electroosmotic slip velocity, which minimizes computational costs. However, there is ambiguity in the selection of the location where the slip boundary condition imposed away from an ion selective membrane transitions to a no-slip boundary condition, imposed on or adjacent to the membrane. Moreover, replacing the boundary layer with the slip boundary condition ignores concentration gradient driven flow [19], [20]. Therefore, there remains a need to reduce the computational cost of representing the boundary layer without oversimplifying the underlying physics.

In this work, we address this need by utilizing an approach used in fluid mechanics – the weak imposition of Dirichlet condition – this is used to efficiently model the no-slip condition (Dirichlet boundary condition) [21]. The weak imposition of Dirichlet condition, also known as Nitsche’s method or symmetric interior penalty Galerkin method (SIPG) if applied to element boundaries, provides a consistent and robust way of enforcing Dirichlet conditions by variational weakening of the no-slip condition into a Neumann type condition, especially in the context of Finite Element (FE) analysis. [22], [23] Such a strategy releases the point-wise no-slip condition imposed at the boundary of the fluid domain, thus minimizing the mesh resolution required to track the steep gradients close to the boundaries. This effect reliably imitates the presence (and effect) of the thin boundary layer [24], [25]. Enforcing Dirichlet boundary conditions weakly allows for an accurate overall flow solution even if the mesh size in the wall-normal direction is relatively large. This approach has substantially benefited efficient simulations of turbulent flow scenarios [26], [27] as well as other multi-physics flow scenarios [28], [29].

The present study develops a FEM framework for the fully coupled Navier–Stokes and Poisson–Nernst–Planck (NS–PNP) equations to simulate electrochemical and electrokinetic systems. To overcome difficulties from thin boundary layers, Dirichlet boundary conditions are weakly enforced [21], [30] in the PNP equations. While usage of the developed framework is not limited to specific application to electrochemical and electrokinetic systems, we demonstrate its efficacy in resolving calculations of electroosmotic flow (EOF) and ion concentration polarization (ICP). EOF and ICP were selected as test cases for two reasons — first, there has been growing interest in these phenomena due to their potential impact in chemical, biomedical, and environmental fields, and second, because these examples include the three fundamental transport mechanisms — convection, diffusion, and migration. Our findings are significant because, despite a much coarser mesh, the results obtained with weak BC showed good agreement to those obtained with strong BC, and furthermore, boundary flux calculations converged much faster to the solution using the weak BC. Collectively, these results demonstrate a significant reduction in computational load while retaining accuracy. Therefore, we expect that this weak BC approach will provide greater stability and accuracy in the simulation of a wide range of electrochemical and electrokinetic systems.

The outline of the rest of the paper is as follows: We begin by revisiting the governing equations for charged species transport followed by the non-dimensional forms of these governing equations in Section 2. Then, the FEM problems are defined with weakly imposed Dirichlet boundary conditions in Section 3. The solving strategy for the numerical method was discussed in Section 4. In Section 5, the developed framework is validated with a manufactured solution and by simulation of EOF, for which an analytical solution exists. Next, 1D and 2D IDZs are simulated with weakly imposed Dirichlet boundary conditions. For 3D applications, the developed platform was tested for electrolyte separation (desalting) in a microchannel. For each example, the calculated boundary flux is compared with that obtained with strongly imposed Dirichlet boundary conditions. We also test the weakly imposed boundary conditions for the simulation of the electrokinetic instabilities near a perm-selective membrane. We conclude in Section 6.

Section snippets

Governing equations

Poisson–Nernst–Planck (PNP): Without loss of generality, we consider a canonical problem of solvent flow and species transport in a (micro)channel configuration. We consider an incompressible fluid (Eq. (10)). This problem encompasses both the electroosmotic and pressure driven regimes. We consider N>1 number of charged species with subscript i indicating the species index. The species flux ji,3

Weak form of the equations

Consider the spatial domain as Ω, with Ω as its boundary, and by Γ the boundary where the weak boundary conditions are enforced. We can define the variational problem as follows.

Definition 1

Let (,) be the standard L2 inner product over the subscript (i.e. either Ω or ΩΓ). We state the variational problem as follows: find u(t,x)L20,T;H01Ω with 7

Strategy for implementation

As specified before, we use a block iterative strategy to solve the set of equations. This approach provides several advantages, including (a) reducing the number of degrees of freedom per solve, (b) mitigating the numerical stiffness that exists between the equations (especially the large body force in the momentum equation), (c) enabling (simplified) weak imposition of Dirichlet boundary conditions by allowing separate treatment for ci and ϕ (see the boundary condition terms in Eqs. (32), (33)

Convergence against manufactured solution

We use the method of manufactured solutions to assess the convergence of our implementation. We select an input “solution”, and substitute it in the full set of governing equations. We then use the residual as a body force on the right-hand side of Eqs. (28), (29), (30). We choose the following “solution” with appropriate body forcing terms: u=cos(2πt)sin(2πx)cos(2πy),v=cos(2πt)cos(2πx)sin(2πy),p=cos(2πt)sin(2πx)cos(2πy),c+=cos(2πt)cos(2πx)sin(2πy),c=cos(2πt)sin(2πx)cos(2πy),ϕ=cos(2πt)cos(2πx

Conclusion

In this study, we simulate electrokinetic systems represented by the Navier–Stokes–Poisson–Nernst–Planck equations, with a key quantity of interest being the current flux at the system boundaries. Accurately computing the current flux is challenging due to the thin boundary layers (small Debye lengths) that require fine mesh to resolve. We address this challenge by using the Weak imposition of Dirichlet conditions. The framework was validated against manufactured solutions and the analytical

CRediT authorship contribution statement

Sungu Kim: Conceptualization, Methodology, Software, Formal analysis, Writing, Visualization. Makrand A. Khanwale: Methodology, Formal analysis, Writing, Visualization. Robbyn K. Anand: Conceptualization, Resources, Supervision, Writing. Baskar Ganapathysubramanian: Conceptualization, Methodology, Analysis, Writing, Resources, Supervision, Project administration.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

We thank Kumar Saurabh from the Ganapathysubramanian group for technical discussions and implementation support, as well as for proof-reading the manuscript.

Funding

The authors acknowledge XSEDE grant number TG-CTS110007 for computing time on TACC Stampede2. BG was funded in part by NSF grant 1935255, 1855902. RKA and SK were funded in part by an NSF CAREER grant awarded by the Chemistry Directorate Chemical Measurement and Imaging Program under award number 1849109.

References (50)

  • BerzinaB. et al.

    Tutorial review: Enrichment and separation of neutral and charged species by ion concentration polarization focusing

    Anal. Chim. Acta

    (2020)
  • KwakR. et al.

    Microscale electrodialysis: Concentration profiling and vortex visualization

    Desalination

    (2013)
  • LiM. et al.

    Recent advancements in ion concentration polarization

    Analyst

    (2016)
  • MavreF. et al.

    Bipolar electrodes: a useful tool for concentration, separation, and detection of analytes in microelectrochemical systems

    (2010)
  • BondarenkoM.P. et al.

    Current-induced ion concentration polarization at a perfect ion-exchange patch in an infinite insulating wall

    ChemElectroChem

    (2020)
  • BerzinaB. et al.

    An electrokinetic separation route to source dialysate from excess fluid in blood

    Anal. Chem.

    (2018)
  • AnandR.K. et al.

    Bipolar electrode focusing: faradaic ion concentration polarization

    Anal. Chem.

    (2011)
  • KimS. et al.

    Concentration enrichment, separation, and cation exchange in nanoliter-scale water-in-oil droplets

    J. Am. Chem. Soc.

    (2020)
  • ProbsteinR.F.

    Physicochemical Hydrodynamics: An Introduction

    (2005)
  • ManiA. et al.

    On the propagation of concentration polarization from microchannel- nanochannel interfaces Part I: analytical model and characteristic analysis

    Langmuir

    (2009)
  • ZangleT.A. et al.

    On the propagation of concentration polarization from microchannel- nanochannel interfaces Part II: numerical and experimental study

    Langmuir

    (2009)
  • DruzgalskiC. et al.

    Direct numerical simulation of electroconvective instability and hydrodynamic chaos near an ion-selective surface

    Phys. Fluids

    (2013)
  • KlerP.A. et al.

    Modeling and high performance simulation of electrophoretic techniques in microfluidic chips

    Microfluid. Nanofluid.

    (2011)
  • HlushkouD. et al.

    Numerical simulation of electrochemical desalination

    J. Phys.: Condens. Matter

    (2016)
  • BoyD.A. et al.

    Simulation tools for lab on a chip research: advantages, challenges, and thoughts for the future

    Lab Chip

    (2008)
  • Cited by (10)

    View all citing articles on Scopus
    1

    Presently at Department of Mechanical Engineering, Stanford University, CA, USA.

    2

    Presently at Center for Turbulence Research, Department of Mechanical Engineering, Stanford University, CA, USA.

    View full text