黏弹性液滴热毛细迁移的对流不稳定1)

2021-05-31章绍能胡开鑫

力学学报 2021年5期

章绍能胡开鑫

(宁波大学机械工程与力学学院力学系,浙江宁波 315211)

引言

当液滴置于温度分布不均的固体壁面上,其表面张力梯度会驱动液滴产生热毛细迁移.这种现象广泛存在于微流控装置[1]、喷墨印刷[2]和材料加工[3]等多种工业应用中.为了对液滴在固壁面上的迁移过程进行精准控制,近几十年来学者们从理论分析[4]、数值模拟[5-7]和实验[8-10]等方面对液滴的热毛细迁移进行了广泛的研究.此外关于液滴的研究还包括液滴撞击固体表面[11]、通过液桥形成液滴[12]和液滴喷射[13]等.

为了得到液滴发生迁移所需要的条件,Ford 和Nadim[14]从多个方面对液滴在固体表面上的热毛细迁移速度进行了理论分析,并得出了液滴迁移速度方程.Smith[15]研究了二维液滴在施加温度梯度的固壁面上的迁移,并利用润滑理论推导了液滴形状方程.Chen 等[16]通过理论与实验两方面对基底上液滴的控制进行研究,得出了液滴迁移速度的影响因素,其结果与Ford 和Nadim[14]的研究理论一致.Brzoska等[17]通过实验研究了液滴在不可润湿固体表面对水平温度梯度的反应,发现只有液滴半径大于临界半径时液滴才会移动.Haj-Hariri 等[18]模拟了恒定温度梯度下可变形液滴的三维热毛细迁移,发现迁移速度与液滴形状有关.Pratap 等[19]通过理论和实验研究了癸烷液滴在PDMS 镀层表面上的热毛细迁移,实验得出了液滴迁移速度随液滴大小和温度梯度的变化关系,与理论预测结果一致.Dai 等[20-23]建立了石蜡油滴在单向和全向热梯度作用下迁移的理论模型并进行实验,还通过实验研究了球板之间液桥的热毛细迁移,并且从基本原理、评估方法和设计操作策略对界面液体流动进行综述.

在液滴热毛细迁移过程中,若雷诺数过大可能会发生流动失稳,使得内部形成湍流.如果液滴足够薄时还会发生破裂[24],甚至产生指纹式不稳定[25].此时的液滴迁移是不可控的,所以对液滴迁移进行流动稳定性分析是必要的.热毛细对流不稳定性的研究已具有多年的历史,Davis[26]以及Schatz 和Neitzel[27]已对此进行了综述.Smith 和Davis[28]采用线性稳定性理论研究了热毛细液层的对流不稳定,发现存在两种不稳定性,分别为静止的纵向转动和不稳定的水热波.他们的结论得到了实验[29]和数值模拟[30]的验证.Burelbach 等[31]对液膜因热毛细不稳定而发生破裂进行了详细讨论.Hu 等[32]对附壁液滴热毛细迁移的流动失稳进行了理论分析,结果表明:临界Marangoni 数总是随着液滴迁移速度的增大而减小.

目前对液滴热毛细迁移的研究主要集中在牛顿流体.然而许多工业应用中的液滴是聚合物流体,既有黏性又有弹性.黏弹性热毛细液层的研究结果表明:流体中的弹性决定了流动失稳的临界数、模态和机制[33].然而目前尚无研究关注弹性对液滴热毛细迁移的影响.本文采用线性稳定性理论研究了附壁黏弹性液滴在热毛细迁移中的对流不稳定性,得到了临界参数,画出了扰动流场图,并进行了能量分析.

1 控制方程

考虑一扁平液滴放置于固壁平面上,如图1 所示.平面水平方向上的温度梯度会导致液滴的表面张力梯度,进而驱动液滴热毛细迁移.其中x,y和z分别代表流向、展向和法向,坐标原点为液滴最左侧接触固壁面的位置.φ 为接触角,ζ 为迁移速度,令液滴厚度d为特征长度.重力作用下液滴毛细长度为κ-1=,ρ,g分别为表面张力、流体密度和重力加速度,液滴厚度和毛细长度的关系为d=2κ-1sin(φ/2)[34].假设液滴厚度d远小于宽度L,则液滴宽度远大于毛细长度,此时液滴可看成薄膜.所以假设液滴扁平是合理的,除边缘外液滴内的流动可近似看作为平行剪切流.

其中u,p,T分别为速度、压强和温度,τ为应力张量.重力因数G,Reynolds 数Re,Marangoni 数Ma和Prandtl 数Pr分别定义为

定义z=0 为固壁面,z=1 为自由面.设固壁面流向温度线性分布,液滴整体以速度ζ 向右迁移,并以液滴整体为参考系,则固壁面边界条件为

自由面边界条件为

