李青 余釗圣 劉朋欣 李婷婷 陳堅強 袁先旭,?

北京大學(xué)學(xué)報(自然科學(xué)版) 第60卷 第1期 2024年1月

Acta Scientiarum Naturalium Universitatis Pekinensis, Vol. 60, No. 1 (Jan. 2024)

10.13209/j.0479-8023.2023.022

國家自然科學(xué)基金(12272398)、中國博士后科學(xué)基金(2021MD703970)和國家重點研發(fā)計劃(2019YFA0405200)資助

2022–12–12;

2023–01–29

高焓流動中的可壓縮顆粒兩相流并行求解器:數(shù)值方法及其驗證

李青1,2余釗圣3劉朋欣1李婷婷4陳堅強1袁先旭1,?

1.空天飛行空氣動力科學(xué)與技術(shù)全國重點實驗室, 中國空氣動力研究與發(fā)展中心, 綿陽 621000; 2.天目山實驗室, 杭州 310000; 3.浙江大學(xué)航空航天學(xué)院, 杭州 310027; 4.西安交通大學(xué)化學(xué)工程與技術(shù)學(xué)院, 西安 710049; ?通信作者, E-mail: yuanxianxu@cardc.cn

在極端力學(xué)環(huán)境下可壓縮點力顆粒兩相流理論方程的基礎(chǔ)上, 提出一種基于動態(tài)鏈表數(shù)組的并行顆粒求解器, 使其與可壓縮流體求解器耦合。與基于歐拉坐標系的攜帶流體相不同, 基于拉格朗日坐標系的求解器采用動態(tài)鏈表數(shù)組對彌散顆粒相進行內(nèi)存分配, 可以解決因使用全局數(shù)組解決拉格朗日/歐拉坐標系轉(zhuǎn)換帶來的彌散顆粒群內(nèi)存利用率低和計算效率低下問題。最后對多物理效應(yīng)的可壓縮顆粒兩相流求解器進行驗證, 并在馬赫數(shù)漸進趨于零的情況下, 對兩個不可壓縮槽道顆粒湍流標模進行驗證。

動態(tài)鏈表數(shù)組; 拉格朗日/歐拉坐標系轉(zhuǎn)換; 可壓縮點力顆粒兩相流直接數(shù)值模擬求解器

有關(guān)顆粒兩相流的系統(tǒng)研究可以追溯至航天火箭發(fā)動機燃燒室固體鋁粉的燃燒問題[1]。隨后, 由于化工流化床的工業(yè)需求推進, 空天和太空工程的工業(yè)需求使得可壓縮顆粒兩相流的相關(guān)研究得到迅速的發(fā)展。高速燒蝕飛行器、航空發(fā)動機、大推力火箭發(fā)動機和超燃沖壓發(fā)動機都包含極端力學(xué)過程[2]: 特征溫度為為(103)～(104)K, 特征速度∞為(103) m/s, 特征馬赫數(shù) Ma 為(10), 特征雷諾數(shù) Re∞為(106), 飛行高度為(10)km 等。航空航天中的顆粒兩相流體動力學(xué)包含燃燒[3]、固體顆?；蛞旱闻鲎伯a(chǎn)生電荷[4–5]、熱輻射和熱對流[6]以及冷熱壁誘導(dǎo)的熱泳力[7]等多個子物理過程。多物理場的各個子物理過程同時同場耦合, 使得航空航天設(shè)備的流體力學(xué)問題變得復(fù)雜。傳統(tǒng)的單相可壓縮流動的理論、計算和實驗方法已經(jīng)無法滿足當前緊迫的工業(yè)需求。

1 研究背景

1.1 可壓縮點力顆粒兩相流理論方程

李青等[8]給出極端力學(xué)條件下, 考慮多物理效應(yīng)的可壓縮點力顆粒兩相流直接數(shù)值模擬求解器的理論方程(point-particle direct numerical simulation, PP-DNS), 如式(1)～(4)所示。式(1)和(2)為多物理效應(yīng)下顆粒動力學(xué)和熱力學(xué)的控制方程, 式(3)和(4)為流體動力學(xué)和熱力學(xué)的顆粒兩相流雙向耦合方程。其中, 假設(shè)顆粒為數(shù)學(xué)上無窮小的沒有體積的點源, 其動力學(xué)和熱力學(xué)過程可以通過不同物理過程的理論或半經(jīng)驗公式表示, 并進行線性疊加。

計算流體力學(xué)的核心是求解離散 Navier-Stokes偏微分方程組, 即離散相容方程組, 而不是 Navier-Stokes 方程本身, 因此 PP-DNS 數(shù)值求解器本質(zhì)上是求解一個數(shù)值離散的偏微分方程組[9–10]來刻畫流體力學(xué)過程, 并耦合求解一系列常微分方程組來刻畫顆粒群的動力學(xué)和熱力學(xué)過程。流體力學(xué)方程屬于雙曲守恒類偏微分方程[11–12], 如果數(shù)值離散帶來的數(shù)值誤差不能隨時間衰減, 則會使得原本沿著特征線傳播的信息不再滿足雙曲守恒律的影響域和依賴域關(guān)系, 最終導(dǎo)致數(shù)值求解器發(fā)散。顆粒動力學(xué)方程與流體力學(xué)方程的反饋耦合等效于在數(shù)值求解域的離散 Euler 點上增加了數(shù)值源項, 會增加流體力學(xué)的離散相容方程組的數(shù)值不穩(wěn)定性和數(shù)值求解器的數(shù)值剛性, 使得滿足數(shù)值穩(wěn)定性條件的最大時間步更加苛刻。為了消除數(shù)值不穩(wěn)定性, Pierce[13]采用 Runge-Kutta 三階的多時間子步方式, 將彌散顆粒相的動力學(xué)和熱力學(xué)信息分別反饋至可壓縮流體力學(xué)的動量方程和能量方程。

