复杂微通道内非混相驱替过程的格子Boltzm ann方法∗

臧晨强 娄钦

(上海理工大学能源与动力工程学院,上海 200093)

复杂微通道内非混相驱替过程的格子Boltzm ann方法∗

臧晨强 娄钦†

(上海理工大学能源与动力工程学院,上海 200093)

(2017年1月10日收到;2017年5月4日收到修改稿)

本文采用改进的基于伪势模型的格子Boltzmann方法研究复杂微通道内的非混相驱替问题.这种方法克服了原始伪势模型中计算结果对网格步长的依赖.首先用Lap lace定律验证模型的正确性,然后用该方法研究壁面润湿性、粗糙结构、黏性比以及距离对非混相驱替过程的影响.模拟结果表明：与壁面粗糙结构和黏性比相比,壁面润湿性的影响是决定性的因素.随着接触角的增加,驱替效率增加,当接触角大于某一值后,驱替效率不再变化;随着黏性比的增加,驱替效率增加;而壁面粗糙性对驱替过程的影响较复杂,只有凸起半圆的半径在一定范围内增加时,驱替效率增加;距离较小时将促进驱替过程.

非混相驱替,格子Boltzmann方法,驱替效率,驱替时间

1 引 言

非混相驱替过程是一个重要的研究课题,广泛存在于日常生活和工业生产中,例如土壤中的非水相液体的运输、聚合物电解质膜燃料电池的气体扩散层中水的输运、石油的采集等.该问题是一个典型的接触线问题,涉及复杂的流体-流体间相互作用力和流体-固体间相互作用力.因此,尽管很多学者关注到了这类问题[1-21],但是目前的研究尚不完善.

一些研究者采用宏观方法和实验方法研究非混相驱替问题.Tecklenburg等[1]采用双连续模型研究水平裂隙岩体中非混相驱替过程.Zhu等[2]采用流体体积(VOF)方法研究聚合物电解质燃料电池气体通道中液态水的驱替过程,他们研究了壁面的润湿性、空气入口速度、注水速度和孔隙的大小对驱替过程的影响.Yang等[3]用实验的方法研究了毛细管力驱动的液-液驱替问题,重点研究了相对黏度对驱替效率的影响.Islam et等[4]研究了亲水性颗粒悬浮液中水的驱替问题.在他们的研究中,主要考察不同体积的水滴对驱替过程的影响.李维仲等[5]将实验方法和格子Boltzmann方法结合,研究了单个气泡在具有三个半圆形喉部的复杂流道内的上升过程.

尽管用宏观方法和实验方法研究非混相驱替过程可以得到流动的宏观参数,如驱替效率等,但得不到流动的细节.为了弥补该不足,一些研究者用微观方法研究驱替问题[6,7].Primkulov[6]利用分子动力学理论研究了微米级球形玻璃表面水滴的驱替问题.在他们的研究中,考察了毛细数、流体黏度和界面特性对驱替过程的影响.数值计算结果表明,分子动力学模型很好地拟合了毛细数大于10-4的实验数据.利用分子动力学方法,Kop lik等[7]研究了低雷诺数下的非混相驱替过程.研究表明,在低雷诺数的情况下,无滑移条件在接触线附近不适用.虽然流动的细节可以通过微观方法获得,但由于微观方法的计算量较大,因此仅适用于较小的计算区域.

