陈上通 吴 笛 王 佳 段 俐 康 琦,

* (中国科学院力学研究所微重力实验室,北京 100190)

† (中国科学院大学工程科学学院,北京 100049)

引言

表面张力驱动流,在此定义为由表面张力引起的附加压力驱动的自发界面流,是航天器贮箱内液体行为的重要组成部分.为了在微重力环境下下高效地管理液体,有必要深入研究表面张力问题.

Washburn[1]首次提出了毛细管内液体爬升的运动方程.他们证实该过程可以用毛细驱动压力、黏滞阻力和重力之间的平衡来描述,推导出了著名的Lucus-Washburn 方程.Concus 和Finn[2]探索了具有内角的容器中的液体平衡界面,并提出了著名的Concus-Finn 条件.Levine 等[3]考虑了表面张力驱动压力、黏滞阻力和对流损失,推导出了圆管内液体爬升高度的二阶微分方程.Stange 等[4]提出了一个更全面的圆管内液体爬升模型,它补充考虑了弯液面的重定位过程、动态接触角和泊肃叶流动的启动过程.Jiang 等[5]在1979 年改进了动态接触的计算方法并被广泛接受.Dreyer 等[6]探究了平行平板之间的表面张力驱动流动并得到了液体爬升高度方程.Weislogel 和Lichter[7]获得了液体沿尖内角处的流动方程,并扩展到由不同润湿性形成的尖内角处[8]和圆内角处[9].Yue 和Wang[10]采用数值方法研究了三维自由表面的动力学问题.Higuera 等[11]研究了润湿液体在两个垂直板之间形成一个小角度的狭窄间隙中的流动过程.Wolf 等[12]提出了一种基Lattice-Boltzmann 方法、通过考虑流体和固壁之间长程相互作用的影响来模拟平行平板之间的液体爬升.Bolleddula 等[13]提出了一种新的沿不同截面形状尖内角处流动的解析解,并探讨了实验模型几何形状的影响.魏月兴等[14]研究了内角处表面张力驱动流动理论并将其应用到板式贮箱的结构设计上.李永等[15-17]利用落塔试验和数值模拟探究了板式贮箱的流体管理性能,得到了丰富的微重力试验数据.Reyssat[18]探究了平板和圆弧板之间形成的液桥界面形貌特征,建立的模型可以准确预测液桥外形.Wu 等[19]研究了弯曲尖内角处的流动过程,通过修正尖内角曲率半径带来的影响,得到了流动距离与曲率半径的关系.Romain 等[20]研究了轴对称几何模型中的表面张力驱动流动.Chen 等[21]利用毛细驱动流动理论对贮箱内PMD 结构进行了优化,并开展了落塔实验和基于VOF 法的仿真分析,检验了其性能.Cheng 等[22]研究了变直径圆管中的流动问题.Chen等[23]研究了微重力下椭圆管内的表面张力驱动流动,并提出了两段式的流动新模型.

本文探究了微重力环境下成一定角度平板间的表面张力驱动流,考虑了动态接触角、对流压力损失、黏滞阻力和液池内弯曲自由面的影响,并利用基于VOF 方法的数值模拟进行验证.此外,微分方程可变换成由作用在板间控制体上一系列合力的方程,通过同时考虑两个主要作用力,将整个爬升过程分成三个阶段.

1 理论推导

研究模型如图1 所示,包含圆柱形液池和自上方插入液池一定深度的平板模型,液池内壁和平板外壁均有防爬挡板,选择直角坐标系来分析该问题.初始液面为水平面,z 向原点位于水平面上.液体爬升高度为h,液体爬升平均速度为.平板水平截面如图2 所示.平板左端的距离是2a,右端的距离是2c,板的宽度是b.平板上特定位置与和x 轴间的距离是d.液池内弯曲界面的曲率半径是用防爬挡板之间的距离i 与液池内气液界面中心线和平板模型对称线之间的距离j 计算.为了计算i 与j 的值,平板之间的梯形入口等效成半径为re的圆形入口,半径re为理论分析的基本假设是

