賓遠(yuǎn)為 武 琦 夏振華 史一蓬,3)

* (北京大學(xué)湍流與復(fù)雜系統(tǒng)國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京 100871)

? (浙江大學(xué)流體工程研究所,杭州 310027)

引言

攜帶粒子的湍流問(wèn)題常見(jiàn)于自然界和工業(yè)工程應(yīng)用中,例如云層中雨滴的形成[1-2],大氣中污染物的擴(kuò)散[3-4],發(fā)動(dòng)機(jī)中的燃燒、霧化和混合等[5-7].由于湍流是一種多尺度的運(yùn)動(dòng)現(xiàn)象,粒子在湍流中的運(yùn)動(dòng)受到不同尺度流動(dòng)的影響.然而,盡管已經(jīng)進(jìn)行了大量理論、實(shí)驗(yàn)和數(shù)值模擬工作,但對(duì)其中物理機(jī)制的認(rèn)識(shí)仍然不夠透徹[8-11].

攜帶粒子的湍流數(shù)值模擬方法可以分為全解析(fully resolved) 模擬和點(diǎn)顆粒 (point-particle) 模擬[8-9].這兩類方法的分類依據(jù)是粒子尺度相對(duì)于流場(chǎng)特征尺度(例如Kolmogorov 尺度或大渦模擬中的濾波尺度) 的大小.全解析模擬要求流場(chǎng)的求解特征尺度(即網(wǎng)格尺寸)小于粒子尺度.在全解析模擬中,粒子的尺度是不可忽略的,因此它們被視為運(yùn)動(dòng)的壁面邊界條件,從而可以直接模擬粒子周圍所有尺度的流動(dòng)[10].湍流中粒子的運(yùn)動(dòng)由周圍流體施加的作用力和其他外力(如重力)共同決定.同時(shí),流體的運(yùn)動(dòng)也受到粒子的反作用力影響.因此,流體與粒子之間的相互作用是一個(gè)雙向耦合的過(guò)程.在某些情況下,當(dāng)粒子之間距離特別小時(shí),也需要考慮粒子之間的相互作用[12].然而,全解析模擬需要用小網(wǎng)格解析粒子及其周圍所有尺度內(nèi)的流動(dòng),因此計(jì)算成本較高,且對(duì)網(wǎng)格質(zhì)量要求較高.目前全解析模擬只能實(shí)現(xiàn)對(duì)單個(gè)粒子或數(shù)千個(gè)粒子的模擬[13-14].而對(duì)于實(shí)際情況中更常見(jiàn)的百萬(wàn)、千萬(wàn)量級(jí)的粒子運(yùn)動(dòng)求解,全解析模擬仍然面臨較大困難.

點(diǎn)粒子模擬是一種基于拉格朗日觀點(diǎn)的方法,其考慮的粒子尺度通常遠(yuǎn)小于流場(chǎng)特征尺度.因此,可以合理地將粒子視為流動(dòng)中的質(zhì)點(diǎn),并利用點(diǎn)粒子運(yùn)動(dòng)模型追蹤每個(gè)點(diǎn)的加速度、速度和位置[11,15].在點(diǎn)粒子模擬中,流體相的求解采用歐拉觀點(diǎn)的場(chǎng)方法,因此該方法也被稱為歐拉?拉格朗日法[15].Tsai[16]和Maxey 等[17]理論上給出了點(diǎn)粒子的受力模型,隨后被廣泛應(yīng)用于模擬湍流中點(diǎn)粒子的運(yùn)動(dòng)[18-20].在歐拉?拉格朗日法中,粒子與流場(chǎng)的耦合方式是關(guān)鍵因素.根據(jù)離散相的體積分?jǐn)?shù)(volume fraction)Φv和質(zhì)量分?jǐn)?shù)(mass fraction) Φm的不同,粒子與流場(chǎng)之間的耦合方式也不同[21].當(dāng)在氣固系統(tǒng)[22]中Φv≤10?6或液固系統(tǒng)[23]中 Φv≤10?4時(shí),離散相對(duì)流體相的作用和離散相之間的相互作用可以被合理忽略.數(shù)值模擬只需考慮離散相粒子在流場(chǎng)作用力下的運(yùn)動(dòng).這被稱為單向耦合(one-way coupling)[8,24].然而,當(dāng) Φv和 Φm相對(duì)較大時(shí),需要同時(shí)考慮流體相對(duì)離散相的作用力和離散相對(duì)流體相的反作用力,即雙向耦合(two-way coupling).隨著 Φv和 Φm的進(jìn)一步增大,還需要考慮離散相之間的相互作用,即四向耦合(four-way coupling).相比于全解析模擬,點(diǎn)粒子模擬對(duì)計(jì)算量的要求更小,因此在大規(guī)模粒子運(yùn)動(dòng)求解方面具有優(yōu)勢(shì).目前,點(diǎn)粒子模擬方法已廣泛應(yīng)用于湍流直接數(shù)值模擬和大渦模擬中[8-9,25].

