A multi-scale model for electrokinetic transport in networks of pores

networks of pores, Part II: computational algorithms and applications ... ‡Center for Turbulence Research, Stanford University, Stanford, California...
0 downloads 0 Views 3MB Size
Subscriber access provided by CORNELL UNIVERSITY LIBRARY

Article

A multi-scale model for electrokinetic transport in networks of pores, Part II: computational algorithms and applications Shima Alizadeh, and Ali Mani Langmuir, Just Accepted Manuscript • Publication Date (Web): 16 May 2017 Downloaded from http://pubs.acs.org on May 17, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Langmuir is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

A multi-scale model for electrokinetic transport in networks of pores, Part II: computational algorithms and applications Shima Alizadeh† and Ali Mani∗,†,‡ Department of Mechanical Engineering, Flow Physics and Computational Engineering, Stanford University, Stanford, California 94305, USA E-mail: [email protected]

Phone: +1650 7251325. Fax: +1650 7253525

Abstract The first part of this two-paper series presented a robust mathematical model for fast and accurate prediction of electrokinetic phenomena in porous networks with complex topologies. In the second part of this series we first present a numerical algorithm that can efficiently solve the model equations. We then demonstrate that the resulting framework is capable of capturing a wide range of transport phenomena in microstructures by considering a hierarchy of canonical problems with increasing complexity. The developed framework is validated against direct numerical simulations of deionization shocks in micro-pore-membrane junctions and concentration polarization in micro- and nano-channel systems. We demonstrate that for thin pores subject to concentration gradients, our model consistently captures correct induced osmotic pressure, which ∗

To whom correspondence should be addressed Department of Mechanical Engineering, Flow Physics and Computational Engineering, Stanford University, Stanford, California 94305, USA ‡ Center for Turbulence Research, Stanford University, Stanford, California 94305, USA †

1

ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

is a macroscopic phenomena originally derived from thermodynamic principles, but here naturally predicted through microscopic electrostatic interactions. Moreover, we show that the developed model captures current rectification phenomenon in a conical nano-pore subject to an axial external electric field. Lastly, we provide discussions on examples involving stationary and moving deionization shocks in micro-pore nano-pore T-junctions, as well as induced flow loops when pores of varying sizes are connected in parallel.

1

Introduction

Many industrial applications involving porous materials are driven by the coupling between fluid flow, ion transport, and electrostatics. Microfluidic lab-on-a-chip systems are fabricated microstructures used for a variety of processes including biological separation, mixing, and bio-molecule sensing. 1–5 Electrochemical cells often consist of porous elements to conduct processes including energy storage 6–9 and conversion, 10,11 pumping, 12 desalination, 24 and decontamination. Other relevant examples include hydraulic fracturing for shale gas extraction, 13,14 or electrokinetic decontamination of soil and porous rocks. 15–17 The porous structures in many of these applications involve complex topologies with large surface to volume ratios, 18–24 which can cause a variety of nonlinear dynamics with impacts on the system performances. Deionization shock 25,27,28 is one example of such dynamics, which allows the passage of overlimiting current 24,31,32 and provides a novel mean for water desalination 24,25,30 and soil decontamination. Additional complexities that can arise in confined structures are internally induced flow circulations, and current rectification by nano-scale pores. 33–39 Having a precise perception of these complex dynamics can shed light on the optimal engineering of synthetic porous media with applications in microfluidics and nanofluidics, desalination, energy storage and conversion. In Part I of this two-paper series, we derived a multi-scale model for the simulation of electrokinetic phenomena in networks of pores. Assuming that the pores of the network 2

ACS Paragon Plus Environment

Page 2 of 37

Page 3 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