嚴格地說, 求解顆粒兩相流基本問題應(yīng)該完全基于 Navier-Stokes 方程。以現(xiàn)有的不可壓縮顆粒解析求解器[14]為例, 其采用比顆粒尺寸更小的數(shù)值離散網(wǎng)格, 獲取有限體積顆粒周圍所有的流體力學(xué)信息, 然后通過定義, 積分獲得顆粒受力。這類方法采用 Lagrange 離散點或 Euler 離散點在顆粒體內(nèi)或表面分布, 通過離散點反饋顆粒與攜帶流體的相間動量交換信息, 并通過積分獲得有限體積的顆粒受力。這類方法屬于全分辨顆粒直接數(shù)值模擬求解方法(particle-resolved direct numerical simulation, PR-DNS)。PR-DNS 的優(yōu)點是不進行任何假設(shè), 直接求解 Navier-Stokes 方程, 獲得彌散顆粒與攜帶流體兩相間的動力學(xué)信息(圖 1)。缺點是需要比 PP-DNS 更小的數(shù)值離散網(wǎng)格, 會導(dǎo)致計算量陡增(圖 2); 另一方面, 在相同顆粒數(shù)目的條件下, PR-DNS比 PP-DNS 帶來更多的數(shù)值源項反饋, 從而帶來比 PP-DNS 更強的數(shù)值不穩(wěn)定性。真實工業(yè)問題的顆粒數(shù)目往往是十億量級, 因此 PR-DNS 方法能且僅能用于簡單標模研究[15–16]。李瑞元等[17]給出單個球形顆粒在可壓縮均勻流中的直接數(shù)值模擬, 并驗證了可壓縮圓球擾流阻力表達式。值得一提的是, 當顆粒尺度比最小流動尺度小時, 流體對顆粒的作用力、顆流體對顆粒的作用力以及顆粒對流體的反饋可以近似地表示為點力模型[18–19], 即 PP-DNS 與PR-DNS 有著相似的精度。目前, 可壓縮條件下的多顆粒 PR-DNS 研究尚罕有報道。

表示一般的顆粒相信息, 如顆粒受力Fp或溫度Tp等

圖2 顆粒解析直接數(shù)值模擬求解器與點力顆粒直接數(shù)值模擬求解器的對比

1.2 基于 PP-DNS 的顆粒湍流和顆粒兩相流研究

采用 PP-DNS 作為直接數(shù)值模擬工具的顆粒兩相湍流基礎(chǔ)研究可分為兩類: 一類是用 PP-DNS 方法進行顆粒湍流機理研究; 另一類是用 PP-DNS 提供直接數(shù)值模擬數(shù)據(jù)庫, 用于構(gòu)建工業(yè)軟件模型, 如 PP-LES (point-particle large eddy simulation)。PP-LES 對攜帶流體相采用 LES 模型, 這類方法的核心問題是如何重構(gòu)被粗網(wǎng)格忽略的亞格子脈動, 從而讓彌散顆粒相能“感受”到原本被 LES 忽略的亞格子流體脈動, 復(fù)現(xiàn)彌散顆粒相的統(tǒng)計過程。LES-LES模型又稱為 Euler-Euler 模型, 它是在 PP-LES 的基礎(chǔ)上, 對彌散顆粒相進行連續(xù)性假設(shè), 將彌散顆粒相的統(tǒng)計動力學(xué)過程?；癁橐粋€基于 Euler 坐標的偏微分方程組。連續(xù)性假設(shè)將彌散顆粒相的控制方程從一系列常微分方程組變成一個偏微分方程組, 大大減少了計算量。以 10 億顆粒的工業(yè)流化床為例, 對 PP-DNS 而言, 彌散顆粒相至少需要求解 10億組 6 個自由度變量的常微分方程組, 總共 60 億個非定常變量; 同時需要考慮 10 億個顆粒對攜帶流體相的插值反饋, 計算量巨大, 無法解決工業(yè)級別規(guī)模的真實航空航天問題[20–21]。對 LES-LES 而言, 基于 Euler 坐標的彌散顆粒相偏微分方程組的非定常變量不超過 10 個。因此, PP-LES和 Euler-Euler 模型往往被用于構(gòu)建工業(yè)軟件內(nèi)核, 而兩個模型的構(gòu)建都需要 PP-DNS 提供參考顆粒湍流數(shù)據(jù)庫。因此, PP-DNS 是開展顆粒兩相流工業(yè)軟件研發(fā)的直接數(shù)值模擬工具和重要手段。

