Numerical simulation of a backward-facing step flow in a microchannel with external electric field

Abstract

A backward-facing step flow in the microchannel with external electric field was investigated numerically by a high-order accuracy upwind compact difference scheme in this work. The Poisson–Boltzmann and Navier–Stokes equations were computed by the high-order scheme, and the results confirmed the ability of the new solver in simulation of micro-scale electric double layer effects. The flow fields were displayed for different Reynolds numbers; the positions of the vortex saddle point of model with external electric field and model without external electric field were compared. The average velocity increases linearly with the electric field intensity; however, the Joule heating effects cannot be neglected when the electric field intensity increases to a certain level.

Introduction

With the development of micro- and nano-technology in recent years, the study on electric double layer (EDL) is changing from theoretical analysis to applied research. The rapid development of micro-fabrication technology has enabled a variety of microfluidic systems, consisting of the microchannel, micro-valves, micro-pumps, and micro-mixer, which have been effectively used in medical, pharmaceutical, defense, and environmental monitoring. The micro- and nano-electroosmotic flow is thus enabled to become a cutting-edge research topic. In microfluidic systems, the electrolyte solution under an applied electric field is driven by the Coulomb force, and the internal liquid molecules start to move; the EDL effect appears on the solid–liquid interface. ¹ The rapid development of microfluidic technologies has got a large number of achievements, such as microelectric mechanics systems (MEMS), Biochips, and biosensor systems.^2–5 This requires people to have profound understanding of physical characteristics of microfluidic systems so as to manufacture new devices and expand the usage.

Due to the restrictions of geometry size, it is difficult to measure electroosmotic transport phenomena of microchannel flows precisely by experimental methods; therefore, efficient and high-order accuracy numerical scheme of the electroosmotic flow is important in numerical simulation. Over the last few decades, many researches have been done on electroosmotic transport phenomena in microchannels. Hunter ⁶ reported that the flow behavior will be greatly affected by the close channel surface because of the existence of EDL. Bier et al. ⁷ used a computer to simulate the electrophoresis of the sample transfer process. In recent years, researchers have employed different numerical methods, such as the finite difference method, ⁸ finite volume method, ⁹ and finite element method, ¹⁰ to simulate the electroosmotic phenomenon in microchannels with different structures. For pressure-driven flows in microchannels, Wang and Wu ¹¹ reported that the periodical flow velocity and the Flow-Induced Electric Field (FIEF) strongly depend on the frequency Reynolds number; Ban et al. ¹² reported that average electric conductivity of electrolyte inside the channel increased significantly as the EDL overlaps. Zhu et al.^13,14 demonstrated a model to investigate the transient heat and liquid moisture transfer through porous fibrous media with consideration of EDL. On the shape of microchannels, Taylor et al. ¹⁵ studied the impact of surface roughness on fluid flow, and pointed out that with the decrease of channel-scale, the roughness on fluid flow becomes more important. However, the influence of protrusions in microchannels needs to be further investigated with detailed examination of electric conductivity of the electrolyte in order to establish a better understanding of the effect of EDL.

By using a high-order accuracy upwind compact difference scheme, the current research focuses on the simulation of the EDL in a microchannel with a backward-facing step. The remaining sections are arranged as follows: section ‘Formulations and modeling’ gives a brief description about the formulas of EDL and the flow field in microchannel. Section ‘Numerical schemes’ describes high-order accuracy numerical schemes. Numerical results and discussions are presented in section ‘Results and discussion’, and conclusions are drawn in section ‘Conclusions’.

Formulations and modeling

EDL effect

In microfluidic systems, the surfaces of the channel walls are polarized to bring electrostatic charge when wall materials are in contact with the electrolyte solution. For example, when a wall made of glass is put in the water, the chemical reaction produces net negative charges. The interaction between ionizing solution and the static charge on the surface of the insulation affects the distribution of ions in the electrolyte solution, and different charged ions gathered near the wall quickly, resulting in the reorganization of the solid–liquid interface near the liquid ion to form a thin liquid with a net charge layer, which is called the EDL.

Figure 1 shows the backward-facing step flow in the microchannel with external electric field. In the microscopic scale, the step can be seen as a small protrusion of the wall. The external electric field is applied to both the upper and lower wall. In microscopic cases, the continuum assumption is still applicable in the liquid flow when electric field is applied.^1,16–20 The electrical potential $\psi$ and the net charge density ${\rho }_e$ is governed by Poisson’s equation for electrostatic field

