2015工程熱物理氣動(dòng)文獻(xiàn)蘇欣榮

1、基金項(xiàng)目：國(guó)家自然科學(xué)基金重點(diǎn)項(xiàng)目（51136003）中國(guó)工程熱物理學(xué)會(huì) 熱機(jī)氣動(dòng)熱力學(xué)學(xué)術(shù)會(huì)議論文 編號(hào)：152100并行計(jì)算環(huán)境中的高效氣熱耦合方法蘇欣榮，袁新清華大學(xué)熱科學(xué)與動(dòng)力工程教育部重點(diǎn)實(shí)驗(yàn)室，北京100084Tel: Email: 摘 要 氣熱耦合計(jì)算有助于高溫透平葉片熱流和溫度分布的準(zhǔn)確預(yù)估。在氣熱耦合計(jì)算中，流體區(qū)域和固體區(qū)域可能采用不同的網(wǎng)格形式和不同的求解器；此外，在DES/LES等大規(guī)模計(jì)算中，流體區(qū)域和固體區(qū)域的耦合需要在并行環(huán)境中進(jìn)行。因此，需要發(fā)展在不同求解區(qū)域支持不同的網(wǎng)格形式和求解器，并且適用于大規(guī)模并行計(jì)算環(huán)境的高效氣熱耦合方法。

2、本文發(fā)展了并行計(jì)算環(huán)境中的高效氣熱耦合方法，其流體域和固體域之間的耦合通過高效的并行數(shù)據(jù)交換來實(shí)現(xiàn)，其對(duì)求解器類型和網(wǎng)格形式以及負(fù)載分配沒有限制。數(shù)值算例驗(yàn)證了本文氣熱耦合方法在并行計(jì)算環(huán)境中的高效性和靈活性。關(guān)鍵詞 氣熱耦合；并行計(jì)算；數(shù)據(jù)交換前言燃?xì)廨啓C(jī)和航空發(fā)動(dòng)機(jī)不斷向著更高的效率發(fā)展，與之對(duì)應(yīng)，透平進(jìn)口溫度不斷提高且已經(jīng)遠(yuǎn)遠(yuǎn)超過材料所能承受的溫度極限。除了材料技術(shù)的不斷進(jìn)步外，還必須依賴于先進(jìn)的冷卻技術(shù)以保證高溫透平葉片的正常運(yùn)行和工作壽命。對(duì)于高溫透平葉片，如果實(shí)際工作溫度高于許可溫度10，其預(yù)測(cè)壽命可能減半1。因此，在透平葉片設(shè)計(jì)中，透平葉片熱流和溫度分布的準(zhǔn)確預(yù)估是必要前提之一。

3、同時(shí)考慮流動(dòng)和固體傳熱的氣熱耦合數(shù)值方法是可行的研究手段之一。Hylton等2首先使用氣熱耦合方法研究了包含冷卻孔的MARK-，C3X葉片，在其研究中，流體域使用基于多塊結(jié)構(gòu)化網(wǎng)格的有限體積程序，固體域使用邊界元方法進(jìn)行求解。Bohn等3在氣熱耦合領(lǐng)域開展了大量研究，在其計(jì)算中，流體域和固體域均使用多塊結(jié)構(gòu)化網(wǎng)格。Duchaine等4發(fā)展了適用于大渦模擬的氣熱耦合方法，在其研究中，流體域和固體域均使用基于非結(jié)構(gòu)網(wǎng)格的三階Taylor-Galerkin方法。周馳等5發(fā)展了基于區(qū)域分解算法的氣熱耦合方法，在其研究中，流體域和固體域均采用基于結(jié)構(gòu)化網(wǎng)格的有限差分方法。郝增榮等6發(fā)展了基于間斷有限元方

