董雷霆 賀雙新

(北京航空航天大學 航空科學與工程學院, 北京 100083)

Abstract: Rotational components undergo complex loads which can easily initiate cracks, resulting in fatigue and fracture failures. For thermal-loaded rotational components with cracks, the finite element method(FEM) was applied to the global structure without cracks, thus exerting the advantage of FEM in numerical analysis of large-scale structure; meanwhile, the symmetric Galerkin boundary element method (SGBEM) was employed to the subdomain around the crack, thereby developing the advantage of SGBEM in fracture analysis. The SGBEM super element, which took the influence of the rotational inertia loading and the thermal loading into consideration, was developed in this paper. The stiffness matrix of the obtained SGBEM super element was symmetric and positive semidefinite, where the rotational inertia loading and the thermal loading only influence the equivalent force vector. Thus, the SGBEM super element can be directly assembled with the stiffness matrix and the equivalent force vector of the FEM, and be coupled with the FEM to solve the thermoelastic fracture problem of rotational components. Stress intensity factors of the cracked rotational disk undergoing the thermal loading were computed, which validates the effectiveness of SGBEM super element and FEM coupling method.

Keywords: symmetric Galerkin boundary element method (SGBEM); super element; thermoelasticity;rotational inertia load; stress intensity factor

渦輪盤和轉(zhuǎn)子葉片等旋轉(zhuǎn)部件的疲勞斷裂是航空發(fā)動機的常見多發(fā)故障[1-3],嚴重影響航空發(fā)動機和飛機的安全性與可靠性。 旋轉(zhuǎn)部件的高效仿真分析,特別是應力強度因子等斷裂參數(shù)的精確快速計算,對飛機和航空發(fā)動機的結(jié)構(gòu)仿真分析和結(jié)構(gòu)完整性評估具有重要意義[4-5]。

有限元法(finite element method,FEM)[6]廣泛應用于大型復雜結(jié)構(gòu)的仿真分析,已成為當今工程技術(shù)領(lǐng)域應用最為廣泛、成效最為顯著的數(shù)值分析方法。 FEM 采用連續(xù)函數(shù)作為形函數(shù),然而在裂紋幾何面處的位移場具有不連續(xù)性,并且在裂紋尖端處的應力場具有奇異性[7],因此FEM在處理裂紋這種強不連續(xù)問題方面具有難度,需要在裂紋尖端的高應力區(qū)劃分致密的網(wǎng)格[8],這增加了問題求解規(guī)模,且限制了斷裂力學參數(shù)的求解精度。

對稱伽遼金邊界單元法(symmetric Galerkin boundary element method,SGBEM)[9-10]在結(jié)構(gòu)斷裂仿真分析方面具有優(yōu)勢。 對于三維斷裂力學問題,由于SGBEM 無需內(nèi)部網(wǎng)格,僅在裂紋面劃分網(wǎng)格,相比有限元,大大降低了裂紋尖端應力集中區(qū)域網(wǎng)格劃分難度;其面力邊界積分方程中的檢驗函數(shù)選取為位移場[11],可以精確計算裂紋張開位移,配以裂尖1/4 結(jié)點單元后[12],SGBEM 可得到很高的應力強度因子計算精度[13];并且SGBEM 的位移和面力邊界積分方程是弱奇異的[14],這降低了數(shù)值積分計算難度。 但由于SGBEM 的系數(shù)矩陣是滿陣[15],隨著自由度數(shù)目的增加,會導致計算規(guī)?？焖僭黾?因此難以用于求解復雜結(jié)構(gòu)。

SGBEM 超單元FEM 耦合法[16]結(jié)合了FEM求解大型復雜結(jié)構(gòu)的優(yōu)勢和SGBEM 求解斷裂問題的優(yōu)勢。 該方法在裂紋附近局部子域使用SGBEM 求解,構(gòu)造了SGBEM 超單元;在全局無裂紋區(qū)域使用FEM 求解;由于SGBEM超單元對稱且半正定,可直接裝配至全局有限元的剛度矩陣中,從而耦合求解含裂紋的復雜結(jié)構(gòu)。 該方法對飛機二維結(jié)構(gòu)如加筋板裂紋黏結(jié)修復、耳片接頭裂紋疲勞失效等實際問題可進行高效的仿真分析。 但該方法沒有考慮旋轉(zhuǎn)慣性載荷和熱載荷對含裂紋局部子域產(chǎn)生的影響,不能直接用于三維旋轉(zhuǎn)部件的熱彈性斷裂分析。

