CN104850749B - 一种圆柱形微流道中电势分布数值的获取方法 - Google Patents

Abstract

现有技术中简单的利用有限差分法对电势分布进行处理求收敛值，而本发明则将电势分布分为大误差区和小误差区，对于不同的区域采用不同的处理方法如局部分区加密有限差分法和变步长有限差分法，通过简单的矩阵运算即可获得圆柱形微流道中双电层电势分布的精确计算结果，减小了小误差区的计算量，提高了大误差区的计算精度。从而有效解决了普通方法求解圆柱形微流道中双电层电势分布时无法获得收敛解的问题，以及微流道及液体的参数在特殊情况下的电势分布的快速、精确的计算问题。

Description

一种圆柱形微流道中电势分布数值的获取方法

技术领域

本发明属于微观技术领域，尤其涉及一种圆柱形微流道中电势分布数值的获取方法。

背景技术

随着微观技术的不断发展，微尺度下的流体流动问题逐渐引起人们的重视。研究发现，在直径小于几百微米的微流道中，由于管壁与流体之间微观作用力的存在，流体的实际流动规律与利用经典流体理论计算得到的结果差异较大。在微米级的管道中，静电力是主要的微观作用力之一，其产生原因是固液两相之间的双电层结构。当固、液相表面接触时，固液两相之间会形成电量相等、电性相反的电荷，即双电层结构。当液相发生流动时，双电层中在液相内的部分电荷会随着液相的流动发生移动，定向移动的电荷形成电流，称之为流动电流，电荷在管道两端积累形成电势，称之为流动电势。流动电势产生阻碍双电层中电荷移动的静电力，影响管道中流体的流动。当微流道尺寸小于几十微米时，由于双电层而形成的静电力变得不可忽略。知道双电层电势在管道中的分布规律对于描述微流道中流体运动规律、流动电流和流动电势都是必不可少的，由于微流道的尺寸限制，直接测量往往难以实现，需要通过数值模拟的方法对相应的问题展开研究。

计算双电层在微流道内的电势分布需要求解Poisson-Boltzmann方程，通常需要对该种非线性微分方程的边值问题进行研究。对于普通的边值问题，传统的打靶法、有限差分法(FDM，Finite Difference Method)和蒙特卡罗模拟常常可作为有效的求解工具，但是在特殊的参数条件下，传统的方法难以快速准确地获得圆柱形微流道中的电势分布。

发明内容

为解决上述问题，本发明提供一种圆柱形微流道中电势分布数值的获取方法，解决了双电层电势分布函数在特殊参数下的求解问题。

本发明的圆柱形微流道中电势分布数值的获取方法，其包括：

步骤1，对表示圆柱形微流道中电势分布的Poisson-Boltzmann方程进行线性化处理获得电势分布的线性微分方程，该线性微分方程表达出电势及其各阶导数在定义域[0,R]上单调，因此将[0,R]均分为n个均分区间，步长h表示为R为圆柱形微流道的半径；

步骤2，利用有限差分法将线性微分方程转换为线性方程组进行求解得到各均分区间电势的近似值，将近似值默认为各均分区间的电势ψ₁～ψ_n；其中，[0,h]对应第1个均分区间，第1个均分区间的电势为ψ₁，((i-1)h,ih]对应第i个均分区间，第i个均分区间的电势为ψ_i；

步骤3，对各均分区间电势ψ₁～ψ_n的近似值求和并计算平均值设定阈值且0.1＜μ＜1；然后将各均分区间电势的近似值与设定的阈值ψ_divide比较，找到r_divide，该r_divide满足当r_i＜r_divide时，ψ_i＜ψ_divide且当r_i≥r_divide时，ψ_i≥ψ_divide，从而将电势分布分为大误差区[r_divide,R]和小误差区[0，r_divide)；r_i为圆柱形微流道点在用柱坐标中的沿半径方向的坐标；

步骤4，逐次加密处理：

