一種基于時域雙向迭代的葉輪機葉片顫振應力預測方法

專利名稱：一種基于時域雙向迭代的葉輪機葉片顫振應力預測方法
技術領域：
本發明涉及一種計算機輔助分析工具領域的方法，具體是涉及一種基于時域雙向 迭代的葉輪機葉片顫振應力預測方法。屬于葉輪機械模擬技術領域。
背景技術：
目前，在葉輪機械中，存在著振動葉片與周圍流體之間復雜的相互作用，彈性葉片 通過其變形影響流場，流場的擾動反之改變葉片的形狀。在這類流固耦合現象中，以顫振對 結構可靠性的影響最大。顫振屬于一種自激振動，葉片表面的非定常氣動力來源于葉片本 身的運動，非定常氣動力和葉片的非定常運動互為因果、循環遞增，使得短時間內葉片振幅 急劇增大，甚至造成葉片斷裂和機器損壞。由于進行顫振試驗難度較高、成本巨大且具有相 當的危險性，所以隨著計算機及相關學科的發展，數值模擬方法因其成本低、效率高，成為 顫振分析的重要工具。在葉輪機械的顫振數值模擬中，通常使用能量法和特征值法，近年來也少量開始 使用時間推進法。現有的顫振數值計算方法至少存在以下四個缺點首先，傳統的能量法和特征值法都假定葉片以給定形式作運動，將流場和葉片運 動互相孤立，忽略了彈性葉片本身變形對流場的影響，同時也極少考慮離心力對葉片變形 和流場的影響，無法實現氣動_結構一體化計算，存在較大的計算誤差。其次，由于流體和結構都存在強烈的非線性，這兩種領域的非線性交替疊加，形成 了更加復雜的耦合非線性。現有的顫振預測方法都引入了大量線性化假設，甚至只采用結 構向流體的單向傳遞，弱化甚至消除了耦合非線性，無法準確地描述非線性現象，降低了模 擬精度。再次，現有的方法不能展示整個顫振發展的時間歷程，也無法得到顫振應力，而顫 振應力恰恰是顫振現象中非常重要的指標。最后，以往的基于時間推進的葉輪機顫振數值模擬方法都采用固定時間步長，無 法在計算中對步長進行有效的調整，造成計算速度和效率較低，甚至會導致收斂困難；而且 大多采用商用的結構有限元計算軟件，在計算中需要反復進行人工干預，不能在軟件內部 針對流固耦合計算的特性對整個計算流程和算法進行優化，增加了操作復雜度，降低了計 算效率，使得現有的基于時間推進的預測方法計算費用較高，經濟型較差。

