（流体力学专业论文）毛细管电泳通道中的弯道效应及其消除.pdf

浙江大学硕士学位论文 聂德明2 0 0 4 2 摘要 首先，本文首次建立了槽流中对流扩散运动的格子b o l t z m a n n 模型，对l b m 中的速度边界条件、浓度边界条件进行了初步的探索为更好地满足速度无滑移 条件，对标准的壁面反弹边界条件进行了改进，提高了计算精度。此外，本文利 用该模型分别对槽流中的纯扩散和对流扩散进行了数值模拟，所得结果- 9 t a y l o r a r is 扩散理论相符 其次，本文对毛细管电泳分离中的主要的流体驱动方式一电渗驱动进行了深 入地研究，建立了完整的电渗流数学模型，并采用涡量一流函数方法求解流场。 对直槽道和18 0 度弯槽道中的电渗流型及影响参数进行了分析 最后，文章综合分析了引起弯道效应的诸多因素，对几种消除弯道效应的方 法进行了比较为了克服这些方法的不足之处，本文提供了一种新的消除弯道效 应的方法，即通过改变弯道内外壁面上的f 电势分布来改变速度分布，通过速度 差异来平衡迁移长度的差异，进而克服弯道效应同时，本文建立了描述弯道效 应的数学模型，将样品的对流扩散运动替代为跟随粒子的对流运动，通过粒子分 布的离散性和对称性双方面来定量描述弯道效应，并在此基础上给出了弯道上 0 5 ) 。 o 2 3i b m 中的边界条件 宏观的流体力学方程组可以由流动的速度、压力、浓度等实际情况，较为方 便地真接给出边界条件。而l b m 需要从微观角度给出分布函数的边界条件，并 且使之表现出的宏观性质要满足宏观流动的边界条件。由于l b m 求解的速度分 布函数疗的自由度通常比宏观量( 如密度、速度等) 大，为了满足边界条件，疋 可以有多种组合，因此就存在怎样恰当准确地给出边界上分布函数的问题。 2 3 1 速度边界条件 在l b m 中为了获得壁面无滑移速度条件通常采用固壁反弹边界条件，具体 是指在壁面格子上的流向壁面的粒子按原方向返回，即离开壁面方向上分布函数 的值赋为其相反方向上分布函数的值。反弹边界条件能够保证壁面的质量守恒和 无滑移速度。然而，近年来的一些研究 2 1 1 2 2 1 证明了反弹边界只能获得一阶精度。 浙江大学硕士学位论文 聂德明2 0 0 4 2 这就降低了l b m 的精度，因为文献【玎】【2 4 】表明，在低m a c h 数下l b m 在内网格 点具有二阶精度。为了获得具有二阶精度的无滑移边界条件，c h c n 等人【2 5 】将 l b m 看作一种特殊的有限差分方法，在固壁边界外再加一层，而这一层的分布 函数利用插值获得：z o u 和h c 2 q 将标准的全反弹边界条件改为只有法向反弹边 界，也改进了l b m 的精度；文献口y l 2 8 1 将壁面置于第一层流体节点和反弹节点层 中间，能够得到二阶精度。 本文对标准的反弹边界条件进行了改进，让壁面节点参与碰撞，并对部分分 布函数f 进行了重分配，使之也达到了二阶精度的要求。 bf c 蕊； 2 兹 固壁 ed 图2 6 壁面反弹边界条件 图2 6 表示槽道的下壁面，对于固壁上的节点a ，当经过流过程后，由a 邻 近的b 、c 、d 、e 、f 点可以分别确定格点a 的5 个方向的分布函数：五，二， ，石，五。 根据壁面反弹的思想可知： 五= 工，正= 石，五= 五( 2 2 6 ) 则点a 所有方向的分布函数可以确定；标准的壁面反弹方法同时还要求壁面节 点不进行碰撞，以保证无滑移条件。然而这种处理方法得到的精度有限。如果让 固壁节点参与碰撞，可以看到反弹条件并不能保证壁面切向速度为零，因为 彳工；根据文献，为了满足无滑移条件，可以对部分分布函数，进行重分配： ( 1 ) 利用反弹边界初步确定所有的分布函数，然后根据( 2 7 ) 计算出p ，u ； ( 2 ) 对具有切向分量的流入流体内部的分布函数，作如下分配： = 一唧m e i ( 2 2 7 ) 其中口为常数，对于d 2 q 9 模型，口= 去；以下壁面为例，可得： 浙江大学硕士学位论文聂德明2 0 0 4 ，2 正= 五一i 1p ( “+ v ) 。 ( 2 2 8 ) 五：五+ 三p ( 一v ) 代入到方程( 2 7 ) 中，可知“= 0 ，v = 0 ，满足无滑移条件。其它壁面情况可作类 似处理。至于进出口条件，既可以采用插值法也可以采用法向反弹法，详见文献 【2 5 】f 2 6 1 。 2 3 2 浓度边界条件 一由于固壁上样品的浓度不为零，不能应用壁面反弹法。根据壁面无渗透条件， 可以采用镜面反射方法f 2 2 确定未知分布函数。如图2 7 所示，以下壁面为例，当 流动过程结束后，l 、3 、4 、7 、8 方向的分布函数己知，通过镜面反射， 石= f 4 ，石= 五，五= 石，可以确定边界上所有方向的分布函数。至于进出1 = 的边 界条件，可以假定样品在无穷远处的浓度为零。因此，进出口处( 工= o ，j = ) 的 分布函数可以全部设为零。 bfc 。 囤壁 e j 1 。 4 。 图2 7 镜面反射边界条件 2 4 算例 2 4 1 直槽道中的纯扩散 如图2 4 所示，当槽流速度u 为零时，样品仅发生扩散运动。描述此过程的 方程为： 箜：d c( 2 2 9 ) 研 在初始条件c ( x = x o ，t = o ) = c d 下有解析解口瑚： 浙江大学硕士学位论文 聂德明2 0 0 4 2 门 c = 兰骗p - ( x - ) 2 “2 ( 2 3 0 ) 2 n - c r 浓度曲线呈高斯分布，其中盯= 2 d - r ，盯是浓度的标准偏差。 计算中取= 1 0 ，f = 0 8 ，c o = 1 0 ，占= 0 0 1 ，r = 1 0 0 ( 格子单位) ：计算 结果如图2 8 所示。由d ：至岩万可得样品扩散系数的理论解见：1 o 口一3 。根 。 据图2 8 及以下两式： = e “ ) d x ( 2 - 3 1 ) 盯2 = e ( x x o ) 2 c 似) d x ( 2 - 3 2 ) 可算得x o = 1 0 ，一= 2 0 3 e 3 ，则n = 1 0 1 5 e 一3 ，可见数值解和理论解符合得很好。 图2 8 样品浓度沿轴向呈高斯分布图2 9 样品浓度曲线随时间的变化 图2 9 给出了不同时刻浓度分布的高斯曲线。随着时间的推移，样品逐渐向 两侧扩散，浓度曲线的带宽( 盯) 也逐渐扩大。在电泳分离的过程中控制其分离 效果实际上就是控制其浓度曲线的带宽，希望曲线越窄小越好。 2 4 2 直槽道中的对流扩散 当平板间的流体速度u 不为零时，样品将发生扩散和对流。首先采用c h e n 侣1 提供的方格d 2 q 9 模型模拟槽道内的p o i s e u i l l e 流，即( 2 1 0 ) 式。速度边界采用 改进的壁面反弹条件。 如图2 1 0 ，抛物型速度剖面与分析解吻合得很好1 3 1 1 ，计算时取 f = 0 8 ，r e = 0 1 。图2 1 1 给出了槽道中轴向的密度分布，平均密度为2 7 ，说明 浙江大学硕士学位论文聂德明2 0 0 4 2 。 图2 】0 标准化槽流速度剖面图2 】1 沿槽道轴向密度分布 压力是线性分布的。因为l b m 在近似不可压假定下有，p = c ，2 p ，t 声速。 是 图2 1 2 不同时刻c o u e u e 流速度剖面图2 1 3 改进速度边界对应的误差收敛精度 图2 1 2 给出了c o u e t t e 剪切流随时间变化的速度剖面，在t = 3 0 0 0 时槽道内 的速度已经成线性分布，说明达到了稳定状态。为了说明改进的速度边界能够使 计算达到二阶精度，文中通过逐步精化网格来观察误差收敛规律。如图2 1 3 ，咖 表示y 方向的格点数，当格点数逐渐增大时，均可以得到误差与格点数的关系： e r r 。c 眵一2 ( 2 3 3 ) 说明达到了二阶精度。从图中还可以发现，对于相同的松弛时间r ，r e 数越大， 精度越低，这也是l b m 的一个缺陷，因为l b m 只有在低m a c h 数才能获得较 好的结果。对于相同的r e 数，当松弛时间f 1 时，f 越大精度越低，这是因为 浙江大学硕士学位论文聂德明2 0 0 4 2 应用c h a p m a n - e n s k o g 展开推导n - s 方程是在假定r 值较小的条件下进行的。r 值 过大，就很可能使l b m 偏离n s 方程的解。 当样品进入已稳定的p o i s e u l l e 流时，在流体的作用下开始运动。p o i s e u l l e 流抛物线型的速度剖面使样品具有类似剖面的浓度分布，即处于流场中间部分的 样品运动的速度大于两旁的样品，随着时间的推移，中间部分的样品与壁面的样 品之间的距离越来越大；但是样品自身还存在着扩散运动，在两者的共同作用下 就产生如图2 1 4 所示的样品浓度分布，图中颜色越深，表示浓度越大3 2 1 。此算 例中所取参数为x o = 5 0 ，r = 0 6 ，c o = 1 0 ，占= o 0 1 a 二二 二二二二二二二二二o 二二 二二二二二二二二瞄。 二二】二二二二二二二二o o ”l , , 1 3 0 0 图2 1 4 样品的对流扩散运动 根据t a y l o r - a r i s t l 9 l 【2 0 】理论，样品扩散系数d 随流体平均速度的平方 。= p 。+ 玄u 蔷2 1 i 2 呈线性变化。= 哦+ 暑蔷，吃为分子扩散系数。定义p e 数， n = 爱，则。= 巩1 1 + 堡2 1 0 ) 1 。计算中参数，c o ，占，f 的选取与2 4 相同，r = 3 。 改变流体速度2 ，并利用2 4 中所采取的方法计算出对应的样品扩散系数d 。 结果如图2 1 5 所示，数值计算的结果比t a y l o r 理论预测的结果偏低，但相对误 差小于2 ，计算结果与理论结果基本符合。 心数表示了对流与扩散作用的比值。当不断增大p e 数时，对流的效应不断 增大，而扩散的作用不断减小。图2 1 6 给出了p e 数从o 2 4 到2 6 4 时对应的样 品浓度分布曲线，可以看出，当p e 增大时，浓度分布曲线更加趋近于对称，且 更加接近于高斯曲线分布，这说明样品浓度分布的变化更小，更易于分离。计算 中取r = 1 0 ，= 0 5 0 ，并通过改变流体速度“。来改变心数：当样品运动到 x l = x 0 + “。t = 0 5 0 6 6 时进行检测，求出样品轴向的平均浓度分布。 塑竖查堂堡主兰丝堡兰 垂堡塑! ! 竺：! 图2 1 5 扩散系数与平均速度成二次关系图2 1 6 样品浓度分布随p e 数的变化 1 9 浙江大学硕士学位论文 聂德明2 0 0 4 2 第三章电渗流的数值模拟 随着微通道装置在生物样品、化学试剂的分离、鉴别及合成等领域的广泛应 用，微通道内流体驱动方式的研究也逐渐成为一研究热点。尤其是近几年生物芯 片技术的进步和“l a b o n a - c h i p ”概念的提出更是迫切要求实现微尺度下流体的自 动、精确的驱动和控制。微流体的驱动与控制和宏观流体的驱动与控制有很大不 同，一这主要是由于当尺度减小时，流体的流动特性发生了变化。这种流动特性的 变化使得宏观流体驱动与控制技术在微流体中的简单移植往往不成功或者效果 不好，微流体的驱动与控制技术更为复杂和多样化【3 3 】。 1 9 9 2 年h a r r i s o n 2 1 等人将样品注射系统毛细管电泳装置集成到平板玻璃芯片 上，他们发现通过施加适当的电压能驱动流体在毛细管内流动。这就是电渗流 ( e l c t r o o s m o t i cf l o w ) 。后来，s e i l e r t 3 4 l 等人对h a r r i s o n 等人的实验技术及装置作 了改进。1 9 9 3 年，h a r r i s o n 等人1 3 5 1 在两条垂直交叉的十字型微通道中成功地实 现了电渗流驱动及样品的电泳分离。次年，f a n 和h a r r i s o n t 3 6 】将这一系统集成到 了玻璃芯片上并再次证明了它的有效性及发展潜质。近几年，许多研究者包括 s e i l e r 等h t ”1 、e 虢n h a u s e r 等人1 3 8 1 、b u 玛g r a f 等人【3 9 l 、j a c o b s o n 等人1 4 0 】及c u l b e r t s o n 等人【4 l 】从实验角度对十字型毛细管电泳装置及电渗技术进行了改进和完善，并 将这一技术广泛应用到微流控芯片、生物芯片等微型分析化学系统中。关于电渗 流的理论研究方面，早在1 9 6 5 年r i c e 和w h i t e h e a d p 2 l 就推导出了描述毛细管中 电渗流的方程。后来j o r g e n s o n 和l u k a c s l 4 3 1 也建立了一维毛细管电泳模型。1 9 9 3 年，a n d r e e v 和l i s i n 4 a 不仅给出了毛细管电泳的一维数学模型，同时研究了电 渗流的速度剖面对分离效果的影响。n e e l e s h 和h o w a r d 4 5 1 对十字型注射及电泳分 离系统进行了数值模拟，研究了不同条件下微通道内的电渗流型。e r i m a k o v 等 人【4 6 1 和q i u 等人【4 7 1 对电渗流也进行了理论分析。 由上观之，毛细管电泳技术的高速发展主要得益于电渗驱动技术的利用，当 然同时也促使了电渗技术的不断完善。应该说，电渗流是目前最成功的微流体驱 动和控制方法之一。与压力等驱动方式相比，电渗流在微通道装置中具有以下优 点【4 5 】【4 8 】1 4 9 】： 浙江大学硕士学位论文聂德明2 0 0 4 2 ( 1 ) 一般而言，电渗流的速度大小与管道或槽道的横向尺寸无关，易于控制； 而压力流的速度除了与压力梯度有关之外，还与管道或槽道的横向尺寸 有关，为了保持一定的流速，需要考虑双重因素。 ( 2 ) 电渗流在管道或槽道中的横向速度剖面几乎是平直的，有“平流泵”之作 用。这样的速度剖面有利于样品的分离，即使在通道内传输很长距离样 品的浓度带宽变化也很小。而压力驱动流将产生抛物线型速度剖面，沿 通道横向的速度梯度比较大，不利于样品的高效分离。 ( 3 ) 电渗流主要通过施加电压驱动带电流体，因此可以控制电压来控制流速。 利用电压的切换可以在微通道的交叉口控制电渗流的方向，实现阀的功 能。优化通道的几何结构，还可以在微流装置的不同部位产生不同的流 速。这在生化分析中，例如液体的混和及多样品的并行处理中很有用处。 对于压力等驱动方式，通常需要安装微泵、微阀等装置，在工艺加工及 维修方面比较困难。 当然，电渗流的驱动与控制也存在着一些缺点【3 3 l ： ( 1 ) 电渗流对管道或槽道壁面材料和被驱动流体的物理化学性质有要求，与 液体接触的表面材料必须能够提供电荷，以形成双电层。因此它只适合 一定范围内的流体和管壁材料。 ( 2 ) 产生电渗流所需要的高电压电源会带来安全、功耗和所占空间大的问题， 这不利于系统的微型化。 ( 3 ) 电渗流尽管适合于驱动和控制狭窄管道或槽道的微量液体，但由于焦耳 热问题，它却不能高速驱动更宽管道或槽道中的流体，而这一能力在许 多微流体应用( 例如样品的预处理) 中是有必要的。 3 1 电渗流产生的原理 电渗现象是一种宏观现象，它是指在电场作用下，管道中或固相多孔物质内， 液体沿固体表面移动的现象【3 3 1 。图3 1 说明了电渗流形成的原理。电渗流产生的 前提是在与电解液接触的管壁上有不动的表面电荷，这种表面电荷来自于离子化 基或是液体中被强力吸附在表面的电荷( 例如，石英毛细管壁上的表面负电荷来 自于硅羟基在水溶液中发生电离所产生的研口一负离子) ，在表面电荷的静电吸附 浙江大学硕士学位论文 聂德明2 0 0 4 2 和分子扩散的作用下，溶液中的抗衡离子就会在固液界面上形成双电层( e l e c t r i c d o u b l el a y e r , e d l ) ，而管道中央液体中的净电荷则几乎为零。如图3 1 ( b ) ，根据 s t e m 模型，双电层由紧密层和扩散层组成。其中紧密层的厚度为1 2 个离子的 厚度。当在管道两端施加适当的电压时，在电场作用下，固液两相就会在紧密层 和扩散层之间的滑动面发生相对运动。由于离子的溶剂化作用及粘滞力的作用， 当形成扩散层的离子发生迁移时，这些离子就会携带着液体一同移动，因此形成 了电渗流。液体随扩散层中的离子移动，从开始到形成稳定的速度轮廓，所需的 时间很短。数值计算表明，这一时间大约为1 0 0 t s 1 m s 。在这段时间之后，电 渗流的流速轮廓是一个平面，就像一个瓶塞，这一点与压差引起的抛物型的流体 溶液，、 g o f = j 0 o o o 0 0 o 紧扩 藿鏊 ( a ) 外加电压作用下的电渗流( b ) 双电层电荷分布 图3 ，l 电渗流形成原理 流速轮廓不同，如图3 2 ，电渗流沿槽道横向的速度梯度很小，除了壁面附近外， 其余各处的速度几乎相等。 ( a ) 压力驱动流抛物型速度剖面( b ) 电渗驱动流瓶塞速度剖面 图3 2 两种不同驱动方式的速度剖面 影响电渗流的主要因素有双电层的结构特征及缓冲溶液的属性 3 3 1 。前者由 槽道壁面的化学属性及样品在内壁面上的吸附能力决定，后者主要由缓冲溶液的 组成、浓度、p h 值及温度决定。在槽道及溶液性质一定的情况下，电渗流主要 浙江大学硕士学位论文 聂德明2 0 0 4 2 与外加电场有关。 3 2 控制方程及边界条件 根据流场k n u d s e n 数大小，流体及其运动特性的描述可以划分成三个层次 ，即 ( 1 ) k n u d s e n 0 o l ，宏观的层次，通过建立在连续介质假定的基础上建立起 来的n a v i e r - s t o k e s 方程组进行描述。 ( 2 ) o 0 1 k n u d s e n 1 0 ，分子层次，从分子动力学模型( m o l e c u l a rd y n a m i c sm o d e l ) 出发进行描述。求解方法有直接m o n t ec a r l o 方法，i p ( i n f o r m a t i o n p r e s e r v a t i o n ) 方法。 对于目前微米微通道内的大多数液体的运动，其k n u d s e n 数通常小于o 0 1 或在其左右，因此仍然属于连续介质区，流体可以当作连续介质考虑。所以本文 仍然采用n - s 方程和无滑移边界条件来描述流场。 3 2 1 控制方程 假定流体不可压，连续性方程和运动方程【4 5 】： v v = 0 ( 3 1 ) p ( 詈+ ( v v ) v ) = 册v 懈 z ， 其中v 是流速，p 是密度，p 是压力，是粘性系数，见是电荷密度，e 是电 场强度。展和e 由下式决定： v 2 中= 钐 ( 3 - 3 ) e = 一v o ( 3 4 ) 式中，占是溶液的介电常数。将( 3 3 ) 、( 3 4 ) 式代入到( 3 2 ) 式，得， p ( 詈+ ( v v ) v = 一即v + 占( v 2 中v o ) s ， 从3 1 知，与流体相接触的通道壁面上带了电荷，使壁面产生定电势，称 之为z e t a 电势( f 电势) a 电渗流主要就是通过施加外加电场使壁面附近的带电 浙江大学硕士学位论文 聂德明2 0 0 4 2 溶质发生迁移，从而带动整个流体运动。通常情况下，通道中的电荷分布由壁面 电势及外加电场两者共同决定。根据文献【4 5 1 ，当d e b y e 厚度 ( 见图3 - 3 ，定义 见3 2 2 ) 很小及壁面电荷为少量时可以认为通道中的电荷分布主要由f 电势决 定，而与外加电场无关。当然，由于带电溶质随流体一起运动，因此电荷分布可 能会受流体的运动特性影响。但是当流体速度很小，即运动方程中惯性项很小或 d e b y e 厚度很小时，流体运动对通道中的电荷分布的影响可以忽略【5 0 】。因此，为 了简化计算模型，将电势。分解为两个部分【4 5 1 ：外加电场产生的电势西及壁面电 荷产生的电势妒，即 中= + 妒 ( 3 6 ) 方程( 3 - 4 ) 相应地分解为两个方程： v 2 西= 0( 3 7 ) v 2 伊一钐 ( 3 - 8 ) 因此，产生电渗流所需的电场强度e 为： e = 一v 西( 3 - 9 ) 根据文献【4 8 1 ，电荷密度织的分布与电势妒之间的关系可以通过平衡态b o t t z m a n n 分布表示， 成= 一2 f z c 。s i n h ( z f 妒r t ) ( 3 - 1 0 ) 式中，f 是f a r a d a y 常数，z 是离子电荷数，c 。是离子体积浓度，r 是气体常数， r 是温度。将上式代入到( 3 8 ) ，得： v 2 伊；2 f z c 3 s i n h ( z f o 胄，) ( 3 。11 ) 将( 3 9 ) ，( 3 一1 1 ) 式代入动量方程( 3 5 ) ，得： p ( 罾- i - ( v v ) v ) 却+ a v 2 v + 孚s i n h ( z f 邮丁，r e ( 3 - 1 2 ) 方程( 3 7 ) ，( 3 - 1 1 ) ，( 3 - 1 ) 及( 3 - 1 2 ) 耦合求解就能得到电渗流场。计算中首 先通过( 3 7 ) ，( 3 1 1 ) 计算出电势分布及电荷密度分布，然后代入到( 3 1 2 ) ， 并联合( 3 1 ) 就可以求出计算区域中的流线和涡量分布，然后根据流函数与速 度的关系就可以确定速度场了。 浙江大学硕士学位论文 聂德明2 0 0 4 2 3 2 2 无量纲化 如图3 3 所示的流场，采用下列无量纲参数： x。v尼 z 2 i y 。言，妒。面妒 妒= 詈，v = 罟，p 2 鲁 伊一f h 一f 一 图3 3 流场示意图 代入到方程( 3 7 ) ，( 3 1 1 ) ，( 3 1 ) 及( 3 1 2 ) 中，并以原始变量表示： v 2 击= 0( 3 - 1 3 ) v 2 p - - a s i n h ( ) ( 3 - 1 4 ) v v = 0( 3 - 1 5 ) 詈+ ( v v ) v + 堆i n h ( q 0 州= 去v 2 v ( 3 - 1 6 ) 式中，口，卢是无量纲参数，定义如下， 蹦2 哪= 等 m 黼t = ( 警) i ，2 ，揣。1 表示。却e 厚她相当于电渗流场的边界层 厚度。 3 2 3 涡量一流函数方法 在不可压流动问题中，虽然方程形式比可压缩情况下的简单，然而其求解过 程并不比后者容易。因为在不可压缩流场的控制方程中，压力梯度是以源项形式 出现在动量方程中的，压力没有独立的方程，而在可压缩流体中压力与密度之间 的关系可以通过状态方程规定。因此，对于不可压问题，求解中最大的困难就是 浙江大学硕士学位论文聂德明2 0 0 4 2 连续性方程v v = 0 的满足。目前二维不可压问题中的三种解法，涡量流函数法、 原始变量法及压力校正法，其立足点都与v v = 0 有关。由于本文课题所关心的 问题并不在于计算区域中的压力分布，而主要在于其中的速度分布，基于此考虑， 本文将采用涡量流函数求解流场，避开压力p o s s i o n 方程【i l 】【5 1 】【5 2 】。 将方程( 3 1 6 ) 两边取旋，得： 詈+ ( v v ) q 一( n v ) v 一历( s i n h ( 矿) ) v = 去v 2 q ( 3 _ 1 8 ) 式中，q = v v 为涡量，定义流函数v , ，v = v v ，显然连续性方程( 3 1 5 ) 自动满足。将叩带入到涡量n = v x v 中，可以得到涡量- 流函数方程 f l ：一v 2 叩( 3 - 1 9 ) 考虑二维、定常情况，则( 3 1 8 ) 简化为： ( v 审) q 一刃( s i n h ( 妒) ) v 2 去v 2 q ( 3 _ 2 0 ) 为了便于后面的数值计算，现在来推导贴体坐标下的方程形式。设x 。表示 笛卡儿坐标，它和贴体坐标一的关系为【5 3 】 x 。= x j ( x 1 ，x 2 ，x 3 ) ( 3 2 1 ) 则笛卡儿坐标x ”和贴体坐标之间的正变换系数屈。可由下式求得， 尼7 = 等 ( 3 2 2 ) 求得正变换系数属7 后，根据下式可求出逆变换系数屈。， 砌叫牡彤捌 ：。， 利用正变化系数属7 和逆变换系数屏。，笛卡儿坐标下的任意矢量分量q 。和 贴体坐标下的逆变矢量分量q 之间就可以通过如下式子相互转化： d = q ”展q 。= d p ，” ( 3 - 2 4 ) 协变度量张量g l ，和逆变度量张量可通过下式计算： 岛= 屈”局g “g p = t = ；三：) c s z s ， 贴体坐标系的第一类c h r i s t o f f e l 符号r 和第二类c h r i s t o f f e l 符号r m 可通过 如下式子计算： 浙江大学硕士学位论文 聂德明2 0 0 4 2 驴甄等+ 等+ 割f ，k - - f q m g r k z s ， 圳，= 等刷_ ( 3 - 2 7 ) 出。 。 及二阶协变倒数d k ： q ：掣删，w 删_ ( 3 _ 2 8 ) 。利用以上笛卡儿坐标系与贴体坐标系之间的转化关系，可以将描述电渗流场 的所有方程包括( 3 - 1 3 ) ，( 3 - 1 4 ) ，( 3 1 9 ) 及( 3 2 0 ) 分别化为便于差分计算张 量形式： 昙( 厨别= 。( 3 - 2 9 ) 去昙( 西9 等 - 戚血( 伊) c ，枷， 击鲁0 x ( 厄“詈) 一 b ， g l ”。撕 去去昙( 玉9 剽圳嚣一肋出c 詈等扩c ，啦， 式中，i = 1 ，2 ，、压= d e t ( g , 1 。方程( 3 中占妒是置换张量，定义如下：, 3 2 ) 占社= f ，女正转 忑 。7 蚍转 一f ，t 逆转 一忑 腿转 0 f ，j ，女有两个相同指标 3 2 4 边界条件 3 2 3 中已经推导出了描述电渗流场的数学方程，但是为了求解该问题，必 须给出相应的初始条件和边界条件。由于是定常问题，所以只存在边界条件。以 图3 4 所示的计算区域为例，下面对电势、流函数及涡量的边界条件逐一进行阐 浙江大学硕士学位论文聂德明 2 0 0 4 2 述。 图3 4 直槽道边界条件 对于外加电压场，壁面采用绝缘条件【4 5 】，即假定计算区域与外界绝缘；对 于进出口条件，则可以直接指定值，表示为： 掣l ：o 儿：破儿：欢 (3-33)o n b r ”“ 对于壁面电荷电势伊，壁面上等于f 电势4 5 1 ，进出口采用充分发展条件，即 札= f 剖。= 。剖。= 。 ( 3 - 3 4 ) 对于流函数v ，在固壁上通常可以设为常数，例如图3 4 所示槽道，下壁面 可以设为0 ，上壁面设为l 。进口条件可以取一个已知的分布，最简单的就是线 性分布，出口可以认为已经稳定，应用充分发展条件【5 l 】。如下式表示： y i 。= 0y i ，- - 1y i ，。= 七y 觌= 涡量q 的边界条件比较复杂，凼为它不能像速度、温度等从物理上提出在边 界上涡量为多少，或涡量导数为多少。因此涡量边界条件是一种数值边界条件。 对于壁面边界b ，利用t a y l o r 展开法将+ 在上展开川： 一+新+锐i11 hn ll c3-360y ， 弘7 b + = 弘0 + = l l+ = 寺i = + l () 】。卵一1 2 h 为网格间距。略去高阶量，且有譬i = ，则： 纠：吾( + 。一一。 ) ( 3 3 7 ) 矿l 。2 万叫8 叫p ) q 弓” 根据涡量的定义： 浙江大学硕士学位论文聂德明2 0 0 4 2 q 。= 一( 锐+ 洲 b s s ， 又因为沿b 界面有y = o ，= 0 ，所以： q s 一紊t ( 3 - 3 9 )， 对于上壁面边界t 可以采用同样的方法得到边界条件， 皿一寺( 一) ( 3 - 4 0 )， 由于进口处的流函数已经给定了，所以可以根据涡量流函数关系式确定进口处 的涡量。至于出口处的涡量，仍可以通过充分发展的条件来获得，即 0( 3 4 1 ) 以上边界条件是针对图3 4 所示的计算区域进行设置的，对于其他计算区域 如1 8 0 度槽道也可以作类似处理。 3 3 离散网格及方程 本章拟对直槽道及1 8 0 度弯槽道中的电渗流场进行模拟。对于直槽道采用均 匀正交网格。1 8 0 度常半径槽道的网格如图3 5 所示，图( a ) 表示了正交的物理网 格，可以直接生成，图( b ) 表示了计算平面上的均匀正交网格，为了计算方便将 间距都设为1 。图3 6 表示了变形弯槽道的网格划分，需要求解微分方程来生成 网格洲，为非正交网格。同样，计算过程在图( b ) 所示的网格内进行。求解微分 方程如下， 翥+ 两2 x l _ p ( 一蚂 ( 3 。z ) 叙1 苏1 蠡2 知2 、。 嘉+ 嚣= o ( x i , x 2 ， 。， 缸”苏。缸2 。 其中一，x 2 为物理平面上的坐标，通常边界点已知：x i , x 2 为贴体坐标，通常定义 为间距等于1 的矩形区域。p ，o 用来控制网格密度。 钾 塑( 毽 浙江大学硕士学位论文 聂德明 2 0 0 4 2 ( a ) 常半径槽道物理网格 ( b ) 常半径槽道计算网格 图3 51 8 0 度常半径槽道网格划分 ( a ) 变形槽道物理网格( b ) 变形槽道计算网格 图3 6 变形槽道网格划分 对于外加电压场方程( 3 2 9 ) ，展开得： 口窘+ a 窘+ c 等+ d 善+ e 嚣= 。 俘。， 式中，各系数口，b ，c ，吐e 定义如下 3 0 浙江大学硕士学位论文 聂德明2 0 0 4 2 口= 厄“6 = 西2 2 一a ( 斛) a ( 玉”) c = i 一+ i ：二 aa a ( 西“) a ( 西2 2 ) 矿+ 矿 p = 2 9 9 1 2 采用中心差分格式离散( 3 4 4 ) 得： v d u 谚，= 6 “谚“+ c “痧一i ，+ d w 一，+ j + p “谚。，一j + s j ， ( 3 4 5 ) 式中各系数为 d u = 2 ( q ，+ 6 ，) b | t 。i = d j + i 1c j 。l 。“2 a i ，厂j q j d u = 6 l ，j + i id l ， 。= b 。一；d l j s 。= 五1e u ( 声i + l j + l - - 唬扎，一一死。，。+ 谚 卜。) 其中s 。是当作源项来处理，值采用上一次迭代的值求出各系数后，离散方 程( 3 4 5 ) 就可以进行迭代求解了。本文采用t d m a 算法将离散方程分解成两 个方向上的三对角方程，分别通过追赶法进行求解。 对于壁面电荷电势方程( 3 - 3 0 ) ，其展开及离散过程与方程( 3 - 2 9 ) 的过程 类似，但是前者方程右端存在非线性项s i i l 1 ( 妒) ，难以进行差分离散，因此有必 要对其进行线性化【4 9 】： s i i l l l ( 伊) = s i l l l ( p ) + 妒c o s h ( p ) 一妒+ c o s h ( 妒) ( 3 4 6 ) 式中，妒表示计算过程中上一次迭代的结果。将上式代入到方程( 3 - 3 0 ) 中，就 可以进行差分离散了。 涡量流函数方程( 3 - 3 1 ) 展开后得： 浙江大学硕士学位论文聂德明2 0 0 4 2 盯争+ 。器+ c 詈+ d 詈+ e 嘉+ 徊= 。( 3 - 4 7 ) 各系数日，b ，c ，西p 同上，厂= g 。 涡量方程( 3 3 2 ) 展开后得， 口窘“窘“等材詈“嘉+ 毋= 。( 3 - 4 8 ) 式中各系数为 d = ab = b一= e c = c r e g 矿1 d = d - r e 以矿2 毋= r e