4、法的氣熱耦合方法，在其研究中流體域和固體域均采用基于非結(jié)構(gòu)網(wǎng)格的二階間斷有限元方法。總結(jié)已有氣熱耦合研究，從求解器的角度來看可以分為兩類。第一類在固體域中使用的求解器通過改動(dòng)流體求解器得到，即在固體域中求解密度為常數(shù)、速度為0的流體控制方程，這類氣熱耦合方法可以直接使用已有流體求解器，因此簡(jiǎn)便易行，同時(shí)由于流體域和固體域求解器采用相同的數(shù)值算法和數(shù)據(jù)結(jié)構(gòu)，因此易于數(shù)據(jù)交換。但在此類方法中固體域計(jì)算量較大，流體域和固體域的求解器緊密耦合，同時(shí)也限制了其靈活性。在第二類氣熱耦合方法中，固體域求解熱傳導(dǎo)方程，流體域和固體域的耦合通過界面上的數(shù)據(jù)傳遞來實(shí)現(xiàn)。此類方法具有更好的靈活性，因此在已有研究中大

5、量采用。但在目前多數(shù)研究中，流體域和固體域求解器均基于相同的網(wǎng)格類型和底層數(shù)據(jù)結(jié)構(gòu)，因此在一定程度上仍限制了其靈活性和對(duì)大規(guī)模復(fù)雜幾何結(jié)構(gòu)氣熱耦合問題的適應(yīng)性。例如，對(duì)于透平內(nèi)流研究，主流通道內(nèi)適宜使用多塊結(jié)構(gòu)化網(wǎng)格以達(dá)到較高的數(shù)值精度，而在盤腔、二次空氣系統(tǒng)、冷卻通道等復(fù)雜幾何流體域和固體域則適宜使用非結(jié)構(gòu)網(wǎng)格。在基于DES、LES等高精度湍流模擬的大規(guī)模計(jì)算中，流體域和固體域的計(jì)算網(wǎng)格通常分解至多個(gè)計(jì)算進(jìn)程中，此時(shí)流體域和固體域之間的數(shù)據(jù)交換需要在并行環(huán)境中高效進(jìn)行，且求解域之間的耦合盡可能不對(duì)流體求解器和固體求解器的負(fù)載平衡產(chǎn)生影響。本文發(fā)展了在并行計(jì)算環(huán)境中的高效氣熱耦合方法，其中流體

6、域和固體域的耦合通過高效的并行數(shù)據(jù)交換來實(shí)現(xiàn)，該方法對(duì)流體域和固體域的求解器和計(jì)算網(wǎng)格類型沒有限制，同時(shí)其對(duì)流體域和固體域的負(fù)載分配也沒有限制，因此具有高效、靈活的特點(diǎn)。數(shù)值方法控制方程對(duì)于流體域，考慮積分形式的三維Navier-Stokes方程 MACROBUTTON AuroraSupport.EditInitialCounterValues ADDIN MACROBUTTON AuroraSupport.NoMacro ADDIN (1)其中守恒變量(2)通量F, G和H分別表征對(duì)流通量；Fv，Gv和Hv分別表征粘性通量。在本文計(jì)算中，均使用理想氣體假設(shè) (3)，其中R為理想氣體常數(shù)。在

7、本文計(jì)算中，分子粘性系數(shù)通過Sutherland公式計(jì)算得到；湍流粘性系數(shù)通過求解湍流控制方程給出。對(duì)于固體域，控制方程為： (4)，其中，cp，和S分別表征密度，比熱，熱傳導(dǎo)系數(shù)和內(nèi)熱源強(qiáng)度?；隈詈掀鞯臍鉄狁詈嫌?jì)算模型對(duì)于式（1）中給出的流體域內(nèi)的Navier-Stokes方程，可以使用有限差分方法、有限體積方法或者間斷有限元等方法求解。與之對(duì)應(yīng)，所使用的網(wǎng)格形式可以是多塊結(jié)構(gòu)化網(wǎng)格，也可以是非結(jié)構(gòu)化網(wǎng)格。對(duì)于式（4）中給出的固體域中的熱傳導(dǎo)方程，通常使用基于非結(jié)構(gòu)網(wǎng)格的有限體積方法或有限元方法求解。同時(shí)，對(duì)于大規(guī)模計(jì)算，如基于DES、LES等的透平復(fù)雜結(jié)構(gòu)氣熱耦合計(jì)算，通常需要借助于大規(guī)

