雷穎昊南，雷三惠，張生浩，王 萍

(蘭州大學(xué) 顆粒湍流研究中心，蘭州 710000)

0 引 言

可壓縮顆粒兩相流常見(jiàn)于超燃發(fā)動(dòng)機(jī)[1-2]、沙塵或雨中飛行器航行[3-4]和混合物沖擊/爆炸[5]等。另一個(gè)典型的例子是高超聲速飛行器再入大氣過(guò)程中，熱防護(hù)材料高溫?zé)g“層裂”過(guò)程中產(chǎn)生顆粒[6-7]，與飛行器周圍氣流相互作用。因?yàn)閱蝹€(gè)顆?？赡茉谙”』蜻B續(xù)氣體中經(jīng)歷亞聲速到高超聲速的流態(tài)，這些類型氣流中的顆粒運(yùn)動(dòng)建模很復(fù)雜。同時(shí)，由于兩相系統(tǒng)的復(fù)雜性和多尺度性，用理論和實(shí)驗(yàn)方法進(jìn)行研究存在諸多困難，數(shù)值模擬可直觀地展示顆粒內(nèi)部的物理現(xiàn)象，因此可發(fā)揮較大作用。

用于模擬顆粒兩相流的方法大體上可以分為歐拉-歐拉、歐拉-拉格朗日兩類。在歐拉-歐拉方法中，流體和顆粒均被視為連續(xù)介質(zhì)，采用Navier-Stokes方程組求解[8]，該方法適用于顆粒體積分?jǐn)?shù)較高的情況；而在歐拉-拉格朗日方法中，氣體方程在歐拉坐標(biāo)系下離散，顆粒在拉格朗日坐標(biāo)系下跟蹤，所受到的力可能包括拖曳力、升力、布朗力、熱泳力等[9]。假設(shè)氣相占據(jù)整個(gè)空間，而顆粒被考慮為具有有限質(zhì)量和動(dòng)量的點(diǎn)源。當(dāng)顆粒體積分?jǐn)?shù) φp較低時(shí)，顆粒相與流體間的動(dòng)量交換可以忽略，稱之為單向耦合；增大φp至顆粒與流體間相互作用不可忽略時(shí)，兩相通過(guò)質(zhì)量（以化學(xué)反應(yīng)的形式）、動(dòng)量和能量方程雙向耦合；更進(jìn)一步增大 φp，在密相流的歐拉-拉格朗日方法中需要考慮顆粒間碰撞，即四向耦合[10]。

歐拉-拉格朗日方法由于可以更為準(zhǔn)確地模擬單個(gè)顆粒行為，最近被廣泛用于可壓縮多相流研究。如Zhang、Dai、Xiao 等采用湍流直接數(shù)值模擬求解器[11-13]，結(jié)合點(diǎn)顆粒運(yùn)動(dòng)模型研究了顆粒的傾向性聚集和擴(kuò)散。Dai、Xia 等研究了湍流受顆粒調(diào)制的規(guī)律[14-15]。Li 等[16]在FLUENT 軟件求解流場(chǎng)的基礎(chǔ)上，開(kāi)發(fā)顆粒求解器代碼，研究超重力作用下亞微米顆粒在超聲速層流邊界層中的運(yùn)動(dòng)。隨后，其在此基礎(chǔ)上進(jìn)一步研究了超聲速繞楔流中的顆粒沉積[9]。Teh 等[17]基于OpenFOAM 模擬流場(chǎng)研究了顆粒對(duì)斜激波與層流邊界層相互作用的影響。Palmer 等[4]在美國(guó)宇航局艾姆斯研究中心開(kāi)發(fā)的計(jì)算流體動(dòng)力學(xué)解算器（DPLR，其中包含有限速率反應(yīng)動(dòng)力學(xué)、熱和化學(xué)非平衡、精確的高溫傳輸系數(shù)和電離流物理的通用模型）的基礎(chǔ)上，根據(jù)火星的地面和大氣觀測(cè)，使用一組耦合的常微分方程計(jì)算沙塵通過(guò)激波層的軌跡，研究沙塵顆粒影響下火星再入飛行器隔熱層侵蝕。Davuluri 等[7]采用歐拉-拉格朗日方法數(shù)值研究燒蝕層裂現(xiàn)象中表面噴射后剝落顆粒的軌跡，該研究使用高超聲速KATS-Kentucky 氣動(dòng)熱力學(xué)和熱響應(yīng)求解器，用于獲得流場(chǎng)。Davuluri 等[18]在此基礎(chǔ)上進(jìn)一步討論了顆?？刂品匠讨械耐弦妨?，利用基于雷諾數(shù)、馬赫數(shù)和克努森數(shù)的一系列經(jīng)驗(yàn)關(guān)系，建立了混合阻力系數(shù)模型。