近年来,作为连接微观方法和宏观方法桥梁的介观方法开始被诸多学者采用,格子Boltzmann方法作为介观方法的一种开始受到广泛关注.该方法的优点是其有着宏观方法的本质和微观方法的特性,因此它能够直接处理流体与流体间的作用力及流体与固壁间的作用力.目前基于格子Boltzm ann方法的非混相驱替问题的研究越来越多.Kang等[12,13]在二维和三维通道里研究了不混溶液滴受重力作用下的驱替流动,包括壁面润湿性、Bond数、液滴大小和密度比对驱替流体的影响.之后,他们又研究了驱替过程中的指进现象并得到了一些重要的结果[14].然而他们的研究主要集中在没有孔隙介质的单通道的驱替问题.Huang等[15]采用基于颜色梯度的多相格子Boltzmann方法,研究了毛细数和黏性比对非均相多孔介质中非混相驱替过程的影响.Dong等[16,17]采用格子Boltzmann方法研究了微通道和多孔介质中毛细数、Bond数和黏性比对非混相驱替过程的影响.李维仲等[18]采用格子Boltzm ann方法研究了水平通道内流体黏度比以及壁面润湿性对非混相驱替过程的影响.李娟等[19]采用多相多组分格子Boltzm ann模型中的伪势模型研究了单通道中壁面上液滴驱替过程.彭本利等[20]采用自由能格子Boltzmann方法研究了在不同蒸汽速度剪切作用下,液滴在具有不同润湿性固体表面上的变形和运动过程.以上研究推动并揭示了一些驱替机理,但对于一些典型的多孔结构,如孔喉结构,其流动细节却不清楚.为了得到孔喉结构中的非混相驱替问题的详细信息,Liang等[21]利用数值模拟方法研究了含腔微通道内的非混相驱替问题.在该工作中,主要研究了壁面润湿性、毛细数和密度比对驱替效率的影响,然而没有考虑微通道的粗糙性.本文在Liang等的工作基础上,研究了含粗糙壁面的微通道内的驱替问题,其中微通道的粗糙结构由一个凸起的半圆来体现,在研究中采用了格子Boltzmann方法,并主要考察了壁面润湿性、粗糙度、黏性比和粗糙表面与含液相腔间的距离对非混相驱替过程的影响.

目前有很多常用的多相格子Boltzmann模型被相继提出,例如颜色梯度模型[22]、伪势模型[23,24]、自由能模型[25,26]、动力学模型[27-29].在这些方法中,伪势模型由于其简单性成为最常用的方法之一.然而它也有一些局限性,例如Yu等[30]发现流体性质对网格步长有依赖性,并提出了一个改进的格子Boltzmann方法来克服原模型的缺陷.本文采用了这种改进的格子Boltzmann方法来研究复杂通道内的非混相驱替问题.

2 数值方法

本文采用了Yu等[30]提出的改进的伪势模型进行数值模拟研究.这种模型具有简单性和可靠性.由于该模型在文献[30]中已经被详细描述,因此在这里我们只给出一个简短的描述.在这个模型中粒子的分布函数fi的演化方程为

其中i=0,1,2,···,b-1,b是离散速度的个数;x和t分别表示是位置和时间;ci表示离散速度;τ是松弛时间,它与运动黏度ν的关系为(δt表示时间步长),是一个模型常数,其中c=δx/δt(δx代表网格步长).是平衡态分布函数,且与局部密度ρ和速度u有关：

其中,wi是权重系数.方程(1)中的Fi是力项,它的表达式为[31]

其中F是总的相互作用力,它包含三部分：流体间相互作用力(Fint)、流固作用力(Fads)和外力(G).在本文使用的模型中,导致相分离的流体间作用力F in t为[24,30]

其中g是一个常数,表示流体间作用力强度.ψ(x)是“有效质量”,它是与局部密度有关的一个函数：

流固作用力Fads的表达式为[32]

其中ρw是壁面密度;s(x)是一个开关函数,当x是固体点时,s(x)取值为1,而当x为流体点时,s(x)取值为0.宏观物理量的表达式为：

本文使用D2Q9模型进行数值模拟研究.其中权重系数wi的参数设置如下：w0=4/9;当i=1-4时,wi=1/9;当i=5-8时,wi=1/36. c i为

3 模型验证

与其他研究者的模型验证方法一样[21,33-36],本文用Laplace定律验证程序的正确性.初始时,在N x×N y的格子区域中心放置一个半径为r、密度为ρl的静止液滴,其余区域充满着蒸汽,密度为ρv.当系统达到平衡状态时,液滴呈圆形.根据Laplace定律,液滴内、外压力差(Δp)与表面张力σ和液滴半径r满足关系式

本文模拟不同的液滴半径下的系统状态以验证Lap lace定律,半径的变化为r=18,20,22,24, 26,28,30,32,35,40,45,50,并考虑了不同黏性比(M)下的系统状态,以确保可信度,其取值分别为M=3.78,6.19,7.57,9.10.表1列出了由Maxwell重构得到的具体参数值.在以上所有的模拟中N x×N y=128×128,x方向和y方向均为周期边界条件.图1给出了不同黏性比M 下的液滴内外压力差(pi-po)和液滴半径1/r的关系.从图中可以得到数值模拟结果和Lap lace定律很符合.同时,根据方程(11)可知,图中直线的斜率即为对应的表面张力,在本文中黏性比M=3.78,6.19,7.57,9.10下的表面张力分别为 0.0176,0.0445,0.0605和0.0782.