分别表示热毛细效应和热流量.假设液滴基本流是平行的,温度在x方向的分布是线性的,则有

其中带有下标0 的变量表示基本流,Tb为垂直方向上的温度分布.在液滴整体参考系中,任意垂直截面上的质量流量为0,可推出基本流

因为黏弹性流体的弹性作用,在τ0中存在一正应力.式中无下标0 的变量表示扰动,σ=σr+iσi,σr和σi分别为增长率和频率,α,β 分别表示在x,y轴上的波数.总波数、波传播角和波速分别用k=θ=tan-1(β/α)和c=-σi/k表示.扰动方程和边界条件在附录A 中给出.σ 可以通过Chebyshev 配点法求解,配置点数Nc通常设置为70～100.表1 中对牛顿流体[32]与η/λ →1 时的Oldroyd-B 流体特征值进行比较,两者结果完全一致.表2 中给出了临界模态不同配置点下特征值的趋势情况.从表中可以看出,本文数值结果是足够准确的.

表1 Pr=0.01,ζ=0.2,Ma=15.5,k=1.49,θ=83.4° 时牛顿流体与Oldroyd-B 流体最不稳定特征值的比较Table 1 A comparison of most unstable eigenvalues for Newtonian fluid and Oldroyd-B fluid at Pr=0.01,ζ=0.2,Ma=15.5,k=1.49,θ=83.4°

表2 Pr=1,ζ=0.1,η/λ=0.1,ε=0.01, Ma=123.61,k=2.582,θ=46.5° 时不同配置点数下的中性模态特征值Table 2 Eigenvalues of neutral modes computed by different numbers of collocation points at Pr=1,ζ=0.1,η/λ=0.1,ε=0.01, Ma=123.61,k=2.582,θ=46.5°

2 数值结果

定义临界Marangoni 数Mac为所有波数下中性Marangoni 数MaN(σr=0)的最小值

第二个方面，益智区游戏材料的投放要根据主题开展，讲究层次性。在幼儿园时期的每个阶段都会有不同的主题设定，每一个主题都有不同的活动内容和重点、难点。所以，在进行游戏材料投放时一定要结合主题，适时投放，让幼儿可以在每个主题下集中的学习，掌握一方面的知识或者能力，得到更好的实践。在投放游戏材料的时候一定要讲究层次性，一方面，要按照由易到难，由浅入深的规律进行投放，做到与幼儿的认知发展相吻合；另一方面，要考虑到每个幼儿之间在能力、兴趣、爱好等方面的差异性。这需要教师根据幼儿的特点提供不同层次的游戏材料，实行因材施教的教育方式，根据幼儿的兴趣和能力进行投放，以此来促进每一个幼儿能力的提高。

为简便起见,只考虑ζ=0.1,η/λ=0.1 的情况,即迁移速度中等且溶剂黏度占主导.在计算中,分析Mac与弹性数ε=λ/Re=之间的关系.因为后者只与流体的性质和流动几何特性有关,而与特征速度U无关.

2.1 临界曲线

根据Prandtl 数的定义可知,Prandtl 数反映流体内热对流和热传导的相对重要性.Prandtl 数较小时失稳由惯性机制引起,Prandtl 数较大时失稳机制由热毛细机制引起.当Prandtl 数在1 附近时,热传导与热对流的重要性相当,此时流体失稳机制会变得更加复杂,所以本文选择了Prandtl 数为0.3,1 和3 进行计算.

图2～图5 为Pr=0.01 时Mac、波数、波传播角和波速随ε 的变化曲线.可以发现曲线存在3种临界模态,当ε <0.002 和ε >0.12 时(图2～图5 中的曲线a和d),临界模态为逆向斜波(θ>90°);0.002 <ε <0.035 时(图2～图5 中的曲线b),临界模态为同向流向波(θ=0°);0.035 <ε <0.12 时(图2～图5 中的曲线c),波的传播方向发生变化,此时的临界模态为逆向流向波(θ=180°).综合起来可以看出,曲线a的各临界值均随ε 的增大而略微增大;曲线b和c对应的流向波的Mac先随ε 的增大而减小,随后有短暂的上升,波数一直随ε 的增大而减小,波速先随ε 的增大而减小直到0,此时模态从同向流向波转为逆向,之后随ε 的增大而增大;曲线d的Mac随ε 的增大而减小,波传播角随ε 的增大而迅速减小,而波数的变化是微小的,波速随ε 的变化并不是单调的.

图2 Pr=0.01 时Mac 随ε 的变化曲线,其中a,d 为逆向斜波,b 为同向流向波,c 为逆向流向波Fig.2 Variation of Mac with ε at Pr=0.01,where a and d show upstream oblique waves,b shows downstream streamwise wave,c shows upstream streamwise wave