图1 研究模型正向剖视图Fig.1 Front view of the model

图2 平板模型的水平截面Fig.2 Horizontal cross section of the plates

(1)气液界面上的剪切力忽略不计;

(2)流动过程是恒温的;

(3)流动是充分发展的泊肃叶流;

(4)流体是牛顿流体,不可压缩且均质;

(5)液体与壁面之间无滑移.

平板上特定位置与和x 轴间的距离为

流动是充分发展的泊肃叶流,液体仅具有z 向速度分量u,流动边界条件为

u 为关于x,y,t 的函数,t 为时间.可以认为板间区域由x 向上的无数个dx 分段组成,在每个dx 分段内u 为关于y,t 的函数,与x 无关.结合式(2)并根据N-S 方程可得每个dx 分段内速度场分布为

结合平板间流量公式

式中 Ω 代表板间水平截面全域.可将式(3)转化为

N-S 方程z 向分量去除零项后为

式中 ρ 为密度,ν 为液体运动学黏性系数.对该式在板间水平截面上对x 和y 进行积分,得到

将式(7)对z 从z=0 到z=h 进行积分,可以得到

采用与文献[14]相似的处理方法,可以得到板间顶端凹液面的压力为

式中 p0为液面附近大气压,σ 为表面张力,αd和 αs分别为动态接触角和静态接触角.对于完全润湿的液体,即静态接触角为0°,其动态接触角有一经验公式

式中μ 为动力学黏性系数.该计算模型受到了广泛运用和验证,文献[4,6,20]均采用该计算模型.

为了计算板间入口处的作用力,在等效圆形入口附近建立半径为re的半球形控制体2,如图3所示.

图3 板间入口处的等效半球形控制体Fig.3 Equivalent circular entrance and the control volume around the inlet

式中 pR为液池内弯曲界面引起的附加压力,为

液池中自由面的曲率半径Rc的近似表达式为

在板间入口处的作用力沿z 向分力为

经半球面R=re进入控制体2 的流量z 向分量为

经板间入口流出控制体2 的流量z 向分量为

经半球面R=re进入控制体2 的加速度流z 向分量为

经入口处流出控制体2 的加速度流z 向分量为

因此,半球形控制体2 内的整体加速度的近似表达式为

他的脚离我的头只有几英寸远，我理应去安慰他，我本应该主动去安慰他才对，因为我从小就是受这种教育长大的。相反，我觉得那样做很恶心。看起来那么强壮的人，不应该表现得这么脆弱。为什么他不能像其余人一样悄悄地哭呢？

结合式(11)～ 式(18)及控制体2 上的动量方程,可以求得板间入口处的作用力z 向分量为

将式(1)、式(5)、式(9)、式(10)和式(20)代入式(8),可以得到液体爬升高度h 的二阶微分方程

该二阶微分方程可利用四阶Runge-Kutta 法并结合初始条件 h(t=0)=h˙(t=0)=0 求解.

2 数值模拟

本文建立的三维计算模型如图4(a)所示.由于液池内气液界面弯曲造成的影响很小,为了建模的方便,选择边长为140 mm 的方形液池来代替理论分析中的圆柱形液池,并在理论计算中使用79 mm的等效半径,该简化所造成的影响可以忽略不计,模型的总高度为120 mm.建立的网格模型如图4(b)所示,网格总数约为100 万,在所有壁面附近建立边界层网格,首层边界层的高度约为0.1 mm,相邻两层网格膨胀比为1.2.由于计算模型简单,网格质量好,调整边界层网格后进行第二次仿真,两次结果差异很小,验证了网格对仿真结果的无关性.同时为了减小仿真可能产生的随机误差,取两次结果的平均值作为最终结果.