步骤41，在小误差区[0，r_divide)中对每个均分区间进行变步长的稀疏处理：从第一个均分区间开始，依次合并相邻的均分区间，当合并到某个均分区间时，合并的均分区间电势之和大于或等于ψ_divide，停止合并；再从下一个均分区间开始，重复合并操作，合并的均分区间组成稀疏区间；稀疏区间的步长为h＝r_b-r_a；a为稀疏区间中起始的均分区间，b为稀疏区间中结束的均分区间；

在大误差区[r_divide,R]中对每个均分区间进行变步长的加密处理：将每个均分区间重新均分为n_j个加密区间，且n_j为大于加密区间电势与ψ_divide比值的最小整数，则加密区间中步长为

步骤42，利用二次拉格朗日插值多项式来表示经步骤41划分的新区间的电势，所述新区间包括：稀疏区间或加密区间；

步骤43，对于大误差区[r_divide,R]中的每个新区间进行逐次加密：当次加密的步长为上一次加密的步长的二分之一，直到各区间电势的残差平方和满足设定的求解精度停止加密，获得逐次加密区间；

对于小误差区[0，r_divide)保留步骤41的分区；

步骤5，在经逐次加密处理后的各区上将所述的微分方程转化为线性方程组并进行矩阵运算，获得的电势计算结果即为电势分布的收敛值。

进一步的，所述步骤3中μ＝0.2。

本发明还提供一种圆柱形微流道中电势分布数值的获取方法，其包括：

步骤1，对表示圆柱形微流道中电势分布的Poisson-Boltzmann方程进行线性化处理获得电势分布的线性微分方程，该线性微分方程表达出电势及其各阶导数在定义域[0,R]上单调，因此将[0,R]均分为n个均分区间，步长h表示为R为圆柱形微流道的半径；

步骤2，利用有限差分法将线性微分方程转换为线性方程组进行求解得到各均分区间电势的近似值，将近似值默认为各均分区间的电势ψ₁～ψ_n；其中，[0,h]对应第1个均分区间，第1个均分区间的电势为ψ₁，((i-1)h,ih]对应第i个均分区间，第i个均分区间的电势为ψ_i；

步骤3，对各均分区间电势ψ₁～ψ_n的近似值求和并计算平均值设定阈值且0.1＜μ＜1；然后将各均分区间电势的近似值与设定的阈值ψ_divide比较，找到r_divide，该r_divide满足当r_i＜r_divide时，ψ_i＜ψ_divide且当r_i≥r_divide时，ψ_i≥ψ_divide，从而将电势分布分为大误差区[r_divide,R]和小误差区[0，r_divide)；r_i为圆柱形微流道点在用柱坐标中的沿半径方向的坐标；

步骤4，进行逐次加密处理：

在小误差区[0，r_divide)中保留步骤2的分区；

在大误差区[r_divide,R]中对每个均分区间进行逐次加密：当次加密的步长为上一次加密的步长的二分之一，直到各区间电势的残差平方和满足求解精度停止加密，获得逐次加密区间；

步骤5，在经逐次加密处理后的各区上将所述的微分方程转换为线性方程组并进行矩阵运算，获得的电势计算结果即为电势分布的收敛值。

其中，所述步骤3中μ＝0.2。

有益效果：

本发明利用模型中电势分布函数的单调性设计LRFDM和VSFDM两种新方法，通过简单的矩阵运算即可获得圆柱形微流道中双电层电势分布的精确计算结果，效率高，耗时少。可以有效解决普通方法求解圆柱形微流道中双电层电势分布时无法获得收敛解的问题，可以得到精确的计算结果。

现有技术中简单的利用有限差分法对电势分布进行处理求收敛值，而本发明则将电势分布分为大误差区和小误差区，对于不同的区域采用不同的处理方法如局部分区加密有限差分法和变步长有限差分法，减小了小误差区的计算量，提高了大误差区的计算精度。从而解决了微流道及液体的参数在特殊情况下的电势分布的快速、精确的计算问题。

附图说明

图1圆柱形微流道结构示意图；

图2(a)-(d)局部分区加密有限差分法步骤示意图；

图3(a)-(d)变步长有限差分法步骤示意图；

图4(a)-(d)FDM，LRFDM和VSFDM残差变化图；

图5圆柱形微流道中电势分布LRFDM算法结果图

