可压缩多介质粘性流体和湍流的大涡模拟*

柏劲松,王 涛,邹立勇,李 平

(中国工程物理研究院流体物理研究所冲击波物理与爆轰物理国防科技重点实验室,四川 绵阳621900)

1 引 言

对于复杂湍流的研究,解析方法几乎无能为力,实验测量和数值模拟是研究复杂湍流的必要手段。湍流的数值模拟主要包括直接数值模拟(direct numerical simulation,DNS)、雷诺平均数值模拟(Reynolds averaged Navier-Stokes,RA NS)和大涡数值模拟(large eddy simulation,LES)3 种。其中LES 是基于网格尺度封闭模型及对大尺度涡进行直接求解NS 方程,可以模拟湍流发展过程。近年来,湍流大涡数值模拟普遍接受的观念是:湍流统计特性能否与DNS 或实验结果符合良好是衡量大涡数值模拟准确性的标准,而不是追求它与DNS 结果的相关性[1-2]。大尺度脉动到小尺度脉动的能量传递是大涡数值模拟的关键量,只要亚格子模型正确模拟能量传递,大尺度脉动的统计量就基本上能够计算正确。

亚格子应力模型是湍流大涡模拟的核心问题,Smagorisky 模型就是参照雷诺平均模式最早提出的唯象涡粘模型,至今仍在使用。涡粘模型的最大优点是简单,调整模型系数能够保证模型的亚格子耗散与实际亚格子耗散一致,但是Smagorisky 模型属于耗散型,在流场中任意一点都是从可解尺度湍流向不可解尺度湍流输送能量,即所谓正转,而不存在相反方向的能量传递,在实际复杂湍流中,已经发现可能存在局部正转。Smagorisky 模型总体上耗散过大,因此后来又出现了动力模型[3]、理性亚格子模型[4]等。另一种涡粘模式亚格子应力模型叫Vreman 模型,由A.W.Vreman[5]提出,该模型的耗散较S magorisky 模型小,在低雷诺数湍流混合层计算中的计算结果与动力模型基本一致,更接近于DNS 的计算值。Vreman 模型保持了Smagorisky 模型的优点,在湍流混合层和槽道湍流大涡模拟中得到推广使用[6]。

本文中以可压缩多介质粘性流体界面不稳定性诱发湍流混合研究为目的,在M FPPM[7](multi-fluid piecew ise parabolic method)和M VPPM[8]基础上,利用Smagorinsky 和V reman 亚格子湍流模型,采用大涡数值模拟方法求解可压缩流体NS 方程,通过算子分裂分步计算,以多介质流体高精度(parabolic piecew ise method)方法计算无粘流场,二阶空间中心差方法和两步Rung-Kutta 时间推进相结合的方法计算分子动力学粘性、湍流粘性以及热流部分对流场的影响,研制了适用于可压缩多介质流体界面不稳定性发展演化至湍流阶段的计算方法和二维计算程序M VFT。文中在2 种亚格子湍流模型下采用大涡模拟方法计算了在LANL 激波管上进行的单气柱实验[9],给出了不同时刻气柱的高度和宽度,流场中SF6上游边界、下游边界和涡边界处的速度以及涡的特征。通过与LAN L 实验和计算结果的比较可见,在单气柱激波管界面不稳定性诱发至湍流混合大涡数值模拟研究中,Vreman 模型和Smagorinsky 模型计算的气柱形状基本一致,而在流场速度计算中Vreman 模型的计算结果与实验更接近一些,比较2 种模型计算涡的分布,V reman 模型耗散小于Smagorinsky 模型。通过比较和分析初步确认MVFT 方法和计算程序可用于对界面不稳定性诱发湍流混合流动的大涡数值模拟。

2 控制方程和计算方法

2.1 基本控制方程

考虑热传导和粘性情况下的大涡模拟控制方程采用笛卡尔坐标系下的NS 方程,经过Faver 滤波得到,略去非线性项后的基本形式为

式中:i、j 分别代表x、y、z 等3 个方向,相同的i、j 表示求和分别为密度、速度和压力,为单位质量的总能量,N 为介质的种类,φs 为第s 种介质的体积分数,在含有N 种不同介质的混合网格中,各体积分数满足。状态方程采用p=(γ-1)ρe+γπ形式,对于气体,γ为比热比(γ＞1),π=0,对于液体和用冲击Hugoniot 曲线描述的弹性密实介质,γ为材料的拟合常数,π为具有粘性张量量纲的材料常数。为有效粘性应力张量,为湍流有效粘性系数为流体动力粘性系数,为亚格子粘性系数为热传导在单位时间单位空间的能量流为流体温度为有效导热系数,Prlam+μsgscp/Prsgs,λlam为流体导热系数,λsgs为亚格子导热系数,cp为比定压热容,Prlam和Prsgs为普朗特数为热导率。

2.2 2 种亚格子应力模型

