Oldroyd-B 黏弹性液滴碰撞过程的数值模拟1)

2022-04-07关新燕富庆飞杨立军

力学学报 2022年3期

关新燕富庆飞,2) 刘虎杨立军

* (北京航空航天大学宇航学院,北京 100191)

† (北京宇航系统控制研究所,北京 100076)

引言

凝胶推进剂因具有安全[1-4]、高比冲[5-6]、推力易调节[7]、钝感与环保[8-9]等优点,在导弹武器系统和快速发射航天运载器动力系统等方面有着广阔的应用前景[10-12].但是复杂的流变特性使其雾化过程存在一定困难.凝胶推进剂一般通过添加聚合物胶凝剂配置而成[13-16],制得的推进剂具有黏弹性效应,在雾化过程中会产生黏弹性液丝和液滴等凝胶碎片.

对于牛顿液滴,在各种撞击条件下进行碰撞实验,得到了合并、反弹、反向分离和拉伸分离的碰撞结果[17].对于非牛顿液滴的碰撞也已经进行了一些实验,Finotello 等[18]通过液滴碰撞实验研究了剪切变稀非牛顿流体黄原胶的碰撞行为,尽管其中存在着复杂的碰撞动力学,仍然基于无量纲韦伯数(We)和冲击因子(B)得到了完整的相图,观察到了圆盘状以及碰撞液滴的振荡行为,这与黏性能量耗散的增加和拉伸作用有关.对于液滴的二次雾化已经有相关的数值研究[19-23].Khare和Yang[24]使用VOF 捕捉了非牛顿液滴的变形和破裂,分析了幂律型液滴破裂的基本物理原理.发现当液滴向下减速时,会在流动方向上进行拉伸,不同的剪切速率导致液滴不同点的黏度不相等,使拉伸变得不对称,最终液滴形成凹坑直到破裂,破裂时出现串珠结构,并且不平衡力产生的扭矩使液滴旋转.G u p t a和Sbragaglia[25]数值研究了牛顿/黏弹性液滴在黏弹性/牛顿流体基质中受到剪切时的变形和破裂.采用Lattice-Boltzmann 模型(LBM)模拟两种具有可变黏度比(液滴与基质的黏度比)的不混溶流体,有限差分法(finite difference,FD)用于黏弹性的建模,并使用FENE-P 本构方程解释了聚合物的动力学.结果发现在无穷大基质中黏弹性的影响很小,而在基质局限的情况下会有较明显的影响.在De(应力松弛时间与观测时间的比值)较小的情况下,与液滴是黏弹性的情况相比,基质的黏弹性改变受限液滴稳定性的现象更为明显.Izbassarov和Muradoglu[26]计算研究了黏弹性对固体表面上的液滴撞击,扩散和反弹的影响,采用FENE-CR 模型.发现更大的Wi(对于特征长度与速度尺度单一的物理问题,Wi与De相一致)的黏弹性更有利于液滴回弹,并且在低韦伯数和高雷诺数下,黏弹性的影响更为明显.此外,增加平衡接触角和雷诺数、减小韦伯数都有助于液滴的回弹.通过改变聚合物黏度研究聚合物浓度的影响,结果是液滴的回弹随着聚合物黏度的增加而被抑制,这主要归因于黏性耗散的增强.在边界层中,界面附近的聚合物会产生拉伸应力,从而推动接触线以增加铺展速率,这是聚合物液滴比牛顿液滴铺展更广的主要原因.

文献[27-28]对幂律液滴相撞进行了数值模拟.Liu 等[27]通过VOF 方法模拟了幂律型凝胶推进剂液滴的正向碰撞,研究了反弹、合并和反向分离液滴的现象,结果表明液滴的反弹取决于韦伯数和流体的黏度.在低韦伯数下,液滴的最小中心厚度出现晚于其最大变形,这与在高韦伯数下的现象相反.而且,复杂的流场表明最大剪切速率出现在流动被重新定向和加速的点.Sun 等[29]采用LBM 方法在二维坐标系中模拟了非牛顿液滴的碰撞,Carreau-Yasuda (CY)模型用于描述非牛顿流变学对液滴的合并、分离和内部混合的影响,We值的范围大约在20 到200 之间.结果表明,剪切稀化黏度对碰撞动力学的影响非常复杂,并且强烈依赖于流体流变学,随着非牛顿效应的增加,剪切稀化流体促进了碰撞液滴的内部混合,变形和分离.对于剪切稠化流体,由于随着剪切增稠而增加的黏性耗散,液滴变形被显着抑制,促进了永久合并.Focke和Bothe[30]采用直接数值模拟(direct numerical simulation,DNS)和VOF 方法,对Motzigemba 等[31]的实验数据进行了研究,结果表明对于牛顿流体,可以通过选择有效黏度产生与剪切稀化流体相同的液滴碰撞动力学.他们还提出了一种液膜稳定方法[32],以证明高We值下的偏心碰撞也存在有效黏度,非牛顿液滴碰撞会在碰撞复合体的中心形成薄的液膜,只有在碰撞薄片保持完好无损时,得到的结果才可靠.Motzigemba等[31]基于流体体积法,对非牛顿液滴的碰撞进行了数值模拟,发现剪切稀化流体的碰撞液滴变形比牛顿液体的更大.液滴碰撞达到最大变形所需的无量纲时间与黏度无关,但与韦伯数有关,从得到的速度场可以看出,在所有碰撞过程中拉伸流占主导.