图6圆柱形微流道中电势分布VSFDM算法结果图

具体实施方式

本文研究的圆柱形微流道的参数定义如图1所示。圆柱形微流道中的点用柱坐标表示(r,θ,z)，0≤r＜+∞,0≤θ≤2π,-∞＜z＜+∞，R为圆柱形微流道的半径。

双电层的电势分布和电荷之间的关系可用Poisson方程描述，如式(1)

其中ψ为双电层在管道中形成的电势，ρ_e为溶液中的电荷密度，ε₀为真空绝对介电常数，ε_r为介质的相对介电常数。式(1)的柱坐标形式如式(2)

根据圆柱形微流道和双电层的对称性，在稳态时，ψ与z和θ无关，则ψ(r,θ,z)＝ψ(r)，且易得故圆柱形微流道结构中的Poisson方程可写为

同时，双电层中的电荷由于电场力的作用，其热运动的分布满足Boltzmann分布律。

其中n_i(r)为微流道某处单位体积内具有we电荷离子的离子浓度，n₀为该离子在溶液中的离子浓度，w为每个离子所带电荷量，e为基本电荷量，k_b为波尔兹曼常数，T为绝对温度。

当液体中的每个正负离子所带的电荷量相等，符号相反，则对于正负电荷分别有

则液体中的某点的电荷密度为

联立(3)和(7)可得用于描述圆柱形微流道中电势分布的Poisson-Boltzmann方程，如(8)

式(8)为非线性微分方程，考虑到该式中的r，e和ψ均为小量，可将非线性微分方程线性化

式(10)边界条件为：

其中ζ为双电层的zeta电势。

对于式(10)的求解，当溶液中的离子浓度n₀很高、每个离子所带电荷量w很大、溶液的相对介电常数ε_r很小或在其他特殊参数条件下时，α的值会很大(α＞10⁶)，此时会导致ψ(r)在局部区域斜率过大，从而引入较大的计算误差。若利用常用的计算方法求解ψ(r)，无法得到精确、收敛的结果。本发明在有限差分法的基础上，针对式(10)的求解问题，设计了局部分区加密有限差分法(LRFDM，Local Refinement Finite Difference Method)和变步长有限差分法(VSFDM，Variable Step Finite Difference Method)，以获得精确、收敛的结果。这两种新算法适用的前提条件是

ψ(r)本身及其各阶导数在定义域(0,R]上单调

ψ(r)的大小与ψ″(r_i),ψ⁽⁴⁾(r_i)的大小具有一致性

(1)ψ(r)本身及其各阶导数在定义域(0,R]上单调，即：

定理1：若ψ(r)满足微分方程式(10)及其边界条件，则ψ(r)本身及其各阶导数在(0,R]上恒为正。

证明：对于式(10)可以令则有

则式(10)可以写成

式(13)的边界条件为：

式(13)为标准的虚宗量Bessel方程，在本发明中v＝0。式(13)的解为

ψ(x)＝c₁I_v(x)+c₂K_v(x) (14)

其中

I_v(x)，K_v(x)分别称为第一类和第二类虚宗量v阶Bessel函数。

式(14)的导数形式为

ψ′(x)＝c₁I_v′(x)+c₂K_v′(x) (18)

其中K_v(x)在x＝0处无界，故

由式(13)的边界条件ψ′(0)＝0易知c₂＝0，因此

由边界条件可知

由式(19)和(20)易得，ψ(x)本身及其各阶导数在上各阶导数为正，同理易得ψ(r)本身及其各阶导数在(0,R]上为正。至此定理1得证。

(2)ψ(r)的大小与ψ″(r_i),ψ⁽⁴⁾(r_i)的大小具有一致性，即：

定理2：若ψ(r)满足微分方程式(10)及其边界条件，且ψ(r₁)＜ψ(r₂)，r₁,r₂∈(0,R]，则有ψ″(r₁)＜ψ″(r₂)，ψ⁽⁴⁾(r₁)＜ψ⁽⁴⁾(r₂)。

