Adaptive mesh simulations of ionic systems in microcapillaries based on the estimation of transport times

Computers and Chemical Engineering 30 (2006) 674–685

Adaptive mesh simulations of ionic systems in microcapillaries based on the estimation of transport times ˇ Michal Pˇribyl ∗ , Dalimil Snita, Milan Kub´ıcˇ ek Department of Chemical Engineering and Centre for Nonlinear Dynamics of Chemical and Biological Systems, Institute of Chemical Technology, Prague, Technick´a 5, 166 28 Prague 6, Czech Republic Received 29 April 2005; received in revised form 10 November 2005; accepted 23 November 2005 Available online 19 January 2006

Abstract We present results of testing an empirical algorithm of mesh adaptation for modeling of spatially one-dimensional reaction-transport systems with initially separated components and large moving gradients. The algorithm combines newly developed procedures of r-refinement and standard FEMLAB procedures. The developed method uses an expansion of a dense mesh in neighborhoods of localized large gradients and estimates the mesh adaptation interval (based on the evaluation of transport times). The size of the neighborhood is dynamically controlled. The adaptation method has been tested on systems of parabolic–elliptic PDEs with extremely large moving gradients of concentrations and electric potential. We have found that CPU time requirements are comparable to other mesh adaptation solvers. The simplicity of the empirical procedure of the mesh adaptation and its easy implementation in standard dynamic solvers represent main advantages of the proposed method. The studied system, called an “electrolyte diode”, is described by four mass balances of the electrolyte components (parabolic PDEs) and by Poisson equation of electrostatics (an elliptic PDE). Dynamics of formation of open and closed modes of the electrolyte diode are described and current–voltage characteristics are explained. The results obtained with the developed solver agree with the results of steady state analysis. Our tests prove that the proposed algorithm of mesh adaptation can be used in modeling of microfluidics application, e.g., DNA chips, capillary electrophoresis and isoelectrical focusing, where formation of extremely large gradients of electric potential and concentrations of the electrolyte components is expected. © 2005 Elsevier Ltd. All rights reserved. Keywords: Microchip; Modeling; Mesh adaptation; Moving mesh; Dynamical simulation; Electrolyte diode

1. Introduction This work is focused on mathematical modeling of reacting ionic systems in microcapillaries with initially separated electrolytes. A microcapillary connecting two reservoirs with well-defined concentration of chemical components (salts, buffers, acids, bases, etc.) is considered (see Fig. 1A). A difference of electric potential can be imposed on the microcapillary boundaries. Similar microcapillary systems are often used, for example, in diagnostic microchips, capillary electrophoresis, isoelectrical focusing, fast PCR, etc. (Fukuba, Yamamoto, Naganuma, & Fujii, 2004; Han & Singh, 2004; Mitnik et al., ˇ 2001; Pˇribyl, Snita, Hasal, & Marek, 2004; Righetti, Gelfi, &

∗

Corresponding author. Tel.: +420 2 2044 3168; fax: +420 2 2044 4320. E-mail address: [email protected] (M. Pˇribyl).

0098-1354/$ – see front matter © 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.compchemeng.2005.11.006

D’Acunto, 2002; Schmalzing, Koutny, Taylor, Nashabeh, & Fusch, 1997). The microcapillary contains pure water at the start. Then ionic components from the reservoirs are transported into the microcapillary by diffusion and/or electromigration mechanisms. The transport is usually accompanied, especially at the start, with large concentration gradients and with the gradient of electric potential. These gradients move to the center of the microcapillary (see Fig. 1B). Some of the ionic chemical components from the opposite compartments can interact within the microcapillary. These reaction zones are characterized by specific properties. For example, the interaction between protons and hydroxyl ions gives rise to a region with low electric conductivity and a sharp change of electric potential. (Lindner, Snita, & Marek, 2002; Merkin, Simon, & Noszticzius, 2000). If an insoluble chemical component is a result of an ionic interaction then regions of supersaturation could be

M. Pˇribyl et al. / Computers and Chemical Engineering 30 (2006) 674–685

Fig. 1. (A) Scheme of the electrolyte microsystem with a connecting microcapillary and the imposed electric field. (B) The gradients of concentrations and/or electric potential move through the microcapillary in time.

ˇ detected inside the microcapillary (Pˇribyl, Snita, & Marek, 2005). Dynamical simulation of transients in such microsystems is often difficult because of extremely large moving gradients of the model variables and the resulting stiffness of the problems. Electrolyte microsystems in one spatial dimension are usually described by a set of parabolic partial differential equations (PDE, in this case mass balance equations) and one elliptic partial differential equation (equation of electrostatics). A system of differential–algebraic equations (DAE) is obtained after the spatial discretization of the domain. There are methods that use a uniform fixed mesh (Lim, Le Lann, & Joulia, 2001) or a piecewise uniform mesh (Liu & Bhatia, 1999, 2001) for solution of problems with moving gradients. These methods are suitable for singularly perturbed problems, e.g., adsorption, with large gradients in thin layers. However a non-equidistant mesh has to be often generated to meet the accuracy requirements. If large gradients move outside the regions where the mesh is dense, numerical problems can arise and the approximation errors grow. Hence, a dynamic adaptation of the spatial mesh is often necessary. In principle, there are three kinds of methods of the dynamic adaptation (Bieniasz, 2000): (i) r-refinement or the moving mesh method (a fixed number of the grid points, the point positions move), (ii) h-refinement or the local refinement method (a number of grid points is not fixed, the original points do not move), and (iii) p-refinement (a fixed mesh is used, the order of spatial discretization changes). A combination of these methods can be used. Almost all methods of dynamic mesh adaptation are based on the equidistribution principle according to which the mesh point localization is selected so that a measure of the solution error is equalized over each subinterval on the mesh (Huang, Ren, & Russell, 1994). The error distribution is monitored by a monitor function, which can be, for example, first or second derivative of a dependent variable or their combinations. The first derivative localizes the large gradients in the system and the second derivative localizes sudden folds on spatial profiles. The second derivative control was successfully used by Kulkarni and Dudukovic (1996) for reactors with moving temperature fronts.

675