图4 仿真模型Fig.4 Numerical model

本文使用KF-96 系列硅油开展研究工作,该系列型号硅油以其运动学黏性系数作为标号,比如2 号硅油(SF 2)、5 号硅油(SF 5)和10 号硅油(SF 10),他们的物性参数如表1 所示.本文选择层流作为流动模式,压力-速度耦合方程由SIMPLIC 算法进行数值求解,压力空间离散式方程选用“Body Force Weighted”算法求解,梯度空间离散化方程选择“Least Square Cell”算法求解,动量空间离散化方程使用二阶迎风格式求解,松弛因子采用默认设置,当方程迭代余量减少到10-6时,计算被认为收敛.在计算过程中,Courant 数(Co)大部分时候小于1,这表明计算过程非常稳定.Co 对于瞬态流动很重要.对于一维网格,它的定义为

表1 硅油物性参数(25 °C)Table 1 Liquid properties (25 °C)

式中,u 代表流速,Δx 代表网格尺寸,Δ t 代表时间步长.

板间流动的仿真侧视图如图5 所示,黄色曲面代表气液界面.在t=0 s 时,液体全部位于液池底部,气液界面为一平面.仿真开始后,液体快速向上爬升并在板间形成一弯曲液面.液体与平板的接触线不是直线,如图6 所示.因此,为了计算液体爬升高度,选择接触线的两端点H1,H3和最高点H2,将其平均值(H1+H2+H3)/3 视为液体爬升高度h.

图5 仿真结果侧视图.液体为10 号硅油,a=2 mm,b=14 mm,c=5 mmFig.5 Front view of the numerical model as the liquid is SF 10 and a=2 mm,b=14 mm,c=5 mm

图6 板上液相分布,红色部分代表液体Fig.6 Liquid distribution in a plate and the red part represents liquid

理论和仿真爬升高度的对比如图7 所示.本文使用了6 种尺寸的平板模型.图中实线代表理论结果,方块代表仿真结果,并用不同颜色代表不同型号硅油的爬升数据.黑色代表2 号硅油的数据,红色代表5 号硅油的数据,蓝色代表10 号硅油的数据.由于模型高度的限制,对于2 号硅油,取其2 s 内的爬升数据;对于5 号硅油,取其4 s 内的爬升数据;对于10 号硅油,取其6 s 内的爬升数据.本文忽略了板间液面从一开始的平面到弯曲平衡界面的发展过程.此外,对于爬升高度的测量方法也将带来一定误差.总体来看,数值结果与理论结果吻合良好.对比图8(a)和图8(b),或者图8(c)和图8(d)可以看出,相同黏性的硅油,在较宽间距板间的爬升速度比较小间距板间的慢,平板之间的距离越宽,板间流速越慢.从理论分析中可知,表面张力的增加有助于增加毛细驱动压力,加快液体流速,而液体黏性的增加会减小流速.

图7 液体爬升高度随时间变化Fig.7 Liquid climbing height h vs.time t

3 机理分析

为了进一步分析,可以将式(14)转化成作用在控制体1 上合力的方程,如下

式中各力的含义为

(1)板间毛细驱动压力

(2)液池内惯性力

(3)板间入口对流压力损失

(4)板间黏滞阻力

(5)板间惯性力

(6)液池内黏滞阻力

(7)液池内附加压力

