带电波纹壁纳米通道中纳米流体的流向势与能量转化效率

The Streaming Potential and Energy Conversion Efficiency of Nanofluid in Charged Corrugated Wall Nanochannels

摘要: 本文通过理论建模与数值模拟，系统研究了正弦波纹壁面微通道中纳米流体的流动特性、流向势及能量转化效率。建立了综合考虑波纹几何参数(波幅δ、波数n、相位差θ)与纳米颗粒属性(体积分数φ)的电动力学流动模型，采用摄动法推导了电势、速度、流向势及能量转化效率的解析解，并利用切比雪夫谱方法获取了数值解并进行验证。研究结果表明：纳米流体的引入虽因有效粘度增加而降低了整体流速，但显著抑制了波谷区域的回流现象，增强了流动稳定性；当上下壁面波纹相位差θ = π时，通道横截面积均匀，中心处流速最大；对称管道中心流速、流向势及能量转化效率均存在极值，在临界波数nπ ≈ 0.1π和较大波幅时达到峰值，且当纳米粒子浓度增加时，能量转化效率提高。本研究揭示了波纹壁面与纳米颗粒的协同作用机制，为优化微流控系统中机械能–电能转化器件的设计提供了理论依据。

Abstract: This study systematically investigates the flow characteristics, flow direction potential, and energy conversion efficiency of nanofiuid in microchannels with sinusoidal corrugated walls through theoretical modeling and numerical simulation. An electrodynamics model was established, comprehensively considering the corrugation geometric parameters (amplitude δ, number of waves n, phase difference θ) and nanoparticle properties (volume fraction φ). The analytical solutions for potential, velocity, flow direction potential, and energy conversion efficiency were derived using the small parameter perturbation method, with numerical solutions validated via Chebyshev spectral analysis. Results demonstrate that while nanofluid introduction reduces overall flow velocity due to increased effective viscosity, it significantly suppresses backflow in the trough regions, enhancing flow stability. When the phase difference θ between upper and lower corrugations equals π, the channel cross-section becomes uniform with maximum velocity at the center. Symmetrical channels exhibit extremum values for central velocity, flow direction potential, and energy conversion efficiency, peaking at critical wave numbers nπ ≈ 0.1π and large amplitudes. Moreover, the energy conversion efficiency increases with increasing nanoparticle concentration. This research elucidates the synergistic mechanism between corrugated walls and nanoparticles, providing theoretical foundations for optimizing mechanical-to-electrical energy conversion devices in microfluidic systems.

1. 引言

近年来，微流体装置在微机电系统[1]、生物和化学传感器[2]、药物输送芯片[3]等方面有很重要的应用。微通道中的流体运动输送特性得到了很好的研究，并且采用了几种常见的机制来驱动流体流动，包括经典的压力梯度[4]、外加电场[5]和电磁场[6]。在决定微尺度流体行为的各种现象中，电渗流(EOF)是最主要的传输机制。通常，当固体表面与电解质溶液接触时，由于界面的化学相互作用，微通道的璧面上会产生表面电荷。这些表面电荷会影响固体表面附近电解质溶液中阳离子和阴离子的分布，从而产生双电层(EDL) [7]。当施加带电表面切线方向的电场时，EDL中的离子在电场力作用下发生定向移动。此时可移动的离子沿电场方向运动，导致经典的电渗流(EOF) [8]。

另一个重要的电动力学现象是在压力驱动的流体中产生流向势，这是一种非电驱动的电动力学现象。当通道两端施加较大的压力差时，双电层内离子的流动会形成与流动方向相同的流向电流。同时，反离子(与壁面极性相反)会聚集到通道的下游，使得下游的电势比上游的高，形成与原流动方向相反的电场，称为流向势[9]-[11]。M.B. Andersen1等人[12]的实验表明流向势与表面电荷，电解质溶液浓度和通道尺寸都有关系。Xie等人[13]分析了纳米通道中两种不混溶流体的流向势和电动能量转化效率。理论结果表明，流向势可以看作是估计电动能量转化效率的标准。

一般来说，以上所提到的研究都是在光滑壁面的假设下进行的。但是，在实践中是不可能得到完全光滑的壁面。由于微通道的尺寸通常很小，因此壁面的微小变化就能对流体流动产生显著影响。杨大勇等人[14]针对具有正弦表面粗糙度的平行板微通道，采用有限元方法求解控制方程，研究了粗糙度单元高度和密度对EOF行为的影响。结果为微通道表面的优化设计提供了有价值的见解，以实现微芯片中的精确流量控制。长龙等人[15]采用摄动展开法，研究了具有正弦波纹的环形微管道中脉冲流动。Xie等人[16]研究了带电对称波纹表面的纳米通道中，电动流动变量的扰动和数值解，并系统地研究了波的振幅和波数等波纹几何形状对电动流动特性的影响。

在优化流动工质方面，研究表明，微通道流动工质中引入纳米颗粒也可提高流体性能。将纳米流体与微通道相组合，可更大程度提升流场性能。张莹等人[17]探讨了多孔介质中偶应力纳米流体的电动能量转化，揭示了纳米颗粒体积分数和电场参数对转化效率的影响。这些工作表明，纳米流体与微通道结构的结合是提升流体流动性能的有效策略。

目前，关于波纹壁效应研究的文献主要集中在牛顿流体领域。然而，微管中具有波纹壁表面的非牛顿流体的流向势现象尚未得到足够的关注。特别是在理论分析和数值模拟这两个关键方面。本文利用摄动法，从理论上研究了波动壁面微通道中纳米流体的流向势和电动能量转化。首先基于Debye-Hückel线性化假设，得到了电势、流速和压降分布，并进一步给出了流向势和电动能量转化效率的解析解，之后利用切比雪夫谱方法获得了数值解并进行比较。研究结果对具有高转化效率的微纳米流体器件的设计具有重要价值。

2. 数学模型

Figure 1. The schematic diagram of electrodynamic flow in the nano-channel of corrugated surface

图1. 波纹表面纳米通道中电动力学流动示意图