$${\nabla }^2\psi =-\frac{{\rho }_e}{\epsilon {\epsilon }_0}$$

(1)

where $\epsilon$ is the permittivity of the solution and ${\epsilon }_0$ is the vacuum permittivity. Assume that the space distribution of the ion concentration in the electrolyte solution meets the Boltzmann distribution, then the net charge density can be write as

$${\rho }_e=e\sum z_i{n}_{i\infty }\exp (-\frac{z_{i}e\psi }{k_{b}T})$$

(2)

where $n_{i\infty }$ presents the ion concentration in electrolyte solution, $z_i$ is the ion valence of ith kind of ion, $e$ denotes the electron charge constant, $k_b$ is Boltzmann constant, and T is the absolute temperature.

Figure 1.

A microchannel with a backward-facing step.

For the symmetric binary electrolyte solution used in this work, the Poisson–Boltzmann equation can be reduced to

$${\nabla }^2\psi =\frac{2e{n}_{i\infty }}{\epsilon {\epsilon }_0}\sinh \left(\frac{e\psi }{k_{b}T}\right)$$

(3)

Flow field in microchannels

When an electric field is applied at the top and bottom of the microchannel, electrostatic charges in the EDL are driven by the electric force and start to move. The movement of charges drives the movement of liquid within the channel because of the viscosity, thus forming the electroosmotic flow. Because the length of the microchannel is usually much larger than its width, the flows in the microchannel can be well approximated by a fully developed laminar flow except at the entrance and the exit of the channel.

Let $\partial \Omega$ be the boundary of a polyhedral domain $\Omega$ . $H^1\left(\Omega \right)$ is the first-order Sobolev space and $L^2\left(\Omega \right)$ is be the space of second power summable functions on $\Omega$ . Driven by the Coulomb force, the Navier–Stokes equation for incompressible viscous fluids in a two-dimensional uniform microchannel is to find $(u,p)\in H^1{\left(\Omega \right)}^3\times L^2\left(\Omega \right)$ such that for any $t\in (0,T)$ , the following set of equations hold

