基于VOF模型的凝胶推进剂的伪流动研究

曹钦柳,封 锋,邓寒玉

(南京理工大学 机械工程学院,江苏 南京 210094)

基于VOF模型的凝胶推进剂的伪流动研究

曹钦柳,封 锋,邓寒玉

(南京理工大学 机械工程学院,江苏 南京 210094)

利用VOF模型计算凝胶推进剂一次雾化的两相流动时,为研究引入表面张力后伪流动的特征,开发了基于同步重叠相加(SOLA)交错网格算法的非牛顿不可压流体求解器。以凝胶推进剂液滴为对象,研究了VOF模型下的伪流动。结果表明:使用本方法的验证算例中数值解和解析解吻合度较高,误差在5%以内;伪流动强度定义中包含表面张力系数,一定程度上抵消了凝胶推进剂表面张力计算误差增大对伪流动强度的影响;4种网格尺度下,伪流动强度的变化趋势及时间平均值基本一致,说明网格尺度对伪流动无影响;伪流动强度随着本构参数n的升高而增强,并且当n>0.6时上升得更快;伪流动强度与本构参数K呈比例系数为0.05的线性关系。

凝胶推进剂;表面张力;流体体积函数;伪流动

凝胶推进剂作为一种新型的推进剂兼具了固体和液体推进剂各自的优点,具有高安全性、易储藏、不易泄漏和可重复启动的特点,因此得到迅速的发展。要使凝胶推进剂实现高效的燃烧,需要通过相关措施使其充分雾化。对于液滴的二次雾化问题,已经提出相关的振荡和破碎模型[1-2]。但是对于凝胶推进剂的一次雾化,目前仍无较为完备的理论,相关研究者提出了一些近似的计算模型,希望能够从雾化到燃烧进行完整的数值模拟,采用VOF、level set和SPH等方法都取得了较好的效果[3-5]。在液体一次雾化中表面张力是不容忽视的,在数值模拟中引入表面张力模型后,由于不能精确计算表面张力和界面的压力跳跃条件,给计算上带来了数值误差。这种数值误差的直接后果便是在模拟中出现伪流动现象,当伪流动足够强时,会引起数值振荡甚至得到错误的结果。Gueyffier等[6]指出,在一些高密度比和高黏性比的情况下,很多VOF方法都不能按照真实的参数进行计算。

对于伪流动现象国内外研究者进行了大量研究工作。Lafaurie等[7]基于VOF方法对气泡的伪流动进行了研究,提出了一种衡量伪流动的参数,结果表明牛顿流体中气泡的表面张力越强时伪流动就越明显。Popinet和Zaleski[8]用一种front-tracking方法细化了网格压力的梯度离散,较好地处理了表面张力和界面压力跳跃平衡,使得伪流动强度大大降低。Aleinov等[9]基于核函数提出了计算表面张力的高阶内核概念,使得更精确地计算界面曲率成为可能。马东军等[10]对VOF方法的伪流动进行了研究,得出其强度几乎为常值,并且不随网格加密而收敛。不过上述文献都是针对牛顿流体进行的研究,目前针对凝胶等非牛顿流体伪流动现象的研究还未见公开报道。

本文采用VOF模型进行界面捕捉,利用CSF模型处理表面张力,开发了求解不可压流体控制方程组的求解器,并以凝胶推进剂液滴为对象,对伪流动现象进行了研究。分析了不同表面张力系数,不同网格尺度和不同幂率模型本构参数对伪流动强度的影响。相关结论可以为采用VOF方法研究凝胶推进剂的一次雾化提供参考。

1 数值方法

1.1 计算模型

本文研究凝胶推进剂液滴在VOF模型下的伪流动,液滴的直径、计算区域大小及其边界类型如图1所示。凝胶液滴周围为空气,四周为自由边界,速度零梯度外推。初始时压力设为环境压力101 325 Pa,速度均为0。

1.2 界面重构方法

在两相共存的网格中,实际上是有界面的。采用施主-受主思想计算时,需要几何重构方法构造界面。研究人员提出了一些几何重构方法,如SLIC、FLAIR和PLIC等[11]。其中PLIC方法是由Youngs首先提出的。它通过单个网格的斜线来近似两相界面,斜线的斜率需要通过周围网格确定,在重构精度上可以达到二阶。因此,本文采用PLIC重构方法进行数值计算。如图2给出了界面法向在第三象限的4种情况。

1.3 黏性模型

