Hamilton體系下的擴展有限元方法

卿光輝，沈 雄

（中國民航大學航空工程學院，天津 300300）

Hamilton體系下的擴展有限元方法

卿光輝，沈 雄

（中國民航大學航空工程學院，天津 300300）

基于Hamilton體系下的彈性力學理論,推導了擴展有限元法的有關(guān)公式.具體計算了含裂紋的等厚度薄板的應(yīng)力強度因子,數(shù)值分析所得結(jié)果與精確解進行了比較,證明了本文方法的正確性.本文為擴展有限元研究提供了一種新的數(shù)值計算方法,拓展了擴展有限元的基礎(chǔ)理論.

擴展有限元；裂紋；Hamilton體系；應(yīng)力強度因子

19世紀80年代，Hamilton計算方法的陸續(xù)出現(xiàn)以及Hamilton體系理論研究的發(fā)展，使其在數(shù)值求解方面上的應(yīng)用得到前所未有的發(fā)展。20世紀末，Hamilton體系又被引入彈性力學中，給彈性力學問題的求解帶來了一種新的思路[1]。

擴展有限元法是在常規(guī)有限元框架內(nèi)求解不連續(xù)問題的數(shù)值計算方法，它是基于單位分解思想[2-3]，以常規(guī)有限元為基本框架，在位移插值函數(shù)中加入了模擬裂紋界面的跳躍函數(shù)和裂尖漸近位移場函數(shù)，求解過程與常規(guī)有限元方法基本相似，離散網(wǎng)格與裂紋的幾何位置無關(guān)，且是相互獨立的，在裂尖處不必加密網(wǎng)格，因此便于分析含有裂紋的不連續(xù)問題[4-6]。

本文的目是為擴展有限元研究提供一種新的數(shù)值計算方法，豐富擴展有限元的基礎(chǔ)理論，以及Hamilton理論體系的工程應(yīng)用范疇。

本文針對裂紋問題，在前人擴展有限元理論的基礎(chǔ)上，引入Hamilton理論，即采用一種強有力的半解析法[7]（在一個方向采用常規(guī)有限元插值，而在另一個方向采用狀態(tài)空間法給出控制方向的解析解答）。首先是將擴展的形函數(shù)代入Hamilton彈性力學理論中，就基于Hamilton彈性力學理論的相關(guān)公式進行了詳細的推導。其次，利用Matlab、Mathematica軟件編寫了相關(guān)程序，計算了含裂紋的等厚度薄板的裂紋尖端的應(yīng)力和位移，從而計算得到應(yīng)力強度因子，數(shù)值結(jié)果與理論值進行了比較，驗證了此方法是可行的。

1 擴展有限元

在常規(guī)有限元的基礎(chǔ)上，擴展有限元位移插值函數(shù)中，加入反映裂紋面的跳躍函數(shù)和裂尖漸近位移場函數(shù)[8]，根據(jù)前人的相關(guān)文獻資料可知，用Hamilton體系表示的彈性力學基本方程中，應(yīng)力和位移是完全獨立的，因此，對于應(yīng)力也可用相同的形函數(shù)進行插值處理，成立

其中：i是離散單元的節(jié)點集合；j是有裂紋穿過單元的節(jié)點集合；k是含裂紋尖端單元的節(jié)點集合，節(jié)點標記如圖1所示；Ni、Nj、Nk分別是對應(yīng)節(jié)點形函數(shù)；位移插值式（1a）中，ui為單元相應(yīng)節(jié)點的位移分別為單元相應(yīng)附加節(jié)點的位移；應(yīng)力插值式（1b）中，σi為單元相應(yīng)節(jié)點的應(yīng)力分別為單元相應(yīng)附加節(jié)點的應(yīng)力；H（x）為跳躍函數(shù)，反映裂紋面位移的不連續(xù)性

圖1 附加函數(shù)的加強節(jié)點Fig.1 Enrich nodes of addition function

如圖2所示，x是高斯點，x*是裂紋面上離x最近一點，n是x*處垂直于裂紋的法向量，F(xiàn)α（x，y）是裂尖漸近位移場的基本函數(shù)，若為四階漸近，在以裂尖為圓心的極坐標下取