本文基于國(guó)產(chǎn)自主研發(fā)的計(jì)算流體力學(xué)（Computational Fluid Dynamics，CFD）軟件PHengLEI（風(fēng)雷）[19]，開(kāi)發(fā)可壓縮邊界層顆粒多相流程序，為國(guó)家數(shù)值風(fēng)洞軟件開(kāi)發(fā)提供技術(shù)儲(chǔ)備和集成模塊?；诘湫退憷龑?duì)所開(kāi)發(fā)的模塊進(jìn)行檢驗(yàn)，并初步研究了層流邊界層中顆粒的運(yùn)動(dòng)行為及其對(duì)流動(dòng)性質(zhì)的影響。

1 數(shù)值模型和方法

風(fēng)雷軟件是中國(guó)空氣動(dòng)力研究與發(fā)展中心研發(fā)的面向流體工程的混合CFD 平臺(tái)。風(fēng)雷軟件的計(jì)算范圍覆蓋低速、亞跨聲速和高超聲速，是一款具有完全自主知識(shí)產(chǎn)權(quán)的面向工程應(yīng)用和學(xué)術(shù)研究的通用CFD 軟件，借鑒了面向?qū)ο蟮拇笮蛙浖O(shè)計(jì)理念，采用C++語(yǔ)言編程。風(fēng)雷軟件的整體框架可參考文獻(xiàn)[19-20]，這里僅介紹與多相流開(kāi)發(fā)相關(guān)的數(shù)值模型和方法。

1.1 流體控制方程

本文考慮了顆粒-流體的雙向耦合，在求解過(guò)程中，通過(guò)在流體相控制方程中添加顆粒的動(dòng)量、能量源項(xiàng)的方式，來(lái)體現(xiàn)顆粒對(duì)流體相的反作用。實(shí)際程序計(jì)算中，通過(guò)定義歐拉變量，將每一時(shí)間步各單元的反饋源項(xiàng)存儲(chǔ)到單元中心處。笛卡爾坐標(biāo)系下，包含顆粒源項(xiàng)的守恒形式流動(dòng)控制方程為：

其中，ui、 ρ、Tf、E分別表示流體的速度、密度、溫度、總能，k為流體傳熱系數(shù)， σij為黏性切應(yīng)力張量，R為氣體常數(shù)，計(jì)算中取值為 287.974 4 J/(kg·K)。流體計(jì)算采用二階有限體積方法，空間上采用Roe 格式，時(shí)間上采用LUSGS 格式。顆粒的動(dòng)量源項(xiàng) φui以及能量源項(xiàng) φei的表達(dá)式分別為：

其中，F(xiàn)p,i是顆粒受到的力，Tp為 顆粒的溫度，Qp,i是顆粒釋放的熱量，Gp,i是 顆粒對(duì)流體做的功，dp是顆粒的粒徑，Nu是努塞特?cái)?shù)，np表示計(jì)算網(wǎng)格單元體積內(nèi)的顆粒數(shù)， ΔV表示網(wǎng)格單元的體積，up,i是顆粒速度。

1.2 計(jì)算相關(guān)參數(shù)

本文算例的流體介質(zhì)為空氣，其可壓縮計(jì)算涉及的黏度由著名的Sutherland 公式給出：

其中，Ts為 薩瑟蘭常數(shù)，計(jì)算中一般取 1 10.4 K，T0、μ0為參考狀態(tài)的溫度與動(dòng)力學(xué)黏度，計(jì)算中分別取228.15 K、1.789 4×10?5kg/(m·s)。計(jì)算中涉及的無(wú)量綱參數(shù)：普朗特?cái)?shù)Pr=cF,pμk，cF,p為流體的比定壓熱容；單位雷諾數(shù)Re?∞=ρ∞U∞/μ∞， ρ∞和 μ∞分別為來(lái)流密度和黏度；來(lái)流馬赫數(shù)Ma∞=U∞/c∞，c∞為來(lái)流聲速。

1.3 顆粒控制方程

本文離散顆粒相計(jì)算的控制方程分別為顆粒位移、顆粒平動(dòng)、顆粒轉(zhuǎn)動(dòng)及顆粒溫度控制方程，見(jiàn)式（10～13）：

其中，Xp、vp、 ωp是顆粒的運(yùn)動(dòng)學(xué)參數(shù)，分別是位移、速度、角速度，Tp為 顆粒溫度，Ip是顆粒的慣性矩，F(xiàn)p和Tc分 別是顆粒受到的力及扭矩，cp,p為顆粒相物質(zhì)的比定壓熱容。

