蘇毅，王生楠，高文

（西北工業(yè)大學(xué)航空學(xué)院，陜西西安 710072）

廣義擴(kuò)展有限元及其在裂紋擴(kuò)展中的應(yīng)用

蘇毅，王生楠，高文

（西北工業(yè)大學(xué)航空學(xué)院，陜西西安 710072）

廣義擴(kuò)展有限元是廣義有限元和擴(kuò)展有限元兩者結(jié)合起來形成的一種新的數(shù)值方法。介紹了廣義擴(kuò)展有限元的基本原理并推導(dǎo)了相應(yīng)的公式，提出了將Westergaard裂紋尖端奇異場的基函數(shù)作為結(jié)點位移插值函數(shù)，探討了數(shù)值積分策略，給出了裂紋尖端應(yīng)力強(qiáng)度因子的計算方法，編寫廣義擴(kuò)展有限元程序。通過典型含裂紋平板的計算，表明廣義擴(kuò)展有限元計算應(yīng)力強(qiáng)度因子精度更高，也不需要劃分過密的網(wǎng)格。

廣義擴(kuò)展有限元；廣義有限元；擴(kuò)展有限元；數(shù)值積分策略；應(yīng)力強(qiáng)度因子

有限元法在解決工程問題中發(fā)揮著越來越重要的作用，隨著其應(yīng)用的深入，對基于變分原理的有限元法的發(fā)展始終沒有停止過。在有限元法中，由于單元的位移多項式插值函數(shù)在描述幾何、材料等不連續(xù)性上的不足，僅要求有限元網(wǎng)格能夠描述區(qū)域的幾何特征、材料變化等，只有足夠細(xì)的網(wǎng)格才能達(dá)到所需的精度，這在很大程度上限制了有限元法的應(yīng)用。為了解決有限元法在處理不連續(xù)性上的不足，人們發(fā)展了一些新的有限元法，其中基于單位分解法的廣義有限元法（GFEM）和擴(kuò)展有限元法（XFEM）最具代表性，它們都是在有限元法思想上的延伸。GFEM通過在結(jié)點處引入廣義自由度，對結(jié)點進(jìn)行再次插值，從而提高數(shù)值逼近精度，滿足對特定問題的特殊逼近要求［1-6］。廣義自由度的物理意義可以理解為以該結(jié)點為中心的支集上的局部逼近函數(shù)的權(quán)重。XFEM通過改進(jìn)單元的位移形狀函數(shù)使之包含問題不連續(xù)性的基本成分，從而放松對網(wǎng)格密度的過分要求。XFEM附加的自由度的物理意義為結(jié)點增強(qiáng)基函數(shù)的權(quán)重。XFEM在處理諸如模擬界面、裂紋及其擴(kuò)展、復(fù)雜流體等不連續(xù)問題時特別有效［7-10］。GFEM和XFEM在劃分用于插值逼近的有限元法網(wǎng)格時，都不需要考慮區(qū)域的內(nèi)部結(jié)構(gòu)細(xì)節(jié)。

文獻(xiàn)［11］結(jié)合擴(kuò)展有限元法和廣義有限元法的特點提出了廣義擴(kuò)展有限元法（GXFEM），并用于裂紋擴(kuò)展分析。在相同網(wǎng)格密度下，GXFEM較XFEM的計算精度高。但是該文的GXFEM繼承了XFEM的缺點，混合單元對解的收斂性有影響，而且對結(jié)點位移插值函數(shù)的形式?jīng)]有進(jìn)一步探討。

本文選用Westergaard裂紋尖端奇異場的基函數(shù)作為結(jié)點位移插值函數(shù)，對裂紋擴(kuò)展問題的剛度矩陣重新分析，采用Cholesky準(zhǔn)則計算剛度矩陣，減少了運算時間。編寫了GXFEM裂紋擴(kuò)展分析程序，通過對典型含裂紋板的數(shù)值分析，表明了本文建立的GXFEM在裂紋擴(kuò)展分析時的有效性和高的計算精度，且Westergaard裂紋尖端奇異場的基函數(shù)作為結(jié)點位移插值函數(shù)較常規(guī)的1階結(jié)點位移插值函數(shù)精度高。

1　GXFEM的基本原理

1.1不連續(xù)的位移場