Wang 等[22]和 Rosa 等[23]研究重力對各向同性均勻顆粒湍流的影響, 給出重力、顆粒慣性和湍流雷諾數(shù)三者之間的歸一化表達式。Ferrante 等[24]研究微小氣泡在壁湍流中的減阻機理, 發(fā)現(xiàn)彌散氣泡相會調(diào)制減弱湍流雷諾應(yīng)力, 并減小壁面摩阻。Lain 等[25]通過研究氣固顆粒槽道湍流, 發(fā)現(xiàn)壁面粗糙度對湍流脈動統(tǒng)計量的影響不可忽略。Marchioli 等[26]發(fā)現(xiàn)湍流相干結(jié)構(gòu)與顆粒渦泳現(xiàn)象密切相關(guān), 并利用統(tǒng)計采樣的顆粒湍流數(shù)據(jù)闡明湍流上掃和下掠過程與顆粒法向輸運的關(guān)系。Mehrabadi 等[15]指出, 有限尺寸誘導(dǎo)的尾流效應(yīng)只有在顆粒慣性大的時候才需要修正, 對于 St 為 1 量級的顆粒湍流, 采用 PP-DNS 就可以獲得與 PR-DNS 近似的精度。Daitche[27]發(fā)現(xiàn), 顆粒歷史力效應(yīng)會對顆粒的統(tǒng)計分布產(chǎn)生影響。Zamansky 等[28]發(fā)現(xiàn), 受到太陽熱輻射的黑體顆?？梢酝ㄟ^與攜帶流體相之間的相間傳熱, 改變顆粒湍流動力學(xué)過程。Zhao 等[29]推導(dǎo)了顆粒湍流能量平衡方程, 利用 PP-DNS 研究顆粒壁湍流的能量平衡, 闡明顆粒湍流相互作用的基本過程。在此基礎(chǔ)上, Zhou 等[30]通過進一步的研究, 發(fā)現(xiàn)顆粒對壁湍流脈動應(yīng)力的調(diào)制作用是非單調(diào)的。Pan 等[31]發(fā)現(xiàn)顆粒對湍流的做功過程是多尺度且各向異性的: 在平均尺度上, 顆粒在壁湍流外區(qū)吸收平均流場的動能, 且在近壁面為平均流場輸入動能; 在脈動尺度上, 顆粒在流向上為湍流場輸入動能, 在法向和展向上則從湍流場從吸收動能。Li 等[32]研究了慣性顆粒在不可壓縮空間發(fā)展邊界層中的動力學(xué)過程, 發(fā)現(xiàn)在其考察的參數(shù)范圍內(nèi), 慣性顆粒會增加層流邊界層和湍流邊界層的摩阻。近年來, PP-DNS 在可壓縮顆粒湍流的研究方興未艾。Dai等[33]研究慣性顆粒在可壓縮各向同性均勻湍流中的動力學(xué), 發(fā)現(xiàn)顆粒減弱了可壓縮性; 同時, 隨著Stokes 數(shù)的變化, 慣性顆粒非單調(diào)地調(diào)制湍動能。Xiao 等[34]采用單向耦合方法, 研究低馬赫數(shù)可壓縮顆粒壁湍流中的顆粒運動和分布特性, 發(fā)現(xiàn)顆粒傾向性聚團機理與不可壓情況類似, 與“上拋”、“下掃”和流動條帶結(jié)構(gòu)密切相關(guān)。

Wang 等[35]和 Fevrier 等[36]用 PP-DNS 研究各向同性均勻顆粒湍流, 在統(tǒng)計意義上給出顆粒脈動和湍流脈動的 Euler-Euler 預(yù)測模型。顆粒兩相流工業(yè)軟件 Neptune 的 Euler-Euler 就是基于 Fevrier 等[36]的工作開發(fā)的。隨后, Yu 等[37]用 Neptune 數(shù)值預(yù)測化工管道中的磨蝕問題。Marchioli[38]指出, 由于 PP-LES 忽略了亞格子流體脈動, 使得顆粒的渦泳統(tǒng)計過程與 PP-DNS 有差異, 并且, 如何在 PP-LES 中重構(gòu)被 LES 忽略的小尺度脈動, 并重新讓顆?！案惺堋钡? 是構(gòu)建 PP-LES 模型的關(guān)鍵。

2 數(shù)值方法

2.1 流體求解器

本文在極端力學(xué)環(huán)境下可壓縮點力顆粒兩相流理論方程的基礎(chǔ)上, 開發(fā)一種可壓縮流動條件下基于動態(tài)鏈表數(shù)組的 PP-DNS 求解器, 使其能夠與可壓縮流體求解器耦合。

Tian 等[39]開發(fā)了可壓縮點力顆粒兩相流直接數(shù)值模擬求解器并進行驗證, 其基本理論方程考慮了顆粒相間阻力和相間阻力做功。本文開發(fā)的顆粒求解器在數(shù)值方法上與之類似, 并在其基礎(chǔ)上考慮了多物理效應(yīng)[3,8]。

本文采用任意階 Lagrange 插值[40]來實現(xiàn) Euler 坐標點到彌散相顆粒 Lagrange 點的插值轉(zhuǎn)換, 如式(5)和(6)以及圖 3(a)所示。

其中,I(,,,)一般地表示 3 個方向的速度分量(I(,,,),I(,,,)和I(,,,))、可壓縮流體的密度I(,,,)、溫度I(,,,)或動力黏度I(,,,),th表示某方向上的 Lagrange 插值階數(shù)。

2.2 顆粒求解器

將顆粒動力學(xué)以及熱力學(xué)方程寫入流體求解器的 Runge Kutta 三階時間步推進循環(huán)(式(1)和(2))中, 相間阻力、相間阻力做功以及熱對流項對流體相的反饋采用與流體求解器 OpenCFD 同步的時間推進[9,41–42]:

紅色點為Lagrange坐標點, 藍色點為Euler坐標點

圖3 從Euler坐標點到Lagrange坐標點(a)和從Lagrange坐標點到Euler坐標點(b)的插值示意圖

Fig. 3 Illustration of interpolation from Euler point to Lagrange point (a) and from Lagrange point to Euler point (b)

紅色部分是任意位置的MPI分塊

其中,表示 Runge-Kutta 三階時間步推進的子循環(huán)步, 先后取=1, 2, 3; conservation 項表示可壓縮流體力學(xué)數(shù)值離散方程的守恒項; source 項表示顆粒相對流體相的反饋。

真實的航空航天可壓縮顆粒兩相流問題往往更復(fù)雜, 復(fù)雜的偏微分方程會導(dǎo)致其數(shù)值離散相容方程變得具有數(shù)值剛性。從數(shù)值穩(wěn)定性的角度講, 滿足穩(wěn)定性的最大時間步變得越來越小, 使得計算量劇增, 導(dǎo)致可以研究的參數(shù)范圍受到極大的限制。針對此問題, Pierce[13]開發(fā)了時間分步推進結(jié)合壓力泊松方程修正的數(shù)值算法, 得到廣泛應(yīng)用[43–44]。壓力修正算法最早源于不可壓縮單相流動的數(shù)值求解[45–46]。由于流體力學(xué)方程組本質(zhì)上是一種雙曲型和拋物型的偏微分方程組, 其不穩(wěn)定性主要來源于雙曲型偏微分方程的非線性慣性對流項[47]。數(shù)值穩(wěn)定性的獲得與高精度解析湍流小尺度脈動是一對矛盾體: 一方面, 如果引入過多的數(shù)值黏性獲得穩(wěn)定性, 會導(dǎo)致湍流高階統(tǒng)計矩無法準確捕捉湍流特征[48]; 另一方面, 如果采用高精度算法, 數(shù)值剛性會導(dǎo)致最大穩(wěn)定時間步太小, 使得計算沒有效率[10,49–50]。本質(zhì)上, 壓力修正方法等效于引入質(zhì)量守恒條件約束。壓力修正算法的意義在于, 在數(shù)值精度和數(shù)值穩(wěn)定性二者之間達到平衡。

2.3 彌散顆粒群在并行分塊中的通信

圖 4 為顆粒求解器的 MPI 分塊間的通信示意圖。任意位置的 MPI 分塊與其相鄰的 26 個 MPI 分塊依次發(fā)送或接收顆粒相信息, 從而依次修改或增刪動態(tài)鏈表數(shù)組信息。動態(tài)鏈表數(shù)組結(jié)合MPI分塊的通信技術(shù), 提高了通信效率, 避免了 MPI 通信過程中的死鎖。

3 標模驗證算例

3.1 單顆粒問題

3.1.1插值精度和并行分塊

表 2 給出流場參數(shù)設(shè)置, 式(8)給出慣性顆粒在Taylor-Green 渦中的控制方程, 采用單向耦合的 PP-DNS 與理論解對比。關(guān)閉 OpenCFD 的可壓縮求解器, 全場定常施加一個不可壓縮的 Taylor-Green 流動, 從而避免構(gòu)造可壓縮流場的繁瑣。理論解如式(9)和(10)所示, 理論速度*采用 Taylor-Green 解析解, PP-DNS 使用的*從施加流場插值獲得。圖 6 表明, Lagrange 插值算法的精度誤差隨著階數(shù)的增加而減小, 隨著網(wǎng)格尺寸的減小而減小。

一個單向耦合慣性顆粒在 Taylor Green 渦中的運動如圖 7 所示?？梢钥闯? PP-DNS 求解器的動態(tài)鏈表數(shù)組通信模塊調(diào)試成功, 顆?？梢栽诓煌?MPI分塊之間準確地收發(fā)通信, 不需要單獨分配緩沖內(nèi)存來交換彌散顆粒相的數(shù)據(jù)。

圖5 可壓縮均勻流的DNS數(shù)值網(wǎng)格示意圖

表1 可壓縮均勻流的DNS流場參數(shù)設(shè)置

表2 單向耦合施加的不可壓縮Taylor-Green流場參數(shù)設(shè)置

圖7 動態(tài)鏈表數(shù)組并行分塊驗證

3.1.2考慮可壓縮效應(yīng)的相間阻力

對于均勻流中的顆粒非定常速度問題采用雙向耦合驗證, DNS 設(shè)置見表 1, 來流速度為*=0.6。顆粒在可壓縮流動中的統(tǒng)計定常阻力公式可以近似表示為滑移馬赫數(shù)和滑移雷諾數(shù)的函數(shù)[51–52](式(11)), 理論解可以表示為式(12)[8]。圖 8 表明, 本文的 PP-DNS 結(jié)果與半經(jīng)驗理論解(式(12))吻合。

(12)

3.1.3 歷史力

3.1.4 重力

重力作用往往在湍流脈動尺度上影響顆粒動力學(xué)[22–23], 有研究發(fā)現(xiàn)在槽道流中重力影響湍流統(tǒng)計矩[55–56]。式(14)為考慮 Stokes 阻力和一般重力場作用下的顆粒動力學(xué)方程, 式(15)為理論解[8]。圖 10表明, 本文 PP-DNS 結(jié)果與半經(jīng)驗理論解(式(15))相吻合。