A class of the r-methods that use curvilinear coordinates instead of the Cartesian coordinates was described and reviewed, for example, by Adjerid and Flaherty (1986) and Thomson (1985). Computation is then done on a fixed grid in the curvilinear space even though the grid points move in the Cartesian space. There is no need for interpolation of the solution on the new mesh. Coimbra, Sereno, and Rodrigues (2004) recently proposed an r-method suitable for reaction–diffusion–convection systems. New positions of the mesh points on spatially 1D and 2D domains is predicted by minimizing the norm of the residual of the approximation in the model PDEs. The method requires introducing of empirical constants for a particular model. Another adaptive moving mesh method based on this idea was reported by Huang, Zheng, and Zhan (2002). The method is based on an alternating procedure where a new mesh distribution is computed from artificial parabolic moving mesh PDEs and then the physical PDEs are integrated on the new mesh. The necessity of computation of the mesh distribution from the additional moving mesh PDEs is the main disadvantage of the method. An h-method for solving parabolic PDEs based on finite volume methods was presented in the work of Roussel, Schneider, Tsigulin, and Bockhorn (2003). They successfully simulate a diffusion–convection PDE with extremely large gradients of velocity. The size ratio of the maximal and the minimal volume up to 211 was used. Comparison of available software codes that employs mesh adaptation for solution of parabolic PDEs was made by Wang, Keast, and Muir (2004a, 2004b). Their software BACOL was most efficient for problems with rapidly moving spatial gradients. BACOL uses an adaptive collocation method and allows a decrease or an increase of the number of subintervals. The spatial error is evaluated after each integration step, which can decelerate the simulation procedure. The software EPDCOL also uses the collocation method (Keast & Muir, 1991). The NAG library source codes contain the algorithm D03PPF (D03PPF), which is an h-refinement method based on a finite difference method. TOMS731 (Blom & Zegeling, 1994) and MOVCOL (Huang & Russell, 1996) are r-refinement codes. They use the finite element method and the cubic Hermite collocation. Bieniasz (2000) has dealt with electrochemistry problems where large gradients of spatially and time dependent variables are presented. He developed an h-refinement method based on finite differences called the patch adaptive strategy (PAS). Each integration step has to be evaluated at least for two different densities of the grid points. A new mesh refinement and reevaluation of the solution has to be carried out in the places with an unacceptable error. This probably substantially slows down the dynamic simulation. However, the PAS method was successfully used in simulating the Nernst–Planck–Poisson equations (Bieniasz, 2004a, 2004b), the Nernst–Planck-electroneutrality equations (Bieniasz, 2004a, 2004b), and other PDEs producing large moving gradients (Drake & Manoranjan, 1996). A combined hp-adaptive strategy was reported by Almeida and Silva (2002). The method allows using a high order discretization method in spatial regions with sharp gradients. Simultaneous computation of the mesh parameters (order

676

M. Pˇribyl et al. / Computers and Chemical Engineering 30 (2006) 674–685

and size of each element) is based on minimizing of error functions. Currently, the r-refinement methods are less developed than the h-refinement ones. On the other hand, the r-refinement method can be easily implemented in any available fixed mesh software and difficulties arising from change in the number of grid points can be avoided (Cao, Huang, & Russell, 2003). All the above-mentioned approaches for solution of the problems with moving gradients required development of spatial and time discretizations of model equations and a complete implementation. In this work, a different strategy has been chosen. The standard FEMLAB routines have been used (based on finite element methods) for time integration of the spatially discretized system of model equations (FEMLAB 2.3 ref. manual). After evaluation of the monitor function at a selected time, new distribution of the mesh is determined. Evaluation of the monitor function and the successive mesh adaptation at each time integration step substantially decelerates the numerical simulations. For that reason, the expansion of the dense mesh in neighborhoods of the large gradients has been applied, which enables to make the mesh adaptations after longer time intervals. The time interval between two mesh adaptations is estimated as the transport time of the chemical components through the region of the dense mesh. Thiele, Falk, and Tobiska (2001) selected similar strategy for prediction of elution profiles in the annular chromatography. However, their problem is described by a set of stationary PDEs. Thus the spatial mesh is fixed but regions with a higher mesh density are a priori determined from the knowledge of transport parameters of the system. In our case, the large gradients move in time and space and it is not possible to predict their spatial position on the entire time domain. In dynamical systems, the transport times (the diffusion time and the electromigration time) of chemical components at large gradients have to be evaluated repeatedly after a time period. This time period can be estimated from the size of the region with the dense mesh and from the transport velocities of ionic components. General features of the algorithm are described in the following section. The algorithm has been tested on two microcapillary systems. A microsystem with precipitation interaction is the first one. The obtained results were recently reported (Pˇribyl et al., 2005). In this paper, dynamical behavior of the system with initially separated components (a strong acid and a strong base) is analyzed. We call this system an “electrolyte diode”. Typical imposed differences per unit length of the pH value and the electric potential were pH 12 mm−1 and
Φ = 10 V mm−1 , respectively. This gives rise to extremely large gradients on the spatial profiles of pH and electric potential (locally ∂Φ/∂x ≈ 106 V m−1 ). Among the above-cited methods, only one was tested on a system with such extreme gradients (Roussel et al., 2003). Capabilities and properties of the proposed algorithm are tested on a transient behavior study of the electrolyte diode system. The proposed dynamical solver helps to determine stability of the steady state solutions that were recently obtained by the use of a steady state solver (Lindner et al., 2002).

2. Reaction-transport model of a microcapillary The mathematical model describes processes in an effectively one-dimensional (capillary) microsystem shown in Fig. 1A. The capillary connects two reservoirs. Compositions of the solutions and the electric potential value in the reservoirs are considered to be constant in time. Transport of components and chemical interactions occur in the microcapillary. If no convection flow is present in the microcapillary, e.g., in a gel-like medium, the distribution of concentration of the i-th component in the microcapillary connecting both reservoirs is described by the mass conservation equation:
 ∂ ∂ci zi Di F ∂Φ ∂ci =− −Di − ci ∂t ∂x ∂x RT ∂x +

n 

υij rj ,

i = 1, . . . , m

(1)

j=1