廣義擴(kuò)展有限元法的基本的思想就是通過增加插值函數(shù)的階次提高計算精度并在單位分解的基礎(chǔ)上通過增加描述不連續(xù)性的附加函數(shù)，模擬裂紋的不連續(xù)性及裂紋尖端的奇異性，避免有限元法分析不連續(xù)問題時重新劃分網(wǎng)格的缺點。為此，構(gòu)造廣義擴(kuò)展有限元法的位移模式如下：

富集單元的結(jié)點位移并不是真實的結(jié)點位移，為了使所有的結(jié)點都是真實的結(jié)點位移，需要對（1）式進(jìn)行轉(zhuǎn)基處理，變換如下：

當(dāng)結(jié)點為常規(guī)單元的結(jié)點時，只有第1項；當(dāng)結(jié)點為被裂紋貫穿單元的結(jié)點時，有第1項和第2項；當(dāng)結(jié)點為裂紋尖端所在單元的結(jié)點時，有第1項和第3項；結(jié)點同時處于裂紋貫穿單元和裂紋尖端所在單元，應(yīng)優(yōu)先屬于裂紋尖端所在單元，加強(qiáng)方式選擇裂紋尖端所在單元結(jié)點加強(qiáng)方式，如圖1所示。

位移模式可表示為：

圖1　裂紋的加強(qiáng)結(jié)點

式中，常規(guī)單元的結(jié)點位移插值函數(shù)和裂紋貫穿單元的結(jié)點位移插值函數(shù)為0階：

對應(yīng)每個結(jié)點有2個自由度。裂紋尖端單元的結(jié)點位移插值函數(shù)為1階：

對應(yīng)每個結(jié)點有6個自由度。結(jié)點位移的插值函數(shù)階次的提高，廣義擴(kuò)展有限元形函數(shù)的階次也相應(yīng)的提高。對于裂紋尖端和裂紋貫穿單元結(jié)點來說，結(jié)點的廣義形函數(shù)為：

1.2Westergaard裂紋尖端奇異場

對于裂紋問題，裂紋尖端附近的漸進(jìn)位移場通?？梢赃\用Westergaard裂紋尖端奇異場：

將Westergaard裂紋尖端奇異場的基函數(shù)作為結(jié)點位移插值函數(shù)，得到F為：

1.3廣義擴(kuò)展有限元法離散方程的建立

將廣義擴(kuò)展有限元法的位移模式代入彈性體虛功方程，經(jīng)過變分運算，可得到離散的線性方程組為：

式中：d是結(jié)點廣義位移，K和f分別為結(jié)構(gòu)總剛度矩陣和外力矢量，K和f按常規(guī)步驟由單元剛度矩陣組集而成。單元層次上的K和f分別為：

單元荷載列陣f：

單元荷載力矢量分量表示如下：

常規(guī)應(yīng)變矩陣：

裂紋貫穿單元附加應(yīng)變矩陣：

裂紋尖端單元附加應(yīng)變矩陣：

采用J積分計算應(yīng)力強(qiáng)度因子，傳統(tǒng)的J積分形式［12］為：

式中：Γ是積分路徑，nj是該路徑上弧元素外法線的方向余弦。J積分為裂紋擴(kuò)展能量釋放的能力，與積分路徑無關(guān)。在復(fù)合加載模式下，J積分與相應(yīng)的應(yīng)力強(qiáng)度因子有如下形式［13］：

本文應(yīng)用最大周向拉應(yīng)力準(zhǔn)則，確定裂紋尖端擴(kuò)展方向θc：

同時采用裂紋增量形式模擬裂紋的擴(kuò)展過程。

2　重新分析GXFEM中裂紋的擴(kuò)展問題

（14）式也可以寫成

雙出口教室的人群疏散建模和模擬···················董力耘 藍(lán)冬愷 段曉茵 (2,217)

從（26）、（27）式可以得到：

隨著裂紋擴(kuò)展，（27）式中子矩陣A12、A22發(fā)生改變，在總剛度矩陣中，與裂紋尖端聯(lián)系的部分矩陣為：

A11在裂紋擴(kuò)展的迭代中不發(fā)生改變可以重復(fù)使用。對于新的貫穿單元結(jié)點，以上算法可以用來計算新的剛度矩陣分量的分解因式，A11的迭代結(jié)果可以附加到新的因式分解上，這樣L21可以得到。裂紋尖端的剛度矩陣可以用來求L22，用Cholesky準(zhǔn)則可以計算出整體剛度矩陣，因為需要計算的部分比較小，節(jié)省了計算時間。