直接數(shù)值模擬(direct numerical simulation,DNS)是一種精確的湍流模擬方法.歐拉?拉格朗日點(diǎn)粒子模型已被成功應(yīng)用于一些顆粒湍流相互作用問(wèn)題的DNS 研究中.例如,Maxey[18]使用單向耦合的點(diǎn)粒子模擬計(jì)算粒子在均勻湍流中的沉降.Squires等[26]研究了不可壓縮各向均勻同性湍流中粒子的慣性聚集,以及粒子對(duì)湍流的反作用[27].Bec 等[28]著重研究了重粒子在各向均勻同性湍流中的加速度概率分布并探究了不同Stokes 數(shù)(S t)對(duì)加速度的影響.Xia 等[29]進(jìn)一步在可壓縮均勻湍流中研究了重粒子對(duì)湍流場(chǎng)的調(diào)制作用.然而,DNS 本身計(jì)算成本高昂且計(jì)算需求苛刻,目前的歐拉?拉格朗日點(diǎn)粒子的DNS 仍受限于低雷諾數(shù)和簡(jiǎn)單流動(dòng)場(chǎng)景[9,30].

在湍流問(wèn)題的數(shù)值研究中,大渦模擬(largeeddy simulation,LES)方法近年來(lái)得到快速發(fā)展,逐漸在學(xué)術(shù)界和工程應(yīng)用中占據(jù)重要地位.LES 通過(guò)直接求解湍流中的大尺度運(yùn)動(dòng),而小尺度運(yùn)動(dòng)對(duì)大尺度運(yùn)動(dòng)的影響則通過(guò)亞格子(sub-grid stress,SGS)模型來(lái)表示.由于不需要解析小尺度運(yùn)動(dòng),LES 的計(jì)算量遠(yuǎn)小于DNS[31].盡管單相湍流的SGS 模型近些年來(lái)得到不斷完善,但是否可以直接應(yīng)用于含顆粒的兩相湍流問(wèn)題仍值得進(jìn)一步研究.在LES 中,由于只能模擬大尺度運(yùn)動(dòng)對(duì)應(yīng)的流場(chǎng),原先存在于DNS流場(chǎng)中的小尺度結(jié)構(gòu)無(wú)法被解析,湍流場(chǎng)中的小尺度運(yùn)動(dòng)與粒子的相互作用也被忽略.然而,這種忽略小尺度運(yùn)動(dòng)的合理性尚未完全驗(yàn)證,因此需要進(jìn)一步探索.