where the total flux of the i-th component is given by Nernst–Planck equation. Di , zi , and ci denote diffusion coefficient, charge number, and concentration of the i-th component, respectively. The symbols F, R, T, and Φ represent Faraday constant, the molar gas constant, temperature, and electric potential, respectively. Components of the electrolyte can interact in n chemical reactions with reaction rates rj . The distribution of electric potential is given by the solution of Poisson equation of electrostatics on the spatially onedimensional domain: ∂2 Φ F  =− z i ci 2 ∂x εr ε0

(2)

i

where the product εr ε0 denotes the permittivity of the electrolyte. The composition of the reservoir solutions and the imposed difference of electric potential define boundary conditions of the system (see Fig. 1A). The Dirichlet boundary conditions with constant values of concentrations and a constant value of the electric potential are considered ci |x=0 = ci,L , ci |x=l = ci,R , Φ|x=0 = ΦL , Φ|x=l = ΦR .

i = 1, . . . , m

(3)

3. Case study Chemical interactions between the solutions of potassium hydroxide (the L reservoir) and hydrochloric acid (the R reservoir) are subject of the case study. The neutralization between the hydrogen and the hydroxyl ions occurs in the microcapillary H+ + OH−
H2 O.

(4)

Kinetics of this chemical reaction can be expressed as: rw = kw (cH+ cOH− − Kw ),

(5)

where kw and Kw are the kinetic water recombination constant and the equilibrium constant, respectively.

M. Pˇribyl et al. / Computers and Chemical Engineering 30 (2006) 674–685

677

The electrolyte solutions in the reservoirs L and R are considered to be in chemical equilibrium and satisfy the condition of local electroneutrality: Kw = cH+ cOH− ,  zi ci = 0.

(6) (7)

i

The boundary concentrations were chosen to be cCl− ,L = 0 kmol m−3 , cK+ ,L = 0.1 kmol m−3 , cCl− ,R = 0.1 kmol m−3 , cK+ ,R = 0 kmol m−3 . The concentrations of the hydroxyl and hydrogen ions on the boundaries are computed from Eqs. (6) and (7). The value of electric potential on the L boundary was kept at ΦL = 0 V. The value ΦR is the studied parameter of the model and its value is specified in the figure captions. We have considered all processes in the microcapillary to occur far from the electrodes. Hence, the interactions at the electrode surface are not included in the model. The obtained mathematical model consists of the four parabolic PDEs (1) for the component: H+ , OH− , K+ , and Cl− , and one elliptic PDE (2) for the evaluation of electric potential. One chemical interaction is considered, Eq. (5). The boundary conditions are expressed by Eq. (3). In dynamical simulations, we have considered a microcapillary filled with pure water at pH 7 as an initial condition of the microsystem. 4. Numerical algorithm The developed numerical algorithm is capable to solve the set of Eqs. (1) and (2) in a spatially one-dimensional system with moving sharp concentration gradients and/or gradients of electric potential. The algorithm combines the standard FEMLAB procedures and the developed empirical procedure for dynamical adaptation of a non-equidistant mesh based on the estimation of transport times. The FEMLAB procedures employ the finite element method for the space discretization. The Lagrange quadratic elements were used in our study. If we fix the number of discrete elements on the studied domain (nx), the shortest elements should be concentrated in the regions of the sharp gradients in order to minimize the error of numerical approximation of the solution (the size of elements is given by a vector h of dimension nx). The problem is that the gradients can move in time. Hence, it is convenient to build a mesh with a high density in a certain neighborhood of the sharp gradient regions to cover a gradient shift in a time interval. The radius of the neighborhood (d) can be initially chosen and updated during the simulation. The time interval
t between two adaptations of the spatial mesh (i.e. adaptation interval) can be then estimated from the knowledge of the rate of mass transport within the sharp spatial gradients. 4.1. The main procedure (Fig. 2) Optional parameters have to be chosen before numerical simulations: the number of finite elements (nx), minimal radius of the neighborhood (dmin ), maximal radius of the neighborhood (dmax ), final time of the dynamical simulation (tend ), maximal

Fig. 2. The main algorithm of the solver.

time interval between two adaptations
tmax and the parameter tstop (Fig. 2). We choose d = dmin for the first adaptation interval. If the sharp gradients vanish during integration, the estimated
t could be higher than
tmax and thus we chose
t =
tmax in such case. Further, if the system is reaching a steady state, the gradients of concentrations and/or electric potential stop moving. The time tstop at which the system is near the steady state can be estimated from the knowledge of the transport times of the chemical components in the system. If the time of dynamical simulation t > tstop then the adaptation time in the next rounds is
t =
tmax and the density of the mesh in any neighborhood of the sharp gradients need not be increased. Hence a good guess of tstop can substantially speed up dynamical simulation. We used dmin = 2 ␮m, dmax = 100 ␮m and nx = 1000 in the case study. Consistent initial conditions and initial mesh have to be selected before starting the integration procedure. The dynamical simulation runs in a cycle, in which the time integration of

678

M. Pˇribyl et al. / Computers and Chemical Engineering 30 (2006) 674–685

Fig. 3. (A) The detail of the initial profile of pH at the L boundary. (B) The detail of the initial profile of pH at the R boundary. (C) The initial profile of the electric potential; ΦR = 10 V.