凝胶推进剂在流动上表现出非牛顿流体的性质,其流变本构方程可写为如下形式[12]:

(1)

(2)

(3)

Δ可以根据剪切速度张量各分量求得:

(4)

由上述公式,可以得到应力张量各分量表达式为

(5)

(6)

(7)

式中:η0为表观黏度,它是衡量非牛顿流体黏性的重要参数,可以表示为

η0=KΔ(n-1)/2

(8)

1.4 表面张力模型

计算表面张力时,本文采用Brackbill提出的连续表面张力模型(CSF),假设表面张力为一种体积力,并且连续分布于两相界面两侧的一定厚度内。其大小与界面的曲率成正比,则网格内两相界面受的表面张力为[13]

fs=σκ(F)

(9)

式中:σ为表面张力系数,κ为两相界面的曲率,F为网格中液体相的体积分数。上式求得的表面张力实际是一种压力对坐标的变化量,并作用于网格内所有的流体,它作为源项添加到动量方程中。

在牛顿流体的VOF伪流动研究中Popinet[8]等引入了一个衡量伪流动强度的参数,其表达式为

Cas=|U|μ/σ

(10)

式中:μ为牛顿流体的黏性系数,U为速度矢量。目前还未提出专门针对凝胶伪流动的衡量参数,因此本文尝试直接用表观黏度η0代替上述黏性系数。

表面张力的引入给计算带来了一定的数值误差。采用重构方法获得的两相界面始终是对真实界面的一种近似,导致表面张力的计算值与实际值存在一定误差。目前的数值方法很难处理表面张力带来的压力跳跃平衡条件,二者带来的误差会驱动流场的非物理变化,这是产生伪流动的根本原因。需要指出的是,伪流动是计算模型本身的固有特性,提高重构精度或者采用更好的数值方法可以减小伪流动的强度。

2 数值验证

本文采用基于交错网格的SOLA算法[14],压力和密度值存储于格心处,速度值则保存在网格线上。压力项和黏性项采用中心差分格式,对流项采用不完全迎风差分格式,即迎风格式与中心差分格式的线性组合,本文的线性系数取为0.5。为验证本文数值方法的有效性,选取如图3并行两相通道流作为验证算例,图中上下为两平行平板。

3 计算结果与分析

本文为深入探讨VOF模型中凝胶液滴的伪流动现象,对不同表面张力系数,网格尺度和幂率模型本构参数下凝胶液滴的伪流动进行了数值计算。空气密度取1.25kg/m3,凝胶推进剂的密度在不同的配比下,其值有所差异,本文直接取为1 000kg/m3。按照幂率模型,当剪切速率为0时,其表观黏度趋向于无穷大,因此数值计算时需要对上界值进行限定,本文取最大表观黏度值ηmax为10Pa·s。

3.1 表面张力系数的影响

Lafaurie等[7]基于VOF方法对牛顿流体气泡的伪流动进行了研究,指出表面张力增大时伪流动越明显。表面张力系数是表面张力中的一个重要参数。为研究伪流动强弱随表面张力系数的变化规律,本文选取0.01,0.02,0.03,0.04,0.05,0.06和0.07 7种不同的表面张力系数值,采用幂率型煤油凝胶(K=13.5,n=0.47),网格数为2002,对其进行数值计算。

图5给出了5ms时4种张力系数下的凝胶液滴表面图。从图中可以看出,随着表面张力的升高,界面扭曲程度增大,说明表面出现了较大的伪流动速度。为定量描述伪流动强度Cas值的变化情况,图6给出了3个时刻伪流动强度随表面张力系数的变化规律。从图中可以看出,表面张力系数增大时,3个时刻的伪流动强度稍有下降或者上升,但是其变化都比较小,基本保持一致。从图5可以看出,表面张力系数较大时表面处伪流动带来的影响明显强于系数较小的情况。虽然表面张力系数大时,造成的界面处伪流动速度大,但同时伪流动强度公式(10)中也包含了表面张力系数的影响,从而一定程度上抵消了表面张力误差对伪流动强度的影响。

3.2 网格尺度的影响

数值上计算两相界面曲率时存在较大误差是造成表面张力误差的原因。如果采用更加精细的网格,细化对界面的描述,就会对界面曲率的计算精度产生影响。为研究凝胶液滴的伪流动强度随网格尺度的变化规律,选用相同的幂率型煤油凝胶,在不同的网格尺度下,对其进行了数值研究。网格的相关参数如表1所示,表中m为网格数量,md为液滴直径上网格点的数量。