為了研究這一問(wèn)題,一個(gè)直接方法是考慮粒子在濾波后DNS 流場(chǎng)(filtered DNS,fDNS)中的運(yùn)動(dòng).在不可壓縮均勻各向同性湍流(incompressible homogenous isotropic turbulence,IHIT)中,已經(jīng)進(jìn)行了相關(guān)研究.例如,Fede 等[32]研究了不同濾波寬度下濾波后的大尺度流動(dòng)中粒子對(duì)的運(yùn)動(dòng),以研究粒子碰撞和聚集性.研究結(jié)果表明,只有當(dāng)濾波寬度足夠大,以至于流體動(dòng)能明顯下降時(shí),粒子的聚集性才會(huì)受到影響.Jin 等[33]進(jìn)一步研究了LES中SGS 模型對(duì)于粒子在湍流中統(tǒng)計(jì)特性的影響.他們的計(jì)算分析表明,小尺度流動(dòng)對(duì)低密度粒子碰撞、聚集的影響是不可忽視的.因此,他們提出對(duì)低密度粒子的LES 必須要考慮SGS 模型,否則可能導(dǎo)致高達(dá)60%的相對(duì)誤差.Ray 等[34]選取了多個(gè)濾波寬度,大致涵蓋了整個(gè)慣性子區(qū)范圍,主要研究粒子徑向分布函數(shù)(RDF)和粒子徑向相對(duì)速度隨濾波寬度的變化.他們發(fā)現(xiàn)低密度粒子的RDF 會(huì)隨著濾波寬度變大而降低,而高密度粒子的趨勢(shì)相反.這一結(jié)果驗(yàn)證了Jin 等[33]的結(jié)論,即不能忽略小尺度流動(dòng)對(duì)粒子在湍流中運(yùn)動(dòng)的影響.Pozorski 等[35]在綜述中對(duì)LES 濾波對(duì)粒子聚集的影響進(jìn)行了回顧,并總結(jié)了不同學(xué)者提出的可提升粒子計(jì)算準(zhǔn)確性的SGS 模型.在對(duì)粒子運(yùn)動(dòng)統(tǒng)計(jì)量的研究方面,Lalescu 等[36]研究了示蹤粒子在濾波流場(chǎng)中與加速度有關(guān)的統(tǒng)計(jì)量.示蹤粒子可以視為密度趨近于零的慣性粒子.他們發(fā)現(xiàn)濾波后的流場(chǎng)會(huì)使粒子加速度的間歇性減弱.這表現(xiàn)為,在濾波后流場(chǎng)中粒子加速度較未濾波的流場(chǎng)中的加速度變化更為平緩.

盡管目前已有許多關(guān)于重粒子在可壓縮湍流中的物理機(jī)理研究[29,37],但關(guān)于LES 中濾波寬度的影響還沒(méi)有得到系統(tǒng)研究.因此,本文以含重粒子的可壓縮均勻各向同性湍流(compressible homogenous isotropic turbulence,CHIT)為研究對(duì)象,通過(guò)先驗(yàn)分析研究不同尺度的湍流結(jié)構(gòu)和粒子的初速度對(duì)粒子的聚集行為以及后續(xù)穩(wěn)態(tài)統(tǒng)計(jì)結(jié)果的影響.

1 數(shù)值細(xì)節(jié)

1.1 流體相控制方程

三維的可壓縮Navier-Stokes 方程為

其中,ui,ρ,p和T分別是無(wú)量綱速度分量、質(zhì)量、壓強(qiáng)和溫度.Fi是施加在流體上的大尺度力,Λ 是大尺度的冷卻函數(shù).黏性應(yīng)力 σij和總能E被分別定義為

其中 θ=?uk/?xk是無(wú)量綱的速度散度.分子黏性 μ 和熱導(dǎo)率 κ 基于常數(shù)普朗特?cái)?shù)假設(shè)[38],由Sutherland's律[39]給定

在上述無(wú)量綱方程中,我們使用了特征物理量L,U,ρ0,ρ0RT0,ρ0U2和c0≡對(duì)長(zhǎng)度、速度、密度、壓強(qiáng)、能量和聲速進(jìn)行無(wú)量綱化.分子黏性μ和熱導(dǎo)率 κ 由 μ0和 κ0無(wú)量綱化.這里,比熱γ ≡Cp/Cv為1.4.Cp,Cv和R≡Cp?Cv分別是定壓比熱容.定容比熱容和氣體常數(shù).流動(dòng)由3 個(gè)相對(duì)參數(shù)控制,它們分別是相對(duì)雷諾數(shù)Re=ρ0U/L,相對(duì)馬赫數(shù)M=U/c0和相對(duì)普朗特?cái)?shù)Pr=μ0Cp/κ0.在本研究中,這3 個(gè)參數(shù)分別為Pr=0.7,Re=200 和M=0.45.參數(shù) α 的定義為 α=PrRe(γ?1)M2.

1.2 粒子控制方程

當(dāng)粒子直徑d遠(yuǎn)小于Kolmogorov 長(zhǎng)度尺度 η,且粒子密度 ρp遠(yuǎn)大于流體密度 ρ 時(shí),粒子位置xp,i和速度vp,i由下式?jīng)Q定[16-17,29]