the model equations is carried out on time interval
t between two mesh adaptations, the radius of neighborhood is estimated (if t < tstop ), the new mesh is created and new time interval between two mesh adaptations is determined and integration is again carried out. The standard FEMLAB solver daspk (FEMLAB 2.3 ref. manual) based on the backward differentiation formulae is used for the time integration of the discretized set of the model equations. The values of relative and absolute tolerances for all quantities have to be carefully selected (we have chosen: Rtol = 10−8 , AtolcH = AtolcOH = 10−10 mol m−3 , AtolcK = AtolcCl = 10−3 mol m−3 , Atol = 10−5 V). 4.2. The problem of consistent initial conditions Although the boundary concentration of chemical components can be chosen arbitrarily, the discontinuity in the concentration profiles at the boundaries of the domain bring problems with evaluation of the initial profile of electric potential. Here, we have used the following strategy: concentrations of chloride and potassium ions are frozen at zero value inside the microcapillary (pure water in the microcapillary) and on the microcapillary boundaries. In the next step, the molar balances of hydrogen and hydroxyl ions, Eq. (1), together with the Poisson equation, Eq. (2), are solved by the Newton method to obtain a steady state distribution of these quantities on an equidistant mesh (usually nx = 2000). Because the electroneutrality condition, Eq. (7), is not satisfied on the boundaries, extremely large gradients of the concentrations of hydrogen and hydroxyl ions and electric potential are observed (Fig. 3). Similar phenomena can be found at surfaces with bounded electric charge (Probstein, 1994). The dimension of this boundary layer corresponds to the dimension of the Debye length observed for electrolytes with the characteristic concentration of ions about ≈10−7 kmol m−3 (Dutta & Beskok, 2001). The pH value of 7 is established in a short distance (few microns) from the microcapillary boundaries, which reflects the presence of pure water in the microcapillary at the start of simulation. The spatial mesh is then reconstructed according to values of the computed monitor functions (see the next paragraph) and the steady state solution is again computed on the refined mesh with nx = 1000. Finally, the potassium and chloride boundary conditions are set to the correct values (see Table 1). The computed steady state

profiles of cH+ , cOH− , and Φ together with zero concentration profiles of potassium and chloride ions (with discontinuity at the boundaries) are successfully used as an initial condition for the dynamical simulation. 4.3. Mesh adaptation (Fig. 4) The value of a reciprocal monitor function eri is computed by means of the chosen monitor functions f j that are evaluated in j the center of each finite element i, i.e. fi , after integration over Table 1 Model parameters Parameter (unit)

Description

Value

DH+ (m2 s−1 ) DOH− (m2 s−1 )

Diffusion coefficient of protons Diffusion coefficient of hydroxyl ions Diffusion coefficient of chloride ions Diffusion coefficient of potassium ions pH value on the R boundary Concentration of chloride ions on the R boundary Concentration of potassium ions on the R boundary pH value on the L boundary Concentration of chloride ions on the L boundary Concentration of potassium ions on the L boundary Electric potential on the L boundary Electric potential on the R boundary Dimension of the system Charge number of protons Charge number of hydroxyl ions Charge number of chloride ions Charge number of potassium ions Faraday constant Molar gas constant Temperature Vacuum permittivity Relative permittivity of water Kinetic constant of water recombination Ionic product of water

9.31 × 10−9 5.28 × 10−9

DCl− (m2 s−1 ) DK+ (m2 s−1 ) pHR cCl− ,R (kmol m−3 ) cK+ ,R (kmol m−3 ) pHL cCl− ,L (kmol m−3 ) cK+ ,L (kmol m−3 ) ΦL (V) ΦR (V) l (m) zH+ zOH− zCl− zK+ F (C mol−1 ) R (J mol−1 K−1 ) T (K) ε0 (F m−1 ) εr kw (m3 mol−1 s−1 ) Kw (mol2 m−6 )

2.04 × 10−9 1.96 × 10−9 1 0.1 0 13 0 0.1 0

1 × 10−3 1 −1 −1 1 9.6487 × 104 8.314 310 8.8542 × 10−12 7.85 × 101 1.3 × 108 1 × 10−8

M. Pˇribyl et al. / Computers and Chemical Engineering 30 (2006) 674–685

679

the adaptation interval
t:   m j  1 |fi |l   = + 0.1 , nx j eri |fs hs | + 0.1 + 0.1 j=1

s=1

i = 1, 2, . . . , nx

(8)

where l is length of the domain and hs is length of the s-th finite element. The small addition constants 0.1 are used to avoid division by zero in the regions without any detectable gradients. We adopted four monitor function f j with the same weight in the case study: ∂Φ ∂2 Φ ∂pH f3 = , f2 = 2 , , ∂x ∂x ∂x ∂pOH f4 = . (9) ∂x The monitor functions detect gradients on the profiles of the model variables ( f 1 , f 3 , f 4 ) and also sudden changes of gradient of electric potential, which contributes to the stability of simulations (Kulkarni & Dudukovic, 1996). A new length of the finite elements hei is estimated by means of Eq. (10) where the equidistribution principle is applied: f1 =

hequi eri l . hei = nx s=1 ers hs

(10)

The symbol hequi denotes the equidistant size of a finite element obtained from hequi = l/nx. Let us note that i hei = l. For that reason, the element sizes are finally normalized (as will be explained later). If t < tstop , the moving gradients of concentration and/or electric potential can be probably detected. Hence we have to increase the element density in neighborhoods of the gradients in the next step. We compute an integer quantity k that is defined as the rounded natural logarithm of the ratio of lengths of the largest and the shortest finite element (the lengths are evaluated from Eq. (10)): 
 max(he ) k = round ln . (11) min(he ) The radius of the neighborhood d is then increased or decreased according to the parameter k because its high value corresponds with the presence of sharp gradients. For example, k = 11 at the start of dynamical simulation in the case study. If k ≤ 3, then we increase the parameter d, otherwise, the radius of neighborhood is decreased. The constant kinc and the kdec (see Fig. 4) are usually chosen equal to 1.1 and 0.5, respectively. Locations and lengths of the finite elements on the original mesh are given by vectors p and h, respectively. Now, a new vector of element sizes he is assigned to the original location vector p. The increase of the element density in the neighborhood d of large gradients is carried out in a cycle for n growing from 1 to k. All element positions pi , for which hei < [max(he )]/exp(n), are found. Then for each such i, all element positions p j for which |p j − pi | ≤ d are found. The length of each such corresponding element hej is modified according to hej = [max(he )]/exp(n), cf. Fig. 4.

Fig. 4. The algorithm of the dynamical mesh adaptation.

Because the sum of lengths of all finite elements differs from the required length of the domain l, e.g., particular elements can overlap, the new mesh has to be normalized, i.e. new vectors p and h have to be found. The normalization is provided by the standard FEMLAB procedure meshinit. The input parameters of this procedure are: the original vector p, the computed vector he and the maximum ratio of lengths of two

M. Pˇribyl et al. / Computers and Chemical Engineering 30 (2006) 674–685

