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

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

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

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

10.13209/j.0479-8023.2023.022

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

2022–12–12;

2023–01–29

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

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

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

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

動態鏈表數組; 拉格朗日/歐拉坐標系轉換; 可壓縮點力顆粒兩相流直接數值模擬求解器

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

1 研究背景

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

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

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

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

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

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

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

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

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

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

2 數值方法

2.1 流體求解器

本文在極端力學環境下可壓縮點力顆粒兩相流理論方程的基礎上, 開發一種可壓縮流動條件下基于動態鏈表數組的 PP-DNS 求解器, 使其能夠與可壓縮流體求解器耦合。

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

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

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

2.2 顆粒求解器

將顆粒動力學以及熱力學方程寫入流體求解器的 Runge Kutta 三階時間步推進循環(式(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 三階時間步推進的子循環步, 先后取=1, 2, 3; conservation 項表示可壓縮流體力學數值離散方程的守恒項; source 項表示顆粒相對流體相的反饋。

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

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

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

3 標模驗證算例

3.1 單顆粒問題

3.1.1插值精度和并行分塊

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

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

圖5 可壓縮均勻流的DNS數值網格示意圖

表1 可壓縮均勻流的DNS流場參數設置

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

圖7 動態鏈表數組并行分塊驗證

3.1.2考慮可壓縮效應的相間阻力

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

(12)

3.1.3 歷史力

3.1.4 重力

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

3.1.5 熱泳力

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

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

3.1.6 電場力

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

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

3.1.7 磁場力

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

表3 可壓縮槽道流的DNS流場參數設置

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

3.1.8 熱對流

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

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

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

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

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

3.2 顆粒群問題

3.2.1數值標模1

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

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

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

3.2.2數值標模2

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

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

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

4 并行效率

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

以表 3 中算例的顆粒數目p=105為基準, 考察Re=180, Ma=0.3 條件下, 不同顆粒數目p=105, 107,108(“超載”)和不同并行核數對并行綜合效率的影響, 結果如圖 26 所示?？梢钥闯? 顆粒兩相流求解器的運行時間隨著運行核數的增加而線性下降, 相對于標準算例p=105的運行時間, “超載”算例的運行時間隨著顆粒數據增加而線性增加。顆粒求解器并行效率的線性特性體現了動態鏈表數組在求解顆粒數目巨大問題上的并行優勢。與全局數組相比, 顆粒數目增加帶來的并行效率是非線性降低的, 甚至會因為占用的內存太多而導致內存溢出, 使得算例無法進行。

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

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

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

5 結論和展望

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

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

[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] 岳朋濤. 超燃沖壓發動機燃燒室若干問題的研究[D]. 北京: 中國科學技術大學, 2002

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

[4] 胡文文. 風沙流中沙粒帶電機理及荷質比研究[D]. 蘭州: 蘭州大學, 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 部分): 考慮多物理效應的點力顆粒兩相流理論方程. 空氣動力學學報, 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] 應隆安, 滕振寰. 雙曲型守恒律方程及其差分方法. 北京: 科學出版社, 1991

[12] 谷超豪. 數學物理方程. 第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] 李瑞元, 陳飛國, 葛蔚, 等. 高馬赫數低雷諾數條件下圓球繞流曳力系數. 空氣動力學學報, 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] 袁先旭, 陳堅強, 杜雁霞, 等. 國家數值風洞(NNW)工程中的CFD基礎科學問題研究進展. 航空學報, 2021, 42(9): 625733–625733

[21] 陳堅強. 國家數值風洞(NNW)工程關鍵技術研究進展. 中國科學: 技術科學, 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] 李新亮, 馬延文, 傅德薰. 可壓槽道湍流的直接數值模擬及標度律分析. 中國科學: A輯, 2001, 31(2): 153–164

[42] 童福林, 周桂宇, 周浩, 等. 激波/湍流邊界層干擾物面剪切應力統計特性. 航空學報, 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] 尚慶, 沈清, 莊逢甘. 超聲速氣流中簡化微小液滴橫向噴射的汽化過程探索. 空氣動力學學報, 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] 董雙嶺, 曹炳陽, 過增元. 作用在粒子上的熱泳升力研究. 工程熱物理學報, 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