对于Smagorinsky 模型,其亚格子粘性系数

对于V reman 模型,其亚格子粘性系数

2.3 计算方法

采用算子分裂技术,将方程组(1)描述的物理过程分解为3 个子过程进行计算,即整个流量的计算分解为无粘流量、粘性流量和热流量3 部分。那么方程组(1)就可以分裂为如下的2 个方程组

对于无粘流量部分包括对冲击波、稀疏波及接触间断的计算,采用PPM 方法求解方程组(4),对于多物质界面,由于不稳定性后期界面混乱甚至破碎,因此无法进行界面重构,采用VOF(volume of fluid)追踪算法来捕捉,根据介质的体积分数区分界面附近的不同材料,采用压力一致、速度相同基本假定确定混合网格中的计算参量,具体推导过程见文献[8]。对于粘性流量和热流量的计算,考虑湍流的亚格子应力、牛顿流体粘性应力和能量流的影响,所以在求解方程组(4)的基础上,还需要求解方程组(5)。由于亚格子应力、牛顿流体粘性应力和能量流仅对流场的动量、能量和界面扩散产生影响,因此方程组(5)中第1 个方程描述的质量保持守恒,离散求解时可以不予考虑,方程组(5)的第2 ～4 个方程采用二阶空间中心差分方法计算,然后再利用两步Rung-Kutta 方法进行时间推进,求解过程见文献[11]。

3 LANL 激波加载气柱不稳定性实验模型大涡模拟

C.A.Zoldi[9]在LANL 激波管上进行激波加载重气柱实验,实验测试段的布局如图1,1.2Ma 的空气冲击波作用于SF6气柱,由于在空气和SF6气柱界面处存在压力和密度梯度,因此在激波作用下的RM 界面不稳定性发展演化产生涡结构,最后形成空气和SF 6 气体的混合。涡通过速度的旋度进行计算,即ω=▽×v。C.A.Zoldi[9]采用平面激光诱导荧光技术PLIF(planar laser induced fluorescence)给出了实验观察到的6 个时刻SF6 气柱形状,如图2(a)所示,运用粒子图像速度仪PIV(particle image velocimetry)技术给出了SF6气柱3 个测试部位的速度,并应用RAGE 计算程序进行数值模拟,给出了对应时刻的计算图形,如图2(b)所示。RAGE 是求解Euler 方程的计算程序,忽略了气体的粘性和热效应,未考虑湍流模型。本文中根据第2 节方法研制了用于求解NS 方程的多介质粘性流体和湍流大涡模拟计算程序M VFT,既考虑了气体的物理粘性和热效应,还加入湍流模型模拟湍流混合过程,从物理上更适用于模拟激波加载重气柱至混合实验模型。

图1 激波管测试段示意图Fig.1 Test section of the shock tube

图2 实验动态图像和RAGE 模拟结果Fig.2 Experimental dynamic images and RAGE simulations at different times

计算采用的空气和SF6气体参数如表1 所示,普朗特数取为0.67 和0.7。在数值模拟中考虑到SF 6 气柱在空气中的耗散,在SF 6 和空气之间设置了一个高斯函数型的耗散界面过渡层ITL(interface transition layer),如图3 所示,过渡层内的密度运用高斯函数计算,其中Rc=0.12 cm,δ=0.275 cm,计算区域网格为600×100。计算采用Smagorinsky 和Vreman 亚格子湍流模型,分别给出了模型参数CS=0.17 和CV=0.07 时SF6 气柱各对应时刻的密度、涡等位图,如图4、5 所示。定性比较2 模型给出的密度和涡图形,在早期区别较小,只是在后期图形出现明显差异,在750 μs时刻V reman 模型给出的SF 6 气柱形状与相同时刻实验图形更接近一些,为了定量比较2 种亚格子应力模型计算结果以及与实验和RAGE 计算结果,图6 中给出了SF 6 气柱在几个时刻的宽度l 和高度h的比较,从图6 可见,RAGE 的计算结果偏离实验最大,M VFT 所用的2 个亚格子模型计算结果基本一致,均较RAGE 接近实验值,但实验值均大于2 程序的计算结果;对于流场中SF6 气柱运动速度的比较,给出了图7 中定义的U E(upstream edge)、DE(dow nstream edge)和VE(vortex edge)3 个位置的位移历史,图7(a)、(b)分别为Smagorinsky 模型和Vreman 模型的计算结果,表2 列出了3 个位置平均运动速度vU E、vDE、vV E的M VFT 计算值、实验值和RAGE 程序计算值。通过比较可见,M VFT 计算中Vreman 模型(VM)给出的速度值大于Smagorinsky 模型(SM),但均小于RAGE 程序的计算值,MVFT 和RAGE 的计算值均远小于实验值。

表1 空气和SF6 计算参数[9]Table 1 Computational parameters for air and SF6[9]

