1.   引　言

微流体装置在微机电系统、生物和化学传感器、药物递送芯片、溶剂分离装置和热控制系统等方面有很重要的应用. 电解质与微通道壁面发生接触时微通道中的化学成分和电解质在表面发生化学反应, 在微通道壁面与电解质之间进行电荷交换. 这些表面电荷通过影响电解质溶液壁面附近的阳离子和阴离子分布情况形成双电层(EDL). 当施加带电表面切线方向的电场时, EDL中的离子在电场力作用下发生定向移动. 此时可移动的离子沿电场方向运动, 由于流体黏性效应, 附近的流体质点随之运动, 电渗流(EOF)由此而形成^[1–4].

电渗流是电动运输中的一种主要现象. 1809年, Reuss^[5]观察水在外电场传输过程中的运动时发现了这一现象. 随后微流体装置中的电渗流问题成为研究的热点问题. 关于微管道内完全发展的牛顿流体的电渗流已经有了大量研究^[6–17].

许多重要流体的流动不能用一个简单的牛顿模型来描述. 非牛顿流体广泛存在于大多数生理流体和不同工业使用的流体中^[18,19]. 然而, 这些流变学模型大多没有考虑到微观结构的影响, 不足以模拟流体微元素的旋转. 而微极性流体表现出耦合应力, 流体的颗粒除了速度矢量外还具有独立的旋转矢量. Hoyt等^[20]在实验中发现, 在牛顿流体中添加极少量的聚合物, 物体在湍流中的摩擦阻力会降低30%—35%. Eringen等^[21]为了解释这一现象, 提出了微极性流体模型, 其与牛顿流体的主要差别在于放弃欧拉柯西应力原理的假设, 并考虑了流体运动中非零偶应力张量存在的情况. Papautsky^[22]基于微极性流体理论的数值模型描述了微通道流体行为. Ali和Hayat ^[23]研究了不可压缩微极性流体在不对称通道中的蠕动运动, 通过数值积分研究了相关参数对各波长压力上升的影响. Siddiqui和Lakhtakia ^[24]求解了在均匀外加电场作用下矩形微通道中微极性流体的稳态、对称和一维电渗流的边界值问题. Wang等^[25]采用有限差分法分析了圆形、圆柱形管中微极性流体的磁流体动力学流动. Misra等^[26]研究两个平行多孔板包围的微通道中微极性流体的电渗流, 说明了微极性参数的增大会导致微旋转幅度的减小. Ding等^[27]研究了不可压缩微极性流体在两个无限延伸的微平行板之间的时间周期性电渗流动, 给出了相关参数对速度和微旋转的影响. 之后Ding等^[28]又分析了流体颗粒微观结构对电动现象的影响, 结果表明, 微极性的引入会显著影响电动效应, 特别是在EDL重叠的情况下. Chaube等^[29]通过考虑在没有体力和耦合效应的情况下微极性流体的控制方程, 分析流体的微极性对通过平行板微通道的电动驱动蠕动泵送引起的完全流动的影响, 结果表明, 施加外部电场会改变蠕动流动. Huang和Huang ^[30]考虑了在一般滑移边界条件影响下狭缝微通道几何形状中微极性液体的电动扩散渗透流, 解得了微旋转、线速度和体积流量的解析解. Rana等^[31]利用纳米粒子的相关性和微极性流体模型, 建立了新的复耦合非线性微分方程集. Zhu^[32]基于微极性流体本构方程, 建立了Hele-Shaw流动的数学模型, 将微极性流体的运动方程扩展到多孔介质中的流动, 得到了流速、压力梯度和流速之间的关系. Karampour等^[33]研究了带疏水壁的微极性流体的电渗流, 并给出多个参数对流动模式的影响. Jaiswal和Yadav ^[34]考虑了微极性流体在充满多孔介质的旋转环空区域的Couette流动, 采用解析方法求解流经环形多孔区域的微极性流体的线速度、微旋转速度、剪切应力和耦合应力; 在自旋条件下, 与无自旋条件下的速度分布相比, 速度分布的值较低. Narla等^[35]分析了弯曲通道中热微极性流体电渗流中的传热, 得到各参数对速度及传热特性的影响. Fatunmbi等^[36]在具有多种滑移特性的指数可拉伸薄片的非定常流动中进行了磁微极性流体的熵生成分析, 微极性参数、Eckert和Prandtl数的增长导致系统中熵生成的增大. Rauf等^[37]研究了三维不可压缩和时间无关的微极性流体在两个平行圆盘之间的层流流动, 通过Runge-Kutta-Fehlberg方法求解非线性方程组, 得到速度和微旋转的变化趋势等.

综上所述, 尽管微通道中微极性流体电渗流已有大量研究, 但大多是基于壁面低Zeta势(即Zeta势远小于25 mV时)的情形, 利用Debye-Hückel(D-H)线性近似求解问题的控制方程. 在实际的工程应用中, 很多情况下Zeta势高于25 mV. 基于此, 本文研究高Zeta势下平行板微通道中微极性流体的时间周期电渗流问题. 利用有限差分法求解非线性Poisson-Boltzmann方程和不可压缩微极性流体的连续性方程、动量方程、角动量方程及本构方程, 并分析相关参数对流体流动速度及微旋转的影响.