粒子控制方程(式(9)和式(10))是基于牛頓不可壓縮流體得到的.能否將這一方程用于可壓縮流動(dòng),需要特別討論.Loth[42]指出當(dāng)由粒子與當(dāng)?shù)亓黧w的相對(duì)速度和當(dāng)?shù)亓黧w聲速定義的粒子相對(duì)馬赫數(shù)時(shí),粒子的控制方程必須要考慮可壓縮性的影響.在本文中,即使對(duì)于密度最大的粒子(S t=10),平均的粒子相對(duì)馬赫數(shù)也僅約等于 0.2.所以本文是可以合理地忽略可壓縮效應(yīng)對(duì)粒子控制方程的影響.

1.3 數(shù)值方法

在本研究的數(shù)值模擬中,計(jì)算域?yàn)殚L(zhǎng)寬高均為2π的立方盒子,3 個(gè)方向均為周期性邊界條件.模擬采用了 2563的均勻網(wǎng)格.其中空間離散格式采用了混合的緊致加權(quán)本質(zhì)無(wú)振蕩 (WENO) 格式[38].該混合格式結(jié)合了光滑區(qū)域的八階緊致有限差分格式[43]和激波區(qū)域的7 階WENO 格式[44].時(shí)間推進(jìn)選取3 階龍格?庫(kù)塔格式.計(jì)算時(shí)間步長(zhǎng)由Δt=CFLΔx/max(uiui)動(dòng)態(tài)決定.在我們的研究中,網(wǎng)格大小 Δx等于2π/256 ≈0.024 5,CFL數(shù)設(shè)置為 0.05.流動(dòng)由在大尺度上施加的剪切力 Fi驅(qū)動(dòng).我們?cè)谟?jì)算中假設(shè)粒子的直徑為d=0.1Δx.在本工作中,Kolmogorov 尺度η ≈Δx,粒子的直徑d≈0.1η,這使得粒子直徑d遠(yuǎn)小于Kolmogorov 長(zhǎng)度尺度 η 的假設(shè)成立.粒子的密度由Stokes 數(shù)的定義給定.另外百萬(wàn)顆粒對(duì)應(yīng)的體積分?jǐn)?shù) Φv約為 3 ×10?5.對(duì)于S t=10 的算例,顆粒的質(zhì)量分?jǐn)?shù)約為0.5.從嚴(yán)格意義上講,本文所涉及的參數(shù)已經(jīng)略微超出單向耦合成立的范圍,顆?？赡軙?huì)對(duì)流場(chǎng)有一定的反饋?zhàn)饔肹24].但由于本文聚焦于流場(chǎng)對(duì)顆粒運(yùn)動(dòng)的影響,因此選擇忽略該反饋?zhàn)饔肹40].后續(xù)將進(jìn)一步研究該反饋?zhàn)饔玫挠绊?粒子位置和速度的時(shí)間推使用二階龍格?庫(kù)塔格式.粒子的時(shí)間步長(zhǎng) Δtp與流動(dòng)的時(shí)間步長(zhǎng) Δt相同.第i個(gè)粒子位置的流場(chǎng)物理屬性值在光滑區(qū)域或激波區(qū)域分別由線性插值或非線性插值得到.

在當(dāng)前新農(nóng)村發(fā)展過(guò)程當(dāng)中，體育相關(guān)產(chǎn)業(yè)的發(fā)展也發(fā)揮著非常重要的影響。以往環(huán)境下，受限于農(nóng)村地區(qū)經(jīng)濟(jì)發(fā)展水平，農(nóng)民群眾的消費(fèi)意識(shí)相對(duì)有限，不具備良好的消費(fèi)習(xí)慣和消費(fèi)意識(shí)。內(nèi)需不足對(duì)產(chǎn)業(yè)發(fā)展是十分不利的。