3　幾點說明

3.1數(shù)值積分策略

在有裂紋這種強(qiáng)不連續(xù)存在時，擴(kuò)展有限元單元不需要重新劃分網(wǎng)格，廣義擴(kuò)展有限元也有這樣的特點。廣義擴(kuò)展有限元進(jìn)行單元分解的目的是數(shù)值積分。在裂紋經(jīng)過單元采用和常規(guī)有限元相同數(shù)目的積分顯然是不能達(dá)到要求的。為了保證積分的精度，采用以下的積分方案：

圖2　單元的分塊積分策略

如圖2所示，對于傳統(tǒng)四結(jié)點等參單元的單元勁度矩陣即沒有結(jié)點加強(qiáng)的單元，其積分方案采用2×2個高斯積分點。由于階次的提高（較傳統(tǒng)四結(jié)點等參單元被積函數(shù)階次至少要高一次），對廣義擴(kuò)展有限元四結(jié)點等參單元計算單元勁度矩陣宜采用9個高斯積分點。

沒有裂紋穿過，但有結(jié)點加強(qiáng)的單元，采用6×6個高斯積分點。

若4個結(jié)點都是Heaviside加強(qiáng)，在裂紋貫穿單元，裂紋將單元分為2個子域，分別由每個子域的角點形成Delaunay三角形，每個三角形內(nèi)采用3×3個高斯積分點。

在裂紋尖端單元，縫尖延長，將單元分成兩個子域，分別由每個子域的角點形成Delaunay三角形，每個三角形內(nèi)采用9×9個高斯點積分。在積分時將積分區(qū)分為這些三角形，并不增加整體的自由度數(shù)，只是積分的需要。

3.2廣義擴(kuò)展有限元與有限元的聯(lián)合運用

為了解決實際問題的計算效率和精度這一問題，可針對具體實際問題的結(jié)構(gòu)形式或所關(guān)心問題的區(qū)域，將廣義擴(kuò)展有限元與擴(kuò)展有限元聯(lián)合運用。因為廣義擴(kuò)展有限元兼具廣義有限元和擴(kuò)展有限元的特點，因此廣義擴(kuò)展有限元中結(jié)點位移的插值多項式階次可以按照需要任意選取，而不影響位移的協(xié)調(diào)性問題，對廣義擴(kuò)展有限元的常規(guī)四結(jié)點單元中按照1階的結(jié)點位移插值函數(shù)建立的程序，每個結(jié)點有6個廣義位移（自由度），只要將某些結(jié)點后4個自由度施加約束，那么就自然退化成傳統(tǒng)有限元的結(jié)點位移，因此只要對約束信息做恰當(dāng)?shù)奶幚恚瑹o需修改程序就可以實現(xiàn)廣義擴(kuò)展有限元與擴(kuò)展有限元的聯(lián)合運用。因為廣義有限元中結(jié)點的自由度為廣義位移，故在求得結(jié)點廣義位移后，還需要通過結(jié)點的位移插值函數(shù)求的結(jié)點最終位移值。

3.3廣義擴(kuò)展有限元線性相關(guān)性的處理

由于廣義擴(kuò)展有限元的局部逼近函數(shù)大多存在線性相關(guān)性，廣義擴(kuò)展有限元方法的剛度矩陣如廣義有限元的剛度矩陣一般具有奇異性。廣義擴(kuò)展有限元方法中，零特征根的個數(shù)一般不知道。因此需要求解以半正定矩陣A為系數(shù)的線性方程組即Ac =b，文獻(xiàn)［15］提出了解決方法。

4　數(shù)值算例

4.1單邊裂紋平板受單向均勻拉伸

長度為a的單邊裂紋平板如圖3所示。

圖3　單邊裂紋平板受拉伸

L=6 m，W=3 m，受單向均布拉應(yīng)力σ=1 MPa，板的彈性模量E=10 MPa，泊松比為0.3。文獻(xiàn)［16］給出了該問題的應(yīng)力強(qiáng)度因子（SIF）的經(jīng)驗公式。有限元網(wǎng)格密度為15×31，采用不同方法計算得到的應(yīng)力強(qiáng)度因子列于表1。

