專利名稱：葉輪機械計入葉間相位角的氣彈穩(wěn)定性流固耦合預(yù)測方法
技術(shù)領(lǐng)域：
本發(fā)明屬于葉輪機械模擬技術(shù)領(lǐng)域，具體涉及基于能量法的一種葉輪機械計入葉 間相位角的氣彈穩(wěn)定性流固耦合預(yù)測方法。
背景技術(shù)：
科學(xué)技術(shù)的進步使得葉輪機械的葉尖切線速度和壓比等不斷提高，導(dǎo)致風(fēng)扇/壓 氣機葉片氣彈穩(wěn)定性問題突顯。隨著葉輪機械向著大功率、高性能的發(fā)展，風(fēng)扇/壓氣機葉 片的工作條件越來越惡劣，葉尖切線速度越來越高，葉片的剛性相對越來越小，葉片發(fā)生氣 彈不穩(wěn)定的可能性越來越大。國內(nèi)外許多發(fā)動機均發(fā)生過風(fēng)扇/壓氣機葉片的氣彈穩(wěn)定性 故障，葉輪機械葉片氣彈穩(wěn)定性問題成為現(xiàn)代航空發(fā)動機設(shè)計所面臨的主要難點之一。20世紀80年代以來，隨著計算機的發(fā)展和計算技術(shù)的進步，出現(xiàn)了數(shù)值方法，而 計算流體力學(xué)的發(fā)展更進一步推動了振蕩葉柵氣彈穩(wěn)定性問題的研究。目前，研究最多的 數(shù)值方法主要為能量法和特征值法。由于葉輪機械內(nèi)部的流場呈現(xiàn)著復(fù)雜的三維非定常特 性，而且分離流動、旋轉(zhuǎn)失速、激波、激波與附面層的相互干擾、激波與激波間的相互干擾等 非線性因素的影響使得氣動載荷很難用準確的函數(shù)來表達，通過特征值法來定量地準確描 述氣彈穩(wěn)定性幾乎不可能實現(xiàn)，因此，采用能量法能夠很好地保證流場計算精度。能量法忽略了流體和固體之間的相互作用而僅僅考慮固體振動對流場的影響，因 此，分別對流體域和固體域單獨求解，并通過發(fā)展的基于有限元形函數(shù)的概念實現(xiàn)流體域 和固體域之間的數(shù)據(jù)傳遞，這也是氣彈穩(wěn)定性數(shù)值預(yù)測技術(shù)的關(guān)鍵。另一方面，由于葉片的振蕩會引起葉柵之間非定常流場的相互作用，對葉片氣彈 穩(wěn)定性有很大的影響。早期對葉輪機械的氣動彈性問題大多假設(shè)葉柵中的所有葉片均以相 同的振幅振蕩，并且不計葉間相位角的影響，但是隨著對氣動彈性研究的深入，發(fā)現(xiàn)葉間相 位角是影響葉片氣動彈性的一個相當關(guān)鍵的參數(shù)，它對葉片氣彈穩(wěn)定性起著重要的作用， 而且實際葉片并不可能滿足上述的假設(shè)。因此，對于葉輪機械葉片氣動彈性這樣的流固耦合問題，有必要發(fā)展和完善適用 于工程設(shè)計需求的氣彈穩(wěn)定性數(shù)值預(yù)測方法，并考慮葉間相位角等關(guān)鍵參數(shù)的影響。其中 最大的困難是如何準確高效地實現(xiàn)流固耦合交界面上的數(shù)據(jù)傳遞，并在振蕩葉柵中計入葉 間相位角。目前還沒有一項較為成熟的技術(shù)來預(yù)測計入葉間相位角影響的葉輪機械葉片的 氣彈穩(wěn)定性。

