可压缩均匀湍流中重粒子运动行为的先验研究1)

宾远为 武 琦 夏振华 史一蓬,3)

* (北京大学湍流与复杂系统国家重点实验室,北京 100871)

† (浙江大学流体工程研究所,杭州 310027)

引言

携带粒子的湍流问题常见于自然界和工业工程应用中,例如云层中雨滴的形成[1-2],大气中污染物的扩散[3-4],发动机中的燃烧、雾化和混合等[5-7].由于湍流是一种多尺度的运动现象,粒子在湍流中的运动受到不同尺度流动的影响.然而,尽管已经进行了大量理论、实验和数值模拟工作,但对其中物理机制的认识仍然不够透彻[8-11].

携带粒子的湍流数值模拟方法可以分为全解析(fully resolved) 模拟和点颗粒 (point-particle) 模拟[8-9].这两类方法的分类依据是粒子尺度相对于流场特征尺度(例如Kolmogorov 尺度或大涡模拟中的滤波尺度) 的大小.全解析模拟要求流场的求解特征尺度(即网格尺寸)小于粒子尺度.在全解析模拟中,粒子的尺度是不可忽略的,因此它们被视为运动的壁面边界条件,从而可以直接模拟粒子周围所有尺度的流动[10].湍流中粒子的运动由周围流体施加的作用力和其他外力(如重力)共同决定.同时,流体的运动也受到粒子的反作用力影响.因此,流体与粒子之间的相互作用是一个双向耦合的过程.在某些情况下,当粒子之间距离特别小时,也需要考虑粒子之间的相互作用[12].然而,全解析模拟需要用小网格解析粒子及其周围所有尺度内的流动,因此计算成本较高,且对网格质量要求较高.目前全解析模拟只能实现对单个粒子或数千个粒子的模拟[13-14].而对于实际情况中更常见的百万、千万量级的粒子运动求解,全解析模拟仍然面临较大困难.

点粒子模拟是一种基于拉格朗日观点的方法,其考虑的粒子尺度通常远小于流场特征尺度.因此,可以合理地将粒子视为流动中的质点,并利用点粒子运动模型追踪每个点的加速度、速度和位置[11,15].在点粒子模拟中,流体相的求解采用欧拉观点的场方法,因此该方法也被称为欧拉−拉格朗日法[15].Tsai[16]和Maxey 等[17]理论上给出了点粒子的受力模型,随后被广泛应用于模拟湍流中点粒子的运动[18-20].在欧拉−拉格朗日法中,粒子与流场的耦合方式是关键因素.根据离散相的体积分数(volume fraction)Φv和质量分数(mass fraction) Φm的不同,粒子与流场之间的耦合方式也不同[21].当在气固系统[22]中Φv≤10−6或液固系统[23]中 Φv≤10−4时,离散相对流体相的作用和离散相之间的相互作用可以被合理忽略.数值模拟只需考虑离散相粒子在流场作用力下的运动.这被称为单向耦合(one-way coupling)[8,24].然而,当 Φv和 Φm相对较大时,需要同时考虑流体相对离散相的作用力和离散相对流体相的反作用力,即双向耦合(two-way coupling).随着 Φv和 Φm的进一步增大,还需要考虑离散相之间的相互作用,即四向耦合(four-way coupling).相比于全解析模拟,点粒子模拟对计算量的要求更小,因此在大规模粒子运动求解方面具有优势.目前,点粒子模拟方法已广泛应用于湍流直接数值模拟和大涡模拟中[8-9,25].

直接数值模拟(direct numerical simulation,DNS)是一种精确的湍流模拟方法.欧拉−拉格朗日点粒子模型已被成功应用于一些颗粒湍流相互作用问题的DNS 研究中.例如,Maxey[18]使用单向耦合的点粒子模拟计算粒子在均匀湍流中的沉降.Squires等[26]研究了不可压缩各向均匀同性湍流中粒子的惯性聚集,以及粒子对湍流的反作用[27].Bec 等[28]着重研究了重粒子在各向均匀同性湍流中的加速度概率分布并探究了不同Stokes 数(S t)对加速度的影响.Xia 等[29]进一步在可压缩均匀湍流中研究了重粒子对湍流场的调制作用.然而,DNS 本身计算成本高昂且计算需求苛刻,目前的欧拉−拉格朗日点粒子的DNS 仍受限于低雷诺数和简单流动场景[9,30].