我们考虑一个具有波纹表面的二维(2D)微通道，如图1所示，建立以通道入口中心为原点的笛卡尔坐标系$\left(x^{*},y^{*}\right)$ ，其中$x^{*}$ 轴与流动方向一致，$y^{*}$ 轴垂直于通道壁向上。微通道的上壁面位于$y^{*}=h^{*}\left(x^{*}\right)=H\left(1+\delta \sin \left(\frac{n\pi x^{*}}{L}\right)\right)$ ，下壁面位于$y^{*}=m^{*}\left(x^{*}\right)=-H\left(1+\delta \sin \left(\frac{n\pi x^{*}}{L}+\theta \right)\right)$ 。其中H表示平均半通道高度，δ为无量纲波幅，L为特征长度，$k=n\pi /L$ 为无量纲波数，$\lambda =1/k$ 为波长，θ为上下壁面之间的相位差。微通道的几何形态主要由这些波纹参数决定。本研究中，纳米通道长度设定为6 L，并假设其长度远大于平均高度，即$6L\gg 2H$ 。流体域是由电解质溶液和Al₂O₃纳米颗粒组成的稳定纳米流体悬浮液。在压力驱动条件下，假设波纹状纳米通道壁带电且具有均匀的ζ电位，电解质溶液沿带电表面流动时，由电流输送的离子会在上游和下游之间产生电势差，此时出现了与流动方向相反的诱导电场，称为流向势。在纳米通道两端连接外部负载，可将压力产生的机械能转化为电能。

3. 控制方程

3.1. 电势场

根据静电学理论[7]，双电层(EDL)电势分布满足一维泊松方程：

$\frac{{\text{d}}^2{\psi }^{*}}{\text{d}{y}^{*}{}^2}=-\frac{{\rho }_e}{{\epsilon }_{nf}}$ , (1)

其中，${\rho }_e$ 是净电荷密度，${\epsilon }_{nf}$ 是纳米流体的介电常数，考虑纳米流体中纳米颗粒分布均匀，纳米颗粒与基液之间无相对滑移速度，其数值取决于纳米颗粒${\epsilon }_s$ 的介电常数和基础流体${\epsilon }_f$ 的介电常数。二者之间的关系为：

$\frac{{\epsilon }_{nf}}{{\epsilon }_f}=1+\frac{3\left({\epsilon }_s-{\epsilon }_f\right)\phi }{\left({\epsilon }_s+2{\epsilon }_f\right)-\left({\epsilon }_s-{\epsilon }_f\right)\phi }$ , (2)

其中$\phi$ 表示电解质溶液中纳米颗粒的体积分数，可以看到，在没有纳米颗粒的情况下(即$\phi =0$ )，纳米流体的介电常数与牛顿流体完全一致。此外，我们假设在对称电解质溶液中，净电荷密度的分布服从玻尔兹曼分布。

${\rho }_e=-2n_{\infty }{e}_0{z}_v\sinh \left(\frac{e_0{z}_v{\psi }^{*}}{k_B{T}_a}\right)$ , (3)

其中$n_{\infty }$ 为离子浓度，$e_0$ 为基本电荷，$z_v$ 为离子化合价，$k_B$ 为玻尔兹曼常数，$T_a$ 为绝对温度。在本研究中，我们关注低zeta电位(通常小于25 mV)并采用Debye-Hückel线性化近似，即$\text{sinh}\left(\frac{e_0{z}_v{\psi }^{*}}{k_B{T}_a}\right)\approx \frac{e_0{z}_v{\psi }^{*}}{k_B{T}_a}$ [7] [18]。将方程(3)代入方程(1)后，我们得到：

$\frac{{\text{d}}^2{\psi }^{*}}{\text{d}{y}^{*2}}=\frac{{\epsilon }_f}{{\epsilon }_{nf}}{\kappa }^2{\psi }^{*}$ , (4)

其中Debye-Hückel参数为$\kappa ={\lambda }^{-1}={\left\{\left(2n_{\infty }{e}_0^2{z}_v^2\right)/\left({\epsilon }_f{k}_B{T}_a\right)\right\}}^{1/2}$ ，λ是双电层厚度。方程(4)在上壁面处的边界条件为[19]：

${\psi }^{*}|{}_{y^{*}=h^{*}\left(x^{*}\right)}={\psi }_w$ , (5)

在下壁面处：

${\psi }^{*}|{}_{y^{*}=m^{*}\left(x^{*}\right)}={\psi }_w$ . (6)

为了得到泊松–玻尔兹曼方程及其边界条件的无量纲形式，引入以下无量纲变量与参数

$x=\frac{x^{*}}{L},y=\frac{y^{*}}{H},\psi =\frac{{\psi }^{*}}{{\psi }_w^{*}},{\psi }_w^{*}=-\frac{k_B{T}_a}{e_0{z}_v},K=\kappa H,{\epsilon }_r=\frac{{\epsilon }_{nf}}{{\epsilon }_f}$ . (7)

将上述无量纲变量和参数代入到方程(4)、(5)、(6)，我们得到如下方程和边界条件：

$\frac{{\text{d}}^2\psi }{\text{d}{y}^2}=\frac{K^2}{{\epsilon }_r}\psi$ , (8)

$\psi |{}_{y=h\left(x\right)}=1$ , (9)

$\psi |{}_{y=m\left(x\right)}=1$ . (10)

为求解方程(8)，采用正则摄动法，将$\delta$ 作为摄动参数(即$\delta \ll 1$ )，展开$\psi$ ：

$\psi ={\psi }_0+\delta {\psi }_1+O\left({\delta }^2\right)$ . (11)

将方程(11)代入方程(8)，我们得到零阶方程和一阶方程：

$\frac{{\text{d}}^2{\psi }_0}{\text{d}{y}^2}=\frac{K^2}{{\epsilon }_r}{\psi }_0,\frac{{\text{d}}^2{\psi }_1}{\text{d}{y}^2}=\frac{K^2}{{\epsilon }_r}{\psi }_1$ . (12)

