基于流固耦合的部分進(jìn)氣渦輪數(shù)值模擬研究

趙瑞勇，陳 暉，劉軍年，毋 杰

（西安航天動(dòng)力研究所，陜西西安710100）

0 引言

部分進(jìn)氣設(shè)計(jì)的燃?xì)鉁u輪機(jī)內(nèi)部流動(dòng)極為復(fù)雜，由于粘性和復(fù)雜幾何條件引起的激波存在相互耦合，造成了流動(dòng)的非定常性和非穩(wěn)定性，其內(nèi)流場(chǎng)氣動(dòng)特性不同于一般燃?xì)鉁u輪。某型液體火箭發(fā)動(dòng)機(jī)渦輪由于其部分進(jìn)氣設(shè)計(jì)和葉輪高速旋轉(zhuǎn)導(dǎo)致葉輪受到強(qiáng)烈的交變力沖擊，對(duì)葉片應(yīng)力分布產(chǎn)生很大影響?？紤]真實(shí)工況的氣動(dòng)、熱載荷進(jìn)行葉輪強(qiáng)度計(jì)算對(duì)渦輪結(jié)構(gòu)設(shè)計(jì)和工程研制有著重要意義。限于整機(jī)試驗(yàn)研究成本，隨著計(jì)算機(jī)技術(shù)和CFD技術(shù)的發(fā)展，對(duì)渦輪進(jìn)行流固耦合數(shù)值仿真成為研究該問題的重要手段。

在流固耦合仿真計(jì)算中，由于流體域和固體域耦合交界面網(wǎng)格疏密不一致，因此要想實(shí)現(xiàn)載荷傳遞與流固數(shù)據(jù)交換，尋求高效率、小誤差的CSD/CFD數(shù)據(jù)交換方法是實(shí)現(xiàn)耦合技術(shù)的關(guān)鍵。無限平板樣條 (IPS)內(nèi)插值方法和其他插值方法發(fā)展到如今已成為處理機(jī)翼氣動(dòng)彈性計(jì)算數(shù)據(jù)交換較為流行的方法。對(duì)于一般模型，有研究者采用最近取樣法和鄰近節(jié)點(diǎn)加權(quán)平均法實(shí)現(xiàn)流固耦合計(jì)算數(shù)據(jù)插值。在流固交界面網(wǎng)格疏密差異不太大且網(wǎng)格節(jié)點(diǎn)足夠密集的情況下，最近取樣法能夠滿足要求；在流固耦合交界面網(wǎng)格節(jié)點(diǎn)分布比較稀疏的情況下，鄰近點(diǎn)加權(quán)平均法能夠滿足要求，但是耗時(shí)較長(zhǎng)。由于在流固耦合計(jì)算中，流體模型與固體模型網(wǎng)格密度往往不一致或者當(dāng)流體網(wǎng)格發(fā)生形變交界面稀疏情況通常不明確，最近取樣法、鄰近點(diǎn)加權(quán)平均法誤差較大，有時(shí)甚至得到不正確的結(jié)果。文獻(xiàn) [1]中提出的六面體Lagrange插值法和四面體體積插值法對(duì)網(wǎng)格單元類型有限制性。文獻(xiàn) [2-3]提出了一種改進(jìn)的常體積轉(zhuǎn)換法用于非線性氣動(dòng)彈性。

本文發(fā)展了全三維線性插值算法來進(jìn)行CSD/CFD耦合界面數(shù)據(jù)交換。采用全三維線性插值算法對(duì)某型液體火箭發(fā)動(dòng)機(jī)高轉(zhuǎn)速部分進(jìn)氣渦輪進(jìn)行了氣/熱/固多學(xué)科耦合仿真計(jì)算，為工程研制提供了重要的參考價(jià)值。

1 研究對(duì)象

部分進(jìn)氣燃?xì)鉁u輪機(jī)通流部分由噴嘴和工作葉片組成。首先高溫、高壓燃?xì)饨?jīng)過噴嘴將氣體的可用焓降轉(zhuǎn)變?yōu)闅怏w的動(dòng)能；隨后高速燃?xì)庖砸欢ń嵌却迪蛉~輪，使葉輪高速旋轉(zhuǎn)，在此過程中完成了燃?xì)鈩?dòng)能到葉輪動(dòng)能的轉(zhuǎn)變。

圖1 部分進(jìn)氣燃?xì)鉁u輪示意圖Fig.1 Schematic diagram of partial admission turbine