在湍流问题的数值研究中,大涡模拟(largeeddy simulation,LES)方法近年来得到快速发展,逐渐在学术界和工程应用中占据重要地位.LES 通过直接求解湍流中的大尺度运动,而小尺度运动对大尺度运动的影响则通过亚格子(sub-grid stress,SGS)模型来表示.由于不需要解析小尺度运动,LES 的计算量远小于DNS[31].尽管单相湍流的SGS 模型近些年来得到不断完善,但是否可以直接应用于含颗粒的两相湍流问题仍值得进一步研究.在LES 中,由于只能模拟大尺度运动对应的流场,原先存在于DNS流场中的小尺度结构无法被解析,湍流场中的小尺度运动与粒子的相互作用也被忽略.然而,这种忽略小尺度运动的合理性尚未完全验证,因此需要进一步探索.

为了研究这一问题,一个直接方法是考虑粒子在滤波后DNS 流场(filtered DNS,fDNS)中的运动.在不可压缩均匀各向同性湍流(incompressible homogenous isotropic turbulence,IHIT)中,已经进行了相关研究.例如,Fede 等[32]研究了不同滤波宽度下滤波后的大尺度流动中粒子对的运动,以研究粒子碰撞和聚集性.研究结果表明,只有当滤波宽度足够大,以至于流体动能明显下降时,粒子的聚集性才会受到影响.Jin 等[33]进一步研究了LES中SGS 模型对于粒子在湍流中统计特性的影响.他们的计算分析表明,小尺度流动对低密度粒子碰撞、聚集的影响是不可忽视的.因此,他们提出对低密度粒子的LES 必须要考虑SGS 模型,否则可能导致高达60%的相对误差.Ray 等[34]选取了多个滤波宽度,大致涵盖了整个惯性子区范围,主要研究粒子径向分布函数(RDF)和粒子径向相对速度随滤波宽度的变化.他们发现低密度粒子的RDF 会随着滤波宽度变大而降低,而高密度粒子的趋势相反.这一结果验证了Jin 等[33]的结论,即不能忽略小尺度流动对粒子在湍流中运动的影响.Pozorski 等[35]在综述中对LES 滤波对粒子聚集的影响进行了回顾,并总结了不同学者提出的可提升粒子计算准确性的SGS 模型.在对粒子运动统计量的研究方面,Lalescu 等[36]研究了示踪粒子在滤波流场中与加速度有关的统计量.示踪粒子可以视为密度趋近于零的惯性粒子.他们发现滤波后的流场会使粒子加速度的间歇性减弱.这表现为,在滤波后流场中粒子加速度较未滤波的流场中的加速度变化更为平缓.

尽管目前已有许多关于重粒子在可压缩湍流中的物理机理研究[29,37],但关于LES 中滤波宽度的影响还没有得到系统研究.因此,本文以含重粒子的可压缩均匀各向同性湍流(compressible homogenous isotropic turbulence,CHIT)为研究对象,通过先验分析研究不同尺度的湍流结构和粒子的初速度对粒子的聚集行为以及后续稳态统计结果的影响.

1 数值细节

1.1 流体相控制方程

三维的可压缩Navier-Stokes 方程为

其中,ui,ρ,p和T分别是无量纲速度分量、质量、压强和温度.Fi是施加在流体上的大尺度力,Λ 是大尺度的冷却函数.黏性应力 σij和总能E被分别定义为

其中 θ=∂uk/∂xk是无量纲的速度散度.分子黏性 µ 和热导率 κ 基于常数普朗特数假设[38],由Sutherland's律[39]给定

在上述无量纲方程中,我们使用了特征物理量L,U,ρ0,ρ0RT0,ρ0U2和c0≡对长度、速度、密度、压强、能量和声速进行无量纲化.分子黏性µ和热导率 κ 由 µ0和 κ0无量纲化.这里,比热γ ≡Cp/Cv为1.4.Cp,Cv和R≡Cp−Cv分别是定压比热容.定容比热容和气体常数.流动由3 个相对参数控制,它们分别是相对雷诺数Re=ρ0U/L,相对马赫数M=U/c0和相对普朗特数Pr=µ0Cp/κ0.在本研究中,这3 个参数分别为Pr=0.7,Re=200 和M=0.45.参数 α 的定义为 α=PrRe(γ−1)M2.

1.2 粒子控制方程

当粒子直径d远小于Kolmogorov 长度尺度 η,且粒子密度 ρp远大于流体密度 ρ 时,粒子位置xp,i和速度vp,i由下式决定[16-17,29]

粒子控制方程(式(9)和式(10))是基于牛顿不可压缩流体得到的.能否将这一方程用于可压缩流动,需要特别讨论.Loth[42]指出当由粒子与当地流体的相对速度和当地流体声速定义的粒子相对马赫数时,粒子的控制方程必须要考虑可压缩性的影响.在本文中,即使对于密度最大的粒子(S t=10),平均的粒子相对马赫数也仅约等于 0.2.所以本文是可以合理地忽略可压缩效应对粒子控制方程的影响.