此外，为精确获取通道壁的边界条件，我们采用泰勒级数展开法将边界条件表示为渐近等价的形式。

通过这种方法，我们可推导出边界条件在$\delta \ll 1$ 时的简化表达式，在上壁面($y=h\left(x\right)$ )和下壁面($y=m\left(x\right)$ )处，函数ψ在$y=1$ 和$y=-1$ 处的泰勒展开式为：

${\psi }_i|{}_{y=h}={\psi }_i|{}_{y=1}+\frac{\partial {\psi }_i}{\partial y}|{}_{y=1}\delta \sin \left(n\pi x\right)+O\left({\delta }^2\right)$ , (13)

${\psi }_i|{}_{y=m}={\psi }_i|{}_{y=-1}-\frac{\partial {\psi }_i}{\partial y}|{}_{y=-1}\delta \sin \left(n\pi x\right)+O\left({\delta }^2\right)$ . (14)

其中$i=0,1$ ，将方程(13)和方程(14)代入方程(11)，结合上述边界条件，我们可以得到：

${\psi }_0|{}_{y=\pm 1}=1,{\psi }_1|{}_{y=1}=-{\frac{\partial {\psi }_0}{\partial y}|}_{y=1}\sin \left(n\pi x\right),{\psi }_1|{}_{y=-1}={\frac{\partial {\psi }_0}{\partial y}|}_{y=-1}\sin \left(n\pi x+\theta \right)$ . (15a~c)

根据方程(12)和方程(15a-c)，我们推导出零阶电势和一阶电势的解析解：

${\psi }_0=\cosh \left(\frac{Ky}{\sqrt{{\epsilon }_r}}\right)\mathrm{sech}\left(\frac{K}{\sqrt{{\epsilon }_r}}\right)$ , (16)

${\psi }_1=\frac{K\left({\text{e}}^{\frac{K-Ky}{\sqrt{{\epsilon }_r}}}-\left({\text{e}}^{\frac{2K}{\sqrt{{\epsilon }_r}}}-{\text{e}}^{\frac{2Ky}{\sqrt{{\epsilon }_r}}}\right)\sin \left(\theta +n\pi x\right)-\left({\text{e}}^{\frac{2K(y+1)}{\sqrt{{\epsilon }_r}}}-1\right)\sin \left(n\pi x\right)\right)}{\sqrt{{\epsilon }_r}{\left({\text{e}}^{\frac{2K}{\sqrt{{\epsilon }_r}}}+1\right)}^2}$ . (17)

当${\epsilon }_r=0$ 且$\theta =0$ 时，其结果与我们之前研究中的波动壁面微通道中无相位差的牛顿流体一致[16]。

3.2. 速度分布

我们假设微通道中的流动是稳态的，压力驱动下的流体流动速度满足连续性方程和动量方程[20]：

$\frac{\partial u^{*}}{\partial x^{*}}+\frac{\partial v^{*}}{\partial y^{*}}=0$ , (18)

${\rho }_{nf}\left(u^{*}\frac{\partial u^{*}}{\partial x^{*}}+v^{*}\frac{\partial u^{*}}{\partial y^{*}}\right)=-\frac{\partial p^{*}}{\partial x^{*}}+\frac{\partial {\tau }_{xx}^{*}}{\partial x^{*}}+\frac{\partial {\tau }_{xy}^{*}}{\partial y^{*}}+{\rho }_e{E}_s^{*}$ , (19)

${\rho }_{nf}\left(u^{*}\frac{\partial v^{*}}{\partial x^{*}}+v^{*}\frac{\partial v^{*}}{\partial y^{*}}\right)=-\frac{\partial p^{*}}{\partial y^{*}}+\frac{\partial {\tau }_{yx}^{*}}{\partial x^{*}}+\frac{\partial {\tau }_{yy}^{*}}{\partial y^{*}}$ , (20)

${\tau }_{xy}^{*}={\mu }_{nf}\left(\frac{\partial u^{*}}{\partial y^{*}}+\frac{\partial v^{*}}{\partial x^{*}}\right)$ . (21)

其中，$u^{*}$ ，$v^{*}$ 分别表示流体在$x^{*}$ ，$y^{*}$ 方向上的速度，${\rho }_{nf}$ 为纳米流体密度，${\rho }_{nf}=\left(1-\phi \right){\rho }_f+\phi {\rho }_p$ [21]。$p^{*}$ 为压力，${\tau }^{*}$ 为剪切应力，$E_s^{*}$ 为流向势，${\mu }_{nf}$ 为纳米流体粘度。这里${\mu }_{nf}=\frac{1}{{\left(1-\phi \right)}^{2.5}}{\mu }_f$ [22]，需要指出的是，流向势仍是一个未知量。引入以下变量与参数，将流动方程与边界条件转化为无量纲形式：

$x=\frac{x^{*}}{L},y=\frac{y^{*}}{H},u=\frac{u^{*}}{u_{HS}},v=\frac{v^{*}L}{u_{HS}H},p=\frac{p^{*}{H}^2}{{\mu }_f{u}_{HS}L},\tau =\frac{{\tau }^{*}}{{\mu }_f{u}_{HS}/H}$ ,

$u_{HS}=-\frac{{\epsilon }_f{k}_B{T}_a{E}_0}{{\mu }_f{e}_0{z}_v},\mathrm{Re}=\frac{{\rho }_f{u}_{HS}H}{{\mu }_f},\beta =\frac{H}{L},{\mu }_r=\frac{{\mu }_{nf}}{{\mu }_f},{\rho }_r=\frac{{\rho }_{nf}}{{\rho }_f}$ (22)

方程(18)~(21)可以表示为：

$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$ , (23)

$\beta \mathrm{Re}{\rho }_r\left(u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}\right)=-\frac{\partial p}{\partial x}+\beta \frac{\partial {\tau }_{xx}}{\partial x}+\frac{\partial {\tau }_{xy}}{\partial y}+K^2\psi E_s$ , (24)

