Bingham 流体双自由面热毛细液层的稳定性分析1)

2023-01-15王胜胡开鑫

力学学报 2022年12期

王胜胡开鑫

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

引言

当液体层表面存在水平温度梯度时,其诱导的表面张力梯度将驱动液体运动.这种流动被称为热毛细对流.由于其在晶体生长[1]中的重要作用,热毛细对流得到了广泛的研究.Smith等[2]采用线性稳定性方法研究了液层的热毛细不稳定.Patne等[3]对倾斜温度梯度液层进行稳定性分析.Chen等[4]研究了在零重力条件下由矩形池的等温侧壁引起的瞬态热毛细对流.

近年来热毛细效应的研究从单一自由面发展到双自由面.美国宇航局宇航员Pettit[5]在国际空间站上进行一系列热毛细流的微重力实验.其中环形水的液膜暴露在不均匀的温度分布中,这些实验表明,具有两个自由表面的液层有可能获得一种新的材料结晶过程[6].Ueno等[7]通过数值和实验方法研究在矩形孔液体薄膜中的热毛细流动.Messmer等[8]研究在热毛细作用力下的的双自由表面薄膜,表明在低Marangoni 数下有两种基本的流动结构.Toshiki等[9]研究具有两个气液界面的自由液体膜中的热毛细流动,发现该流动表现出从二维稳流态向三维振荡态的转变.Hu等[10]用线性稳定性分析方法研究双自由面热毛细液层的失稳机理.赵诚卓等[11]研究双自由面溶质-热毛细液层的不稳定性,发现在大Pr时,溶质毛细力的出现使流动更加稳定;在其他参数下,溶质毛细力会减弱稳定性.

非牛顿流体的热毛细对流广泛存在于镀膜[12]、薄膜干燥[13]、脱湿[14]、光刻[15]、喷墨打印[16]和微重力聚合物加工[17]等应用中.非牛顿流体与牛顿流体在流动特性上有很大的差异.黏弹性热毛细液层已被许多作者[18-20]研究过,结果表明:弹性导致的正应力会显著影响流动稳定性.Hu等[21]对剪切稀化流体的热毛细液层进行了线性稳定性分析.

黏塑性流体出现在许多工业应用和自然环境中,如聚合物加工[22]、熔浆[23]等.黏塑性流体的主要特征是其屈服应力,在应力较高时表现出类液体行为,在应力较低时表现出类固体行为.Bingham 流体是黏塑性流体的一种理想化模型.当剪切力未达到屈服值之前,Bingham 流体表现出弹性固体的性质而没有流动性,当超过屈服值之后,才发生流动且流体的黏度保持不变,类似牛顿流体[24].

近几十年来,对Bingham 流体进行了许多的研究.Nouar等[25]从特征值敏感性、最小缺陷和转捩标度律方面研究Bingham 流体在通道中流动的稳定性.Hu等[26]研究无限大液层中Bingham 流体热毛细对流的不稳定性.王世芳等[27]研究Bingham 流体在低渗透多孔介质中球向渗流的分形模型.韩建国等[28]使用同轴双圆柱流变仪测试Bingham 流体流变参数.Baioumy等[29]研究圆形管道入口区的Bingham 流体的流动.Fusi[30]提出一种基于有限差分的数值格式来研究Bingham 流体在通道中的非定常运动(平面泊肃叶流).Maurya等[31]研究Bingham 流体在垂直T 通道的子通道中的流动反转情况.Esmaeili等[32]分析Bingham 流体在非平行板间的挤压流动.Zhang等[33]利用射频加热法对稠油油藏开发中的Bingham 流体流动进行数值研究.

近年来,中国科学家开始对聚合物空间3D 打印技术进行试验.Dou等[34]提出在微重力条件下设计和制造高精度陶瓷组件的方法,其中制备的陶瓷浆料的流变特性类似Bingham 流体.本文以陶瓷浆料的空间3D 打印为应用背景,研究Bingham 流体在双自由面液体层中的热毛细对流及其稳定性.本文的组织结构如下.在第1 节中,给出了该问题的物理模型和数值描述.推导出了基本流的解析解和控制方程.然后在第2 节中,得到了不同Bingham 数(B)和垂直温差(Q)时的临界参数;显示了临界模态的流场图,研究了能量机制.最后在第3 节进行总结.

1 控制方程