对于计算和实验差异,C.A.Zoldi[9]认为实验中空气冲击波马赫数略大于1.2 是造成实验值偏大的原因之一,另外如果实验所用的SF6 气柱密度小于计算所取的SF6 气柱密度也会导致实验值偏大;由于MVFT 考虑了气体的物理粘性、热效应以及湍流模型,因此其计算值小于RAG E 的计算值是合理的;对于MV FT 所用的2 个亚格子模型计算的差异,则是由于V reman 模型的耗散小于Smagorinsky模型所引起,通过图8 比较470 μs 时刻经过最大涡中心沿x 方向的分布图可以进一步证实,前者的涡峰大于后者,最左边的负峰值和正峰值则分别对应界面过渡层中上游边界和下游边界位置。

表2 运动气柱速度Table 2 Velocities of the evolving cylinder

图3 高斯函数型耗散界面过渡层Fig.3 Gaussian interfacial transition layer

图4 不同时刻MVFT 计算的密度等位图Fig.4 Density contour images by M VFT at different times

图5 不同时刻M VFT 计算的涡等位图Fig.5 Vortex contour images by M VFT at different times

图6 SF 6 气柱宽度和高度的实验和计算结果Fig.6 Width and height measurements and calculations of the SF 6 evloving cylinder

图7 SF6 气柱上游、下游和涡边界的位置Fig.7 Location and displacement of the upstream edge(UE),the downstream edge(DE),and the vortex edge(VE)on the SF6 evolving cylinder

图8 470 μs 时刻过最大涡中心沿x 方向分布图Fig.8 Vorticity profile alone the x axis through the maximum vorticity at t=470 μs

4 结 论

(1)将Smagorinsky 和V reman 亚格子模型引入到可压缩多介质粘性流体动力学高精度计算程序M VPPM 中,实现了RM 界面不稳定性后期湍流混合的大涡数值模拟,研制了适用于可压缩多介质流体界面不稳定性发展演化至湍流阶段的计算方法和二维计算程序M VFT;(2)分析了Smagorinsky 和V reman 涡粘型亚格子模型对能量的耗散,通过对LA NL 激波管单气柱RM不稳定性实验中气柱上下游边界、涡边界速度的定量计算及与实验和RAGE 程序计算结果比较可见,V reman亚格子模型对能量的耗散小于Smagorinsky 模型;(3)RM 界面不稳定性后期湍流混合大涡数值模拟还存在一些不确定因素,Smagorinsky 模型和Vreman 模型不包含转捩机制,模型选取和模型参数选择具有不唯一性等。

[1] Jimenez J,Moser R D.Large-eddy simulation:Where are w e and w hat can we expect?[J].AIAA Journal,2000,38(4):605-612.

[2] Pope S.Ten questions concerning the large eddy simulation of turbulent flow s[J].New Journal of Physics,2004,6(1):1-24.

[3] 张兆顺,崔桂香,许春晓.湍流理论与模拟[M].北京:清华大学出版社,2005.

[4] Cui G X, Zhou H B,Zhang Z S.A new dynamic subgrid eddy viscosity model w ith application to tubulent channel flow[J].Physics of Fluids,2004,16(8):2835-2942.

[5] Vreman A W.An eddy-viscosity subgrid-scale model for turbulent shear flow:Algebraic theory and applications[J].Physics of Fluids, 2004,16(10):3670-3681.

[6] Park N, Lee S, Lee J, et al.A dynamic subgrid-scale eddy viscosity model w ith a global model coefficient[J].Physics of Fluids,2006,18(2):125109.

[7] Bai J S, Li P,Zhang Z J,et al.Application of the high-resolution Godunov method to the multifluid flow calculations[J].Chinese Physics, 2004,13(12):1992-1997.

[8] 柏劲松,李平,王涛,等.可压缩多介质粘性流体的数值计算[J].爆炸与冲击,2007,27(6):515-521.BAI Jing-song, LI Ping,WANG Tao,et al.Computation of compressible multi-viscosity-fluid flow s[J].Explosion and Shock Waves, 2007,27(6):515-521.

[9] Zoldi C A.A numerical and experimental study of a shock accelerated heavy gas cylinder[D].New York:Applied Mathematics and Statistics, State University of New York at Stony Brook,2002.

[10] Lilly D K.The representation of small-scale turbulence in numerical simulation experiments[C]∥Goldstine H H.Proceedings of IBM Scientific Com puting Symposium on Environmental Sciences.Yorktown Heights,New York,1967:195-210.

[11] 柏劲松,李平,邹立勇,等.界面不稳定性引起混合过程的二维数值计算[J].力学学报,2008,40(4):464-471.BAI Jing-song, LI Ping,ZOU Li-yong,et al.Two-dimensional simulation of the mixing induced by interface instabilities[J].Chinese Journal of Theoretical and Applied Mechanics,2008,40(4):464-471.