${\beta }^3\mathrm{Re}{\rho }_r\left(u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}\right)=-\frac{\partial p}{\partial y}+{\beta }^2\frac{\partial {\tau }_{yx}}{\partial x}+\beta \frac{\partial {\tau }_{yy}}{\partial y}$ , (25)

${\tau }_{xy}={\mu }_r\left(\frac{\partial u}{\partial y}+{\beta }^2\frac{\partial v}{\partial x}\right)$ , (26)

我们假设$\beta \ll 1$ 。且$Re \ll 1$ [23]，因此等式(23)~(26)可以简化为：

$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$ , (27)

$-\frac{\partial p}{\partial x}+{\mu }_r\frac{{\partial }^{2}u}{\partial y^2}+K^2\psi E_s=0$ , (28)

$-\frac{\partial p}{\partial y}=0$ . (29)

接下来，我们分别给出无滑移边界条件与无渗透边界条件，以及进出口压力条件：

$u|{}_{y=h(x)}=0,u|{}_{y=m(x)}=0,v|{}_{y=h(x)}=0,v|{}_{y=m(x)}=0,p|{}_{x=0}=p_{\text{in}}\text{ and }p|{}_{x=6}=p_{\text{out}}$ . (30a~f)

方程(27)~(29)具有显著的非线性特征。为获得流场的渐近解析解，我们采用摄动分析法。即将δ作为摄动参数。参照3.1的方法，我们将速度和压力变量展开为摄动级数。

$u=u_0+\delta u_1+O\left({\delta }^2\right),v=v_0+\delta v_1+O\left({\delta }^2\right),p=p_0+\delta p_1+O\left({\delta }^2\right)$ . (31)

将方程(31)代入方程(27)~(29)后，我们得到了零阶方程和一阶方程：

$\frac{\partial u_i}{\partial x}+\frac{\partial v_i}{\partial y}=0$ , (32)

$-\frac{\partial p_i}{\partial x}+{\mu }_r\frac{{\partial }^2{u}_i}{\partial y^2}+K^2{\psi }_i{E}_s=0$ , (33)

$-\frac{\partial p_i}{\partial y}=0$ . (34)

其中i = 0，1。此外，通过泰勒级数展开边界条件，我们可以得到：

$u_0|{}_{y=\pm 1}=0,v_0|{}_{y=\pm 1}=0,p_0|{}_{x=0}=p_{\text{in}},p_0|{}_{x=6}=p_{\text{out}}$ . (35a~d)

$u_1|{}_{y=1}=-\frac{\partial u_0}{\partial y}|{}_{y=1}\sin \left(n\pi x\right),u_1|{}_{y=-1}=\frac{\partial u_0}{\partial y}|{}_{y=-1}\sin \left(n\pi x+\theta \right),v_1|{}_{y=\pm 1}=0,p_1|{}_{x=0}=0,p_1|{}_{x=6}=0$ . (36a~e)

在考虑波纹几何结构时，压力梯度呈现非均匀分布，且流向势仍为未知常数[19]。因此，我们需要同时确定二维速度、压力及流向势。我们假设初始压力$p_0$ 与y方向无关，即$\frac{\partial p_0}{\partial x}=M_1$ 。根据该假设可得：

$p_0=M_{1}x+M_2$ . (37)

其中$M_1=\frac{1}{6}\left(p_{\text{out}}-p_{\text{in}}\right)$ ，$M_2=p_{\text{in}}$ 。现在可以得到零阶方程的解析解：

$u_0=c_1+c_{2}y+\frac{-\cosh \left(\frac{Ky}{\sqrt{{\epsilon }_r}}\right)\text{Sech}\left(\frac{K}{\sqrt{{\epsilon }_r}}\right)E_s{\epsilon }_r+\frac{1}{2}{y}^2{M}_1}{{\mu }_r}$ . (38)

根据边界条件(35a)：

$c_1=-\frac{M_1-2E_{\text{s}}{\epsilon }_r}{2{\mu }_r}$ , (39)

$c_2=0$ . (40)

我们可以看到，零阶的解对应于平行板微通道内发生的电动力学流动现象。而一阶解可表示为

$\begin{array}{c}{u}_1=c_3+c_{4}y+\frac{1}{{\mu }_r}(\frac{{\text{e}}^{\frac{K}{\sqrt{{\epsilon }_r}}+\frac{Ky}{\sqrt{{\epsilon }_r}}}K\left({\text{e}}^{\frac{2K}{\sqrt{{\epsilon }_r}}}\sin \left(n\pi x\right)-\sin \left(n\pi x+\theta \right)\right)E_s\sqrt{{\epsilon }_r}}{{\left(1+{\text{e}}^{\frac{2K}{\sqrt{{\epsilon }_r}}}\right)}^2}\\ +\frac{{\text{e}}^{\frac{K}{\sqrt{{\epsilon }_r}}-\frac{Ky}{\sqrt{{\epsilon }_r}}}K\left(-\sin \left(n\pi x\right)+{\text{e}}^{\frac{2K}{\sqrt{{\epsilon }_r}}}\sin \left(n\pi x+\theta \right)\right)E_s\sqrt{{\epsilon }_r}}{{\left(1+{\text{e}}^{\frac{2K}{\sqrt{{\epsilon }_r}}}\right)}^2}+\frac{1}{2}{y}^2{p^{\prime }}_1\left(x\right))\end{array}$ (41)

我们发现一阶速度$u_1$ 包含了未知的一阶压力。因此，我们必须先确定一阶压力。首先，利用连续性方程，将方程(41)代入方程(32)，可以得到一阶速度$v_1$ 的解析解，然后在纳米通道壁上考虑无渗透边界条件(36c~36e)。从而可以得到一阶压力的解析表达式：

$p_1=\frac{2n\pi \left(c_5+c_{6}x\right)+3\left(\cos \left(n\pi x\right)+\cos \left(n\pi x+\theta \right)\right)M_1-6\cos \left(n\pi x+\frac{\theta }{2}\right)\cos \left(\frac{\theta }{2}\right)E_s{\epsilon }_r\tanh {\left(\frac{K}{\sqrt{{\epsilon }_r}}\right)}^2}{2n\pi }$ (42)