當(dāng)流動(dòng)達(dá)到統(tǒng)計(jì)穩(wěn)態(tài)后,粒子St數(shù)為1 時(shí),在圖1 中繪制了某一時(shí)刻流場(chǎng)的無(wú)量綱散度 θ/θrms,無(wú)量綱渦量大小 |ω|/|ω|rms的等值圖和粒子位置.由于本文中的瞬時(shí)湍流馬赫數(shù)Mt=urms/c0在0.987～0.989 之間,可以在圖1(a) 觀察到大量的小激波(Shocklets,θ/θrms＜?3).這些激波對(duì)粒子的運(yùn)動(dòng)會(huì)起到明顯影響[40].圖1(b) 顯示了無(wú)量綱渦量大小|ω|/|ω|rms和粒子的對(duì)應(yīng)位置.圖中可以看出粒子局部聚集于低渦量的區(qū)域內(nèi).

為了驗(yàn)證時(shí)間步長(zhǎng)度的收斂性,我們將 CFL 數(shù)分別調(diào)整為0.05,0.025 和0.012 5.因此,在每種對(duì)應(yīng)的情況下,時(shí)間步長(zhǎng)變?yōu)?Δt,Δt/2 和 Δt/4.驗(yàn)證算例包含了不同的S t數(shù),即St=1,5,10,并使用I.C.4作為初始條件.之后,隨機(jī)選擇一個(gè)粒子并跟蹤它的路徑線,如圖2 所示.很明顯,對(duì)于相同的S t數(shù),不同時(shí)間步長(zhǎng)的路徑線是相同的.這意味著當(dāng)CFL=0.05時(shí),時(shí)間步長(zhǎng)已經(jīng)收斂.

圖2 一個(gè)隨機(jī)選擇的粒子在流動(dòng)中的路徑.初始速度采用I.C.4Fig.2 Path of a randomly selected particle.Here,the initial velocity is I.C.4

1.4 初始條件

本節(jié)介紹模擬中的初始條件的生成過(guò)程.在無(wú)粒子情況下,首先通過(guò)DNS 模擬得到了達(dá)到統(tǒng)計(jì)定常狀態(tài)的可壓縮均勻各向同性湍流.為了確保湍流達(dá)到統(tǒng)計(jì)定常狀態(tài),流場(chǎng)演化了5 個(gè)大渦翻轉(zhuǎn)時(shí)間(TE)以保證湍流達(dá)到統(tǒng)計(jì)定常狀態(tài).然后,該流場(chǎng)被保存,并作為模擬粒子運(yùn)動(dòng)的初始流場(chǎng).總共Np=1.0×106個(gè)粒子以不同的給定初始速度隨機(jī)撒落在計(jì)算域中.為了避免粒子初始位置對(duì)結(jié)果的影響,這些粒子的位置也被保存用于后續(xù)的模擬.因此,在本研究中的所有算例中,都采用相同的初始流場(chǎng)和粒子位置.本文中粒子有5 種給定的初速度,具體見(jiàn)表1.具體來(lái)說(shuō),對(duì)于I.C.1,粒子被靜態(tài)地放置在流場(chǎng)中.對(duì)于I.C.2～I(xiàn).C.4,分別給每個(gè)粒子設(shè)置對(duì)應(yīng)的給定速度作為初始速度.在I.C.5 中,粒子的初始速度設(shè)定為當(dāng)?shù)氐牧鲃?dòng)速度.

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

1.5 粒子聚集行為表征方法

本文使用了兩種方法對(duì)粒子聚集行為進(jìn)行表征,分別為概率密度函數(shù)(PDF)方法和香農(nóng)熵(Shannon entropy)方法.其中,PDF 方法適用于對(duì)粒子相關(guān)物理量進(jìn)行觀察,尤其在統(tǒng)計(jì)穩(wěn)態(tài)時(shí)較為有效,但無(wú)法追蹤粒子聚集程度隨時(shí)間的演化.而香農(nóng)熵方法則能夠追蹤粒子聚集時(shí)間的演化.香農(nóng)熵方法將整個(gè)計(jì)算域統(tǒng)一劃分為Nbin=1283個(gè)小盒子.在本文中,計(jì)算域?yàn)檫呴L(zhǎng)為 2π 的正立方體,所以,經(jīng)過(guò)劃分,會(huì)在計(jì)算域中得到 1283個(gè)大小相同的邊長(zhǎng)為2π/128的正立方體.香農(nóng)熵S(t)由以下方程定義