图1 不同黏性比M下的(pi-po)和1/r的关系Fig.1.Relationship between(pi-po)and 1/r at different viscosity ratios M.

表1 不同黏性比下的参数值Tab le 1.Param eter values of different viscosity ratios.

4 模拟结果与分析

图2给出了微通道结构示意图,灰色区域代表宽度为R的固体,由于其宏观参数不发生变化因此不参与计算;白色区域为流体流动的区域,数值模拟时只需要模拟白色区域.模型包含一个圆心O1半径为R1的凸起半圆(obstacle A),和一个圆心O2半径为R2润湿性为θ的半圆腔(cavity B),两者之间的距离为D,凸起A代表壁面的粗糙度.初始时刻,密度为ρl的液体放在腔B中,其他计算区域为密度为ρv的气体.气体从左侧进口流入,右侧出口流出.在模拟中,整个系统的网格为N×L,其中N=R+W是竖直方向上的网格数,L为水平方向上的网格数,在流体流动的x方向上还施加有一个质量力G.不失一般性,本文忽略了重力的影响,因为重力在垂直方向上的作用很小,而且流体是在毛细力和黏性力的作用下流动的[21].固体壁面上选用无滑移边界条件,进出口选用周期性边界条件以保证质量守恒.

在驱替过程中,流体流动的无量纲参数设置如下：毛细数Ca=Uρvν/σ,其中U是在质量力作用下流体的最大速度,表达式为U=W2G/(8ν).在模拟中,主要是通过两个参数衡量驱替效果：驱替效率De和驱替时间t.其中De是被驱出腔B的液体质量与原始腔B液体总质量的比值;t是液体从初始时刻开始到被驱到微通道出口的时间,这里取W/U为特征时间,在文中出现的时间为无量纲时间.

由于不同的接触角对应不同的润湿性,在这里计算了接触角以满足后续研究的需要.润湿性有三种形式：亲水性、疏水性、和中性.当液滴在材料表面接触角大于90°时,材料表面具有疏水性;当液滴在材料表面接触角小于90°时,材料表面具有亲水性;当液滴在材料表面接触角等于90°时,材料表面具有中性.在本文使用的模型中,通过改变壁面密度ρw来控制接触角的大小[37,38].

我们模拟了一个液滴在壁面上的展开情况来得到壁面密度和接触角的关系.初始时,一个半径为r密度为ρl的液滴放在一个矩形计算区域底部的中心,其他区域充满着密度为ρv的气体.在所有的模拟中,计算区域网格数为Nx×Ny=150×300.液体密度ρl=3.0404,气体密度ρv=0.3342, g=-10.5.x方向为周期性边界条件,y方向为无滑移边界条件.图3给出了不同壁面密度(ρw)下分别对应的接触角(θ).从图中可以看出当壁面密度从ρl减小到ρv时,对应的接触角从0°增加到180°.

图2 微通道结构示意图Fig.2.Schem atic illustration of the channel structure.

图3 (网刊彩色)不同壁面密度下的接触角Fig.3.(color on line)Contact angle at d iff erent su rface densities.

4.1 润湿性的影响

固体壁面的润湿性是驱替过程的关键因素之一,本节研究了润湿性对驱替过程的影响.众所周知,固体壁面润湿性可以量化为气液固三相的接触角,因此本节研究了不同的接触角θ对驱替过程的影响.在研究中,分别考虑了以下接触角θ=20°, 40°,50°,60°,75°,90°,105°,130°和160°.在数值模拟中,其他参数设置如下：N×L=100×1000, ρl=3.0404,ρv=0.3342,g=-10.5,ν=0.1667, R=W=R2=2R1=2D=50,Ca=0.05,圆心O2的坐标为(200,50).

