旋轉(zhuǎn)部件SGBEM-FEM 耦合熱彈性斷裂分析

董雷霆 賀雙新

(北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院, 北京 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ā)動(dòng)機(jī)的常見(jiàn)多發(fā)故障[1-3],嚴(yán)重影響航空發(fā)動(dòng)機(jī)和飛機(jī)的安全性與可靠性。 旋轉(zhuǎn)部件的高效仿真分析,特別是應(yīng)力強(qiáng)度因子等斷裂參數(shù)的精確快速計(jì)算,對(duì)飛機(jī)和航空發(fā)動(dòng)機(jī)的結(jié)構(gòu)仿真分析和結(jié)構(gòu)完整性評(píng)估具有重要意義[4-5]。

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

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

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

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

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

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

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

圖1 熱彈性體Fig.1 Thermoelastic body

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

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

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

邊界上位移可離散為

1.3 含裂紋SGBEM 超單元

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

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

2 SGBEM 超單元FEM 耦合法

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

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

位移有限元的表達(dá)格式可寫(xiě)為

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

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

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

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

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

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

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

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

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

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

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

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

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

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

如圖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(見(jiàn)圖3);在圓盤內(nèi)環(huán)面中心處施加剛體位移約束。

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

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

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

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

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

4 結(jié) 論

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

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

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

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

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

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

附錄A：

附錄B：

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