發(fā)明內(nèi)容
本發(fā)明的目的是提供葉輪機械計入葉間相位角的氣彈穩(wěn)定性流固耦合預(yù)測方法， 它解決了現(xiàn)有技術(shù)的不足。本發(fā)明選擇精度和效率平衡的能量法作為模型，它滿足葉片小幅振動的假設(shè)，并 且忽略流體對振動系統(tǒng)的影響，而僅僅考慮固體域?qū)α黧w域的影響，能夠降低計算成本。本發(fā)明葉輪機械計入葉間相位角的氣彈穩(wěn)定性流固耦合預(yù)測方法，其具體步驟如下步驟一建立單扇區(qū)有限元模型；針對給定的實體模型在有限元軟件MARC中建立葉輪機械的單扇區(qū)有限元模型；步驟二 對建立的有限元模型進行含預(yù)應(yīng)力的模態(tài)分析；首先，將MARC中建立的有限元模型導(dǎo)入到有限元軟件ANSYS中，定義相應(yīng)的材料 參數(shù)，并給定轉(zhuǎn)速以及位移約束條件，通過靜力分析得到離心力引起的預(yù)應(yīng)力；然后，采用振型歸一化的模態(tài)分析，并計入離心力引起的預(yù)應(yīng)力，獲得葉片的固有 振動特性，提取葉片振動固有頻率以及模態(tài)振型，得到固體域葉片表面的節(jié)點振幅；步驟三建立計算流體力學(xué)模型；首先，通過給定的葉輪機械葉型數(shù)據(jù)，在計算流體力學(xué)仿真工具CFX的TurboGrid 模塊中建立葉柵的單通道流場模型，其中葉片表面附近為0型網(wǎng)格，其他區(qū)域為H型網(wǎng)格；然后，將建立的單扇區(qū)流場模型導(dǎo)入到ICEM CFD中，通過繞旋轉(zhuǎn)軸旋轉(zhuǎn)、復(fù)制 (N-I)份網(wǎng)格得到全環(huán)葉柵模型，并輸出流體域網(wǎng)格點和實體單元信息以及流體域各葉片 表面網(wǎng)格點和表面單元信息；步驟四采用數(shù)據(jù)傳遞方法獲得插值到流體域葉片表面的振動位移；首先，基于恒定位置矢量差方法，采用有限單元法中局部坐標形式的形函數(shù)進行 三維線性插值，將固體域中模態(tài)分析得到的第P階模態(tài)對應(yīng)的節(jié)點振動位移插值到流固耦 合交界面上的表面單元網(wǎng)格點上，獲得流體域中葉片表面網(wǎng)格點的振幅；然后，通過流體域中各葉片扇區(qū)之間的角度關(guān)系，即可依次逐個得到全環(huán)各葉片 表面網(wǎng)格點的振幅；步驟五采用動網(wǎng)格技術(shù)獲得每一個葉片扇區(qū)的可動域網(wǎng)格點振幅；根據(jù)流體域中網(wǎng)格點編號順序，對每一個葉片沿0型域網(wǎng)格的法向逐層遞推搜 索，得到與該葉片表面網(wǎng)格點編號順序相同的各層網(wǎng)格點，并按照各層網(wǎng)格點與葉片表面 網(wǎng)格點之間法向距離的比值大小來線性插值，得到該葉片可動域的網(wǎng)格點振幅；步驟六獲得非定常計算流體力學(xué)(CFD)分析所需的考慮葉間相位角的網(wǎng)格文 件；首先，將葉片一個振動周期分為K步，葉片以簡諧規(guī)律振動，給定各葉片相對于第 一個葉片的初始相位角，通過初始坐標以及各時間步的振動位移，可以得到一個周期內(nèi)每 一時間步上各葉片扇區(qū)可動域的坐標值；然后，通過固定域坐標值和每一時間步上各葉片扇區(qū)可動域的坐標值，可以得到 一個周期內(nèi)每一時間步的流體網(wǎng)格點坐標組成的網(wǎng)格文件；步驟七采用CFX中Junction Box模塊進行非定常分析；流場分析通過計算流體力學(xué)仿真工具CFX實現(xiàn)；首先，將步驟三建立的流場模型的網(wǎng)格點坐標值換成步驟六中獲得的第一個時間 步的網(wǎng)格點坐標值，對流場模型加載邊界條件和初始條件；通道入口給定總溫和總壓，出口 給定平均靜壓，輪緣、輪轂給定無滑移、光滑壁面邊界條件，葉片表面給定動網(wǎng)格壁面邊界 條件；以葉片表面給定無滑移、光滑壁面邊界條件的定常解作為非定常分析的初始條件；然后，指定葉片的振動模態(tài)以及振動幅值，其運動周期為固有頻率的倒數(shù)，在每一 個振動時間步上求解采用k- ε湍流模型封閉的Reynolds平均Navier-Stokes方程，并寫出每一個求解時間步的瞬時結(jié)果文件；步驟八獲得各葉間相位角對應(yīng)的模態(tài)氣動阻尼比，預(yù)測振蕩葉柵的氣彈穩(wěn)定 性；首先，在不同葉間相位角情況下，根據(jù)一個穩(wěn)定周期內(nèi)各時間步的結(jié)果文件，基 于等效粘性阻尼推導(dǎo)的模態(tài)氣動阻尼比獲得葉片一個振動周期的氣動功和模態(tài)氣動阻尼 比；其中氣動功是流場各網(wǎng)格點氣動力在一個穩(wěn)定振動周期內(nèi)所做的總功； 基于等效粘性阻尼推導(dǎo)的模態(tài)氣動阻尼比公式為么㈤=-(aero、2，其中W為非定 常氣動力在一個周期內(nèi)所做的氣動功，qaer°為葉片的振動幅值；然后，通過不同葉間相位角情況下的模態(tài)氣動阻尼比，判斷某階模態(tài)下某一工況 的氣彈穩(wěn)定性，在各葉間相位角下，只要有某一個葉間相位角對應(yīng)的模態(tài)氣動阻尼比為負， 則該階模態(tài)下，該工況氣彈失穩(wěn)。這種基于能量法的模擬方法適用于葉輪機械全環(huán)葉柵考慮葉間相位角影響的氣 彈穩(wěn)定性分析，在保證了一定計算精度的基礎(chǔ)上，能夠很好地提高計算效率，節(jié)約計算資源 和減少時間成本。這種方法中所用到的基于有限元形函數(shù)的數(shù)據(jù)傳遞方法準確并且效率 高。而為了保證計算效率，將流體域分為固定域和每個葉片表面附近的可動域，并在每一個 可動域采用動網(wǎng)格技術(shù)，摒棄CFX本身只能運動表面網(wǎng)格節(jié)點的約束，使得葉片能夠在更 大的振幅下振動，同時通過節(jié)徑數(shù)給出全環(huán)葉柵各葉片的葉間相位角，并通過不同葉間相 位角下的模態(tài)氣動阻尼比判斷葉輪機械葉片的氣彈穩(wěn)定性。下面分別介紹該數(shù)值預(yù)測方法中用到的數(shù)據(jù)傳遞方法和考慮葉間相位角的網(wǎng)格 文件生成方法。數(shù)據(jù)傳遞方法的具體實施步驟如下步驟一獲得與第一個葉片表面的流體網(wǎng)格點F距離最近的固體節(jié)點S ；對流體域中第一個葉片表面的任意網(wǎng)格點F，在固體域中找到與其距離最近的固 體表面節(jié)點S，與節(jié)點S相關(guān)的單元至多有4個，分別為Si、S2, S3和S4,并從固體域中獲得 各單元對應(yīng)的結(jié)構(gòu)單元的表面編號；步驟二 獲得流體網(wǎng)格點F在不同固體單元中的局部坐標值；流體網(wǎng)格點F的坐標為(X。，y0, Z0)，固體單元八個節(jié)點的坐標分別為(Xi，yi Zi)， (i = 1，...，8)，每個節(jié)點的形函數(shù)在局部坐標系下可以表示為隊(ε，η，ζ)，(i = 1，..， 8)，則通過流體網(wǎng)格點坐標和固體單元節(jié)點坐標之間的關(guān)系
可以得到流體網(wǎng)格點F在與節(jié)點S相關(guān)的不同單元的局部坐標值 其中局部坐標系下的形函數(shù)可以表示為 Ν( ε，η，ζ) = (1士 ε ) (1士 η) (1士 ζ)
步驟三獲得流體網(wǎng)格點F對應(yīng)固體域的投影單元M ；分別對比不同單元對應(yīng)的局部坐標值，找出各局部坐標值的平方和最小值對應(yīng)的 單元M以及該單元的表面編號，即為流體網(wǎng)格點F在固體域中的投影單元和單元面，將網(wǎng)格 點F投影到單元M的表面得到F'，由于流體域網(wǎng)格和固體域網(wǎng)格對原始幾何構(gòu)型的不同精 度導(dǎo)致網(wǎng)格點F和投影點F'之間存在坐標差，采用恒定位置矢量差法，使流體網(wǎng)格點F的 函數(shù)值始終等于其投影點F'的函數(shù)值；步驟四插值得到流體網(wǎng)格點F的位移；首先，獲得流體網(wǎng)格點F在投影單元M上的局部坐標值，并進一步得到單元M的八 個形函數(shù)隊(￡，η，ζ) = (1士 O (1士 η) (1士 ζ)，(i = 1，···，8)；然后，通過固體域中實體單元的八個節(jié)點振動位移(Ui, Vi，Wi) (i = 1，. . .，8)可以
插值獲得流固耦合交界面上流體網(wǎng)格點F的位移(Utl，v0, Wtl)，具體的位移插值為 8 步驟五獲得其他葉片的表面網(wǎng)格點振動位移；首先，獲得與第一個葉片表面網(wǎng)格點編號相對應(yīng)的其他各葉片的表面網(wǎng)格點編 號；在旋轉(zhuǎn)坐標軸為z的情況下，各葉片逆時針編號，第i個葉片表面網(wǎng)格點坐標(Xi， Yi, Zi)與第一個葉片表面網(wǎng)格點坐標(X1, Y1, Z1)之間滿足如下關(guān)系
;當?shù)谝粋€葉片在其他象其中當?shù)谝粋€葉片在第一象限時，
限時，按照對應(yīng)的角度關(guān)系給出Q1 ； 2π其中
N為全環(huán)葉片個數(shù)；然后，將第一個葉片表面網(wǎng)格點的振幅分別賦值給其他各葉片的對應(yīng)網(wǎng)格點，獲 得全環(huán)所有葉片表面網(wǎng)格點的振幅?？紤]葉間相位角的網(wǎng)格文件生成方法具體實施步驟如下步驟一將葉片的一個振動周期分為K步，其中一個周期對應(yīng)的頻率為葉片在某 階模態(tài)下的固有頻率ω，給定第一個葉片正弦振動的初始相位角為零，則第i個葉片在第k 個時間步的相位角為
其中，1為節(jié)徑數(shù)，N為全環(huán)葉片個數(shù)；則該葉片扇區(qū)的每一個網(wǎng)格點在該時間步的振動位移為
rUlk =UlSiiKplk
vik = V1 sin 灼 ^lk =^ sin^4其中，Ui, Vi, Wi分別為第i個葉片扇區(qū)的某一個網(wǎng)格點在x，y，ζ三個方向的振幅；步驟二 通過動網(wǎng)格技術(shù)搜索獲得第i個葉片對應(yīng)可動域的網(wǎng)格點編號，并從全 域網(wǎng)格點坐標值可以得到該葉片扇區(qū)的可動域網(wǎng)格點初始坐標值(Xi，Yi, Zi)，而該葉片扇 區(qū)的每一個網(wǎng)格點在第k個時間步的X，y，ζ三個方向振動位移分別為uik，vik, wik,則該葉 片扇區(qū)的每一個網(wǎng)格點在第k個時間步的坐標值為