其中 max[H(t)]=lnNbin,p(k,t) 是在時(shí)間t時(shí),在第k個(gè)小盒子中發(fā)現(xiàn)一個(gè)粒子的頻率,Np(k,t) 是在時(shí)間t時(shí),第k個(gè)小盒子中的粒子數(shù).當(dāng)粒子足夠均勻分布時(shí),S(t)趨近于1.而當(dāng)所有粒子都集中在一個(gè)小盒子中時(shí),S(t)等于0.在本文的算例中,在初始時(shí)刻S(0)約等于0.93.根據(jù)該定義,可以使用任一瞬時(shí)流場(chǎng)來(lái)計(jì)算香農(nóng)熵S(t).

2 不同濾波寬度的影響

在本節(jié)中,我們研究不同濾波寬度的影響.由于湍流中存在非線性效應(yīng)和黏性效應(yīng),湍流中渦的尺度從大尺度L到Kolmogorov 長(zhǎng)度尺度 η 是連續(xù)分布的.同時(shí),流體動(dòng)能從最大的渦傳輸?shù)阶钚〉臏u.這個(gè)過(guò)程被稱為能量級(jí)串.在 CHIT 中,能量譜E(k)與波數(shù)k的? 5/3 冪律已被大量文獻(xiàn)[38,45]證明.

圖3 中實(shí)線顯示了達(dá)到統(tǒng)計(jì)定常狀態(tài)的能譜E(k).由于網(wǎng)格分辨率的限制,我們的DNS 中最大可識(shí)別波數(shù)是kmax=128.為了對(duì)流場(chǎng)進(jìn)行空間濾波,本文采用了經(jīng)典的譜截?cái)酁V波器,其數(shù)學(xué)形式為

圖3 流動(dòng)能量譜 E(k) (實(shí)線)和5 個(gè)不同的截?cái)嗖〝?shù)(虛線).從左到右分別是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

當(dāng)CHIT 達(dá)到統(tǒng)計(jì)定常狀態(tài)時(shí),粒子被均勻隨機(jī)地撒落到計(jì)算域中.在湍流中運(yùn)動(dòng)的粒子最終會(huì)達(dá)到統(tǒng)計(jì)上的穩(wěn)定狀態(tài),并會(huì)顯示出局部聚集的空間分布,這被稱為慣性聚集現(xiàn)象[26].在上述過(guò)程中,存在兩個(gè)目前尚不清楚的問(wèn)題.首先,從初始時(shí)刻到粒子達(dá)到統(tǒng)計(jì)穩(wěn)定狀態(tài)的時(shí)刻,不同尺度的渦對(duì)粒子聚集的影響尚不明確.其次,當(dāng)達(dá)到統(tǒng)計(jì)穩(wěn)定狀態(tài)時(shí),粒子聚集的程度受到kc和S t的影響程度也不明確.

圖4 顯示了粒子以I.C.2 為初始速度條件時(shí),隨時(shí)間演化的香農(nóng)熵S(t) 在不同的S t和kc情況下的變化.通過(guò)圖4,可以得出以下4 點(diǎn)結(jié)論: (1) 圖4 中的S(t) 首先減小并趨于常數(shù).這表明粒子撒落在計(jì)算域后會(huì)逐漸聚集,直到達(dá)到統(tǒng)計(jì)定常狀態(tài).(2) 通過(guò)比較粒子直接由DNS 流場(chǎng)驅(qū)動(dòng)的結(jié)果,我們觀察到在統(tǒng)計(jì)定常狀態(tài)時(shí),低密度粒子(S t=1)的聚集程度(S(∞)≈0.88)比高密度粒子(St=5,10)的聚集程度(S(∞)≈0.91,0.92) 更強(qiáng).(3) 當(dāng)kc從kmax/kc=2 變化到kmax/kc=8 時(shí),S(t) 并沒(méi)有表現(xiàn)出明顯的差異,這表明,當(dāng)kc選取在合適的范圍內(nèi),使用LES 是能夠得到可信的粒子聚集性統(tǒng)計(jì)量的.(4) 如果kmax/kc大于8,隨著kmax/kc增加,S(t) 下降的速率更快.此外,當(dāng)kmax/kc增加,在統(tǒng)計(jì)定常狀態(tài)時(shí),低密度粒子的聚集性降低,而高密度粒子的聚集性增加.因此,小尺度渦有利于低密度粒子的聚集.相反,當(dāng)S t較大時(shí),例如St＞5,小尺度渦對(duì)高密度粒子聚集性的影響不大.該結(jié)論與文獻(xiàn)[32]中的研究結(jié)果一致.