2.   数学模型与公式推导

本文研究在高Zeta势下一类不可压缩微极性流体通过平行板微通道的时间周期电渗流. 在平行板微通道中, 假设通道长度和宽度远远大于其高度, 如图1所示, 在直角坐标系$(X, Y, Z)$下, 通道的高度为$2 H$($H = 100 {\text{ μm}}$), 该通道具有平行于X轴分布均匀强度为${{\boldsymbol{{E}}}_{{\text{app}}}}$的外加时变电场$({{\boldsymbol{E}}_{{\text{app}}}} = ({E_x}, {\text{ }}0, {\text{ }}0),$ ${E_x} = {E_0}\cos (\omega t))$, 在此外加电场作用下微极性流体沿X轴流动. 假设电荷在X轴方向分布均匀, 则电势仅是关于Y的函数.

2.1   电势分布

根据静电学理论, 单位体积的静电荷密度$\rho _{\text{e}}^{{*}}$与EDL的电势$\bar \psi$由下式表示^[26]:

$$\Delta \bar \psi = - \frac{{\rho _{\text{e}}^{\text{*}}}}{\varepsilon },$$

(1)

$$\rho _{\text{e}}^{\text{*}} = - 2{n_0}{z_{{\nu }}}e\sinh \left(\frac{{{z_{{\nu }}}{\text{e}}\bar \psi }}{{{k_{\text{b}}}T}}\right),$$

(2)

其中, $\varepsilon$为电解质液体的介电常数, ${n_0}$为电解质溶液的离子浓度, ${z_{{\nu }}}$为离子化合价, e为电子电荷, ${k_{\text{b}}}$为Boltzmann常数, T为绝对温度.

如果平行板的两个板是由不同材料制成, 那么上下Zeta势可能会有所不同. Zeta势可以写成如下形式^[26]:

$$\bar \psi = {\bar \psi _{\text{u}}},{\text{ }}Y = H,$$

(3)

$$\bar \psi = {\bar \psi _{1}},{\text{ }}Y = - H{.}$$

(4)

类似于文献[26], 电势$\bar \psi$不受外加电场时间波动的影响, 联立(1)式和(2)式, 并简化可以得到:

$$\frac{{{{\text{d}}^2}\bar \psi }}{{{\text{d}}{Y^2}}} = \frac{{2{n_0}{z_{{\nu }}}e}}{\varepsilon }\sinh \left(\frac{{{z_{{\nu }}}e\bar \psi }}{{{k_{\text{b}}}T}}\right).$$

(5)

定义以下无量纲量:

$$\left\{ \begin{aligned} &\psi = \frac{{{z_{{\nu }}}e\bar \psi }}{{{k_{\text{b}}}T}},~~{\text{ }}y = \frac{Y}{H},~~{\text{ }}m = \kappa H, \\ &\kappa = {\left(\frac{{2{n_0}z_{{\nu }}^{2}{e^2}}}{{\varepsilon {k_{\text{b}}}T}}\right)^{{1/2}}},~~{\text{ }}\beta = \frac{{{{\bar \psi }_{\text{u}}}}}{{{{\bar \psi }_1}}}, \\ &{\psi _u} = \frac{{{z_{{\nu }}}e{{\bar \psi }_u}}}{{{k_{\text{b}}}T}},~~{\text{ }}{\psi _1} = \frac{{{z_{{\nu }}}e{{\bar \psi }_1}}}{{{k_{\text{b}}}T}}, \end{aligned} \right.$$

(6)

其中$\kappa$称为Debye-Hückel参数, $1/\kappa$代表EDL的厚度, 称为Debye长度, $m$为无量纲的电动宽度, 表示微通道的半宽度与Debye长度的比值.

利用无量纲量, 方程(5)和边界条件(3)式和(4)式可化为

$$\frac{{{{\text{d}}^2}\psi }}{{{\text{d}}{y^2}}} = {m^2}\sinh (\psi ),$$

(7)

$$\psi = \beta {\psi _1},{\text{ }}y = 1,$$

(8)

$$\psi = {\psi _1},~~ {\text{ }}y = - 1.$$

(9)

单位体积的静电荷密度$\rho _{\text{e}}^{\text{*}}$可以化为

$$\rho _{\text{e}}^{\text{*}} = - 2{n_0}{z_{{\nu }}}e\sinh (\psi ) = - \frac{{\varepsilon {k_{\text{b}}}T}}{{{z_{{\nu }}}{\text{e}}{H^2}}}{m^2}\sinh (\psi ).$$

(10)

2.2   运动控制方程

微极性流体是具有微观结构的流体, 属于一类具有非对称应力张量的流体. 微极性流体表现出耦合应力, Eringen^[21]提出的微极性流体模型中流体的颗粒除了速度矢量外还具有独立的旋转矢量, 即以旋转张量${\nu _{ij}}$绕体积元的质心进行旋转. 考虑到偏对称条件, ${\nu _{ij}}$的独立个数为3, 因此引入一个新的向量${{\boldsymbol{\nu }}_k}$为^[26]

$${{\boldsymbol{\nu }}_k} = \frac{1}{2}{\varepsilon _{kij}}{\nu _{ij}},$$

其中${\varepsilon _{ijk}}$是交替向量, 轴矢量${{\boldsymbol{\nu }}_k}$是旋转矢量.

对于本文所考虑不可压缩微极性流体, 其连续性方程、动量方程和角动量方程表示为^[26]

$$\nabla \cdot {\boldsymbol{V}} = 0,$$

(11)

$$\rho \frac{{{\text{D}}{\boldsymbol{V}}}}{{{\text{D}}t}} = f + \nabla {\boldsymbol{T}},$$

(12)

$$\rho j\frac{{{\text{D}}{\boldsymbol{\nu}} }}{{{\text{D}}t}} = {\boldsymbol{l}} + \nabla {\boldsymbol{C}} + {{\boldsymbol{T}}^*},$$

(13)

其中${\boldsymbol{V}}$为流体速度矢量, $\rho$为密度, ${\boldsymbol{f }}$为单位体积的体力, T为应力张量, $j$为微惯性系数, ${\boldsymbol{\nu}}$为微旋转矢量, ${\boldsymbol{l}}$为单位体积的偶体力, ${\boldsymbol{C}}$为偶应力张量, ${\boldsymbol{T}}_k^*$为应力张量的矢量形式, 定义为

$${\boldsymbol{T}}_k^* = {\varepsilon _{kij}}{{\boldsymbol{T}}_{ij}}.$$

本文所考虑微极性流体是线性且各向同性的^[26], 其本构方程为

$$\begin{split} {\boldsymbol{T}} =\;& ( - p + \lambda \nabla {\boldsymbol{V}})I + \mu (\nabla {\boldsymbol{V}} + \nabla {{\boldsymbol{V}}^{\text{T}}})\\ &+ \chi (\nabla {\boldsymbol{V}} - \nabla {{\boldsymbol{V}}^{\text{T}}}) - 2\chi {\sigma ^*},\end{split}$$

(14)

$${\boldsymbol{C}} = {\alpha _\nu }(\nabla \cdot {\boldsymbol{\nu }})I + {\beta _\nu }(\nabla {\boldsymbol{\nu }}+ \nabla {{\boldsymbol{\nu}} ^{\text{T}}}) + \gamma (\nabla {\boldsymbol{\nu}} - \nabla {{\boldsymbol{\nu}} ^{\text{T}}}),$$

(15)

其中$\sigma _{ij}^* = {\nu _{ij}}$为旋转张量, $\mu$为牛顿剪切黏度系数, $\chi$是涡旋黏度系数, $\lambda$为二阶黏度系数, ${\alpha _\nu }$, ${\beta _\nu }$和$\gamma$为自旋梯度黏度系数^[26].

将本构方程(14)和方程(15)代入动量方程(12)和角动量方程(13), 消去偶体力$l$^[22], 得到:

$$\begin{split}&\rho \left(\frac{{\partial {\boldsymbol{V}}}}{{\partial t}} + ({\boldsymbol{V}} \cdot \nabla ){\boldsymbol{V}}\right)\\ =\;& - \nabla p + (\mu + \chi )\Delta {\boldsymbol{V}} + \chi \nabla \times {\boldsymbol{\nu }}+ {\boldsymbol{f}}, \end{split}$$

(16)

$$\begin{split}&~ \rho j\left(\frac{{\partial {\boldsymbol{\nu}} }}{{\partial t}} + ({\boldsymbol{V}} \cdot \nabla ){\boldsymbol{\nu }}\right) = \\ & \gamma \nabla {\boldsymbol{\nu}} {+} \chi \nabla \times {\boldsymbol{V}} {+} ({\alpha _\nu } + {\beta _\nu })\nabla (\nabla \cdot {\boldsymbol{\nu}} ) - 2\chi {\boldsymbol{\nu}} ,\end{split}$$

(17)

其中${\boldsymbol{f}} = \rho _{\text{e}}^{\text{*}}{{\boldsymbol{E}}_{{\text{app}}}}$为单位体积的电场力, $j$为微惯性系数, 定义为

$$j = \frac{{2\gamma }}{{2\mu + \chi }}.$$

壁面无滑移边界条件为

$${\left. {\boldsymbol{V}} \right|_{Y = \pm H}} = 0.$$

(18)

施加于微旋转的边界条件为

$${\boldsymbol{\nu}} + s\nabla \times{\boldsymbol{ V}} = 0,~~ {\text{ }}Y = \pm H,$$

(19)

其中, 常数$s \in [ - 1, 0]$. 为简单起见, 本文选择$s = 0$, 考虑没有压力驱动的情况即令(16)式中$p = 0$.

流体的速度和微旋转矢量表示为^[26]

$${\boldsymbol{ V}} = (U(Y,{\text{ }}t),{\text{ }}0,{\text{ }}0),{\text{ }}{\boldsymbol{\nu}} = (0,{\text{ }}0,{\text{ }}\varPhi (Y,{\text{ }}t)){.}$$

则连续性方程(11)自动满足, 方程(16)和方程(17)简化为

$$\rho \frac{{\partial U}}{{\partial t}} = (\mu + \chi )\frac{{{\partial ^2}U}}{{\partial {Y^2}}} + \chi \frac{{\partial \varPhi }}{{\partial Y}} + \rho _{\text{e}}^{\text{*}}{E_x},$$

(20)

$$\rho j\frac{{\partial \varPhi }}{{\partial t}} = \gamma \frac{{{\partial ^2}\varPhi }}{{\partial {Y^2}}} - \chi \frac{{\partial U}}{{\partial Y}} - 2\chi \varPhi ,$$

(21)

其中${{\boldsymbol{E}}_{{\text{app}}}} = {E_x}{{\boldsymbol{e}}_x}$.

假设周期性电渗流的外加电场、速度和微旋度表示为^[26]

$$\left\{ \begin{aligned} &{E_x} = {\text{Re}}({E_0}{{\text{e}}^{{\text{i}}\omega t}}), \\ &U = {\text{Re}}(U_0^*{{\text{e}}^{{\text{i}}\omega t}}), \\ &\varPhi = {\text{Re}}(\varPhi _0^*{{\text{e}}^{{\text{i}}\omega t}}), \end{aligned} \right.$$

(22)

其中${E_0}$, $U_0^*$和$\varPhi _0^*$分别为外加电场E的幅值、电渗透速度U和微旋转$\varPhi$的复幅值, $\omega$为外加电场的振荡角频率.

引入以下无量纲量:

$$\left\{ \begin{aligned} & y = {Y}/{H},~~{\text{ }}u = {U}/{V},~~{\text{ }}{u_0} = {{U_0^*}}/{V}, \\ &\phi = \frac{H}{V}\varPhi ,{\text{ }}{\phi _0} = \frac{H}{V}\varPhi _0^*,~~ {\rho _{\text{e}}} = \frac{{{z_{{\nu }}}e{H^2}}}{{\varepsilon {k_{\text{b}}}T}}\rho _{\text{e}}^{\text{*}},\\ & {\text{ }}\varOmega = \frac{{\omega \rho {H^2}}}{\mu }, ~~~ V = - \frac{{\varepsilon {k_{\text{b}}}T{E_0}}}{{\mu {z_{{\nu }}}e}}, \end{aligned} \right.$$

(23)

其中, $\varOmega$为施加电场频率$\omega$与流体系统特征频率${\omega ^*} = \mu /(\rho {H^2})$之比.

结合(22)式和(23)式, (20)式和(21)式的无量纲形式可写为

$$\frac{{{{\text{d}}^2}{u_0}}}{{{\text{d}}{y^2}}} + {k_1}\frac{{{\text{d}}{\phi _0}}}{{{\text{d}}y}} - {\text{i}}Re{u_0} + {m^2}(1 - {k_1})\sinh (\psi ) = 0,$$

(24)

$$\frac{{{{\text{d}}^2}{\phi _0}}}{{{\text{d}}{y^2}}} - {k_2}\frac{{{\text{d}}{u_0}}}{{{\text{d}}y}} - (2{k_2} + {\text{i}}R){\phi _0} = 0,$$

(25)

其中,

$$\begin{aligned} & {k_1} = \frac{\chi }{{\mu + \chi }},{\text{ }} ~~{k_2} = \frac{{\chi {H^2}}}{\gamma },~~ {\text{ }}Re = \frac{{\omega \rho {H^2}}}{{\mu + \chi }},{\text{ }}\\ &R = \frac{{\omega \rho j{H^2}}}{\gamma } = \frac{{2\omega \rho {H^2}}}{{2\mu + \chi }}, \end{aligned}$$

(26)

式中${k_1}$耦合了两个黏度系数, 表示微极性参数$(0 \leqslant {k_1} < 1)$, ${k_2}$为耦合应力参数^[28], R表示微旋转雷诺数, Re表示微极性流体的电振荡雷诺数.

方便起见, Re和R在无量纲频率$\varOmega$和微极性参数${k_1}$下重新定义为

$$Re = \varOmega (1 - {k_1}),{\text{ }} ~~~ R = \varOmega \frac{{2(1 - {k_1})}}{{2 - {k_1}}}.$$

(27)

无量纲边界条件为

$${u_0} = 0,{\text{ }}y = \pm 1,$$

(28)

$${\varphi _0} = 0,{\text{ }}y = \pm 1.$$

(29)

2.3   体积流量

为了研究相关无量纲参数对微极性流体体积流量的影响, 利用速度的表达式, 可以求出微通道单位宽度的无量纲体积流量$Q$为

$$Q = \int_{ - 1}^1 {u(y,t){\text{d}}y} = {{\mathrm{Re}}} ({Q_0}{{\text{e}}^{{\text{i}}\omega t}}),{\text{ }}{Q_0} = \int_{ - 1}^1 {{u_0}(y){\text{d}}y} .$$

(30)

2.4   微旋转强度

为了研究相关无量纲参数对微极性流体微旋转的影响, 利用微旋转的表达式, 可以求出微通道单位面积的无量纲微旋转强度$W$:

$${ W = \bigg|\displaystyle\int_0^1 {\phi (y,t){\text{d}}y} \bigg| = \big|{{{\mathrm{Re}}} ({W_0}{{\text{e}}^{{\text{i}}\omega t}})} \big|,\; {W_0} = \bigg|\displaystyle\int_0^1 {{\phi _0}(y){\text{d}}y} \bigg|.}$$

(31)

2.5   壁面剪切应力

不可压缩微极性流体的壁面剪切应力无量纲化后可写成如下形式:

$${\sigma _{12}} = {\left. {\frac{{\partial u}}{{\partial y}}} \right|_{y = \pm 1}},~~{\text{ }}{\sigma _{21}} = {\left. {\frac{1}{{1 - {k_1}}}\frac{{\partial u}}{{\partial y}}} \right|_{y = \pm 1}},$$

(32)

其中, ${\sigma _{12}}$和${\sigma _{21}}$表示$\sigma _{12}^*$和$\sigma _{21}^*$的无量纲形式:

$${\sigma _{12}} = \frac{H}{{\mu V}}\sigma _{12}^*,~~ {\text{ }}{\sigma _{21}} = \frac{H}{{\mu V}}\sigma _{21}^*.$$

3.   问题求解

本文求解平行板微通道中一类不可压缩微极性流体在高Zeta势下的时间周期电渗流, 利用有限差分法求解电势分布、速度分布、微旋转分布及剪切应力分布所满足的方程, 进而利用数值积分求解体积流量和微旋转强度, 并利用Matlab软件画出数值解图形, 分析各相关参数对流体速度、体积流量、微旋转等的影响. 本文中几个相关量均取以下数值: $\rho = 1.2 \times {10^3} \;{\text{kg}} {\cdot} {{\text{m}}^{{{ - 3}}}}$, $\mu = 3 \times {10^{ - 2}} \;{\text{kg}} {\cdot} {{\text{m}}^{{{ - 1}}}} {\cdot} {{\text{s}}^{{{ - 1}}}}$, $\gamma = 4.8 \times {10^{ - 10}}\; {\text{kg}} {\cdot} {\text{m}} {\cdot} {{\text{s}}^{{{ - 1}}}}$, 此时有${k_2} = \dfrac{\mu {k_1}{H^2}}{\gamma (1 - {k_1})}$.

3.1   电势分布

3.1.1   有限差分法

根据边界条件的无量纲化, $y$轴方向区间为$[-1,1]$, 对(7)式—(9)式进行中心差分得到如下格式:

$$({{{\psi _{i + 1}} - {\psi _i} + {\psi _{i - 1}}}})/{{{h^2}}} = {m^2}\sinh ({\psi _i}),$$

(33)

$${\psi _0} = {\psi _1}, {\text{ }}{\psi _N} = \beta {\psi _1},$$

(34)

其中, 步长$h = 1/N, \;i = 1—(N - 1)$.

3.1.2   结果比较

利用有限差分法和Matlab自带的BVP4c求解低Zeta势下电势分布所满足的方程(7)—(9), 将结果与使用D-H线性近似所得解析解对比. 如图2所示, BVP4c方法、有限差分法与解析解均有较好的吻合性. 因此, 考虑到有限差分法求解此类问题的有效性和便捷性, 本文将其推广用于求解高Zeta势的情形.

3.2   速度及微旋转分布

3.2.1   有限差分法

根据边界条件的无量纲化, $y$轴方向区间为$[-1,1]$, 对方程组(24)和(25)进行中心差分和向后差分得到如下格式(后续${u_0}$, ${\phi _0}$下标均省略):

$$\begin{split}& \frac{{{u_{j + 1}} - 2{u_j} + {u_{j - 1}}}}{{{h^2}}} + {k_1}\frac{{{\phi _j} - {\phi _{j - 1}}}}{h} - iRe{u_j} \\ =\;& - (1 - {k_1}){m^2}{\psi _j}, \end{split}$$

(35)

$$\frac{{{\phi _{j + 1}} - 2{\phi _j} + {\phi _{j - 1}}}}{{{h^2}}} - {k_2}\frac{{{u_j} - {u_{j - 1}}}}{h} - (2{k_2} + iR){\phi _j} = 0,$$

(36)

其中, 步长$h = 1/N,\; i = 1—(N - 1)$.

相应的边界条件由(28)式、(29)式经过离散得到:

$${u_0} = {u_N} = 0,$$

(37)

$${\phi _0} = {\phi _N} = 0.$$

(38)

3.2.2   结果比较

利用有限差分法求解低Zeta势下速度及微旋转分布所满足的方程(24)和(25)及边界条件(28)式和(29)式, 将得到的结果与使用D-H线性近似所得解析解对比. 如图3和图4所示, 结果非常吻合. 因此, 将有限差分法推广到速度及微旋转分布的求解是可行的, 由此可以将有限差分法推广到求解高Zeta势下的情形.

3.3   体积流量

根据边界条件的无量纲化, $y$轴方向区间为[–1, 1], 将(30)式得到的Q在$y$方向划分为${N_1}$个小区间, 对(30)式利用复化梯形法即可得体积流量近似值为

$${Q_0} = {h_1}[{u_{0,0}} + 2\sum\limits_{i = 1}^{{N_1} - 1} {{u_{0,i}} + {u_{0,{N_1}}}} ]/2,$$

(39)

其中, $h = 1/N, \;i = 1—(N - 1)$.

3.4   微旋转强度

根据边界条件的无量纲化, $y$轴方向区间为$[-1,1]$, 将方程(31)得到的$W$在$y$方向划分为${N_1}$个小区间. 对(31)式利用复化梯形法即可得微旋转强度近似值为

$${{W_0} = \bigg| {{h_1} \bigg[{\phi _{0,{N_1}/2}} + 2\sum\limits_{i = {N_1}/2 + 1}^{{N_1} - 1} {{\phi _{0,i}} + {\phi _{0,{N_1}}}} \bigg]\Big/2} \bigg|,}$$

(40)

其中${h_1} = 1/{N_1}, {\text{ }}i = {N_1}/2 + 1—({N_1} - 1)$.

3.5   壁面剪切应力

将方程(32)中关于速度的导数分别用差商近似代替, 即可得壁面$y = - 1$处的剪切应力${\sigma _{12}}$和${\sigma _{21}}$的近似值分别为

$${\sigma _{12}} = \frac{{{u_{0,1}} - {u_{0,0}}}}{{{h_1}}},{\text{ }}{\sigma _{21}} = \frac{{{u_{0,1}} - {u_{0,0}}}}{{{h_1}(1 - {k_1})}},$$

(41)

其中, ${h_1} = 1/{N_1}$.

4.   结果与讨论

本文利用有限差分法研究高Zeta势下平行板微通道中一类不可压缩微极性流体的时间周期电渗流. 下面讨论不可压缩微极性流体在高Zeta势情形下电动宽度m、无量纲电振荡频率$\varOmega$、微极性参数${k_1}$对流体速度$u$、微旋转$\phi$、体积流量$Q$、微旋转强度$W$以及剪切应力${\sigma _{12}}$和${\sigma _{21}}$等物理量的影响. 典型的参数取值范围为$0 \leqslant {k_1} < 1$和$5 \leqslant m \leqslant 100$(以下均用${u_0}$, ${\phi _0}$, ${Q_0}$, ${W_0}$表示复函数$u$, $\phi$, $Q$, $W$的幅值).

图5分析了不同Zeta势和微极性参数${k_1}$对无量纲速度的影响. 从图5可以看出, EOF速度的振幅随着Zeta势的增大而增大, 这是由于双电层内的电荷密度会随着Zeta势的增大而增大, 导致电渗力增大, 进而加快电渗流速度(图5(a)), 说明与低Zeta势相比, 高Zeta势对微极性流体的速度有显著的促进作用. 微极性参数${k_1}$的增大使速度减小, 这是因为${k_1}$是涡旋黏度系数$\chi$与涡旋黏度系数$\chi$和牛顿剪切黏度系数$\mu$之和的比值, 因此${k_1}$增大说明$\mu$减小或者$\chi$增大, 意味着涡旋黏度系数发挥越来越重要的作用, 此时越来越多的流体动量去支配流体粒子的微旋转, 因此速度就会减小并趋于零(图5(b)).

图6分析了高Zeta势下, 对于不同的$m$和$\varOmega$, 无量纲速度${u_0}$的变化趋势. 从图6(a), (b)可以看出, 当电振荡频率$\varOmega$较低时, 随着电动宽度$m$的增大, 电渗流速度分布趋于稳定的速度加快, 且呈现塞子型, 这是因为在远离壁面处, 流体溶液中正负离子数几乎相等, 导致电势趋于零, 所以速度趋于稳定; 而当电振荡频率$\varOmega$较高时, 速度分布主要集中在双电层内, 双电层外的速度明显减小, 特别是位于通道中心的流体几乎不流动(图6(c), (d)).

图7分析了不同Zeta势和微极性参数${k_1}$对无量纲微旋转的影响. 从图7(a)可以看出, 微旋转的振幅随着Zeta势的增大而增大, 微通道壁面和中心处没有发生微旋转, 双电层内微旋转较为明显, 说明与低Zeta势相比, 高Zeta势对微极性流体微旋转有显著的促进作用. 从图7(b)可以看出, 随着微极性参数${k_1}$的增大, 微旋转幅度是增大的; 但是当它增大到${k_1} \to 1$时, 微旋转幅度又开始减小, 这是由于此时涡旋黏度系数$\chi$远远大于牛顿剪切黏度系数$\mu$, 由其引起的能量耗散占主导地位, 因此微旋转效应变弱.

图8描述了高Zeta势下, 对于不同的m和$\varOmega$, 无量纲微旋转幅值${\phi _0}$的变化趋势. 由图8可以看出, 随着$\varOmega$的增大微旋转幅值${\phi _0}$减小; 电振荡频率$\varOmega$较低时, m的值越小, 微旋转幅值${\phi _0}$越大(图8(a), (b)); 当$\varOmega$较高时, 变化趋势则相反(图8(c), (d)). 此外, 随着$\varOmega$的增大, 微旋转幅值${\phi _0}$的振荡特性变得明显, 这是因为, 从图6(a)可以看出当电振荡频率$\varOmega$较低、电动宽度m较大时, 除壁面附近微通道中其他区域的速度几乎没有变化, 微旋转主要集中在EDL区域内(图8(a), (b)), 而图6(d)显示了较高电振荡频率下微通道中心的速度平缓, 此时微旋转振幅减小, 且远离壁面处的微旋转逐渐趋于零(图8(c), (d)). 图9分析了不同Zeta势下体积流量${Q_0}$的变化趋势. 可以看出电势越大, 体积流量越大, 原因是Zeta势增大, 速度加快, 使得单位宽度的体积流量增多; 说明与低Zeta势相比, 高Zeta势对微极性流体体积流量有显著的促进作用.

图10显示了高Zeta势下体积流量${Q_0}$随$m$和${k_1}$的变化趋势. 从图10(a)可以看出, 对于给定的微极性参数${k_1}$和较小的电振荡频率$\varOmega$, 当$m$增大时, ${Q_0}$的数值增大达到最大值后有略微的减小, ${k_1}$只是会改变最大值, 不会改变这一趋势; 这是因为在较小的电振荡频率$\varOmega$下, 速度呈现塞子型. 从图10(a)还可以看出, 对于给定的电动宽度$m$和较小的电振荡频率$\varOmega$, 当${k_1}$增大时, ${Q_0}$减小并趋于零, 这是由于速度随着${k_1}$的增大而减小并趋于零(图5(b)). 从图10(b)可以看出, 对于给定的微极性参数${k_1}$和电动宽度$m$, 电振荡频率$\varOmega$有一个临界值${\varOmega _0} = 1$, 当$\varOmega < 1$时体积流量${Q_0}$不变, 当$\varOmega > 1$时体积流量${Q_0}$减小.

图11分析了不同Zeta势下微旋强度${W_0}$的变化趋势. 从图11可以看出, 微旋强度幅值与Zeta势是正相关的, 随着Zeta势的增大, 微旋强度整体是增大的, 而且关于${k_1}$和$\varOmega$的趋势不变, 说明与低Zeta势相比, 高Zeta势对微极性流体微旋强度有显著的促进作用. 这是因为, 随着Zeta势的增大, 微旋转幅度增大(图7).

图12描述了高Zeta势下微旋强度${W_0}$在不同${k_1}$和$\varOmega$下的变化趋势. 从图12(a)可以看出, 当微极性参数${k_1} \to 0$或${k_1} \to 1$时, ${W_0}$的数值趋于零; 这是因为, 当微极性参数${k_1} \to 0$时, 流体几乎不发生微旋转, 且当微极性参数增大到${k_1} \to 1$时, 微旋转幅值又减小. 定义微旋的特征穿透深度^[22]${\delta _m} \approx O(\sqrt {{\chi \mathord{\left/ {\vphantom {\chi {\rho \omega }}} \right. } {\rho \omega }}} )$, 由于${k_1}$是关于$\chi$的单调递增函数, ${k_1}$从0增大到1, 意味着$\chi$是增大的, 穿透深度随之增大, 因此微旋强度增大. 微旋强度不仅与穿透深度有关, 还取决于微旋转振幅的定义. 从图12(b)可以看出, 对于给定的微极性参数${k_1}$, 当电振荡频率$\varOmega$很低($\varOmega < 1$)时, ${W_0}$是一个固定常数值, 而且电动宽度$m$的值越小${W_0}$的值越大. 对于较大的$m$, 在$\varOmega$达到中等数值处${W_0}$有最大值. 随着电振荡频率$\varOmega$的进一步增大微旋强度振幅趋于零, 且随着电动宽度$m$的增大而增大; 这是由于在频率非常高的情况下, 扩散时间尺度远大于振荡时间周期. 因此, 流动角动量没有足够的时间从壁面双电层向主体区域扩散很远, 穿透深度较小, 所以微旋强度减弱.

图13分析了不同Zeta势下剪切应力${\sigma _{21}}$和${\sigma _{12}}$的变化趋势. 由图13可以看出, 随着Zeta势的增大, 两个剪切应力均增大. 这是因为随着Zeta势的增大, 壁面附近流体速度是增大的, 从而由(27)式可知剪切应力${\sigma _{21}}$和${\sigma _{12}}$增大. 说明与低Zeta势相比, 高Zeta势对微极性流体的剪切应力有显著的促进作用.

图14显示了高Zeta势下在壁面$y = - 1$处剪切应力${\sigma _{21}}$在不同$m$, ${k_1}$和$\varOmega$下的变化趋势. 从图14可以看出, 随着电动宽度$m$的增大, 剪切应力${\sigma _{21}}$增大; 这是因为电动宽度$m$增大, 速度增大. 当电振荡频率$\varOmega$较小($\varOmega < 1$)时, 微极性参数${k_1}$的取值不影响壁面剪切应力的幅值, 这是由于较低的电振荡频率$\varOmega$下, 微极性参数${k_1}$几乎不影响速度的变化率; 当电振荡频率$\varOmega$较大($\varOmega > 1$)时, 速度随着微极性参数${k_1}$的增大而减小, 因此壁面剪切应力${\sigma _{21}}$随着微极性参数${k_1}$的增大而减小. 当电振荡频率$\varOmega$较小($\varOmega < 1$)时, 速度不随$\varOmega$的变化而变化, 因此剪切应力${\sigma _{21}}$不随$\varOmega$的变化而变化; 当电振荡频率$\varOmega$较大($\varOmega > 1$)时, 速度随着电振荡频率$\varOmega$的增大而减小, 因此壁面剪切应力${\sigma _{21}}$随着电振荡频率$\varOmega$的增大而减小.

图15分析了壁面剪切应力${\sigma _{12}}$在不同$m$, ${k_1}$和$\varOmega$下的变化趋势. 从图15(a)可以看出, 对于给定的较小的电振荡频率$\varOmega$($\varOmega < 1$), ${\sigma _{12}}$随着微极性参数${k_1}$的增大而线性减小, 且当${k_1} \to 1$时, ${\sigma _{12}}$的振幅趋于零; ${\sigma _{12}}$随着电动宽度m的增大而增大; 这是因为, 当电振荡频率$\varOmega$较小($\varOmega < 1$)时, 壁面附近速度随着微极性参数${k_1}$的增大而快速减小, 而随着电动宽度$m$的增大而增大. 从图15(b)可以看出, 对于给定的电动宽度m, 当电振荡频率$\varOmega$较小($\varOmega < 1$)时, ${\sigma _{12}}$不随$\varOmega$的增大而变化; 当电振荡频率$\varOmega$较大($\varOmega > 1$)时, ${\sigma _{12}}$随$\varOmega$的增大而减小. 这是由于当电振荡频率$\varOmega$较小($\varOmega < 1$)时, 对于给定的微极性参数${k_1}$, 速度斜率不变, 于是壁面剪切应力${\sigma _{12}}$处于恒定值; 当电振荡频率$\varOmega$较大($\varOmega > 1$)时, 速度幅值减小, 因而壁面剪切应力${\sigma _{12}}$随$\varOmega$的增大而减小.

5.   结　论

本文研究了平行板微通道中一类不可压缩微极性流体在高Zeta势下的时间周期电渗流, 利用有限差分法求解了高Zeta势下电势所满足的非线性Poisson-Boltzmann方程以及不可压缩微极性流体的连续性方程、动量方程、角动量方程和本构方程, 分析了高Zeta势下电动宽度m、电振荡频率$\varOmega$及微极性参数${k_1}$等无量纲参数对不可压缩微极性流体的速度和微旋转效应的影响. 结果表明: 1)随着Zeta势的增大, 微极性流体的速度、微旋转、体积流量、微旋强度以及剪切应力增大, 说明与低Zeta势相比, 高Zeta势对微极性流体电渗流有显著的促进作用. 2)在高Zeta势下, 随着微极性参数${k_1}$的增大, 速度减小, 微旋转增强, 相应的体积流量减小, 微旋强度增大; 但是当微极性参数增大到${k_1} \to 1$时, 微旋转效应变弱, 微旋强度大大降低. 说明在高Zeta势下, 随着微极性参数的增大, 微极性流体的速度减小, 但是对微旋转效应的影响呈现先增强后减弱的趋势. 3)在高Zeta势下, 当电振荡频率$\varOmega$较低(小于1)时, 随着电动宽度m的增加, 电渗流速度分布趋于稳定的速度越来越快, 影响了速度的变化率, 而微旋转幅值减小; 当电振荡频率$\varOmega$较高(大于1)时, 随着电动宽度m的增大, 速度幅值及微旋转幅值减小, 微旋强度减弱, 但体积流量快速增大并趋于恒定. 4)在高Zeta势下, 随着电振荡频率$\varOmega$的增大, 速度和微旋转的振荡变得明显, 具体表现为: 在低电振荡频率$\varOmega$情况下速度和微旋转的峰值、体积流量及微旋强度均不变; 在高电振荡频率$\varOmega$情况下速度和微旋转的幅值减小, 体积流量及微旋强度减小直至趋于零, 说明电振荡频率对流体的流动特性有阻碍作用. 5)在高Zeta势下, 壁面剪切应力${\sigma _{21}}$及${\sigma _{12}}$的幅值随电动宽度m的增大而增大; 当电振荡频率$\varOmega$较小($\varOmega < 1$)时, 剪切应力${\sigma _{21}}$与${\sigma _{12}}$均取恒定值, 且微极性参数${k_1}$的取值不影响壁面剪切应力${\sigma _{21}}$的幅值; 当电振荡频率$\varOmega$较大($\varOmega > 1$)时, 壁面剪切应力${\sigma _{21}}$, ${\sigma _{12}}$的幅值随着电振荡频率$\varOmega$的增大而减小, 且壁面剪切应力${\sigma _{21}}$的幅值随着微极性参数${k_1}$的增大而减小, 而壁面剪切应力${\sigma _{12}}$的振幅随着微极性参数${k_1}$的增大而线性减小.