基于風雷軟件的可壓縮氣固顆粒兩相流計算

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

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

0 引 言

可壓縮顆粒兩相流常見于超燃發動機[1-2]、沙塵或雨中飛行器航行[3-4]和混合物沖擊/爆炸[5]等。另一個典型的例子是高超聲速飛行器再入大氣過程中，熱防護材料高溫燒蝕“層裂”過程中產生顆粒[6-7]，與飛行器周圍氣流相互作用。因為單個顆粒可能在稀薄或連續氣體中經歷亞聲速到高超聲速的流態，這些類型氣流中的顆粒運動建模很復雜。同時，由于兩相系統的復雜性和多尺度性，用理論和實驗方法進行研究存在諸多困難，數值模擬可直觀地展示顆粒內部的物理現象，因此可發揮較大作用。

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

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

本文基于國產自主研發的計算流體力學（Computational Fluid Dynamics，CFD）軟件PHengLEI（風雷）[19]，開發可壓縮邊界層顆粒多相流程序，為國家數值風洞軟件開發提供技術儲備和集成模塊?；诘湫退憷龑λ_發的模塊進行檢驗，并初步研究了層流邊界層中顆粒的運動行為及其對流動性質的影響。

1 數值模型和方法

風雷軟件是中國空氣動力研究與發展中心研發的面向流體工程的混合CFD 平臺。風雷軟件的計算范圍覆蓋低速、亞跨聲速和高超聲速，是一款具有完全自主知識產權的面向工程應用和學術研究的通用CFD 軟件，借鑒了面向對象的大型軟件設計理念，采用C++語言編程。風雷軟件的整體框架可參考文獻[19-20]，這里僅介紹與多相流開發相關的數值模型和方法。

1.1 流體控制方程

本文考慮了顆粒-流體的雙向耦合，在求解過程中，通過在流體相控制方程中添加顆粒的動量、能量源項的方式，來體現顆粒對流體相的反作用。實際程序計算中，通過定義歐拉變量，將每一時間步各單元的反饋源項存儲到單元中心處。笛卡爾坐標系下，包含顆粒源項的守恒形式流動控制方程為：

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

其中，Fp,i是顆粒受到的力，Tp為 顆粒的溫度，Qp,i是顆粒釋放的熱量，Gp,i是 顆粒對流體做的功，dp是顆粒的粒徑，Nu是努塞特數，np表示計算網格單元體積內的顆粒數， ΔV表示網格單元的體積，up,i是顆粒速度。

1.2 計算相關參數

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

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

1.3 顆?？刂品匠?/h3>

本文離散顆粒相計算的控制方程分別為顆粒位移、顆粒平動、顆粒轉動及顆粒溫度控制方程，見式（10～13）：

其中，Xp、vp、 ωp是顆粒的運動學參數，分別是位移、速度、角速度，Tp為 顆粒溫度，Ip是顆粒的慣性矩，Fp和Tc分 別是顆粒受到的力及扭矩，cp,p為顆粒相物質的比定壓熱容。

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

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

當顆粒位于流場強剪切區域時，會受到Saffman（薩夫曼）升力。Saffman[30]基于無界線性剪切流情況，對單個球體受力問題進行研究，通過理論推導發現了這種垂向流體誘導的剪切升力，這一升力因此被命名為Saffman 升力：

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

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

式中：Aw、A0、Hw及H0均 為克努森數的函數，k? 是顆粒和流體的導熱系數之比。布朗力可以被考慮為一個高斯白噪聲隨機過程，其譜強度Snij的表達式為：

式中：Tf是流體溫度， δij是 克羅內克符號，kB是玻爾茲曼常數，Cc是斯托克斯-坎寧安修正項。最終，布朗力表示為：

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

1.4 流場速度插值方法

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

1.5 顆粒邊界條件

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

圖1 計算網格點內顆粒壁面碰撞和反彈示意圖Fig.1 Schematic of particle-wall collision and rebound in the calculation cell