3.1.5 熱泳力

由于極端力學(xué)條件下的壁流邊界層存在溫度梯度, 細小的顆粒會受到溫度梯度誘導(dǎo)的熱泳力[8]影響。式(16)為考慮 Stokes 阻力和熱泳力的方程, 式(17)為理論解。圖 11 顯示, 本文 PP-DNS 結(jié)果與理論解(式(17))吻合。

圖11 受到熱泳力的慣性顆粒在均勻流中的非定常速度與理論解的驗證

3.1.6 電場力

在極端力學(xué)環(huán)境下, 顆粒群間的碰撞會產(chǎn)生局部電場[4,58], 因此帶電顆粒動力學(xué)不可忽略。本文不考慮顆粒群相互碰撞產(chǎn)生的動態(tài)電場, 只考慮電場是常矢量的情況。式(18)為考慮 Stokes 阻力和電場力的方程, 式(19)為理論解[8]。圖 12 顯示, 本文PP-DNS 結(jié)果與理論解(式(19))吻合。

圖13 受到磁場力的慣性顆粒在均勻流中的非定常速度與理論解的驗證

3.1.7 磁場力

根據(jù)經(jīng)典電磁理論, 帶電顆粒群運動會感生動態(tài)電磁場。當前暫時不考慮顆粒群相互碰撞產(chǎn)生的動態(tài)電磁場, 只考慮磁場是常矢量的情況。式(20)為考慮 Stokes 阻力和磁場力的方程, 式(21)為理論解[8]。圖 13 顯示, 本文 PP-DNS 結(jié)果與理論解(式 (21))吻合。

表3 可壓縮槽道流的DNS流場參數(shù)設(shè)置

圖15 單相可壓縮槽道湍流Reτ=150的湍流一階矩和二階矩驗證

3.1.8 熱對流

圓球在可壓縮流動中的對流熱交換系數(shù)會被可壓縮效應(yīng)修正[3,8,52]。式(22)為考慮熱對流效應(yīng)的顆粒熱力學(xué)方程, 式(23)為理論解[8]。圖 14 顯示, 本文 PP-DNS 結(jié)果與理論解(式(23))吻合。

圖16 彌散顆粒相St+=25的一階和二階統(tǒng)計矩驗證

圖17 彌散顆粒相St+=5的一階和二階統(tǒng)計矩驗證

綜上所述, 單顆粒驗證問題本質(zhì)上分為兩個方面。1)部分模塊采用單向耦合。任何形式的 Bug 和錯誤都與流體求解器無關(guān), 因此可以很好地鎖定由編程錯誤帶來 Bug 的位置, 提高編程效率。2)相間阻力和熱對流模塊采用雙向耦合。這是因為根據(jù)理論方程, 盡管多物理效應(yīng)復(fù)雜, 但其反饋到攜帶流體相的路徑是唯一確定的[8,39], 即顆粒動力學(xué)方程只能通過相間阻力反饋到流體動量方程中, 多物理效應(yīng)只能通過改變相間阻力, 間接地調(diào)制流體動量。同理, 顆粒熱力學(xué)方程只能通過熱對流反饋到流體能量方程中, 熱輻射效應(yīng)只能通過改變熱對流, 間接地調(diào)制流體能量。從算例結(jié)果可知, 本文可壓縮顆粒兩相流求解器可以很好地數(shù)值模擬點力顆粒在極端力學(xué)環(huán)境中的動力學(xué)和熱力學(xué)過程, 并能實現(xiàn)兩相耦合。

圖18 彌散顆粒相St+=1的一階和二階統(tǒng)計矩驗證

3.2 顆粒群問題

3.2.1數(shù)值標模1

表4 可壓縮槽道流Ma=0.3, Reτ=180的DNS流場參數(shù)設(shè)置

圖19 單相可壓縮槽道湍流Reτ=180的湍流一階矩、二階矩驗證和湍流能量平衡驗證

圖20 槽道顆粒湍流Reτ=180, =0.50脈動速度驗證

3.2.2數(shù)值標模2

圖21 槽道顆粒湍流Reτ=180, =0.50湍動能驗證

圖22 槽道顆粒湍流Reτ=180, =0.75脈動速度驗證

圖23 槽道顆粒湍流Reτ=180, =0.75湍動能驗證

4 并行效率

以表 3 的參數(shù)設(shè)置為例, 在 Re=180, Ma=0.3,p=105的條件下, 考察不同并行核數(shù)和不同量級的顆粒數(shù)目對顆粒兩相流求解器并行綜合效率的影響。圖 24 對比顆粒兩相流求解器與單獨攜帶流體相求解器的并行運行時間和加速度比?？梢钥闯? 本文開發(fā)的顆粒兩相流求解器的絕對運行時間和加速比與流體求解器相差不大。圖 25 對比兩者的并行效率[39]和顆粒求解器的計算增量?？梢钥闯? 兩者的并行效率基本上相同, 且顆粒兩相流求解器導(dǎo)致的計算增量為 30%。