圖4 S(t) 隨著時(shí)間的演化.初始速度采用I.C.2Fig.4 S(t) as a function of time.Initial velocity is I.C.2

從粒子空間數(shù)密度的概率密度函數(shù)(PDF)可以得到相同的結(jié)論.數(shù)密度的獲取首先可以將計(jì)算域劃分為若干個(gè)立方體,在粒子達(dá)到統(tǒng)計(jì)定常狀態(tài)后,可以得到每個(gè)立方體中的平均粒子數(shù)n,稱為平均數(shù)密度[26].平均數(shù)密度n為零的概率越大,代表粒子的聚集性越強(qiáng).本文將整個(gè)計(jì)算域劃分為 2563的立方體,并統(tǒng)計(jì)達(dá)到穩(wěn)態(tài)后關(guān)于n的PDF.圖5 繪制了在圖4 中所有算例的PDF,并證明可以得到與之前相同的結(jié)論.具體地,圖5 顯示出達(dá)到穩(wěn)態(tài)后,在低密度粒子情況下,增加kmax/kc會(huì)減弱粒子的聚集性.而在高密度粒子情況下,增加kmax/kc會(huì)增加粒子的聚集性.

圖5 平均數(shù)密度 n 的PDF.初始速度采用I.C.2Fig.5 PDF of the average particle number density n.Initial velocity is I.C.2

圖6 與圖5 相同,但是為粒子絕對(duì)加速度的PDFFig.6 Same as Fig.5,but for the magnitudes of the particles’acceleration

圖7 與圖5 相同,但是為粒子絕對(duì)速度的PDFFig.7 Same as Fig.5,but for the magnitudes of the particles’ velocity

3 初始速度影響

在不同的流動(dòng)場(chǎng)景模擬中,流體中粒子的初速度由不同的情況給出.例如,初始時(shí)刻靜止的顆粒沉降[46](I.C.1),含顆粒的射流對(duì)應(yīng)[47](I.C.2～I(xiàn).C.4),和示蹤粒子(I.C.5).關(guān)于這些初始條件的詳細(xì)信息已經(jīng)在表1 中列出.雖然許多研究已經(jīng)使用了這些初始條件,但在作者了解的范圍里,初始條件的影響還沒(méi)有得到充分的研究.在本節(jié)中,將使用香農(nóng)熵和平均速度來(lái)分析這些初始條件對(duì)粒子聚集行為的影響.為了驗(yàn)證LES 的適用性,本文選擇了一個(gè)合適的截?cái)嗖〝?shù),即kmax/kc=4,該值在 CHIT 的LES 模擬中被廣泛使用.

粒子平均速度的變化表明了粒子的運(yùn)動(dòng)學(xué)過(guò)程.如圖8 所示,對(duì)于各種初始條件,平均速度在經(jīng)過(guò)初始的弛豫時(shí)間后會(huì)趨于相同的值.且該值會(huì)隨著S t的增加而增加.此外,與傳統(tǒng)觀點(diǎn)相反,我們發(fā)現(xiàn)在粒子初始時(shí)刻等于當(dāng)?shù)亓鲃?dòng)速度的情況下(I.C.5),其平均速度需要較長(zhǎng)的時(shí)間才能收斂到統(tǒng)計(jì)定常值.這里的統(tǒng)計(jì)定常值為不同初速度的算例達(dá)到同一平均速度的對(duì)應(yīng)速度值.在圖8 中可以看I.C.5 對(duì)應(yīng)的平均速度是最晚與其他算例一致的.值得一提的是,速度達(dá)到統(tǒng)計(jì)定常值并不意味著粒子運(yùn)動(dòng)已經(jīng)達(dá)到統(tǒng)計(jì)定常狀態(tài).例如,圖8 (a) 所示,當(dāng)t/TE=0.5,S t=1時(shí),對(duì)于所有初始條件,其平均速度值已經(jīng)收斂,但香農(nóng)熵仍然還未收斂.同樣的現(xiàn)象也可以在S t=1 和 5 時(shí)觀察到.