本文將SGBEM 超單元FEM 耦合法推廣至三維旋轉(zhuǎn)部件的熱彈性斷裂問題。 推導了三維情況下考慮熱載荷和旋轉(zhuǎn)慣性載荷的SGBEM 超單元,將SGBEM 超單元應用于含裂紋局部子域,并給出SGBEM 超單元FEM 耦合的過程,計算算例驗證了該方法。

1 旋轉(zhuǎn)部件熱彈性SGBEM 超單元

1.1 位移和面力邊界積分方程

如圖1 所示,設(shè)彈性體Ω內(nèi)存在體力f和溫度場T,源點ξ和場點x均位于彈性體Ω的邊界?Ω上。 伽遼金式位移邊界積分方程[17-18]為

圖1 熱彈性體Fig.1 Thermoelastic body

伽遼金式位移邊界積分方程式(1)是將傳統(tǒng)邊界元法中的位移邊界積分方程寫為弱形式,并對核函數(shù)使用亥姆霍茲(Helmholtz)分解進行弱奇異化后得到的;伽遼金式面力邊界積分方程式(2)是將強奇異的應力邊界積分方程寫為弱形式,并對核函數(shù)使用亥姆霍茲分解進行弱奇異化后得到的。 將伽遼金式位移和面力邊界積分方程在邊界網(wǎng)格上離散,再求解相應的線性方程組即可得到邊界上的位移和面力,從而可以求解裂紋的張開位移,并使用位移外推法獲得裂紋應力強度因子。

1.2 SGBEM 超單元構(gòu)造格式

將伽遼金式位移邊界積分方程式(1)和伽遼金式面力邊界積分方程式(2)離散。 為簡化推導過程,本文不列出一個邊界單元系數(shù)矩陣組裝到對稱伽遼金邊界單元法中系數(shù)矩陣的過程。 不失一般性,假設(shè)邊界上位移和面力均為未知的。

邊界上位移可離散為

1.3 含裂紋SGBEM 超單元

添加此約束條件后,無需將2 個裂紋面分別劃分網(wǎng)格,僅需在一個裂紋面劃分網(wǎng)格即可,可直接計算裂紋面張開位移。

式(13)代入SGBEM 超單元的構(gòu)造格式(9),即得含裂紋彈性體的SGBEM 超單元構(gòu)造格式。

2 SGBEM 超單元FEM 耦合法

2.1 SGBEM 超單元FEM 耦合構(gòu)造格式

如圖1 所示,局部子域和全局有限元區(qū)域相交面上有重合的面單元和重合的結(jié)點。

位移有限元的表達格式可寫為

分別裝配式(14)、式(15)的剛度矩陣和等效結(jié)點力向量,可得SGBEM 超單元FEM 耦合構(gòu)造格式:

注意到,有限元區(qū)域與局部子域相交面上存在著相交結(jié)點,有限元在相交結(jié)點上的面力tI與SGBEM 超單元在相交結(jié)點上的面力ˉtI大小相等,方向相反,即tI=- ˉtI,故式(16)不顯示相交結(jié)點上的面力。

2.2 SGBEM 超單元FEM 耦合求解流程

SGBEM 超單元FEM 耦合求解流程采用自動化的軟件完成,用戶僅需要獨立劃分有限元網(wǎng)格與裂紋面網(wǎng)格。 如圖2 所示,詳細耦合求解過程如下。

圖2 SGBEM 超單元FEM 耦合求解流程Fig.2 Process of SGBEM super element and FEM coupling solution

1) 搜索局部子域。 當旋轉(zhuǎn)部件劃分好有限元網(wǎng)格后,根據(jù)旋轉(zhuǎn)部件中裂紋位置,找到裂紋附近有限元單元,這些有限元單元構(gòu)成了局部子域有限元網(wǎng)格。