整理可得：

$u_1=-\frac{1}{4{\mu }_r}\left[\begin{array}{l}\frac{1}{{\left(1+{\text{e}}^{\frac{2K}{\sqrt{{\epsilon }_r}}}\right)}^{2}n\pi }(\frac{1}{2}{\left(1+{\text{e}}^{\frac{2K}{\sqrt{{\epsilon }_r}}}\right)}^2\left(1+y\right)\left(-1+y\right)\left(-1+\cos \left(6n\pi \right)-\cos \theta +\cos \left(6n\pi +\theta \right)+2n\pi \left(-1+3y\right)\sin \left(n\pi x\right)\right)\\ +2n\pi \left(-1+y\right)\left(1+3y\right)\sin \left(n\pi x+\theta \right)M_1-4{\text{e}}^{\frac{K-Ky}{\sqrt{{\epsilon }_r}}}Kn\pi \left(\left(-1+{\text{e}}^{\frac{2K(1+y)}{\sqrt{{\epsilon }_r}}}\right)\sin \left(n\pi x\right)+\left({\text{e}}^{\frac{2K}{\sqrt{{\epsilon }_r}}}-{\text{e}}^{\frac{2Ky}{\sqrt{{\epsilon }_r}}}\right)\sin \left(n\pi x+\theta \right)\right)\\ e_s\sqrt{{\epsilon }_r}+{\left(-1+{\text{e}}^{\frac{2K}{\sqrt{{\epsilon }_r}}}\right)}^2\left(-1+y^2\right)\cos \left(\frac{\theta }{2}\right)\left(-\cos \left(6n\pi +\frac{\theta }{2}\right)+\cos \left(\frac{\theta }{2}\right)-6n\pi \sin \left(n\pi x+\frac{\theta }{2}\right)\right)e_s{\epsilon }_r)\end{array}\right]$ (43)

3.3. 流向势

在无外加电场作用时，电解液溶液会在流动过程中产生一个正向的流向电流和反向的电导流，进而形成流向势。在稳态条件下，流向电流I_s与电导流I_c的总和为零。

$I_s+I_c=0$ , (44)

此外，流向电流与电导流表示如下：

$I_s={\displaystyle {\int }_0^{6L}{\displaystyle {\int }_{m^{*}(x)}^{h^{*}(x)}{e}_0{z}_v\left(n_{+}-n_{-}\right)u^{*}\text{d}{y}^{*}\text{d}{x}^{*}}}$ , (45)

$I_c={\displaystyle {\int }_0^{6L}{\displaystyle {\int }_{m^{*}(x)}^{h^{*}(x)}{e}_0{z}_v\left(n_{+}+n_{-}\right)\frac{e_0{z}_v{E}_s^{*}}{f}\text{d}{y}^{*}\text{d}{x}^{*}}}$ , (46)

其中$n_{+}$ 和$n_{-}$ 分别为阳离子和阴离子的离子浓度，f为离子摩擦系数。阳离子$n_{+}$ 和阴离子$n_{-}$ 的离子数密度服从玻尔兹曼分布，并近似为Debye-Hückel线性化模型：

$n_{\pm }=n_0\exp \left(\mp \left(e_0{z}_v{\psi }^{*}\right)/k_{B}T\right)\approx n_0\left(1\mp \left(e_0{z}_v{\psi }^{*}\right)/\left(k_{B}T\right)\right)$ (47)

综合上述方程，我们得出：

${\displaystyle {\int }_0^{6L}{\displaystyle {\int }_{m^{*}(x)}^{h^{*}(x)}\frac{-2e_0^2{z}_v^2{\psi }^{*}{n}_0}{k_B{T}_a}{u}^{*}\text{d}{y}^{*}\text{d}{x}^{*}+}}{\displaystyle {\int }_0^{6L}{\displaystyle {\int }_{m^{*}(x)}^{h^{*}(x)}\frac{2e_0^2{z}_v^2{n}_0{E}_s^{*}}{f}\text{d}{y}^{*}\text{d}{x}^{*}}}=0$ . (48)

第一项代表流向电流，第二项表示电导流。使用方程(7)和方程(22)中给出的无量纲变量，方程(48)可以写成如下形式：

${\displaystyle {\int }_0^6{\displaystyle {\int }_{m(x)}^{h(x)}u\psi \text{d}y\text{d}x=R{E}_s}}{\displaystyle {\int }_0^6\left(h\left(x\right)-m\left(x\right)\right)\text{d}x}$ . (49)

这里$R=\left(\mu e_0^2{z}_v^2\right)/\left(\epsilon f{k}_B{T}_a\right)$ 是Péclet离子数。我们重新表述流速的表达式，提取流向势Es的项，使得$u=u_0+\delta u_1=u_p+E_s{u}_e$ 。结合速度表达式，我们推导出流向势：

$E_s=-\frac{{\displaystyle {\int }_0^6{\displaystyle {\int }_m^h{u}_p\psi \text{d}y\text{d}x}}}{R{\displaystyle {\int }_0^6\left(h\left(x\right)-m\left(x\right)\right)\text{d}x}+{\displaystyle {\int }_0^6{\displaystyle {\int }_m^h{u}_e\psi \text{d}y\text{d}x}}}$ . (50)

3.4. 能量转化效率

在微通道两端施加压力为流体提供的机械能可以被转化为电能。从数学角度而言，转化效率可表述为输出功率与输入功率的比值，即：

$\eta =\frac{P_{out}^{*}}{P_{in}^{*}}$ , (51)

$P_{out}^{*}=\frac{1}{4}\left|I_s^{*}{E}_s^{*}\right|$ , (52)

$P_{in}^{*}=\left|Q_{E_s^{*}=0}{\left(-\frac{\text{d}{p}^{*}}{\text{d}{x}^{*}}\right)}_m\right|$ , (53)

其中，$Q_{E_s^{*}=0}$ 为忽略流向势的流量，${\left(-\frac{\text{d}{p}^{*}}{\text{d}{x}^{*}}\right)}_m$ 为忽略流向势的平均压力梯度。经过计算，可得到无量纲能量转化效率：