步驟三通過固定域坐標值和每一時間步上各葉片扇區(qū)可動域的坐標值，可以得 到一個周期內(nèi)每一時間步的流體網(wǎng)格點坐標，并將每一個時間步的網(wǎng)格點坐標值輸出為網(wǎng) 格文件。本發(fā)明的優(yōu)點在于(1)基于能量法的流固耦合預(yù)測方法充分借助當前計算結(jié)構(gòu)動力學(xué)和計算流體 動力學(xué)的前沿技術(shù)來實現(xiàn)各部分的單獨計算，并通過設(shè)計的數(shù)據(jù)傳遞方法來進行流體域和 固體域的數(shù)據(jù)傳遞，從而使得耦合計算大大簡化，并在每個子系統(tǒng)內(nèi)均能保持計算的高精 度；(2)基于恒定位置矢量差方法，采用有限元法中局部坐標形式的形函數(shù)進行三維 線性插值，方法簡單，雖然會在流體網(wǎng)格造成一定的橫向剪切變形，但是計算效率較高；(3)將全環(huán)葉柵模型分為葉片表面附近的可動域和其他區(qū)域的固定域能夠很好地 保證計算效率，而在通過節(jié)徑數(shù)給出的葉間相位角情況下計算全范圍內(nèi)的模態(tài)氣動阻尼比 能夠更準確地用于判斷葉輪機械葉片的氣彈穩(wěn)定性。