以表 3 中算例的顆粒數(shù)目p=105為基準, 考察Re=180, Ma=0.3 條件下, 不同顆粒數(shù)目p=105, 107,108(“超載”)和不同并行核數(shù)對并行綜合效率的影響, 結(jié)果如圖 26 所示。可以看出, 顆粒兩相流求解器的運行時間隨著運行核數(shù)的增加而線性下降, 相對于標準算例p=105的運行時間, “超載”算例的運行時間隨著顆粒數(shù)據(jù)增加而線性增加。顆粒求解器并行效率的線性特性體現(xiàn)了動態(tài)鏈表數(shù)組在求解顆粒數(shù)目巨大問題上的并行優(yōu)勢。與全局數(shù)組相比, 顆粒數(shù)目增加帶來的并行效率是非線性降低的, 甚至?xí)驗檎加玫膬?nèi)存太多而導(dǎo)致內(nèi)存溢出, 使得算例無法進行。

圖24 可壓縮槽道顆粒湍流的并行計算時間和加速比與分塊數(shù)目的關(guān)系

圖25 可壓縮槽道顆粒湍流的并行計算效率η和顆粒求解器的計算增量與分塊數(shù)目的關(guān)系

圖26 可壓縮槽道顆粒湍流的并行計算時間和加速比與分塊數(shù)目的關(guān)系

5 結(jié)論和展望

本研究開發(fā)了一種基于點力模型的含顆?？蓧嚎s流直接數(shù)值模擬求解器, 并與當前工業(yè)需求密切相關(guān)的多物理效應(yīng)模塊進行驗證。采用動態(tài)鏈表數(shù)組對彌散相顆粒群進行內(nèi)存分配, 使得每個 MPI 并行分塊使用的內(nèi)存是全局數(shù)組分配方式的 1/。并行分塊數(shù)目越大, 基于動態(tài)鏈表數(shù)組的顆粒求解器相對于全局數(shù)組的內(nèi)存利用效率越高。理論分析表明, 對于顆粒碰撞問題, 動態(tài)鏈表數(shù)組的顆粒求解器的計算量是全局數(shù)組的 1/2。本文檢驗了不同MPI 分塊間的動態(tài)鏈表數(shù)組通信, 并對多物理效應(yīng)的可壓縮顆粒兩相流求解器進行驗證。在流動馬赫數(shù)漸進趨于零的情況下, 對不可壓顆粒湍流標模進行湍流各階統(tǒng)計矩的驗證。顆粒并行求解器采用動態(tài)鏈表技術(shù), 對彌散顆粒相進行分塊內(nèi)存分配以及邏輯層次封裝技術(shù), 提高了計算效率并節(jié)約了內(nèi)存, 使得顆粒并行求解器具備向 GPU 版本拓展的可 能性。

可壓縮顆粒解析的直接數(shù)值模型算法涉及大量拉格朗日點的信息反饋到局部歐拉點[16], 導(dǎo)致局部的歐拉點“接收”到大量不同來源的數(shù)值間斷, 從而增加可壓縮顆粒兩相流求解器的數(shù)值剛性。因此, 在未來的研究中, 一方面需要增加壓力修正的數(shù)值算法模塊來穩(wěn)定可壓縮顆粒兩相流求解器, 另一方面, 需要增加多重網(wǎng)格算法對求解器進行加速來解決由于數(shù)值剛性帶來的計算量大的問題。

[1] Crowe C T, Schwarzkopf J D, Sommerfield M, et al. Multiphase flows with droplets and particles. Multi-phase flows with droplets and particles. Boca Raton: CRC Press, 2011

[2] 岳朋濤. 超燃沖壓發(fā)動機燃燒室若干問題的研究[D]. 北京: 中國科學(xué)技術(shù)大學(xué), 2002

[3] 李婷婷, 劉朋欣, 袁先旭, 等. 高焓可壓縮流動熱輻射顆粒兩相流理論分析. 北京大學(xué)學(xué)報(自然科學(xué)版), 2022, 58(6): 989–998

[4] 胡文文. 風沙流中沙粒帶電機理及荷質(zhì)比研究[D]. 蘭州: 蘭州大學(xué), 2012

[5] Zheng Xiaojing, He Lihong, Zhou Youhe. Theoretical model of the electric field produced by charged parti-cles in windblown sand flux. Journal of Geophysical Research-Atmospheres, 2004, 109: D15208

[6] Zamansky R, Coletti F, Massot M, et al. Radiation induces turbulence in particle-laden fluids. Physics of Fluids, 2014, 26(7): 071701

[7] Healy D P, Young J B. An experimental and theoretical study of particle deposition due to thermophoresis and turbulence in an annular flow. International Journal of Multiphase Flow, 2010, 36(11/12): 870–881

[8] 李青, 涂國華, 李婷婷, 等. 高焓流動中的可壓縮顆粒求解器(第 1 部分): 考慮多物理效應(yīng)的點力顆粒兩相流理論方程. 空氣動力學(xué)學(xué)報, 2023, 41(8): 71–86

[9] Li Xinliang, Fu Dexun, Ma Yanwen. Direct numerical simulation of a spatially evolving supersonic turbulent boundary layer at Ma=6. Chinese Physics Letters, 2006, 23(6): 1519–1522

[10] Hirsch C. Numerical computation of internal and ex-ternal flows. New York: John Wiley & Sons Inc, 1988

[11] 應(yīng)隆安, 滕振寰. 雙曲型守恒律方程及其差分方法. 北京: 科學(xué)出版社, 1991

[12] 谷超豪. 數(shù)學(xué)物理方程. 第3版. 北京: 高等教育出版社, 2012

