- Email: [email protected]
Open-source toolbox for electromigrative separations
Accepted Manuscript Open-source toolbox for electromigrative separations Santiago Márquez Damián, Federico Schaumburg, Pablo A. Kler
PII: DOI: Reference:
S0010-4655(18)30407-7 https://doi.org/10.1016/j.cpc.2018.11.015 COMPHY 6673
To appear in:
Computer Physics Communications
Received date : 21 April 2018 Revised date : 30 October 2018 Accepted date : 19 November 2018 Please cite this article as: S.M. Damián, F. Schaumburg and P.A. Kler, Open-source toolbox for electromigrative separations, Computer Physics Communications (2018), https://doi.org/10.1016/j.cpc.2018.11.015 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Open-source toolbox for electromigrative separations Santiago M´arquez Dami´ana,b,, Federico Schaumburgc , Pablo A. Klera,d a Centro
de Investigaci´on de M´etodos Computacionales (CIMEC, UNL–CONICET) de Ingenier´ıa Mec´anica, Facultad Regional Santa Fe, Universidad Tecnol´ogica Nacional c Instituto de Desarrollo Tecnol´ ogico para la Industria Qu´ımica (INTEC, UNL–CONICET) d Departamento de Ingenier´ıa en Sistemas de Informaci´ on, Facultad Regional Santa Fe, Universidad Tecnol´ogica Nacional b Departamento
Abstract We present a novel simulation toolbox for electromigrative problems. The formulation presented includes all the fields, interaction and effects that configure both electrokinetic and electroosmotic phenomena for studying capillary and chip electromigrative separation methods. This simulation toolbox was developed by using the Finite Volume Method (FVM) in the platform OpenFOAM R . This implementation defines a new scenario for numerical simulations of electromigrative separations in terms of four main novel and outstanding characteristics offered by OpenFOAM R : (i) native 3D support, (ii) electroosmotic fluid flow solution by using Navier-Stokes equation, (iii) automatic parallel and supercomputing support, and (iv) GNU-GPL license. In this way, we put into consideration of the electromigrative separation community, for the first time, a simulation tool for a wide range of experimental methods such as capillary zone electrophoresis, isotachophoresis, and isoelectric focusing, among others. After presenting the main characteristics of the implementation, different validation and application examples are provided in order to illustrate the capabilities and potentialities of the library, but also trying to motivate a discussion about its future. Keywords: Finite volume method, OpenFOAM R , High performance parallel simulations, Electroosmotic flow, Electrophoresis, Isoelectric focusing, Isotachophoresis. Program summary
1. Introduction
Program title: electroMicroTransport Program Files doi: http://dx.doi.org/10.17632/9wpvypzj9y.1 Electrophoresis can be defined as the differential Licensing provisions: GNU General Public License v3 movement of charged species in solution by attraction External routines: OpenFOAM R v17.12, v18.06 or repulsion under the action of an electric field [1]. Nature of problem: 3D multiphysics problems involving fluid Electrophoretic separations are widely used in many dynamic and electromigrative phenomena in capillaries and fields as the academy, industry or clinical laboratories microfluidic chips. due to their high separation efficiency, low reactive and Solution method: Navier-Stokes, charge conservation and energy consumption, portability and integrability with transport equations solved via the Finite Volume Method on-chip technologies, and its proven capability to interunning in shared and distributed memory parallel platforms. grate analytical systems with standard and non-standard
∗ Corresponding author. E-mail address: [email protected]
Preprint submitted to Computer Physics Communications
detection techniques [2]. Moreover, the main operation modes for electromigrative separations, such as capillary zone electrophoresis (CZE), isoelectric focusing (IEF) and isotachophoresis (ITP), allow processing samples with different bio-physico-chemical characteristics, while obtaining from each method complementary information [3]. Modeling and computer simulation of such techniques is mandatory for developers of analytical methods in order to avoid unnecessary experimentation while consuming time, reactants and expensive equipment operated by highly qualified personNovember 29, 2018
nel [4]. Software allowing 1D simulation of the most common electromigrative separation techniques has been available for years [5], being the most relevant GENTRANS [6], SIMUL 5 [7] and SPRESSO [8]. These tools solve a set of non-linear partial differential equations describing the charge and mass conservation principles. Thus, by setting the boundary value problem, specific 1D cases are possible to simulate. These tools were reported to successfully simulate 1D CZE, IEF and ITP problems. Moreover, they arrive to equivalent results when the same problem is solved [9]. Nonetheless, although some of them are free of charge, easily available and simple to install, they do not allow 2D or 3D modeling, and then, several relevant applications are excluded of their simulation scope. Applications that involve electroosmotic flow (EOF) with non-uniform velocity, which is the most common case found under practical laboratory conditions, represent the most significant case for this limitation. Other applications requiring 2D or 3D modeling include free flow methods, like free flow electrophoresis (FFE) and free flow IEF (FFIEF), or lab-on-a-chip (LOC) devices, where microfluidic channel turnings and crossings effects have to be accounted for. Moreover, electrophoretic separations are gaining relevance among the paper–based microfluidic community, where ITP was already implemented [10, 11]. In the past years, several improvements to named conventional models where reported, regarding more complex physicochemical models for particular experimental conditions [12, 13, 14], and for complex models of electrically-driven flows [15]. Moreover, some interesting applications like IEF performed on nanochannels [16], 2D IEF [17, 18], FFIEF [19, 20], and even FFIEF coupled to CZE in a 2D LOC [21, 22] were successfully simulated using different computational tools based in the finite element or the finite volume methods. However, none of this software is widely available, nor is easy to install/use, preventing their extensively use by the electromigrative separation community. Recently, a commercial tool was also reported to have a good performance when compared to conventional software [23]. Although 2D and 3D modeling is supported, closed packages may not be flexible enough for method development or investigation of the new formulations required by emerging applications, besides the fact that their license price can be prohibitive for developing countries. Thus, computational tools for the more challenging applications, such as 2D and 3D separations or LOC
electrophoresis, where computer simulation is intensively required, currently are not widely accessible, to the best of our knowledge. In this context, we consider the importance of the development of a computational toolbox that satisfies the following requirements: i) able to handle 2D and 3D problems, ii) capable of solving the multiscale and multiphysics phenomena involved in the electromigrative separation techniques, iii) open source and flexible enough to enable the inclusion of new elements and features that will be necessary in the future, and iv) capable of parallelization as an answer to large problems involving multiple species and complex geometries. The OpenFOAM R (Open Field Operation and Manipulation) platform is a free open source project for the solution of multiphysic problems based on the Finite Volume Method (FVM). Besides from the wellknown advantages of object oriented programming which makes code reuse an easy task, OpenFOAM R offers three main outstanding characteristics. First, native 3D support is given. Naturally, the platform offers the possibility of running 2D and 1D problems also, but it is optimized to solve numerical problems in complex 3D domains. Second, automatic parallel and supercomputing support is offered. For any user, independently of his/her programing skills, OpenFOAM R is ready to run in parallel for both shared and distributed memory computational facilities. Third, OpenFOAM R is licensed under the GNU–GPL, enabling the community to work following collaborative paradigms in order to improve and develop the library, including new physiochemical models as well as more efficient numerical approaches [24, 25]. In this work, the electroMicroTransport toolbox, based on OpenFOAM R is presented. It includes a solver application, new boundary conditions and tutorials. The toolbox is capable of dealing with electromigrative separation problems, comprising both electrokinetic and electroosmotic phenomena, involving fluid flow and arbitrary number of chemical species. Tutorials are given in the form of working examples. Three of them also operate as validation cases, and the fourth as an application example. 2. Mathematical model The solution of electromigrative separation problems implies finding the velocity, pressure, electric, and concentration fields. Then, the mathematical model involves equations for the fluid dynamics, electric fields, mass transport and the coupling of all these. Isothermal conditions are assumed, since temperature increase is 2
negligible provided that Joule heating is fully dissipated by the capillary electrophoresis or LOC device [21].
Due to the fact that the channel walls are supposed to be perfectly insulating, there are no components of the applied electric field normal to such surface. Therefore, near the interface region, Φ varies only tangentially to the wall plane and the total electric potential can be thought as the superposition of ζ and Φ, whose derivatives, i.e. the associated electric fields, are strictly orthogonal. This superposition is valid if the applied electric field is small compared to the EDL field, as often happens in practice. Introducing the electric potentials into Eq. (3) leads to the following equation:
2.1. Fluid dynamics In the framework of continuum fluid mechanics, fluid velocity u and pressure p are governed by the Navier–Stokes equations for newtonian incompressible flows [26]: ∇·u=0 (1) ! ∂u + u · ∇u = ∇ · (−pI + µ(∇u + ∇uT )) + ρel E. ρ ∂t (2)
∇2 ζ + ∇2 Φ = −
Eq. (1) is the incompressibility constraint, while Eq. (2) represents the conservation of momentum. In this case, the gravitational forces have been neglected due to the small length scale of capillaries and microdevices, and the consequently small Bond number. Here, ρ is the fluid density, µ the viscosity, and I the identity tensor. The last term on the right hand side of Eq. (2) is the contribution of the electrical forces to the momentum balance, where E is the electric field strength and ρel is the electric charge density of the electrolyte solution, obtained as the net charge summation over all ions present in the solution. In the FVM formulation presented in this work, this last term is neglected due to the reasons presented next.
(4)
The second term on the left hand side of Eq. (4) is nonzero in electromigrative problems because of the variations on ∇Φ along the channel due to the presence of different ions generating zones with different electrical conductivities. However, ∇2 Φ is several orders of magnitude smaller than ∇2 ζ (see also [29, 30, 31]), which allows splitting the computation of the electric field into two parts as explained below. 2.2.1. Electric double layer According to the previous analysis, by neglecting ∇Φ in Eq. (4), the EDL potential is governed by, ∇2 ζ = −
2.2. Electric field Charge density distribution for an electrolyte solution with solvent permittivity , is related to the electric potential ψ as [27]: ∇ · (−∇ψ) = ρel .
ρel .
ρel .
(5)
Considering that ζ decreases to zero within 1– 10 [nm] from the wall, while the cross-sectional channel or capillary dimensions are typically 10–500 [µm], computational requirements are excessively large when an entire electrophoretic device is modeled. In that sense, we simplify the EDL model and EOF computation by introducing the so-called thin EDL approximation [32, 33]. Consequently, EOF is regarded as an electrically induced slip velocity (the EOF velocity ueo ) parallel to the direction of the applied electric field (∇Φ). This velocity is determined by the Helmholtz–Smoluchowski equation,
(3)
Modeling electromigrative problems requires special considerations for the electric potential, since it involves different length scales depending on which part of the domain is analyzed: the bulk of the electrolyte solution or the walls of the capillary/microchannel [28]. The first contribution to the electric field (and consequently the electric potential) comes from the presence of electrostatic charges at the solid–liquid interfaces conforming the well-known electric double layer (EDL) [27]. The interfacial charge has an associated electric potential ζ that decreases exponentially in the direction normal to the wall due to the screening produced by counterions in the electrolyte solution. There is also another electric potential Φ, which arises when a potential difference ∆Φ is externally applied to drive the electromigration and/or induce EOF.
ueo =
ζ∇Φ . µ
(6)
This approximation also implies that ρel ≈ 0 in the fluid outside the EDL, meaning that the last term on the right-hand side of Eq. (2) is negligible, and enabling the splitting of Eq. (4) in two independent equations. Thus, the electroosmotic velocity is included in the fluid prob3
lem as a boundary condition, reducing significantly the computational demand. Finally, for the calculation of the electrokinetic potential ζ value at the interfase, we consider those frequent cases found in practice such as interfaces containing weak acid groups, such as silanol in fused silica capillaries and carboxyl in synthetic polymer materials. For such interfaces the following relationship is appropriate for symmetric monovalent background electrolytes (BGE) [34]: Fζ √ −en s 8RT I sinh = (pK 2RT 1 + 10 S −pH) e−Fζ/(RT )
where D j is the diffusion coefficient of the j−species, and σ is the electrical conductivity calculated as:
pH = 3 − log10 ([H ]) +
2.3. Mass transport and chemistry The non-stationary mass transport of a weakly concentrated sample molecule or a buffer electrolyte component j is governed by,
(7)
∂c j + ∇ · (−z j Ω j ∇Φc j + uc j − D j ∇c j ) = 0. ∂t
(8) (9)
(10)
where N is the total amount of ionic species considered, c j the concentration, and z j the effective charge of every j−species present in the electrolyte solution. Therefore, if the parameters that characterize the interface (nS , KS ) are known, the ζ-potential can be locally predicted for particular values of pH and I.
zj =
3 X
nαnj .
(14)
n=−3
where αnj is the mass fraction of the n charge state. After this consideration, z2j must be redefined to properly calculate Eqs. (10) and (12) as:
2.2.2. Bulk electrolyte Based on the considerations above (ρe ≈ 0, and ∇ζ ∇Φ), the electric potential related to the external source, Φ, can be computed by using the steady state charge conservation equation [27], N X ∇ · − σ∇Φ − F z j D j ∇c j = 0
(13)
Here, a linear superposition of three transport mechanisms is considered, namely the electromigrative, advective, and diffusive phenomena. Diverse electrolytes regarding its acid-base behavior, and particularly protons coming from water dissociation, have to be considered. A suitable approximation to model the acid-base reactions is to adopt chemical equilibrium constants, since in electrolyte chemistry the association and dissociation processes are much faster than transport processes [35]. Then, on a general form, we consider seven different possible ionic states for every electrolyte1 , one neutral (z j = 0), three positives (one, two and three positive charges, i.e. z j = 1, z j = 2, and z j = 3, respectively), and three negative (one, two and three negative charges, i.e. z j = −1, z j = −2, and z j = −3, respectively). Then, the effective charge z j of any electrolyte is calculated as the linear superposition of every possible charge state affected by its relative mass fraction, i.e.:
N
1X c jz j2 2 j=1
(12)
where Ω j is the electrophoretic mobility of the j−species. The terms in brackets in Eq. (11) represent electric current density contributions, which account for the ion fluxes due to electrical forces and diffusion.
where [H+ ] represents the proton concentration in mol m−3 . Finally, I represents the ionic strength of the electrolyte solution in the surroundings of the wall, calculated as: I=
z j2Ω jc j
j=1
where, R is the ideal gas constant, T represents temperature, F the Faraday constant, and e the fundamental charge. Simultaneously, pK s and n s represent the relevant physicochemical characteristics of the wall by means of the dissociation point of the representative charged groups, and its surface density, respectively. Here, both pK s and pH are considered in the classical way for physical chemistry, i.e.: pK s = 3 − log10 (K s )
N X
σ=F
z2j =
3 X
(nαnj )2 .
(15)
n=−3
(11)
1 More ionic states can be considered, but they can be though as exceptional cases, justifying an ad-hoc modification of the code.
j=1
4
zero for Forward Euler scheme, one for Backward Euler Scheme and a half for the Crank-Nicolson scheme. The sum over f implies the sum of fluxes over all the faces of the faceted cell j, while S f is a vector normal to face f , with a magnitude proportional to its area.
Clearly, the effective charge of every electrolyte depends on the pH of the surrounding solution and the different pK characteristic of every charge transition. Then, the relative amount of every charged state is calculated as: αkj = P3
[H+ ]k+3 Π−k j=−3 K j
i=−3
[H+ ]i+3 Π−i j=−3 K j
(16)
3.2. Discretization Schemes The diffusive, advective, and source terms require further discretization. For the first one, central difference is used with non-orthogonal correction. For advection, TVD methods [39] are used to ensure high order resolution along bounded and non-oscillatory solutions. Here is worth to note that advective velocities can be due to fluid flow as well as electromigration. All the schemes available in OpenFOAM R are ready to use. In the application examples provided within the toolbox the Minmod [40] scheme was selected. Finally, source terms include constant and reactive terms which are treated semi-implicitly by linearization. The whole discretization theoretically reaches second-order convergence in space and time for a single and totally uncoupled scalar [37]. For the solution of fluid flow a special method is applied to couple pressure and velocity, since pressure has no evolutive equation for incompressible flows. Therefore, the Pressure Implicit Split of Operators method (PISO) [41] is applied to solve Eq. (2) by means of nCorrectors PISO loops. Regarding Eqs. (7) and (17), these are solved by using standard NewtonRaphson loops with analytical jacobians. The implementation of such methods are available in the source code, more precisely in SmoluchowskyWallVelocityFvPatchVectorField.C and ampholyteSet.C, respectively. The linear systems resulting from the discretization procedure are solved by using a Pre-conditioned BiConjugate Gradient (PBiCG) method coupled with a Diagonal Incomplete-LU (DILU) pre-conditioner, for the non-symmetric matrices that result from the discretization of fluid flow and mass transport equations (Eqs. (1), (2), and (13)). In the case of the symmetric matrix that results from the discretization of the charge conservation equation (Eq. (11)), a Preconditioned Conjugate Gradient (PCG) method coupled with a Diagonal Incomplete-Cholesky (DIC) preconditioner, was used. The absolute tolerance for the solvers is typically set at 10−12 for Eq. (13) and 10−9 for the other cases. More detailed information about solvers and discretization schemes can be found at the source code of the tutorials on the corresponding ./system/fvSolution and ./system/fvSchemes files, respectively.
where Ki represents the dissociation constant for the icharge state of the ion. Finally, in order to calculate the pH of the solution, we rely on the already mentioned hypothesis that in the bulk computational domain the electroneutrality verifies (ρe ≈ 0). Thus, a local sum over all the charges present has to be null. Then a non-linear equation can be solved in order to determine the pH on every domain cell, i.e.: N
[H+ ] −
X Kw z jc j = 0 + + [H ] j=1
where Kw is the water Kw =1 × 10−8 mol m−3 .
dissociation
(17) constant
3. Numerical modeling 3.1. Finite volume method formulation The model equations presented in the previous sections are solved by means of the CellCentered Finite Volume Method [36, 37, 38] over the OpenFOAM R platform [24], allowing for rapid code development, re-usability and collaborative work due its GNU-GPL license. This platform was designed to solve general transport problems from diverse background physics, using a native three-dimensional treatment of the problems and ready-to-use parallelization algorithms. Therefore Eqs. (2), (11) and (13) are treated as advection–diffusion equations with source terms, leading the discrete equation for cell j, cn+1 − cnj j ∆t
+
X f
γ∇cn+θ + v cn+θ · S f = S V j + b j
(18)
where this equation corresponds to the discretization of an scalar problem for the quantity c, with diffusivity γ, advective velocity v and source S . Additionally b is the source term due to boundary conditions, ∆t is the timestep and θ is a temporal implicitness parameter being 5
changes rapidly and then stays constant during time until next time-step. Then two time-step are set, the main one ∆t, and a secondary one ∆t s :
3.3. Numerical details Next, a series of details about numerical methods are explained, particularly those which are considered of key importance to obtain successful and efficient simulations.
∆t s =
3.3.1. Electric double layer boundary condition Since it is not computationally feasible to capture the complete physics of the EDL, this effect is included in the model by means of a boundary condition for the velocity. It is set as Dirichlet boundary condition where velocity is computed by the Helmholtz-Smoluchowsky equation as explained in Section 2.2.1.
∆Φ ∆l
Φa = Φa − Φc = ∆Φ =
(19)
where nUSC is the number of subcycles in which the main time-step is divided. Therefore the hydrodynamic flow is solved a given number of ∆t s time-steps (usually nUSC /1000) until reaching steady-state and then all species are transported with this new flow using ∆t. 3.3.5. Courant number calculation and adjustable time step The maximum Courant number is computed for both effective flux φ and hydrodynamic flux φh . The first one is the maximum of all species. The calculation is implemented as follows:
3.3.2. Constant current boundary condition Physics in section 2.2 and OpenFOAM R are described and implemented in terms of the electric potential Φ. Nevertheless, many practical problems are described in terms of electric current density i = −σ∇Φ. Thus, through the use of Φ, a Dirichlet boundary condition for prescribing i was implemented for 1D cases. The uniform1DCurrentDensity boundary condition averages σ values in every cell, for every time step, allowing to fix i according to: i = σav
∆t nUSC
Comax
)) ( P 1 f φf ∆t = max max n j 2 Vi (
(20)
where a maximum over all cells n, and species j, is implied. The global time-step ∆t is adjusted in order to fullfil the Courant number restrictions for the momentum-equations and each ampholyte transport equation. Therefore an adaptive time-step control is coded based on the combination of fluid and electrodriven velocities and a maximum allowed Courant number maxCo, based on the effective flux φ.
i∆l σav
where ∆l is the distance between the anode and cathode, Φa is the potential at the anode and Φc , the potential in the cathode, is assumed to be zero.
F∆t,max =
3.3.3. Efective flux and mixed boundary conditions The standard inletOutlet boundary condition is used as usual but redefining the hydrodynamic flux φ, in order to account not only for the fluid flow flux φh = u · S f , but also for the electromigrative flux φe = (−z j Ω j ∇Φ) · S f , being φ = φh + φe and S f an outward vector normal to each boundary face. With this new definition, if φ ≥ 0, a zero gradient is imposed on that face; on the other hand, if φ < 0 the boundary is recognized as an inlet and a c f value is prescribed, giving a mixed Dirichlet-Neumann boundary condition.
maxCo max{Comax,φ ,Comax,φh }
F∆t = min min F∆t,max , 1 + 0.1 F∆t,max , 1.2
(21)
∆t = min {F∆t ∆t, ∆tmax }
3.3.6. Solution algorithm The set of equations involved in the problem are solved sequentially in a segregated fashion: 1. Start program and load physical and numerical parameters and do the following steps the required number of time-steps
3.3.4. Dual time-stepping Time-scales for hydrodynamic and electromigrative processes are completely different. Species evolve in long time-steps meanwhile the hydrodynamic flow
2. Correct H + by means of Eq. (17) 3. Solve for Φ, Eq. (11) 6
The next level in code organization is the class ampholyteSet, whose main objective is to solve Eq. (13), for each ampholyte and provide group properties as the ionic strenght I and the global conductivity σ required to solve Eq. (11). Due to the modularity of OpenFOAM R the boundary conditions are selected at run-time allowing great flexibility at user level. Therefore, the SmoluchowskyWallVelocity boundary condition (see section 3.3.1) is implemented in the SmoluchowskyWallVelocityFvPatchVectorField class while the uniform1DCurrentDensity is implemented in the homonymous class. Finally, the solver electroMicroTransport is coded in the main loop, solving Eqs. (11), (2) and (13) sequentially in a segregated fashion as explained in section 3.3.6.
4. Use time sub-looping in p-u solving, dividing the global time-step by nUSubCycles and solving the coupling by nEndUSubCycles sub-steps: 5. Start the p-u coupling loop iterating nOuterCorrectors (PIMPLE loop) times a. Compute prediction for u, Eq. (2) b. Solve the PISO loop nCorrectors times coupling p and u, Eq. (1) and (2) 6. End PIMPLE loop 7. End p-u sub-looping 8. Solve sequentially Eq. (13) for all the c j species using global ∆t. 4. electroMicroTransport toolbox
4.2. Boundary conditions
As stated above, the electroMicroTransport toolbox based on OpenFOAM R comprises of the solver itself, two boundary conditions and four tutorials, as schematized in Fig.1.
Two new boundary conditions were implemented from scratch. The first modifies the fluid dynamic field, i.e. the SmoluchowskyWallVelocity boundary condition, which prescribes a velocity due to the electroosmotic phenomenon as explained in section 3.3.1. In order to setup this boundary condition no further parameters are required, since the SmoluchowskyWallVelocityFvPatchVectorField class obtains the local velocity values to be imposed as Dirichlet conditions from the local values of ∇Φ, I, and pH fields. Second, the uniform1DCurrentDensity boundary condition for the electromigrative phenomena, explained in section 3.3.2, was implemented. It is only valid for 1D problems and it prescribes a given current density along the domain. This should be indicated in the 0/Phi file, together with a characteristic length L, which is the distance between the electrodes. 4.3. Tutorials
Figure 1: Components of the electroMicroTransport toolbox.
Four tutorials are provided to exemplify the electroMicroTransport usage. These tutorials are working examples which will be explained next in sections 5 and 6. The first three tutorials address wellknown 1D validation cases, while the fourth is a 2D application example specially chosen to show the tool potential.
4.1. Code organization The code is organized in three classes plus the code for the solver itself. The base class is ampholyte which reads the initial concentration distribution for each ampholyte j as long as its physical properties such as diffusivity D j , electrophoretic mobility Ω j and pK j values. This class provides methods to return these parameters, and also the quantities z j , z2j and ∂z j /∂H + required to solve the associated transport equations.
4.4. Download, installation and basic usage Since electroMicroTransport is an add-on, the first step is to download the OpenFOAM R source code from the developer’s site, and compile it following 7
Table 1: Compound concentrations in mM for each electrolyte solution for ITP example
the instructions provided [25] in a user-writable folder. Then electroMicroTransport has to be downloaded from its repository [42] to a destination within the user’s home folder. Next, for compilation and assuming that the user runs a Linux plataform, it is required to run the following command line:
Compound Acetic acid Lithium Sodium hydroxide Potassium
$ . / Allwclean $ . / Allwmake It has to be noted that the first instruction involves the overwriting of the native finiteVolume library and thus it takes a few minutes. OpenFOAM R runs on Linux and Windows platforms, then electroMicroTransport should theoretically work on both operative systems but was only tested on Ubuntu 16.04 LTS, Fedora 20 and 24, Debian 9 and CentOS 6.5 with OpenFOAM R versions 18.06 and 17.12. From then, electroMicroTransport can be used as any other OpenFOAM R application, being the easiest way to start, to copy and run one of the given tutorials, by running the following command lines:
TE [mM] 10.0 10.0 0.0 0.0
Sample [mM] 10.0 0.0 10.0 0.0
LE [mM] 10.0 0.0 0.0 10.0
Profiles obtained show the characteristic triangular shape expected due to electrodispersion and are practically identical to those obtained with other simulation tools [9, 23]. 5.2. Isotachophoresis A typical example for isotachophoretic separation of cations is shown in this section. ITP consists of sandwiching a sample between a leading electrolyte (LE) and a terminating electrolyte (TE). The electrolyte system is designed in order to fulfill the condition that the LE coion has higher mobility, while the TE coion has lower mobility, than the sample. When an electric field is applied, a configuration of consecutive zones ordered by mobilities is formed. The whole electrolyte system migrates with the same velocity [43]. Computationally, ITP represents a challenge due to the large gradients in conductivity, particularly in the first stages of the separation, i.e. the first 1-3 seconds after applying electric field. Table 1 gathers the data of the electrolyte system studied in this section. Here, a constant 2000 A m−2 electrical current density is applied on a 1 cm capillary. The sample is initially placed between 0.05% and 0.9% of the total capillary length. The 1D grid consists of 8000 segments. The numerical concentration profiles of all the chemical species obtained at 4.8 s after applying the electric current is depicted in figure 3a, where the three zones (TE, sample and LE) can be clearly distinguished. It can also be seen that the model handles properly the abrupt transition between the different zones without oscillations in the case of the electrolyte concentrations and conductivity. In contrast, the pH shows small oscillatory transitions between the zones meaning that diffusion current in the transport equations is correctly included in the formulation, but it is not consistent with the stabilization terms. We can conclude this, by carefully taking into account all the considerations present in [44], as far as the
$ . / i n i c o n d . sh $ . / electroMicroTransport
5. Validation examples The three tutorials provided within the electroMicroTransport toolbox serve also as validation examples, since they constitute well-known cases previously used as benchmarks for the evaluation of other 1D models [23]. 5.1. Capillary Zone Electrophoresis In CZE charged chemical species are separated based on differences in their electrophoretic mobilities, in a uniform background electrolyte (BGE) placed in a capillary tube [1]. In this case, the BGE is a buffer composed of acetic acid 20 mM, and Tris 12 mM in a 20 cm capillary, where a 2500 A m−2 electrical current density is applied. A 5 mm length sample consisting of 1 mM pyridine and 1 mM aniline is placed 5 mm away from the anode position as seen in Fig. 2a. At the boundaries, concentration fields are managed with the inletOutlet condition. Regarding the grid, the capillary tube was divided into 16000 linear segments. The simulation was run for 240 s giving the concentration profiles for pyridine and aniline shown in Fig. 2c. 8
Figure 2: Validation results for tutorial example 5.1. a) Initial aniline and pyridine concentration profiles and separation achieved by CZE simulation after b)120 s and c) 240 s. Results obtained with GENTRANS are represented with dotted lines [23].
6. Application example: Double-L injection device
oscillatory behavior for the pH vanishes for finer grids (see supplemetary material). Precisely, it was demonstrated that this oscillatory shape is related to numerical conditions, while a hollow-like shape for the pH corresponds to those reported cases where diffusion current was neglected in the modeling.
To show the ability of the electroMicroTransport toolbox to solve more complex cases, the full operation of a microfluidic Double-L device, for pinched injection and CZE, as shown in Fig. 5, was simulated [46]. Proper representation of the problem requires EOF, electric potential and transport to be handled in 2D. A double-L injection device allows controlling plug size and sample leakage, improving detection efficiency. Different modes of operation are possible according to potentials applied to the four inlets/reservoirs [47]. The operation mode simulated involves two sequential steps known as injection and a separation. In the first step, a sample composed of glycine (Gly), alanine (Ala) and phenylalanine (Phe) 1 mM each is inserted via the inlet 2 in a 2.3 M acetic acid buffer. In the same inlet a 100 V potential is applied, while Outlet 2 is grounded. The sample is transported due to the EOF and electromigration until reaching the separation channel. In the next step, a 1000 V separation potential is applied at the inlet 1 while outlet 1 is grounded. Thus, the sample plug migrates through the main channel where electrophoretic separation occurs. Boundary conditions used are summarized in table 2. A total of 181250 cells was used, where the width of the channel was discretized with 50 elements and the grid was refined at the inlets, outlets and channel intersections. Simulation results are shown in Fig. 5. In particular, Fig. 5e shows the time evolution of the sample (c sample = cAla + cGly + cPhe ) in a virtual detector located at 3.0 cm from the injection arm. This signal can be correlated with the electropherogram obtained with such detector. Here, it can be seen that the peaks appear in the correct order: i) Glycine
5.3. Isoelectric focusing Amphoteric species are chemical compounds which can act as an acid or a base according to the pH of the electrolyte solution. The pH at which the molecule exhibits no charge is called the isoelectric point (pI). IEF consists on separating amphoteric species based on their pI, through the use of a pH gradient usually generated with a mixture of carrier ampholytes. Thus, IEF experiments require a large number of components [45]. The IEF problem solved in this paper consists of a 5 cm capillary where 35 hypothetical carrier ampholytes are uniformly distributed at 1 mM concentration each. The carrier ampholytes have pI values ranging from 3.1 to 9.9 with ∆pI = 0.2, pK+1 = pI −1 and pK−1 = pI +1. A constant 300 V cm−1 electric field is applied to the boundaries, which are closed, avoiding any kind of mass flow. The grid used for simulation consisted of 10000 linear segments. Fig. 4 shows the concentration profiles, the pH and conductivity along the capillary after 10 min. It can be seen how the carrier ampholytes order themselves with their pI increasing from anode to cathode. Also, it can be appreciated how the pH increases monotonically along the capillary. Results obtained are very similar to those obtained with other solver methods [23]. 9
Figure 3: Validation results for tutorial example 5.2. a) Cations concentration b) conductivity and pH profiles along the capillary tube, after a 4.8 s ITP simulation of the system described in 5.2. Results obtained with GENTRANS are plotted with dotted lines [23].
(tGly = 121.6 s), ii) Alanine (tAla = 127.9 s)iii) Phenylalanine (tPhe = 156.8 s). Moreover, the position of the peaks coincides with an error lower than 4% with the position calculated analytically [48], taking into consideration the position of the detector, the applied field, the length of the main channel, the effective charge of each species and the electroosmotic flow(110 µm s−1 ) at the working pH (pH = 2.2). The simulation performed allows, among other things, the assessment of the peak longitudinal dispersion, and the magnitude of the sample leakage produced by the injection method. It is worth to mention here that the 2D simulation enables the user to see the tailing effect for the peaks. This effect, that is typical from these injection procedures made with microfluicic chips [49], is clearly captured by the simulation, emphasizing the difference between the presented toolbox and the other 1D tools previously developed for electrophoretic separations.
i) a boundary condition implementing the HelmoltzSmoluchowsky equation Eq. (6) for EOF modeling via the thin EDL approximation, ii) a boundary condition for Φ for fixing the current density i on single column problems, iii) the use of the native inletOutlet boundary condition accounting not only for the fluid flow flux φh , but also for the electromigrative flux φe , iv) a dual time-stepping algorithm to deal with the different hydrodynamic and electromigrative time-scales, and v) Courant number calculation involving the effective φ and hydrodynamic φh fluxes for every cell and species. Moreover, three examples are given for both tutorial and validation purposes. Furthermore, a final tutorial example of a complex 2D real problem is included to demonstrate the toolbox capabilities. Finally, the flexibility of the presented tool endows it with the possibility to easily expand it to new applications, like, for example, paper electrophoresis.
7. Concluding remarks
Acknowledgments We want to acknowledge Wolfgang Thormann for his valuable contribution to this work. This research was supported by CONICET, ANPCyT (Grant PICT 20160640), and UTN (Grant PID ASUTNFE0004475), Argentina.
A toolbox for the solution of 3D multiphysic problems involving fluid dynamic and electromigrative phenomena in capillaries and microfluidic chips was developed. This toolbox was implemented by using OpenFOAM R meaning that native 3D problem handling, parallel computing and GNU-GPL license are features supported from scratch. Also, the electroMicroTransport toolbox implements a number of novel characteristics that makes it particularly interesting for the scientific community, namely:
Supplementary material Supplementary material associated to different numerical tests performed with the toolbox and a technical 10
Figure 4: Validation results for tutorial example 5.3. a) Ampholyte concentration, b) conductivity and pH profiles along the capillary tube, achieved with a 10 min simulation of a IEF separation. The gray lines are the pH and conductivity counterparts obtained with GENTRANS [23].
discussion about numerical convergence can be found in the on-line version of the article.
[12] Q. Mao, J. Pawliszyn, W. Thormann, Dynamics of capillary isoelectric focusing in the absence of fluid flow: High-resolution computer simulation and experimental validation with whole column optical imaging, Anal. Chem. 72 (21) (2000) 5493– 5502. [13] W. Thormann, J. Caslavska, R. Mosher, Modeling of electroosmotic and electrophoretic mobilization in capillary and microchip isoelectric focusing, J. Chromatogr. A 1155 (2) (2007) 154–163. [14] R. A. Mosher, W. Thormann, High-resolution computer simulation of the dynamics of isoelectric focusing: In quest of more realistic input parameters for carrier ampholytes, Electrophoresis 29 (5) (2008) 1036–1047. [15] F. Pimenta, M. A. Alves, Numerical simulation of electricallydriven flows using openfoam, arXiv preprint arXiv:1802.02843. [16] W.-L. Hsu, D. W. Inglis, M. A. Startsev, E. M. Goldys, M. R. Davidson, D. J. Harvie, Isoelectric focusing in a silica nanofluidic channel: effects of electromigration and electroosmosis, Anal. Chem. 86 (17) (2014) 8711–8718. [17] A. Chatterjee, Generalized numerical formulations for multiphysics microfluidics-type applications, J. Micromech. Microeng. 13 (5) (2003) 758. [18] J. Shim, P. Dutta, C. F. Ivory, Modeling and simulation of ief in 2-d microgeometries, Electrophoresis 28 (4) (2007) 572–586. [19] K. Yoo, J. Shim, J. Liu, P. Dutta, Mathematical and numerical model to study two-dimensional free flow isoelectric focusing, Biomicrofluidics 8 (3) (2014) 034111. [20] K. Yoo, J. Shim, J. Liu, P. Dutta, Efficient algorithm for simulation of isoelectric focusing, Electrophoresis 35 (5) (2014) 638– 645. [21] P. Kler, C. Berli, F. Guarnieri, Modeling and high performance simulation of electrophoretic techniques in microfluidic chips, Microfluid. Nanofluid. 10 (1) (2011) 187. [22] P. Kler, L. Dalcin, R. Paz, T. Tezduyar, Supg and discontinuitycapturing methods for coupled fluid mechanics and electrochemical transport problems, Comput. Mech. 51 (2) (2013) 171. [23] S. Mikkonen, H. Ekstr¨om, W. Thormann, High-resolution dynamic computer simulation of electrophoresis using a multi-
References [1] J. P. Landers, Handbook of capillary and microchip electrophoresis and associated microtechniques, CRC press, 2007. [2] P. A. Kler, D. Sydes, C. Huhn, Column–coupling strategies for multidimensional electrophoretic separation techniques, Anal. Bioanal. Chem. 407 (1) (2015) 119–138. [3] M. Ramos-Pay´an, J. A. Oca˜na-Gonzalez, R. M. Fern´andez´ Bello-L´opez, Recent trends in capTorres, A. Llobera, M. Angel illary electrophoresis for complex samples analysis: A review, Electrophoresis (2017) 111–125. [4] W. Thormann, J. Caslavska, M. C. Breadmore, R. A. Mosher, Dynamic computer simulations of electrophoresis: Three decades of active research, Electrophoresis 30 (2009) S16–S26. [5] M. Bier, O. Palusinski, R. Mosher, D. Saville, Electrophoresis: Mathematical modeling and computer simulation, Science 219 (4590) (1983) 1281–1287. [6] W. Thormann, J. Caslavska, R. Mosher, Modeling of electroosmotic and electrophoretic mobilization in capillary and microchip isoelectric focusing, J. Chromatogr. A 1155 (2) (2007) 154–163. [7] V. Hruˇska, M. Jaroˇs, B. Gaˇs, Simul 5–free dynamic simulator of electrophoresis, Electrophoresis 27 (5-6) (2006) 984–991. [8] M. Bercovici, S. K. Lele, J. G. Santiago, Open source simulation tool for electrophoretic stacking, focusing, and separation, J. Chromatogr. A 1216 (6) (2009) 1008–1018. [9] R. A. Mosher, M. C. Breadmore, W. Thormann, High-resolution electrophoretic simulations: Performance characteristics of onedimensional simulators, Electrophoresis 32 (5) (2011) 532–541. [10] B. Y. Moghadam, K. T. Connelly, J. D. Posner, Isotachophoretic preconcenetration on paper-based microfluidic devices, Anal. Chem. 86 (12) (2014) 5829–5837. [11] X. Li, L. Luo, R. M. Crooks, Low-voltage paper isotachophoresis device for dna focusing, Lab Chip 15 (20) (2015) 4090–4098.
11
Figure 5: Simulation of a Double-L injection device for electrophoretic separations. Dimensions and name of the boundaries are indicated in the main figure, while sample profile (c sample = cAla + cGly + cPhe ) at different instants is shown in inset a)–d) respectively. a) Sample entering the injection arm when t =22 s. b) Sample already arrived channel intersections and is headed to Outlet 2, 60 s after the experience began. c) and d) separation obtained in the main arm, when t =120 s, during the separation step. e) Time evolution of c sample . The first peak is due to Glycine (zΩ =−21.83 × 10−9 m2 V−1 s−1 ), the second due to Alanine (zΩ =−18.39 × 10−9 m2 V−1 s−1 ) and the third to Phenylalanine zΩ =−12.02 × 10−9 m2 V−1 s−1 ).
[24] [25] [26] [27] [28]
[29] [30]
physics software platform, J. Chromatogr. A 1532 (2018) 216– 222. H. G. Weller, G. Tabor, H. Jasak, C. Fureby, A tensorial approach to computational continuum mechanics using objectoriented techniques, Comput. Phys. 12(6) (1998) 620–631. OpenFOAM the open source CFD tool, https://www. openfoam.com/, accessed: 2018-03-19. J. Lyklema, Fundamentals of interface and colloid science: soft colloids, Vol. 5, Academic press, 2005. R. Probstein, Physicochemical Hydrodynamics. An Introduction, 2nd Edition, Wiley-Interscience, 2003. P. A. Kler, L. D. Dalcin, R. R. Paz, T. E. Tezduyar, Supg and discontinuity-capturing methods for coupled fluid mechanics and electrochemical transport problems, Computational Mechanics 51 (2) (2013) 171–185. J. M. MacInnes, Computation of reacting electrokinetic flow in microchannel geometries, Chem. Eng. Sci. 57 (21) (2002) 4539–4558. T. L. Sounart, J. C. Baygents, Lubrication theory for electroosmotic flow in a non-uniform electrolyte, J. Fluid Mech. 576 (2007) 139–172.
[31] T. J. Craven, J. M. Rees, W. B. Zimmerman, On slip velocity boundary conditions for electroosmotic flow near sharp corners, Phys. Fluid. 20 (4) (2008) 043603. [32] E. Brunet, A. Adjari, Generalized Onsager relations for electrokinetic effects in anisotropic and heterogeneous geometries, Phys. Rev. E 69 (1) (2004) 016306. [33] C. L. A. Berli, Equivalent circuit modeling of electrokinetically driven analytical microsystems, Microfluid. Nanofluid. 4 (5) (2008) 391–399. [34] C. L. A. Berli, M. Piaggio, J. Deiber, Modeling the zeta potential of silica capillaries in relation to the background electrolyte composition, Electrophoresis 24 (2003) 1587–1595. [35] I. Arnaud, J. Josserand, J. Rossier, H. Girault, Finite element simulation of off-gel buffering, Electrophoresis 23 (2002) 3253– 3261. [36] J. Ferziger, M. Peric, Computational Methods for Fluid Dynamics, Vol. I, Springer-Verlag, 2002. [37] H. Jasak, Error analysis and estimation for the finite volume method with applications to fluid flows, Ph.D. thesis, Department of Mechanical Engineering Imperial College of Science, Technology and Medicine (1996).
12
[38] S. M´arquez Dami´an, An extended mixture model for the simultaneous treatment of short and long scale interfaces, Ph.D. thesis, FICH, Universidad Nacional del Litoral, Santa Fe, Argentina (2013). [39] M. Darwish, F. Moukalled, TVD schemes for unstructured grids, Int. J. Heat Mass Tran. 46 (4) (2003) 599–611. [40] P. L. Roe, Characteristic-based schemes for the Euler equations, Annual Review of Fluid Mechanics 18 (1) (1986) 337–365. [41] R. Issa, Solution of implicitly discretised fluid flow equations by operator splitting, Journal of Computational Physics 62 (1986) 40–65. [42] electoMicroTransport repository, https://gitlab.com/ santiagomarquezd/electroMicroTransport/, accessed: 2018-03-19. [43] P. A. Kler, C. Huhn, Non-aqueous electrolytes for isotachophoresis of weak bases and its application to the comprehensive preconcentration of the 20 proteinogenic amino acids in column-coupling itp/ce–ms, Anal. Bioanal. Chem. 406 (28) (2014) 7163–7174. [44] T. Sounart, J. Baygents, Simulation of electrophoretic separations: effect of numerical and molecular diffusion on ph calculations in poorly buffered systems, Electrophoresis 21 (12) (2000) 2095–2103. [45] R. A. Mosher, W. Thormann, High-resolution computer simulation of the dynamics of isoelectric focusing: In quest of more realistic input parameters for carrier ampholytes, Electrophoresis 29 (5) (2008) 1036–1047. [46] J. M. Karlinsey, Sample introduction techniques for microchip electrophoresis: a review, Anal. Chim. Acta 725 (2012) 1–13. [47] C.-H. Lin, R.-J. Yang, C.-H. Tai, C.-Y. Lee, L.-M. Fu, Double-l injection technique for high performance capillary electrophoresis detection in microfluidic chips, J. Micromech. Microeng. 14 (4) (2004) 639. [48] P. Coufal, J. Zuska, T. van de Goor, V. Smith, B. Gaˇs, Separation of twenty underivatized essential amino acids by capillary zone electrophoresis with contactless conductivity detection, Electrophoresis 24 (4) (2003) 671–677. [49] D. Sydes, P. A. Kler, M. Hermans, C. Huhn, Zero-deadvolume interfaces for two–dimensional electrophoretic separations, Electrophoresis 37 (22) (2016) 3020–3024.
Table 2: Set of boundary conditions used to simulate the double-L injection device presented in section 6.
Boundary
Condition type
Equation: Transport Walls zeroGradient inletOutlet Inlet 1 Inlet 2
inletOutlet
Outlet 1
inletOutlet
Outlet 2
inletOutlet
Value
–c j =0 mM –cAcet =2300 mM –cAla =0.01 mM –cGly =0.01 mM –cPhe =0.01 mM –cAcet =2300 mM –c j =0 mM –cAcet =2300 mM –c j =0 mM –cAcet =2300 mM
Equation: Charge conservation (injection step) Walls zeroGradient Inlet 1 zeroGradient Inlet 2 fixedValue –Φ =100 V Outlet 1 zeroGradient Outlet 2 fixedValue –Φ =0 V Equation: Charge conservation (separation step) Walls zeroGradient Inlet 1 fixedValue –Φ =1000 V Inlet 2 zeroGradient fixedValue –Φ =0 V Outlet 1 zeroGradient Outlet 2 zeroGradient Equation: Fluid-Flow Walls SmoluchowskyWallVelocity zeroGradient –p Inlet 1 zeroGradient –u fixedValue –p =0 Pa zeroGradient –u Inlet 2 fixedValue –p =0 Pa Outlet 1 zeroGradient –u fixedValue –p =0 Pa Outlet 2 zeroGradient –u fixedValue –p =0 Pa
13
PII: DOI: Reference:
S0010-4655(18)30407-7 https://doi.org/10.1016/j.cpc.2018.11.015 COMPHY 6673
To appear in:
Computer Physics Communications
Received date : 21 April 2018 Revised date : 30 October 2018 Accepted date : 19 November 2018 Please cite this article as: S.M. Damián, F. Schaumburg and P.A. Kler, Open-source toolbox for electromigrative separations, Computer Physics Communications (2018), https://doi.org/10.1016/j.cpc.2018.11.015 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Open-source toolbox for electromigrative separations Santiago M´arquez Dami´ana,b,, Federico Schaumburgc , Pablo A. Klera,d a Centro
de Investigaci´on de M´etodos Computacionales (CIMEC, UNL–CONICET) de Ingenier´ıa Mec´anica, Facultad Regional Santa Fe, Universidad Tecnol´ogica Nacional c Instituto de Desarrollo Tecnol´ ogico para la Industria Qu´ımica (INTEC, UNL–CONICET) d Departamento de Ingenier´ıa en Sistemas de Informaci´ on, Facultad Regional Santa Fe, Universidad Tecnol´ogica Nacional b Departamento
Abstract We present a novel simulation toolbox for electromigrative problems. The formulation presented includes all the fields, interaction and effects that configure both electrokinetic and electroosmotic phenomena for studying capillary and chip electromigrative separation methods. This simulation toolbox was developed by using the Finite Volume Method (FVM) in the platform OpenFOAM R . This implementation defines a new scenario for numerical simulations of electromigrative separations in terms of four main novel and outstanding characteristics offered by OpenFOAM R : (i) native 3D support, (ii) electroosmotic fluid flow solution by using Navier-Stokes equation, (iii) automatic parallel and supercomputing support, and (iv) GNU-GPL license. In this way, we put into consideration of the electromigrative separation community, for the first time, a simulation tool for a wide range of experimental methods such as capillary zone electrophoresis, isotachophoresis, and isoelectric focusing, among others. After presenting the main characteristics of the implementation, different validation and application examples are provided in order to illustrate the capabilities and potentialities of the library, but also trying to motivate a discussion about its future. Keywords: Finite volume method, OpenFOAM R , High performance parallel simulations, Electroosmotic flow, Electrophoresis, Isoelectric focusing, Isotachophoresis. Program summary
1. Introduction
Program title: electroMicroTransport Program Files doi: http://dx.doi.org/10.17632/9wpvypzj9y.1 Electrophoresis can be defined as the differential Licensing provisions: GNU General Public License v3 movement of charged species in solution by attraction External routines: OpenFOAM R v17.12, v18.06 or repulsion under the action of an electric field [1]. Nature of problem: 3D multiphysics problems involving fluid Electrophoretic separations are widely used in many dynamic and electromigrative phenomena in capillaries and fields as the academy, industry or clinical laboratories microfluidic chips. due to their high separation efficiency, low reactive and Solution method: Navier-Stokes, charge conservation and energy consumption, portability and integrability with transport equations solved via the Finite Volume Method on-chip technologies, and its proven capability to interunning in shared and distributed memory parallel platforms. grate analytical systems with standard and non-standard
∗ Corresponding author. E-mail address: [email protected]
Preprint submitted to Computer Physics Communications
detection techniques [2]. Moreover, the main operation modes for electromigrative separations, such as capillary zone electrophoresis (CZE), isoelectric focusing (IEF) and isotachophoresis (ITP), allow processing samples with different bio-physico-chemical characteristics, while obtaining from each method complementary information [3]. Modeling and computer simulation of such techniques is mandatory for developers of analytical methods in order to avoid unnecessary experimentation while consuming time, reactants and expensive equipment operated by highly qualified personNovember 29, 2018
nel [4]. Software allowing 1D simulation of the most common electromigrative separation techniques has been available for years [5], being the most relevant GENTRANS [6], SIMUL 5 [7] and SPRESSO [8]. These tools solve a set of non-linear partial differential equations describing the charge and mass conservation principles. Thus, by setting the boundary value problem, specific 1D cases are possible to simulate. These tools were reported to successfully simulate 1D CZE, IEF and ITP problems. Moreover, they arrive to equivalent results when the same problem is solved [9]. Nonetheless, although some of them are free of charge, easily available and simple to install, they do not allow 2D or 3D modeling, and then, several relevant applications are excluded of their simulation scope. Applications that involve electroosmotic flow (EOF) with non-uniform velocity, which is the most common case found under practical laboratory conditions, represent the most significant case for this limitation. Other applications requiring 2D or 3D modeling include free flow methods, like free flow electrophoresis (FFE) and free flow IEF (FFIEF), or lab-on-a-chip (LOC) devices, where microfluidic channel turnings and crossings effects have to be accounted for. Moreover, electrophoretic separations are gaining relevance among the paper–based microfluidic community, where ITP was already implemented [10, 11]. In the past years, several improvements to named conventional models where reported, regarding more complex physicochemical models for particular experimental conditions [12, 13, 14], and for complex models of electrically-driven flows [15]. Moreover, some interesting applications like IEF performed on nanochannels [16], 2D IEF [17, 18], FFIEF [19, 20], and even FFIEF coupled to CZE in a 2D LOC [21, 22] were successfully simulated using different computational tools based in the finite element or the finite volume methods. However, none of this software is widely available, nor is easy to install/use, preventing their extensively use by the electromigrative separation community. Recently, a commercial tool was also reported to have a good performance when compared to conventional software [23]. Although 2D and 3D modeling is supported, closed packages may not be flexible enough for method development or investigation of the new formulations required by emerging applications, besides the fact that their license price can be prohibitive for developing countries. Thus, computational tools for the more challenging applications, such as 2D and 3D separations or LOC
electrophoresis, where computer simulation is intensively required, currently are not widely accessible, to the best of our knowledge. In this context, we consider the importance of the development of a computational toolbox that satisfies the following requirements: i) able to handle 2D and 3D problems, ii) capable of solving the multiscale and multiphysics phenomena involved in the electromigrative separation techniques, iii) open source and flexible enough to enable the inclusion of new elements and features that will be necessary in the future, and iv) capable of parallelization as an answer to large problems involving multiple species and complex geometries. The OpenFOAM R (Open Field Operation and Manipulation) platform is a free open source project for the solution of multiphysic problems based on the Finite Volume Method (FVM). Besides from the wellknown advantages of object oriented programming which makes code reuse an easy task, OpenFOAM R offers three main outstanding characteristics. First, native 3D support is given. Naturally, the platform offers the possibility of running 2D and 1D problems also, but it is optimized to solve numerical problems in complex 3D domains. Second, automatic parallel and supercomputing support is offered. For any user, independently of his/her programing skills, OpenFOAM R is ready to run in parallel for both shared and distributed memory computational facilities. Third, OpenFOAM R is licensed under the GNU–GPL, enabling the community to work following collaborative paradigms in order to improve and develop the library, including new physiochemical models as well as more efficient numerical approaches [24, 25]. In this work, the electroMicroTransport toolbox, based on OpenFOAM R is presented. It includes a solver application, new boundary conditions and tutorials. The toolbox is capable of dealing with electromigrative separation problems, comprising both electrokinetic and electroosmotic phenomena, involving fluid flow and arbitrary number of chemical species. Tutorials are given in the form of working examples. Three of them also operate as validation cases, and the fourth as an application example. 2. Mathematical model The solution of electromigrative separation problems implies finding the velocity, pressure, electric, and concentration fields. Then, the mathematical model involves equations for the fluid dynamics, electric fields, mass transport and the coupling of all these. Isothermal conditions are assumed, since temperature increase is 2
negligible provided that Joule heating is fully dissipated by the capillary electrophoresis or LOC device [21].
Due to the fact that the channel walls are supposed to be perfectly insulating, there are no components of the applied electric field normal to such surface. Therefore, near the interface region, Φ varies only tangentially to the wall plane and the total electric potential can be thought as the superposition of ζ and Φ, whose derivatives, i.e. the associated electric fields, are strictly orthogonal. This superposition is valid if the applied electric field is small compared to the EDL field, as often happens in practice. Introducing the electric potentials into Eq. (3) leads to the following equation:
2.1. Fluid dynamics In the framework of continuum fluid mechanics, fluid velocity u and pressure p are governed by the Navier–Stokes equations for newtonian incompressible flows [26]: ∇·u=0 (1) ! ∂u + u · ∇u = ∇ · (−pI + µ(∇u + ∇uT )) + ρel E. ρ ∂t (2)
∇2 ζ + ∇2 Φ = −
Eq. (1) is the incompressibility constraint, while Eq. (2) represents the conservation of momentum. In this case, the gravitational forces have been neglected due to the small length scale of capillaries and microdevices, and the consequently small Bond number. Here, ρ is the fluid density, µ the viscosity, and I the identity tensor. The last term on the right hand side of Eq. (2) is the contribution of the electrical forces to the momentum balance, where E is the electric field strength and ρel is the electric charge density of the electrolyte solution, obtained as the net charge summation over all ions present in the solution. In the FVM formulation presented in this work, this last term is neglected due to the reasons presented next.
(4)
The second term on the left hand side of Eq. (4) is nonzero in electromigrative problems because of the variations on ∇Φ along the channel due to the presence of different ions generating zones with different electrical conductivities. However, ∇2 Φ is several orders of magnitude smaller than ∇2 ζ (see also [29, 30, 31]), which allows splitting the computation of the electric field into two parts as explained below. 2.2.1. Electric double layer According to the previous analysis, by neglecting ∇Φ in Eq. (4), the EDL potential is governed by, ∇2 ζ = −
2.2. Electric field Charge density distribution for an electrolyte solution with solvent permittivity , is related to the electric potential ψ as [27]: ∇ · (−∇ψ) = ρel .
ρel .
ρel .
(5)
Considering that ζ decreases to zero within 1– 10 [nm] from the wall, while the cross-sectional channel or capillary dimensions are typically 10–500 [µm], computational requirements are excessively large when an entire electrophoretic device is modeled. In that sense, we simplify the EDL model and EOF computation by introducing the so-called thin EDL approximation [32, 33]. Consequently, EOF is regarded as an electrically induced slip velocity (the EOF velocity ueo ) parallel to the direction of the applied electric field (∇Φ). This velocity is determined by the Helmholtz–Smoluchowski equation,
(3)
Modeling electromigrative problems requires special considerations for the electric potential, since it involves different length scales depending on which part of the domain is analyzed: the bulk of the electrolyte solution or the walls of the capillary/microchannel [28]. The first contribution to the electric field (and consequently the electric potential) comes from the presence of electrostatic charges at the solid–liquid interfaces conforming the well-known electric double layer (EDL) [27]. The interfacial charge has an associated electric potential ζ that decreases exponentially in the direction normal to the wall due to the screening produced by counterions in the electrolyte solution. There is also another electric potential Φ, which arises when a potential difference ∆Φ is externally applied to drive the electromigration and/or induce EOF.
ueo =
ζ∇Φ . µ
(6)
This approximation also implies that ρel ≈ 0 in the fluid outside the EDL, meaning that the last term on the right-hand side of Eq. (2) is negligible, and enabling the splitting of Eq. (4) in two independent equations. Thus, the electroosmotic velocity is included in the fluid prob3
lem as a boundary condition, reducing significantly the computational demand. Finally, for the calculation of the electrokinetic potential ζ value at the interfase, we consider those frequent cases found in practice such as interfaces containing weak acid groups, such as silanol in fused silica capillaries and carboxyl in synthetic polymer materials. For such interfaces the following relationship is appropriate for symmetric monovalent background electrolytes (BGE) [34]: Fζ √ −en s 8RT I sinh = (pK 2RT 1 + 10 S −pH) e−Fζ/(RT )
where D j is the diffusion coefficient of the j−species, and σ is the electrical conductivity calculated as:
pH = 3 − log10 ([H ]) +
2.3. Mass transport and chemistry The non-stationary mass transport of a weakly concentrated sample molecule or a buffer electrolyte component j is governed by,
(7)
∂c j + ∇ · (−z j Ω j ∇Φc j + uc j − D j ∇c j ) = 0. ∂t
(8) (9)
(10)
where N is the total amount of ionic species considered, c j the concentration, and z j the effective charge of every j−species present in the electrolyte solution. Therefore, if the parameters that characterize the interface (nS , KS ) are known, the ζ-potential can be locally predicted for particular values of pH and I.
zj =
3 X
nαnj .
(14)
n=−3
where αnj is the mass fraction of the n charge state. After this consideration, z2j must be redefined to properly calculate Eqs. (10) and (12) as:
2.2.2. Bulk electrolyte Based on the considerations above (ρe ≈ 0, and ∇ζ ∇Φ), the electric potential related to the external source, Φ, can be computed by using the steady state charge conservation equation [27], N X ∇ · − σ∇Φ − F z j D j ∇c j = 0
(13)
Here, a linear superposition of three transport mechanisms is considered, namely the electromigrative, advective, and diffusive phenomena. Diverse electrolytes regarding its acid-base behavior, and particularly protons coming from water dissociation, have to be considered. A suitable approximation to model the acid-base reactions is to adopt chemical equilibrium constants, since in electrolyte chemistry the association and dissociation processes are much faster than transport processes [35]. Then, on a general form, we consider seven different possible ionic states for every electrolyte1 , one neutral (z j = 0), three positives (one, two and three positive charges, i.e. z j = 1, z j = 2, and z j = 3, respectively), and three negative (one, two and three negative charges, i.e. z j = −1, z j = −2, and z j = −3, respectively). Then, the effective charge z j of any electrolyte is calculated as the linear superposition of every possible charge state affected by its relative mass fraction, i.e.:
N
1X c jz j2 2 j=1
(12)
where Ω j is the electrophoretic mobility of the j−species. The terms in brackets in Eq. (11) represent electric current density contributions, which account for the ion fluxes due to electrical forces and diffusion.
where [H+ ] represents the proton concentration in mol m−3 . Finally, I represents the ionic strength of the electrolyte solution in the surroundings of the wall, calculated as: I=
z j2Ω jc j
j=1
where, R is the ideal gas constant, T represents temperature, F the Faraday constant, and e the fundamental charge. Simultaneously, pK s and n s represent the relevant physicochemical characteristics of the wall by means of the dissociation point of the representative charged groups, and its surface density, respectively. Here, both pK s and pH are considered in the classical way for physical chemistry, i.e.: pK s = 3 − log10 (K s )
N X
σ=F
z2j =
3 X
(nαnj )2 .
(15)
n=−3
(11)
1 More ionic states can be considered, but they can be though as exceptional cases, justifying an ad-hoc modification of the code.
j=1
4
zero for Forward Euler scheme, one for Backward Euler Scheme and a half for the Crank-Nicolson scheme. The sum over f implies the sum of fluxes over all the faces of the faceted cell j, while S f is a vector normal to face f , with a magnitude proportional to its area.
Clearly, the effective charge of every electrolyte depends on the pH of the surrounding solution and the different pK characteristic of every charge transition. Then, the relative amount of every charged state is calculated as: αkj = P3
[H+ ]k+3 Π−k j=−3 K j
i=−3
[H+ ]i+3 Π−i j=−3 K j
(16)
3.2. Discretization Schemes The diffusive, advective, and source terms require further discretization. For the first one, central difference is used with non-orthogonal correction. For advection, TVD methods [39] are used to ensure high order resolution along bounded and non-oscillatory solutions. Here is worth to note that advective velocities can be due to fluid flow as well as electromigration. All the schemes available in OpenFOAM R are ready to use. In the application examples provided within the toolbox the Minmod [40] scheme was selected. Finally, source terms include constant and reactive terms which are treated semi-implicitly by linearization. The whole discretization theoretically reaches second-order convergence in space and time for a single and totally uncoupled scalar [37]. For the solution of fluid flow a special method is applied to couple pressure and velocity, since pressure has no evolutive equation for incompressible flows. Therefore, the Pressure Implicit Split of Operators method (PISO) [41] is applied to solve Eq. (2) by means of nCorrectors PISO loops. Regarding Eqs. (7) and (17), these are solved by using standard NewtonRaphson loops with analytical jacobians. The implementation of such methods are available in the source code, more precisely in SmoluchowskyWallVelocityFvPatchVectorField.C and ampholyteSet.C, respectively. The linear systems resulting from the discretization procedure are solved by using a Pre-conditioned BiConjugate Gradient (PBiCG) method coupled with a Diagonal Incomplete-LU (DILU) pre-conditioner, for the non-symmetric matrices that result from the discretization of fluid flow and mass transport equations (Eqs. (1), (2), and (13)). In the case of the symmetric matrix that results from the discretization of the charge conservation equation (Eq. (11)), a Preconditioned Conjugate Gradient (PCG) method coupled with a Diagonal Incomplete-Cholesky (DIC) preconditioner, was used. The absolute tolerance for the solvers is typically set at 10−12 for Eq. (13) and 10−9 for the other cases. More detailed information about solvers and discretization schemes can be found at the source code of the tutorials on the corresponding ./system/fvSolution and ./system/fvSchemes files, respectively.
where Ki represents the dissociation constant for the icharge state of the ion. Finally, in order to calculate the pH of the solution, we rely on the already mentioned hypothesis that in the bulk computational domain the electroneutrality verifies (ρe ≈ 0). Thus, a local sum over all the charges present has to be null. Then a non-linear equation can be solved in order to determine the pH on every domain cell, i.e.: N
[H+ ] −
X Kw z jc j = 0 + + [H ] j=1
where Kw is the water Kw =1 × 10−8 mol m−3 .
dissociation
(17) constant
3. Numerical modeling 3.1. Finite volume method formulation The model equations presented in the previous sections are solved by means of the CellCentered Finite Volume Method [36, 37, 38] over the OpenFOAM R platform [24], allowing for rapid code development, re-usability and collaborative work due its GNU-GPL license. This platform was designed to solve general transport problems from diverse background physics, using a native three-dimensional treatment of the problems and ready-to-use parallelization algorithms. Therefore Eqs. (2), (11) and (13) are treated as advection–diffusion equations with source terms, leading the discrete equation for cell j, cn+1 − cnj j ∆t
+
X f
γ∇cn+θ + v cn+θ · S f = S V j + b j
(18)
where this equation corresponds to the discretization of an scalar problem for the quantity c, with diffusivity γ, advective velocity v and source S . Additionally b is the source term due to boundary conditions, ∆t is the timestep and θ is a temporal implicitness parameter being 5
changes rapidly and then stays constant during time until next time-step. Then two time-step are set, the main one ∆t, and a secondary one ∆t s :
3.3. Numerical details Next, a series of details about numerical methods are explained, particularly those which are considered of key importance to obtain successful and efficient simulations.
∆t s =
3.3.1. Electric double layer boundary condition Since it is not computationally feasible to capture the complete physics of the EDL, this effect is included in the model by means of a boundary condition for the velocity. It is set as Dirichlet boundary condition where velocity is computed by the Helmholtz-Smoluchowsky equation as explained in Section 2.2.1.
∆Φ ∆l
Φa = Φa − Φc = ∆Φ =
(19)
where nUSC is the number of subcycles in which the main time-step is divided. Therefore the hydrodynamic flow is solved a given number of ∆t s time-steps (usually nUSC /1000) until reaching steady-state and then all species are transported with this new flow using ∆t. 3.3.5. Courant number calculation and adjustable time step The maximum Courant number is computed for both effective flux φ and hydrodynamic flux φh . The first one is the maximum of all species. The calculation is implemented as follows:
3.3.2. Constant current boundary condition Physics in section 2.2 and OpenFOAM R are described and implemented in terms of the electric potential Φ. Nevertheless, many practical problems are described in terms of electric current density i = −σ∇Φ. Thus, through the use of Φ, a Dirichlet boundary condition for prescribing i was implemented for 1D cases. The uniform1DCurrentDensity boundary condition averages σ values in every cell, for every time step, allowing to fix i according to: i = σav
∆t nUSC
Comax
)) ( P 1 f φf ∆t = max max n j 2 Vi (
(20)
where a maximum over all cells n, and species j, is implied. The global time-step ∆t is adjusted in order to fullfil the Courant number restrictions for the momentum-equations and each ampholyte transport equation. Therefore an adaptive time-step control is coded based on the combination of fluid and electrodriven velocities and a maximum allowed Courant number maxCo, based on the effective flux φ.
i∆l σav
where ∆l is the distance between the anode and cathode, Φa is the potential at the anode and Φc , the potential in the cathode, is assumed to be zero.
F∆t,max =
3.3.3. Efective flux and mixed boundary conditions The standard inletOutlet boundary condition is used as usual but redefining the hydrodynamic flux φ, in order to account not only for the fluid flow flux φh = u · S f , but also for the electromigrative flux φe = (−z j Ω j ∇Φ) · S f , being φ = φh + φe and S f an outward vector normal to each boundary face. With this new definition, if φ ≥ 0, a zero gradient is imposed on that face; on the other hand, if φ < 0 the boundary is recognized as an inlet and a c f value is prescribed, giving a mixed Dirichlet-Neumann boundary condition.
maxCo max{Comax,φ ,Comax,φh }
F∆t = min min F∆t,max , 1 + 0.1 F∆t,max , 1.2
(21)
∆t = min {F∆t ∆t, ∆tmax }
3.3.6. Solution algorithm The set of equations involved in the problem are solved sequentially in a segregated fashion: 1. Start program and load physical and numerical parameters and do the following steps the required number of time-steps
3.3.4. Dual time-stepping Time-scales for hydrodynamic and electromigrative processes are completely different. Species evolve in long time-steps meanwhile the hydrodynamic flow
2. Correct H + by means of Eq. (17) 3. Solve for Φ, Eq. (11) 6
The next level in code organization is the class ampholyteSet, whose main objective is to solve Eq. (13), for each ampholyte and provide group properties as the ionic strenght I and the global conductivity σ required to solve Eq. (11). Due to the modularity of OpenFOAM R the boundary conditions are selected at run-time allowing great flexibility at user level. Therefore, the SmoluchowskyWallVelocity boundary condition (see section 3.3.1) is implemented in the SmoluchowskyWallVelocityFvPatchVectorField class while the uniform1DCurrentDensity is implemented in the homonymous class. Finally, the solver electroMicroTransport is coded in the main loop, solving Eqs. (11), (2) and (13) sequentially in a segregated fashion as explained in section 3.3.6.
4. Use time sub-looping in p-u solving, dividing the global time-step by nUSubCycles and solving the coupling by nEndUSubCycles sub-steps: 5. Start the p-u coupling loop iterating nOuterCorrectors (PIMPLE loop) times a. Compute prediction for u, Eq. (2) b. Solve the PISO loop nCorrectors times coupling p and u, Eq. (1) and (2) 6. End PIMPLE loop 7. End p-u sub-looping 8. Solve sequentially Eq. (13) for all the c j species using global ∆t. 4. electroMicroTransport toolbox
4.2. Boundary conditions
As stated above, the electroMicroTransport toolbox based on OpenFOAM R comprises of the solver itself, two boundary conditions and four tutorials, as schematized in Fig.1.
Two new boundary conditions were implemented from scratch. The first modifies the fluid dynamic field, i.e. the SmoluchowskyWallVelocity boundary condition, which prescribes a velocity due to the electroosmotic phenomenon as explained in section 3.3.1. In order to setup this boundary condition no further parameters are required, since the SmoluchowskyWallVelocityFvPatchVectorField class obtains the local velocity values to be imposed as Dirichlet conditions from the local values of ∇Φ, I, and pH fields. Second, the uniform1DCurrentDensity boundary condition for the electromigrative phenomena, explained in section 3.3.2, was implemented. It is only valid for 1D problems and it prescribes a given current density along the domain. This should be indicated in the 0/Phi file, together with a characteristic length L, which is the distance between the electrodes. 4.3. Tutorials
Figure 1: Components of the electroMicroTransport toolbox.
Four tutorials are provided to exemplify the electroMicroTransport usage. These tutorials are working examples which will be explained next in sections 5 and 6. The first three tutorials address wellknown 1D validation cases, while the fourth is a 2D application example specially chosen to show the tool potential.
4.1. Code organization The code is organized in three classes plus the code for the solver itself. The base class is ampholyte which reads the initial concentration distribution for each ampholyte j as long as its physical properties such as diffusivity D j , electrophoretic mobility Ω j and pK j values. This class provides methods to return these parameters, and also the quantities z j , z2j and ∂z j /∂H + required to solve the associated transport equations.
4.4. Download, installation and basic usage Since electroMicroTransport is an add-on, the first step is to download the OpenFOAM R source code from the developer’s site, and compile it following 7
Table 1: Compound concentrations in mM for each electrolyte solution for ITP example
the instructions provided [25] in a user-writable folder. Then electroMicroTransport has to be downloaded from its repository [42] to a destination within the user’s home folder. Next, for compilation and assuming that the user runs a Linux plataform, it is required to run the following command line:
Compound Acetic acid Lithium Sodium hydroxide Potassium
$ . / Allwclean $ . / Allwmake It has to be noted that the first instruction involves the overwriting of the native finiteVolume library and thus it takes a few minutes. OpenFOAM R runs on Linux and Windows platforms, then electroMicroTransport should theoretically work on both operative systems but was only tested on Ubuntu 16.04 LTS, Fedora 20 and 24, Debian 9 and CentOS 6.5 with OpenFOAM R versions 18.06 and 17.12. From then, electroMicroTransport can be used as any other OpenFOAM R application, being the easiest way to start, to copy and run one of the given tutorials, by running the following command lines:
TE [mM] 10.0 10.0 0.0 0.0
Sample [mM] 10.0 0.0 10.0 0.0
LE [mM] 10.0 0.0 0.0 10.0
Profiles obtained show the characteristic triangular shape expected due to electrodispersion and are practically identical to those obtained with other simulation tools [9, 23]. 5.2. Isotachophoresis A typical example for isotachophoretic separation of cations is shown in this section. ITP consists of sandwiching a sample between a leading electrolyte (LE) and a terminating electrolyte (TE). The electrolyte system is designed in order to fulfill the condition that the LE coion has higher mobility, while the TE coion has lower mobility, than the sample. When an electric field is applied, a configuration of consecutive zones ordered by mobilities is formed. The whole electrolyte system migrates with the same velocity [43]. Computationally, ITP represents a challenge due to the large gradients in conductivity, particularly in the first stages of the separation, i.e. the first 1-3 seconds after applying electric field. Table 1 gathers the data of the electrolyte system studied in this section. Here, a constant 2000 A m−2 electrical current density is applied on a 1 cm capillary. The sample is initially placed between 0.05% and 0.9% of the total capillary length. The 1D grid consists of 8000 segments. The numerical concentration profiles of all the chemical species obtained at 4.8 s after applying the electric current is depicted in figure 3a, where the three zones (TE, sample and LE) can be clearly distinguished. It can also be seen that the model handles properly the abrupt transition between the different zones without oscillations in the case of the electrolyte concentrations and conductivity. In contrast, the pH shows small oscillatory transitions between the zones meaning that diffusion current in the transport equations is correctly included in the formulation, but it is not consistent with the stabilization terms. We can conclude this, by carefully taking into account all the considerations present in [44], as far as the
$ . / i n i c o n d . sh $ . / electroMicroTransport
5. Validation examples The three tutorials provided within the electroMicroTransport toolbox serve also as validation examples, since they constitute well-known cases previously used as benchmarks for the evaluation of other 1D models [23]. 5.1. Capillary Zone Electrophoresis In CZE charged chemical species are separated based on differences in their electrophoretic mobilities, in a uniform background electrolyte (BGE) placed in a capillary tube [1]. In this case, the BGE is a buffer composed of acetic acid 20 mM, and Tris 12 mM in a 20 cm capillary, where a 2500 A m−2 electrical current density is applied. A 5 mm length sample consisting of 1 mM pyridine and 1 mM aniline is placed 5 mm away from the anode position as seen in Fig. 2a. At the boundaries, concentration fields are managed with the inletOutlet condition. Regarding the grid, the capillary tube was divided into 16000 linear segments. The simulation was run for 240 s giving the concentration profiles for pyridine and aniline shown in Fig. 2c. 8
Figure 2: Validation results for tutorial example 5.1. a) Initial aniline and pyridine concentration profiles and separation achieved by CZE simulation after b)120 s and c) 240 s. Results obtained with GENTRANS are represented with dotted lines [23].
6. Application example: Double-L injection device
oscillatory behavior for the pH vanishes for finer grids (see supplemetary material). Precisely, it was demonstrated that this oscillatory shape is related to numerical conditions, while a hollow-like shape for the pH corresponds to those reported cases where diffusion current was neglected in the modeling.
To show the ability of the electroMicroTransport toolbox to solve more complex cases, the full operation of a microfluidic Double-L device, for pinched injection and CZE, as shown in Fig. 5, was simulated [46]. Proper representation of the problem requires EOF, electric potential and transport to be handled in 2D. A double-L injection device allows controlling plug size and sample leakage, improving detection efficiency. Different modes of operation are possible according to potentials applied to the four inlets/reservoirs [47]. The operation mode simulated involves two sequential steps known as injection and a separation. In the first step, a sample composed of glycine (Gly), alanine (Ala) and phenylalanine (Phe) 1 mM each is inserted via the inlet 2 in a 2.3 M acetic acid buffer. In the same inlet a 100 V potential is applied, while Outlet 2 is grounded. The sample is transported due to the EOF and electromigration until reaching the separation channel. In the next step, a 1000 V separation potential is applied at the inlet 1 while outlet 1 is grounded. Thus, the sample plug migrates through the main channel where electrophoretic separation occurs. Boundary conditions used are summarized in table 2. A total of 181250 cells was used, where the width of the channel was discretized with 50 elements and the grid was refined at the inlets, outlets and channel intersections. Simulation results are shown in Fig. 5. In particular, Fig. 5e shows the time evolution of the sample (c sample = cAla + cGly + cPhe ) in a virtual detector located at 3.0 cm from the injection arm. This signal can be correlated with the electropherogram obtained with such detector. Here, it can be seen that the peaks appear in the correct order: i) Glycine
5.3. Isoelectric focusing Amphoteric species are chemical compounds which can act as an acid or a base according to the pH of the electrolyte solution. The pH at which the molecule exhibits no charge is called the isoelectric point (pI). IEF consists on separating amphoteric species based on their pI, through the use of a pH gradient usually generated with a mixture of carrier ampholytes. Thus, IEF experiments require a large number of components [45]. The IEF problem solved in this paper consists of a 5 cm capillary where 35 hypothetical carrier ampholytes are uniformly distributed at 1 mM concentration each. The carrier ampholytes have pI values ranging from 3.1 to 9.9 with ∆pI = 0.2, pK+1 = pI −1 and pK−1 = pI +1. A constant 300 V cm−1 electric field is applied to the boundaries, which are closed, avoiding any kind of mass flow. The grid used for simulation consisted of 10000 linear segments. Fig. 4 shows the concentration profiles, the pH and conductivity along the capillary after 10 min. It can be seen how the carrier ampholytes order themselves with their pI increasing from anode to cathode. Also, it can be appreciated how the pH increases monotonically along the capillary. Results obtained are very similar to those obtained with other solver methods [23]. 9
Figure 3: Validation results for tutorial example 5.2. a) Cations concentration b) conductivity and pH profiles along the capillary tube, after a 4.8 s ITP simulation of the system described in 5.2. Results obtained with GENTRANS are plotted with dotted lines [23].
(tGly = 121.6 s), ii) Alanine (tAla = 127.9 s)iii) Phenylalanine (tPhe = 156.8 s). Moreover, the position of the peaks coincides with an error lower than 4% with the position calculated analytically [48], taking into consideration the position of the detector, the applied field, the length of the main channel, the effective charge of each species and the electroosmotic flow(110 µm s−1 ) at the working pH (pH = 2.2). The simulation performed allows, among other things, the assessment of the peak longitudinal dispersion, and the magnitude of the sample leakage produced by the injection method. It is worth to mention here that the 2D simulation enables the user to see the tailing effect for the peaks. This effect, that is typical from these injection procedures made with microfluicic chips [49], is clearly captured by the simulation, emphasizing the difference between the presented toolbox and the other 1D tools previously developed for electrophoretic separations.
i) a boundary condition implementing the HelmoltzSmoluchowsky equation Eq. (6) for EOF modeling via the thin EDL approximation, ii) a boundary condition for Φ for fixing the current density i on single column problems, iii) the use of the native inletOutlet boundary condition accounting not only for the fluid flow flux φh , but also for the electromigrative flux φe , iv) a dual time-stepping algorithm to deal with the different hydrodynamic and electromigrative time-scales, and v) Courant number calculation involving the effective φ and hydrodynamic φh fluxes for every cell and species. Moreover, three examples are given for both tutorial and validation purposes. Furthermore, a final tutorial example of a complex 2D real problem is included to demonstrate the toolbox capabilities. Finally, the flexibility of the presented tool endows it with the possibility to easily expand it to new applications, like, for example, paper electrophoresis.
7. Concluding remarks
Acknowledgments We want to acknowledge Wolfgang Thormann for his valuable contribution to this work. This research was supported by CONICET, ANPCyT (Grant PICT 20160640), and UTN (Grant PID ASUTNFE0004475), Argentina.
A toolbox for the solution of 3D multiphysic problems involving fluid dynamic and electromigrative phenomena in capillaries and microfluidic chips was developed. This toolbox was implemented by using OpenFOAM R meaning that native 3D problem handling, parallel computing and GNU-GPL license are features supported from scratch. Also, the electroMicroTransport toolbox implements a number of novel characteristics that makes it particularly interesting for the scientific community, namely:
Supplementary material Supplementary material associated to different numerical tests performed with the toolbox and a technical 10
Figure 4: Validation results for tutorial example 5.3. a) Ampholyte concentration, b) conductivity and pH profiles along the capillary tube, achieved with a 10 min simulation of a IEF separation. The gray lines are the pH and conductivity counterparts obtained with GENTRANS [23].
discussion about numerical convergence can be found in the on-line version of the article.
[12] Q. Mao, J. Pawliszyn, W. Thormann, Dynamics of capillary isoelectric focusing in the absence of fluid flow: High-resolution computer simulation and experimental validation with whole column optical imaging, Anal. Chem. 72 (21) (2000) 5493– 5502. [13] W. Thormann, J. Caslavska, R. Mosher, Modeling of electroosmotic and electrophoretic mobilization in capillary and microchip isoelectric focusing, J. Chromatogr. A 1155 (2) (2007) 154–163. [14] R. A. Mosher, W. Thormann, High-resolution computer simulation of the dynamics of isoelectric focusing: In quest of more realistic input parameters for carrier ampholytes, Electrophoresis 29 (5) (2008) 1036–1047. [15] F. Pimenta, M. A. Alves, Numerical simulation of electricallydriven flows using openfoam, arXiv preprint arXiv:1802.02843. [16] W.-L. Hsu, D. W. Inglis, M. A. Startsev, E. M. Goldys, M. R. Davidson, D. J. Harvie, Isoelectric focusing in a silica nanofluidic channel: effects of electromigration and electroosmosis, Anal. Chem. 86 (17) (2014) 8711–8718. [17] A. Chatterjee, Generalized numerical formulations for multiphysics microfluidics-type applications, J. Micromech. Microeng. 13 (5) (2003) 758. [18] J. Shim, P. Dutta, C. F. Ivory, Modeling and simulation of ief in 2-d microgeometries, Electrophoresis 28 (4) (2007) 572–586. [19] K. Yoo, J. Shim, J. Liu, P. Dutta, Mathematical and numerical model to study two-dimensional free flow isoelectric focusing, Biomicrofluidics 8 (3) (2014) 034111. [20] K. Yoo, J. Shim, J. Liu, P. Dutta, Efficient algorithm for simulation of isoelectric focusing, Electrophoresis 35 (5) (2014) 638– 645. [21] P. Kler, C. Berli, F. Guarnieri, Modeling and high performance simulation of electrophoretic techniques in microfluidic chips, Microfluid. Nanofluid. 10 (1) (2011) 187. [22] P. Kler, L. Dalcin, R. Paz, T. Tezduyar, Supg and discontinuitycapturing methods for coupled fluid mechanics and electrochemical transport problems, Comput. Mech. 51 (2) (2013) 171. [23] S. Mikkonen, H. Ekstr¨om, W. Thormann, High-resolution dynamic computer simulation of electrophoresis using a multi-
References [1] J. P. Landers, Handbook of capillary and microchip electrophoresis and associated microtechniques, CRC press, 2007. [2] P. A. Kler, D. Sydes, C. Huhn, Column–coupling strategies for multidimensional electrophoretic separation techniques, Anal. Bioanal. Chem. 407 (1) (2015) 119–138. [3] M. Ramos-Pay´an, J. A. Oca˜na-Gonzalez, R. M. Fern´andez´ Bello-L´opez, Recent trends in capTorres, A. Llobera, M. Angel illary electrophoresis for complex samples analysis: A review, Electrophoresis (2017) 111–125. [4] W. Thormann, J. Caslavska, M. C. Breadmore, R. A. Mosher, Dynamic computer simulations of electrophoresis: Three decades of active research, Electrophoresis 30 (2009) S16–S26. [5] M. Bier, O. Palusinski, R. Mosher, D. Saville, Electrophoresis: Mathematical modeling and computer simulation, Science 219 (4590) (1983) 1281–1287. [6] W. Thormann, J. Caslavska, R. Mosher, Modeling of electroosmotic and electrophoretic mobilization in capillary and microchip isoelectric focusing, J. Chromatogr. A 1155 (2) (2007) 154–163. [7] V. Hruˇska, M. Jaroˇs, B. Gaˇs, Simul 5–free dynamic simulator of electrophoresis, Electrophoresis 27 (5-6) (2006) 984–991. [8] M. Bercovici, S. K. Lele, J. G. Santiago, Open source simulation tool for electrophoretic stacking, focusing, and separation, J. Chromatogr. A 1216 (6) (2009) 1008–1018. [9] R. A. Mosher, M. C. Breadmore, W. Thormann, High-resolution electrophoretic simulations: Performance characteristics of onedimensional simulators, Electrophoresis 32 (5) (2011) 532–541. [10] B. Y. Moghadam, K. T. Connelly, J. D. Posner, Isotachophoretic preconcenetration on paper-based microfluidic devices, Anal. Chem. 86 (12) (2014) 5829–5837. [11] X. Li, L. Luo, R. M. Crooks, Low-voltage paper isotachophoresis device for dna focusing, Lab Chip 15 (20) (2015) 4090–4098.
11
Figure 5: Simulation of a Double-L injection device for electrophoretic separations. Dimensions and name of the boundaries are indicated in the main figure, while sample profile (c sample = cAla + cGly + cPhe ) at different instants is shown in inset a)–d) respectively. a) Sample entering the injection arm when t =22 s. b) Sample already arrived channel intersections and is headed to Outlet 2, 60 s after the experience began. c) and d) separation obtained in the main arm, when t =120 s, during the separation step. e) Time evolution of c sample . The first peak is due to Glycine (zΩ =−21.83 × 10−9 m2 V−1 s−1 ), the second due to Alanine (zΩ =−18.39 × 10−9 m2 V−1 s−1 ) and the third to Phenylalanine zΩ =−12.02 × 10−9 m2 V−1 s−1 ).
[24] [25] [26] [27] [28]
[29] [30]
physics software platform, J. Chromatogr. A 1532 (2018) 216– 222. H. G. Weller, G. Tabor, H. Jasak, C. Fureby, A tensorial approach to computational continuum mechanics using objectoriented techniques, Comput. Phys. 12(6) (1998) 620–631. OpenFOAM the open source CFD tool, https://www. openfoam.com/, accessed: 2018-03-19. J. Lyklema, Fundamentals of interface and colloid science: soft colloids, Vol. 5, Academic press, 2005. R. Probstein, Physicochemical Hydrodynamics. An Introduction, 2nd Edition, Wiley-Interscience, 2003. P. A. Kler, L. D. Dalcin, R. R. Paz, T. E. Tezduyar, Supg and discontinuity-capturing methods for coupled fluid mechanics and electrochemical transport problems, Computational Mechanics 51 (2) (2013) 171–185. J. M. MacInnes, Computation of reacting electrokinetic flow in microchannel geometries, Chem. Eng. Sci. 57 (21) (2002) 4539–4558. T. L. Sounart, J. C. Baygents, Lubrication theory for electroosmotic flow in a non-uniform electrolyte, J. Fluid Mech. 576 (2007) 139–172.
[31] T. J. Craven, J. M. Rees, W. B. Zimmerman, On slip velocity boundary conditions for electroosmotic flow near sharp corners, Phys. Fluid. 20 (4) (2008) 043603. [32] E. Brunet, A. Adjari, Generalized Onsager relations for electrokinetic effects in anisotropic and heterogeneous geometries, Phys. Rev. E 69 (1) (2004) 016306. [33] C. L. A. Berli, Equivalent circuit modeling of electrokinetically driven analytical microsystems, Microfluid. Nanofluid. 4 (5) (2008) 391–399. [34] C. L. A. Berli, M. Piaggio, J. Deiber, Modeling the zeta potential of silica capillaries in relation to the background electrolyte composition, Electrophoresis 24 (2003) 1587–1595. [35] I. Arnaud, J. Josserand, J. Rossier, H. Girault, Finite element simulation of off-gel buffering, Electrophoresis 23 (2002) 3253– 3261. [36] J. Ferziger, M. Peric, Computational Methods for Fluid Dynamics, Vol. I, Springer-Verlag, 2002. [37] H. Jasak, Error analysis and estimation for the finite volume method with applications to fluid flows, Ph.D. thesis, Department of Mechanical Engineering Imperial College of Science, Technology and Medicine (1996).
12
[38] S. M´arquez Dami´an, An extended mixture model for the simultaneous treatment of short and long scale interfaces, Ph.D. thesis, FICH, Universidad Nacional del Litoral, Santa Fe, Argentina (2013). [39] M. Darwish, F. Moukalled, TVD schemes for unstructured grids, Int. J. Heat Mass Tran. 46 (4) (2003) 599–611. [40] P. L. Roe, Characteristic-based schemes for the Euler equations, Annual Review of Fluid Mechanics 18 (1) (1986) 337–365. [41] R. Issa, Solution of implicitly discretised fluid flow equations by operator splitting, Journal of Computational Physics 62 (1986) 40–65. [42] electoMicroTransport repository, https://gitlab.com/ santiagomarquezd/electroMicroTransport/, accessed: 2018-03-19. [43] P. A. Kler, C. Huhn, Non-aqueous electrolytes for isotachophoresis of weak bases and its application to the comprehensive preconcentration of the 20 proteinogenic amino acids in column-coupling itp/ce–ms, Anal. Bioanal. Chem. 406 (28) (2014) 7163–7174. [44] T. Sounart, J. Baygents, Simulation of electrophoretic separations: effect of numerical and molecular diffusion on ph calculations in poorly buffered systems, Electrophoresis 21 (12) (2000) 2095–2103. [45] R. A. Mosher, W. Thormann, High-resolution computer simulation of the dynamics of isoelectric focusing: In quest of more realistic input parameters for carrier ampholytes, Electrophoresis 29 (5) (2008) 1036–1047. [46] J. M. Karlinsey, Sample introduction techniques for microchip electrophoresis: a review, Anal. Chim. Acta 725 (2012) 1–13. [47] C.-H. Lin, R.-J. Yang, C.-H. Tai, C.-Y. Lee, L.-M. Fu, Double-l injection technique for high performance capillary electrophoresis detection in microfluidic chips, J. Micromech. Microeng. 14 (4) (2004) 639. [48] P. Coufal, J. Zuska, T. van de Goor, V. Smith, B. Gaˇs, Separation of twenty underivatized essential amino acids by capillary zone electrophoresis with contactless conductivity detection, Electrophoresis 24 (4) (2003) 671–677. [49] D. Sydes, P. A. Kler, M. Hermans, C. Huhn, Zero-deadvolume interfaces for two–dimensional electrophoretic separations, Electrophoresis 37 (22) (2016) 3020–3024.
Table 2: Set of boundary conditions used to simulate the double-L injection device presented in section 6.
Boundary
Condition type
Equation: Transport Walls zeroGradient inletOutlet Inlet 1 Inlet 2
inletOutlet
Outlet 1
inletOutlet
Outlet 2
inletOutlet
Value
–c j =0 mM –cAcet =2300 mM –cAla =0.01 mM –cGly =0.01 mM –cPhe =0.01 mM –cAcet =2300 mM –c j =0 mM –cAcet =2300 mM –c j =0 mM –cAcet =2300 mM
Equation: Charge conservation (injection step) Walls zeroGradient Inlet 1 zeroGradient Inlet 2 fixedValue –Φ =100 V Outlet 1 zeroGradient Outlet 2 fixedValue –Φ =0 V Equation: Charge conservation (separation step) Walls zeroGradient Inlet 1 fixedValue –Φ =1000 V Inlet 2 zeroGradient fixedValue –Φ =0 V Outlet 1 zeroGradient Outlet 2 zeroGradient Equation: Fluid-Flow Walls SmoluchowskyWallVelocity zeroGradient –p Inlet 1 zeroGradient –u fixedValue –p =0 Pa zeroGradient –u Inlet 2 fixedValue –p =0 Pa Outlet 1 zeroGradient –u fixedValue –p =0 Pa Outlet 2 zeroGradient –u fixedValue –p =0 Pa
13