圖1是本發(fā)明計入葉間相位角的氣彈穩(wěn)定性流固耦合預(yù)測方法流程圖；圖2是現(xiàn)有MARC中的NASA67葉片有限元模型；圖3是現(xiàn)有NASA67葉片第一階彎曲模態(tài)；圖4是現(xiàn)有NASA67葉片全環(huán)葉柵流場模型；圖5是本發(fā)明數(shù)據(jù)傳遞方法流程圖；圖6是現(xiàn)有NASA67葉片單扇區(qū)葉柵的可動域和固定域；圖7是本發(fā)明動網(wǎng)格技術(shù)示意圖8是本發(fā)明模態(tài)氣動阻尼比與葉間相位角的關(guān)系圖； 具體實施例方式下面將結(jié)合附圖和實施例對本發(fā)明作進一步的詳細說明。本發(fā)明是葉輪機械計入 葉間相位角的氣彈穩(wěn)定性流固耦合預(yù)測方法，方法的流程如圖1所示。步驟一建立單扇區(qū)有限元模型；針對給定的實體模型在有限元軟件MARC中建立葉輪機械的單扇區(qū)有限元模型；步驟二 對建立的有限元模型進行含預(yù)應(yīng)力的模態(tài)分析；首先，將MARC中建立的有限元模型導(dǎo)入到有限元軟件ANSYS中，定義相應(yīng)的材料 參數(shù)，并給定轉(zhuǎn)速以及位移約束條件，通過靜力分析得到離心力引起的預(yù)應(yīng)力；然后，采用振型歸一化的模態(tài)分析，并計入離心力引起的預(yù)應(yīng)力，獲得葉片的固有 振動特性，提取葉片振動固有頻率以及模態(tài)振型，得到固體域葉片表面的節(jié)點振幅；步驟三建立計算流體力學(xué)模型；首先，通過給定的葉輪機械葉型數(shù)據(jù)，在計算流體力學(xué)仿真工具CFX的TurboGrid 模塊中建立葉柵的單通道流場模型，其中葉片表面附近為0型網(wǎng)格，其他區(qū)域為H型網(wǎng)格；然后，將建立的單扇區(qū)流場模型導(dǎo)入到ICEM CFD中，通過繞旋轉(zhuǎn)軸旋轉(zhuǎn)、復(fù)制 (N-I)份網(wǎng)格得到全環(huán)葉柵模型，并輸出流體域網(wǎng)格點和實體單元信息以及流體域各葉片 表面網(wǎng)格點和表面單元信息；步驟四采用數(shù)據(jù)傳遞方法獲得插值到流體域葉片表面的振動位移；首先，基于恒定位置矢量差方法，采用有限單元法中局部坐標形式的形函數(shù)進行 三維線性插值，將固體域中模態(tài)分析得到的第P階模態(tài)對應(yīng)的節(jié)點振動位移插值到流固耦 合交界面上的表面單元網(wǎng)格點上，獲得流體域中葉片表面網(wǎng)格點的振幅；然后，通過流體域中各葉片扇區(qū)之間的角度關(guān)系，即可依次逐個得到全環(huán)各葉片 表面網(wǎng)格點的振幅；步驟五采用動網(wǎng)格技術(shù)獲得每一個葉片扇區(qū)的可動域網(wǎng)格點振幅；根據(jù)流體域中網(wǎng)格點編號順序，對每一個葉片沿0型域網(wǎng)格的法向逐層遞推搜 索，得到與該葉片表面網(wǎng)格點編號順序相同的各層網(wǎng)格點，并按照各層網(wǎng)格點與葉片表面 網(wǎng)格點之間法向距離的比值大小來線性插值，得到該葉片可動域的網(wǎng)格點振幅；步驟六獲得非定常計算流體力學(xué)(CFD)分析所需的考慮葉間相位角的網(wǎng)格文 件；首先，將葉片一個振動周期分為K步，葉片以簡諧規(guī)律振動，給定各葉片相對于第 一個葉片的初始相位角，通過初始坐標以及各時間步的振動位移，可以得到一個周期內(nèi)每 一時間步上各葉片扇區(qū)可動域的坐標值；然后，通過固定域坐標值和每一時間步上各葉片扇區(qū)可動域的坐標值，可以得到 一個周期內(nèi)每一時間步的流體網(wǎng)格點坐標組成的網(wǎng)格文件；步驟七采用CFX中Junction Box模塊進行非定常分析；流場分析通過計算流體力學(xué)仿真工具CFX實現(xiàn)；首先，將步驟三建立的流場模型的網(wǎng)格點坐標值換成步驟六中獲得的第一個時間 步的網(wǎng)格點坐標值，對流場模型加載邊界條件和初始條件；通道入口給定總溫和總壓，出口給定平均靜壓，輪緣、輪轂給定無滑移、光滑壁面邊界條件，葉片表面給定動網(wǎng)格壁面邊界 條件；以葉片表面給定無滑移、光滑壁面邊界條件的定常解作為非定常分析的初始條件；然后，指定葉片的振動模態(tài)以及振動幅值，其運動周期為固有頻率的倒數(shù)，在每一 個振動時間步上求解采用k- ε湍流模型封閉的Reynolds平均Navier-Stokes方程，并寫 出每一個求解時間步的瞬時結(jié)果文件；步驟八獲得各葉間相位角對應(yīng)的模態(tài)氣動阻尼比，預(yù)測振蕩葉柵的氣彈穩(wěn)定 性；首先，在不同葉間相位角情況下，根據(jù)一個穩(wěn)定周期內(nèi)各時間步的結(jié)果文件，基 于等效粘性阻尼推導(dǎo)的模態(tài)氣動阻尼比獲得葉片一個振動周期的氣動功和模態(tài)氣動阻尼 比；其中氣動功是流場各網(wǎng)格點氣動力在一個穩(wěn)定振動周期內(nèi)所做的總功；
-W基于等效粘性阻尼推導(dǎo)的模態(tài)氣動阻尼比公式為《 。=7^^7，其中W為非定
2徹 )
常氣動力在一個周期內(nèi)所做的氣動功，qaer°為葉片的振動幅值；然后，通過不同葉間相位角情況下的模態(tài)氣動阻尼比，判斷某階模態(tài)下某一工況 的氣彈穩(wěn)定性，在各葉間相位角下，只要有某一個葉間相位角對應(yīng)的模態(tài)氣動阻尼比為負， 則該階模態(tài)下，該工況氣彈失穩(wěn)。這種基于能量法的模擬方法適用于葉輪機械全環(huán)葉柵考慮葉間相位角影響的氣 彈穩(wěn)定性分析，在保證了一定計算精度的基礎(chǔ)上，能夠很好地提高計算效率，節(jié)約計算資源 和減少時間成本。這種方法中所用到的基于有限元形函數(shù)的數(shù)據(jù)傳遞方法準確且高效。而 為了保證計算效率，將流體域分為固定域和葉片表面附近的可動域，并在每一個可動域保 證葉片能夠在更大的振幅下振動，同時通過節(jié)徑數(shù)給出全環(huán)葉柵各葉片的葉間相位角，并 通過不同葉間相位角下的模態(tài)氣動阻尼比判斷葉輪機械葉片的氣彈穩(wěn)定性。下面通過一個實例來說明此數(shù)值方法。選用NASA67轉(zhuǎn)子葉片，其基本的設(shè)計氣動 參數(shù)見表1?；诮o定的葉型數(shù)據(jù)在有限元軟件MARC中建立單個葉片的有限元模型如圖 2，其中單元選取八節(jié)點六面體實體單元。表INASA 67轉(zhuǎn)子葉片基本設(shè)計參數(shù) 將MARC中建立的有限元模型導(dǎo)入到有限元軟件ANSYS中進行含預(yù)應(yīng)力的模態(tài)分 析。其中材料選取鈦合金，邊界條件給定葉根固支。首先，通過靜力分析得到離心力引起的 預(yù)應(yīng)力，然后，進行含預(yù)應(yīng)力的模態(tài)分析，得到葉片的各階動頻和模態(tài)，如圖3為葉片的第 一階彎曲模態(tài)。
通過給定的葉型數(shù)據(jù)在CFX TurboGrid中建立流體域單通道模型，并通過旋轉(zhuǎn)復(fù) 制得到全環(huán)葉柵模型，如圖4，其中，模型忽略了葉尖間隙的作用。采用恒定位置矢量差法，基于有限元局部坐標系下形函數(shù)的概念發(fā)展了數(shù)據(jù)傳遞 方法，將單個葉片的第一階彎曲模態(tài)插值到流體域的流固耦合交界面的葉片表面網(wǎng)格點 上，其數(shù)據(jù)傳遞的流程如圖5，具體實施步驟如下步驟一獲得與第一個葉片表面的流體網(wǎng)格點F距離最近的固體節(jié)點S ；對流體域中第一個葉片表面的任意網(wǎng)格點F，在固體域中找到與其距離最近的固 體表面節(jié)點S，與節(jié)點S相關(guān)的單元至多有4個，分別為Sp S2, S3和S4,并從固體域中獲得 各單元對應(yīng)的結(jié)構(gòu)單元的表面編號；步驟二 獲得流體網(wǎng)格點F在不同固體單元中的局部坐標值；流體網(wǎng)格點F的坐標為(X(1，y0, z0)，固體單元八個節(jié)點的坐標分別為(Xi，yi Zi)， (i = 1，...，8)，每個節(jié)點的形函數(shù)在局部坐標系下可以表示為隊(ε，η，ζ)，(i = 1，..， 8)，則通過流體網(wǎng)格點坐標和固體單元節(jié)點坐標之間的關(guān)系
可以得到流體網(wǎng)格點F在與節(jié)點S相關(guān)的不同單元的局部坐標值； 其中局部坐標系下的形函數(shù)可以表示為 Ν( ε，η，ζ) = (1士 ε ) (1士 η) (1士 ζ)步驟三獲得流體網(wǎng)格點F對應(yīng)固體域的投影單元M ；分別對比不同單元對應(yīng)的局部坐標值，找出各局部坐標值的平方和最小值對應(yīng)的 單元M以及該單元的表面編號，即為流體網(wǎng)格點F在固體域中的投影單元和單元面，將網(wǎng)格 點F投影到單元M的表面得到F'，由于流體域網(wǎng)格和固體域網(wǎng)格對原始幾何構(gòu)型的不同精 度導(dǎo)致網(wǎng)格點F和投影點F'之間存在坐標差，采用恒定位置矢量差法，使流體網(wǎng)格點F的 函數(shù)值始終等于其投影點F'的函數(shù)值；步驟四插值得到流體網(wǎng)格點F的位移；首先，獲得流體網(wǎng)格點F在投影單元M上的局部坐標值，并進一步得到單元M的八 個形函數(shù)隊(￡，η，ζ) = (1士 O (1士 η) (1士 ζ)，(i = 1，···，8)；然后，通過固體域中實體單元的八個節(jié)點振動位移(Ui, Vi，Wi) (i = 1，...，8)可以 插值獲得流固耦合交界面上流體網(wǎng)格點F的位移Ov v0, Wtl)，具體的位移插值為 步驟五獲得其他葉片的表面網(wǎng)格點振動位移
首先，獲得與第一個葉片表面網(wǎng)格點編號相對應(yīng)的其他各葉片的表面網(wǎng)格點編 號；在旋轉(zhuǎn)坐標軸為ζ的情況下，各葉片逆時針編號，第i個葉片表面網(wǎng)格點坐標(Xi， Yi, Zi)與第一個葉片表面網(wǎng)格點坐標(Xl，Y1, Z1)之間滿足如下關(guān)系
Xj = ^jxf + y\ sin(0j + Oi) 其中當?shù)谝粋€葉片在第一象限時，
當?shù)谝粋€葉片在其他象
限時，按照對應(yīng)的角度關(guān)系給出θ工； 2π其中，Θ, = ^r( -1)，N為全環(huán)葉片個數(shù)；然后，將第一個葉片表面網(wǎng)格點的振幅分別賦值給其他各葉片的對應(yīng)網(wǎng)格點，獲 得全環(huán)所有葉片表面網(wǎng)格點的振幅。在振蕩葉柵作用下的非定常流場分析時，僅將振蕩葉片附近的流體域指定為可動 域，用0型網(wǎng)格來實現(xiàn)，0型可動域的網(wǎng)格點坐標在每一求解時刻均隨著葉片的振蕩實時 更新，遠離葉片的其他區(qū)域均被定義為固定域，網(wǎng)格點坐標始終保持不變，用H型網(wǎng)格來實 現(xiàn)。在振蕩葉柵作用下的流場分析中，將每一個葉柵通道的流體域網(wǎng)格劃分為葉片表面附 近的0型網(wǎng)格和其他區(qū)域的H型網(wǎng)格如圖6。從網(wǎng)格點的編號來看，每個葉片扇區(qū)的0型可動域的各層網(wǎng)格點編號規(guī)律相同， 而且對應(yīng)網(wǎng)格點位于同一直線上，如圖7，當葉片振動時，流體域中固體邊界附近的網(wǎng)格點 坐標均會按法向距離的大小比例變化，且始終能夠保證流體域表面一層的網(wǎng)格尺寸大小滿 足精度要求。當獲得每個葉片附近可動域的振幅后，生成考慮葉間相位角的網(wǎng)格文件，生成 方法的具體實施步驟如下步驟一將葉片的一個振動周期分為K步，其中一個周期對應(yīng)的頻率為葉片在某 階模態(tài)下的固有頻率ω，給定第一個葉片正弦振動的初始相位角為零，則第i個葉片在第k 個時間步的相位角為 其中，1為節(jié)徑數(shù)，N為全環(huán)葉片個數(shù)；則該葉片扇區(qū)的每一個網(wǎng)格點在該時間步的振動位移為 其中，Ui, Vi, Wi分別為第i個葉片扇區(qū)的某一個網(wǎng)格點在x，y，ζ三個方向的振幅；步驟二 通過動網(wǎng)格技術(shù)搜索獲得第i個葉片對應(yīng)可動域的網(wǎng)格點編號，并從全域網(wǎng)格點坐標值可以得到該葉片扇區(qū)的可動域網(wǎng)格點初始坐標值(Xi，Ii, Zi)，而該葉片扇 區(qū)的每一個網(wǎng)格點在第k個時間步的X，y，ζ三個方向振動位移分別為uik，vik, wik,則該葉 片扇區(qū)的每一個網(wǎng)格點在第k個時間步的坐標值為
Xik = Xi + Uik< yik=y,+vlk
z,k = z, + w,k步驟三通過固定域坐標值和每一時間步上各葉片扇區(qū)可動域的坐標值，可以得 到一個周期內(nèi)每一時間步的流體網(wǎng)格點坐標，并將每一個時間步的網(wǎng)格點坐標值輸出為網(wǎng) 格文件。得到計入葉間相位角影響的各時間步網(wǎng)格文件后，在CFX的Junction Box動網(wǎng)格 模塊中進行非定常分析，其中以葉片不振蕩、考慮葉間相位角的全環(huán)葉柵模型在設(shè)計狀態(tài) 的邊界條件下計算得到的定常解作為初始條件，并設(shè)置設(shè)計狀態(tài)的非定常邊界條件，在每 一個求解時間步讀取網(wǎng)格文件，實現(xiàn)葉片的振蕩作用。CFX非定常分析后，讀取一個穩(wěn)定周期內(nèi)各時間步的葉片表面氣動力以及相應(yīng)的 振蕩位移，進而求得一個周期內(nèi)的非定常氣動功，并得到給定葉間相位角情況下基于等效 粘性阻尼推導(dǎo)的模態(tài)氣動阻尼比，進而曲線擬合得到不同葉間相位角對應(yīng)的模態(tài)氣動阻尼 比，如圖8。圖中模態(tài)氣動阻尼比與葉間相位角近似呈簡諧規(guī)律，且在接近0°時模態(tài)氣動 阻尼比達到最小值，在接近180°時模態(tài)氣動阻尼比達到最大值，而且在整個葉間相位角范 圍內(nèi)，模態(tài)氣動阻尼比均為正，表示NASA67在葉片的一階彎曲模態(tài)下不會出現(xiàn)氣彈失穩(wěn)的 現(xiàn)象，這與實驗值是吻合的。對于其他的葉輪機械葉片的氣彈穩(wěn)定性判斷，可以通過以上的方法，在葉片的某 一階模態(tài)下，得到不同葉間相位角情況下的模態(tài)氣動阻尼比，只有當整個葉間相位角范圍 內(nèi)模態(tài)氣動阻尼比均為正時，葉輪機械葉片才氣彈穩(wěn)定。
權(quán)利要求
葉輪機械計入葉間相位角的氣彈穩(wěn)定性流固耦合預(yù)測方法，其特征在于，該方法包括如下步驟步驟一建立單扇區(qū)有限元模型；針對給定的實體模型在有限元軟件MARC中建立葉輪機械的單扇區(qū)有限元模型；步驟二對建立的有限元模型進行含預(yù)應(yīng)力的模態(tài)分析；首先，將MARC中建立的有限元模型導(dǎo)入到有限元軟件ANSYS中，定義相應(yīng)的材料參數(shù)，并給定轉(zhuǎn)速以及位移約束條件，通過靜力分析得到離心力引起的預(yù)應(yīng)力；然后，采用振型歸一化的模態(tài)分析，并計入離心力引起的預(yù)應(yīng)力，獲得葉片的固有振動特性，提取葉片振動固有頻率以及模態(tài)振型，得到固體域葉片表面的節(jié)點振幅；步驟三建立計算流體力學(xué)模型；首先，通過給定的葉輪機械葉型數(shù)據(jù)，在計算流體力學(xué)仿真工具CFX的TurboGrid模塊中建立葉柵的單通道流場模型，其中葉片表面附近為O型網(wǎng)格，其他區(qū)域為H型網(wǎng)格；然后，將建立的單扇區(qū)流場模型導(dǎo)入到ICEM CFD中，通過繞旋轉(zhuǎn)軸旋轉(zhuǎn)、復(fù)制(N-1)份網(wǎng)格得到全環(huán)葉柵模型，并輸出流體域網(wǎng)格點和實體單元信息以及流體域各葉片表面網(wǎng)格點和表面單元信息；步驟四采用數(shù)據(jù)傳遞方法獲得插值到流體域葉片表面的振動位移；首先，基于恒定位置矢量差方法，采用有限單元法中局部坐標形式的形函數(shù)進行三維線性插值，將固體域中模態(tài)分析得到的第P階模態(tài)對應(yīng)的節(jié)點振動位移插值到流固耦合交界面上的表面單元網(wǎng)格點上，獲得流體域中葉片表面網(wǎng)格點的振幅；然后，通過流體域中各葉片扇區(qū)之間的角度關(guān)系，即可依次逐個得到全環(huán)各葉片表面網(wǎng)格點的振幅；步驟五采用動網(wǎng)格技術(shù)獲得每一個葉片扇區(qū)的可動域網(wǎng)格點振幅；根據(jù)流體域中網(wǎng)格點編號順序，對每一個葉片沿O型域網(wǎng)格的法向逐層遞推搜索，得到與該葉片表面網(wǎng)格點編號順序相同的各層網(wǎng)格點，并按照各層網(wǎng)格點與葉片表面網(wǎng)格點之間法向距離的比值大小來線性插值，得到該葉片可動域的網(wǎng)格點振幅；步驟六獲得非定常計算流體力學(xué)即CFD分析所需的考慮葉間相位角的網(wǎng)格文件；首先，將葉片一個振動周期分為K步，葉片以簡諧規(guī)律振動，給定各葉片相對于第一個葉片的初始相位角，通過初始坐標以及各時間步的振動位移，可以得到一個周期內(nèi)每一時間步上各葉片扇區(qū)可動域的坐標值；然后，通過固定域坐標值和每一時間步上各葉片扇區(qū)可動域的坐標值，可以得到一個周期內(nèi)每一時間步的流體網(wǎng)格點坐標組成的網(wǎng)格文件；步驟七采用CFX中Junction Box模塊進行非定常分析；流場分析通過計算流體力學(xué)仿真工具CFX實現(xiàn)；首先，將步驟三建立的流場模型的網(wǎng)格點坐標值換成步驟六中獲得的第一個時間步的網(wǎng)格點坐標值，對流場模型加載邊界條件和初始條件；通道入口給定總溫和總壓，出口給定平均靜壓，輪緣、輪轂給定無滑移、光滑壁面邊界條件，葉片表面給定動網(wǎng)格壁面邊界條件；以葉片表面給定無滑移、光滑壁面邊界條件的定常解作為非定常分析的初始條件；然后，指定葉片的振動模態(tài)以及振動幅值，其運動周期為固有頻率的倒數(shù)，在每一個振動時間步上求解采用k-ε湍流模型封閉的Reynolds平均Navier-Stokes方程，并寫出每一個求解時間步的瞬時結(jié)果文件；步驟八獲得各葉間相位角對應(yīng)的模態(tài)氣動阻尼比，預(yù)測振蕩葉柵的氣彈穩(wěn)定性；首先，在不同葉間相位角情況下，根據(jù)一個穩(wěn)定周期內(nèi)各時間步的結(jié)果文件，基于等效粘性阻尼推導(dǎo)的模態(tài)氣動阻尼比獲得葉片一個振動周期的氣動功和模態(tài)氣動阻尼比；其中氣動功是流場各網(wǎng)格點氣動力在一個穩(wěn)定振動周期內(nèi)所做的總功；基于等效粘性阻尼推導(dǎo)的模態(tài)氣動阻尼比公式為其中W為非定常氣動力在一個周期內(nèi)所做的氣動功，qaero為葉片的振動幅值；然后，通過不同葉間相位角情況下的模態(tài)氣動阻尼比，判斷某階模態(tài)下某一工況的氣彈穩(wěn)定性，在各葉間相位角下，只要有某一個葉間相位角對應(yīng)的模態(tài)氣動阻尼比為負，則該階模態(tài)下，該工況氣彈失穩(wěn)。FSA00000166970900021.tif
2.根據(jù)權(quán)利要求1所述的葉輪機械計入葉間相位角的氣彈穩(wěn)定性流固耦合預(yù)測方法， 其特征在于步驟四中所述的數(shù)據(jù)傳遞方法的具體實施步驟如下步驟一獲得與第一個葉片表面的流體網(wǎng)格點F距離最近的固體節(jié)點S ； 對流體域中第一個葉片表面的任意網(wǎng)格點F，在固體域中找到與其距離最近的固體表 面節(jié)點S，與節(jié)點S相關(guān)的單元至多有4個，分別為Sp S2、S3和S4,并從固體域中獲得各單 元對應(yīng)的結(jié)構(gòu)單元的表面編號；步驟二 獲得流體網(wǎng)格點F在不同固體單元中的局部坐標值；流體網(wǎng)格點F的坐標為(Xq，y0, z0)，固體單元八個節(jié)點的坐標分別為(Xi，Zi)，(i =·1，...，8)，每個節(jié)點的形函數(shù)在局部坐標系下可以表示為隊(ε，η, ζ)，(i = 1，...，8)，則通過流體網(wǎng)格點坐標和固體單元節(jié)點坐標之間的關(guān)系 ‘ 8 可以得到流體網(wǎng)格點F在與節(jié)點S相關(guān)的不同單元的局部坐標值； 其中局部坐標系下的形函數(shù)可以表示為 Ν( ε，η，ζ) = (1士 ε ) (1士 η) (1士 ζ) 步驟三獲得流體網(wǎng)格點F對應(yīng)固體域的投影單元M ；分別對比不同單元對應(yīng)的局部坐標值，找出各局部坐標值的平方和最小值對應(yīng)的單元 M以及該單元的表面編號，即為流體網(wǎng)格點F在固體域中的投影單元和單元面，將網(wǎng)格點F 投影到單元M的表面得到F'，由于流體域網(wǎng)格和固體域網(wǎng)格對原始幾何構(gòu)型的不同精度 導(dǎo)致網(wǎng)格點F和投影點F'之間存在坐標差，采用恒定位置矢量差法，使流體網(wǎng)格點F的函 數(shù)值始終等于其投影點F'的函數(shù)值；步驟四插值得到流體網(wǎng)格點F的位移；首先，獲得流體網(wǎng)格點F在投影單元M上的局部坐標值，并進一步得到單元M的八個形 函數(shù) NiG，η, ζ) = (1士 ε) (1 士 η) (1士 ζ)，(i = 1，···，8)；然后，通過固體域中實體單元的八個節(jié)點振動位移(Ui,Vi，Wi) (i = 1，...，8)可以插值獲得流固耦合交界面上流體網(wǎng)格點F的位移(Utl，v0, Wtl)，具體的位移插值為步驟五獲得其他葉片的表面網(wǎng)格點振動位移；首先，獲得與第一個葉片表面網(wǎng)格點編號相對應(yīng)的其他各葉片的表面網(wǎng)格點編號； 在旋轉(zhuǎn)坐標軸為Z的情況下，各葉片逆時針編號，第i個葉片表面網(wǎng)格點坐標(Xi，Yi, Zi)與第一個葉片表面網(wǎng)格點坐標(Xl，Y1, Z1)之間滿足如下關(guān)系 其中當?shù)谝粋€葉片在第一象限時， 當?shù)谝粋€葉片在其他象限時，按照對應(yīng)的角度關(guān)系給出θ工； 其中， N為全環(huán)葉片個數(shù)；然后，將第一個葉片表面網(wǎng)格點的振幅分別賦值給其他各葉片的對應(yīng)網(wǎng)格點，獲得全 環(huán)所有葉片表面網(wǎng)格點的振幅；
3.根據(jù)權(quán)利要求1所述的葉輪機械計入葉間相位角的氣彈穩(wěn)定性流固耦合預(yù)測方法， 其特征在于步驟六中所述生成考慮葉間相位角的網(wǎng)格文件實施步驟如下步驟一將葉片的一個振動周期分為K步，其中一個周期對應(yīng)的頻率為葉片在某階模 態(tài)下的固有頻率ω，給定第一個葉片正弦振動的初始相位角為零，則第i個葉片在第k個時 間步的相位角為 其中，1為節(jié)徑數(shù)，N為全環(huán)葉片個數(shù)； 則該葉片扇區(qū)的每一個網(wǎng)格點在該時間步的振動位移為 其中，Ui, Vi, Wi分別為第i個葉片扇區(qū)的某一個網(wǎng)格點在X，y，ζ三個方向的振幅； 步驟二 通過動網(wǎng)格技術(shù)搜索獲得第i個葉片對應(yīng)可動域的網(wǎng)格點編號，并從全域網(wǎng) 格點坐標值可以得到該葉片扇區(qū)的可動域網(wǎng)格點初始坐標值(Xi，Ii, Zi)，而該葉片扇區(qū)的 每一個網(wǎng)格點在第k個時間步的X，y，ζ三個方向振動位移分別為uik，vik, wik,則該葉片扇 區(qū)的每一個網(wǎng)格點在第k個時間步的坐標值為Xik ~ Xi + Uik -yik =兄 + v,k = zi + ^lk步驟三通過固定域坐標值和每一時間步上各葉片扇區(qū)可動域的坐標值，可以得到一 個周期內(nèi)每一時間步的流體網(wǎng)格點坐標，并將每一個時間步的網(wǎng)格點坐標值輸出為網(wǎng)格文 件。
全文摘要
葉輪機械計入葉間相位角的氣彈穩(wěn)定性流固耦合預(yù)測方法，采用能量法，通過不同葉間相位角的模態(tài)氣動阻尼比預(yù)測葉輪機械的氣彈穩(wěn)定性。該方法對建立的單扇區(qū)有限元模型采用含預(yù)應(yīng)力的模態(tài)分析，通過數(shù)據(jù)傳遞方法實現(xiàn)流固耦合交界面上的振動位移傳遞，并采用動網(wǎng)格技術(shù)得到考慮葉間相位角的網(wǎng)格文件，用于CFX中Junction Box模塊的非定常計算流體力學(xué)分析，進而獲得各葉間相位角的模態(tài)氣動阻尼比，用于預(yù)測氣彈穩(wěn)定性。該基于流固耦合方法，充分利用了計算結(jié)構(gòu)動力學(xué)和計算流體動力學(xué)的前沿技術(shù)，在每個子系統(tǒng)內(nèi)保證了計算精度，而固定域和可動域的劃分提高了計算效率。它在葉輪機械模擬技術(shù)領(lǐng)域具有良好的實用價值和廣闊的應(yīng)用前景。
文檔編號G06F17/50GK101882177SQ20101020979
公開日2010年11月10日 申請日期2010年6月18日 優(yōu)先權(quán)日2010年6月18日
發(fā)明者張小偉, 王延榮 申請人:北京航空航天大學(xué)