表1 网格参数

图7为4个算例10ms时的凝胶液滴表面示意图。从图中可以看出,网格数目为502的算例发生了明显的扭曲,但是界面比较平缓且褶皱较少。随着网格数目的增加,两相界面都显著偏离了初始的圆形界面,并且界面的褶皱逐渐增多,界面的锐利性逐渐增强。网格数目的增加并未有效延缓伪流动的发生,并且网格数目增加时界面将在更小的范围内发生扭曲,反而使得褶皱逐渐增加。因此,网格数增加可以使得界面在更小的范围内变化,但是这种变化终会逐渐发展到可观察的程度。

马东军[10]使用VOF模型研究了牛顿流体伪流动现象随网格尺度的变化,结果表明Cas值不随网格尺度的变化而变化,说明伪流动强度不具备网格收敛性。图8给出了4种网格尺度Cas值在10ms内随时间的变化图。从图中可以看出,4种情况Cas值都随时间逐渐升高,伪流动都随时间逐渐增强,并且其变化规律基本一致,其时间平均值分别为0.233,0.221 8,0.222和0.221,时间平均值的标准方差为0.004 2。所以,加密网格后计算量大了,但是对伪流动的时间增长率和强度都未产生多大影响,本文采用的数值方法在凝胶的伪流动中同样不具备网格收敛性。因此,通过加密网格来控制伪流动,对提高含表面张力计算的数值精度是无效的。

3.3 本构参数的影响

在针对牛顿流体的伪流动研究中,一些无量纲参数是影响伪流动强度的重要因素,如凝胶液滴与空气的密度比ρl/ρg和黏性比μl/μg,以及衡量表面张力和黏性力相互关系的Ohnesorge数[8],但是,对于凝胶推进剂,其表观黏度随剪切速率而变化,流场中每个网格的表观黏度都有所差别。本文算例中空气和凝胶液滴的密度均为常数,密度比为定值,而黏性比和Ohnesorge数没有统一的数值,所以很难用这两个参数来研究黏度对伪流动的影响。式(1)的本构参数K和n体现了黏度的变化,因此本文选取不同的K和n值,如表2,网格数为2002,对凝胶液滴的伪流动进行了计算。

表2 凝胶推进剂的本构参数

为定量描述伪流动强度的变化,图9给出了Cas值随n的变化规律。从图中可以看出,随着n的增大,伪流动的强度逐渐增大,最大值5.18比最小值0.043要高出2个数量级,并且当n>0.6时,伪流动强度的变化对n很敏感。这说明,虽然在目前的算法中伪流动是不可避免的数值误差,但是对于不同的流体,其强度是有差别的。非牛顿流体稀化能力越强,相同时间内由伪流动带来的误差就越小。所以,对于稀化能力较弱的凝胶推进剂,采用包含表面张力的VOF模型进行计算时,其结果需要谨慎对待。

Lafaurie等[7]采用基于VOF的新型算法,得出牛顿流体伪流动强度随黏性的增大而降低。文献[10]采用VOF和levelset模型进行计算,结果表明Cas值几乎为常数,与黏性无关。K值体现了凝胶推进剂的黏度平均水平,图10为Cas值随K值的变化规律。Cas随K值几乎成系数为0.05的比例增加。凝胶推进剂具有较高的黏度,对Cas值贡献较高。流场速度小,其差异未能体现出来。所以Cas值与凝胶推进剂的平均黏度呈现出了明显的相关性。

4 结束语

表面张力系数增大时,两相界面处的伪流动速度增大,界面扭曲程度增大,但是伪流动强度值变化不大,因为伪流动强度定义中同样包含了表面张力系数的影响,从而抵消表面张力系数的影响。

VOF模型下4种网格精度的Cas值变化趋势和时间平均值基本一致,其时间平均值的标准方差仅为0.004 2,网格精度对伪流动的强度无影响。采用加密网格不能减小小伪流动带来的数值误差。

伪流动强度随着凝胶推进剂本构参数n值的增加而增强,并且当n>0.6时增加的梯度更大。伪流动强度最大值5.18比0.043要高出2个数量级。

伪流动强度随凝胶推进剂的平均黏度度量K值的增加而增强,并且几乎成系数为0.05的比例增加。