從表1的計算結(jié)果可以看出，GXFEM較XFEM的計算結(jié)果更接近解析解，且使用（12）式的結(jié)點位移插值函數(shù)的計算結(jié)果較（9）式作為結(jié)點位移插值函數(shù)更接近解析解。

表1　單邊裂紋拉伸平板SIF計算結(jié)果比較

4.2中心斜裂紋平板受單向均勻拉伸

如圖4所示，一矩形板內(nèi)含有一條中心斜裂紋，矩形板的上端和下端受到豎直方向的單向拉伸載荷σ=1 MPa，矩形板的高和寬為W=10 m，裂紋長度為2a=1 m。材料彈性模量為E=71.7 GPa，泊松比為μ=0.33。該模型的應(yīng)力強(qiáng)度因子解析解取自文獻(xiàn)［16］。

圖4　中心斜裂紋板

選取網(wǎng)格密度為100×101。采用不同方法計算得到的應(yīng)力強(qiáng)度因子KⅠ和KⅡ隨裂紋角的變化曲線如圖5所示。

圖5　應(yīng)力強(qiáng)度因子隨裂紋角的變化曲線

從圖5可以看出，GXFEM比XFEM計算結(jié)果更接近解析解。但是因為網(wǎng)格密度較密，各種方法的精度差距不大，GXFEM的優(yōu)勢主要體現(xiàn)在網(wǎng)格密度不大的時候，仍能精確的得到應(yīng)力強(qiáng)度因子，而且本例是單向拉伸，高階位移插值函數(shù)的效果有限。

選取半裂紋長度0.5m，角度為30°，應(yīng)用廣義擴(kuò)展有限元模擬裂紋擴(kuò)展。分別采用4種網(wǎng)格密度，網(wǎng)格1為50×51，網(wǎng)格2為80×81，網(wǎng)格3為100 ×101，網(wǎng)格4為120×121，圖6給出了不同網(wǎng)格對應(yīng)的裂紋擴(kuò)展路徑。

圖6　不同網(wǎng)格對應(yīng)的裂紋擴(kuò)展路徑

從圖6可以看出，隨著網(wǎng)格密度的增加，裂紋擴(kuò)展的路徑幾乎重合，表明當(dāng)網(wǎng)格密度增加到一定程度時，可以很好的模擬裂紋擴(kuò)展。

5　結(jié) 論

從廣義擴(kuò)展有限元的位移場構(gòu)造和驗證算例兩方面可以看出，廣義擴(kuò)展有限元法的精確明顯高于擴(kuò)展有限元法，而且選用Westergaard裂尖奇異場的基函數(shù)作為結(jié)點位移插值函數(shù)（即（12）式）比采用一階結(jié)點位移插值函數(shù)（即（9）式）的精度要高。廣義擴(kuò)展有限元結(jié)點位移插值多項式可以任意選取，且不影響位移的協(xié)調(diào)性問題，其繼承了廣義有限元的優(yōu)點，但是對于其剛度奇異性問題還需要做更進(jìn)一步的研究。

廣義擴(kuò)展有限元可以退化為擴(kuò)展有限元、廣義有限元甚至常規(guī)有限元。在考慮不連續(xù)問題時，只需要對不連續(xù)的區(qū)域進(jìn)行結(jié)點自由度的廣義化和增加富集函數(shù)，對于其他的區(qū)域可以采用常規(guī)有限元。

［1］ Babushka I，Osborn.Generalized Finite Element Methods：Their Performance and Their Relation to Mixed Methods［J］.SIAM Journal of Numerical Analysis，1983，20（3）：510-535

［2］ Duarte C A，Babuska I，Oden J T.Generalized Finite Element Methods for Three-Dimensional Structural Mechanics Problems ［J］.Computer&Structures，2000，77：215-232

［3］ Strouboulis T，Copps K，Babuska I.The Generalized Finite Element Method［J］.Computer Methods in Applied Mechanics and Engineering，2001，190：4081-4193

［4］ Strouboulis T，Zhang L，Babuska I.Generalized Finite Element Method Using Mesh-Based Handbooks：Application to Problems in Domains with Many Voids［J］.Computer Methods in Applied Mechanics and Engineering，2003，192（28/29/30）：3109-1361

［5］ Duarte C A，Babuska I.Mesh-Independent Porthotropic Enrichment Using the Genernalized Finite Element Method［J］.International Journal for Numerical Methods in Engineering，2001，55（12）：1477-1492