本文研究的某型液體火箭發(fā)動(dòng)機(jī)采用軸流式部分進(jìn)氣單級(jí)沖動(dòng)式渦輪機(jī)，葉片安裝圍帶以減小葉尖漏氣損失，同時(shí)考慮渦輪與氧化劑泵之間的間隙密封（直徑間隙0.06 mm）泄漏。該型渦輪具有高入口溫度、高渦輪壓比、高轉(zhuǎn)速的特點(diǎn)，因此由于其部分進(jìn)氣特點(diǎn)及高氣動(dòng)設(shè)計(jì)參數(shù)所帶來的結(jié)構(gòu)強(qiáng)度問題，對(duì)工程型號(hào)研制至關(guān)重要。傳統(tǒng)工程計(jì)算中，渦輪轉(zhuǎn)子強(qiáng)度計(jì)算氣動(dòng)、熱載荷邊界條件通常采用經(jīng)驗(yàn)預(yù)估或近似給出，無法準(zhǔn)確考慮氣動(dòng)、熱載荷影響。對(duì)某型部分進(jìn)氣渦輪開展氣/熱/固多學(xué)科耦合仿真研究具有重要的工程價(jià)值。

2 CFD/CSD數(shù)據(jù)交換

由于流場(chǎng)計(jì)算網(wǎng)格模型與強(qiáng)度計(jì)算網(wǎng)格模型采用分別建模分網(wǎng)，因此耦合交界面網(wǎng)格疏密程度不一致，且流體網(wǎng)格類型和固體網(wǎng)格類型往往并不統(tǒng)一。發(fā)展了一種三維線性插值算法用于CSD/CFD流固耦合仿真數(shù)據(jù)交換。算法主要原理如下：每個(gè)FE（有限元）網(wǎng)格點(diǎn)選距離最近的10個(gè)CFD（流體）網(wǎng)格點(diǎn)，引入面積控制因子選出其中不小于其值的最近的4個(gè)點(diǎn)，這樣4個(gè)CFD網(wǎng)格點(diǎn)每3個(gè)就可以組成1個(gè)平面三角，依次選取每個(gè)平面三角作為插值平面△A1B1C1，在點(diǎn)A1，B1及C1上分別作△A1B1C1的法線，使得A1A2，B1B2及C1C2分別為點(diǎn)A1，B1及C1的流體載荷值（壓力、溫度）。過F（FE節(jié)點(diǎn)）點(diǎn)做△A1B1C1的法線，F(xiàn)E網(wǎng)格點(diǎn)沿選定平面法向量投影。在△A1B1C1和△A2B2C2上的交點(diǎn)分別為F1，F(xiàn)2，則F1F2為F（FE節(jié)點(diǎn)）點(diǎn)在該插值平面的流體載荷值（壓力、溫度），如圖2所示。對(duì)每一個(gè)F（FE節(jié)點(diǎn)）點(diǎn)在4個(gè)插值平面上求得的4個(gè)F1F2流體載荷值選擇合適的權(quán)函數(shù)取加權(quán)平均值即為該F（FE節(jié)點(diǎn)）點(diǎn)的插值結(jié)果。

圖2 三維線性插值幾何示意圖Fig.2 Geometric schematic diagram of three-dimensional linear interpolation

程序?qū)λ憷哪承蜏u輪葉輪數(shù)值對(duì)比，其結(jié)構(gòu)場(chǎng)有限元壓力、溫度載荷準(zhǔn)確的傳遞了CFD結(jié)果，本文發(fā)展的三維線性插值算法程序，對(duì)于渦輪機(jī)流固耦合計(jì)算過程中的數(shù)據(jù)交換，沒有網(wǎng)格類型限制，插值誤差較小，從而減小了耦合計(jì)算中由于數(shù)據(jù)傳遞帶來的結(jié)果誤差，使計(jì)算結(jié)果可信度高。

3 算例及結(jié)果討論

3.1 某型渦輪定常流場(chǎng)計(jì)算及結(jié)果分析