顆粒中心與壁面距離小于顆粒半徑時，在碰撞點所在壁面網格單元平面上將速度矢量進行反射變換，更新顆粒的位置信息與速度信息，并繼承其他物理參數，反射過程無動量損失。

2 顆粒求解可靠性檢驗

本節將用3 個算例對代碼進行驗證測試，測試算例分別為：1）泊肅葉流動中顆粒運動算例，用于檢驗顆粒插值方法；2）超聲速絕熱平板邊界層中單顆顆粒運動算例，用于檢驗不同受力模型對顆粒運行軌跡的影響；3）等溫平板邊界層顆粒群對熱流的影響算例，驗證代碼雙向耦合部分的準確性。

2.1 泊肅葉流動中顆粒運動

計算二維穩態泊肅葉流動中單個顆粒在流場中的運動速度和軌跡的模擬示意圖見圖2。圖中，T∞為來流總溫和等溫壁的壁面溫度，U∞為來流速度，pin和pout分別為入口和出口總壓強。

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

根據解析解，該流動中速度均具有拋物線形式解：

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

計算域總長度設置為 0 mm ≤Lx≤20 mm，高度?0.2 mm ≤Ly≤0.2 mm。在流向上采用均勻網格，垂向網格采用沿壁面指數加密并關于x軸對稱。流向200 個網格節點，垂向21 個網格節點。來流溫度T∞=334.5 K ， 來流速度U∞=11.7 m/s，來流馬赫數Ma∞=0.031 9 ， 入口壓強pin=101.625 kPa，普朗特數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，該位置處流動不受入口影響，充分發展，平均壓力梯度為 dp/dx=?10.149 kPa。圖3 給出了顆粒軌跡和顆粒速度隨時間變化的數值解和解析解的對比，其中圖3（a）的縱坐標表示無量綱的顆粒速度，即up/U∞和vp/U∞；圖3（b）的縱坐標表示為無量綱的顆粒位置坐標，即xp/L和yp/L，從圖中可以看出顆粒解析解和數值解驗證良好。

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

2.2 超聲速絕熱平板邊界層中顆粒運動

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

表1 平板邊界層流動參數Table 1 Parameters of the flat-plate boundary layer flow

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

圖4 流場速度廓線計算值與解析解對比Fig.4 Comparison between calculated values and analytical solutions of velocity profile

直接數值模擬平板邊界層達到穩定后，在入口邊界層發展前緣處釋放單顆顆粒（見圖4 中插圖）。釋放顆粒的具體參數見表2。其中，顆粒的Stokes 數表示為S t=τp/τf， τp表示顆粒松弛時間。這里以流體來流黏度 μ∞和 δout表示的有量綱顆粒松弛時間尺度為τp,out=ρpdp2/(18μ∞)，對 應 的 顆 粒Stokes 數 表 示 為S tout=τp,out/τf,out，顆粒在入口邊界初始釋放高度為yp0=0.15 mm，顆粒釋放的初始速度為當地流場速度。

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

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

圖5 平板邊界層中運動顆粒y 坐標隨時間變化Fig.5 Temporal variation of coordinate values of moving particles along wall-normal direction in boundary layer flow

同時，這一流動情況下其他作用力對軌跡的影響見圖6?？梢詮膱D中看出，在來流馬赫數為2.05 時，在拖曳力和升力的基礎上，進一步考慮熱泳力以及布朗力，對粒徑為1 um 的顆粒的運動軌跡并無顯著影響，因此，熱泳力以及布朗力可以忽略。這也與文獻[9]中的結論相吻合。

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

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

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

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

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

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

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

3.1 顆粒對流體的調制研究

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

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

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

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

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

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

3.2 流體對顆粒輸運的影響研究

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

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