图3 Pr=0.01 时波数随ε 的变化曲线,其中a,d 为逆向斜波,b 为同向流向波,c 为逆向流向波Fig.3 Variation of wave number with ε at Pr=0.01,where a and d show upstream oblique waves,b shows downstream streamwise wave,c shows upstream streamwise wave

图4 Pr=0.01 时波传播角随ε 的变化曲线,其中a,d 为逆向斜波Fig.4 Variation of wave propagation angle with ε at Pr=0.01,where a and d show upstream oblique waves

图5 Pr=0.01 时波速随ε 的变化曲线,其中a,d 为逆向斜波,b 为同向流向波,c 为逆向流向波Fig.5 Variation of wave speed with ε at Pr=0.01,where a and d show upstream oblique waves,b shows downstream streamwise wave,c shows upstream streamwise wave

图6～图9 为Pr=0.3,1 和3 时Mac、波数、波传播角和波速随ε 的变化曲线.可以发现每个Prandtl数均存在3 种相同的临界模态,且Mac随Pr的增大而增大.对于Pr=1,当ε <0.017,ε >0.35 时(图6～图9 中的曲线e和h),临界模态为逆向斜波(θ>90°);0.017<ε<0.028 时(图6～图9 中的曲线f),临界模态为同向斜波(θ <90°);0.028 <ε <0.35 时(图6～图9 中的曲线g),临界模态为展向稳态模态(θ=90°,c=0).综合起来可以看出,对于Mac,曲线e随ε 的增大略有上升,而其余曲线均随ε 的增大而减小;对于波数,曲线h开始随ε 的增大有短暂的上升随后减小,而其余曲线均随ε 的增大而增大,其中曲线g的波数开始有明显变化,随后只有微小变化;对于传播角,曲线e和f均随ε 的增大而增大,曲线h则无明显变化;对于波速,曲线e的后半段存在轻微变化,曲线f随ε 的增大而减小,曲线h随ε 的增大有明显增大.其余Prandtl 数的曲线特征是类似的.

图6 Pr=0.3,1 和3 时Mac 随ε 的变化曲线,其中a,d,e,h,i,l曲线对应逆向斜波;b, f, j 曲线对应同向斜波;c,g,k 曲线对应展向稳态模态Fig.6 Variation of Mac with ε at Pr=0.3,1 and 3,where a,d,e,h,i,l curves correspond to upstream oblique wave;b, f, j curves correspond to downstream oblique wave;c,g,k curves correspond to spanwise stationary mode

图7 Pr=0.3,1 和3 时波数随ε 的变化曲线,其中a,d,e,h,i,l曲线对应逆向斜波;b, f, j 曲线对应同向斜波;c,g,k 曲线对应展向稳态模态Fig.7 Variation of wave number with ε at Pr=0.3,1 and 3,where a,d,e,h,i,l curves correspond to upstream oblique wave;b, f, j curves correspond to downstream oblique wave;c,g,k curves correspond to spanwise stationary mode

图8 Pr=0.3,1 和3 时波传播角随ε 的变化曲线,其中a,d,e,h,i,l曲线对应逆向斜波;b, f, j 曲线对应同向斜波Fig.8 Variation of wave propagation angle with ε at Pr=0.3,1 and 3,where a,d,e,h,i,l curves correspond to upstream oblique wave;b, f, j curves correspond to downstream oblique wave

图9 Pr=0.3,1 和3 时(I)波速随ε 的变化曲线,其中a,d,e,h,i,l曲线对应逆向斜波;b, f, j 曲线对应同向斜波Fig.9 Variation of wave speed with ε at Pr=0.3,1 and 3,where a,d,e,h,i,l curves correspond to upstream oblique wave;b, f, j curves correspond to downstream oblique wave

图10～图13 为Pr=100 时Mac、波数、波传播角和波速随ε 的变化曲线.可以发现曲线也存在3 种与Pr=1 类似的临界模态,当ε <0.12 时(图2～图5 中的曲线a),临界模态为逆向斜波(θ >90°);0.12<ε<0.55 时(图2～图5 中的曲线b),临界模态为同向斜波(θ <90°);ε >0.55 时(图2～图5 中的曲线c),临界模态为展向稳态模态(θ=90°,c=0).第一种模态的Mac随ε 的增大而增大,而后两种模态的Mac随ε 的增大而减小.综合各曲线可以看出,曲线a的波数和传播角均随ε 的增大而减小,而波速随ε 的增大而增大;曲线b的波数和波速均随ε 的增大而减小,而传播角无明显的变化;曲线c的波数随ε 的增大而增大.

图10 Pr=100 时Mac 随ε 的变化曲线,其中a 对应逆向斜波,b 对应同向斜波,c 对应展向稳态模态Fig.10 Variation of Mac with ε at Pr=100,where a shows upstream oblique waves,b shows downstream oblique wave,c shows upstream streamwise wave