表2给出了驱替过程在不同接触角下的驱替效率(D e)和驱替时间(t).从表中可以清晰地得到壁面的润湿性对驱替过程有很重要的影响：当接触角小于某一角度时(这里是75°),驱替效率D e随着θ的增加而增加,当接触角增加到某值后,De不再变化;当接触角很小时,例如θ=20°,驱替效率D e为0.6552,在这种情况下腔B液体有大部分未被驱出,随着接触角的增加,驱替效率增加;当θ=60°, De增加到0.8951,这说明腔B中液体大部分已经被驱出;当θ大于75°时,D e等于1.0,即腔B液体全部被驱出.上述数值模拟的结果与理论预测的一致.事实上,固体壁面的润湿性受黏附力(adhesive forces)和凝聚力(cohesive forces)的影响.接触角较小对应着较大的黏附力,液体不容易与壁面分离,因此较难从腔B中被驱出;随着接触角的增加,黏附力变小,部分液体容易被驱出腔B,因而De增加;当凝聚力大于黏性力时,腔B中的液体全部被驱出.从表2还可以看出,接触角越大,驱替时间(t)越短.综上所述,增加接触角对驱替过程有利.

为了更直观地揭示接触角对驱替过程的影响,图4给出了不同接触角下的瞬时流动情况.当驱替过程在接触角为20°和60°时,初始时一部分液滴流出腔B,并铺展在壁面上,这部分液体上受到来自进口流体的拖拽力,使得液体在壁面上缓慢地向出口流动.同时在固体壁面上的液体和腔内剩余液体之间的连接处,表面张力和拖拽力同时施加在此处液体上.随着越来越多的液体流出,拖拽力会克服表面张力,使得液体在此处断裂.随后壁面上的液体流向出口,腔内液体依然留在腔内.同时接触角20°驱替过程的驱替时间大于接触角为60°驱替过程的驱替时间.另外,接触角大于75°时会略有不同,主要体现在液体不会发生断裂,并全部被驱出微通道,这是由于接触角增大,腔B中液体受到的黏附力变小,因此作用在液体上的驱力更容易将液体驱出.

表2 不同接触角下的驱替效率和时间Tab le 2.D isp lacem ent efficiency and times at d iff erent contact angle.

图4 (网刊彩色)不同接触角下的流动过程 (a)θ=20°,从左到右t=0,20,40,50;(b)θ=60°,从左到右t=0,20,34,44;(c)θ=90°,从左到右t=0,20,24,36Fig.4.(color on line)Flow process at d iff erent contact angle:(a)θ=20°,from left to right t=0,20,40, 50;(b)θ=60°,from left to right t=0,20,34,44;(c)θ=90°,from left to right t=0,20,24,36.

4.2 壁面粗糙度的影响

粗糙的壁面将影响流场的变化,因而本节研究壁面粗糙度对驱替过程的影响.在本节中,通过改变凸起A的尺寸(R1),研究不同粗糙度的壁面对驱替过程的影响.其中凸起A半径的参数设置如下：R1=0.0R2,0.1R2,0.2R2,0.3R2,0.4R2,0.5R2, 0.6R2,0.7R2,0.8R2和0.9R2.其他参数和边界条件与4.1节中设置相同.在上一节中的研究已经说明固体壁面的润湿性对驱替过程有重要的影响,因此为了模拟结果的一般性,在这里也考虑了不同接触角下(θ=40°,50°,60°,75°和90°)的驱替过程.

不同尺寸的凸起A在不同接触角下对驱替效率(D e)和驱替时间(t)的影响如图5所示,其中对应θ=40°和θ=50°的垂直向上的折线代表t=∞.从图中可以看出凸起A对驱替过程有着不可忽略的影响,并且R1对驱替过程的影响依赖于壁面润湿性(θ).当接触角(θ)小于90°并大于60°时,D e随着R1的增加而增加,当De增加到1.0后便保持不变;然而当接触角等于或小于50°时,D e首先随着R1的增加而增加,当R1大于某一值后De急剧减小.当θ大于等于90°,由于效率均为1.0,曲线重合因而未给出.情况下所有的液体都可以从腔B中驱出,而且驱替时间(t)几乎不随着凸起A的尺寸(R1)改变而改变.这说明对于液体可以全部被驱出的情况,凸起A的尺寸R1对驱替过程的影响并不大,从这个意义上说,与粗糙度的大小相比,润湿性对驱替过程影响更大.相应地,t的变化也是一个复杂的过程.对于所有的接触角来说,当R1小于0.7R2时,t随着R1的增加而减小;当R1大于0.7R2时,t随着R1的增加而增加,液体不能被驱出腔B,驱替时间无穷大.这种现象出现的原因可以从驱替过程的流场分析,如图6所示,随着R1的增加,腔B中x方向的速度减小,y方向的速度增加,这使得液体在腔B中更容易被驱出：驱替效率增加,驱替时间减小.但当R1大于某一值后,通过凸起的速度会变得很小,因此会使驱替过程变难：驱替时间增加(如θ=60°和θ=75°的情况);或者液体不能从腔B中驱出,使得驱替时间无穷大(如θ=50°和θ=40°的情况).