证明：由定理1可知ψ′(r)＞0在(0,R]上恒成立，故当ψ(r₁)＜ψ(r₂)，则有r₁＜r₂。又由定理1可知，ψ⁽³⁾(r)＞0,ψ⁽⁵⁾(r)＞0在(0,R]上恒成立，故当r₁＜r₂，ψ″(r₁)＜ψ″(r₂)，ψ⁽⁴⁾(r₁)＜ψ⁽⁴⁾(r₂)必成立。至此定理2得证。

基于以上两个前提条件，设计局部分区加密有限差分法(LRFDM)和变步长有限差分法(VSFDM)。

LRFDM方法步骤

Step1：将[0,R]平均分为n份，步长为利用有限差分法进行预处理，利用差商代替导数，将线性微分方程转换为线性方程组进行求解，其中ψ_i表示ψ(r_i)的近似解：

将式(21)和式(22)代入式(10)，化简可得

对于边界条件的处理：

●终值条件ψ(R)＝ζ

令ψ_n+1＝ζ

●初值条件ψ′(0)＝0

采用前向有限差商公式。

式(23)、(24)和式(25)共含n个线性无关的方程组，联立可得矩阵方程如(26)

利用矩阵运算即可得到ψ₁～ψ_n的粗略近似值。如图2(a)，图中各变量用归一化的结果表示。

Step2：设定阈值ψ_divide，并将各点电势的近似值ψ_i与设定的阈值ψ_divide比较，可以找到r_divide，使得当r_i＜r_divide时，有ψ_i＜ψ_divide，当r_i≥r_divide时，有ψ_i≥ψ_divide，将区间[r_divide,R]作为加密区域，如图2(b)。

Step3：在加密区域[r_divide,R]中利用逐次分区加密的有限差分方法，即步长选择并利用式(21)～(26)所述的有限差分法进行处理，直到第k步与k+1计算得到各点电势的残差平方和满足求解精度，此时在[0,R]上总分区数为如图2(c)(d)。

VSFDM方法步骤

Step1：将[0,R]平均分为n份，步长为利用式(21)～(26)所述的有限差分法进行预处理，得到各点电势的粗略近似值，如图3(a)，图中各变量用归一化的结果表示。

Step2：对各点电势的近似值ψ_i求和，并计算平均值设定阈值0.1＜μ＜1。将各点电势的近似值ψ_i与设定的阈值ψ_divide比较，可以找到r_divide，使得当r_i＜r_divide时，有ψ_i＜ψ_divide，当r_i≥r_divide时，有ψ_i≥ψ_divide，如图3(b)。

Step3：在[0，r_divide)区间中对差分区间进行变步长的稀疏处理，当时将[r_a,r_b]重新划分为一个新的区间，步长为h^*＝r_b-r_a，如图3(c)。

Step4：在[r_divide,R]区间中对差分区间进行变步长的加密处理，对于每个ψ_i对应的区间[r_i-1,r_i)重新划分为n_i个区间，n_i为大于的最小整数，此时[r_i-1,r_i)区间中步长为对于两个相邻区间步长不相等的情形，利用二次拉格朗日插值多项式来表示电势，如式(27)。

则电势的一阶导数和二阶导数分别可表示为：

将式(28)和(29)代入(10)易知

其边界条件为：

利用式(30)～(32)并通过矩阵运算得到变步长情形下的有限差分计算结果。对每个区间[r_i-1,r_i)进行逐次加密直到第k步与k+1计算得到各点电势的残差平方和满足求解精度，此时在[0,R]上的总分区数为如图3(d)。

图1为本实例中针对的圆柱形微流道的结构示意图，在微流道中有稳定流动的液体。微流道为圆柱形，半径R＝50μm，长度L＝2mm，环境温度为20℃，即T＝293.15K，微流道中流动的液体为丙酮液体，相应的物理参数为zeta电势ζ＝5mV，电荷浓度分别为n₀＝6.022×10²¹/m³，n₀＝6.022×10²⁰/m³，n₀＝6.022×10¹⁹/m³，n₀＝6.022×10¹⁸/m³，液体的相对介电常数ε_r＝20.7(按照丙酮液体在20℃时的参数设定)。其他常用的物理参数如下：真空介电常数ε₀＝8.854×10^-12F/m，元电荷e＝1.6×10^-19C，波尔兹曼常数k_b＝1_.38×10^-23J/K。局部分区加密有限差分法LRFDM和变步长有限差分法VSFDM在Step1中的预分区数n＝500，残差用归一化后的结果表示，即获得最后结果的条件是err^*＜10^-5。