在顆粒受力方面，Magnaudet 等[21]和Ling 等[22-23]對(duì)不可壓縮和可壓縮流動(dòng)中顆粒受力及其發(fā)展進(jìn)行了全面概述。然而，Thomas[24]、Tedeschi[25-26]和Hughes[27]等的研究發(fā)現(xiàn)，對(duì)于粒徑在0.1～100.0 μm范圍內(nèi)的顆粒，非定常力并不顯著。因此，由于本文研究的顆粒-流體密度比例較大且顆粒尺度較低，非定常力對(duì)結(jié)果影響較小。這里考慮的顆粒受力模型參考了Li[9]的工作，包括拖曳力FD、 升力FSaff、布朗力FBi、熱泳力Ft等。顆粒受力Fp表達(dá)式見(jiàn)式（14）：

本文研究的問(wèn)題中，流動(dòng)范圍涉及亞聲速到高超聲速、連續(xù)流到過(guò)渡流區(qū)，顆粒從低雷諾數(shù)到高雷諾數(shù)，因此，所采用的拖曳力FD模型需要適用較大的顆粒雷諾數(shù)Rep（以顆粒直徑為參考長(zhǎng)度的雷諾數(shù)）、克努森數(shù)Kn和 馬赫數(shù)Map（ 定義為 |vr|/af,af為顆粒位置處流體的聲速，vr為流體與顆粒之間的相對(duì)速度）范圍。而顆粒常見(jiàn)的Stokes[28]均勻不可壓縮蠕流條件下單顆顆粒的阻力公式（FD=?6πμavr，a=dp/2為顆粒半徑， μ是流體相介質(zhì)的動(dòng)力黏度）已不適用。因此，參考Tedeschi等[29]的工作，本文計(jì)算所用的阻力公式（15、16）考慮了可壓縮性影響和稀薄效應(yīng)：

當(dāng)顆粒位于流場(chǎng)強(qiáng)剪切區(qū)域時(shí)，會(huì)受到Saffman（薩夫曼）升力。Saffman[30]基于無(wú)界線性剪切流情況，對(duì)單個(gè)球體受力問(wèn)題進(jìn)行研究，通過(guò)理論推導(dǎo)發(fā)現(xiàn)了這種垂向流體誘導(dǎo)的剪切升力，這一升力因此被命名為Saffman 升力：

其中，升力系數(shù)J=1， ω為顆粒所在位置處流體渦量。

Li 等[9]指出，熱泳力Ft在弓形激波附近和溫度梯度較高的熱邊界層中可能很重要, 對(duì)于亞微米顆粒，布朗運(yùn)動(dòng)及布朗力FBi也可能不可忽視。本文參考Li 等的工作，在顆粒運(yùn)動(dòng)模型中也增加了熱泳力和布朗力，其中，熱泳力采用了Yamamoto 等[31]的無(wú)量綱形式的公式：

式中：Aw、A0、Hw及H0均 為克努森數(shù)的函數(shù)，k? 是顆粒和流體的導(dǎo)熱系數(shù)之比。布朗力可以被考慮為一個(gè)高斯白噪聲隨機(jī)過(guò)程，其譜強(qiáng)度Snij的表達(dá)式為：

式中：Tf是流體溫度， δij是 克羅內(nèi)克符號(hào)，kB是玻爾茲曼常數(shù)，Cc是斯托克斯-坎寧安修正項(xiàng)。最終，布朗力表示為：

其中， ξi是零均值的單位方差無(wú)關(guān)的高斯隨機(jī)數(shù)[17]。注意到，在以往的歐拉-拉格朗日模擬中，大部分研究認(rèn)為拖曳力是占主導(dǎo)地位的，因此，往往在模擬顆粒軌跡時(shí)也僅考慮了拖曳力[20]。Li 等[9]在平板邊界層層流顆粒軌跡的模擬中發(fā)現(xiàn)，熱泳力在其模擬的參數(shù)范圍內(nèi)（自由來(lái)流馬赫數(shù)小于3.01）是不重要的。同時(shí)，他們的模擬結(jié)果表明小粒徑（dp=0.05μm）顆粒的布朗運(yùn)動(dòng)比大粒徑（dp=1.0μm）顆粒的布朗運(yùn)動(dòng)嚴(yán)重得多，且隨著主流馬赫數(shù)的增大，布朗運(yùn)動(dòng)將變得更加劇烈。但無(wú)論如何，升力的影響非常重要，當(dāng)顆粒進(jìn)入邊界層時(shí)，考慮薩夫曼升力的顆粒軌跡顯著低于未考慮薩夫曼升力的顆粒軌跡。

1.4 流場(chǎng)速度插值方法