本文研究的流动模型如图1 所示.从狭缝中射出的聚合物流体液层,其表面的水平温度梯度会导致内部的热毛细对流.d是液体层深度的一半,τ13是表面的剪应力,U0是速度场.在该参考系中,流体流向的平均流速为0.具有两个自由表面的液层在其上下表面上受到恒定的温度梯度b=-dT/dx>0的影响.其中x,y,z为流向、展向和法向.其中表面张力σˆ随温度T的变化满足 σ ˆ=σˆ0-γ(T-T0) 的关系.由于层内部的剪切速率小于表面附近的值,因此液层中间会有一个栓塞区域.该流动由3 个区域组成,III 是栓塞区,I和II 是屈服区域,其范围分别是(-z0,z0),[-1,-z0],[z0,1].

图1 Bingham 流体双自由面热毛细液层示意图Fig.1 The thermocapillary liquid layer with two free surfaces for a Bingham fluid

流动的无量纲控制方程组由下式给出,分别为连续性方程,动量方程和能量方程

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

采用Bingham 流体模型来反映聚合物的本构关系.其本构方程为

z=±1的界面为液层自由面,边界条件为

其中速度场只有x方向的分量,Tb为垂直方向上的温度分布.在液层整体参考系中,任意垂直截面上的质量流量为0,可推出基本流,详见附录A.图2是基本流的速度分布和垂直温度分布,从图2(a)可知,B越大,栓塞区域越大;从图2(b)可知,随着B的增加,垂直温度梯度越来越小.

图2 基本流的(a)速度分布和(b)垂直温度分布(Ma=100,Q=0)Fig.2 (a) The velocity distribution and(b) vertical temperature distribution at Ma=100for the basic flow

式(13)中Q是上、下自由面垂直温度差

设当Q=0时,Tb(±1)=0,根据边界条件可得到Q1,Q2和Q的关系

2 数值结果

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

本章讨论了Q=0和Q＞0两种情况,得到了不同参数下的临界曲线,画出了对应的扰动流场图,并分析了其中的能量机制.

2.1 临界曲线

2.1.1Q=0

本节画出了Q=0,Bi=0时不同B下的临界曲线.其中Mac随Pr的变化如图3(a)所示.可以发现,Mac随着B的增加而增大;当B=0时,流体为牛顿流体,图中显示有斜波和流向波两种临界模态,在较小Pr,最不稳定的模态是斜波,在大Pr下,流向波比斜波更不稳定,因此临界模态发生变化.临界曲线的交叉处是两种波中性MaN(σr=0)相等的情况.当B=0.01 时,作为Bingham 流体的临界参数发生明显突变,且在大Pr下,模态由同向流向波(θ=0°)变成了逆向斜波 θ ∈(90°,180°) ;在中小Pr下的临界模态是逆向斜波,而在大Pr下则为逆向流向波(θ=180°);当B=0.6 时,临界模态存在逆向斜波、同向斜波θ ∈(0°,90°)以及同向流向波3 种.结合图3(b)综合来看,k随着B的增加而增大.在中小Pr下,牛顿流体的波数变化很小,而Bingham流体的波数随着Pr的增加而明显增大.图3(c)为波传播角 θ随Pr变化的曲线,可以看出牛顿流体在Pr=1.2 的时候,θ 降为0°,之后随着B的增加而变大.在大Pr的情况下,B=0.2 时,θ 变成180°,而对于B=0.6,θ 则变为0°.

图3 Q=0,Bi=0时,不同B 下,(a) Mac随Pr 的变化曲线以及(a)所对应的(b)波数和(c)波传播角.曲线对应:a,c,d,e,g,h为斜波,b,f,i,j为流向波Fig.3 (a) The variation of Macwith Pr under different B at when Q=0and Bi=0,(b) wave number and(c) wave propagation angle.The curves correspond to:a,c,d,e,g,h oblique wave andb,f,i,j streamwise wave

此外可以发现,Bi的增大并不会影响临界模态的种类,而只会增强流动的稳定性.

2.1.2Q＞0

本小节考虑上、下自由面有温度差时的临界模态.图4为B=0.2,Ma=100时垂直方向上的温度分布,可以看出,当Q=0时,Tb是关于z=0对称,Q的增加导致了Tb的不对称.

图4 B=0.2,Ma=100时垂直方向上的温度分布Fig.4 Vertical temperature distribution at B=0.2 and Ma=100