(1)算法精度及收敛速度验证

对于本专利所述的微流道，利用传统的有限差分法(FDM)和文中所述的局部分区加密有限差分法(LRFDM)及变步长有限差分法(VSFDM)进行处理，对比它们的精度和收敛速度。残差变化示意图如图4(a)-(d)所示。图4a至4d分别为n₀＝6.022×10²¹/m³、n₀＝6.022×10²⁰/m³、n₀＝6.022×10¹⁹/m³、n₀＝6.022×10¹⁸/m³下的FDM，LRFDM和VSFDM残差变化图。

易知传统的FDM因为计算量随迭代次数的增长迅速增大，在第三次迭代时即出现内存溢出的情况。当溶液离子浓度较大，即α较大时，FDM算法可以达到的计算精度较低，在n₀＝6.022×10²¹/m³时，err^*只能达到10^-4，随着溶液的离子浓度减小，即当α不断减小时，FDM可以达到的计算精度逐渐提高，在n₀＝6.022×10¹⁸/m³时，err^*可以达到10^-8。这种现象出现的原因为当α较大时，ψ(r)函数在R附近的斜率很大，利用FDM计算时此处的误差较大，而FDM通过逐次分半加密的方法来提高精度，在误差较小的区域浪费了计算量，因此导致内存溢出，而当α较小时，ψ(r)函数在R附近的斜率减小，因而FDM计算的误差减小，可以达到较高的精度。因此可知，FDM在α较大时计算式(10)难以达到较高的精度。

LRFDM在各离子浓度的条件下计算时err^*都能达到10^-9，并且计算量得到了很好的控制，不会出现内存溢出的情况。因此无论α的大小如何，LRFDM可以在很少的迭代次数中获得较高精度的解，因此LRFDM的收敛速度，算法精度以及算法的鲁棒性都要明显优于FDM。主要是因为LRFDM通过划分误差较大的区域进行分区加密，并没有在误差较小的地方浪费计算量，因而控制了整个算法的计算量，并使得大误差区域有更高的计算精度，进而使得整个算法的精度提高。

VSFDM从第一次迭代时就可以达到较高的精度，通过很少次数的迭代，各离子浓度的条件下计算时err^*都能优于10^-9，可见VSFDM在收敛速度和算法精度上都要优于LRFDM和FDM。主要是因为VSFDM比LRFDM有更好的分区方法，变步长的方法使得大误差区域可以得到更多的分区，所以可以更有效地提高了计算精度。

(2)算法计算速度验证

算法的计算速度不仅与迭代次数有关，还与每一次迭代的计算量有关，通过软件仿真，FDM迭代3次时即会发成内存溢出的情况，计算量随迭代次数显著增加，效率低，速度慢。此处对FDM，LRFDM和VSFDM计算速度进行对比。并以离子浓度为n₀＝6.022×10²¹/m³时的计算速度为例，算法完成每个步骤的耗时如表1所示(数据保留四位有效数字)。

当然，本发明还可有其他多种实施例，在不背离本发明精神及其实质的情况下，熟悉本领域的技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

Claims (4)

1.一种圆柱形微流道中电势分布数值的获取方法，其特征在于，包括：

步骤1，对表示圆柱形微流道中电势分布的Poisson-Boltzmann方程进行线性化处理获得电势分布的线性微分方程，该线性微分方程表达出电势及其各阶导数在定义域[0,R]上单调，因此将[0,R]均分为n个均分区间，步长h表示为R为圆柱形微流道的半径；

步骤2，利用有限差分法将线性微分方程转换为线性方程组进行求解得到各均分区间电势的近似值，将近似值默认为各均分区间的电势ψ₁～ψ_n；其中，[0,h]对应第1个均分区间，第1个均分区间的电势为ψ₁，((i-1)h,ih]对应第i个均分区间，第i个均分区间的电势为ψ_i；