圖2 光滑裂紋的坐標示意圖Fig.2 Smooth crack coordinates

2 Hamilton理論下的公式推導

對一個任意劃分的等厚度板，如圖3所示，用Hamiltonian等參元[7，9]離散x-y平面，場函數(shù)和形函數(shù)分別為

圖3 線性單元Fig.3 Linear element

其中

因此，Hamiltonian單元的控制微分方程[7，9-10]為

其中：Ci（i=1，2，…，6）是與材料彈性常數(shù)Cij（i，j=1，2，…，6）相關(guān)的組合項[7，10]；α=?/?x，β=?/?y，γ=?/?z為微分運算符為板側(cè)面所受的應(yīng)力邊界值。

各單元的形函數(shù)如下：

1）常規(guī)單元控制微分方程系數(shù)矩陣按上式進行計算，該單元形函數(shù)如式（5）所示

簡寫為

2）含加強節(jié)點單元的控制微分方程系數(shù)矩陣計算須擴充形函數(shù)，加強節(jié)點的形函數(shù)為

i）有跳躍函數(shù)加強的節(jié)點i

ⅱ）有裂尖漸近位移場函數(shù)加強的節(jié)點i

其中：Fαj（j=1，2，3，4）為式（3）Fα（x，y）的第j項，該單元形函數(shù)為

通過上面所求單元形函數(shù)，代入式（6）所有離散單元按常規(guī)有限元方法進行組裝，可求得單層板的整體控制微分方程為

其中：下標m表示第m層，上式的精確解為

其中

計算Tm和Sm的方法很多，考慮到計算結(jié)果的精度，采用精細積分方法計算Tm（z）[11]，采用五點拋物線規(guī)則來計算Sm，即

3 單元高斯積分方法

由于離散x-y平面的區(qū)域中，如圖1所示，如果部分單元的積分計算直接進行高斯積分，則積分結(jié)果將不夠精確；對于有裂紋存在的區(qū)域，需進行分塊積分，即將這些單元再分成幾個子單元，以裂紋為邊界劃分成幾個子單元，再進行高斯積分，這樣就能提高計算的精度。

有裂紋穿過的單元：分成4個子單元，如圖4所示；含裂紋尖端的單元：分成6個子單元，如圖5所示。本文具體單元相關(guān)的系數(shù)矩陣的積分計算方法，采用如下方式的高斯積分[12-13]：

圖4 完全貫穿單元Fig.4 Crack split element

圖5 裂尖單元Fig.5 Crack tip element

1）常規(guī)單元（無加強節(jié)點），積分點采用2×2的高斯點；

2）過渡單元（無裂紋貫穿，有加強節(jié)點），積分點采用4×4的高斯點；

3）貫穿單元（裂紋完全穿過），用Delaunay三角剖分，每個子三角形采用3個高斯點，即12個高斯點；

4）裂尖單元用Delaunay三角剖分，將裂紋延長，劃分成6個子三角形，每個子三角形采用13個高斯點，即78個高斯點。

4 應(yīng)力強度因子的計算

應(yīng)力強度因子是判斷裂紋開裂的重要依據(jù)，在擴展有限單元法中常采用相互作用積分來求解應(yīng)力強度因子。

相互作用積分是一個能量積分，積分路徑如圖6所示，裂紋尖端在回路內(nèi)，其定義式為[14-15]

其中：E*為材料常數(shù)E（彈性模量）和μ（泊松比）的組合項：E*=E（平面應(yīng)力），E*=E/（1-μ2）（平面應(yīng)變），從式（10）可知，若附加場滿足則用這種相互作用積分求解真實應(yīng)力-變形場的Ⅰ型的應(yīng)力強度因子同樣道理，若求解真實應(yīng)力-變形場的Ⅱ型應(yīng)力強度因子只要選擇滿足的附加場即可。