[13] Pierce C D. Progress-variable approach for large-eddy simulation of turbulent combustion [D]. Stanford: Stan-ford University, 2001

[14] Yu Zhaosheng, Shao Xueming. A direct-forcing fic-titious domain method for particulate flows. Journal of Computational Physics, 2007, 227(1): 292–314

[15] Mehrabadi M, Horwitz J, Subramaniam S, et al. A direct comparison of particle-resolved and point-particle methods in decaying turbulence. Journal of Fluid Mechanics, 2018, 850: 336–369

[16] Yu Zhaosheng, Xia Yan, Guo Yu, et al. Modulation of turbulence intensity by heavy finite-size particles in upward channel flow. Journal of Fluid Mechanics, 2021, 913: A3

[17] 李瑞元, 陳飛國, 葛蔚, 等. 高馬赫數(shù)低雷諾數(shù)條件下圓球繞流曳力系數(shù). 空氣動力學(xué)學(xué)報, 2021, 39(3): 201–208

[18] Maxey M R. Equation of motion for a small rigid sphere in a nonuniform flow. Physics of Fluids, 1983, 26(4): 883–889

[19] Balachandar S, Eaton J K. Turbulent dispersed multi-phase flow. Annual Review of Fluid Mechanics, 2010, 42: 111–133

[20] 袁先旭, 陳堅強, 杜雁霞, 等. 國家數(shù)值風洞(NNW)工程中的CFD基礎(chǔ)科學(xué)問題研究進展. 航空學(xué)報, 2021, 42(9): 625733–625733

[21] 陳堅強. 國家數(shù)值風洞(NNW)工程關(guān)鍵技術(shù)研究進展. 中國科學(xué): 技術(shù)科學(xué), 2021, 51(11): 1326–1347

[22] Wang Lianping, Maxey M. Settling velocity and con-centration distribution of heavy particles in homoge-neous isotropic turbulence. J Fluid Mech, 1993, 256 (1): 27–68

[23] Rosa B, Parishani H, Ayala O, et al. Settling velocity of small inertial particles in homogeneous isotropic turbulence from high-resolution DNS. International Journal of Multiphase Flow, 2016, 83: 217–231

[24] Ferrante A, Elghobashi S. Reynolds number effect on drag reduction in a microbubble-laden spatially deve-loping turbulent boundary layer. Journal of Fluid Me-chanics, 2005, 543(1): 93–106

[25] Lain S, Sommerfield M. Euler/Lagrange computations of pneumatic conveying in a horizontal channel with different wall roughness. Powder Technology, 2008, 184(1): 76–88

[26] Marchioli C, Soldati A. Mechanisms for particle trans-fer and segregation in a turbulent boundary layer. Jour-nal of Fluid Mechanics, 2002, 468: 283–315

[27] Daitche A. On the role of the history force for inertial particles in turbulence. Journal of Fluid Mechanics, 2015, 782: 567–593

[28] Zamansky R, Coletti F, Massot M, et al. Turbulent thermal convection driven by heated inertial particles. Journal of Fluid Mechanics, 2016, 809: 390–437

[29] Zhao Lihao, Andersson H I, Gillissen J. Interphasial energy transfer and particle dissipation in particle-laden wall turbulence. Journal of Fluid Mechanics, 2013, 715(1): 32–59

[30] Zhou Tian, Zhao Lihao, Huang Weixi, et al. Non-monotonic effect of mass loading on turbulence modulations in particle-laden channel flow. Physics of Fluids, 2020, 32(4): 043304

[31] Pan Qingqing, Xiang Hong, Wang Ze, et al. Kinetic energy balance in turbulent particle-laden channel flow. Physics of Fluids, 2020, 32(7): 073307

[32] Li Dong, Luo Kun, Fan Jianren. Modulation of tur-bulence by dispersed solid particles in a spatially de-veloping flat-plate boundary layer. Journal of Fluid Mechanics, 2016, 802: 359–394

[33] Dai Qi, Luo Kun, Jin Tai, et al. Direct numerical simulation of turbulence modulation by particles in compressible isotropic turbulence. Journal of Fluid Mechanics, 2017, 832: 438–482

[34] Xiao Wei, JinTai, Luo Kun, et al. Eulerian-Lagrangian direct numerical simulation of preferential accumu-lation of inertial particles in a compressible turbulent boundary layer. Journal of Fluid Mechanics, 2020, 903: A19

[35] Wang Lianping, Wexler A S, Zhou Y. Statistical me-chanical description and modelling of turbulent colli-sion of inertial particles. Journal of Fluid Mechanics, 2000, 415: 117–153

[36] Fevrier P, Simonin O, Squires K D. Partitioning of particle velocities in gas-solid turbulent flows into a continuous field and a spatially uncorrelated random distribution: theoretical formalism and numerical stu-dy. Journal of Fluid Mechanics, 2005, 533: 1–46

[37] Yu Wenchao, Fede P, Yazdanpanah M, et al. Gas-solid fluidized bed simulations using the filtered approach: validation against pilot-scale experiments. Chemical Engineering Science, 2020, 217: 115472

[38] Marchioli C. Large-eddy simulation of turbulent dis-persed flows: a review of modelling approaches. Acta Mechanica, 2017, 228(3): 741–771

[39] Tian Baolin, Zeng Junsheng, MengBaoqing, et al. Com-pressible multiphase particle-in-cell method (CMP-PIC) for full pattern flows of gas-particle system. Journal of Computational Physics, 2020, 418: 109602