發明內容
1、目的本發明的目的是克服常規的葉輪機顫振數值仿真方法的不足，提供一種 基于時域雙向迭代的葉輪機葉片顫振應力預測方法，使之能夠準確描述流體與葉片之間的 相互作用和耦合非線性現象，并對葉片的顫振應力進行預測。2、技術方案本發明一種基于時域雙向迭代的葉輪機葉片顫振應力預測方法，它 是在計算機上采用數值模擬的方法對葉輪機葉片的顫振應力進行預測。該方法具體步驟如 下
步驟一在計算機中設定以下六個模塊結構計算模塊、流體計算模塊、數據轉換 模塊、顫振應力輸出模塊、初值計算模塊和雙向迭代模塊。步驟二 調用初值計算模塊通過非線性迭代獲得葉片的初始靜變形和穩態定常流 場作為以后雙向迭代計算的初值，用以加快整個系統的收斂。步驟三調用雙向迭代模塊，在時間上推進由葉片和周圍流場組成的流固耦合系 統。步驟四通過顫振應力輸出模塊，將結構的瞬態位移響應轉化為瞬態顫振應力，并 輸出為文件，供后處理軟件讀取顯示。其中，步驟一所述的六個模塊的具體結構如下模塊一結構計算模塊結構計算模塊的作用是求解三維結構動力方程，獲得葉片的瞬態變形。該模塊包 含如下步驟步驟1. 1 將葉片表面壓力和結構有限元節點單元信息放到輸入文件里，并在輸 入文件中對葉片的工況進行設置。步驟1. 2 調用Fortran語言編寫的結構有限元計算程序求解三維結構動力方程
Md + Cd + (K-Kc + S)a = f + fa + fc(1)其中M，C和K分別是總體質量矩陣、阻尼矩陣和剛度矩陣。K。是旋轉軟化矩陣。S 是應力剛化矩陣(或稱幾何剛化矩陣)。f, fa和f。是節點外力列陣、氣動力列陣和離心力 列陣。a是待求的有限元節點位移列向量。步驟1. 3 計算結束后，輸出葉片表面有限元節點的位移。模塊二 流體計算模塊在流體計算模塊中，提供了各種流體計算軟件的接口，可以后臺調用流體計算軟 件獲得每個時間步的三維非定常流場。該模塊中根據用戶需要調用各種流體動力學計算軟 件，如Fluent，CFX等，在計算結束后輸出葉片表面壓力。模塊三數據轉換模塊數據轉換模塊包括三個子模塊，一是節點配對子模塊，作用是找到耦合界面上每 個流體節點所對應的結構單元；二是載荷傳遞子模塊，該子模塊通過三維形函數插值，將葉 片表面的定常或非定常壓力轉換為葉片結構有限元模型上的節點力；三是變形傳遞子模 塊，該子模塊將葉片的變形轉換為葉片周圍流體域的邊界運動。下面分別對這三個子模塊 進行介紹節點配對子模塊包含下列步驟步驟3. 1. 1 對耦合界面上每個流體節點f，遍歷所有耦合界面上的結構節點，找 到距離f最近的結構節點S。步驟3. 1. 2 找到所有包含結構節點s的結構單元e2,…，em。步驟3. 1.3 對每個結構單元力(1 = 1，2，…，m)，用擬牛頓法求出f在該單元內
的局部坐標(ri; Si,、)，并求出f與該單元中心的相對距離式=擬+sf+tf。步驟3. 1. 4 找出與f相對距離最小的單元et，即dt = π η^^，.··，《)，該單元 即為流體節點f所對應的結構單元。
變形傳遞子模塊包含下列步驟步驟3. 2. 1 設流體節點f在其所對應結構單元e上的投影點為f’，f’在整體 坐標系下的坐標為(xf，，yf，，zf，)，忽略流體節點與其投影點之間的距離，近似有xf, ^xf, yf, ^ yf，Zf, ^ zf，投影點f’在結構單元e中的局部坐標(r，s, t)通過求解下列方程組得 到 其中η為結構單元e中節點的數目，Xi，yi，Zi(i = 1,2,…，8)為結構單元e中8 個節點在整體坐標系下的坐標。投影點處的形函數NiG = 1,2,…，8)為N1 = (1-s) (1-t) (l-r)/8, N2 = (1+s) (l_t) (1-r) /8N3 = (1+s) (1+t) (1-r)/8, N4 = (1-s) (1+t) (1-r) /8
(3)N5 = (1-s) (1-t) (l+r)/8, N6 = (1+s) (1-t) (1+r) /8N7 = (1+s) (1+t) (1+r)/8, N8 = (1-s) (1+t) (1+r) /8步驟3. 2. 2 在求出投影點f’在結構單元e中的局部坐標(r，s, t)后，如果投影 點在面r= 1上，忽略流體節點與其對應投影點之間的距離，強制令r= 1 ；如果投影點在面 s = 1上，忽略流體節點與其對應投影點之間的距離，強制令s = 1 ；如果投影點在面t = 1 上，忽略流體節點與其對應投影點之間的距離，強制令t = 1。步驟3. 2. 3 假設投影點f’的位移為uf,，Ui為結構有限元節點i的位移，通過下 列公式求解流體節點f的位移Uf 載荷傳遞子模塊包含下列步驟步驟3. 3. 1 設流體表面網格c的中心點為f。，計算流體壓力沿網格c的積分 Ffc = IlPdS^Pfc-Ac，其中Ω為網格c的表面，ρ為耦合界面上的流體壓力，&為中心點f。 上的壓力，Α。為網格c的表面積。步驟3. 3. 2 設流體表面網格中心點f。對應的結構單元為e，通過計算尺-況/^ (i =1,2,…，8)可以得到結構單元e上8個節點的節點力，其中Ni為&在結構單元e上投 影點處的形函數。步驟3. 3. 3 將結構單元e上8個節點的節點力集合入結構有限元節點載荷列陣。模塊四顫振應力輸出模塊顫振應力輸出模塊通過葉片的瞬態位移響應計算出葉片的瞬態顫振應力，并輸出 各種供顫振分析使用的后處理文件。該模塊包括如下步驟
步驟4. 1 通過葉片的瞬態位移求出葉片的瞬態顫振應力，即O=De (5)其中ε為應變列陣，σ為應力列陣，D為彈性矩陣。步驟4. 2 將葉片隨時間變化的位移、顫振應力和壓力輸出為可被Origin軟件讀 取的格式。步驟4. 3 將轉子隨時間變化的質量流量和壓比輸出為可被Origin軟件讀取的格式。步驟4. 4 將葉片瞬間的位移云圖和顫振應力云圖輸出為可被Tecplot軟件讀取 的格式。步驟4. 5 如果整個計算結束，輸出該級轉子葉片在壓氣機特性圖上的顫振邊界。模塊五初值計算模塊初值計算模塊通過反復的非線性迭代獲得葉片的初始靜變形和穩態定常流場，作 為以后雙向迭代計算的初值，用以加快整個系統的收斂。該模塊包含如下步驟步驟5. 1 首先通過三維模型設計軟件UG建立葉片的三維實體模型；然后將實體 模型導入網格劃分軟件HyperMesh中，得到葉片的結構有限元網格；將實體模型導入網格 劃分軟件TurboGrid中，得到葉輪機內流場的流體計算網格。步驟5. 2 調用結構計算模塊得到葉片在離心力作用下的初始位移其中Ktl 是結構的小位移剛度矩陣，Fc為葉片在轉動中的離心力列陣。步驟5. 3 調用數據轉換模塊將葉片變形轉化為流體域的邊界運動。步驟5. 4 通過流體計算模塊調用流體動力學計算軟件，進行網格更新，然后計算 新的定常流場。步驟5. 5 調用數據轉換模塊將葉片表面的定常壓力轉換為結構有限元模型上的 節點力。步驟5.6 調用結構計算模塊獲得葉片的非線性變形，并對結構的剛度矩陣進行 更新。計算葉片非線性變形的具體步驟如下步驟5. 6. 1 根據步驟5. 2求得的初始位移％，求得應力剛化矩陣Ks、大變形剛度 矩陣Kl和反力Fm。步驟5. 6. 2 根據(Kq-Kc+Kq +KL) (Ia1 = Fc-Fnr 求出增量位移(! 。步驟5. 6. 3 計算總位移 B1 = a0+dalo步驟5. 6.4:如果求得的an滿足收斂條件I |dan| |2< ε ||an||2，其中ε取0.001， 則非線性迭代收斂，計算結束，指向步驟5. 7 ；否則用ai代替％，前往步驟5. 3。步驟5. 7 求得新的結構剛度矩陣K = Kci-KfKfKp常規的有限元計算軟件無法 同時考慮幾何非線性和應力剛化的影響，本發明通過步驟5. 7即時對結構剛度矩陣進行修 正，從而克服了這一缺陷。此時的葉片靜變形和定常流場就作為以后雙向迭代計算的初始 值。模塊六雙向迭代模塊雙向迭代模塊反復交替調用結構求解器和流體求解器，同時采用數據轉換模塊在 結構求解器和流體求解器之間傳遞邊界數據，通過流體與結構雙向的交互作用，在時間上 推進整個由葉片及其周圍流場組成的流固耦合系統。雙向迭代模塊包含下列步驟
步驟6. 1 通過數據轉換系統將葉片的變形轉換為流體計算域的邊界運動。步驟6. 2 通過流體計算模塊調用流體動力學計算軟件，進行網格更新，然后計算 新的非定常流場。步驟6. 3 通過數據轉換模塊將葉片表面的非定常壓力轉換為結構有限元模型上 的節點力。步驟6. 4 通過結構計算模塊求得新的葉片變形。步驟6. 5 對雙向迭代的計算時間步進行調整。本發明提供了一種新的流固耦合 自適應時間步技術，通過當前的收斂速度控制下一步計算的時間步長，隨時對后續計算的 時間步長進行調整，從而有效地縮短收斂時間，提高計算效率。該流固耦合自適應時間步技 術包括下列步驟步驟6. 5. 1 計算在[(奶)2+(u,-u/ )2+…+Uru1/ )2]/("/+"/+···+"/)，其中
Ui為當前時間步計算的節點i在某個自由度的位移，Ui’為上個時間步計算的節點i在某個 自由度的位移，η為有限元節點的數目。步驟6. 5. 2 取各個自由度的d的平均值為Cltl，用以表征當前時間步的收斂速度。步驟6.5.3 用e代表收斂速度的閾值，用t代表當前的時間步長，如果Cltl <0.5e， 則新的時間步長tnew = 2t ；如果dQ>2e，則新的時間步長tnew = 0. 5t ；如果0. 5e<d0<2e, 則時間步長不變，tnew = t。步驟6. 5. 4 在獲得流固耦合計算的物理時間步長后，對結構動力方程瞬態求 解的時間步長和流體非定常求解的時間步長作相應調整。步驟6. 6 通過計算所得的葉片瞬態位移響應，觀察顫振的發展進程。如果響應曲 線衰減，則不會發生顫振，計算結束；如果響應曲線增長，發生顫振，指向步驟6. 1 ；如果發 生顫振且計算達到預先指定時間點，計算結束，指向步驟6. 7輸出顫振應力。步驟6. 7 通過顫振應力輸出模塊，將結構的瞬態位移響應轉化為瞬態顫振應力， 并輸出為文件，供后處理軟件讀取顯示。3、本發明的優點在于(1)在計算開始前，加入一個初值計算的過程。初值計算中引入了非設計狀態時葉 片離心力對結構剛度矩陣和流場的影響，通過反復的非線性迭代獲得初值，從而加快了整 個計算的收斂速度；同時，常規的有限元計算軟件無法同時考慮幾何非線性和應力剛化的 影響，本發明通過非線性迭代隨時對結構剛度矩陣進行修改，克服了這一缺陷。(2)常規方法都假定葉片以指定形式作運動，將流場和葉片運動互相孤立，忽略了 彈性葉片本身變形對流場的影響，本發明將流場和葉片運動結合起來，成為一個整體的流 固耦合系統，采用一個雙向迭代的流程，實現了氣動_結構一體化計算，減小了計算誤差。(3)通過雙向迭代可以準確地描述結構非線性、流體非線性和流固耦合非線性，避 免了單向傳遞或大量的線性化假設對模擬結果精度的影響。(4)本發明通過時間推進提供整個顫振發展的時間歷程圖，并得到各個時間點的 顫振應力。一方面可以驗證實際中顫振現象的發生發展過程，為顫振機理的闡釋提供依據； 另一方面可以為顫振試驗提供指導和依據，預測達到預期顫振應力所需的時間，指導試驗 及時停車，防止由于葉片損毀甚至飛脫碰撞造成嚴重的試驗事故。(5)以往的基于時域的葉輪機顫振數值模擬方法都采用固定時間步長，不僅造成收斂困難，而且計算速度和效率較低，本發明中發展了一種流固耦合自適應時間步技術，通 過當前的收斂速度控制下一步計算的時間步長，隨時對后續計算的時間步長進行調整，從 而有效地縮短收斂時間，提高計算效率。同時，采用全自動化的操作流程，通過簡單的輸入 文件進行控制，無需在計算過程中進行大量人工干預，避免了繁瑣的操作。另外，通過自編 的結構有限元求解程序，在程序內部針對流固耦合計算的特性對整個計算流程和算法進行 優化，大大提高了計算效率和經濟性。


圖1是本發明的實現流程圖；圖2是本發明的系統模塊圖；圖3是本發明的初值計算模塊的流程圖；圖4是本發明所用某航空發動機壓氣機葉片及流場模型；圖5是本發明的雙向迭代與現有方法的區別示意圖；圖6是本發明的流固耦合自適應時間步技術的示意圖；圖7是本發明葉根某點顫振應力的時間歷程圖；圖8是本發明預測壓氣機顫振特性圖；圖9是本發明葉根某點靜壓的時間歷程圖；圖10是本發明葉片瞬間顫振應力分布圖。
具體實施例方式下面將結合附圖和實施例對本發明作進一步的詳細說明。本發明是基于時域雙向 迭代的葉輪機葉片顫振應力預測方法，方法的實現流程如圖1所示。步驟一在計算機中設定以下六個模塊結構計算模塊、流體計算模塊、數據轉換 模塊、顫振應力輸出模塊、初值計算模塊和雙向迭代模塊。步驟二 調用初值計算模塊通過非線性迭代獲得葉片的初始靜變形和穩態定常流 場作為以后雙向迭代計算的初值，用以加快整個系統的收斂。步驟三調用雙向迭代模塊，在時間上推進由葉片和周圍流場組成的流固耦合系 統。步驟四通過顫振應力輸出模塊，將結構的瞬態位移響應轉化為瞬態顫振應力，并 輸出為文件，供后處理軟件讀取顯示。步驟一中提到的六個模塊包括結構計算模塊、流體計算模塊、數據轉換模塊、顫振 應力輸出模塊、初值計算模塊和雙向迭代模塊，如圖2所示。這六個模塊并非互相孤立，而 是互相聯系，在使用其中一個模塊時，可能會在該模塊中調用其他模塊。六個模塊的結構如 下模塊一結構計算模塊結構計算模塊的作用是求解三維結構動力方程，獲得葉片的瞬態變形。該模塊按 照下列步驟進行運算求解步驟1. 1 將葉片表面壓力和結構有限元節點單元信息放到輸入文件里，并在輸 入文件中對葉片的工況進行設置。
步驟1. 2 調用Fortran語言編寫的結構有限元計算程序求解三維結構動力方程 獲得葉片的瞬態變形。步驟1. 3 計算結束后，輸出葉片表面有限元節點的位移。模塊二 流體計算模塊在流體計算模塊中，提供了各種流體計算軟件的接口，可以后臺調用流體計算軟 件獲得每個時間步的三維非定常流場。該模塊包含如下步驟步驟2. 1 按照需要選擇求解使用的流體動力學計算軟件，如Fluent，CFX等。步驟2. 2 在后臺調用所選擇的流體動力學計算軟件進行求解。步驟2.3 計算結束后，輸出葉片表面的壓力。模塊三數據轉換模塊數據轉換模塊包括三個子模塊，一是節點配對子模塊，作用是找到耦合界面上每 個流體節點所對應的結構單元；二是載荷傳遞子模塊，該子模塊通過三維形函數插值，將葉 片表面的定常或非定常壓力轉換為葉片結構有限元模型上的節點力；三是變形傳遞子模 塊，該子模塊將葉片的變形轉換為葉片周圍流體域的邊界運動。模塊四顫振應力輸出模塊顫振應力輸出模塊通過葉片的瞬態位移響應計算出葉片的瞬態顫振應力，并輸出 各種供顫振分析使用的后處理文件。該模塊包括如下步驟步驟4. 1 通過葉片的瞬態位移求出葉片的瞬態顫振應力。步驟4. 2 將葉片隨時間變化的位移、顫振應力和壓力輸出為可被Origin軟件讀 取的格式。步驟4. 3 將轉子隨時間變化的質量流量和壓比輸出為可被Origin軟件讀取的格式。步驟4.4 將葉片瞬間的位移云圖和顫振應力云圖輸出為可被Tecplot軟件讀取 的格式。步驟4. 5 如果整個計算結束，輸出該級轉子葉片在壓氣機特性圖上的顫振邊界。模塊五初值計算模塊初值計算模塊通過反復的非線性迭代獲得葉片的初始靜變形和穩態定常流場，作 為以后雙向迭代計算的初值，用以加快整個系統的收斂。該模塊包含如下步驟步驟5. 1 首先通過三維模型設計軟件UG建立葉片的三維實體模型；然后將實體 模型導入網格劃分軟件HyperMesh中，得到葉片的結構有限元網格；將實體模型導入網格 劃分軟件TurboGrid中，得到葉輪機內流場的流體計算網格。步驟5. 2 調用結構計算模塊得到葉片在離心力作用下的初始位移列陣a^AT尺， 其中Ktl是結構的小位移剛度矩陣，Fc為葉片在轉動中的離心力列陣。步驟5. 3 調用數據轉換模塊將葉片變形轉化為流體域的邊界運動。步驟5. 4 通過流體計算模塊調用流體動力學計算軟件，進行網格更新，然后計算 新的定常流場。步驟5. 5 調用數據轉換模塊將葉片表面的定常壓力轉換為結構有限元模型上的 節點力。步驟5.6 調用結構計算模塊獲得葉片的非線性變形，并對結構的剛度矩陣進行更新。如果非線性迭代收斂，計算結束，指向步驟5. 7;否則用新的葉片位移列陣取代舊的 葉片位移列陣，前往步驟5. 3。步驟5. 7 求得新的結構剛度矩陣K = Ktl-KfKdKp常規的有限元計算軟件無法 同時考慮幾何非線性和應力剛化的影響，本發明通過步驟5. 7即時對結構剛度矩陣進行修 正，從而克服了這一缺陷。此時的葉片靜變形和定常流場就作為以后雙向迭代計算的初始值。模塊六雙向迭代模塊雙向迭代模塊反復交替調用結構求解器和流體求解器，同時采用數據轉換模塊在 結構求解器和流體求解器之間傳遞邊界數據，通過流體與結構雙向的交互作用，在時間上 推進整個由葉片及其周圍流場組成的流固耦合系統。雙向迭代模塊包含下列步驟步驟6. 1 通過數據轉換系統將葉片的變形轉換為流體計算域的邊界運動。步驟6. 2 通過流體計算模塊調用流體動力學計算軟件，進行網格更新，然后計算 新的非定常流場。步驟6. 3 通過數據轉換模塊將葉片表面的非定常壓力轉換為結構有限元模型上 的節點力。步驟6. 4 通過結構計算模塊求得新的葉片變形。步驟6. 5 對雙向迭代的計算時間步進行調整。本發明提供了一種新的流固耦合 自適應時間步技術，通過當前的收斂速度控制下一步計算的時間步長，隨時對后續計算的 時間步長進行調整，從而有效地縮短收斂時間，提高計算效率。該流固耦合自適應時間步技 術包括下列步驟步驟6. 5. 1 計算oK("廠V+(U,-U/ )、..+ (" -"/ )2]/("/+"/+·..+"/)，其中
Ui為當前時間步計算的節點i在某個自由度的位移，Ui’為上個時間步計算的節點i在某個 自由度的位移，η為有限元節點的數目。步驟6. 5. 2 取各個自由度的d的平均值為Cltl，用以表征當前時間步的收斂速度。步驟6. 5. 3 用e代表收斂速度的閾值，用t代表當前的時間步長，如果Cltl <0. 5e, 則新的時間步長tnew = 2t ；如果dQ>2e，則新的時間步長tnew = 0. 5t ；如果0. 5e<d0<2e, 則時間步長不變，tnew = t。步驟6. 5. 4 在獲得流固耦合計算的物理時間步長tnew后，對結構動力方程瞬態求 解的時間步長和流體非定常求解的時間步長作相應調整。步驟6. 6 通過計算所得的葉片瞬態位移響應，觀察顫振的發展進程。如果響應曲 線衰減，則不會發生顫振，計算結束；如果響應曲線增長，發生顫振，指向步驟6. 1 ；如果發 生顫振且計算達到預先指定時間點，計算結束，指向步驟6. 7輸出顫振應力。步驟6. 7 通過顫振應力輸出模塊，將結構的瞬態位移響應轉化為瞬態顫振應力， 并輸出為文件，供后處理軟件讀取顯示。現采用某航空發動機壓氣機葉片為實例，對本發明作進一步的詳細說明。本例 中使用的設備為一臺Intel內核的微型計算機，主頻為2. 8GHz,內存為2G，操作系統為 Microsoft Windows XP Professional, Service Pack 2。步驟1 初值計算調用初值計算模塊以加快收斂，其流程如圖3所示，具體步驟如下
步驟1. 1 首先通過三維模型設計軟件UG建立葉片的三維實體模型；然后將實體 模型導入網格劃分軟件HyperMesh中，得到葉片的結構有限元網格；將實體模型導入網格 劃分軟件TurboGrid中，得到葉輪機內流場的流體計算網格；葉片的結構有限元網格和流 體網格見圖4，為清晰起見，流體只顯示了輪轂和葉片表面的網格；步驟1. 2 通過結構計算模塊調用Fortran語言編寫的結構有限元分析程序，得到 葉片在離心力作用下的變形，將葉片表面有限元節點的位移寫入文件bladenode. txt ；步驟1. 3 將文件bladenode. txt輸入數據轉換模塊，將葉片變形轉化為流體域的 邊界運動，生成葉片表面流體網格點的位移，寫入文件blad印oint. txt ；步驟1.4:通過流體計算模塊調用流體動力學計算軟件Fluent，讀入文件 bladepoint. txt，進行動網格計算，獲得更新后的流場網格，并計算新的定常流場，將葉片 表面流體網格點的壓力寫入文件blad印res. txt ；步驟1. 5 將文件blacbpres. txt輸入數據轉換模塊，將葉片表面的定常壓力轉換 為結構有限元模型上的節點力，寫入文件bladeforce. txt ；步驟1.6 通過結構計算模塊調用由標準Fortran語言編寫的結構有限元分析程 序，讀入節點力文件bladeforce. txt，進行非線性迭代獲得葉片的變形和新的結構剛度矩 陣，將葉片表面有限元節點的變形寫入文件bladenode. txt，將更新后的結構剛度矩陣寫入 t^f牛 restartek. dat ；步驟1. 7 如果非線性迭代收斂，計算結束，所得葉片靜變形和定常流場就作為以 后雙向迭代計算的初始值；否則如果非線性迭代不收斂，指向步驟1. 3。步驟2:雙向迭代調用雙向迭代模塊，在時間上推進整個流固耦合系統。與常規的單向迭代不同，本 發明通過流體與結構之間的雙向迭代實現了氣動-結構一體化計算，常規方法與本發明所 采用方法的區別見圖5。雙向迭代的具體步驟如下步驟2. 1 將文件bladenode. txt輸入數據轉換模塊，采用三維形函數插值方 法，將葉片的變形轉換為流體計算域的邊界運動，將葉片表面流體網格點的位移寫入文件 bladepoint. txt 步驟2. 2:通過流體計算模塊調用流體動力學計算軟件Fluent，讀入文件 bladepoint. txt，進行動網格計算，獲得更新后的流場網格，并計算新的三維非定常流場， 將葉片表面流體網格點的壓力寫入文件blacbpres. txt ；步驟2. 3 將文件blacbpres. txt輸入數據轉換模塊，采用三維形函數插值方法， 將葉片表面的非定常壓力轉換為結構有限元模型上的節點力，寫入文件bladeforce. txt ；步驟2. 4 通過結構計算模塊調用Fortran語言編寫的結構有限元分析程序，讀入 節點力文件bladeforce. txt,通過隱式Newmark法求解三維結構動力方程，獲得葉片的瞬 態變形，并將葉片表面有限元節點的位移寫入文件bladenode. txt ；步驟2.5 在雙向迭代中，運用本發明提供的流固耦合自適應時間步技術，通過當 前的收斂速度控制下一步計算的時間步長，隨時對后續計算的時間步長進行調整，從而有 效地縮短收斂時間，提高計算效率。本發明采用的流固耦合自適應時間步技術的示意圖見 圖6。步驟2. 6 通過計算所得的葉片瞬態位移響應，觀察顫振的發展進程。如果響應曲線衰減，則不會發生顫振，計算結束；如果響應曲線增長，則顫振發生，指向步驟2. 1 ；如果 顫振發生且計算達到預先指定的時間點，計算結束，指向步驟3輸出顫振應力。以葉根某點 為例，在圖7中給出了顫振應力的時間歷程圖。由圖中可以看出，在極短的時間內葉片振幅 急劇增大，嚴重威脅到了發動機的運行安全，這符合以往觀察到的顫振的特征。由于本發明 可以展示整個顫振發展的時間歷程，并得到顫振應力，這突破了現有的顫振數值模擬中的 局限。其具體意義如下通過本發明提供的模擬顫振發展歷程圖，可以觀察整個顫振發生和 發展的細節，從而通過數值模擬方法驗證了實際中的顫振現象，為顫振機理的闡釋提供了 依據；通過顫振發展歷程圖和計算所得顫振應力，可以為顫振試驗提供指導和依據，預測達 到預期顫振應力所需的時間，指導試驗及時停車，防止由于葉片損毀甚至飛脫碰撞造成嚴 重的試驗事故。步驟3:結果輸出通過顫振應力輸出模塊，將結構的瞬態位移響應轉化為瞬態顫振應力，同時輸出 一系列供顫振分析使用的后處理文件。所有輸出的數據文件包括步驟3. 1 輸出該級轉子葉片在壓氣機特性圖上的顫振邊界，如圖8所示，同時在 圖中附上了試驗所得的顫振邊界以作比較。由圖中可以看出，數值預測的結果與試驗的結 果比較吻合；步驟3.2 輸出葉片節點位移、葉片節點顫振應力、葉片表面靜壓、壓氣機質量流 量、壓氣機壓比的時間歷程圖；以葉根某點靜壓為例，其時間歷程圖見圖9 ；步驟3.3 輸出每個時間步葉片表面的位移云圖和顫振應力云圖，并將所有時間 步的云圖連接起來形成動畫；以某時刻葉片表面顫振應力為例，其云圖見圖10。通過應力 云圖可以獲得顫振時葉片應力集中的危險位置，為后續分析提供依據。
權利要求
一種基于時域雙向迭代的葉輪機葉片顫振應力預測方法，其特征在于，該方法具體步驟如下步驟一在計算機中設定以下六個模塊結構計算模塊、流體計算模塊、數據轉換模塊、顫振應力輸出模塊、初值計算模塊和雙向迭代模塊；步驟二調用初值計算模塊通過非線性迭代獲得葉片的初始靜變形和穩態定常流場作為以后雙向迭代計算的初值，用以加快整個系統的收斂；步驟三調用雙向迭代模塊，在時間上推進由葉片和周圍流場組成的流固耦合系統；步驟四通過顫振應力輸出模塊，將結構的瞬態位移響應轉化為瞬態顫振應力，并輸出為文件，供后處理軟件讀取顯示。
2.根據權利要求1所述的一種基于時域雙向迭代的葉輪機葉片顫振應力預測方法，其 特征在于步驟一中所述的六個模塊的具體結構如下模塊一結構計算模塊結構計算模塊的作用是求解三維結構動力方程，獲得葉片的瞬態變形；該模塊包含如 下步驟步驟1. 1 將葉片表面壓力和結構有限元節點單元信息放到輸入文件里，并在輸入文 件中對葉片的工況進行設置；步驟1. 2 調用Fortran語言編寫的結構有限元計算程序求解三維結構動力方程 Md + Ca + (K-Kc + S)a = f + fa + fc⑴其中M，C和K分別是總體質量矩陣、阻尼矩陣和剛度矩陣；K。是旋轉軟化矩陣；S是應 力剛化矩陣即幾何剛化矩陣；f，fa和f。是節點外力列陣、氣動力列陣和離心力列陣，a是待 求的有限元節點位移列向量；步驟1. 3 計算結束后，輸出葉片表面有限元節點的位移； 模塊二 流體計算模塊在流體計算模塊中，提供了各種流體計算軟件的接口，可以后臺調用流體計算軟件 獲得每個時間步的三維非定常流場；該模塊中根據用戶需要調用流體動力學計算軟件 Fluent, CFX，在計算結束后輸出葉片表面壓力； 模塊三數據轉換模塊數據轉換模塊包括三個子模塊，一是節點配對子模塊，作用是找到耦合界面上每個流 體節點所對應的結構單元；二是載荷傳遞子模塊，該子模塊通過三維形函數插值，將葉片表 面的定常或非定常壓力轉換為葉片結構有限元模型上的節點力；三是變形傳遞子模塊，該 子模塊將葉片的變形轉換為葉片周圍流體域的邊界運動； 節點配對子模塊包含下列步驟步驟3. 1. 1 對耦合界面上每個流體節點f，遍歷所有耦合界面上的結構節點，找到距 離f最近的結構節點S;步驟3. 1. 2 找到所有包含結構節點s的結構單元e2,…，effl ；步驟3. 1.3 對每個結構單元&(1 = 1,2,…，m)，用擬牛頓法求出f在該單元內的局部坐標(ri; Si,、)，并求出f與該單元中心的相對距離式=仏2矛；步驟3. 1.4:找出與f相對距離最小的單元et，即dt = min ((IljCl2,…，dm)，該單元即為流體節點f所對應的結構單元； 變形傳遞子模塊包含下列步驟步驟3. 2. 1 設流體節點f在其所對應結構單元e上的投影點為f’，f'在整體坐標系 下的坐標為(Xf,，yf,，zf,)，忽略流體節點與其投影點之間的距離，近似有xf, xf，yf, ^yf, Zf, ^ zf，投影點f’在結構單元e中的局部坐標(r，s, t)通過求解下列方程組得到 ‘ 8 步驟3. 2. 2 在求出投影點f’在結構單元e中的局部坐標(r，s, t)后，如果投影點在 面r= 1上，忽略流體節點與其對應投影點之間的距離，強制令r= 1 ；如果投影點在面S = 1上，忽略流體節點與其對應投影點之間的距離，強制令s = 1 ；如果投影點在面t = 1上， 忽略流體節點與其對應投影點之間的距離，強制令t = 1 ；步驟3. 2. 3 假設投影點f’的位移為uf，，ui為結構有限元節點i的位移，通過下列公 式求解流體節點f的位移％ /=1載荷傳遞子模塊包含下列步驟步驟3. 3. 1 設流體表面網格c的中心點為f。，計算流體壓力沿網格c的積分Ffc = I1PdS^pfc-Ac，其中Ω為網格c的表面，ρ為耦合界面上的流體壓力，&為中心點f。上的壓力，Α。為網格c的表面積；步驟3. 3. 2 設流體表面網格中心點f。對應的結構單元為e，通過計算/(i = 1， 2，…，8)可以得到結構單元e上8個節點的節點力，其中Ni為f。在結構單元e上投影點 處的形函數；步驟3. 3. 3 將結構單元e上8個節點的節點力集合入結構有限元節點載荷列陣； 模塊四顫振應力輸出模塊顫振應力輸出模塊通過葉片的瞬態位移響應計算出葉片的瞬態顫振應力，并輸出各種 供顫振分析使用的后處理文件；該模塊包括如下步驟步驟4. 1 通過葉片的瞬態位移求出葉片的瞬態顫振應力，即3δ =De (5)其中ε為應變列陣，σ為應力列陣，D為彈性矩陣；步驟4. 2 將葉片隨時間變化的位移、顫振應力和壓力輸出為可被Origin軟件讀取的 格式；步驟4. 3 將轉子隨時間變化的質量流量和壓比輸出為可被Origin軟件讀取的格式； 步驟4. 4:將葉片瞬間的位移云圖和顫振應力云圖輸出為可被Tecplot軟件讀取的格式；步驟4. 5 如果整個計算結束，輸出該級轉子葉片在壓氣機特性圖上的顫振邊界； 模塊五初值計算模塊初值計算模塊通過反復的非線性迭代獲得葉片的初始靜變形和穩態定常流場，作為以 后雙向迭代計算的初值，用以加快整個系統的收斂；該模塊包含如下步驟步驟5. 1 首先通過三維模型設計軟件UG建立葉片的三維實體模型；然后將實體模型 導入網格劃分軟件HyperMesh中，得到葉片的結構有限元網格；將實體模型導入網格劃分 軟件TurboGrid中，得到葉輪機內流場的流體計算網格；步驟5. 2 調用結構計算模塊得到葉片在離心力作用下的初始位移S^ATFc,其中Ktl是 結構的小位移剛度矩陣，fc為葉片在轉動中的離心力列陣；步驟5. 3 調用數據轉換模塊將葉片變形轉化為流體域的邊界運動； 步驟5. 4 通過流體計算模塊調用流體動力學計算軟件，進行網格更新，然后計算新的 定常流場；步驟5. 5 調用數據轉換模塊將葉片表面的定常壓力轉換為結構有限元模型上的節點力；步驟5. 6 調用結構計算模塊獲得葉片的非線性變形，并對結構的剛度矩陣進行更新； 計算葉片非線性變形的具體步驟如下步驟5. 6. 1 根據步驟5. 2求得的初始位移％，求得應力剛化矩陣Ks、大變形剛度矩陣 Kl和反力Fm;步驟 5. 6. 2 根據(KQ-KC+K。+IQda1 = Fc-Fnr 求出增量位移(Ia1 ； 步驟5. 6. 3 計算總位移a: = a0+dai ；步驟5. 6. 4:如果求得的an滿足收斂條件I |dan| |2 < ￡||an||2，其中ε取0.001，則 非線性迭代收斂，計算結束，指向步驟5. 7 ；否則用ai代替％，前往步驟5. 3 ；步驟5. 7 求得新的結構剛度矩陣K = K0-KC+KS+KL ；此時的葉片靜變形和定常流場就作 為以后雙向迭代計算的初始值； 模塊六雙向迭代模塊雙向迭代模塊反復交替調用結構求解器和流體求解器，同時采用數據轉換模塊在結構 求解器和流體求解器之間傳遞邊界數據，通過流體與結構雙向的交互作用，在時間上推進 整個由葉片及其周圍流場組成的流固耦合系統；雙向迭代模塊包含下列步驟 步驟6. 1 通過數據轉換系統將葉片的變形轉換為流體計算域的邊界運動； 步驟6. 2 通過流體計算模塊調用流體動力學計算軟件，進行網格更新，然后計算新的 非定常流場；步驟6. 3 通過數據轉換模塊將葉片表面的非定常壓力轉換為結構有限元模型上的節點力步驟6. 4 通過結構計算模塊求得新的葉片變形；步驟6.5 對雙向迭代的計算時間步進行調整；該流固耦合自適應時間步技術包括下 列步驟步驟 6. 5. 1 計算廠的，)2+(uru/ Y+…+(M1TU1/ )2]/(—/+.“+"力，其中 Ui 為當前時間步計算的節點i在某個自由度的位移，u/為上個時間步計算的節點i在某個自由 度的位移，η為有限元節點的數目；步驟6. 5. 2 取各個自由度的d的平均值為Cltl，用以表征當前時間步的收斂速度； 步驟6. 5. 3 用e代表收斂速度的閾值，用t代表當前的時間步長，如果Cltl < 0. 5e，則 新的時間步長tnew = 2t ；如果dQ > 2e，則新的時間步長tnew = 0. 5t ；如果0. 5e < d0 < 2e, 則時間步長不變，tnew = t ；步驟6. 5. 4 在獲得流固耦合計算的物理時間步長后，對結構動力方程瞬態求解的 時間步長和流體非定常求解的時間步長作相應調整；步驟6.6 通過計算所得的葉片瞬態位移響應，觀察顫振的發展進程；如果響應曲線衰 減，則不會發生顫振，計算結束；如果響應曲線增長，發生顫振，指向步驟6. 1 ；如果發生顫 振且計算達到預先指定時間點，計算結束，指向步驟6. 7輸出顫振應力；步驟6. 7 通過顫振應力輸出模塊，將結構的瞬態位移響應轉化為瞬態顫振應力，并輸 出為文件，供后處理軟件讀取顯示。
全文摘要
本發明公開了一種基于時域雙向迭代的葉輪機葉片顫振應力預測方法，其特征在于，將葉片及其周圍流場視為一個三維流固耦合系統，設計了一套時域的雙向迭代方法，通過交替求解葉片變形和非定常流場得到葉片的顫振應力。該方法在計算機中設定以下模塊結構計算模塊、流體計算模塊、數據轉換模塊、顫振應力輸出模塊、初值計算模塊和雙向迭代模塊。通過非線性迭代獲得葉片靜變形和穩態流場作為初值；交替調用結構計算模塊和流體計算模塊在時間上推進整個系統；通過數據轉換模塊傳遞流固邊界信息；輸出時間歷程上的顫振應力。本發明實現了葉片與流場的一體化計算，考慮了耦合系統的非線性，能觀察整個顫振發展進程，預測葉片的顫振應力。
文檔編號G06F17/50GK101908088SQ20101023721
公開日2010年12月8日 申請日期2010年7月22日 優先權日2010年7月22日
發明者徐可寧, 王延榮 申請人:北京航空航天大學