现有针对黏弹性液滴碰撞的研究大多将其看作纯剪切变稀流体,并以幂律模型等对其剪切变稀的流变性质进行拟合,较少考虑流体具有的弹性的影响.因此本文采用Oldroyd-B 本构描述液滴的黏弹性,对两个等体积液滴的碰撞过程进行了直接数值模拟,研究了松弛时间、黏度比、韦伯数和碰撞因子等参数对黏弹性液滴碰撞过程的影响规律,以期为凝胶推进剂雾化机理的理解和雾化性能的提升提供参考数据.

1 数值计算方法与模型

1.1 数值算法

数值模拟采用Basilik 开源代码[33],将网格自适应与流体体积法、连续表面力模型(continuum surface force,CSF)相结合,得以有效追踪界面的运动变形和精确计算表面张力.

Basilisk 使用投影法求解流体控制方程(Navier-Stokes 方程)

式中,u=(u,v,w) 为流体的速度,ρ≡ρ(x,t) 为流体的密度,μ≡μ(x,t) 为动力黏度,D为应变张量,定义为Dij≡(∂iuj+∂jui),狄拉克分布函数 δs表示表面张力项集中在界面上,σ 为表面张力系数,κ和n是曲率,并垂直于界面.

对于两相流,引入第一相的质量分数c(x,t),定义密度和黏度为

式中,ρ1,ρ2和μ1,μ2分别是第一相和第二相的密度和黏度.c˜ 等于c,或者通过对c应用一个平滑的空间滤波器构造.

密度的平流方程可以用体积分数的等效平流方程来代替

Basilisk 在空间上使用基于四叉树(二维,三维时为八叉树)的正交笛卡尔网格对计算域进行离散.

对于两相黏弹性流动的数值求解,考虑到聚合物的记忆效应,在N-S 方程中引入聚合物应力张量τp

式中,μs是溶液黏度,聚合物张量 τp通常是构象张量A的函数fs(·)

式中,μp为聚合物黏度,λ 为松弛时间.构象张量A可以看作是测量聚合物链分子变形的内部状态变量,假设它对称且正定,则有

上对流导数算符 ∇ 由下式给出

应变函数fS(A)和松弛函数fR(A)可以通过本构方程推导得到,对于Oldroyd-B 本构方程有

求解时使用构象张量[34-35]的对数 Ψ=lgA代替黏弹性应力张量,来解决黏弹性流体数值求解时常遇到的高维森伯格数问题[36],即强弹性流动产生高应力和精细特征的区域时出现数值不稳定性.

1.2 计算模型

图1 为黏弹性液滴碰撞的仿真模型,初始时刻两个液滴的直径为D0,圆心在x方向相距1.05D0,在y方向相距X,以大小相等、方向相反的速度u进行碰撞,碰撞因子B=X/D0.计算域为6D0×6D0的正四边形,初始速度的方向为x方向,与初始速度垂直的方向为y方向.采用Oldroyd-B 模型表示液滴的黏弹性

图1 黏弹性液滴碰撞的计算模型Fig.1 Computational model of viscoelastic droplet collision

黏度比定义式为β=μs/μ0,其中μ0为零剪切黏度,它是溶液黏度与聚合物黏度之和(μ0=μs+μp).当黏弹性液滴进行正撞,即B=0 时.如图2 所示,t=0.01 ms 为液滴初始运动时的网格,t=0.05 ms 为液滴即将碰撞时刻的自适应网格.对于两液滴间隙的捕捉,由于本文没有使用额外的计算模型,因此会使仿真现象产生一定的误差,未来如果进一步加入液滴气缝的解析模型[37],可能得到更加准确的结果.

图2 不同时刻下液滴撞击的自适应网格Fig.2 Adaptive mesh of droplet collision at different time