Gimenez 等[32]提出了顆粒位置處的插值方法，將流場(chǎng)離散化，并將單元中的流場(chǎng)變量 ?存儲(chǔ)在單元中心xc處 ，即 ?c=?(xc)。 考 慮點(diǎn)xc周 圍 任 意 位 置 處 ?的Taylor 級(jí)數(shù)展開(kāi)，并略去高階項(xiàng)，得到任意位置xp處?(xp)的值：

1.5 顆粒邊界條件

所開(kāi)發(fā)代碼中，顆粒邊界條件包括入口、出口邊界條件，周期邊界條件和固壁完全反彈邊界條件三類。顆粒壁面碰撞和反彈示意圖見(jiàn)圖1。顆粒的入口條件包括顆粒的位置坐標(biāo)和速度信息。當(dāng)顆粒中心超出出口邊界或者顆粒中心距出口邊界距離小于半徑時(shí)，從存儲(chǔ)顆粒的數(shù)組中刪除該顆粒信息，如果將顆粒出口速度信息保留并重置其位置坐標(biāo)，可實(shí)現(xiàn)周期邊界條件。注意后者僅適用于具有周期邊界條件的流場(chǎng)。

圖1 計(jì)算網(wǎng)格點(diǎn)內(nèi)顆粒壁面碰撞和反彈示意圖Fig.1 Schematic of particle-wall collision and rebound in the calculation cell

顆粒中心與壁面距離小于顆粒半徑時(shí)，在碰撞點(diǎn)所在壁面網(wǎng)格單元平面上將速度矢量進(jìn)行反射變換，更新顆粒的位置信息與速度信息，并繼承其他物理參數(shù)，反射過(guò)程無(wú)動(dòng)量損失。

2 顆粒求解可靠性檢驗(yàn)

本節(jié)將用3 個(gè)算例對(duì)代碼進(jìn)行驗(yàn)證測(cè)試，測(cè)試算例分別為：1）泊肅葉流動(dòng)中顆粒運(yùn)動(dòng)算例，用于檢驗(yàn)顆粒插值方法；2）超聲速絕熱平板邊界層中單顆顆粒運(yùn)動(dòng)算例，用于檢驗(yàn)不同受力模型對(duì)顆粒運(yùn)行軌跡的影響；3）等溫平板邊界層顆粒群對(duì)熱流的影響算例，驗(yàn)證代碼雙向耦合部分的準(zhǔn)確性。

2.1 泊肅葉流動(dòng)中顆粒運(yùn)動(dòng)

計(jì)算二維穩(wěn)態(tài)泊肅葉流動(dòng)中單個(gè)顆粒在流場(chǎng)中的運(yùn)動(dòng)速度和軌跡的模擬示意圖見(jiàn)圖2。圖中，T∞為來(lái)流總溫和等溫壁的壁面溫度，U∞為來(lái)流速度，pin和pout分別為入口和出口總壓強(qiáng)。

圖2 計(jì)算域與邊界條件示意圖Fig.2 Schematic diagram of computational domain and boundary conditions

根據(jù)解析解，該流動(dòng)中速度均具有拋物線形式解：

其中Ly為半槽道高度。當(dāng)慣性顆粒（運(yùn)動(dòng)方程中僅考慮Stokes 拖曳力）以流向初速度up0=0和垂向初速度vp0≠0在 初始位置處 (xp0,yp0)處釋放進(jìn)入流動(dòng)中時(shí)，其速度 (up(t),vp(t))和 位置坐標(biāo) (xp(t),yp(t))的理論 解如下：

計(jì)算域總長(zhǎng)度設(shè)置為 0 mm ≤Lx≤20 mm，高度?0.2 mm ≤Ly≤0.2 mm。在流向上采用均勻網(wǎng)格，垂向網(wǎng)格采用沿壁面指數(shù)加密并關(guān)于x軸對(duì)稱。流向200 個(gè)網(wǎng)格節(jié)點(diǎn)，垂向21 個(gè)網(wǎng)格節(jié)點(diǎn)。來(lái)流溫度T∞=334.5 K ， 來(lái)流速度U∞=11.7 m/s，來(lái)流馬赫數(shù)Ma∞=0.031 9 ， 入口壓強(qiáng)pin=101.625 kPa，普朗特?cái)?shù)Pr=0.72 ， 比 熱比 γ=1.4 。 顆 粒 密度 ρp=2 000 kg/m3，粒徑dp=1×10?5m。 顆粒釋放的初始速度為up0=0，vp0=1.0 m/s， 初始位置為xp0=11 mm，yp0=0.01 mm，該位置處流動(dòng)不受入口影響，充分發(fā)展，平均壓力梯度為 dp/dx=?10.149 kPa。圖3 給出了顆粒軌跡和顆粒速度隨時(shí)間變化的數(shù)值解和解析解的對(duì)比，其中圖3（a）的縱坐標(biāo)表示無(wú)量綱的顆粒速度，即up/U∞和vp/U∞；圖3（b）的縱坐標(biāo)表示為無(wú)量綱的顆粒位置坐標(biāo)，即xp/L和yp/L，從圖中可以看出顆粒解析解和數(shù)值解驗(yàn)證良好。