在實際應(yīng)用中為便于計算，把式（9）曲線積分變成曲面積分，如圖6所示，在回路Γ的外圍附近，設(shè)有一回路Γ0，兩回路之間的區(qū)域為A，則式（9）變?yōu)?/p>

圖6 相互作用積分示意圖Fig.6 Interaction integration

其中：C=Γ+C-+Γ0+C+；mj為圍線C上的外法線的單位向量，由格林（Green）公式將式（11）曲線積分轉(zhuǎn)化為曲面積分為

式（13）中定義q為一平滑的權(quán)函數(shù)。對于積分區(qū)域A，采用文獻[12]中的方法，以裂紋尖端為圓心、R為半徑作圓，其穿過的所有單元（如圖7所示，陰影部分）即為區(qū)域A，其中為積分區(qū)域因子，h為含裂紋尖端單元的面積，節(jié)點在圓內(nèi)；對應(yīng)的權(quán)函數(shù)q=1；節(jié)點在圓外，對應(yīng)的權(quán)函數(shù)q=0，單元內(nèi)某處的權(quán)函數(shù)值由單元節(jié)點權(quán)函數(shù)值進行插值計算

圖7 積分區(qū)域Fig.7 Integration domain

5 算例分析

如圖8所示，一等厚度薄板模型，尺寸寬W=1.0m，高L=2.0 m，厚度h=0.05 m，平板上邊界受均布拉應(yīng)力σ=1 MPa，下邊界在板的中間位置靠左側(cè)有一水平靜止裂紋，其長度a=0.5 m，平板彈性模量E=210 GPa，泊松比μ=0.3，這種含裂紋的平板模型的應(yīng)力強度因子精確解可通過如下表達式求得[16]

由式（14）和式（15）計算得：Kexact=2.877 2 MPa·采用本文的XFEM方法，將表面離散區(qū)域（如圖9所示）網(wǎng)格數(shù)分別劃分為45、102、119、133、147、161，厚度分兩層，取厚度為處的應(yīng)力和位移，按平面應(yīng)變問題來近似計算平板應(yīng)力強度因子，計算不同積分區(qū)域因子（rk=1.0～4.0）時的應(yīng)力強度因子數(shù)值解Knum。

從圖10及表1分析可知：

1）當r=1.0時，計算結(jié)果波動大且誤差較大（誤差>10%）；當r=2.0或2.5時，計算結(jié)果隨網(wǎng)格加密而趨于穩(wěn)定，但誤差大；當r=2.8、3.0和3.2時，計算結(jié)果隨網(wǎng)格加密而趨于穩(wěn)定，且誤差小（誤差<10%）。但當r=3.2且網(wǎng)格為45時，精度不夠理想。

2）引入Hamilton理論體系后，離散網(wǎng)格少的情況下就能得到較精確的數(shù)值解，誤差小于1%，而文獻[17]中的位移法程序計算的結(jié)果，網(wǎng)格非常密的情況下才能得到誤差小于2%的數(shù)值解。

圖8 幾何模型圖Fig.8 Geometric model

圖9 x-y表面離散網(wǎng)格Fig.9 Discrete mesh of x-y surface

圖10 應(yīng)力強度因子KΙFig.10 Stress intensity factor KΙ

表1 不同網(wǎng)格及積分半徑r下應(yīng)力強度因子K（MPa·m1/2）的數(shù)值解及誤差Tab.1 SIF numerical solution and error of K（MPa·m1/2）under different meshes and integration radius

6 結(jié)語

本文將Hamilton理論引入到擴展有限元來求解斷裂問題，采用了擴展有限元的位移模式，引入Hamilton理論后，用相同方式處理應(yīng)力，通過算例分析可以得到：

1）研究了積分半徑和網(wǎng)格密度對應(yīng)力強度因子計算結(jié)果的影響。引入Hamilton理論后在網(wǎng)格較稀疏的情況下就能得到較精確的計算結(jié)果，大大節(jié)省了有限元的前處理過程；積分半徑r在文獻[17]雖然已經(jīng)研究過，但同樣對本文方法的計算結(jié)果也有較大的影響。通過數(shù)值分析比較，本文在網(wǎng)格劃分為45～161的范圍內(nèi)，r=3.0時，數(shù)值結(jié)果都比較好。