图3 为黏弹性液滴的正撞过程,此时液滴的初始速度为1 m/s,黏度比β=0.1,松弛时间λ=1 ms,图4 为发生正撞后液滴尺寸随时间的变化,这里采用单个液滴的最大宽度Dmax来表示.可以看出液滴从初始时刻的直径为1 mm 开始,在碰撞挤压的过程中扩散,Dmax会变大,在将近1 ms 时Dmax达到最大值,此时初始撞击动能耗尽,除了克服间隙气体流动耗散掉的能量,动能部分转化为表面能,部分存储为弹性能.最后速度方向改变,液滴内的流体开始向反方向运动,直径又变小,此处称这一过程为回缩,注意此处的回缩液滴并没有分离,后面呈现出小幅度振荡的过程.将液滴的初始速度减小到0.5 m/s,液滴在碰撞后并没有融合,而是出现反弹的现象.图5 为黏弹性液滴碰撞反弹的仿真结果.由于在两液滴靠近后,其动能较小导致不足以将缝隙中的高压气体排除,也没有等到足够的时间让气体自行排出缝隙,就在互相接触之前失去了撞击的惯性,随后液滴在表面张力的作用下试图恢复圆形表面,此时变形了的部分使液滴最终向反方向弹开.

图3 黏弹性液滴的碰撞融合过程(u=1 m/s)Fig.3 Collision and coalescence process of viscoelastic droplet(u=1 m/s)

图4 液滴尺寸随时间的变化曲线Fig.4 Curve of droplet size over time

图5 黏弹性液滴的碰撞反弹过程(u=0.5 m/s)Fig.5 Collision and bounce process of viscoelastic droplet(u=0.5 m/s)

2 结果与分析

2.1 松弛时间对黏弹性液滴正撞过程的影响

图6 为不同松弛时间下液滴最大宽度随时间的变化.可以看出松弛时间越大,液滴的变形程度越大.在液滴靠近与碰撞的过程中,λ越大时,Dmax的峰值越大,即液滴被挤压地越扁.这是由于松弛时间增大时,黏弹性应力增大,而黏弹性有利于液滴的扩散.在液滴回缩的过程中,λ越大时,Dmax的最小值也越小,即在较大的松弛时间下液滴的回缩程度较大,在λ=4 ms 时,Dmax可以减小到初始直径D0以下.可见不同松弛时间下的液滴回缩阶段差异也很大,此时在挤压阶段存储的弹性能释放出来,有助于回缩.在液滴振荡的过程中,λ越大时,振荡的幅值更大、相应的振荡周期越长.

图6 不同松弛时间λ 下液滴最大宽度Dmax 随时间的变化Fig.6 The curve of droplet maximum width Dmax with time at different relaxation time λ

采用如下公式计算液滴的动能KE

式中,ρ为密度,u和v分别为x方向和y方向的速度.在计算模型中通过计算每个网格的动能,然后将所有网格的动能相加来实现系统动能的计算.

液滴表面能SE的计算公式为

这是表面能的二维类比形式.式中σ为表面张力系数,l为气液界面的总长度.

根据式(14)和式(15)计算得到不同松弛时间下液滴碰撞过程中动能和表面能随时间的变化.图7为不同松弛时间下,系统动能随时间的变化曲线.从图中可以看出在系统的初始动能都相同的情况下,松弛时间越大,动能耗散的速率越慢.之后在不同的松弛时间下,动能的变化曲线有相位差,松弛时间增大时,动能变化曲线的拐点向后延迟,振荡幅值与周期也越大.动能在一开始系统迅速下降到最小值后,又开始上升,松弛时间越大时动能上升的越多,从局部放大图可以看到在随后的几次振荡中动能的变化也遵循这一规律.动能峰值随着振荡的次数递减,由于松弛时间越大时动能耗散的更快,在图中λ=0.5 ms 时动能只振荡1 次,λ=1 ms 时动能振荡2 次.

图7 不同松弛时间λ 下动能KE 随时间的变化Fig.7 Kinetic energy KE variation curve with time at different relaxation time λ

图8 为不同松弛时间下液滴表面能随时间的变化曲线.在图中可以看出液滴在挤压过程中界面总长度增加,表面能迅速上升,松弛时间越大时,表面能的峰值越大,随后在液滴回缩过程中表面能快速减小,由于此时初始的两个液滴已经合并为一个大液滴,因此其表面能减小到小于初始时刻表面能,在后续的液滴振荡过程中随着液滴的振荡,相应时刻的表面能也有较小的起伏,对于较大的松弛时间,表面能达到的平衡值相对更小.

图8 不同松弛时间λ 下表面能SE 随时间的变化Fig.8 Surface energy SE variation curve with time at different relaxation time λ

Focke和Bothe[30]进行了剪切变稀流体的液滴碰撞仿真研究,发现由于黏性力只对碰撞最初始阶段产生影响,因此剪切变稀流体的液滴碰撞过程可以被具有相等有效黏度的牛顿流体复现.对于不同黏度的牛顿液滴碰撞合并过程,Sun 等[38]的仿真结果表明,当Oh增长时,动能的振荡频率随之增大,而振幅相差较小,并且碰撞过程中动能的耗散速度与第一次振荡的动能变化曲线在不同Oh下没有显著区别,这与本节得到的改变松弛时间带来的变化有所不同.