步骤3，对各均分区间电势ψ₁～ψ_n的近似值求和并计算平均值设定阈值且0.1＜μ＜1；然后将各均分区间电势的近似值与设定的阈值ψ_divide比较，找到r_divide，该r_divide满足当r_i＜r_divide时，ψ_i＜ψ_divide且当r_i≥r_divide时，ψ_i≥ψ_divide，从而将电势分布分为大误差区[r_divide,R]和小误差区[0，r_divide)；r_i为圆柱形微流道点在用柱坐标中的沿半径方向的坐标；

步骤4，逐次加密处理：

步骤41，在小误差区[0，r_divide)中对每个均分区间进行变步长的稀疏处理：从第一个均分区间开始，依次合并相邻的均分区间，当合并到某个均分区间时，合并的均分区间电势之和大于或等于ψ_divide，停止合并；再从下一个均分区间开始，重复合并操作，合并的均分区间组成稀疏区间；稀疏区间的步长为h＝r_b-r_a；a为稀疏区间中起始的均分区间，b为稀疏区间中结束的均分区间；

在大误差区[r_divide,R]中对每个均分区间进行变步长的加密处理：将每个均分区间重新均分为n_j个加密区间，且n_j为大于加密区间电势与ψ_divide比值的最小整数，则加密区间中步长为

步骤42，利用二次拉格朗日插值多项式来表示经步骤41划分的新区间的电势，所述新区间包括：稀疏区间或加密区间；

步骤43，对于大误差区[r_divide,R]中的每个新区间进行逐次加密：当次加密的步长为上一次加密的步长的二分之一，直到各区间电势的残差平方和满足设定的求解精度停止加密，获得逐次加密区间；

对于小误差区[0，r_divide)保留步骤41的分区；

步骤5，在经逐次加密处理后的各区上将所述的微分方程转化为线性方程组并进行矩阵运算，获得的电势计算结果即为电势分布的收敛值。

2.如权利要求1所述的圆柱形微流道中电势分布数值的获取方法，其特征在于，所述步骤3中μ＝0.2。

3.一种圆柱形微流道中电势分布数值的获取方法，其特征在于，包括：

步骤1，对表示圆柱形微流道中电势分布的Poisson-Boltzmann方程进行线性化处理获得电势分布的线性微分方程，该线性微分方程表达出电势及其各阶导数在定义域[0,R]上单调，因此将[0,R]均分为n个均分区间，步长h表示为R为圆柱形微流道的半径；

步骤2，利用有限差分法将线性微分方程转换为线性方程组进行求解得到各均分区间电势的近似值，将近似值默认为各均分区间的电势ψ₁～ψ_n；其中，[0,h]对应第1个均分区间，第1个均分区间的电势为ψ₁，((i-1)h,ih]对应第i个均分区间，第i个均分区间的电势为ψ_i；

步骤3，对各均分区间电势ψ₁～ψ_n的近似值求和并计算平均值设定阈值且0.1＜μ＜1；然后将各均分区间电势的近似值与设定的阈值ψ_divide比较，找到r_divide，该r_divide满足当r_i＜r_divide时，ψ_i＜ψ_divide且当r_i≥r_divide时，ψ_i≥ψ_divide，从而将电势分布分为大误差区[r_divide,R]和小误差区[0，r_divide)；r_i为圆柱形微流道点在用柱坐标中的沿半径方向的坐标；

步骤4，进行逐次加密处理：

在小误差区[0，r_divide)中保留步骤2的分区；

在大误差区[r_divide,R]中对每个均分区间进行逐次加密：当次加密的步长为上一次加密的步长的二分之一，直到各区间电势的残差平方和满足求解精度停止加密，获得逐次加密区间；

步骤5，在经逐次加密处理后的各区上将所述的微分方程转换为线性方程组并进行矩阵运算，获得的电势计算结果即为电势分布的收敛值。

4.如权利要求3所述的圆柱形微流道中电势分布数值的获取方法，其特征在于，所述步骤3中μ＝0.2。