圖3 數(shù)值模擬的顆粒速度和軌跡與理論結(jié)果對(duì)比Fig.3 Comparison of particle velocities and trajectories from numerical simulations with theoretical results

2.2 超聲速絕熱平板邊界層中顆粒運(yùn)動(dòng)

Li 等[9]對(duì)超聲速平板層流邊界層的亞微米顆粒在不同受力模型下的運(yùn)動(dòng)進(jìn)行了分析，并詳細(xì)討論了顆粒初始位置、流動(dòng)馬赫數(shù)、顆粒尺寸和密度對(duì)于軌跡的影響。本文在與文獻(xiàn)[9]相同的條件下模擬顆粒運(yùn)動(dòng)，具體流動(dòng)參數(shù)見(jiàn)表1。其中，以單位長(zhǎng)度表示的有量綱流體松弛時(shí)間為 τ?f,∞=1/U∞，以單相流邊界層流動(dòng)出口處邊界層厚度 δout表示的有量綱流體松弛時(shí)間為 τf,out=δout/U∞=δoutτ?f,∞， 以 δout表示的流體雷諾數(shù)為Reout=ρ∞U∞δout/μ∞=δoutRe?∞。

表1 平板邊界層流動(dòng)參數(shù)Table 1 Parameters of the flat-plate boundary layer flow

二維超聲速平板層流邊界層流動(dòng)計(jì)算域大小為50 mm×50 mm，在x、y方向均采用330 個(gè)網(wǎng)格節(jié)點(diǎn)，并在近壁處以及平板前緣激波位置處進(jìn)行加密。與文獻(xiàn)[9]不同的是，模擬采用理想氣體狀態(tài)方程，而非Redlich–Kwong 方程。圖4 給出了模擬流場(chǎng)速度廓線與Crocco[33]解析解的對(duì)比，其中橫坐標(biāo)表示為無(wú)量綱高度，縱坐標(biāo)為無(wú)量綱的流體速度，可以看到，數(shù)值模擬結(jié)果與理論解吻合得很好。

圖4 流場(chǎng)速度廓線計(jì)算值與解析解對(duì)比Fig.4 Comparison between calculated values and analytical solutions of velocity profile

直接數(shù)值模擬平板邊界層達(dá)到穩(wěn)定后，在入口邊界層發(fā)展前緣處釋放單顆顆粒（見(jiàn)圖4 中插圖）。釋放顆粒的具體參數(shù)見(jiàn)表2。其中，顆粒的Stokes 數(shù)表示為S t=τp/τf， τp表示顆粒松弛時(shí)間。這里以流體來(lái)流黏度 μ∞和 δout表示的有量綱顆粒松弛時(shí)間尺度為τp,out=ρpdp2/(18μ∞)，對(duì) 應(yīng) 的 顆 粒Stokes 數(shù) 表 示 為S tout=τp,out/τf,out，顆粒在入口邊界初始釋放高度為yp0=0.15 mm，顆粒釋放的初始速度為當(dāng)?shù)亓鲌?chǎng)速度。

表2 模擬的顆粒參數(shù)Table 2 Parameters of a single particle released on plate

顆粒受拖曳力、薩夫曼力、熱泳力和布朗力的共同作用。模擬得到的顆粒運(yùn)動(dòng)軌跡見(jiàn)圖5，可以看到，考慮薩夫曼升力時(shí)，顆粒運(yùn)動(dòng)軌跡呈向下的趨勢(shì)，這與Li 等[9]的結(jié)果一致。注意到，由于氣體狀態(tài)方程的不同，顆粒軌跡會(huì)和參考文獻(xiàn)存在一定區(qū)別，但從定性上分析，顆粒在受到不同力的情況下，總體運(yùn)動(dòng)趨勢(shì)是一致的。

圖5 平板邊界層中運(yùn)動(dòng)顆粒y 坐標(biāo)隨時(shí)間變化Fig.5 Temporal variation of coordinate values of moving particles along wall-normal direction in boundary layer flow

同時(shí)，這一流動(dòng)情況下其他作用力對(duì)軌跡的影響見(jiàn)圖6?？梢詮膱D中看出，在來(lái)流馬赫數(shù)為2.05 時(shí)，在拖曳力和升力的基礎(chǔ)上，進(jìn)一步考慮熱泳力以及布朗力，對(duì)粒徑為1 um 的顆粒的運(yùn)動(dòng)軌跡并無(wú)顯著影響，因此，熱泳力以及布朗力可以忽略。這也與文獻(xiàn)[9]中的結(jié)論相吻合。