图11 Pr=100 时波数随ε 的变化曲线,其中a 对应逆向斜波,b 对应同向斜波,c 对应展向稳态模态Fig.11 Variation of wave number with ε at Pr=100,where a shows upstream oblique waves,b shows downstream oblique wave,c shows upstream streamwise wave

图12 Pr=100 时波传播角随ε 的变化曲线,其中a 对应逆向斜波,b 对应同向斜波Fig.12 Variation of wave propagation angle with ε at Pr=100,where a shows upstream oblique waves,b shows downstream oblique wave

图13 Pr=100 时波速随ε 的变化曲线,其中a 对应逆向斜波,b 对应同向斜波Fig.13 Variation of wave speed with ε at Pr=100,where a shows upstream oblique waves,b shows downstream oblique wave

2.2 扰动流场

为了能更加清晰地了解流动失稳的机制,绘制了不同临界模态下的等温线图和流线图,并将最大扰动温度进行归一化.为了比较不同参数下的流场,将扰动的波长固定为2π,此时原点表示扰动初始位置,横轴数值表示扰动的相位,横轴方向为波的传播方向.

图14～图16 分别显示了Pr为0.01,1 和100 时,不同临界模态所对应的扰动流场.从图中可以看出,当Pr=0.01 时,各模态温度场的热点均分布在自由表面.在图14(b)流场中,一个大涡内还存在两个小涡,扰动速度在垂直方向上存在多次振荡,并且扰动基本分布在表面.当Pr=1 时,图15(c)与图16(c)几乎一致,而且除了逆向斜波的温度场热点分布在流场中间区域外,其余模态的热点均分布在自由表面.当Pr=100 时,斜波对应的温度场的热点均存在于中间区域,且各模态的流线分布几乎呈对称状态.

2.3 能量分析

扰动能量的变化率可以写成以下形式[15]

图14 Pr=0.01 时不同临界模态所对应的扰动流场Fig.14 Perturbation flow field of the different preferred modes at Pr=0.01

图15 Pr=1 时不同临界模态所对应的扰动流场Fig.15 Perturbation flow field of the different preferred modes at Pr=1

式中N为扰动应力做的功,M为Marangoni 力在表面做的功,I为扰动流与基本流之间的相互作用.这里将扰动动能进行归一化处理,即∫|u|2d3r=1.表3 中给出了不同Pr数下各扰动能量变化项的值.可以发现,随着弹性数增大,I由正变负.Pr=100 时,I几乎可以忽略,N>0 代表耗散,说明扰动能量的来源是Marangoni 力在表面做的功.Pr=1 时,M和N均大于0,而I既可能耗散能量又可能提供能量.Pr=0.01 时,各项值均可变化符号,并且N和I是主要的耗散或能量来源.由式(5) 可知在牛顿流体[32]中,τ=S,N始终大于0,对应黏性耗散.而对于黏弹性流体,弹性引起了τ与S之间的相位差,在一些情况下可能对扰动做正功.

图17 分别为Pr=0.01 时各临界模态扰动应力做功和扰动流与基本流之间的相互作用在垂直方向上的分布.各模态对应的弹性数ε 分别为0.001,0.01和0.1.其中

图16 Pr=100 时不同临界模态所对应的扰动流场Fig.16 Perturbation flow field of the different preferred modes at Pr=100

表3 不同Pr 数下各扰动能量变化项的值Table 3 Values of perturbation energy variation terms at different Pr

从图17 中可以看出,在液滴下半区域(z<0.5)各临界模态的PN和PI变化均较为平缓.而在液滴上半区域(z>0.5)中,同向流向波的PN和PI存在多次的振荡.对于逆向流向波,PN和PI的最值均在z=0.8 附近.

图17 Pr=0.01 时各临界模态扰动应力做功(a)和扰动流与基本流之间的相互作用(b)在垂直方向上的分布.其中,UOW 代表逆向斜波,DSW 代表同向流向波,USW 代表逆向流向波Fig.17 Distribution of(a)the work done by perturbation stress and(b)the interaction between the perturbation flow and the basic flow in vertical direction of the different preferred modes at Pr=0.01.Here,UOW stands for upstream oblique wave,DSW stands for downstream streamwise wave,USW stands for upstream streamwise wave

2.4 与热毛细液层对比

将液滴热毛细迁移与黏弹性热毛细液层[33]进行对比可以发现:在小Prandtl 数下,液层中存在展向稳态模态,而液滴中只有斜波和流向波;当弹性数足够大时,液层的临界模态均为展向稳态模态,而液滴在小和中Prandtl 数下的临界模态为斜波;对于展向稳态模态,液层的等温线几乎是垂直的,而在液滴中扰动温度在自由表面处存在振幅.