图5 不同尺寸的凸起A在不同接触角下对驱替效率和驱替时间的影响Fig.5.Dependency of D e and t on the d iff erent size of obstacle A at different contact angle.

图6 接触角为50°时的流动结果:从上到下:R1/R2= 0.1,0.6,0.8Fig.6.F low resu lts of d iff erent sizes of obstacle at contact angle 50°:from top to bottom:R1/R2=0.1, 0.6,0.8.

4.3 黏性比的影响

在本节中我们研究驱替过程中另一个重要的参数黏性比(M).研究中考虑了七组不同的黏性比,分别是M=3.78,4.23,4.94,5.68,6.19,7.57, 9.10(具体的参数如表1所列).Ca=0.1,其他参数值和边界条件与上节相同.

图7给出了不同黏性比(M)对驱替效率(D e)和驱替时间(t)的影响,图中虚线代表t=∞.从图中可以很清楚地看出：如图中虚线所示,在小黏性比时,液体不能被驱出微通道,驱替时间无穷大.随着黏性比(M)的增加,驱替效率De增加,驱替时间t减小,即在黏性比较大时,腔B液体驱出质量较多并且驱替过程更快.为了使结果更直观,图8展示了不同黏性比(M)下的驱替结果.如图所示,小黏性比的情况下,液体不能从腔B中被驱出;增加黏性比会使得腔B中液体被驱出,同时被驱出液体的质量随着黏性比的增加而增加.从液体铺在壁面上的长度可以看出,黏性比越大,长度越长,说明大黏性比下的液体受到的驱力较大.这也同时解释了黏性比(M)越大,驱替效率(D e)越高,驱替时间(t)越小.

图7 不同黏性比(M)对驱替效率(D e)和驱替时间(t)的影响Fig.7.Dependency of D e and t on the d iff erent viscosity ratio(M).

图8 (网刊彩色)不同黏性比(M)下的驱替结果Fig.8.(color on line)Disp lacem ent results at different density ratios(M).

4.4 距离的影响

本节研究凸起A和腔B间的距离(D)对驱替过程的影响. 在数值模拟中D分别取 0.0R2, 0.25R2,0.5R2,0.75R2,1.0R2,1.5R2,2.0R2,2.5R2, 3.0R2,3.5R2,4.0R2,4.5R2,5.0R2,6.0R2,7.0R2,其中R2=50.数值计算中用到的其他参数如下：腔B的圆心为(500,50),R1=0.2R2;Ca=0.08; ρl=3.0404,ρv=0.3342,g=-10.5;θ=60°.

图9给出了在不同距离(D)下的驱替效率(D e)和驱替时间(t).如图所示,当D=0R2即凸起A和腔B的距离最小时,驱替效率最大,此时驱替时间也最短.随着距离D的增大,驱替效率D e减小,同时驱替时间t增加.当D增加到0.75R2时驱替效率达到最小值,对应的驱替时间t达到最大值.然后随着D继续增大,驱替效率增大继而保持不变,对应的驱替时间减小并保持不变.进一步研究发现,最后不再变化的驱替效率和驱替时间对应没有凸起结构时的驱替效率和驱替时间.此模拟结果说明有凸起结构的情况与没有凸起结构的情况相比,当凸起结构和腔B之间的距离较小时该凸起结构对驱替结果是有利的,当凸起结构和腔B之间的距离很大时,凸起结构对驱替没有影响,当凸起结构和腔B之间的距离在两个极值之间时,凸起结构将对驱替有阻碍作用.

图9 不同距离(D)对驱替效率(D e)和驱替时间(t)的影响Fig.9.Dependency of D e and t on the d iff erent distances(D).

5 结 论