圖6 不同受力情況下顆粒軌跡對(duì)比Fig.6 Comparison of particle trajectories under different force conditions

2.3 等溫平板邊界層中顆粒群對(duì)熱流的影響

Ching 等[35]在采用歐拉-拉格朗日方法進(jìn)行雙向耦合模擬時(shí)，與Wang 等[34]進(jìn)行了對(duì)比。在本研究中，參考Ching 等[35]的方法，模擬Wang 等[34]文獻(xiàn)中的大滑移區(qū)域，并與理論結(jié)果進(jìn)行比較，模擬參數(shù)見(jiàn)表3。其中，Ts為薩瑟蘭常數(shù)，顆粒與流體的質(zhì)量比β 定 義 為ρˉd,∞=mˉdnd是來(lái)流顆粒相的體密度，nd是 單位體積的顆粒數(shù)，mˉd是顆粒的平均質(zhì)量。為了與理論值進(jìn)行對(duì)比，顆粒僅受斯托克斯阻力作用。二維計(jì)算域總流向長(zhǎng)度設(shè)置為 [?0.1, 1.0]mm，高度[0, 0.2] mm 。 在流向上采用靠近x=0的指數(shù)加密網(wǎng)格，在垂向上采用靠近y=0的指數(shù)加密網(wǎng)格，網(wǎng)格節(jié)點(diǎn)數(shù)Nx×Ny=91×56。本算例顆粒初始釋放參考Ching 等[35]方式，見(jiàn)圖7。初始顆粒從x=0、y=0～δout處均勻釋放，假設(shè)初始顆粒速度與無(wú)限遠(yuǎn)來(lái)流速度相同并按照來(lái)流的質(zhì)量比 β控制入口處釋放的顆粒數(shù)。每隔 ΔT=2.71×10?9s重復(fù)釋放；當(dāng)顆粒超出出口時(shí)將顆粒移除。

表3 雙向耦合平板邊界參數(shù)Table 3 Parameters of flow and particles in two-way coupling

圖7 顆粒釋放示意圖Fig.7 Schematic of particle release

圖8 單相與兩相流中無(wú)量綱壁面熱通量Fig.8 Non-dimensionl wall heat flux in particle-free and particle-laden flows

3 顆粒與層流邊界層之間的相互影響

3.1 顆粒對(duì)流體的調(diào)制研究

基于以上的測(cè)試結(jié)果，進(jìn)一步分析了層流邊界層中的顆粒-流場(chǎng)相互作用。由2.2 節(jié)可以得知薩夫曼升力在超聲速平板邊界層的顆粒運(yùn)動(dòng)中是不可忽略的（薩夫曼升力對(duì)于顆粒的運(yùn)動(dòng)軌跡具有一定影響），因此，在本算例中，我們將對(duì)2.3 節(jié)的雙向耦合進(jìn)行驗(yàn)證，具體流動(dòng)參數(shù)見(jiàn)表3。在驗(yàn)證雙向耦合的基礎(chǔ)上將流向計(jì)算域延長(zhǎng)至 2.5 mm，給顆粒附加薩夫曼升力，考慮在更真實(shí)顆粒受力情況下的雙向耦合調(diào)制規(guī)律及其沿程信息。

圖9 給出了不同情況下模擬得到的無(wú)量綱平板壁面熱通量沿程變化?？梢钥闯?，兩相流中顆粒導(dǎo)致壁面熱通量顯著增加。注意到，圖中僅考慮拖曳力的顆粒兩相流相比純氣體的無(wú)量綱熱流已經(jīng)有顯著增大。而相比不考慮薩夫曼升力的理論結(jié)果，薩夫曼升力對(duì)壁面熱流有進(jìn)一步的增強(qiáng)作用，并且增強(qiáng)效應(yīng)沿流向不斷增長(zhǎng)。在x?=0.35處，不含薩夫曼升力的算例相對(duì)于純流體壁面，壁面熱通量增強(qiáng)達(dá)到了52.8%，附加薩夫曼升力后，相對(duì)于不加薩夫曼升力的情況又增加了 9.3%。

圖9 不同情況中無(wú)量綱壁面熱通量沿平板流向變化Fig.9 Variation of non-dimensional wall heat flux along flow direction under different conditions