$\zeta =\frac{3}{2}{K}^{2}R{E}_s^2\left|\frac{{\displaystyle {\int }_0^6{\displaystyle {\int }_m^h\text{d}y\text{d}x}}}{{\displaystyle {\int }_0^6{\displaystyle {\int }_m^h{u}_p\text{d}y\text{d}x{\displaystyle {\int }_0^6\frac{\text{d}p}{\text{d}x}\text{d}x}}}}\right|$ . (54)

其中$u_p$ 和$\text{d}p/\text{d}x$ 分别表示忽略流向势的流速和压力梯度。

3.5. 数值方法求解

我们已经采用小参数展开法推导了波纹壁微通道中电势、流速、压力及流向势的渐近解析解。这意味着该假设仅适用于振幅较小的情况。现在，我们将运用数值方法求解上述问题，得到上述问题的数值解，并与前文推导的解析解进行对比。根据方程(8)~(10)，利用“Mathematica”数学软件，我们得到了电势的数值解：

$\psi =\left({\text{e}}^{\frac{K(1-y+\delta \sin (n\pi x)}{\sqrt{{\epsilon }_r}}}\left(1+{\text{e}}^{\frac{K(2y-\delta \sin (n\pi x)+\delta \sin (n\pi x+\theta )}{\sqrt{{\epsilon }_r}}}\right)\right)/\left(1+{\text{e}}^{\frac{K(2+\delta (\sin (n\pi x)+\sin (n\pi x+\theta )))}{\sqrt{{\epsilon }_r}}}\right)$ (55)

将方程(55)代入方程(33)，(34)，同时结合方程(30a~b)中的边界条件，可确定x方向上的速度为：

$u=c_7+c_{8}y+\frac{1}{2{\mu }_r}\frac{\left(-\left(2\left({\text{e}}^{\frac{K(1-y+\delta \sin (n\pi x)}{\sqrt{{\epsilon }_r}}}+{\text{e}}^{\frac{K(1+y+\delta \sin (n\pi x+\theta )}{\sqrt{{\epsilon }_r}}}\right){\epsilon }_r{E}_s\right)\right)}{1+{\text{e}}^{\frac{K(2\delta (\sin (n\pi x)+\sin (n\pi x+\theta )))}{\sqrt{{\epsilon }_r}}}}+y^2{p}^{\prime }\left(x\right)$ (56)

其中：

$c_7=-\frac{1}{2{\mu }_r}\left(-2{\epsilon }_{r}Es+p^{\prime }\left(x\right)+\delta \sin \left(n\pi x\right)p^{\prime }\left(x\right)+\delta \sin \left(\theta +n\pi x\right)p^{\prime }\left(x\right)+{\delta }^2\sin \left(n\pi x\right)\sin \left(\theta +n\pi x\right)p^{\prime }\left(x\right)\right)$

$c_8=-\frac{1}{2{\mu }_r}\left(\delta \sin \left(n\pi x\right)p^{\prime }\left(x\right)-\delta \sin \left(\theta +n\pi x\right)p^{\prime }\left(x\right)\right)$

通过连续性方程，我们可以推导出y方向的速度v，再根据方程$v|{}_{y=h(x)}-v|{}_{y=m(x)}=0$ ，可以得到一个关于压力p的常微分方程：

$\begin{array}{l}\frac{-2}{{\mu }_r}(\delta {\epsilon }_r{E}_{s}n\pi \cos \left(\frac{\theta }{2}\right)\cos \left(\frac{\theta }{2}+n\pi x\right){\left(\tanh \left(\frac{K\left(2+\delta \left(\sin \left(n\pi x\right)+\sin \left(n\pi x+\theta \right)\right)\right)}{2\sqrt{{\epsilon }_r}}\right)\right)}^2\\ +\frac{1}{2{\mu }_r}\delta n\pi \cos \left(\frac{\theta }{2}\right)\cos \left(\frac{\theta }{2}+n\pi x\right){\left(2+\delta \left(\sin \left(n\pi x\right)+\sin \left(n\pi x+\theta \right)\right)\right)}^2{p}^{\prime }\left(x\right)\\ {+\frac{1}{12{\mu }_r}\left(2+\delta (\sin \left(n\pi x\right)+\sin \left(n\pi x+\theta \right)\right))}^3{p}^{″}\left(x\right))=0\end{array}$ . (57)

传统解析方法求解该常微分方程比较困难。为此，我们采用切比雪夫谱方法获得压力的数值解[24]。我们首先计算流向势。在u的数值解中，我们将涉及压力项的表达式单独写出，并代入流向势的相应表达式。随后采用牛顿迭代法求解流向势：首先设定一个初始值，当误差足够小时即视为收敛。此时即可得到流向势的值。将流向势代入压力和速度的表达式后，我们成功求解了波纹壁微通道内的流动特性。此外，将所得数值解与摄动解进行对比，结果显示出良好的一致性。

4. 模型验证

Figure 2. Comparison of analytical and numerical solutions for K = 2, δ = 0.01, n = 3/π, x = 0.25, Es = −0.02, θ = 0,φ = 0, where lines represent analytical solutions and dots represent numerical solutions. (a) Potential, (b) Velocity, (c) Pressure, (d) Pressure gradient

图2. 当K = 2，δ = 0.01，n = 3/π，x = 0.25，Es = −0.02，θ = 0，φ = 0时解析解与数值解的比较，其中线代表解析解，点代表数值解。(a) 电势 (b) 流速 (c) 压力 (d) 压力梯度

在图2中，系统地比较了电势、流速、压力和压力梯度的解析解与数值解。图2(a)中的电势解析解是通过方程(16)和(17)计算得出的，而数值解则通过方程(55)得出。两种计算方法在电势上表现出显著的一致性。但是由于解析解是在边界处渐近展开，而数值解的计算在边界处电势为±1，所以图2(a)在上下壁面处无法吻合，但是不影响整体的一致性。随后，将从方程(55)得出的速度数值解与方程(38)和(41)得出的解析解进行比较，接着比较压力和压力梯度，可以看出，这些比较结果都吻合，证明了本研究计算结果的正确性。需要注意的是，为简化计算，在图2(b)~(d)中假设流向势为常数，这可以避免流向势的复杂迭代过程。在后续计算中，我们包含了迭代过程以获得真实的电动力学流动特性。

5. 结果和讨论

为展示纳米流体速度分布的具体特征，需要对公式中相关物理参数取值范围进行定义。根据参考文献[18] [23]-[25]，颗粒均匀分散、无团聚的稀悬浮纳米流体中纳米颗粒的体积分数范围φ~0~0.3，电解质溶液介电常数ɛ_f = 6.995 × 10⁻¹⁰ CV⁻¹m⁻¹，Al₂O₃纳米颗粒介电常数ɛ_s = 80 CV⁻¹m⁻¹，玻尔兹曼常数k_B = 1.381 × 10⁻²³ J·m⁻¹K⁻¹，基液粘度μ_f = 0.93 × 10⁻³ kg/(m·s)⁻¹，密度ρ = 10³ kg/m³，电子电荷量e = 1.6 × 10⁻¹⁹ C，绝对温度T = 296.15 K，电解质溶液的离子化合价z = 1，离子摩擦系数f = 10⁻¹² N·s/m，平均半通道高度H = 50 nm，纳米通道长度为6 mm (即特征长度L = 1 mm)，电场强度E₀ = 10⁴ V/m，壁面电势为25 mV。雷诺数与长径比均足够小。无量纲双电层厚度K~1~15，R~0~1，无量纲波幅δ~0~0.5。

Figure 3. After adding nanoparticles, the velocity distribution at different cross-sectional positions of the pipeline (red indicates the wave crest, blue indicates the cross-section of the center of the corrugated surface, green indicates the wave trough). Triangles represent Newtonian fluid, while circles represent nano-fluid with φ = 0.25 (θ = 0, δ = 0.5, R = 0.1, K = 5, n = 3/π)

图3. 添加纳米颗粒后，管道不同截面位置的速度分布(红色表示波峰，蓝色表示波中，绿色表示波谷)。三角形代表牛顿流体，圆形代表φ = 0.25的纳米流体(θ = 0, δ = 0.5, R = 0.1, K = 5, n = 3/π)

在图3中可以看到，波谷处牛顿流体出现了回流现象，最大无量纲负速度为−0.0016，回流区域长度在距上下壁面10 nm处，大约占半通道的1/5。而添加φ = 0.25的纳米流体后速度方向没有发生改变，回流抑制效率达100%。根据Cho等人[5]的数值研究，纯压力驱动流在波谷处容易发生回流现象，因为通道两端压力差产生的流动动量可能不足以覆盖波纹表面的波谷区域，此时抛物线速度分布可能会被破坏。但如图3所示，当加入纳米粒子后，速度分布仍保持抛物线形态，这表明在通道横截面上未发生流动回流。在我们之前对牛顿流体的研究中发现，对于波纹型管道中的流动，波谷区域因管道截面积相对较大，流动易出现动量不足的问题，在反向的流向势等因素的影响下，局部流体易形成回流现象，此时速度剖面偏离经典的抛物线形态，在壁面附近出现速度反向的特征。然而，从图3中波谷处的速度分布曲线可以看出，纳米流体的速度剖面接近于抛物线形态，未出现明显的回流现象，这说明加入纳米粒子对流体的回流有抑制作用，这一抑制效应的出现是因为流体的黏度增大减少了波谷中间流体与近壁面流体的速度差，降低了局部流动分离的可能性；同时，纳米颗粒在流动过程中产生的附加剪切应力能够有效抑制近壁区域的反向动量积聚，破坏了回流涡流的形成条件，使得流体在波谷这一易产生回流的区域保持了稳定的流动状态。

图4(a)~(d)展示了不同壁面相位差下通道不同截面处的速度分布情况。在图4(a)中，可以看到对称波纹(反相)壁面的通道中不同横截面的速度分布，我们分别画了θ = π/2、π以及2/3π时的图像，即图4(b)~(d)。我们发现不同相位差导致了固定横截面速度分布的差异。图4(a)中，波峰处流速较大，而波谷处流速相对较小。这是因为当流体从波谷向波峰移动时，通道横截面逐渐减小。为满足质量守恒定律，波峰处的流速会比波谷处更大。同样的，在图4(b)~(d)中，当相位差发生变化时，管道的横截面积发生变化，从而使速度分布发生变化。因此，若某截面处面积因相位差而减小，该区域的流速将增加；反之，截面积增大的区域则伴随流速降低。这种由相位差引起的几何差异，通过改变流道的面积而影响速度。

Figure 4. Velocity at different cross-sectional positions of the pipe (Red indicates x = 0.5, blue indicates x = 1, and green indicates x =1.5). φ = 0, δ = 0.5, R = 0.1, K = 5, n = 3/π. (a) θ = 0, (b) θ = π/2, (c) θ = π, (d) θ = 2/3π

图4. 管道不同截面位置的流速(红色表示x = 0.5，蓝色表示x = 1，绿色表示x =1.5)。φ = 0，δ = 0.5，R= 0.1，K= 5，n= 3/π。(a) θ = 0，(b) θ = π/2，(c) θ = π，(d) θ = 2/3π

然而，考虑不同相位差微通道中心处流体的流速发现，当相位差θ = π时，管道中心处流体的速度高于其他形状微通道，最大速度达0.057。相较于(图4(a))的对称结构，最大速度增加约14%。这是由于θ = π时的壁面形状既避免了部分截面过度扩张导致的流速衰减，也克服了截面变化不均匀引发的局部流动扰动。这种管道设计具备优势：速度的提升意味着更强烈的动量交换，能够增强纳米流体中纳米颗粒的分散均匀性，减少颗粒聚集现象，从而充分发挥纳米颗粒对流体传热、传质性能的提升作用。

图5(a)展示了管道中心处速度随波数的变化规律。研究发现，速度极值出现在波数nπ = 0.1π附近。由于波长等于1/nπ，与波数为倒数关系，即在波长为3.2 mm左右时，速度最大。速度随波数的变化并非单调递增，而是存在一个或多个极值点。当波数从较小值逐渐增大时，速度通常先升至峰值，随后随着波数继续增大，速度会逐渐下降。这是因为当波数较小时，壁面波动的影响过于微弱，流体几乎不受波动影响，导致速度难以提升；而当波数较大时，壁面波动过于频繁，流体受到的阻力更大。在波数区间两端，速度可能趋近某个极限值，但通常不会趋近于零或无穷大。当波长与微通道或流体系统的固有特征尺度相匹配时，可能出现共振现象。在共振时，速度会出现明显峰值，其数值远大于其他波长下的速度。这是因为共振过程中，壁面波动与流体的相互作用达到最大强度，波动能量能最有效地转化为流体的动能，从而显著提升流速。

Figure 5. The velocity of the flow in the center of the pipe is related to the wave number (a) and the wave amplitude (b). (φ = 0.1, R = 0.1, K = 5)

图5. 对称管道中心处流速随 (a) 波数 (b) 波幅的变化关系。(φ = 0.1, R = 0.1, K = 5)

图5(b)展示了管道中心处速度与振幅之间的关系。从图中可以观察到，在振幅范围0~0.5范围内，随着振幅增大，管道中心处的流速整体呈现上升趋势。这表明振幅与管道流速之间存在正相关关系，即振幅增加会导致管道中心流速上升。此外，在相同振幅条件下，管道中心处速度随波数先增大后减小，我们发现，蓝色曲线的流速值高于红色和绿色曲线。综上，当波长为3.2 mm且振幅较大时，管道中心流速达到最大值。

Figure 6. (a) The relationship between flow potential and wavenumber (red indicates δ = 0.1, blue indicates δ = 0.25, green indicates δ = 0.5). (b) The relationship between energy conversion efficiency and wavenumber (δ = 0.1). (φ = 0, θ = 0, R = 0.1, K = 5)

图6. (a) 流向势随波数的变化关系(红色表示δ = 0.1，蓝色表示δ = 0.25，绿色表示δ = 0.5)。(b) 能量转化效率随波数的变化关系(红色表示δ = 0.1，蓝色表示δ = 0) (φ = 0, θ = 0, R = 0.1, K = 5)

波数作为壁面波动的关键参数，会显著影响流体与壁面的相互作用，从而改变流向势。如图6(a)所示，随着波数增大，流体流向势先增大后减小，最终趋于稳定状态。图6(b)中能量转化效率随波数的变化存在若干极值点：当波数较小时，能量转化效率较低；随着波数增大，效率逐渐提升；达到峰值后，继续增加波数会导致效率下降。这一现象进一步验证了我们在图5中得出的结论，所以，当波长为3.2 mm 且振幅较大时，管道中心处速度和流向势较大，且能量转化效率达到最大值0.0246，比平板通道效率提高约3.9%。

Figure 7. (a) The relationship between flow potential and nanoparticle concentration. (b) The relationship between energy conversion efficiency and nanoparticle concentration (δ = 0.5, n = 0.1, θ = 0, R = 0.1, K = 5)

图7. (a) 流向势随纳米粒子浓度的变化。(b) 能量转化效率随纳米粒子浓度的变化(δ = 0.5, n = 0.1, θ = 0, R = 0.1, K = 5)

图7(a)展示了流向势随纳米颗粒体积分数的单调变化规律。当φ增加时，流向势持续增大，纳米颗粒体积分数的提升会增大纳米流体的有效粘度，改变了流体介电特性，增强了双电层结构的极化效应，使通道上下游的电势差进一步扩大。在图7(b)可以看出，能量转化效率随纳米粒子浓度的增加而增加，加入纳米颗粒后，颗粒对回流的抑制效应降低了局部流动的能量耗散，减少了压力驱动机械能的损耗，使机械能的转化效率提升。

6. 总结

本研究通过理论推导与数值模拟，系统探讨了正弦波纹壁面微通道中纳米流体的流动速度、流向势及能量转化效率。结果表明，波纹壁面的几何参数(波幅δ、波数nπ、相位差θ)与纳米颗粒特性(体积分数φ)共同影响流场结构与能量转化。采用摄动法与数值求解相结合，首先推导了电势、速度、流向势及能量转化效率的解析解，并通过切比雪夫谱方法得到数值解，并验证了解的一致性。

• 纳米颗粒对回流的抑制作用：纳米流体的有效粘度高于基液，导致整体流速降低，但显著抑制了波谷区域的回流现象，增强了流动稳定性(图3)。

• 相位差的优化效应：当上下壁面波纹相位差θ = π (同相)时，通道横截面积均匀分布，管道中心处流速最大；非对称结构(如θ = π/2，2π/3)通道横截面积越大，流速越小，反之，通道横截面积越小，流速越大(图4)。

• 波数与波幅的协同影响：存在临界波长(λ = 3.2mm)使管道中心流速、流向势和能量转化效率达到峰值0.0246，比平板通道效率提高约3.9%。波幅δ增大(0~0.5范围内)可提升流速，但需兼顾阻力损失(图5，图6)。

• 能量转化机制：流向势Es由流向电流与传导电流决定，其极值出现在n = 0.1附近且波幅较大时；能量转化效率η随纳米粒子体积分数增加(图7)。

综上，纳米流体的引入可减少回流产生，增强流动稳定性且正弦波纹管道的几何设计(如较大波幅，临界波长、同相布局)可提升管道中心处速度并优化能量转化效率，为微流控系统中机械能–电能转化器件的设计提供了理论依据。

基金项目

本研究由国家自然科学基金地区项目(12262026, 12562027)、内蒙古自治区自然科学基金(2021MS01007, 2025MS01004)以及内蒙古自治区高等学校创新团队项目(NMGIRT2323)资助。

NOTES

^*通讯作者。