2) 消除內(nèi)部結(jié)點。 注意到,裂紋面網(wǎng)格是另外插入的,裂紋面網(wǎng)格結(jié)點與局部子域有限元網(wǎng)格結(jié)點沒有任何關(guān)系。 對于內(nèi)部裂紋,需要將局部子域有限元網(wǎng)格中內(nèi)部結(jié)點消去,僅保留邊界面單元,這些邊界面單元和裂紋面單元共同構(gòu)成了局部子域邊界網(wǎng)格,如圖3 所示,為清楚顯示內(nèi)部,該圖隱藏了部分邊界單元。

3) 重構(gòu)局部子域。 對于面裂紋,以圖2 中④所示半圓形裂紋為例,其在局部子域中的位置如圖3所示。 裂紋前沿(即半圓曲線)上的結(jié)點在局部子域內(nèi)部,無需處理;但該裂紋的開口線(即半圓直徑線)上的結(jié)點既在裂紋面上,也在局部子域邊界面上。 因此,需要重新劃分局部子域邊界面(見圖2 中⑤部分)。 注意,局部子域和全局有限元相交面上的邊界面單元不可以重構(gòu),因為相交面上的結(jié)點是局部子域和全局有限元共有的。

圖3 局部子域邊界網(wǎng)格和裂紋網(wǎng)格Fig.3 Subdomain boundary mesh and crack mesh

4) 裝配剛度矩陣。 在獲得局部子域邊界元網(wǎng)格后,分別求解局部子域SGBEM 超單元和全局有限元,獲取SGBEM 超單元的剛度矩陣和有限元的剛度矩陣,再通過結(jié)點編號直接裝配。

5) 施加載荷。 根據(jù)式(9)等號左端第2 項求解局部子域旋轉(zhuǎn)慣性載荷和熱載荷產(chǎn)生的等效載荷列向量,并求解全局有限元中的載荷列向量。再按照式(16)等號左端第3 項將局部子域等效載荷列向量裝配至全局有限元中的等效載荷列向量。

6) 耦合求解。 步驟4)和步驟5)給出局部子域SGBEM 超單元和全局有限元的剛度矩陣和等效載荷列向量的整體裝配形式。 求解該線性方程組,可獲得局部子域SGBEM 超單元和全局有限元所有結(jié)點的位移解,包括裂紋張開位移,從而可以使用位移外推法求解裂紋應力強度因子。

3 含裂紋受熱旋轉(zhuǎn)圓盤算例

本節(jié)使用SGBEM 超單元FEM 耦合法求解圖2中的受熱旋轉(zhuǎn)圓盤,以展示其優(yōu)點。

如圖4 所示,圓盤內(nèi)徑Ri=0.1 m,外徑Ro=0.2 m,厚t=0.02 m,彈性模量E=1 000 MPa,泊松比ν=0.3,線膨脹系數(shù)α=2.17 ×10-5/℃。 繞圓盤中心軸的轉(zhuǎn)速為ω,轉(zhuǎn)速ω的取值如表1 所示。 沿徑向的溫度分布為

圖4 受熱旋轉(zhuǎn)圓盤Fig.4 Rotational disk undergoing thermal loading

式中:系數(shù)k的取值如表1。 在圓盤內(nèi)環(huán)面中間存在一半圓形面裂紋,半徑為r= 0. 004 m(見圖3);在圓盤內(nèi)環(huán)面中心處施加剛體位移約束。

該圓盤的SGBEM 超單元FEM 耦合法求解流程見2.2 節(jié)。 分別使用SGBEM 超單元FEM 耦合法和商業(yè)有限元軟件ANSYS Workbench 計算半圓形裂紋圓弧中點處的I 型裂紋應力強度因子。 ANSYS Workbench 所使用的有限元網(wǎng)格見(見圖5(a)),該網(wǎng)格比圖2 中⑥部分的耦合法所用的網(wǎng)格更加密集;且ANSYS Workbench 在裂紋附近使用了多種體單元,有五面體棱柱單元、六面體單元、五面體金字塔單元以及四面體單元(見圖5(b))。 可見,耦合方法使用了更加稀疏的網(wǎng)格,在裂紋尖端處僅使用了四邊形單元,網(wǎng)格也更加簡單(見圖2 中④部分)。