[40] Martin J E, Meiburg E. The accumulation and disper-sion of heavy particles in forced two-dimensional mixing layers. I. The fundamental and subharmonic cases. Physics of Fluids, 1994, 6(3): 1116–1132

[41] 李新亮, 馬延文, 傅德薰. 可壓槽道湍流的直接數(shù)值模擬及標度律分析. 中國科學(xué): A輯, 2001, 31(2): 153–164

[42] 童福林, 周桂宇, 周浩, 等. 激波/湍流邊界層干擾物面剪切應(yīng)力統(tǒng)計特性. 航空學(xué)報, 2019, 40(5): 122504–122504

[43] Mirjalili S, Mani A. Consistent, energy-conserving momentum transport for simulations of two-phase flows using the phase field equations. Journal of Computational Physics, 2019, 426: 109918

[44] Zhang Bo, Boyd B, Ling Y. Direct numerical simu-lation of compressible interfacial multiphase flows using a mass-momentum-energy consistent volume-of-fluid method. Computers & Fluids, 2022, 236: 105267

[45] Patankar S V. Numerical heat transfer and fluid flow. Boca Raton: CRC Press, 1980

[46] Kim J, Moin P. Turbulence statistics in fully developed channel flow at low Reynolds number. J Fluid Mech, 1987, 177: 133–166

[47] Lawrence C E. Partial differential equations. Rhode Island: American Mathematical Society, 1998

[48] Lele S K. Compact finite difference schemes with spectral-like resolution. Journal of Computational Phy-sics, 1992, 103(1) : 16–42

[49] Burden R L, Faires J D. Numerical analysis. 9th ed. Boston: Richard Stratton, 2011

[50] Bassenne M, Fu Lin, Mani A. Time-accurate and highly-stable explicit operators for stiff differential equations. Journal of Computational Physics, 2021, 424: 109847

[51] Loth E. Compressibility and rarefaction effects on drag of a spherical particle. AIAA Journal, 2008, 46(9): 2219–2228

[52] 尚慶, 沈清, 莊逢甘. 超聲速氣流中簡化微小液滴橫向噴射的汽化過程探索. 空氣動力學(xué)學(xué)報, 2011, 29(1): 1–9

[53] Basset A B. A treatise on Hydrodynamics. Nature, 1888, 40: 412–423

[54] Mei Renwei, Adrian R. Flow past a sphere with an oscillation in the free-stream velocity and unsteady drag at finite Reynolds number. Journal of Fluid Me-chanics, 1992, 237(1): 323–341

[55] Righetti M, Romano G P. Particle–fluid interactions in a plane near-wall turbulent flow. Journal of Fluid Mechanics, 2004, 505: 93–121

[56] Hout R V. Time-resolved PIV measurements of the in-teraction of polystyrene beads with near-wall-coherent structures in a turbulent channel flow. International Journal of Multiphase Flow, 2011, 37(4): 346–357

[57] 董雙嶺, 曹炳陽, 過增元. 作用在粒子上的熱泳升力研究. 工程熱物理學(xué)報, 2015, 36(5): 1063–1066

[58] Zheng Xiaojing. Electrification of wind-blown sand: recent advances and key issues. European Physical Journal E Soft Matter, 2013, 36(12): no. 138

[59] Marchioli C, Soldati A, Kuerten J, et al. Statistics of particle dispersion in direct numerical simulations of wall-bounded turbulence: results of an international collaborative benchmark test. International Journal of Multiphase Flow, 2008, 34,(9): 879–893

[60] Vreman A W, Kuerten J. Comparison of direct nu-merical simulation databases of turbulent channel flow at Re=180. Physics of Fluids, 2014, 26(1): 015102

[61] Hoyas S, Jiménez J. Reynolds number effects on the Reynolds-stress budgets in turbulent channels. Physics of Fluids, 2008, 20(10): 101511

MPI Solver of Particle-Laden Compressible High Enthalpy Flow: Numerical Method and Validation

LI Qing1,2, YU Zhaosheng3, LIU Pengxin1, LI Tingting4, CHEN Jianqiang1, YUAN Xianxu1,?

1. State Key Laboratory of Aerodynamics, China Aerodynamics Research and Development Center, Mianyang 621000; 2. Tianmushan Laboratory, Hangzhou 310000; 3. School of Aeronautics and Astronautics, Zhejiang University, Hangzhou 310027; 4. School of Chemical Engineering and Technology, Xi’an Jiaotong University, Xi’an 710049; ? Corresponding author, E-mail: yuanxianxu@cardc.cn

Based on the theoretical equations of compressible particle-laden flow in multi-physical conditions, a MPI-Particle Solver based on dynamic linked list array is developed, which is coupled with a compressible flow solver of carried phase. The dispersed particle phase is different from the carried fluid phase, the previous one is based on the Lagrange coordinate system, while the later one is based on the Euler coordinate system. The traditional way is to use the global array to assign the memory to the dispersed particle phase to achieve the transfer between two-phase, but at expense of low computational efficiency. This paper uses dynamic linked list array to allocate memory to dispersed particle phase, achieve both problems. A series of DNS validations with references case have been done, in terms of multi-physical effects and two canonical cases at incompressible limit.

dynamic linked list array; Lagrange/Euler coordinate system conversion; particle solver of particle-laden compressible flow