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

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

（西安航天動力研究所，陜西西安710100）

0 引言

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

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

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

1 研究對象

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[3]徐敏,陳士櫓.CFD/CSD耦合計算研究[J].應(yīng)用力學(xué)學(xué)報,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].航空動力學(xué)報,2007,22(12)：2101-2104.

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

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

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