Computational framework for resolving boundary layers in electrochemical systems using weak imposition of Dirichlet boundary conditions
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 ( ) 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 m to 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 number of charged species with subscript indicating the species index. The species flux ,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 inner product over the subscript (i.e. either or ). We state the variational problem as follows: find 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 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:
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)
- et al.
Direct seawater desalination by ion concentration polarization
Nature Nanotechnol.
(2010) - et al.
Characterization of chaotic electroconvection near flat inert electrodes under oscillatory voltages
Micromachines
(2019) - et al.
Weak imposition of Dirichlet boundary conditions in fluid mechanics
Comput. & Fluids
(2007) - et al.
Weak Dirichlet boundary conditions for wall-bounded turbulent flows
Comput. Methods Appl. Mech. Engrg.
(2007) - et al.
Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows
Comput. Methods Appl. Mech. Engrg.
(2007) - et al.
A residual-based variational multiscale method with weak imposition of boundary conditions for buoyancy-driven flows
Comput. Methods Appl. Mech. Engrg.
(2019) - et al.
Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations
Comput. Methods Appl. Mech. Engrg.
(1982) - et al.
Incompressible flow computations with stabilized bilinear and linear equal-order-interpolation velocity-pressure elements
Comput. Methods Appl. Mech. Engrg.
(1992) - et al.
A better consistency for low-order stabilized finite element methods
Comput. Methods Appl. Mech. Engrg.
(1999) - et al.
Galerkin/least-squares finite element methods for the reduced wave equation with non-reflecting boundary conditions in unbounded domains
Comput. Methods Appl. Mech. Engrg.
(1992)
Tutorial review: Enrichment and separation of neutral and charged species by ion concentration polarization focusing
Anal. Chim. Acta
Microscale electrodialysis: Concentration profiling and vortex visualization
Desalination
Recent advancements in ion concentration polarization
Analyst
Bipolar electrodes: a useful tool for concentration, separation, and detection of analytes in microelectrochemical systems
Current-induced ion concentration polarization at a perfect ion-exchange patch in an infinite insulating wall
ChemElectroChem
An electrokinetic separation route to source dialysate from excess fluid in blood
Anal. Chem.
Bipolar electrode focusing: faradaic ion concentration polarization
Anal. Chem.
Concentration enrichment, separation, and cation exchange in nanoliter-scale water-in-oil droplets
J. Am. Chem. Soc.
Physicochemical Hydrodynamics: An Introduction
On the propagation of concentration polarization from microchannel- nanochannel interfaces Part I: analytical model and characteristic analysis
Langmuir
On the propagation of concentration polarization from microchannel- nanochannel interfaces Part II: numerical and experimental study
Langmuir
Direct numerical simulation of electroconvective instability and hydrodynamic chaos near an ion-selective surface
Phys. Fluids
Modeling and high performance simulation of electrophoretic techniques in microfluidic chips
Microfluid. Nanofluid.
Numerical simulation of electrochemical desalination
J. Phys.: Condens. Matter
Simulation tools for lab on a chip research: advantages, challenges, and thoughts for the future
Lab Chip
Cited by (10)
Direct numerical simulation of electrokinetic transport phenomena in fluids: Variational multi-scale stabilization and octree-based mesh refinement
2024, Journal of Computational PhysicsNew mixed finite element methods for the coupled Stokes and Poisson Nernst Planck equations in Banach spaces
2023, ESAIM: Mathematical Modelling and Numerical Analysis
- 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.