680

Fig. 5. Histograms of distribution of the finite elements on the spatial domain when only one region of sharp gradients is detected. z is number of the finite elements per 25 ␮m interval of the spatial domain. (A) Distribution of the elements without expansion of the dense mesh in the neighborhood of the sharp gradient (t > tstop ). (B) Distribution of the elements with the expansion of the dense mesh in the neighborhood of the sharp gradient (t < tstop , d = 20 ␮m). nx = 1000.

D H+ , d2

DH+ F | max {
Φid } | . (RTd 2 )

adjacent elements (chosen 1.1 in the studied case). The procedure meshinit constructs continuous mesh on the computation domain with the element lengths approximately equal to he in new positions given by a new vector p. If the number of finite elements is not equal to nx (common case), a uniform relative shrinking and swelling of the finite elements he is used and the procedure meshinit is again applied. New mesh h has nx finite elements. The method of linear interpolation is finally employed to fit the solution on the new mesh. An example of a non-equidistant mesh constructed on the basis of the monitor function (9) is plotted in Fig. 5. Here, only one region with large gradients is observed at the position x ≈ 0.4 mm. If the time of simulation exceeds tstop (close to the steady state) then no expansion of the dense mesh in the neighborhood d is applied. More than 60% of the finite elements are concentrated in the region 25 ␮m wide close to the location of the sharp gradients (see Fig. 5A). The other elements are almost uniformly distributed on the rest of the domain. When the gradients move (t < tstop ), the above described expansion of the dense mesh in a neighborhood of the sharp gradients (d = 20 ␮m chosen in this example) is applied. About 50% of the elements are concentrated in the 25 ␮m region with the sharpest gradients (Fig. 5B), which is a bit less than in the case without the expansion. However, more than 40% of the elements can be found in the two adjacent regions. Less than 10% of the elements are distributed on the rest of the domain. Hence, the expansion of the dense mesh can facilitate the simulation when the detected gradients move in time.

trD =

4.4. Estimation of the adaptation time interval

From the numerical point of view, dynamical simulation of the electrolyte diode in the closed mode is the most difficult problem because the gradients of the model variables are large and always presented. Hence, we have studied capabilities of the mesh adaptation during dynamical simulation of the closed mode. Spatial profiles of pH are plotted in Fig. 6A for selected times. A profile just after the initial condition, a profile several seconds after the start of simulation and a profile close to the steady state are shown by solid, dashed and dotted line, respectively.

The length of a new adaptation interval
t is estimated from the characteristic transport times of an electrolyte component according to the following equation:
t =

1 , max{trD , trE }

(12)

where the reciprocal diffusion and electromigration transport times, respectively, can be estimated from:

trE =

(13)

We have chosen hydrogen ions as the representative of the electrolyte in the case study. The symbol
Φid is a difference of electric potential on a segment i of the equidistant spatial mesh. Each segment i has the length d. The difference of electric potential is then computed by means of linear interpolation when the computed spatial profile of electric potential is interpolated on an equidistant mesh with the spatial step d. 5. Results and discussion The electrolyte diode can work in an open or in a closed mode. The open mode is characterized by a high electrolyte conductivity in the microcapillary and a negative value of electric potential imposed on the R boundary (see Fig. 1A). Steep concentration and electric potential gradients are observed in the closed mode, which is characterized by a positive value of the electric potential on the R boundary. We have divided the results section in several subsections. First, we discuss capabilities of the proposed algorithm to adapt the spatial mesh of the finite elements. We also focus on the effect of number of the finite elements. Then examples of dynamical behavior of the electrolyte diode in the both open and closed modes will be shown. Finally, response of the electrolyte diode to fast changes of electric potential on the R boundary will are presented. 5.1. Mesh characteristics

M. Pˇribyl et al. / Computers and Chemical Engineering 30 (2006) 674–685

681

Fig. 6. (A) Spatial profiles of pH value at selected times of the dynamical simulation (solid line t = 3 × 10−5 s, dashed line t = 1 s, dotted line t = 100.65 s). (B) Dependences of the logarithmic ratio of the length of spatial element and the length of the shortest element on the spatial index of the element. (C) The spatial profiles of the relative error of the pH value for the same times as in A. Number of elements nx = 500. (D) The spatial profiles of the relative error of the pH value for the same times as in A. Number of elements nx = 1000; ΦR = 10 V.

Typically 1000 finite elements were used for the simulations and finite elements are indexed from the L to the R boundary. Dependence of the relative logarithmic element size (the ratio of an element size and size of the shortest element) on its index is plotted in Fig. 6B. We have found that this logarithmic ratio almost attains the value of 10 for the central indexes of the finite elements (the solid line in Fig. 6B). This maximum corresponds with the flat central region. Most of the elements of the small size are located almost at the boundaries where the large gradients are observed (Fig. 6A). Let us note that the largest finite element is more than 104 times longer than the shortest element. Two regions with small element sizes are detected on the dashed line of indexes. These regions are located in place of the highest pH gradients (Fig. 6A). Near the steady state (dotted line), most of the short elements lie at the place of the sharp pH drop and the longest elements are located at the boundaries. 5.2. Effect of number of finite elements We have tested the effect of the number of finite elements on the evolution of relative error on selected pH profiles (Fig. 6A). We have computed these profiles for nx = 500, 1000, and 2000. Profiles of the relative error computed for 500 and 1000 elements are plotted in Fig. 6C and D, respectively. The solution computed on 2000 elements is considered to be the correct solution. The largest errors are found at the locations of sharp pH gradients but

they are of the order of 10−4 . Hence, 500 elements are enough to get the numerical solution with a high precision. We have used 1000 elements to avoid possible instabilities of the numerical solver for computation of transients. Consistent initial profiles with a high accuracy were not found for the number of finite elements equal to 250 and less. 5.3. Open mode (Fig. 7) The potassium and chloride ions are transported by the electromigration mechanism into the microcapillary in the open mode (Fig. 7B). The concentration gradients (Fig. 7A and B) and the gradients of electric potential (Fig. 7C and D) quickly vanish as the electric conductivity grows. Monotonous profiles of the model variables without extremely large gradients are observed near the steady state. The observed intensity of the electric field is of the order of 102 –103 V m−1 . These values are smaller than those used in the capillary electrophoresis or biochip applications, e.g., Dodge, Fluri, Verpoorte, and de Rooij (2001). The parameter k (see Section 4) rapidly decreases in the course of dynamical simulation, which results in an increase of both the radius of the neighborhood and the adaptation interval. The mesh expansion in the neighborhood of gradients can be stopped at time 1 s after the start of the simulation (for the given set of parameters). Typical computer time of such a simulation with the chosen relative and absolute tolerances is about