本文采用改进的伪势格子Boltzmann方法研究复杂微通道内的两相流驱替问题,主要研究了壁面润湿性、粗糙度、黏性比以及距离对驱替过程的影响.研究表明：1)随着接触角的增加,驱替效率增加,驱替时间减小;当接触角增加使得驱替效率等于1.0以后,驱替效率保持不变,驱替时间仍然随着接触角的增加而减小;2)壁面的粗糙度对驱替过程也有很重要的影响,本文用一个凸起的半圆凸起A代表壁面的粗糙性,其半径的大小代表粗糙度的大小,数值模拟结果表明,随着凸起半径A的增加,驱替效率增加,驱替时间减小;然而当其半径增加到某一值后,腔B液体不能被驱出;3)黏性比同样影响驱替过程,增加流体的黏性比会使得驱替效率增加,驱替时间减小;4)距离D为0时驱替效率最大,随着D的增加,驱替效率减小然后增加最后趋于常数.

[1]Teck lenburg J,Neuweiler I,Dentz M,Carrera J,Geiger S,Ab ram ow ski C,Silva O 2013 Adv.Water Res.62 475

[2]Zhu X,Sui P C,D jilali N 2008 J.Power Sources 181 101

[3]Yang D,K rasowska M,Priest C,Ralston J 2014 Phys. Chem.Chem.Phys.16 24473

[4]Islam S F,Sundara R V,W hitehouse S,Althaus T O, Palzer S,Hounslow M J,Salm an A D 2016 Chem.Eng. Res.Des.110 160

[5]Li W Z,Sun H M,Dong B 2013 Chin.J.Com putat. M ech.1 106(in Chinese)[李维仲,孙红梅,董波2013计算力学学报1 106]

[6]Prim ku lov B K,Lin F,Xu Z 2016 Co lloids Surf.A: Physicochem.Eng.Aspects 497 336

[7]Kop lik J,Banavar J R,W illem sen J F 1988 Phys.Rev. Lett.60 1282

[8]Zhou G,Chen Z,Ge W,Li J 2010 Chem.Eng.Sci.65 3363

[9]Jam aloei B Y,K harrat R 2010 Transp.Porous M ed.81 1

[10]Pram anik S,M ishra M 2016 Phys.Rev.E 94 043106

[11]Yang K,Guo Z 2016 Com put.F luids 124 157

[12]K ang Q,Zhang D,Chen S 2002 Phys.F luids 14 3203

[13]Kang Q,Zhang D,Chen S 2005 J.Fluid M ech.545 41

[14]K ang Q,Zhang D,Chen S 2004 Adv.W ater Res.27 13

[15]Huang H,Huang J J,Lu X Y 2014 Com put.F luids 93 164

[16]Dong B,Yan Y Y,LiW,Song Y 2010 Com put.F luids 39 768

[17]Dong B,Yan Y Y,LiW Z 2011 Transp.Porous.M ed. 88 293

[18]LiW Z,Dong B,Song Y C 2012 J.Dalian Univ.Techno l.3 343(in Chinese)[李维仲,董波,宋永臣2012大连理工大学学报3 343]

[19]Li J,Song Y C,LiW Z 2009 J.Therm al Sci.Technol. 4 284(in Chinese)[李娟,宋永臣,李维仲2009热科学与技术4 284]

[20]Peng B L,Xu W,W en R F,Lan Z,Bai T,M a X H 2015 J.Eng.Therm ophys.4 820(in Chinese)[彭本利,徐威,温荣福,兰忠,白涛,马学虎2015工程热物理学报4 820]

[21]Liang H,Chai Z,Shi B,Guo Z,Li Q 2015 In t.J.M od. Phys.C 26 1550074

[22]Gunstensen A K,Rothm an D H,Zaleski S,Zanetti G 1991 Phys.Rev.A 43 4320

[23]Shan X,Chen H 1993 Phys.Rev.E 47 1815

[24]Shan X,Chen H 1994 Phys.Rev.E 49 2941

[25]Sw ift M R,Osborn W R,Yeom ans JM 1995 Phys.Rev. Lett.75 830

[26]Sw ift M R,O rland ini E,Osborn W R,Yeom ans J M 1996 Phys.Rev.E 54 5041

[27]Luo L S 1998 Phys.Rev.Lett.81 1618

[28]He X,Chen S,Zhang R 1999 J.Com put.Phys.152 642

[29]Guo Z,Zhao T S 2005 Phys.Rev.E 71 026701