are thin and long, our model solves a set of 1D transport equations in pores’ longitudinal directions. We account for the effects of surface charge and electric double layers (EDLs) by incorporating area-averaged coefficients in the flux terms associated with fluid flow, electric current, and flux of ionic species. The developed model is discretely conservative, does not involve any singularity in the limit of infinitely thick EDLs, preserves equilibrium condition in the absence of external forces, and allows simulation of general networks of pores with random connectivities. The objective of Part II of this two-paper series is to present a validation of the developed model and to demonstrate its capability in capturing a wide range of transport phenomena in porous structures. We begin by summarizing the model equations governing the fluid, ion, and charge transport in a general porous structure. We then describe our computational algorithms used to generate tables of area-averaged transport coefficients, and describe the methodology for temporal and spatial discretization of the transport equations. In Section 4, we present the results of our numerical simulations by considering a wide range of canonical settings involving electrokinetic transport in porous structures. We first present an investigation on the dependency of the induced osmotic pressure gradient in a pore with finite surface charge subject to an axial concentration gradient. We demonstrate that in the limit of small pore radius, our model captures the ideal pressure gradient, consistent with analytical predictions of osmotic pressure derived from thermodynamic principles. However, our model is capable of quantifying transition to the non-ideal limit as the pore size increases, since it utilizes direct forces at the micro-scale. We demonstrate additional validation test cases by comparing our model predictions with those obtained from direct numerical simulations of a micropore-membrane junction, and a microchannel-nanochannel-microchannel system exposed to an external pressure gradient. Furthermore, we consider current rectification by a conical nano-pore, where we illustrate the importance of advective flux on the rectification properties of the nano-porous systems. Finally, we discuss the temporal and spatial evolution of electrokinetic phenomena emerging in networks of micro-pores and nano-pores in parallel

3

ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 37

loops and H-junction configurations.

2

Mathematical model

In Part I of this two-paper series we derived the following area-integrated transport equation for coion species along pores’ longitudinal directions:

S(x)

∂C − ∂ + {S(x)F − } = 0, ∂t ∂x

(1)

where x coordinate is in pore’s axial direction and S(x) is the local pore area cross-section. C − and F − respectively represent the concentration and flux of coions, and the over-bar symbol represents cross-sectional averaging. In Part I we presented a comprehensive derivation of F − in terms of local C − and pore properties. Specifically, we identified two dimensionless local parameters, σ ∗ and λ∗ that determine the pre-factors in front of the flux terms arising from cross-sectional averaging procedure. σ ∗ represents the dimensionless local surface charge, and λ∗ represents the ratio of local Debye length, based on C − , to the local effective pore size, hp . Additionally, we leveraged the cross-sectional invariance of three driving potentials, 40,41 describing the fluxes induced by pressure-driven flow, electro-osmosis, and diffusio-osmosis. These potentials are virtual total pressure, P0 , virtual counterion electrochemical potential, µ+ , and virtual ion concentration C0 . The advantage of these choices is that their magnitude remains bounded when concentration becomes very small. Before describing the final model for F − , we first present a more compact notation indicating effects of advection, diffusion, and electromigration respectively: +

dC0 dµ F¯ − = uC − − 2¯ g + C− , dx dx

(2)

where g¯ is one of the aforementioned pre-factors that solely depends on known local σ ∗ g . uC − is the cross-sectional averaged and λ∗ . Given g¯, C0 can be obtained as C0 = C − /¯ 4

ACS Paragon Plus Environment

Page 5 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

advective flux which itself depends on the combination of pressure-driven, electroosmotic, and diffusio-osmotic flows. As described in Part I, we can re-arrange and write F − with respect to the local gradients of the primary driving potentials multiplied with pre-factors as follows:

F− =