682

M. Pˇribyl et al. / Computers and Chemical Engineering 30 (2006) 674–685

Fig. 7. The electrolyte diode in the open mode. (A) Spatial profiles of pH. (B) Spatial profiles of concentration of potassium ions. (C) Spatial profile of electric potential. (D) Spatial profiles of intensity of the external electric field, solid line t = 0.1071 s, dashed line t = 0.849 s, dotted line t = 19.65 s, dash-dotted line t = 300.7 s; ΦR = −1 V.

20 min, which is comparable those found by Huang et al. (2002) and Wang et al. (2004a, 2004b). 5.4. Closed mode (Fig. 8) Hydrogen and hydroxyl ions migrate from the reservoirs in the microcapillary while potassium and chloride ions are kept in the reservoirs in the closed mode (see Fig. 8A). Thus a thin layer of low conductivity is formed in the place of water recombination reaction. Sharp gradients of the model variables are observed on the entire interval from the start of simulation to the steady state. Moreover, the gradient of electric potential in the place of low conductivity grows in time and moves to the steady state position (Fig. 8D). The observed gradient of electric potential almost reaches the value 106 V m−1 for the given set of parameters. These gradients exceed gradients of electric potential normally used in microfluidics but they can be found, e.g., in cell membranes (a typical value 1.2 × 107 V m−1 , e.g., Alberts et al., 2003). The adaptation parameter k usually remains higher than 3. Hence, radius of the neighborhood is short and the mesh adaptation has to be carried out with a higher frequency. It slows down the simulation procedure and the computational time can exceed several hours. On the other hand, the proposed algorithm is able to successfully simulate dynamical systems with such extreme gradients of electric potential. Dynamical simulations of the both closed and open mode are in agreement the steady state results received by Lindner et al. (2002).

5.5. Transients (Fig. 9) We have also tested the proposed algorithm for simulation of transient behavior of the electrolyte diode. Electric potential in the R reservoir is changed with a period 100 s. The values of electric potential are 0, 0.1, 1, 10, 1, 0.1, 0, −0.1, −1, −10, −1, −0.1, 0 V in this order (see Fig. 9A). To avoid sharp jumps, we have considered a fast linear change of electric potential on the R boundary after each 100 s period dΦR ΦR,a − ΦR,b , = dt
ttr

(14)

on the time interval
ttr . There the symbols ΦR,a and ΦR,b denote the value of electric potential after and before the change, respectively, and
ttr is the time duration of the jumps. We use
ttr = 0.1 s in the case study. As the electric potential is growing on the R boundary, the electrolyte diode is becoming closed. Electric current density does not decrease below ≈−300 A m−2 (see Fig. 9B). Only a short-time decrease of the current density is observed immediately after a change of electric potential, e.g., at time 200 s when ΦR is switched from 0.1 to 1 V. This is because the electrolyte diode is still not fully closed under such conditions. An increase of the imposed potential difference causes both an increase of electric current density according to the Ohm’s law and a change in the intensity of ionic transport. But a thin layer with low conductivity is quickly restored and electric current density returns near the value before the change of electric potential. In the open mode (t ≥ 700 s), the increase of the electric current density is

Fig. 8. The electrolyte diode in the closed mode. (A) Spatial profiles of pH. (B) Spatial profiles of concentration of potassium ions. (C) Spatial profile of electric potential. (D) Spatial profiles of intensity of the external electric field, solid line t = 9.453 × 10−2 s, dashed line t = 0.8451 s, dotted line t = 4.9381 s, dash-dotted line t = 149.6 s; ΦR = 1 V.

M. Pˇribyl et al. / Computers and Chemical Engineering 30 (2006) 674–685

683

Fig. 9. Non-stationary operation of the electrolyte diode. (A) The time course of the electric potential on the R boundary. (B) The time course of the electric current density. (C) The time course of the position of max(∂pH/∂x) in the closed mode.

proportional to the imposed potential difference according to the Ohm’s law because a high electric conductivity remains in the entire microchannel. The obtained volt–ampere characteristics of the electrolyte diode is in agreement with the work reported earlier (Lindner et al., 2002). We have also studied the effect of imposed potential difference in the closed mode on the position of the layer with low conductivity (Fig. 9C). The position of the layer is given by different mobilities of the ions, cf. Eq. (1). Because diffusivities of potassium and chloride ions are almost the same, the change in the layer position is given by the difference between diffusion coefficient of hydrogen and hydroxyl ions. We can observe that the position of the layer moves closer to the L compartment in the closed mode. When an electric field of a given orientation is applied, significance of the electromigration term in Eq. (1) grows. Because the mobility of hydrogen ions is higher than that of hydroxyl ions, the hydrogen ions migrate at a longer distance from the R reservoir. 5.6. Commentary on the method Different positive values of electric potential (ΦR , closed mode, nx = 1000) were imposed on the R boundary. It was found that the dynamical simulation can be carried out if ΦR is not higher than ≈300 V. For the case of 300 V, the maximal gradients ∂pH/∂x and ∂Φ/∂x were approximately 2 × 108 m−1 and 1.2 × 107 V m−1 at the beginning (the initial conditions), respectively, and 1 × 107 m−1 and 1.0 × 107 V m−1 at the end (close to the steady state), respectively. Moving fronts like a flame are usually characterized by extremely large spatial gradients of modeled variables. From this point of view, the proposed method should be also useful for dynamical modeling of such problems. The method parameters (nx, dmin , dmax , tstop ,
tmax ) have to be adjusted properly. However, a particular problem may require very small radius d, which finally results in frequent refinement of the mesh. The computation time then can substantially increase. In such cases, a method that introduces curvilinear coordinates, e.g., Adjerid and Flaherty (1986) and Thomson (1985), can be more favorable. There are two problems during the initial stages of simulations. Consistent steady state profiles of hydrogen ions, hydroxyl ions and electric potential have to be initially obtained on an unrefined mesh. It can lead to divergence problems of the New-