2.2 黏度比对黏弹性液滴正撞过程的影响

图9 为不同黏度比β时液滴最大宽度随时间的变化曲线,可以看出β越大时,液滴的碰撞挤压速率越慢,液滴被挤压的程度越小,其达到的最大宽度值越小.随着β的增大,液滴变形的速率越慢,振荡幅度更小,更早地达到平衡,这说明在相同的总黏度下,当聚合物黏度的占比降低,即体系中的黏弹性物质较少时,回缩与振荡的趋势越小.在β=0.2,0.3,0.4 时Dmax的值在回缩后达到平衡,液滴不再有明显的振荡.

图9 不同黏度比β 下液滴最大宽度Dmax 随时间的变化Fig.9 The curve of droplet maximum width Dmax with time at different viscosity ratio β

图10和图11 分别为不同黏度比β时动能和表面能随时间的变化曲线.从图9可以看到对于β越大的液滴,聚合物占比越小,其动能耗散速率越快,后续动能的振幅更小.而从图10可以看出液滴的表面能在β的改变下没有明显的变化规律,但可以看出整体上β更大的液滴具有更大的表面能.

图10 不同黏度比β 下动能KE 随时间的变化Fig.10 Kinetic energy KE variation curve with time at different viscosity ratio β

图11 不同黏度比β 下表面能SE 随时间的变化Fig.11 Surface energy SE variation curve with time at different viscosity ratio β

2.3 韦伯数对黏弹性液滴正撞过程的影响

图12 为不同韦伯数下液滴的碰撞过程.在液滴碰撞挤压阶段,韦伯数的增长有助于液滴的挤压,因此达到的Dmax峰值更大,但由于不同韦伯数下液滴的挤压速率是相等的,所以大韦伯数时Dmax达到峰值的时间延后.由于随着韦伯数的减小,更多的能量被储存为表面能,因此在收缩阶段,更多的表面能释放,导致液滴收缩得更快,在图中可以看到更小的韦伯数更早完成收缩.在振荡阶段,相比于小韦伯数(We=30,50) 的结果,大韦伯数(We=100,200,1000)时融合体的振幅更小.

图12 不同We 下液滴最大宽度Dmax 随时间的变化Fig.12 The curve of droplet maximum width Dmax with time at different Weber number

图13 为不同韦伯数下动能随时间的变化曲线,在动能下降阶段,韦伯数越小时耗散速率越大,因此液滴的初始碰撞动能更早耗尽,开始回缩.在振荡阶段,小韦伯数时的动能更大,因此液滴的振幅越大.图14 为表面能的变化曲线,可以看出韦伯数更大时,表面能越大,其挤压和回缩过程的相对幅值也更大.

图13 不同We 下动能KE 随时间的变化Fig.13 Kinetic energy KE variation curve with time at different Weber number

图14 不同We 下表面能SE 随时间的变化Fig.14 Surface energy SE variation curve with time at different Weber number

2.4 碰撞因子对黏弹性液滴碰撞过程的影响

改变碰撞因子B,液滴偏心碰撞时会得到不同的结果.如图15 是不同碰撞因子B下两液滴的碰撞过程,可以看出当增大B时,液滴在碰撞后会发生拉伸与转动,B越大时拉伸得越长.当B=0和B=0.2 时,液滴在碰撞挤压后融合,并在挤压拉伸过程后较快地恢复成球形;当B=0.4和B=0.6 时,液滴在挤压后继续向着运动方向拉伸,可以拉伸到更远地距离;当B=0.8 时,液滴在挤压后并没有发生融合,而是继续向相反的方向运动,最终依然分离为两个液滴.在偏心碰撞过程中,碰撞惯性力的法向分量用来克服液滴之间的气膜,而切向分量转换成液滴融合体的旋转运动,当增大B时,更强的切向力使融合体被拉伸.

图15 不同碰撞因子B 下液滴的碰撞过程Fig.15 The collision process of droplet under different impact factor B

图16 为不同B下动能随时间的变化,从中观察到B越小时动能耗散的速率越快,并且耗散到更低的动能.当B=0～ 0.6 时,液滴在其碰撞过程中耗散了绝大部分的动能,但B=0.8 时,液滴在互相接近的过程中只耗散了不到60%的动能,这说明液滴的融合过程会耗散更多的动能.图17 为对应的表面能变化过程,由于B=0.4和B=0.6 时的液滴融合体在计算时间内一直处于旋转拉伸状态,所以其表面能曲线还在上升过程中.B=0.8 时液滴在挤压过程中表面能达到最大值,随后在分离后的球形化过程中降低,并逐渐降到平衡值.