1.3 数值方法

在本研究的数值模拟中,计算域为长宽高均为2π的立方盒子,3 个方向均为周期性边界条件.模拟采用了 2563的均匀网格.其中空间离散格式采用了混合的紧致加权本质无振荡 (WENO) 格式[38].该混合格式结合了光滑区域的八阶紧致有限差分格式[43]和激波区域的7 阶WENO 格式[44].时间推进选取3 阶龙格−库塔格式.计算时间步长由Δt=CFLΔx/max(uiui)动态决定.在我们的研究中,网格大小 Δx等于2π/256 ≈0.024 5,CFL数设置为 0.05.流动由在大尺度上施加的剪切力 Fi驱动.我们在计算中假设粒子的直径为d=0.1Δx.在本工作中,Kolmogorov 尺度η ≈Δx,粒子的直径d≈0.1η,这使得粒子直径d远小于Kolmogorov 长度尺度 η 的假设成立.粒子的密度由Stokes 数的定义给定.另外百万颗粒对应的体积分数 Φv约为 3 ×10−5.对于S t=10 的算例,颗粒的质量分数约为0.5.从严格意义上讲,本文所涉及的参数已经略微超出单向耦合成立的范围,颗粒可能会对流场有一定的反馈作用[24].但由于本文聚焦于流场对颗粒运动的影响,因此选择忽略该反馈作用[40].后续将进一步研究该反馈作用的影响.粒子位置和速度的时间推使用二阶龙格−库塔格式.粒子的时间步长 Δtp与流动的时间步长 Δt相同.第i个粒子位置的流场物理属性值在光滑区域或激波区域分别由线性插值或非线性插值得到.

在当前新农村发展过程当中，体育相关产业的发展也发挥着非常重要的影响。以往环境下，受限于农村地区经济发展水平，农民群众的消费意识相对有限，不具备良好的消费习惯和消费意识。内需不足对产业发展是十分不利的。

当流动达到统计稳态后,粒子St数为1 时,在图1 中绘制了某一时刻流场的无量纲散度 θ/θrms,无量纲涡量大小 |ω|/|ω|rms的等值图和粒子位置.由于本文中的瞬时湍流马赫数Mt=urms/c0在0.987～0.989 之间,可以在图1(a) 观察到大量的小激波(Shocklets,θ/θrms＜−3).这些激波对粒子的运动会起到明显影响[40].图1(b) 显示了无量纲涡量大小|ω|/|ω|rms和粒子的对应位置.图中可以看出粒子局部聚集于低涡量的区域内.

为了验证时间步长度的收敛性,我们将 CFL 数分别调整为0.05,0.025 和0.012 5.因此,在每种对应的情况下,时间步长变为 Δt,Δt/2 和 Δt/4.验证算例包含了不同的S t数,即St=1,5,10,并使用I.C.4作为初始条件.之后,随机选择一个粒子并跟踪它的路径线,如图2 所示.很明显,对于相同的S t数,不同时间步长的路径线是相同的.这意味着当CFL=0.05时,时间步长已经收敛.

图2 一个随机选择的粒子在流动中的路径.初始速度采用I.C.4Fig.2 Path of a randomly selected particle.Here,the initial velocity is I.C.4

1.4 初始条件

本节介绍模拟中的初始条件的生成过程.在无粒子情况下,首先通过DNS 模拟得到了达到统计定常状态的可压缩均匀各向同性湍流.为了确保湍流达到统计定常状态,流场演化了5 个大涡翻转时间(TE)以保证湍流达到统计定常状态.然后,该流场被保存,并作为模拟粒子运动的初始流场.总共Np=1.0×106个粒子以不同的给定初始速度随机撒落在计算域中.为了避免粒子初始位置对结果的影响,这些粒子的位置也被保存用于后续的模拟.因此,在本研究中的所有算例中,都采用相同的初始流场和粒子位置.本文中粒子有5 种给定的初速度,具体见表1.具体来说,对于I.C.1,粒子被静态地放置在流场中.对于I.C.2～I.C.4,分别给每个粒子设置对应的给定速度作为初始速度.在I.C.5 中,粒子的初始速度设定为当地的流动速度.

表1 粒子的初始速度条件Table 1 Initial velocity of particles

1.5 粒子聚集行为表征方法