[1] 刘晨.复杂燃烧流场数值模拟方法研究[D].南京:南京航空航天大学,2009. LIU Chen.Numerical methods for complex combustion flowfields[D].Nanjing:Nanjing University of Aeronautics and Astronautics,2009.(in Chinese)[2] 文华.基于CFD的柴油机喷雾混合过程的多维数值模拟[D].武汉:华中科技大学,2004. WEN Hua.Multi-dimensional numerical modeling of spray mixing process in diesel engines based on CFD[D].Wuhan:Huazhong University of Science and Technology,2004.(in Chinese)

[3] MA D J.Atomization patterns and breakup characteristics of liquid sheets formed by two impinging jets:AIAA 2011-97[R].USA:AIAA,2011.

[4] ARIENTI M,LI X,SOTEROIU M C,et al.Coupled level-set/volume-of-fluid method for simulation of injector atomization[J].Journal of Propulsion & Power,2013,29(1):147-157.

[5] 强洪夫,刘虎,陈福振,等.基于SPH方法的射流撞击仿真[J].推进技术,2012,33(3):424-429. QIANG Hong-fu,LIU Hu,CHEN Fu-zhen,et al.Simulation on jet impingement based on SPH method[J].Journal of Propulsion Technology,2012,33(3):424-429.(in Chinese)

[6] GUEYFFIER D,LI J,NADIM A,et al.Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows[J].Journal of Computational Physics,1999,152(2):423-456.

[7] LAFAURIE B,NARDONE C,SCARDOVELLI R,et al.Modelling merging and fragmentation in multiphase flows with surfer[J].Journal of Computational Physics,1994,113(1):134-147.

[8] POPINET S,ZALESKI S.A front-tracking algorithm for accurate representation of surface tension[J].International Journal for Numerical Methods in Fluids,1999,30(6):775-793.

[9] ALEINOV I,PUCKETT E G.Computing surface tension with high-order kernels[C]//Proceeding of the sixth international symposium on computational fluid dynamics.Lake Tahoe:NV,1995:13-18.

[10] 马东军.可压缩/不可压缩流体交界面高精度数值方法的研究[D].合肥:中国科学技术大学,2002. MA Dong-jun.Numerical investigation of the high resolution interfacial simulation methods applied to compressible and incompressible flows[D].Hefei:University of Science and Technology of China,2002.(in Chinese)

[11] 赵大勇,李维仲.VOF方法中几种界面重构技术的比较[J].热科学与技术,2003,2(4):318-323. ZHAO Da-yong,LI Wei-zhong.Comparison of three kinds of interface reconstruction methods applied to VOF[J].Journal of Thermal Science and Technology,2003,2(4):318-323.(in Chinese)

[12] RAHIMI S,NATAN B.Numerical solution of the flow of power-law gel propellants in converging injectors[J].Propellants Explosives Pyrotechnics,2000,25(4):203-212.

[13] RAESSI M.On modeling surface tension-dominant,large density ratio,two-phase flows[D].Toronto:University of Toronto,2008.

[14] 张德良.计算流体力学教程[M].北京:高等教育出版社,2010:367-374. ZHANG De-liang.Computational fluid dynamics tutorial[M].Beijing:High Education Press,2010:367-374.(in Chinese)

Research on Spurious Currents of Gel Propellant Based on VOF Model

CAO Qin-liu,FENG Feng,DENG Han-yu

(School of Mechanical Engineering,Nanjing University of Science and Technology,Nanjing 210094,China)

In order to study the characteristics of spurious currents considering surface tension in the two-phase flow for gel propellant primary-atomization,an incompressible fluid solver was developed based on synchronous overlap-add(SOLA)staggered mesh.Based on the gel propellant droplet,the spurious currents of VOF method was investigated.The results indicate that the numerical solution using this method agrees well with the analytical solution,and the relevant error is within 5%.Surface tension coefficient is taken into account in the definition of the strength of spurious currents,which partly eliminates the effect of the errors of surface tension of gel propellant on spurious currents values.The four mesh scales has nearly similar trends and time-average values of spurious currents,which reveals little effect of mesh scale on spurious currents.Spurious currents improves with the increase of the indexn,and the growth rate is greater whenn>0.6.The values of spurious currents linearly varies with the indexK,and the proportionality coefficient was 0.05.

gel propellant;surface tension;VOF;spurious currents

2016-12-22

航天科技创新基金(CASC03-02);中央高校基本科研业务费专项基金(30920140112001);江苏省普通高校学术学位研究生科研创新计划(KYLX16_0474)

曹钦柳(1993- ),男,博士研究生,研究方向为航空宇航推进理论与工程。E-mail:115101000190@njust.edu.cn。

V513

A

1004-499X(2017)02-0065-05