圖5 ANSYS Workbench 中圓盤和裂紋的有限元網(wǎng)格Fig.5 Disk and crack FEM mesh used by ANSYS Workbench

受熱旋轉(zhuǎn)圓盤I 型裂紋應力強度因子計算結(jié)果如表1 所示,不同轉(zhuǎn)速和溫度下,SGBEM 超單元FEM 耦合法與ANSYS Workbench 的計算結(jié)果相比,誤差均小于2%,可見SGBEM 超單元FEM耦合法可有效用于旋轉(zhuǎn)部件的熱彈性斷裂分析。

表1 不同轉(zhuǎn)速和溫度下圓盤裂紋應力強度因子Table 1 Stress intensity factors of crack in disk for different rotational speeds and temperature

該算例中本文方法在不含裂紋的全局區(qū)域進行稀疏的有限元剖分;在裂紋附近僅使用二維單元進行邊界元的計算,僅需要對裂紋面和包絡裂紋的局部子域表面進行簡單的二維網(wǎng)格劃分,避免了裂紋附近細致的三維有限元網(wǎng)格剖分,可以顯著減少單元數(shù)量。 關(guān)于SGBEM-FEM 法模擬斷裂與疲勞問題的在網(wǎng)格剖分和計算效率方面的優(yōu)勢,請見文獻[19-20]中給出的詳細論述。

4 結(jié) 論

本文給出了三維旋轉(zhuǎn)部件SGBEM 超單元的具體構(gòu)造格式。 本文將SGBEM 超單元與有限元耦合,計算了旋轉(zhuǎn)圓盤在熱載荷作用下半圓形裂紋的應力強度因子,并與商業(yè)有限元軟件ANSYS Workbench 結(jié)果進行對比,對比結(jié)果表明了該SGBEM 超單元FEM 耦合法可有效用于旋轉(zhuǎn)部件的熱彈性斷裂分析。 通過SGBEM 超單元FEM 耦合求解流程可知:

1) 有限元網(wǎng)格獨立劃分,不受裂紋干擾。 由于裂紋面不使用有限元單元,有限元網(wǎng)格可自由按需劃分,無需考慮裂紋尖端應力集中區(qū)有限元網(wǎng)格加密的困難。

2) 裂紋面網(wǎng)格獨立劃分,難度降低。 裂紋僅使用邊界單元,可自由按需劃分,且內(nèi)部裂紋無需重構(gòu)網(wǎng)格。 雖然面裂紋存在網(wǎng)格重構(gòu)問題,但僅需在局部子域邊界曲面上劃分二維網(wǎng)格,相比傳統(tǒng)有限元需在裂紋附近區(qū)域劃分致密的三維網(wǎng)格,難度大大降低。

3) 耦合過程簡單。 由于SGBEM 超單元僅存在位移自由度,不存在面力自由度,從SGBEM 超單元FEM 耦合構(gòu)造格式可知,其與有限元的裝配耦合過程與有限元中一個單元的剛度矩陣裝配至整體剛度矩陣的過程完全一致。

4) 裂紋尖端網(wǎng)格簡單。 SGBEM 超單元FEM,耦合法在裂紋前沿使用1/4 結(jié)點面單元。而有限元需要使用多種形狀的體單元劃分不規(guī)則形狀的裂紋前沿。

5) 有助于快速仿真裂紋擴展。 當進行裂紋擴展仿真時,每當裂紋擴展一段距離,對于內(nèi)部裂紋,耦合方法僅需在裂紋前沿增加一層邊界面單元即可繼續(xù)仿真;對于面裂紋,僅需再次重構(gòu)局部子域裂紋附近處的邊界面網(wǎng)格。 但對有限元而言,需要重新劃分裂紋附近區(qū)域體單元網(wǎng)格。

附錄A：

附錄B：

當彈性體作用有旋轉(zhuǎn)慣性載荷時,如果使用散度定理將旋轉(zhuǎn)慣性載荷引起的域積分轉(zhuǎn)換為邊界積分,則伽遼金式位移和面力邊界積分方程式(1)和式(2)中關(guān)于旋轉(zhuǎn)慣性載荷的函數(shù)為