8、模并行計(jì)算機(jī)集群，因此流體域和固體域之間的耦合需要在并行環(huán)境中進(jìn)行且流體求解器和固體求解器的負(fù)載平衡特性需要盡可能被保持?；谝陨峡紤]，對(duì)氣熱耦合方法有如下幾點(diǎn)要求：氣熱耦合方法對(duì)流體域和固體域的求解器類型和網(wǎng)格形式?jīng)]有限制；氣熱耦合方法對(duì)流體域和固體域的負(fù)載平衡特性，即流體和固體計(jì)算網(wǎng)格在計(jì)算節(jié)點(diǎn)上的分布沒有限制。上面第一條要求決定了需要在流體域和固體域使用不同的求解器，且求解器之間的耦合通過數(shù)據(jù)交換來實(shí)現(xiàn)。對(duì)于不同類型的求解器，其最優(yōu)的數(shù)據(jù)交換方式也不同。對(duì)于基于網(wǎng)格點(diǎn)的求解器，如有限差分方法或基于格點(diǎn)的有限體積方法，如圖1所示，其最優(yōu)的數(shù)據(jù)交換方式是流體域向固體域傳遞熱流，固體域向流體域

9、傳遞溫度。如圖2所示，如果求解器基于格心，即使用基于結(jié)構(gòu)網(wǎng)格的有限體積方法或格心型非結(jié)構(gòu)有限體積方法，最佳的耦合方法是流體和固體均向壁面?zhèn)鬟f溫度，根據(jù)熱流連續(xù)條件確定壁面溫度。圖1 使用格點(diǎn)型求解器的氣熱耦合計(jì)算中最佳數(shù)據(jù)交換方向示意圖圖2 使用格心型求解器的氣熱耦合計(jì)算中最佳數(shù)據(jù)交換方向示意圖總結(jié)圖1和圖2中不同的數(shù)據(jù)交換方式，可以得到同時(shí)適用于不同類型求解器和網(wǎng)格形式的具有普適性的氣熱耦合方法。如圖3所示，在該耦合方法中，除流體域求解器和固體域固體域求解器外，還包括一個(gè)額外的組成部分：負(fù)責(zé)數(shù)據(jù)交換的耦合器。在此流熱耦合計(jì)算模型中，首先流體域求解器和固體域求解器分別向耦合器發(fā)送數(shù)據(jù)；其次耦合

10、器對(duì)接收到的數(shù)據(jù)進(jìn)行處理；最后流體域求解器和固體域求解器分別從耦合器接收數(shù)據(jù)。借助于耦合器，流體域求解器和固體域求解器之間不再直接進(jìn)行數(shù)據(jù)交換，因此極大的降低了對(duì)兩個(gè)計(jì)算域中求解器和網(wǎng)格類型的限制。圖3 基于流體域求解器、耦合器和固體域求解器的氣熱耦合模型示意圖耦合器的并行實(shí)現(xiàn)圖3中給出的氣熱耦合模型為串行程序，即流體求解器、耦合器和固體求解器均在單個(gè)計(jì)算進(jìn)程中運(yùn)行。透平真實(shí)復(fù)雜幾何計(jì)算，尤其是基于DES、LES等高精度湍流模擬方法的大規(guī)模計(jì)算，通常需要借助于大規(guī)模并行計(jì)算來實(shí)現(xiàn)。在此類計(jì)算中，流體域和固體域通常使用各自的并行負(fù)載平衡策略把計(jì)算載荷盡可能均勻的分配到各個(gè)計(jì)算進(jìn)程中。由于流體域和