2）在劃分合適的網(wǎng)格及積分半徑時，采用本文方法計算的結(jié)果誤差接近1%，說明本文方法是完全可行的。

3）本文的工作拓展了擴展有限元法的研究思路，并且充實了擴展有限元法的基礎(chǔ)理論。

[1]鐘萬勰.彈性力學求解新體系[M].大連：大連理工大學出版社，1995：62-176.

[2]MELENKJM，BABUSKAI.Thepartitionofunityfiniteelement method：basic theory and application[J].Computer Methods in Applied Mechanics and Engineering，1996，139：289-314.

[3]CHESSA J，WANG H，BELYTSCHKO T.On the construction of blending elements for local partition of unity enrichment finite elements[J]. International Journal for Numerical Methods in Engineering，2003，57（7）：1015-1038.

[4]李錄賢，王鐵軍.擴展有限元法（XFEM）及其應(yīng)用[J].力學進展，2005，35（1）：1-20.

[5]余天堂.擴展有限元法的數(shù)值方面[J].巖土力學，2007，28（增刊）：305-310.

[6]郭歷倫，陳忠富，羅景潤，等.擴展有限元法及應(yīng)用綜述[J].力學季刊，2011，32（4）：612-625.

[7]唐立民，褚致中，鄒貴平，等.混合狀態(tài)Hamilton元的半解析解和疊層板的計算[J].計算結(jié)構(gòu)力學及應(yīng)用，1992，9（4）：347-360.

[8]BELYTSCHKO T，BLACK T.Elastic crack growth in finite elements with minimal remeshing[J].International Journal for Numerical Methods in Engineering，1999，45：601–620.

[9]唐立民.彈性力學的混合方程和Hamilton正則方程[J].計算結(jié)構(gòu)力學及應(yīng)用，1991，8（4）：343-350.

[10]卿光輝，邱家俊，塔 娜.彈性體的Hamilton正則方程和整體加筋半的固有頻率分析[J].力學學報，2004，36（6）：749-756.

[11]MOLER C，VAN L C.Nineteen dubious ways to compute the exponential of a matrix.SIAM Review[J].Society for Industrial and Applied Mathematics，1978，20（4）：801-836.

[12]茹忠亮，朱傳銳，張友良，等.斷裂問題的擴展有限元法研究[J].巖土力學，2007，32（7）：2171-2176.

[13]喻葭臨，張其光，于玉貞，等.擴展有限元中非連續(xù)區(qū)域的一種積分方案[J].清華大學學報（自然科學版），2009，49（3）：352-355.

[14]MO?S N，DOLBOW J，BELYTSEHKO T.A finite element method for crack growth without remeshing[J].International Journal for Numerical Methods in Engineering，1999，46（1）：131-150.

[15]莊 茁，柳占立，成斌斌，等.擴展有限單元法[M].北京：清華大學出版社，2012：20-21.

[16]中國航空研究院.應(yīng)力強度因子手冊[M].北京：科學出版社，1981.

[17]王 敏.基于擴展有限元法的平板裂紋模型裂紋擴展研究[D].大連：大連理工大學，2012.

（責任編輯：楊媛媛）

Extended finite element method based on Hamilton theory

QING Guang-hui，SHEN Xiong
（College of Aeronautical Engineering，CAUC，Tianjin 300300，China）

Based on the elasticity theory of Hamilton system，the corresponding formulae of extended finite element method are derived，and the stress intensity factor of the equal thickness thin plate with crack is calculated.The numerical result is compared with the exact solution，verifying the effectiveness of this method.A new numerical approach to study the extended finite element is offered，and the basic theory of the extended finite element is expanded.

extended finite element method；crack；Hamilton system；stress intensity factor

O326

:A

:1674-5590(2015)01-0059-06

2013-10-18；

:2013-12-09

:中國民航大學科研基金項目(2012kye07)

卿光輝(1968—)，男，湖南新化人，教授，博士，研究方向為復合材料力學.