κ dP0 dµ+ − g p + C g p− } − g e + C g e− ) + C − } {C + {κ(C s s 2λ2D dx dx +{

where λD (x) = and κ =

ε V2 µD T

1 hp

q

εkB T ˜ref z2 e2 2C

κ dC0 − g c + C g c− ) − 2¯ (C , g } s 2λ2D dx

(3)

is the ratio of reference Debye length to local effective pore size,

is electrohydrodynamic coupling constant, which is solely dependent on the

choice of electrolyte. Before describing the other coefficients introduced in (3), we discuss the solutions to the fields P0 and µ+ . To fully close equation (3), one can compute P0 and µ+ fields by imposing the conservation of fluid mass,

d {S(x)¯ u} dx

= 0, and the conservation of electric charge,

d {S(x)¯i} dx

= 0 inside

each pore in a porous structure. These equations respectively correspond to equations (59) and (60) in Part I of this two-paper series. u¯ and ¯i are the area-averaged fluid velocity and electric current density respectively, which can be written in the same format of equation (3): dµ+ dC0 κ κ p dP0 e + κg + 2 gc , u¯ = 2 g 2λD dx dx 2λD dx

(4)

+ ¯i = Cs κ {g p + g p+ − g p− } dP0 + {−(2C − + Cs ) + Cs κ(ge + g e+ − g e− )} dµ 2λ2D dx dx

+{Cs

κ dC0 c + g c+ − g c− ) + 2¯ (g . g } 2λ2D dx

(5)

In equations (3) to (5), g¯, g e , g c , g p− , g e− , g c− , g p+ , g e+ , and g c+ are the area-averaged coefficients that generally vary with space and time, but depend solely on local σ ∗ and λ∗ . These coefficients are obtained from interpolations of the data that are tabulated in terms of σ ∗ and λ∗ . Cs = 2|σ ∗ |λ2D represents the excess couterions needed to neutralized surface 5

ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 37

charge at each local cross-section. Therefore, the counterion concentration can be obtained as C + = C − + Cs . One may also represent equations (3) to (5) compactly by combining all pre-factors in front of each gradient term. However, the combined pre-factors depend on λD and κ, as well as σ ∗ and λ∗ , and thus their tabulation requires a high-dimensional table.

3

Computational Algorithms

For the test cases presented in this study we considered two types of pore cross-sectional geometries: cylindrical pores and slits. The latter represents wide and thin rectangular cross-sections where the effects of sidewalls can be ignored. For these pore geometries, the Poisson-Boltzmann equation reduces to a 1D ODE. To solve this equation in the transverse directions, we employed a finite volume approach utilizing uniformly spaced mesh points either in the wall normal direction of a slit or in the radial direction for a circular pore. Using Poisson-Boltzmann equation, we first solve for electrostatic potential (ψ) and

C± , C0

which are

set at computational cell centers. The Laplacian operator in the cross-sectional direction is discretized with second order accuracy. The details of the numerical algorithm leading to solutions to the Poisson-Boltzmann system, and computation of tabulated coefficients, are given in the supporting information. The area-averaged coefficients are computed over a wide range of σ ∗ and λ∗ values, and are tabulated with respect to these parameters. Given the smooth dependence on σ ∗ and λ∗ , these input variables are spaced logarithmically, such that every factor of 10 in σ ∗ and λ∗ is resolved by 10 points. To time-advance equation (1) and compute temporal evolution of ion concentration field, we used uniformly spaced mesh points along each pore longitudinal axis. Note that since the pore area cross-section can vary along the axis, the cell volumes may not be equal. All area-averaged coefficients and the driving potential fields, P0 , µ+ , C0 , and C − are defined at cell centers, whereas we compute the area-averaged fluxes at cell faces. we use second

6

ACS Paragon Plus Environment

Page 7 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

order differentiation and linear interpolation in order to approximate the derivatives and the area-averaged quantities at cell faces. For each internal and terminal node, we consider one computational element whose volume, σ ∗ , and λD are set based on the largest adjacent cell. The inlet and outlet surface areas of each interface with intersecting pores are equal to the cross-section areas of the pores at their connecting points with the interface. This choice provides reasonable numerical robustness when solving for C − in systems involving pore intersections. For the time integration of equation (1), we implemented a second order iterative implicit scheme to enhance the numerical stability of the developed model. The implicit time discretization of equation (1) is: 3C −

(n+1)

(n)

− 4C − 2∆t

+ C−

(n−1)

=

−1 ∂ (n+1) {SFx− } + O(∆t2 ), S ∂x

(6)

where superscript (n + 1) refers to the quantities in the next time step. (n) and (n − 1) correspond to the quantities at the current time step and the previous time step respectively. Fx−

(n+1)

is the area-averaged flux of negative ions at time (n+1), which includes the advection,

diffusion, and electromigration effects as described in equation (1). The existence of stiff terms such as the diffusion term necessitates the use of implicit time integration. With an implicit scheme, one can choose an appropriate time step that resolves the dynamics of unsteady physics, while bypassing the short time scale of quasi-steady processes associated with the stiff terms. To this end, we devised a fast implicit procedure, which is explained in detail in the supporting information. For each problem, the initial condition can be constructed by assuming that fields are globally in equilibrium prior to the start of applied external driving forces. Therefore, we often choose initial conditions for C − such that they correspond to uniform C0 . To demonstrate the capabilities of the developed model in predicting a wide range of electrokinetic phenomena, we present series of canonical and engineering problems in the next section.

7

ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

4 4.1

Page 8 of 37

Numerical results Osmosis phenomenon and its connection to pore size and surface charge

Using our tabulated coefficients, we aim to predict the induced pressure difference across a circular pore and identify the conditions over which the osmotic pressure is recovered across the pore. 41 Figure (1a) shows the schematic of a single pore whose ends lead to two closed large reservoirs with known and unequal salt concentrations. The left reservoir contains higher salt concentration compared to the right reservoir. Given the dead-end flow condition, one expects the imposed diffusio-osmosis flow to be blocked by an induced pressure gradient. However, the system is also coupled to an induced voltage gradient via charge advection by both pressure-driven and diffusio-osmotic flows. Another required condition is that the net current through the system must be zero. Therefore, to predict the induced pressure difference, one needs to set equations (4) and (5) to zero (given zero flow and current at steady state) and solve them in a coupled fashion. As shown in Part I of this two-paper series, equations (4) and (5) can be written in the following compact format for a differential element inside the pore, as shown in Figure (1a):

g11 (x)

dP0 dµ+ dC0 + g12 (x) + g13 (x) = 0, dx dx dx

(7)

g21 (x)

dP0 dµ+ dC0 + g22 (x) + g23 (x) = 0, dx dx dx

(8)

where g11 , g12 , and g13 , are the pre-factors associated with 3 gradient terms in equation (4) and g21 , g22 , and g23 are the coefficients present in equation (5). Eliminating equations above, one can obtain

dP0 dx

in terms of

from the

dC0 : dx

dP0 g23 g12 − g13 g22 dC0 = . dx g11 g22 − g21 g12 dx

8

dµ+ dx

ACS Paragon Plus Environment

(9)

Page 9 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

As discussed in Section 3.2 in Part I of this two-paper series, we can write the virtual total pressure (P0 ) as the superposition of virtual hydrodynamic pressure (Ph0 ) and virtual osmotic pressure (2C0 ). The following relation is eventually obtained between the virtual hydrodynamic pressure gradient and the virtual osmotic pressure gradient: dPh0 1 g23 g12 − g13 g22 d(2C0 ) =( + 1) , dx 2 g11 g22 − g21 g12 dx

(10)

where the pre-factor is only a function of electro-hydrodynamic coupling parameter (κ) and tabulated area-averaged coefficients. Here we considered κ = 0.5, which is relevant to aqueous electrolytes. Figure (1b) depicts the variation of this pre-factor with respect to λ∗ and different values of σ ∗ for the differential element located in an arbitrary longitudinal position inside the pore. The trend indicates the osmotic pressure is recovered as the pore EDLs become highly overlapped (thinner pore).Note that here we only present the ratio of hydrodynamic pressure gradient to the osmotic pressure locally. However, if the local pressure gradient is equal to the osmotic pressure for all the axial positions in the pore (i.e. dPh0 dx

=

d(2C0 ) ), dx

one can conclude that the pressure difference between the end reservoirs is

also equal to the osmotic pressure difference.

9

ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 10 of 37

Page 11 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

but also quantify the non-ideal behavior when the pore size is finite. Additionally, the present result indicates a tradeoff between pore size and pore charge, which can be utilized in designing perm-selective material. For example, Figure (1b) depicts that for highly charged pores, one can achieve roughly the same performance as in a low-charge pore, by choosing a pore size that is about 10 times thicker. This analysis provides quantitative guidelines for designing new membrane material that have lower viscous friction while maintaining ideal osmotic pressure. Figure (1b) reveals the existence of negative induced pressure gradients for intermediate values of λ∗ and some values of σ ∗ . This anomalous phenomenon is not predictable by the van’t Hoff osmotic theory and has been reported in previous literatures 46–52 and denoted as negative osmosis. The existence of electrostatic effect in charged membranes is one primary element of producing negative osmosis. 47,50 It is demonstrated that negative osmosis can be found if an electrolyte diffuse through a leaky membrane, and thus pore sizes must not be in the semi-permeabe limit (ideal membrane). 49,50 This is consistent with our model prediction since the negative osmosis occurs for some intermediate values of λ∗ < 1.

4.2

Deionization shock propagation through a micro-pore

To validate our numerical solution, we compared the numerical results from our reduced model with multi-dimensional direct numerical simulations of single micro-pore implemented by Nielsen et al. 32 As shown in Figure (2a), the geometry of interest is a circular pore with normalized axial length. At the right boundary the pore is blocked by an ideal cationselective membrane, which imposes zero area-averaged flux of negative ion species as well as zero net flow, u¯ = 0. Under this assumption the electro-chemical potential is known at the right end of the pore and is equal to the applied potential at the right boundary (−V0 ), however, the pressure is to be determined such that it results in the zero net flow through the system. At the left end, the pore is connected to a large reservoir with known uniform salt concentration, pressure, and electro-chemical potential. As shown in Figure (2a) both 11

ACS Paragon Plus Environment

Langmuir

pressure and electrochemical potential are set to e zero at the left boundary. ()*+,-'./0/1*+2/' 3/345)-/

7

8$!: ! ! !" & ! ! ! ! (! & ! ! "/./52,+5 !!! ! ! ! ! ! ! 6 7$& ! ! ;
" #"

" !" ./ "

-#&+"

" !"

./ "

" !" , "

" #"

!"#

%$$%$&

&'

(b)

(a)

Figure 8: (a) A network with two pores in parallel forming a loop subject to an external voltage. (b) contours of area-averaged anion concentration and the flow direction at steady state.

4.6

H-Junction

The last model problem is an H-junction microfluidic network. This type of geometry has been used in numerous lab-on-a-chip experiments 5 as a preconcentration device for detection of biological and chemical samples. As shown in Figure (9a), two microchannels with different lengths are connected to a cation-selective element. In practice, this element may be an ion-selective membrane 5 or a nano-fluidic element. 57 Here, we model this element as many parallel replicas of a nano-scale pore. Previous experimental results show that a stationary interface (similar to deionization

23

ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

shocks) is formed in one of the microchannels indicated as Channel 1 in Figure (9a), while the connecting microchannel indicated as Channel 2 experiences the full propagation of a depletion front. Prior to the present study, no mathematical analysis of this system has been performed in the context of deionization shocks. The question that still remains open is whether the experimentally observed interface is a deionization shock, 25 and if so, why this shock does not propagate away from the intersection, as observed in similar analyses. 25,27 We show that our reduced order model captures these effects as observed in the experimental reports and confirms that the stationary interface is indeed a deionization shock. Secondly, inspired by our network model, we present an even simpler phenomenological model to explain why the deionization shock in this system remains stationary. We model this system as a network of five pores, with two internal nodes and four terminal nodes. The two internal nodes are on the two sides of the cation-selective element. Each microchannel is modeled as a single pore connecting one of the terminal nodes to the internal node as shown in Figure (9a). All microchannels have the same surface charge density, thickness, and width, which are σm = −13mC/m2 , h = 1.4µm, and w = 50µm respectively. As shown in Figure (9a), Channel 1 and Channel 2 are both 6.5mm long, while Channel 3 and Channel 4 have lengths equal to 8.7mm. The cation-selective element is modeled as many identical nano-pores with length 90µm, uniformly distributed over the 1µm × 50µm cross-section. Furthermore, we consider that the permeable cross-section area is 40% of the total cross-section. Each individual nano-pore has diameter dn = 10nm and surface charge density, σn = −200mC/m2 . All terminal nodes represent reservoirs containing the same salt concentration of 1mM and are open to atmospheric pressure (i.e. P0 = 0 for all terminal nodes). The terminal nodes that are connected to Channel 3 and 4 have zero potential, while different potential values are applied in the terminal nodes located at the two ends of Channel 1 and 2, whose potentials are indicated by VL and VR in Figure (9a). Since all nodes representing intersection elements are physical components of microchannels, their thicknesses and surface

24

ACS Paragon Plus Environment

Page 24 of 37

Page 25 of 37

charge densities match the quantities specified for microchannels’ computational cells. Furthermore, the computational volumes of two internal nodes are set to their largest adjacent grid cell volumes, which are inside Channel 1 and Channel 3 respectively. We considered 3200 computational cells for each microchannel, and 30 elements for the cation-selective element. ∆t = 10−6 was employed to compute the results presented below. Figure (9b)

!" .. #$ 1 $ 2;