ton method if the number of mesh points at the localized gradients is low. The time solver can fail during the initial stage of time integration due to minimal time step underflow. This difficulty can be for example overcome by means of a smaller value of the parameter d and/or a higher value of parameter nx (number of mesh points). We typically used 1000 mesh points in our case study. Approximately 300 mesh points are necessary to get solution. The use of 250 mesh points leads to the solver failure. The monitor functions should be selected to monitor all significant gradients in the system. For example, first or second derivatives of concentration of potassium and chloride ions can be also selected as the monitor functions in our study. Their inclusion is redundant because their behavior copies the behavior of the other used monitor functions. Selection of the monitor functions always depends on a particular system. The computation time demands of the method are compared for various arrangements of the electrolyte diode (see Table 2). The solver parameters are nx = 1000, dmin = 2 ␮m, dmax = 100 ␮m, tend = 18 s, tstop = 0.5 s,
tmax = 2 s. All simulations were carried out on AMD 2.2 GHz system. Resulting CPU time demands are comparable with the results of an extensive comparative study made by Wang et al. (2004a, 2004b). The problem 5 in their study contains a system of four parabolic PDEs with reaction, convection, and diffusion terms. The studied system is characterized by formation of moving large gradients. The authors employed eight solvers. They observed CPU time demands in order 102 –104 s for the simulation time tend = 18 s when a high precision of computations was used. Bieniasz (2004a, 2004b) reported the time demands of a developed adaptive technique to compute the electric potential drop across a membrane system (Example 3 of the report). The mathemati-

Table 2 CPU demands of the method Diode mode

ΦR (V)

tend (s)

CPU demands (s)

Number of refinements

Open Open Neutral Closed Closed

−10 −1 0 1 10

18 18 18 18 18

1784 1302 2778 3217 4316

38 20 20 20 47

684

M. Pˇribyl et al. / Computers and Chemical Engineering 30 (2006) 674–685

cal model consists of two parabolic PDEs and one elliptic PDE. The problem was computed within 8427 s. A slow computer was used in the study (Intel 0.2 GHz) but the used tolerances and the observed gradients were significantly weaker than in our study. The model system studied in our work consists of four parabolic PDEs and one elliptic PDE. Although the used relative tolerance of the time solver is 10−8 and the absolute tolerances for hydroxyl and hydrogen ions are 10−10 mol m−3 , the CPU time demands are still comparable with the results by Bieniasz (2004a, 2004b) and Wang et al. (2004a, 2004b). 6. Conclusions We have proposed an empirical algorithm for mesh adaptation on a spatially one-dimensional domain. The algorithm combines the developed procedures of r-refinement and a standard FEMLAB procedure for the time integration. The developed method uses an expansion of the dense mesh in neighborhoods of the localized sharp gradients of the model variables. The adaptation method has been tested on a system of parabolic–elliptic PDEs that describes the electrolyte diode. Formation of extremely large and moving gradients of concentrations and electric potential is the main feature of the electrolyte diode. We have found that the developed solver is able to solve this problem on a wide range of imposed potential differences. The CPU time requirements are fully comparable with other mesh adaptation solvers. Simple empirical procedures of mesh adaptation and their easy implementation in standard dynamic solvers are important advantages of the proposed method. Dynamical simulations of the electrolyte diode help to understand processes that result in formation of a layer with low conductivity. Dynamics of development of both open and closed mode has revealed differences in the formation of qualitatively different patterns of the spatial profiles. The current–voltage characteristics of the electrolyte diode have been computed and the corresponding transient behavior has been described. The effect of the imposed difference of electric potential on the position of the layer with low electric conductivity in the closed mode has been explained by means of various mobilities of the electrolyte ions. We have verified stability of the solution obtained by the steady state modeling in Lindner et al. (2002). There are many applications in microfluidics where electrolytes and ions are transported by means of electric field, among others addressing in DNA chips, separation based on capillary electrophoresis or isoelectrical focusing, etc. Formation of large gradients of electric potential and concentrations of the electrolyte components is often possible (e.g., at the microchannel boundaries). We believe that the proposed algorithm of the mesh adaptation can be useful in research and development of such applications. Acknowledgement The authors thank to the Ministry of Education, Youth, and Sport of the Czech Republic (grant nos. 1K03003 and MSM 6046137306) for supporting this project.