不同工况下各力发展情况如图8 所示,横坐标是时间t,纵坐标是力F.深红色、绿色、淡蓝色、深蓝色和洋红色实线分别代表平板间毛细驱动压力Fcp、液池中的惯性力Fir、板间入口处的对流压力损失Fpl、两块平板上的黏滞阻力Ffp和液池内的附加压力Fcr,黑色虚线代表平板间的惯性力Fip.毛细驱动压力从最大值开始变化,迅速减小至最小值后又开始缓慢增加,该变化是由于动态接触角的变化.从上述四图中可以看出,无论是高Oh数(Oh=还是低Oh数工况,液池内惯性力、板间入口压力损失和板上黏滞阻力均先后发挥了主要作用,这点与椭圆管和同心圆管内的流动特征大不相同.液池内惯性力从最大值开始变化,迅速减小,直至被对流压力损失超越;入口处对流压力损失从0 开始变化,迅速增加到最大值后开始缓慢减小.最终在某一时刻,板上黏滞阻力超越入口压力损失发挥主要作用.因此,板间的表面张力驱动流动均可以划分成三个阶段.在第一阶段中,令毛细驱动压力等于液池内惯性力,可以得到

图8 各力随时间变化情况Fig.8 Forces’ development vs time

该方程的初始条件为 h˙(t=0)=0 .结合初始条件可求得第一阶段的爬升高度为

为了计算第二阶段的速度方程,令毛细驱动压力等于入口处压力损失,可以得到

令第一、第二阶段的速度方程相等,可以求出划分这两个阶段的时刻t1

结合第二阶段的初始条件h1(t1),可得到第二阶段爬升高度为

为了计算第三阶段的速度方程,令毛细驱动压力等于板上黏滞阻力,可以得到

令第二、第三阶段的速度方程相等,可以求出划分这两个阶段的时刻t2为

结合第三阶段的初始条件h2(t2),可得第三阶段的爬升高度为

不同工况的计算参数如表2 所示,利用Oh 数区分不同工况.

相比于高Oh 数工况,低Oh 数工况下第二阶段持续时间更长.对比图8(a)和图8(b),或者图8(c)和图8(d),均可看出在低Oh 数时入口对流压力损失发挥主要作用的时间更长,这意味着第二阶段持续时间更长.例如,在第四个仿真工况中,其流动Oh 数为0.004 50,第二阶段持续0.647 s;而在第六个仿真工况中,其流动Oh 数为0.022 1,第二阶段仅持续0.039 s,远小于低Oh 数工况下的持续时间.楔形平板间的表面张力驱动流可以分成三个阶段,如图9(a)所示,x 轴是时间t,y 轴是液体爬升高度h,蓝线代表爬升高度与时间的关系,图中标注了三个阶段的范围.第一阶段中爬升曲线为抛物线型,液体爬升高度与t2成正比;第二阶段中爬升曲线为直线型,液体爬升高度与t 成正比;第三阶段液体爬升高度与t1/2成正比.理论和仿真爬升高度的对比如图9(b)所示,其中黑色实现代表理论结果,方块代表仿真结果,划分不同流动阶段的时刻t1和t2均标注在图中,为了更直观地观察爬升高度的变化规律,建立了一条从时刻t1开始的虚直线.仿真结果和理论结果吻合良好,爬升高度的变化规律也与理论分析一致,充分验证了理论分析的正确性.

4 结论

本文通过理论分析,推导出了微重力环境下成一定夹角平板间表面张力驱动流爬升高度的二阶微分方程,并利用基于VOF 方法的数值模拟开展了验证工作.本文采用了6 种研究模型和3 种型号硅油,理论结果与仿真结果吻合良好.

本文通过进一步的理论分析,还发现同时考虑两个主要作用力,可将平板间的表面张力驱动流动分为三个阶段,在这三个阶段中液体爬升高度h 分别与t2,t 和t1/2成正比.在第一阶段,液池中的惯性力和毛细驱动压力发挥主要作用.但是液池中的惯性力迅速减小,流动进入第二阶段,由对流压力损失和毛细驱动压力发挥主要作用.并且,在高Oh 数情况下,第二阶段的持续时间很短暂.第三阶段时由板上黏滞阻力和毛细驱动压力发挥主要作用.相较于第三阶段,第一、第二阶段持续时间都很短暂,整个流动过程以第三阶段为主.本章研究结果可为板式贮箱中导流板和蓄液叶片等的设计和空间流体管理提供理论依据.