$$\{\begin{array}{c}\multicolumn{1}{c}{\rho (\frac{\partial \stackrel{\rightharpoonup }{u}}{\partial t}+\stackrel{\rightharpoonup }{u}·\nabla \stackrel{\rightharpoonup }{u})=-\nabla P+\mu {\nabla }^2\stackrel{\rightharpoonup }{u}+\stackrel{\rightharpoonup }{E}{\rho }_e}\\ \multicolumn{1}{c}{\nabla ·\stackrel{\rightharpoonup }{u}=0}\end{array}\text{in}\Omega \times (0,T)$$

(4)

where $\stackrel{\rightharpoonup }{u}$ is the velocity, $\mu$ is the kinematic viscosity, $\rho$ is the density of the solution, $E$ is the strength of the electric field applied to the microchannel, and it is a constant in this work.

Let $\zeta$ denote the vorticity and $\varphi$ denote the stream function, the vortex-stream form of equation (4) can be written as

$$\{\begin{array}{c}\multicolumn{1}{c}{\frac{\partial \zeta }{\partial t}+u\frac{\partial \zeta }{\partial x}+v\frac{\partial \zeta }{\partial y}=\mu (\frac{{\partial }^2\zeta }{\partial x^2}+\frac{{\partial }^2\zeta }{\partial y^2})+\frac{E}{\rho }\left(\frac{\partial {\rho }_e}{\partial x}\right)}\\ \multicolumn{1}{c}{\frac{{\partial }^2\varphi }{\partial x^2}+\frac{{\partial }^2\varphi }{\partial y^2}=-\zeta }\end{array}$$

(5)

where $u=\partial \varphi /\partial y$ and $v=-\partial \varphi /\partial x$ are the velocity in $x$ and $y$ directions, respectively.

Let $H$ be the characteristic length, $U$ be the characteristic velocity, and $T$ be the characteristic time, then the vortex-stream form governing equations can be nondimensionalized by

$$\overline{x}=\frac{x}{H},\overline{y}=\frac{y}{H},\overline{\psi }=\frac{e\psi }{k_{b}T},K=\frac{2n_{i\infty }{k}_{b}T}{\rho U^2},\overline{E}=\frac{H\text{Re}}{k_{b}T/e}E,\overline{\varphi }=\frac{\varphi }{\mu },\overline{\zeta }=\frac{\zeta }{\mu },$$

(6)

where

$$\kappa =\sqrt{\frac{2e^2{n}_{i\infty }}{\epsilon {\epsilon }_0{k}_{b}T}}$$

(7)

is the Debye–Huckel coefficient and $\text{Re}$ is the Reynolds number.

The dimensionless form of equation (3) and equation (5) is

$$\{\begin{array}{c}\multicolumn{1}{c}{\frac{{\partial }^2\overline{\psi }}{\partial {\overline{x}}^2}+\frac{{\partial }^2\overline{\psi }}{\partial {\overline{y}}^2}={\kappa }^2{H}^2\sinh \overline{\psi }}\\ \multicolumn{1}{c}{\frac{\partial \overline{\zeta }}{\partial t}+\overline{u}\frac{\partial \overline{\zeta }}{\partial \overline{x}}+\overline{v}\frac{\partial \overline{\zeta }}{\partial \overline{y}}=\upsilon (\frac{{\partial }^2\overline{\zeta }}{\partial {\overline{x}}^2}+\frac{{\partial }^2\overline{\zeta }}{\partial {\overline{y}}^2})+\frac{\overline{E}}{\rho }\left(\frac{\partial {\rho }_e}{\partial \overline{x}}\right)}\\ \multicolumn{1}{c}{\frac{{\partial }^2\overline{\varphi }}{\partial {\overline{x}}^2}+\frac{{\partial }^2\overline{\varphi }}{\partial {\overline{y}}^2}=-\overline{\zeta }}\end{array}$$

(8)

where $\overline{u}=\partial \overline{\varphi }/\partial \overline{y}$ and $\overline{v}=-\partial \overline{\varphi }/\partial \overline{x}$ .

The non-slip boundary conditions are

$$\overline{\psi }=\overline{\xi }=\frac{e\xi }{k_{b}T}$$

(9)

where $\xi$ is the potential of the upper and lower wall applied on all the end-walls.

Numerical schemes

A fifth-order discretization of the first derivatives

As is known, the general compact differential expression is

$$\sum _l{\alpha }_l{F}_{j+l}=\sum _l{\beta }_l{u}_{j+l}$$

(10)

In this work, a differential approximation of five points is considered ²¹

$${\alpha }_1{F}_{j+1}+{\alpha }_0{F}_j+{\alpha }_{-1}{F}_{j-1}={\beta }_{-2}{u}_{j-2}+{\beta }_{-1}{u}_{j-1}+{\beta }_0{u}_j+{\beta }_1{u}_{j+1}+{\beta }_2{u}_{j+2}$$

(11)

Let $F/\Delta x=\partial u/\partial x$ be the difference approximation of first derivate. Set value for the free variables in equation (11) and the coefficients of $\alpha$ and $\beta$ in equation (11) are determined as follows

$$\{\begin{array}{c}\multicolumn{1}{c}{{\alpha }_0=\frac{3}{5},{\alpha }_{-1}=\frac{2}{5},{\alpha }_1=0}\\ \multicolumn{1}{c}{{\beta }_{-2}=-\frac{1}{20},{\beta }_{-1}=-\frac{44}{60},{\beta }_0=\frac{36}{60},}\\ \multicolumn{1}{c}{{\beta }_1=\frac{12}{60},{\beta }_2=-\frac{1}{60}}\end{array}$$

(12)

Then equation (11) can be written as

$$\frac{3}{5}{F}_j+\frac{2}{5}{F}_{j-1}=\frac{1}{60}{\delta }_x^{-}(3u_{j-1}+47u_j+11u_{j+1}-u_{j+2})$$

(13)

where ${\delta }_x^{-}{u}_j=u_j-u_{j-1}$ .

The Taylor series of a discrete function ${\text{T}}_{j+1}$ at point $j$ is

$${\text{T}}_{j+1}={\text{T}}_j+\sum _{n=1}^{\infty }\frac{{\left(\Delta x\right)}^n}{n!}{\left(\frac{{\partial }^n\text{T}}{\partial x^n}\right)}_j$$

(14)

Applying equation (14) to both sides of equation (13) yields

$$\frac{F_j}{\Delta x}={\left(\frac{\partial f}{\partial x}\right)}_j+O{\left(\Delta x\right)}^5$$

(15)

As can be seen, $F_j/\Delta x$ is a fifth-order accuracy discretization of the first derivative $f_x$ .

A fourth-order discretization of the second derivatives

By using Taylor expansion in equation (14), a fourth-order accuracy discretization method for the second derivatives can be written as

$${\delta }_x^2{u}_j=u_{j+1}-2u_j+u_{j-1}=\Delta x^2{\left(\frac{{\partial }^{2}u}{\partial x^2}\right)}_j+\frac{1}{12}\Delta x^4{\left(\frac{{\partial }^{4}u}{\partial x^4}\right)}_j+O\left(\Delta x^6\right)$$

(16)

Assuming $S{\prime }_j=\Delta x^2{({\partial }^{2}u/\partial x^2)}_j$ , then

$$\begin{array}{c}\multicolumn{1}{c}{\Delta x^4{\left(\frac{{\partial }^{4}u}{\partial x^4}\right)}_j=\Delta x^2{\left[\frac{{\partial }^2}{\partial x^2}\left(\Delta x^2\frac{{\partial }^{2}u}{\partial x^2}\right)\right]}_j}\\ \multicolumn{1}{c}{=\Delta x^2{\left(\frac{{\partial }^{2}s{\prime }_j}{\partial x^2}\right)}_j}\\ \multicolumn{1}{c}{=S{\prime }_{j+1}-2S{\prime }_j+S{\prime }_{j-1}+O\left(\Delta x^6\right)}\end{array}$$

(17)

Substituting (17) into (16) yields

$${\delta }_x^2{u}_j=\frac{1}{12}S{\prime }_{j+1}+\frac{5}{6}S{\prime }_j+\frac{1}{12}S{\prime }_{j-1}$$

(18)

It is noted that $S{\prime }_j/\Delta x^2$ is a fourth-order accuracy discretization method of the second derivatives $u_{\mathrm{xx}}$ , as can be seen from equation (17).

Discretization of the Poisson’s equation

The stream function equations, pressure equation of governing equations are all Poisson’s equations, where $t$ serves as coefficient and the equations are for steady problems; therefore, the following equations are considered

$$\{\begin{array}{l}\begin{array}{cc}\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}=f\left(x,y\right)& \text{on}D\end{array}\\ \begin{array}{cc}u\left(x,y\right)=g\left(x,y\right)& \text{on}\partial D\end{array}\end{array}$$

(19)

where $D$ is a bounded region in the plane, with piecewise smooth boundaries $\partial D$ .

Expand $u_{i+1,j},u_{i,j+1},u_{i-1,j},u_{i,j-1},u_{i+1,j+1},u_{i-1{,}_{j+1}},u_{i-1,j-1},u_{i+1,j-1}$ at point $(i,j)$ with the Taylor expansion and let

$$\begin{array}{c}\multicolumn{1}{c}{T_1=u_{i+1,j}+u_{i,j+1}+u_{i-1,j}+u_{i,j-1}}\\ \multicolumn{1}{c}{=4u_{i,j}+h^2(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2})+\frac{h^4}{12}(\frac{{\partial }^{4}u}{\partial x^4}+\frac{{\partial }^{4}u}{\partial y^4})+O\left(h^4\right)}\end{array}$$

(20)

and

$$\begin{array}{c}\multicolumn{1}{c}{T_2=u_{i+1,j+1}+u_{i-1,j+1}+u_{i-1,j-1}+u_{i+1,j-1}}\\ \multicolumn{1}{c}{=4u_{i,j}+2h^2(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2})+\frac{h^4}{6}(\frac{{\partial }^{4}u}{\partial x^4}+\frac{{\partial }^{4}u}{\partial y^4})}\\ \multicolumn{1}{c}{+h^4\frac{{\partial }^{4}u}{\partial x^2\partial y^2}+O\left(h^4\right)}\end{array}$$

(21)

then we have

$$\begin{array}{c}\multicolumn{1}{c}{\frac{a{T}_1+b{T}_2+c{u}_{i,j}}{h^2}=\frac{(4a+4b+c)u_{i,j}}{h^2}}\\ \multicolumn{1}{c}{+(a+2b)(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2})+\frac{h^2}{12}(a+2b)(\frac{{\partial }^{4}u}{\partial x^4}+\frac{{\partial }^{4}u}{\partial y^4})+h^{2}b\frac{{\partial }^{4}u}{\partial x^2\partial y^2}+O\left(h^4\right)}\end{array}$$

(22)

From equation (19), ${\partial }^{4}u/\partial x^4$ and ${\partial }^{4}u/\partial y^4$ can be obtained. By substituting them into equation (22), it can be seen that to keep a fourth-order accuracy, the $u_{i,j}$ terms should vanish and the ${\partial }^{4}u/\partial x^2\partial y^2$ term should be the cross derivative term. That is to say, the following equations

$$\{\begin{array}{c}\multicolumn{1}{c}{4a+4b+c=0}\\ \multicolumn{1}{c}{-2a+8b=0}\end{array}$$

(23)

should be satisfied. Using one of the solutions (e.g. $a=4$ , $b=1$ , $c=-20$ ), the fourth-order compact difference scheme can be expressed as

$$\begin{array}{c}\multicolumn{1}{c}{4(u_{i+1,j}+u_{i,j+1}+u_{i-1,j}+u_{i,j-1})}\\ \multicolumn{1}{c}{+(u_{i+1,j+1}+u_{i-1,j+1}+u_{i-1,j-1}+u_{i+1,j-1})-20u_{i,j}}\\ \multicolumn{1}{c}{}& =6h^2{f}_{i,j}+\frac{h^4}{2}(\frac{f_{i+1,j}-2f_{i,j}+f_{i-1,j}}{\Delta x^2}+\frac{f_{i,j+1}-2f_{i,j}+f_{i,j-1}}{\Delta y^2})\hfill & \hfill \end{array}$$

(24)

A relaxation iteration method is employed to solve equation (24) by

$$\begin{array}{cc}\multicolumn{1}{c}{}& u_{i,j}^{k+1}=\frac{1}{5}(u_{i+1,j}^k+u_{i,j+1}^k+u_{i-1,j}^{k+1}+u_{i,j-1}^{k+1})\hfill \\ \multicolumn{1}{c}{+\frac{1}{20}(u_{i+1,j+1}^k+u_{i-1,j+1}^k+u_{i-1,j-1}^{k+1}+u_{i+1,j-1}^k)}\\ \multicolumn{1}{c}{}& -\frac{3}{10}{h}^2{f}_{i,j}^{k+1}-\frac{h^4}{40}(\frac{f_{i+1,j}^k-2f_{i,j}^{k+1}+f_{i-1,j}^{k+1}}{\Delta x^2}+\frac{f_{i,j+1}^k-2f_{i,j}^{k+1}+f_{i,j-1}^{k+1}}{\Delta y^2})\hfill \end{array}$$

(25)

Let $u_{i,j}^{k+1}=u_{i,j}^k+w\Delta u_{i,j}$ , where

$$\begin{array}{cc}\multicolumn{1}{c}{}& \Delta u_{i,j}=\frac{1}{5}(u_{i+1,j}^k+u_{i,j+1}^k+u_{i-1,j}^{k+1}+u_{i,j-1}^{k+1})\hfill \\ \multicolumn{1}{c}{+\frac{1}{20}(u_{i+1,j+1}^k+u_{i-1,j+1}^k+u_{i-1,j-1}^{k+1}+u_{i+1,j-1}^k)}\\ \multicolumn{1}{c}{}& -\frac{3}{10}{h}^2{f}_{i,j}^{k+1}-\frac{h^4}{40}(\frac{f_{i+1,j}^k-f_{i-1,j}^{k+1}}{2\Delta x^2}+\frac{f_{i,j+1}^k+f_{i,j-1}^{k+1}}{2\Delta y^2})-u_{i,j}^k\hfill \end{array}$$

(26)

In this work, the relaxation factor $w=1.2$ .

Time discretization scheme

A third-order Runge–Kutta method is employed to discretize the temporal space

$$\{\begin{array}{c}\multicolumn{1}{c}{{\zeta }_{i,j}^1={\zeta }_{i,j}^n+\Delta t*H_{i,j}}\\ \multicolumn{1}{c}{{\zeta }_{i,j}^2=\frac{3}{4}{\zeta }_{i,j}^n+\frac{1}{4}({\zeta }_{i,j}^1+\Delta t*H_{i,j})}\\ \multicolumn{1}{c}{{\zeta }_{i,j}^{n+1}=\frac{1}{3}{\zeta }_{i,j}^n+\frac{2}{3}({\zeta }_{i,j}^2+\Delta t*H_{i,j})}\end{array}$$

(27)

where

$$H_{i,j}={[-(u\frac{\partial \zeta }{\partial x}+\upsilon \frac{\partial \zeta }{\partial y})+v(\frac{{\partial }^2\zeta }{\partial x^2}+\frac{{\partial }^2\zeta }{\partial y^2})+\frac{E}{\rho }\left(\frac{\partial {\rho }_e}{\partial x}\right)]}_{i,j}$$

(28)

In order to keep the convergence of the iteration, $\Delta t$ should satisfy the following inequality

$$\Delta t_{i,j}\le \frac{2}{5}max(\frac{\left|u\right|}{\Delta x}+\frac{\left|v\right|}{\Delta y}+c\sqrt{\frac{1}{\Delta x^2}+\frac{1}{\Delta y^2}})$$

(29)

Results and discussion

In this work, KCl solution is used as the electrolyte solution in the microchannel, and the coefficients are given by Table 1.

Table 1.

The material properties.

Parameters	Values
$n_{i\infty }$	1.0 × 10−4 mol/L
$\rho$	1.0 × 103 kg/m3
$\epsilon$	80 C2·J−1·m−1
${\epsilon }_0$	8.854 × 10−12 C2·J−1·m−1
$k_b$	1.38 × 10−23 J·K−1
T	287 K
$e$	1.6 × 10−19 C

The influence of with and without EDL to a backward-facing step model under the high-precision numerical methods was first studied. The strength of electric field applied to the microchannel was 100 kV/m and the potential was 150 mV

By comparing the left and right diagrams of Figures 2 and 3, it can be seen that when Re = 10 and the EDL is applied, the saddle point position of the vortex moves to the left; the location of the vortices tends the backward-facing step and reattachment point of the vortex is closer to the backward-facing step. The length of the recirculation zone is smaller.

Figure 2.

Streamline of Re = 10 (left: no EDL; right: with EDL).

Figure 3.

Velocity vectors of Re = 10 (left: no EDL; right: with EDL).

Second, a higher Reynolds number was applied. As can be seen from Figures 4 and 5, when Re = 15 is applied, the trend of fluid motion is more clear; the vortex produced by the backward-facing step is bigger, and the location of the vortex is further from the overall steps. When the EDL is applied, the position of the saddle point moves to the left. Compared with the case of no EDL effects (the left diagram), both the location of the vortices and the reattachment point of the vortex are closer to the backward-facing step. The length of the recirculation zone is smaller.

Figure 4.

Streamline of Re = 15 (left: no EDL; right: with EDL).

Figure 5.

Velocity vectors of Re = 15 (left: no EDL; right: with EDL).

In order to study the influence of the electric field intensity, the external electric field applied to the microchannel was set to 50, 100, 150, and 200 kV/m; the potential was 150 mV and the Reynolds number was set to 15. Streamlines are given in Figure 6.

Figure 6.

Streamlines under different electric field intensities.

In Figure 6, it can be seen that the vortex moves to the left gradually with the increase of the electric field intensity, and the position of the saddle point becomes lower and closer to the step. The flow field moves toward the bottom, and the length of the recirculation zone gets smaller. To further confirm the influence of the external electric field, the relationship between the velocity and the strength of the electric filed is a concern, as is shown in Figure 7.

Figure 7.

The average velocity at the outlet versus the electric field intensity.

Figure 7 shows the relationship between the electric field intensity and the average velocity at the outlet. It can be seen with the increase of the electric field intensity that the speed increases linearly. It also should be noted that the electric field strength cannot increase indefinitely because once the field strength increases to a certain value, the Joule heating effects will too obviously be neglected.

Conclusion

Based on a high-accuracy upwind compact difference method, the Navier–Stokes equations of vortex-stream form and the Poisson–Boltzmann equation were solved in this work. A backward-facing step flow in the microchannel with external electric field was investigated numerically by the scheme; the results confirmed the ability of the new solver in simulation of micro-scale EDL effects. Several characters of the backward-facing step flow in the microchannel with external electric field are discovered:

1. In a backward-facing step flow in microchannels, EDL appears and the position of the vortex saddle point moves to the left in the presence of the external electric field. Compared with the case with no EDL effects, both the location of the vortices and the reattachment point of the vortex are closer to the backward-facing step, and the length of the recirculation zone is smaller. With an increase of the Reynolds number, the discrepancy in fluid motion becomes clearer.

2. The velocity increases linearly with the electric field intensity. However, the electric field intensity cannot increase indefinitely, as the Joule heating effects will too obviously be neglected when the electric field intensity increases to a certain level.

The computation process and the results show the effectiveness of the method and confirm the ability of the new code to simulate the EDL effect within microchannels. By using the high-order upwind compact scheme, the simulation is more accurate, and to confirm this is a future job for us.