[30]Yu Z,Fan L S 2009 J.Com put.Phys.228 6456

[31]Guo Z,Zheng C,Shi B 2002 Phys.Rev.E 65 046308

[32]M artys N S,Chen H 1996 Phys.Rev.E 53 743

[33]Zhang R,He X,Chen S 2000 Com pu t.Phys.Comm un. 129 121

[34]Fakhari A,Rahim ian M H 2009 Comm un.Non linear Sci.Num er.Sim u lat.14 3046

[35]Zheng H W,Shu C,Chew Y T 2006 J.Com put.Phys. 218 353

[36]Hao L,Cheng P 2000 J.Power Sources 190 435

[37]Huang H,Lu X 2009 Phys.Fluids 21 092104

[38]Li Q,Luo K H,Kang Q J,Chen Q 2014 Phys.Rev.E 90 053301

(Received 10 January 2017;revised manuscript received 4 May 2017)

Lattice Boltzmann simulation of immiscible displacement in the complex micro-channel∗

Zang Chen-Qiang Lou Qin†

(School of Energy and Power Engineering,University of Shanghai for Science and Technology,Shanghai 200093,China)

The imm iscible disp lacem ent process inmicro-channel,which w idely existes in daily life and industrialp roduction,is an im portant research sub ject.This sub ject is a typical contact line prob lem involving com p licated fl uid-fluid interactions and fl uid-solid interactions which have attracted the interest ofmany scholars.Although the imm iscib le disp lacement in micro-channels has been studied by some researches,the problem is still not fully understood because them echanism of the imm iscib le disp lacement is very com p lex.In order to further explain the physical mechanism of imm iscible disp lacement process in micro-channels,detailed numerical simu lations are carried out in a com p lex micro-channel containinga semicircular cavity and a sem icircular by bulge using an im p roved pseudo-potential lattice Boltzm ann method(LBM).Thismodel overcomes the drawback of the dependence of the fl uid properties on the grid size,which exists in the original pseudo-potential LBM.Initially,the cavity is fi lled with the liquid and the rest of the area is fi lled with its vapour.The sem icircular bulge rep resents the roughness of the micro-channel.The approach is fi rst validated by the Lap lace law.The results show that the numerical resu lts are in good agreement with the theoretical predictions.Then themodel is em p loyed to study the imm iscible disp lacem ent process in themicro-channel.The effects of the surface wettability,the surface roughness,the viscosity ratio between the liquid phase and the gas phase,and the distance between the sem icircular cavity and the sem icircu lar bu lge are studied.The simulation resu lts show that the in fluence of the surface wettability on the disp lacem ent process is a decisive factor com pared with other factors. W ith the increase of the contact angle,the disp lacem ent efficiency increases and the disp lacem ent time decreases.W hen the contact angle is larger than a certain value,all of the liquid can be disp laced from the cavity.At that time,the disp lacement efficiency is equal to 1.The above resu lts are consistent with the theoretical prediction that with the increase of the contact angle,the liquid is easily driven out of the cavity because the adhesion force of the liquid in the cavity decreases.On the other hand,the influence of the surface roughness on the disp lacem ent process ism ore com p lex. The disp lacem ent efficiency increases with the radius of the sem icircle bulge increasing in a certain range.W hen the radius is larger than a certain value,the liquid cannot be ejected from the cavity due to the velocity around the cavity is too sm all.Furtherm ore,the liquid cannot be disp laced from the cavity at a sm all viscosity ratio.As the viscosity ratio increases,the disp lacem ent efficiency increases and the disp lacem ent time decreases.As for the distance between the sem icircular bulge and the sem icircular cavity,it promotes the disp lacement process at an early stage.W hen the distance exceeds a certain value,it has little effect on the disp lacem ent process.

imm iscible disp lacement,lattice Boltzmann method,disp lacement efficiency,disp lacem ent time

PACS:47.11.-j,47.55.Ca,47.56.+r DO I:10.7498/aps.66.134701

∗国家自然科学基金(批准号:51406120)资助的课题.

†通信作者.E-m ail:louqin560916@163.com

PACS:47.11.-j,47.55.Ca,47.56.+r DO I:10.7498/aps.66.134701

*Pro ject supported by the National Natural Science Foundation of China(G rant No.51406120).

†Corresponding author.E-m ail:louqin560916@163.com