某型液體火箭發(fā)動(dòng)機(jī)渦輪泵由于其高比功率設(shè)計(jì)，渦輪機(jī)采用單噴嘴部分進(jìn)氣結(jié)構(gòu)，同時(shí)渦輪腔設(shè)計(jì)小間隙（直徑間隙0.06 mm）浮動(dòng)環(huán)密封結(jié)構(gòu)減少工作中的燃?xì)庑孤?。?duì)該型渦輪機(jī)的數(shù)值模擬采用全周模型計(jì)算。限于復(fù)雜結(jié)構(gòu)，仿真計(jì)算中發(fā)展了分塊建模技術(shù)，將計(jì)算域分為7個(gè)塊，最后將各個(gè)計(jì)算域分塊網(wǎng)格采用非一致網(wǎng)格技術(shù)合并為整個(gè)計(jì)算域，實(shí)現(xiàn)了部分進(jìn)氣渦輪考慮微小間隙密封泄漏特性仿真計(jì)算。某型渦輪CFD網(wǎng)格模型如圖3所示。計(jì)算模型中進(jìn)口域、渦輪轉(zhuǎn)子域、排氣管域及密封泄漏腔域5塊采用非結(jié)構(gòu)四面體網(wǎng)格，軸向密封間隙采用6面體網(wǎng)格。全域1 532 807個(gè)單元節(jié)點(diǎn)。

氣動(dòng)、熱力計(jì)算采用FLUENT商用軟件求解3D定常粘性雷諾平均Navier-Stokes方程組。對(duì)控制方程的求解采用基于網(wǎng)格單元中心有限體積法，耦合隱式格式的時(shí)間推進(jìn)算法；對(duì)控制方程對(duì)流項(xiàng)的離散采用2階迎風(fēng)格式。湍流模型選擇RNG k-ε湍流模型。采用高階精度格式對(duì)方程進(jìn)行離散求解，平均殘差小于1×10-6作為收斂判別。動(dòng)/靜交界面采用Mix Plane。

圖3 某型渦輪CFD計(jì)算模型Fig.3 CFD calculation model of one special turbine

圖4分別給出了仿真計(jì)算渦輪馬赫數(shù)、葉片表面壓力及溫度分布。

從結(jié)果看出，總體上葉輪入口燃?xì)庀鄬?duì)速度要比出口大。葉輪入口燃?xì)怦R赫數(shù)最高可達(dá)2.4，葉輪流道入口為噴嘴出來的一斜切口，超音速燃?xì)庠谶@里受到阻滯產(chǎn)生正激波；燃?xì)馀龅饺~片前緣，分開向吸力面和壓力面流動(dòng)，在吸力面氣流繞過前緣以后先有一段加速過程，通道中間處燃?xì)饩S持亞音速流動(dòng)。渦輪轉(zhuǎn)子在工作過程中近噴嘴處轉(zhuǎn)子葉片要承受由于局部進(jìn)氣帶來的熱、氣動(dòng)力不均勻帶來的交變載荷的影響，這在結(jié)構(gòu)設(shè)計(jì)上對(duì)渦輪轉(zhuǎn)子葉片提出了更高的要求。仿真計(jì)算與實(shí)驗(yàn)結(jié)果對(duì)比如表1所示。

圖4 某型渦輪CFD仿真計(jì)算結(jié)果Fig.4 CFD simulation result of a certain turbine

表1 某型渦輪仿真計(jì)算結(jié)果對(duì)比Tab.1 Comparison of simulation results of acertain turbine

3.2 某型渦輪氣、熱、固耦合強(qiáng)度計(jì)算及結(jié)果分析

采用ANSYS有限元程序?qū)δ承蜏u輪轉(zhuǎn)子進(jìn)行氣、熱、固單向耦合強(qiáng)度計(jì)算。渦輪機(jī)采用部分進(jìn)氣方式造成葉片氣動(dòng)、熱載荷分布不均勻，對(duì)葉輪轉(zhuǎn)子進(jìn)行全周數(shù)值仿真。葉輪采用10節(jié)點(diǎn)四面體單元?jiǎng)澐志W(wǎng)格，共213 151個(gè)節(jié)點(diǎn)?？紤]離心力，轉(zhuǎn)速為54 000 r/min。材料為GH4169。有限元計(jì)算網(wǎng)格模型如圖5所示。

圖5 某型渦輪有限元計(jì)算網(wǎng)格Fig.5 FE calculation mesh of a certain turbine

壓力、溫度載荷邊界條件為流場(chǎng)計(jì)算的插值結(jié)果，插值結(jié)果如圖6所示。

圖6 渦輪轉(zhuǎn)子載荷邊界條件（FE）Fig.6 Load boundary condition of turbo rotor