11、固體域的計(jì)算網(wǎng)格數(shù)目和計(jì)算量通常差異較大，分配到每個(gè)計(jì)算進(jìn)程上的流體域和固體域通常是不完全匹配的，也即在大規(guī)模并行計(jì)算環(huán)境中，在耦合器內(nèi)部進(jìn)行的數(shù)據(jù)插值和傳遞等也需要并行的進(jìn)行。圖4 基于流體域求解器、耦合器和固體域求解器的并行氣熱耦合模型示意圖本文發(fā)展的基于耦合器數(shù)據(jù)交換的氣熱耦合方法在并行計(jì)算環(huán)境中的示意圖如圖4所示。除在每個(gè)CPU計(jì)算進(jìn)程內(nèi)部的數(shù)據(jù)交換外，不同CPU進(jìn)程上的流體求解器之間需要進(jìn)行并行數(shù)據(jù)交換；同樣的在不同CPU進(jìn)程上的固體求解器之間也需要進(jìn)行并行數(shù)據(jù)交換；此外，根據(jù)流體域和固體域計(jì)算網(wǎng)格的實(shí)際并行分配情況，不同CPU進(jìn)程上的耦合器之間同樣需要進(jìn)行并行數(shù)據(jù)交換。在本文發(fā)展的

12、并行氣熱耦合數(shù)值工具中，這些并行數(shù)據(jù)交換均基于Message Passing Interface（MPI）標(biāo)準(zhǔn)函數(shù)庫(kù)進(jìn)行。通過對(duì)這些數(shù)據(jù)交換的先后順序進(jìn)行不同的組合，即可實(shí)現(xiàn)流熱強(qiáng)耦合、流熱弱耦合等不同的耦合計(jì)算方法。由圖3和圖4可知，在氣熱耦合計(jì)算的前處理階段，耦合器需要獲取固體壁面上流體側(cè)和固體側(cè)的計(jì)算網(wǎng)格。然后對(duì)每個(gè)流體側(cè)網(wǎng)格單元，耦合器搜尋固體側(cè)對(duì)應(yīng)的插值模板點(diǎn)以便在求解過程中進(jìn)行插值以獲取流體側(cè)所需要的數(shù)據(jù)。同樣對(duì)固體側(cè)網(wǎng)格單元需要執(zhí)行相似的操作。為了加速插值模板點(diǎn)的搜索過程，在本文所發(fā)展的并行耦合器中使用了基于二叉樹的空間遞歸劃分，相對(duì)于遍歷型的搜索插值模板點(diǎn)方法，使用該方法可以得

13、到極大的加速。計(jì)算結(jié)果計(jì)算對(duì)象使用本文所發(fā)展的并行氣熱耦合數(shù)值系統(tǒng)對(duì)某跨音透平進(jìn)行了數(shù)值模擬。在計(jì)算中，流體域使用本課題組獨(dú)立發(fā)展的高精度流體求解器7，湍流模型為Spalart-Allmaras一方程模型8，離散方程使用隱式時(shí)間推進(jìn)方法求解并配合多重網(wǎng)格方法以加速收斂；固體域中使用本課題組獨(dú)立發(fā)展的基于非結(jié)構(gòu)網(wǎng)格的有限體積求解器，離散方程的求解同樣使用高效隱式時(shí)間推進(jìn)方法。在本算例中，流體域進(jìn)口給定總溫、總壓和氣流角分布；流體域出口給定背壓。固體域的所有邊界均為氣熱耦合邊界。在本算例中使用了4個(gè)CPU進(jìn)程進(jìn)行并行計(jì)算，在該算例中的流體域和固體域計(jì)算網(wǎng)格并行負(fù)載分配情況如圖5所示。在計(jì)算的前處理