图5为Q=0.05,Bi=0时不同B下的临界曲线.可以看出,牛顿流体存在逆向斜波、同向流向波、展向稳态模态3 种临界模态,而Bingham 流体虽然在B=0.01 时存在与牛顿流体相同的临界模态,但在B=0.2时,临界模态只有逆向斜波和展向稳态模态两种.对于斜波和流向波模态,Mac基本随着Pr的增加而增加,而展向稳态模态则与之相反.此外对于同种模态,B的增加对 θ 的影响较小.

图5 Q=0.05 时,(a) Mac随Pr 的变化曲线以及(a)所对应的(b)波数和(c)波传播角.曲线对应:a,d,g为斜波,b,e为流向波,c,h,f为展向稳态模态Fig.5 (a) The variation of Macwith Pr under Q=0.05,(b) wave number and(c) wave propagation angle.The curves correspond to:a,d,g oblique wave,b,e streamwise wave,c,h,f spanwise stationary mode

2.2 扰动流场

本节画出了不同临界模态下的等温线和流线图,并将最大扰动温度进行归一化.为了比较不同参数下的流场,将扰动的波长固定为 2π,此时横轴数值表示扰动的相位,横轴方向为波的传播方向.

图6 显示了Q=0,Bi=0时临界模态所对应的扰动流场.从图6(b) 图可以发现,在B=0.2,Pr=0.01 时,扰动温度分布在屈服区和栓塞区,且呈对称分布;在图6(d)中可知,B=0.4,Pr=20时的扰动温度呈反对称分布,且扰动温度只存在于屈服区.

图6 Q=0,Bi=0时不同临界模态所对应的扰动流场:逆向斜波(B=0.2,Pr=0.01)的(a)上屈服面扰动,(b)温度和流线;逆向斜波(B=0.4,Pr=20)的(c)上屈服面扰动,(d)温度和流线Fig.6 The perturbation flow field of the different preferred modes at Q=0and Bi=0:(a) upper yield surface perturbation,(b) temperature and streamlines of the upstream oblique wave(B=0.2,Pr=0.01);(c) upper yield surface perturbation,(d) temperature and streamlines of the upstream oblique wave(B=0.4,Pr=20)

图7 显示了Q=0,Bi=5 时的扰动流场.从图中可以看出,在小Pr下,扰动温度呈对称分布,随着B的增加,扰动温度在栓塞区呈截断的状态.

图7 Q=0,Bi=5 时不同临界模态所对应的扰动流场:逆向斜波(B=0.4,Pr=0.01)的(a)上屈服面扰动,(b)温度和流线;逆向斜波(B=0.6,Pr=0.01)的(c)上屈服面扰动,(d) 温度和流线Fig.7 The perturbation flow field of the different preferred modes at Q=0and Bi=5:(a) upper yield surface perturbation,(b) temperature and streamlines of the upstream oblique wave(B=0.4,Pr=0.01);(c) upper yield surface perturbation,(d) temperature and streamlines of the upstream oblique wave(B=0.6,Pr=0.01)

图8 显示了Q=0.05,Bi=0时的扰动流场.从图8(a)中可以看出牛顿流体的扰动温度在整个流场中分布不规则,且主要集中在下半部分流场区域;而从图8(c)中可知,Bingham 流体的扰动温度分布在下半部分流场区域内.此外由于栓塞区的存在,Bingham 流体中无法找到如图8(a)所示的贯穿液层的涡胞结构.

图8 Q=0.05,Bi=0时不同临界模态所对应的扰动流场:同向流向波(B=0,Pr=10)的(a)温度和流线;同向斜波(B=0.2,Pr=0.1)的(b)上屈服面扰动,(c)温度和流线Fig.8 The perturbation flow field of the different preferred modes at Q=0.05 and Bi=0:(a) temperature and streamlines of the downstream streamwise wave(B=0,Pr=10);(b) upper yield surface perturbation;(c) temperature and streamlines of the downstream oblique wave(B=0.2,Pr=0.1)

2.3 能量分析

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

式中N为扰动应力做的功,M为Marangoni 力在表面做的功,I为扰动流与基本流之间的相互作用.这里将扰动动能进行归一化处理,即表1中给出了不同参数下各扰动能量变化项的值.可以发现在Q=0,Bi=0时,Bingham 流体的I在小Pr数下比牛顿流体的值大很多,因此不能忽略.然而Pr=100时,I则可以忽略不计,N＞0代表耗散,说明扰动能量的主要来源是Marangoni 力在表面做的功.