圖8 粒子數(shù)平均速度隨時(shí)間的演化.驅(qū)動(dòng)流場(chǎng)采用kmax/kc=4Fig.8 Averaged particle velocity as a function of time.Flow field is filtered by kmax/kc=4

從圖9 所示,可以觀察到不同初始速度條件導(dǎo)致聚集性的不同演化.對(duì)于所有不同的S t,粒子初速度為I.C.1～I(xiàn).C.3 時(shí)得到的S(t) 的演化基本相同.然而,當(dāng)初速度過(guò)大時(shí)(I.C.4),S(t) 的演化與初速度較小的情況會(huì)有一些差異.此外,當(dāng)粒子初始速度為當(dāng)?shù)亓鲃?dòng)速度時(shí)(I.C.5),粒子較快達(dá)到統(tǒng)計(jì)穩(wěn)定狀態(tài).最后,可以得出結(jié)論,不同的初始條件不影響統(tǒng)計(jì)定常狀態(tài)下粒子的聚集程度.

圖9 S(t) 隨時(shí)間的演化.驅(qū)動(dòng)流場(chǎng)采用kmax/kc=4Fig.9 S(t) as a function of time.Flow field is filtered bykmax/kc=4

通過(guò)前面的討論,本文可以得出以下結(jié)論,當(dāng)kmax/kc=16且粒子初速度采用 I.C.1 時(shí)的情況下,S(t)的演化沒(méi)有明顯差異.但是,對(duì)于其他初始條件,該結(jié)論是否適用仍然是一個(gè)問(wèn)題.圖10 展示了當(dāng)kmax/kc=16 時(shí) 對(duì)應(yīng)I.C.1～I(xiàn).C.5 的S(t) 的演化.與圖9相比,對(duì)于我們給定的所有初始條件,粒子的聚集性演化過(guò)程沒(méi)有明顯區(qū)別.

圖10 與圖9 相同,但驅(qū)動(dòng)流場(chǎng)采用kmax/kc=16Fig.10 Same as Fig.9,but the flow field is filtered by kmax/kc=16

4 結(jié)論

在本文中,我們使用歐拉?拉格朗日方法對(duì)重粒子在可壓縮各向同性湍流中的運(yùn)動(dòng)行為開(kāi)展了先驗(yàn)研究.重點(diǎn)關(guān)注了不同濾波寬度和初始速度對(duì)粒子運(yùn)動(dòng)的影響.為了研究粒子聚集性的演化過(guò)程,本文引入了香農(nóng)熵方法.根據(jù)本文所展示的結(jié)果,得出以下主要結(jié)論: 首先,小尺度流動(dòng)結(jié)構(gòu)對(duì)低S t數(shù)(S t＜5)的粒子聚集性有促進(jìn)作用,而對(duì)高S t(St≥5),的粒子聚集性有抑制作用.粒子的速度和加速度的PDF隨著S t的增加和截?cái)嗖〝?shù)kc的減小而變窄.因此,為了準(zhǔn)確刻畫重粒子的運(yùn)動(dòng),需要考慮小尺度流動(dòng)結(jié)構(gòu)的影響.特別是對(duì)于低 Stokes 數(shù)的粒子,需要考慮小尺度的流動(dòng)結(jié)構(gòu)對(duì)聚集性的促進(jìn)作用.其次,粒子的初始速度的不同對(duì)統(tǒng)計(jì)結(jié)果沒(méi)有明顯影響.即使在不同初始速度條件下,粒子最終的聚集性在統(tǒng)計(jì)定常狀態(tài)下保持一致.

綜上所述,本文的研究為重粒子在可壓縮均勻各向同性湍流中的運(yùn)動(dòng)行為提供了先驗(yàn)分析,突出了小尺度流動(dòng)結(jié)構(gòu)和初始速度對(duì)粒子運(yùn)動(dòng)的重要性,并對(duì)粒子的聚集性演化進(jìn)行了深入探討.這些結(jié)果對(duì)于深入理解重粒子在湍流中的運(yùn)動(dòng)行為和進(jìn)行相應(yīng)的工程應(yīng)用具有重要的參考價(jià)值.