14、過程中分別對(duì)流體域和固體域進(jìn)行負(fù)載平衡，因此在每個(gè)CPU進(jìn)程中，固體壁面上的流體側(cè)網(wǎng)格單元和固體側(cè)網(wǎng)格單元并不完全匹配，插值模板點(diǎn)的搜尋過程同樣需要并行的完成。圖5 流體域和固體域計(jì)算網(wǎng)格并行負(fù)載分配示意圖計(jì)算結(jié)果及討論圖6中給出了本算例中流體域和固體域的收斂歷史。從圖6可以看出：使用本文所圖6 氣熱耦合算例中流體域和固體域的收斂歷史發(fā)展的并行氣熱耦合方法，經(jīng)過約30000次迭代后流體域和固體域中的計(jì)算均可穩(wěn)定的收斂至機(jī)器零。此外需要說明的是：在已有部分氣熱耦合研究中，需要在氣熱耦合界面上對(duì)熱流密度和溫度進(jìn)行亞松弛以保證穩(wěn)定性。在本文算例中，不需要使用亞松弛系數(shù)即可保證計(jì)算的穩(wěn)定收斂，即本文氣

15、熱耦合方法具有良好的穩(wěn)定性。從圖6中的收斂特性可知：當(dāng)計(jì)算推進(jìn)約5000步后，流體域和固體域的殘差均已下降至少6個(gè)數(shù)量級(jí)，已經(jīng)達(dá)到工程計(jì)算許可的收斂標(biāo)準(zhǔn)，說明本文氣熱耦合方法具有良好的收斂性。圖7 氣熱耦合算例中流體域馬赫數(shù)分布圖圖8 氣熱耦合算例中固體壁面上的溫度分布圖9 氣熱耦合算例中固體壁面上的熱流分布本算例中流體域計(jì)算結(jié)果如圖7所示。氣熱耦合固體壁面上的溫度分布如圖8所示，可見固體域內(nèi)的實(shí)際溫差很??；固體壁面上的熱流密度分布如圖9所示（由流體域向固體域傳遞的熱量定義為正），可見雖然實(shí)際溫差很小，但壁面上的熱流密度不可忽略，尤其是在尾緣附近熱流密度達(dá)到最大值且變化劇烈，這些現(xiàn)象在常用的絕

16、熱壁面計(jì)算是無法考慮的。結(jié) 論本文發(fā)展了并行計(jì)算環(huán)境中基于數(shù)據(jù)交換的高效氣熱耦合計(jì)算方法。該耦合方法中包含流體域求解器、固體域求解器和耦合器等三個(gè)組成部分。流體域和固體域之間的耦合以及數(shù)據(jù)傳遞通過并行耦合器實(shí)現(xiàn)。本文發(fā)展的并行氣熱耦合計(jì)算方法對(duì)流體域和固體域所分別使用的求解器和計(jì)算網(wǎng)格類型沒有限制，同時(shí)對(duì)并行環(huán)境中流體域和固體域的負(fù)載分配也沒有限制，因而有高效、靈活的特點(diǎn)。數(shù)值算例驗(yàn)證了本文并行氣熱耦合計(jì)算方法的穩(wěn)定性和收斂性。參考文獻(xiàn)Han J C, Sandip D, Srinath E. Gas turbine heat transfer and cooling technologyB.

17、 CRC Press, 2012.Hylton L D, Milhec M S, Turner E R. Analytical and experimental evaluation of the heat transfer distribution over the surface of turbine vanesB, NASA-CR-168015, 1983.Bohn D, Bonhoff B, Schonenborn H, Wilhelmi H. Validation of a numerical model for the coupled simulation of fluid flo

18、w and diabatic walls with application to film-cooled turbine blades. VDI-BerichtJ, 1995, 1186(1): 272-295.Duchaine F, Corpron A, Moureau V, Nicoud F, Poinsot T. Development and assessment of a coupled strategy for conjugate heat transfer with Large Eddy Simulation. International Journal of Heat and Fluid FlowJ, 2009, 30(6): 1129-1141.周馳，馮國(guó)泰，王松濤，顧中華. 渦輪葉柵氣熱耦合數(shù)值模擬. 工程熱物理學(xué)報(bào)J, 2003, 24(2): 224-226.Hao Z, Gu C, Ren X. The application of Discontinuous Galerkin Methods in conjugate heat transfer simulations of gas turbines. EnergiesJ, 2014, 7(12): 7857-7877