接下來(lái)通過(guò)邊界層厚度的變化來(lái)分析熱流和摩阻變化的原因。圖11 給出了不同模擬中的邊界層位移厚度，其中，位移厚度 δ′的 無(wú)量綱形式為δ?=δ′/(λ∞Re1∞/2)。從圖11 可以看出，由于顆粒以來(lái)流速度進(jìn)入邊界層，邊界層內(nèi)的流體通過(guò)和顆粒的相互作用從顆粒相獲得動(dòng)量，進(jìn)而加速，導(dǎo)致邊界層厚度減小，從而進(jìn)一步導(dǎo)致邊界層內(nèi)的流體速度梯度和溫度梯度增加，這是導(dǎo)致圖9 和圖10 含顆粒流相較于純流體摩阻與熱流劇烈增加的主要原因。

圖10 不同情況中壁面摩阻系數(shù)沿平板流向變化Fig.10 Variation of skin friction coefficient along flow direction under different conditions

圖11 不同情況下邊界層位移厚度對(duì)比Fig.11 Variation of displacement thickness along flow direction under different conditions

3.2 流體對(duì)顆粒輸運(yùn)的影響研究

首先，討論薩夫曼升力對(duì)顆粒運(yùn)動(dòng)的影響，由式（22）可知，薩夫曼升力的方向主要由顆粒-流體間的滑移速度以及渦量確定。因此，圖12 給出了典型單個(gè)顆粒在進(jìn)入邊界層后，顆粒位置處流向上的滑移速度 ΔU?及 流體渦量 ω?x隨 時(shí)間的變化（其中 ΔU?和 ω?x均用來(lái)流速度無(wú)量綱化）?？梢钥吹筋w粒以初始來(lái)流速度進(jìn)入邊界層后，會(huì)受到一個(gè)向下的薩夫曼升力，這是由于顆粒與流體流向滑移速度為負(fù)導(dǎo)致的，如圖12（a）。

圖12 單個(gè)顆粒流向滑移速度和渦量隨時(shí)間步變化Fig.12 Variation of individual particle slip velocity and vorticity with time steps

其次，為分析流場(chǎng)對(duì)顆粒運(yùn)動(dòng)的影響，進(jìn)一步統(tǒng)計(jì)了流場(chǎng)中顆粒流向、垂向速度分布uˉp、vˉp和顆粒濃度 ?p。由于顆粒對(duì)流體影響的沿程變化，且在大滑移區(qū)含顆粒與不含顆粒情況下的流動(dòng)性質(zhì)差別不顯著，為簡(jiǎn)化分析，首先沿流向?qū)⒂?jì)算域劃分為四個(gè)部分，分 別 為x*= [0, 0.087 5],x*= [0.087 5, 0.175 0],x*=[0.175 0, 0.262 5],x*= [0.262 5, 0.350 0]，對(duì)這 四個(gè)區(qū)域的顆粒分別進(jìn)行統(tǒng)計(jì)平均，uˉp、vˉp和 ?p的統(tǒng)計(jì)結(jié)果分別見(jiàn)圖13、圖14 和圖15。圖中，縱坐標(biāo)y均以單相流邊界層流動(dòng)出口處邊界層最大厚度 δout進(jìn)行無(wú)量綱化，而速度uˉp和vˉp分 別采用U∞無(wú) 量綱化，表示為uˉ?p、vˉ?p。為了更直觀地觀察顆粒與流體的動(dòng)量交換情況，將相同流向段流體的速度uˉ?f也同時(shí)繪制在圖中，并用虛線表示。

圖13 不同流向段顆粒和流體平均流向速度分布廓線Fig.13 Profiles of time-averaged streamwise veloicity of particle and fluid in different streamwise segments

圖14 不同流向段顆粒和流體平均垂向速度分布廓線Fig.14 Profiles of time-averaged vertical veloicity of particle and fluid in different streamwise segments

圖15 不同流向段顆粒濃度沿高度分布Fig.15 Distribution of particle concentration along wall-normal direction in different streamwise segments

圖13 給出了不同流向段顆粒流向速度分布廓線。從圖中可以觀察到，顆粒速度沿著流向不斷衰減，并且越靠近壁面，衰減程度越大，這一速度變化和流場(chǎng)的速度廓線相類似。在計(jì)算域的四個(gè)流向段中，顆粒速度均大于流體速度，因此當(dāng)顆粒進(jìn)入邊界層時(shí)，由于速度高于流體速度，會(huì)產(chǎn)生一個(gè)方向向下的薩夫曼升力，可能會(huì)使顆粒向壁面運(yùn)動(dòng)，這一點(diǎn)也與圖12 觀察到的結(jié)果一致。