本文使用了两种方法对粒子聚集行为进行表征,分别为概率密度函数(PDF)方法和香农熵(Shannon entropy)方法.其中,PDF 方法适用于对粒子相关物理量进行观察,尤其在统计稳态时较为有效,但无法追踪粒子聚集程度随时间的演化.而香农熵方法则能够追踪粒子聚集时间的演化.香农熵方法将整个计算域统一划分为Nbin=1283个小盒子.在本文中,计算域为边长为 2π 的正立方体,所以,经过划分,会在计算域中得到 1283个大小相同的边长为2π/128的正立方体.香农熵S(t)由以下方程定义

其中 max[H(t)]=lnNbin,p(k,t) 是在时间t时,在第k个小盒子中发现一个粒子的频率,Np(k,t) 是在时间t时,第k个小盒子中的粒子数.当粒子足够均匀分布时,S(t)趋近于1.而当所有粒子都集中在一个小盒子中时,S(t)等于0.在本文的算例中,在初始时刻S(0)约等于0.93.根据该定义,可以使用任一瞬时流场来计算香农熵S(t).

2 不同滤波宽度的影响

在本节中,我们研究不同滤波宽度的影响.由于湍流中存在非线性效应和黏性效应,湍流中涡的尺度从大尺度L到Kolmogorov 长度尺度 η 是连续分布的.同时,流体动能从最大的涡传输到最小的涡.这个过程被称为能量级串.在 CHIT 中,能量谱E(k)与波数k的− 5/3 幂律已被大量文献[38,45]证明.

图3 中实线显示了达到统计定常状态的能谱E(k).由于网格分辨率的限制,我们的DNS 中最大可识别波数是kmax=128.为了对流场进行空间滤波,本文采用了经典的谱截断滤波器,其数学形式为

图3 流动能量谱 E(k) (实线)和5 个不同的截断波数(虚线).从左到右分别是kmax/kc=2,4,8,16,32Fig.3 Flow energy spectrum E(k) (solid line) and five different cut-off wavenumbers (dashed line).From left to right,they are kmax/kc=2,4,8,16,32,respectively

当CHIT 达到统计定常状态时,粒子被均匀随机地撒落到计算域中.在湍流中运动的粒子最终会达到统计上的稳定状态,并会显示出局部聚集的空间分布,这被称为惯性聚集现象[26].在上述过程中,存在两个目前尚不清楚的问题.首先,从初始时刻到粒子达到统计稳定状态的时刻,不同尺度的涡对粒子聚集的影响尚不明确.其次,当达到统计稳定状态时,粒子聚集的程度受到kc和S t的影响程度也不明确.

图4 显示了粒子以I.C.2 为初始速度条件时,随时间演化的香农熵S(t) 在不同的S t和kc情况下的变化.通过图4,可以得出以下4 点结论: (1) 图4 中的S(t) 首先减小并趋于常数.这表明粒子撒落在计算域后会逐渐聚集,直到达到统计定常状态.(2) 通过比较粒子直接由DNS 流场驱动的结果,我们观察到在统计定常状态时,低密度粒子(S t=1)的聚集程度(S(∞)≈0.88)比高密度粒子(St=5,10)的聚集程度(S(∞)≈0.91,0.92) 更强.(3) 当kc从kmax/kc=2 变化到kmax/kc=8 时,S(t) 并没有表现出明显的差异,这表明,当kc选取在合适的范围内,使用LES 是能够得到可信的粒子聚集性统计量的.(4) 如果kmax/kc大于8,随着kmax/kc增加,S(t) 下降的速率更快.此外,当kmax/kc增加,在统计定常状态时,低密度粒子的聚集性降低,而高密度粒子的聚集性增加.因此,小尺度涡有利于低密度粒子的聚集.相反,当S t较大时,例如St＞5,小尺度涡对高密度粒子聚集性的影响不大.该结论与文献[32]中的研究结果一致.

图4 S(t) 随着时间的演化.初始速度采用I.C.2Fig.4 S(t) as a function of time.Initial velocity is I.C.2

从粒子空间数密度的概率密度函数(PDF)可以得到相同的结论.数密度的获取首先可以将计算域划分为若干个立方体,在粒子达到统计定常状态后,可以得到每个立方体中的平均粒子数n,称为平均数密度[26].平均数密度n为零的概率越大,代表粒子的聚集性越强.本文将整个计算域划分为 2563的立方体,并统计达到稳态后关于n的PDF.图5 绘制了在图4 中所有算例的PDF,并证明可以得到与之前相同的结论.具体地,图5 显示出达到稳态后,在低密度粒子情况下,增加kmax/kc会减弱粒子的聚集性.而在高密度粒子情况下,增加kmax/kc会增加粒子的聚集性.

图5 平均数密度 n 的PDF.初始速度采用I.C.2Fig.5 PDF of the average particle number density n.Initial velocity is I.C.2