其次，為分析流場對顆粒運動的影響，進一步統計了流場中顆粒流向、垂向速度分布uˉp、vˉp和顆粒濃度 ?p。由于顆粒對流體影響的沿程變化，且在大滑移區含顆粒與不含顆粒情況下的流動性質差別不顯著，為簡化分析，首先沿流向將計算域劃分為四個部分，分 別 為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]，對這 四個區域的顆粒分別進行統計平均，uˉp、vˉp和 ?p的統計結果分別見圖13、圖14 和圖15。圖中，縱坐標y均以單相流邊界層流動出口處邊界層最大厚度 δout進行無量綱化，而速度uˉp和vˉp分 別采用U∞無 量綱化，表示為uˉ?p、vˉ?p。為了更直觀地觀察顆粒與流體的動量交換情況，將相同流向段流體的速度uˉ?f也同時繪制在圖中，并用虛線表示。

圖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 給出了不同流向段顆粒流向速度分布廓線。從圖中可以觀察到，顆粒速度沿著流向不斷衰減，并且越靠近壁面，衰減程度越大，這一速度變化和流場的速度廓線相類似。在計算域的四個流向段中，顆粒速度均大于流體速度，因此當顆粒進入邊界層時，由于速度高于流體速度，會產生一個方向向下的薩夫曼升力，可能會使顆粒向壁面運動，這一點也與圖12 觀察到的結果一致。

同時，流體的垂向速度也可能對顆粒的垂向輸運造成影響。圖14 給出了不同流向段顆粒相和流體相的平均垂向速度分布廓線?？梢钥吹?，在靠近邊界層前緣段（第一段），由于存在激波，流場垂向速度沿高度增加；在第二段，由于遠離了激波區域，流場的垂向速度驟減并且隨著沿程逐漸降低。根據圖14 流場的垂向沿程速度分布，可以觀察到流場的垂向速度對顆粒垂向速度起著加速作用。同時，由于在近壁處顆粒會受到更強的升力作用，顆粒在近壁區域的垂向速度會存在負值。

對圖15 中顆粒沿程濃度進行分析可以發現，在近壁區域，顆粒會存在壁面的累積效應，這是由于近壁處邊界層內流體速度梯度使得顆粒產生靠近壁面的垂向速度導致的（見圖14），并且這種累積效應會沿流向逐漸加?。ㄒ妶D15）。在遠離壁面區域（ y ＞0.5δout）可以觀察到，顆粒垂向速度場相比于近壁區域（y＜0.5δout）沿高度并沒有明顯的變化，因此在垂向統計域內不會出現類似于近壁區域的明顯聚集現象。其次，分析顆粒沿流向的分布，可以觀察到，顆粒平均濃度沿流向逐漸增加，這是由于在顆粒進入邊界層后，邊界層內的低速流體在流向上不斷將顆粒減速，導致沿程顆粒平均流向速度降低（見圖13），從而造成沿程顆粒沿流向累積。

4 結 論

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

在同時考慮拖曳力和升力的情況下，雙向耦合模擬了馬赫數為1.5 的等溫無滑移壁面上的層流邊界層兩相流。結果表明，顆粒和邊界層流動存在顯著的相互作用，具體表現為當顆粒以來流速度和溫度進入邊界層時，通過和流體的能量交換，顆粒的存在會增強邊界層的壁面熱通量，通過和邊界層的動量交換，極大提高邊界層的壁面摩阻。顆粒的存在還會減少邊界層的位移厚度。同時，超聲速層流邊界層流場沿流向變化也會影響顆粒沿程分布。邊界層內顆粒速度大于流體速度而產生向下的剪切升力，使顆粒垂向速度為負，向壁面堆積，這種堆積效應沿流向不斷增強，導致近壁顆粒濃度增加。

后續，將進一步持續開展顆粒求解器功能的擴充、開發和優化。同時，基于求解器研究可壓縮邊界層兩相流的氣動特性和顆粒相運動、分布規律。討論在不同的流動參數（馬赫數、雷諾數、普朗特數、壁溫等）和顆粒參數（體積分數、混合粒徑、密度比、釋放方式、重力等）情況下顆粒-流體的相互作用，為高超聲速飛行器燒蝕顆粒邊界層特性等問題提供相應參考。