分別對(duì)比圖 4（b）、圖 4（c）、圖 6（a）及圖6（b），可見應(yīng)用三維線性插值程序?qū)δ承蜏u輪CFD/CSD數(shù)據(jù)傳遞都得到了很好的結(jié)果。

圖7給出了考慮氣動(dòng)、溫度、離心力載荷工況下的葉輪von-mises應(yīng)力分布。

圖7 渦輪轉(zhuǎn)子應(yīng)力分布Fig.7 Stress distribution of turbine rotor

由圖7可以看出，某型渦輪轉(zhuǎn)子在設(shè)計(jì)轉(zhuǎn)速、氣動(dòng)參數(shù)工況下，最大應(yīng)力為495 MPa，出現(xiàn)在渦輪轉(zhuǎn)子近噴嘴部分葉片前緣葉根處。渦輪轉(zhuǎn)子前緣與后緣應(yīng)力較為平衡，整體應(yīng)力水平低于材料許用應(yīng)力895 MPa，滿足使用要求。

4 結(jié)論

發(fā)展了三維線性插值算法用于CFD/CSD數(shù)據(jù)交換，對(duì)某型液體火箭發(fā)動(dòng)機(jī)高轉(zhuǎn)速部分進(jìn)氣渦輪進(jìn)行了氣動(dòng)、熱、強(qiáng)度耦合仿真，得出以下結(jié)論：發(fā)展的三維線性插值算法對(duì)網(wǎng)格限制性小，結(jié)果精確，為多學(xué)科耦合仿真計(jì)算提供了數(shù)據(jù)交換工具，采用該方法對(duì)渦輪機(jī)進(jìn)行氣動(dòng)、熱、固耦合數(shù)值仿真具有較高的工程參考價(jià)值。對(duì)某型液體火箭發(fā)動(dòng)機(jī)部分進(jìn)氣渦輪機(jī)計(jì)算分析認(rèn)為，葉輪入口燃?xì)庀鄬?duì)馬赫數(shù)超音最高可達(dá)2.4，超音速燃?xì)庠谌~輪流道入口斜切口受到阻滯，產(chǎn)生正激波；近噴嘴處轉(zhuǎn)子葉片要承受由于局部進(jìn)氣帶來的熱、氣動(dòng)力不均勻帶來的交變載荷的影響，在設(shè)計(jì)轉(zhuǎn)速、氣動(dòng)參數(shù)工況下，葉輪最大應(yīng)力為495 MPa，出現(xiàn)在渦輪轉(zhuǎn)子近噴嘴部分葉片前緣葉根處，滿足強(qiáng)度使用要求，為工程研制提供了仿真依據(jù)。

[1]汪學(xué)鋒等.流固耦合網(wǎng)格插值方法研究 [J].船舶力學(xué),2009,13(4)：571-578.

[2]徐敏,史忠軍,陳士櫓.一種流體-結(jié)構(gòu)耦合計(jì)算問題的網(wǎng)格數(shù)據(jù)交換方法[J].西北工業(yè)大學(xué)學(xué)報(bào),2003，21(5)：533-535.

[3]徐敏,陳士櫓.CFD/CSD耦合計(jì)算研究[J].應(yīng)用力學(xué)學(xué)報(bào),2004,21(2)：33-36.

[4]SMITH M J,HODGES D H.Evaluation of computational algorithms suitable for fluid-structure interactions[J].Journal of Aircraft,2000,37(2)：282-294.

[5]MARTENSSON Hans,LAUMERT Bjorn,FRANSSON T H.Aeromechanical aspects on unsteady flow in turbines[C]//Proceedings of 33rd AIAA Fluid Dynamics Conference and Exhibit.Oriando,Florida：AIAA,2003：3997-4002.

[6]李立州.基于網(wǎng)格變形技術(shù)的渦輪葉片變形傳遞[J].航空動(dòng)力學(xué)報(bào),2007,22(12)：2101-2104.

[7]趙瑞勇,王延榮.葉片顫數(shù)值模擬方法研究[D].北京：北京航空航天大學(xué),2010.

[8]胡運(yùn)聰,周新海.振動(dòng)葉柵非定常流動(dòng)數(shù)值模擬與葉片顫振分析[D].西安：西北工業(yè)大學(xué),2003.

[9]伊進(jìn)寶.部分進(jìn)氣燃?xì)鉁u輪機(jī)葉輪流場(chǎng)數(shù)值模擬[J].魚雷技術(shù),2010,18(6)：456-460.