图6 与图5 相同,但是为粒子绝对加速度的PDFFig.6 Same as Fig.5,but for the magnitudes of the particles’acceleration

图7 与图5 相同,但是为粒子绝对速度的PDFFig.7 Same as Fig.5,but for the magnitudes of the particles’ velocity

3 初始速度影响

在不同的流动场景模拟中,流体中粒子的初速度由不同的情况给出.例如,初始时刻静止的颗粒沉降[46](I.C.1),含颗粒的射流对应[47](I.C.2～I.C.4),和示踪粒子(I.C.5).关于这些初始条件的详细信息已经在表1 中列出.虽然许多研究已经使用了这些初始条件,但在作者了解的范围里,初始条件的影响还没有得到充分的研究.在本节中,将使用香农熵和平均速度来分析这些初始条件对粒子聚集行为的影响.为了验证LES 的适用性,本文选择了一个合适的截断波数,即kmax/kc=4,该值在 CHIT 的LES 模拟中被广泛使用.

粒子平均速度的变化表明了粒子的运动学过程.如图8 所示,对于各种初始条件,平均速度在经过初始的弛豫时间后会趋于相同的值.且该值会随着S t的增加而增加.此外,与传统观点相反,我们发现在粒子初始时刻等于当地流动速度的情况下(I.C.5),其平均速度需要较长的时间才能收敛到统计定常值.这里的统计定常值为不同初速度的算例达到同一平均速度的对应速度值.在图8 中可以看I.C.5 对应的平均速度是最晚与其他算例一致的.值得一提的是,速度达到统计定常值并不意味着粒子运动已经达到统计定常状态.例如,图8 (a) 所示,当t/TE=0.5,S t=1时,对于所有初始条件,其平均速度值已经收敛,但香农熵仍然还未收敛.同样的现象也可以在S t=1 和 5 时观察到.

图8 粒子数平均速度随时间的演化.驱动流场采用kmax/kc=4Fig.8 Averaged particle velocity as a function of time.Flow field is filtered by kmax/kc=4

从图9 所示,可以观察到不同初始速度条件导致聚集性的不同演化.对于所有不同的S t,粒子初速度为I.C.1～I.C.3 时得到的S(t) 的演化基本相同.然而,当初速度过大时(I.C.4),S(t) 的演化与初速度较小的情况会有一些差异.此外,当粒子初始速度为当地流动速度时(I.C.5),粒子较快达到统计稳定状态.最后,可以得出结论,不同的初始条件不影响统计定常状态下粒子的聚集程度.

图9 S(t) 随时间的演化.驱动流场采用kmax/kc=4Fig.9 S(t) as a function of time.Flow field is filtered bykmax/kc=4

通过前面的讨论,本文可以得出以下结论,当kmax/kc=16且粒子初速度采用 I.C.1 时的情况下,S(t)的演化没有明显差异.但是,对于其他初始条件,该结论是否适用仍然是一个问题.图10 展示了当kmax/kc=16 时 对应I.C.1～I.C.5 的S(t) 的演化.与图9相比,对于我们给定的所有初始条件,粒子的聚集性演化过程没有明显区别.

图10 与图9 相同,但驱动流场采用kmax/kc=16Fig.10 Same as Fig.9,but the flow field is filtered by kmax/kc=16

4 结论

在本文中,我们使用欧拉−拉格朗日方法对重粒子在可压缩各向同性湍流中的运动行为开展了先验研究.重点关注了不同滤波宽度和初始速度对粒子运动的影响.为了研究粒子聚集性的演化过程,本文引入了香农熵方法.根据本文所展示的结果,得出以下主要结论: 首先,小尺度流动结构对低S t数(S t＜5)的粒子聚集性有促进作用,而对高S t(St≥5),的粒子聚集性有抑制作用.粒子的速度和加速度的PDF随着S t的增加和截断波数kc的减小而变窄.因此,为了准确刻画重粒子的运动,需要考虑小尺度流动结构的影响.特别是对于低 Stokes 数的粒子,需要考虑小尺度的流动结构对聚集性的促进作用.其次,粒子的初始速度的不同对统计结果没有明显影响.即使在不同初始速度条件下,粒子最终的聚集性在统计定常状态下保持一致.

综上所述,本文的研究为重粒子在可压缩均匀各向同性湍流中的运动行为提供了先验分析,突出了小尺度流动结构和初始速度对粒子运动的重要性,并对粒子的聚集性演化进行了深入探讨.这些结果对于深入理解重粒子在湍流中的运动行为和进行相应的工程应用具有重要的参考价值.