同時(shí)，流體的垂向速度也可能對(duì)顆粒的垂向輸運(yùn)造成影響。圖14 給出了不同流向段顆粒相和流體相的平均垂向速度分布廓線。可以看到，在靠近邊界層前緣段（第一段），由于存在激波，流場(chǎng)垂向速度沿高度增加；在第二段，由于遠(yuǎn)離了激波區(qū)域，流場(chǎng)的垂向速度驟減并且隨著沿程逐漸降低。根據(jù)圖14 流場(chǎng)的垂向沿程速度分布，可以觀察到流場(chǎng)的垂向速度對(duì)顆粒垂向速度起著加速作用。同時(shí)，由于在近壁處顆粒會(huì)受到更強(qiáng)的升力作用，顆粒在近壁區(qū)域的垂向速度會(huì)存在負(fù)值。

對(duì)圖15 中顆粒沿程濃度進(jìn)行分析可以發(fā)現(xiàn)，在近壁區(qū)域，顆粒會(huì)存在壁面的累積效應(yīng)，這是由于近壁處邊界層內(nèi)流體速度梯度使得顆粒產(chǎn)生靠近壁面的垂向速度導(dǎo)致的（見(jiàn)圖14），并且這種累積效應(yīng)會(huì)沿流向逐漸加劇（見(jiàn)圖15）。在遠(yuǎn)離壁面區(qū)域（ y ＞0.5δout）可以觀察到，顆粒垂向速度場(chǎng)相比于近壁區(qū)域（y＜0.5δout）沿高度并沒(méi)有明顯的變化，因此在垂向統(tǒng)計(jì)域內(nèi)不會(huì)出現(xiàn)類似于近壁區(qū)域的明顯聚集現(xiàn)象。其次，分析顆粒沿流向的分布，可以觀察到，顆粒平均濃度沿流向逐漸增加，這是由于在顆粒進(jìn)入邊界層后，邊界層內(nèi)的低速流體在流向上不斷將顆粒減速，導(dǎo)致沿程顆粒平均流向速度降低（見(jiàn)圖13），從而造成沿程顆粒沿流向累積。

4 結(jié) 論

基于PHengLEI 軟件開(kāi)源框架，采用緊耦合的模式開(kāi)發(fā)了點(diǎn)顆粒模擬求解器。文中介紹了點(diǎn)顆粒模型的具體求解方法，包括顆?？刂品匠毯褪芰δＰ?、顆粒追蹤算法和邊界條件、顆粒位置處流動(dòng)信息插值等。分別采用泊肅葉流動(dòng)中顆粒運(yùn)動(dòng)軌跡和速度、超聲速絕熱平板邊界層中顆粒的運(yùn)動(dòng)和顆粒對(duì)等溫平板邊界層熱流的影響三個(gè)算例，對(duì)所設(shè)計(jì)的求解器進(jìn)行驗(yàn)證。計(jì)算表明，顆粒求解器對(duì)可壓縮流動(dòng)中顆粒的運(yùn)動(dòng)軌跡、雙向耦合兩相流中壁面熱流的預(yù)測(cè)與文獻(xiàn)中的計(jì)算和理論結(jié)果吻合，證明了求解器的可靠性以及對(duì)于可壓縮兩相流的模擬能力。

在同時(shí)考慮拖曳力和升力的情況下，雙向耦合模擬了馬赫數(shù)為1.5 的等溫?zé)o滑移壁面上的層流邊界層兩相流。結(jié)果表明，顆粒和邊界層流動(dòng)存在顯著的相互作用，具體表現(xiàn)為當(dāng)顆粒以來(lái)流速度和溫度進(jìn)入邊界層時(shí)，通過(guò)和流體的能量交換，顆粒的存在會(huì)增強(qiáng)邊界層的壁面熱通量，通過(guò)和邊界層的動(dòng)量交換，極大提高邊界層的壁面摩阻。顆粒的存在還會(huì)減少邊界層的位移厚度。同時(shí)，超聲速層流邊界層流場(chǎng)沿流向變化也會(huì)影響顆粒沿程分布。邊界層內(nèi)顆粒速度大于流體速度而產(chǎn)生向下的剪切升力，使顆粒垂向速度為負(fù)，向壁面堆積，這種堆積效應(yīng)沿流向不斷增強(qiáng)，導(dǎo)致近壁顆粒濃度增加。

后續(xù)，將進(jìn)一步持續(xù)開(kāi)展顆粒求解器功能的擴(kuò)充、開(kāi)發(fā)和優(yōu)化。同時(shí)，基于求解器研究可壓縮邊界層兩相流的氣動(dòng)特性和顆粒相運(yùn)動(dòng)、分布規(guī)律。討論在不同的流動(dòng)參數(shù)（馬赫數(shù)、雷諾數(shù)、普朗特?cái)?shù)、壁溫等）和顆粒參數(shù)（體積分?jǐn)?shù)、混合粒徑、密度比、釋放方式、重力等）情況下顆粒-流體的相互作用，為高超聲速飛行器燒蝕顆粒邊界層特性等問(wèn)題提供相應(yīng)參考。