References Adjerid, S., & Flaherty, J. E. (1986). A moving-mesh finite element method with local refinement for parabolic partial differential equations. Computer Methods in Applied Mechanics and Engineering, 55, 3. Alberts, B., Johnson, A., Lewis, J., Raff, M., Bray, D., Hopkin, K., Roberts, K., & Walter, P. (2003). Essential cell biology. New York: Garland Science/Taylor & Francis Group. Almeida, R. C., & Silva, R. S. (2002). An hp-adaptive strategy in a Petrov–Galerkin method for convection–diffusion problems. Communications in Numerical Methods in Engineering, 18, 203. Bieniasz, L. K. (2000). Use of dynamically adaptive grid techniques for the solution of electrochemical kinetic equations. Part 5. A finitedifference, adaptive space/time grid strategy based on a patch-type local uniform spatial grid refinement, for kinetic models in onedimensional space geometry. Journal of Electroanalytical Chemistry, 481, 115. Bieniasz, L. K. (2004a). Use of dynamically adaptive grid techniques for the solution of electrochemical kinetic equations. Part 14: Extension of the patch-adaptive strategy to time dependent models involving migration–diffusion transport in one-dimensional space geometry, and its application to example transient experiments described by Nernst–Planck–Poisson equations. Journal of Electroanalytical Chemistry, 565, 251. Bieniasz, L. K. (2004b). Use of dynamically adaptive grid techniques for the solution of electrochemical kinetic equations. Part 15: Patch-adaptive simulation of example transient experiments described by Nernst–Planckelectroneutrality equations in one-dimensional space geometry. Journal of Electroanalytical Chemistry, 565, 273. Blom, J. G., & Zegeling, P. A. (1994). Algorithm 731: A moving grid interface for systems of one-dimensional time-dependent partial differential equations. ACM Transactions on Mathematical Software, 20, 194. Cao, W., Huang, W., & Russell, R. D. (2003). Approaches for generating moving adaptive meshes: Location versus velocity. Applied Numerical Mathematics, 47, 121. Coimbra, M. C., Sereno, C., & Rodrigues, A. (2004). Moving finite element method: Applications to science and engineering problems. Computers and Chemical Engineering, 28, 597. D03PPF NAG Fortran library routine document. http://www.nag.co.uk/ numeric/fl/manual19/pdf/D03/d03ppf fl19.pdf. Dodge, A., Fluri, K., Verpoorte, E., & de Rooij, N. F. (2001). Electrokinetically driven microfluidic chips with surface-modified chambers for heterogeneous immunoassays. Analytical Chemistry, 73, 3400. Drake, R., & Manoranjan, V. S. (1996). A method of dynamic mesh adaptation. International Journal for Numerical Methods in Engineering, 39, 939. Dutta, P., & Beskok, A. (2001). Analytical solution if combined electroosmotic/pressure driven flow in two-dimensional straight channels: Finite Debye layer effects. Analytical Chemistry, 73, 1979. FEMLAB 2.3 ref. manual FEMLAB 2.3, reference manual. Fukuba, T., Yamamoto, T., Naganuma, T., & Fujii, T. (2004). Microfabricated flow-through device for DNA amplification—Towards in situ gene analysis. Chemical Engineering Journal, 101, 151. Han, J., & Singh, A. K. (2004). Rapid protein separations in ultrashort microchannels: Microchip sodium dodecyl-polyacrylamide gel electrophoresis and isoelectrical focusing. Journal of Chromatography A, 1049, 205. Huang, W., Ren, Y., & Russell, R. D. (1994). Moving mesh partial differential equations (MMPDES) based on the equidistribution principles. SIAM Journal on Numerical Analysis, 31, 709. Huang, W., & Russell, R. D. (1996). A moving collocation method for solving time dependent partial differential equation. Applied Numerical Mathematics, 20, 101. Huang, W., Zheng, L., & Zhan, X. (2002). Adaptive moving mesh methods for simulating one-dimensional groundwater problems with sharp moving fronts. International Journal for Numerical Methods in Engineering, 54, 1579.

M. Pˇribyl et al. / Computers and Chemical Engineering 30 (2006) 674–685 Keast, P., & Muir, P. (1991). Algorithm 688. EPDCOL: A more efficient PDECOL code. ACM Transactions on Mathematical Software, 17, 153. Kulkarni, M. S., & Dudukovic, M. P. (1996). A robust algorithm for fixedbed reactors with steep moving temperature and reaction fronts. Chemical Engineering Science, 51, 571. Lim, Y. I., Le Lann, J. M., & Joulia, X. (2001). Accuracy, temporal performance and stability comparisons of discretization methods for the numerical solution of partial differential equations (PDEs) in the presence of steep moving fronts. Computers and Chemical Engineering, 25, 1483. Lindner, J., Snita, D., & Marek, M. (2002). Modelling of ionic systems with a narrow acid base boundary. Physical Chemistry Chemical Physics, 4, 1348. Liu, F., & Bhatia, S. K. (1999). Computationally efficient solution techniques for adsorption problems involving steep gradients in bidisperse particles. Computers and Chemical Engineering, 23, 933. Liu, F., & Bhatia, S. K. (2001). Solution techniques for transport problems involving steep concentration gradients: Application to noncatalytic fluid solid reactions. Computers and Chemical Engineering, 25, 1159. Merkin, J. H., Simon, P. L., & Noszticzius, Z. (2000). Analysis of the electrolyte diode. Electro-diffusion and chemical reaction within a hydrogel reactor. Journal of Mathematical Chemistry, 28, 43. Mitnik, L., Novotny, M., Felten, C., Buonocore, S., Koutny, L., & Schmalzing, D. (2001). Recent advances in DNA sequencing by capillary and microdevice electrophoresis. Electrophoresis, 22, 4104. Probstein, R. F. (1994). Physicochemical hydrodynamics: An Introduction. New York: Wiley and Sons.

685

ˇ Pˇribyl, M., Snita, D., Hasal, P., & Marek, M. (2004). Modeling of electricfield driven transport processes in microdevices for immunoassay. Chemical Engineering Journal, 101, 303. ˇ Pˇribyl, M., Snita, D., & Marek, M. (2005). Nonlinear phenomena and qualitative evaluation of risk of clogging in a capillary microreactor under imposed electric field. Chemical Engineering Journal, 105, 99. Righetti, P. G., Gelfi, C., & D’Acunto, M. R. (2002). Recent progress in DNA analysis by capillary electrophoresis. Electrophoresis, 23, 1361. Roussel, O., Schneider, K., Tsigulin, A., & Bockhorn, H. (2003). A conservative fully adaptive multiresolution algorithm for parabolic PDEs. Journal of Computational Physics, 188, 493. Schmalzing, D., Koutny, L. B., Taylor, T. A., Nashabeh, W., & Fusch, M. (1997). Immunoassay for thyroxine (T4) in serum using capillary electrophoresis and micromachinated devices. Journal of Chromatography B, 697, 175. Thiele, A., Falk, T., Tobiska, L., & Seidel-Morgenstern, A. (2001). Prediction of elution profiles in annular chromatography. Computers and Chemical Engineering, 25, 1089. Thomson, J. F. (1985). A survey of dynamically-adaptive grids in the numerical solution of partial differential equations. Applied Numerical Mathematics, 1, 3. Wang, R., Keast, P., & Muir, P. (2004a). A comparison of adaptive software for 1D parabolic PDEs. Journal of Computational and Applied Mathematics, 169, 127. Wang, R., Keast, P., & Muir, P. (2004b). A high-order global spatially adaptive collocation method for 1-D parabolic PDEs. Applied Numerical Mathematics, 50, 239.