［6］ Babuska I，Banerjee U，Osborn J E.Survey of Meshless and Generalized Finite Element Methods：a Unified Approach［J］.Acta Numerica，2003，12：1-125

［7］ Sukumar N，Chopp D L，Moran B.Extended Finite Element Method and Fast Marching Method for Three-Dimensional Fatigue Crack Propagation［J］.Engineering Fracture Mechanics，2003，70：29-48

［8］ 陳金龍，戰(zhàn)楠，張曉川.基于擴(kuò)展有限元法的裂尖場精度研究［J］.計算力學(xué)學(xué)報，2014，31（4）：425-430

Chen Jinlong，Zhan Nan，Zhang Xiaochuan.Numerical Study on the Accuracy of Crack Tip Field by Extended Finite Element Method［J］.Chinese Journal of Computational Mechanics，2014，31（4）：425-430（in Chinese）

［9］ 宋娜，周儲偉.擴(kuò)展有限元裂尖場精度研究［J］.計算力學(xué)學(xué)報，2009，26（4）：544-547

Song Na，Zhou Chuwei.Accuracy Study of Crack Tip Field in Extended Finite Element Method［J］.Chinese Journal of Computational Mechanics，2009，26（4）：544-547（in Chinese）

［10］Tarancon J E，Vercher A，Giner E.Enhanced Blending Elements for XFEM Applied to Linear Elastic Fracture Mechanics［J］. International Journal for Numerical Methods in Engineering，2009，77（1）：126-148

［11］Zhang Qing，Liu Kuan，Xia Xiaozhou，Yang Jing.Generalized Extended Finite Element Method and Its Application in Crack Growth Analysis［J］.Chinese Journal of Computational Mechanics，2012，29（3）：427-432

［12］Moran B，Shih C F.Crack Tip and Associated Domain Integrals from Momentum and Energy Balance［J］.Engineering Fracture Mechanics，1987，27（6）：615-642

［13］Yau F J，Wang S S，Corten H T.A Mixed-Mode Crack Analysis of Isotropic Solids Using Conservation Laws of Elasticity［J］. Journal of Applied Mechanics，1980，47：335-341

［14］Pais M J，Kim N H，Davis T.Reanalysis of the Extended Finite Element Method for Crack Initiation and Propagation［R］. AIAA-2010-2536

［15］李錄賢，劉書靜，張慧華，等.廣義有限元方法研究進(jìn)展［J］.應(yīng)用力學(xué)學(xué)報，2009，26（1）：96-108

Li Luxian，Liu Shujing，Zhang Huihua，et al.Researching Progress of Generalized Finite Element Method［J］.Chinese Journal of Applied Mechanics，2009，26（1）：96-108（in Chinese）

［16］中國航空研究院.應(yīng)力強(qiáng)度因子手冊［M］（增訂版）.北京：科學(xué)出版社，1993，12：251-252；348-349

China Aviation Institute.Manual of Stress Intensity Factor［M］（Revised Edition）.Beijing：Science Press，1993，12：251-252；348-349（in Chinese）

Generalized Extended Finite Element Method and Its Application in Crack Growth

Su Yi，Wang Shengnan，Gao Wen
（College of Aeronautics，Northwestern Polytechnical University，Xi′an 710072，China）

The generalized extended finite element method is a new numerical method that combines both the generalized finite element and the extended finite element.This paper presents its basic principles and derives its formulas.It proposes that the base function of the Westergaard crack tip field should be used as the node displacement interpolation function and then discusses the numerical integration strategy.It uses the generalized extended finite element method to calculate the stress intensity factor of a crack tip and then develops the crack propagation analysis programming.Finally it gives a numerical example of the propagation of a structure with an edge crack，and the numerical results show that the stress intensity factor calculated with the generalized extended finite element has a higher precision and that there is no need to divide too dense meshes.

crack tips；calculations；crack propagation；degrees of freedom（mechanics）；efficiency；errors；finite element method；integration；matrix algebra；stress intensity factors；stiffness matrix；vectors；generalized extended finite element method；generalized finite element；extended finite element；numerical integration strategy

O346.1

A

1000-2758（2015）06-0921-07

2015-04-23

蘇毅（1983—），女，西北工業(yè)大學(xué)博士研究生，主要從事飛機(jī)結(jié)構(gòu)疲勞斷裂可靠性及損傷容限的研究。