光滑擴(kuò)展有限元方法及其特點(diǎn)研究

徐志民， 謝 偉， 竇鵬鵬

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

1 引 言

隨著有限元法在工程中的廣泛應(yīng)用，其局限性漸漸凸顯。首先，常規(guī)有限元法的形狀函數(shù)是連續(xù)的，導(dǎo)致單元內(nèi)部材料參數(shù)不能階躍，因此分析斷裂力學(xué)的能力很弱。其次，由于常規(guī)有限元法在數(shù)值積分過(guò)程中需要進(jìn)行坐標(biāo)映射，對(duì)單元形狀有著較高的要求，因此有限單元法在處理畸變網(wǎng)格時(shí)表現(xiàn)無(wú)力。

為了解決真實(shí)結(jié)構(gòu)中存在的疲勞裂紋擴(kuò)展問(wèn)題，針對(duì)斷裂問(wèn)題的數(shù)值分析方法得到了廣泛的研究。但是采用常規(guī)有限元法進(jìn)行斷裂力學(xué)分析時(shí)，由于其裂紋尖端存在奇異性，計(jì)算易不收斂，必須對(duì)裂尖網(wǎng)格進(jìn)行細(xì)化加密處理。另外，在裂紋擴(kuò)展計(jì)算中，裂紋必須沿著已有的網(wǎng)格邊線(xiàn)進(jìn)行擴(kuò)展，每一步的擴(kuò)展都需要進(jìn)行網(wǎng)格重構(gòu)，增加了計(jì)算的負(fù)擔(dān)。所以，為了降低對(duì)網(wǎng)格的依賴(lài)性，引入了放寬網(wǎng)格負(fù)擔(dān)的離散化技術(shù)，擴(kuò)展有限元法XFEM[1-6]應(yīng)運(yùn)而生。XFEM基于單元分解法PUM[7,8]，通過(guò)在位移場(chǎng)的基礎(chǔ)上添加擴(kuò)充函數(shù)來(lái)反映裂尖區(qū)域的不連續(xù)性。XFEM的不連續(xù)界面獨(dú)立于網(wǎng)格，可以沿著任意角度擴(kuò)展且無(wú)需網(wǎng)格加密和重構(gòu)，有效彌補(bǔ)了常規(guī)有限元法在處理斷裂力學(xué)方面的不足，但XFEM仍然需要進(jìn)行坐標(biāo)映射。Liu等[9]提出了S -FEM，二維 S -FEM 將單元內(nèi)部的面積分轉(zhuǎn)換為沿著光滑域的線(xiàn)積分，避免了單元的坐標(biāo)映射(等參轉(zhuǎn)換)和雅可比矩陣計(jì)算等問(wèn)題，降低了常規(guī)有限元法對(duì)網(wǎng)格質(zhì)量的依賴(lài)性。隨后根據(jù)光滑應(yīng)變技術(shù)的不同演化出基于節(jié)點(diǎn)光滑有限元法NS -FEM、光滑子域有限元法CS -FEM、光滑邊域有限元法ES -FEM和光滑面域有限元法FS -FEM。

Bordas等[10]將光滑應(yīng)變技術(shù)CS -FEM應(yīng)用于XFEM中，通過(guò)圍線(xiàn)積分，簡(jiǎn)化了不連續(xù)近似積分，避免了在近似中引入Westergaard solution時(shí)對(duì)奇異函數(shù)進(jìn)行數(shù)值積分的需要，其穩(wěn)定性、魯棒性和計(jì)算成本得到大幅提高。Chen等[11,12]提出將ES -FEM引入到XFEM計(jì)算中，將內(nèi)部積分轉(zhuǎn)化為邊界積分，消除了對(duì)被不連續(xù)切割單元進(jìn)行細(xì)分的需要，簡(jiǎn)化了數(shù)值積分過(guò)程，在收斂性、精度和計(jì)算效率方面優(yōu)于標(biāo)準(zhǔn)的XFEM。Vu-Bac等[13]提出NS -XFEM，與ES -XFEM和標(biāo)準(zhǔn)XFEM-T3對(duì)比發(fā)現(xiàn)，NS -XFEM可以產(chǎn)生超收斂解。Jiang等[14]將ES -XFEM應(yīng)用于正交異性復(fù)合材料的斷裂機(jī)理分析中，再次證明ES -XFEM比XFEM更準(zhǔn)確。Wu等[15]證明了ES -XFEM在二維彈性體動(dòng)態(tài)斷裂問(wèn)題分析中的有效性和效率。Surendran等[5]針對(duì)非多項(xiàng)式擴(kuò)充函數(shù)提出了線(xiàn)性光滑技術(shù)(Linear smoothing)來(lái)近似模擬高階多項(xiàng)式場(chǎng)。Wan等[16]針對(duì)軸對(duì)稱(chēng)弱不連續(xù)問(wèn)題，提出了一種完全光滑擴(kuò)展有限元方法，并在彈性靜力學(xué)和彈性動(dòng)力學(xué)中證明了方法的可行性。

上述文獻(xiàn)只是單獨(dú)將某一種光滑有限元技術(shù)引入到XFEM中進(jìn)行計(jì)算，本文在研究XFEM的基礎(chǔ)上結(jié)合光滑子域和光滑邊域兩種方法的優(yōu)點(diǎn)，提出了一種新的有限元算法，在單元與擴(kuò)充結(jié)點(diǎn)的選取上采用ES -FEM的光滑域劃分思想，在數(shù)值積分計(jì)算剛度矩陣時(shí)采用基于三角形子域的 CS -FEM 積分思路。通過(guò)拉伸載荷下的純Ⅰ型邊界直裂紋模型，確定S -XFEM中裂尖光滑單元子胞數(shù)目、高斯點(diǎn)個(gè)數(shù)、交互積分區(qū)域放大系數(shù)及擴(kuò)充方式等關(guān)鍵參數(shù)的選取。通過(guò)剪切載荷下純II型邊界直裂紋模型和拉伸載荷下混合邊界直裂紋模型，對(duì)S -XFEM的實(shí)用性和精確性進(jìn)行研究。

2 S -XFEM基本理論

2.1 光滑應(yīng)變技術(shù)

由常規(guī)有限元法可以得到單元內(nèi)相容應(yīng)變場(chǎng)的計(jì)算公式為[17]

(1)

(2)

利用光滑應(yīng)變技術(shù)創(chuàng)建光滑應(yīng)變場(chǎng)，光滑應(yīng)變技術(shù)將光滑域內(nèi)的應(yīng)變場(chǎng)設(shè)定為常量，其大小可通過(guò)修正相容應(yīng)變場(chǎng)獲得，也可通過(guò)直接對(duì)位移場(chǎng)求導(dǎo)得

(3)

(4)

并且光滑函數(shù)必須滿(mǎn)足

(5)

2.2 基于CS -FEM的S -XFEM剛度矩陣

(6)

式(6)擴(kuò)充項(xiàng)的位移值可通過(guò)位移擴(kuò)充函數(shù)使位移在結(jié)點(diǎn)處的值等于結(jié)點(diǎn)位移，本文選用一種經(jīng)典的修正方案,

(7)

式中 第一個(gè)擴(kuò)充項(xiàng)含有一個(gè)階躍函數(shù)(Heaviside jump function)H(x)，在裂紋的上方取+1，在裂紋的下方取-1，其數(shù)學(xué)表達(dá)式為

(8)

式中x為問(wèn)題點(diǎn)(高斯點(diǎn))，x*為裂紋線(xiàn)上離x最近的一個(gè)點(diǎn)，n為裂紋在x*點(diǎn)處的單位外法線(xiàn)向量。

(9)

式中 (r,θ)為裂尖處的局部極坐標(biāo)系。

(10)

(10)

(11)

(11)

(12)

(12)

所以系統(tǒng)總體剛度矩陣可以表達(dá)為

(13)

式中Ns為光滑域個(gè)數(shù)。

(14)

2.3 基于ES -FEM的S -XFEM擴(kuò)充結(jié)點(diǎn)選擇

如圖1(a)所示，含有裂尖的光滑域用網(wǎng)格線(xiàn)標(biāo)識(shí)，稱(chēng)為裂尖光滑域，該光滑域所在的兩個(gè)單元的所有結(jié)點(diǎn)為裂尖結(jié)點(diǎn)Nes -asympt，用圓圈標(biāo)識(shí)。裂紋貫穿的光滑域的支持結(jié)點(diǎn)為貫穿結(jié)點(diǎn)Nes -disc，用三角標(biāo)識(shí)?？梢园l(fā)現(xiàn)，相比于XFEM的結(jié)點(diǎn)類(lèi)型選擇，基于ES -FEM的選擇方案可以增加裂尖結(jié)點(diǎn)的個(gè)數(shù)。

2.4 S -XFEM的積分策略

CS -FEM的光滑子域完全在單元內(nèi)部，所以其單元?jiǎng)偠染仃嚪e分域與FEM中所使用的積分域一致，避免對(duì)ES -FEM的一個(gè)光滑域中分屬不同單元的光滑子域的形函數(shù)的單獨(dú)計(jì)算，確保一個(gè)光滑域內(nèi)所有高斯點(diǎn)的形函數(shù)可以在同一個(gè)單元中計(jì)算。所以S -XFEM與ES -FEM的最大區(qū)別在于進(jìn)行計(jì)算的光滑域的不同，如圖1所示。圖1(b)為基于S -XFEM的光滑單元的劃分方案。與XFEM相似，S -XFEM也有5種積分區(qū)域，即裂尖光滑單元、貫穿光滑單元、裂尖混合光滑單元、貫穿混合光滑單元和常規(guī)光滑單元[11,12]。

裂尖光滑單元，由于被積函數(shù)為非多項(xiàng)式，簡(jiǎn)單地將光滑單元?jiǎng)澐譃槎噙呅巫佑蜻M(jìn)行計(jì)算無(wú)法達(dá)到數(shù)值精度要求，所以需要在裂尖附近使用更多的積分點(diǎn)。通過(guò)以下步驟實(shí)現(xiàn)。(1) 使用XFEM中Delaunay三角化將光滑單元?jiǎng)澐譃槿切巫佑騭ub -sd1～sub -sd5； (2) 將三角形子域進(jìn)一步劃分為ns c個(gè)子胞，圖2(a)顯示了將sub -sd1和sub -sd2劃分為3個(gè)子胞的情況，sub -sd1和sub -sd2分別劃分為sc1～sc3和sc4～sc6； (3) 在每個(gè)子胞邊界上進(jìn)行積分。當(dāng)在共用一條裂紋的兩個(gè)子胞(sc3和sc5)上布置高斯點(diǎn)進(jìn)行積分時(shí)，階躍函數(shù)和裂尖擴(kuò)充函數(shù)雖然在兩個(gè)子胞上的計(jì)算值相同，但在裂紋處的位移實(shí)際上卻是不連續(xù)的。所以在計(jì)算時(shí)不使用高斯點(diǎn)坐標(biāo)，而使用子胞形心坐標(biāo)。與常規(guī)光滑有限元相比，由于裂尖附近的梯度很大，所以需要在邊界上使用更多的高斯積分點(diǎn)，圖示每個(gè)子胞邊界上使用5個(gè)高斯積分點(diǎn)。

圖1 基于ES -FEM和S -XFEM的擴(kuò)充光滑域和結(jié)點(diǎn)選擇方案

貫穿光滑單元，將裂紋貫穿的光滑單元沿著裂紋線(xiàn)劃分為若干個(gè)子域，使得每個(gè)子域上被積函數(shù)連續(xù)可微。與XFEM類(lèi)似，在S -XFEM中同樣需要將非三角形多邊形使用Delaunay三角化為三角形子域，如圖2(b)所示，貫穿單元?jiǎng)澐譃閟ub -sd1～sub -sd3三部分。線(xiàn)性形函數(shù)與階躍函數(shù)的乘積NiH在外部邊界和內(nèi)部邊線(xiàn)上是線(xiàn)性的，所以在貫穿光滑單元的光滑子域邊界上使用一個(gè)高斯點(diǎn)即可。

圖2 裂尖單元和貫穿單元的子域劃分方式

裂尖混合光滑單元和貫穿混合光滑單元使用的積分策略和文獻(xiàn)[11,12]的一樣，本文不再?gòu)?fù)述。

3 數(shù)值分析結(jié)果

3.1 拉伸載荷下的邊界直裂紋

采用一個(gè)承受拉伸載荷作用的邊界直裂紋算例[11,12]研究上述算法。材料參數(shù)為楊氏模量E=30 GPa，泊松比υ=0.3，假設(shè)為平面應(yīng)變問(wèn)題。

平板的尺寸為1 m×2 m，邊界直裂紋長(zhǎng)度為a，在上邊界施加σ=1.0 Pa的拉伸應(yīng)力。平板底部邊界施加y方向上的位移約束，底部左端點(diǎn)同時(shí)施加x和y方向的位移約束。模型的幾何、載荷和邊界條件如圖3(a)所示，T3和Q4網(wǎng)格劃分如圖3(b，c)所示，圖示網(wǎng)格布種個(gè)數(shù)為20×40。

文獻(xiàn)[18]給出了該模型應(yīng)力強(qiáng)度因子的精確解：

(15)

式中C為幾何修正系數(shù)：

C=1.12-0.231(a/b)+10.55(a/b)2-

21.72(a/b)3+30.39(a/b)4

(16)

3.1.1 S -XFEM裂尖光滑單元光滑子胞數(shù)目

和高斯點(diǎn)數(shù)目的影響研究

關(guān)于高斯點(diǎn)個(gè)數(shù)對(duì)應(yīng)力強(qiáng)度因子的影響已經(jīng)在文獻(xiàn)[11,12]給出，因此，本文不再進(jìn)行討論，直接引用文獻(xiàn)的結(jié)論，選取光滑子胞邊界高斯點(diǎn)數(shù)目ngau=5。

CS -FEM在計(jì)算靜力學(xué)問(wèn)題時(shí)其精度會(huì)隨著子胞個(gè)數(shù)的改變而發(fā)生較大的變化，當(dāng)子胞個(gè)數(shù)ns c=4時(shí)計(jì)算結(jié)果最接近于解析解[19]。由于本文擴(kuò)充函數(shù)的引入，子胞個(gè)數(shù)對(duì)于結(jié)果的影響不能直接引用已有結(jié)論，下面對(duì)光滑子胞個(gè)數(shù)對(duì)裂尖光滑單元積分結(jié)果的影響進(jìn)行討論。

圖3 拉伸載荷下邊界直裂紋模型及網(wǎng)格劃分

采用本節(jié)計(jì)算模型，取a=0.6 m，劃分70×140個(gè)結(jié)點(diǎn)，裂尖采用拓?fù)鋽U(kuò)充，交互積分所使用的半徑放大系數(shù)rk=5。圖4為Q4裂尖單元內(nèi)使用不同子胞個(gè)數(shù)時(shí)的劃分結(jié)果。表1為S -XFEM使用T3和Q4單元計(jì)算的應(yīng)力強(qiáng)度因子隨子胞個(gè)數(shù)的變化情況。由表1可知，KI的計(jì)算值始終小于解析解，隨著ns c的增加，KI越來(lái)越小并逐漸收斂于固定值。由于ns c=1時(shí),KI與最終收斂的固定值相差較大，計(jì)算結(jié)果存在不穩(wěn)定性。所以綜合考慮，本文認(rèn)為ns c=2時(shí)，S -XFEM可以取得較為穩(wěn)定且精確的數(shù)值解。

因此，本文如不特殊說(shuō)明，選取ngau=5，ns c=2。

3.1.2 交互積分區(qū)域半徑放大系數(shù)的影響研究

使用交互積分法計(jì)算應(yīng)力強(qiáng)度因子，需要確定一個(gè)特征長(zhǎng)度的放大系數(shù)rk。采用3.1.1的模型，計(jì)算結(jié)果列入表2。當(dāng)rk很小時(shí)計(jì)算結(jié)果波動(dòng)很大，這是因?yàn)榻换シe分區(qū)域靠近裂尖，裂尖擴(kuò)充結(jié)點(diǎn)使計(jì)算結(jié)果分散性很大，隨著rk的逐漸增大，積分區(qū)域離裂尖越來(lái)越遠(yuǎn)，計(jì)算結(jié)果趨于穩(wěn)定，當(dāng)rk=5時(shí)，4種方法的結(jié)果均收斂于穩(wěn)定值。因此，本文如不特殊說(shuō)明，取交互積分區(qū)域半徑放大系數(shù)rk=5。

表1 不同ns c下的KI計(jì)算值和解析解

圖4 裂尖Q4單元子胞劃分

3.1.3 擴(kuò)充方式的影響研究

基于二維幾何的裂尖擴(kuò)充，以裂尖為圓心定義一個(gè)半徑為r的圓，在此區(qū)域內(nèi)的結(jié)點(diǎn)都用裂尖擴(kuò)充函數(shù)擴(kuò)充。半徑r的定義為r=rehlocal，re為幾何擴(kuò)充半徑的放大系數(shù)。拓?fù)鋽U(kuò)充只擴(kuò)充含有裂尖的單元結(jié)點(diǎn)。

采用3.1.1的模型，圖5給出了Q4使用不同擴(kuò)充方式下的結(jié)點(diǎn)擴(kuò)充示意圖，圖中方框?yàn)榱鸭鈹U(kuò)充結(jié)點(diǎn)，叉形為貫穿擴(kuò)充結(jié)點(diǎn)。不同擴(kuò)充方式下應(yīng)力強(qiáng)度因子的計(jì)算結(jié)果列入表3。可以看出，re=1時(shí)拓?fù)鋽U(kuò)充與幾何擴(kuò)充的計(jì)算結(jié)果相同，因?yàn)檫@時(shí)幾何擴(kuò)充半徑內(nèi)的結(jié)點(diǎn)與拓?fù)鋽U(kuò)充一致。隨著re的增大，KI計(jì)算值增大并趨近于解析解，這是因?yàn)閞e的增大使得擴(kuò)充結(jié)點(diǎn)增加，裂尖擴(kuò)充函數(shù)作用更加明顯，相應(yīng)的計(jì)算結(jié)果也就更加精確。當(dāng)re進(jìn)一步增加時(shí)，計(jì)算誤差突然增大，這是因?yàn)樵诮换シe分半徑放大系數(shù)(rk=5)一定的情況下，裂尖擴(kuò)充結(jié)點(diǎn)范圍的增大使得交互積分區(qū)域內(nèi)包含了裂尖結(jié)點(diǎn)，從而降低了KI的計(jì)算精度。另外，由于ES -FEM結(jié)點(diǎn)的支持域比XFEM大，所以相比于XFEM，其對(duì)幾何擴(kuò)充半徑放大系數(shù)re更為敏感。因此面對(duì)不同需求時(shí)可以選擇不同的單元擴(kuò)充方式。

圖5 不同擴(kuò)充方式下Q4結(jié)點(diǎn)擴(kuò)充

3.1.4 精確性與收斂性研究

對(duì)比S -XFEM和XFEM在分別使用T3和Q4單元時(shí)應(yīng)力強(qiáng)度因子的計(jì)算結(jié)果，以研究不同算法的精確度與收斂性。采用3.1.1的模型，裂紋長(zhǎng)度選為0.3 m和0.6 m。圖6和圖7分別給出了4種方法應(yīng)力強(qiáng)度因子計(jì)算值隨著結(jié)點(diǎn)數(shù)增加的變化趨勢(shì)和其相對(duì)誤差隨結(jié)點(diǎn)數(shù)增加的變化趨勢(shì)。

表2 不同交互積分半徑放大系數(shù)下的KI計(jì)算值Tab.2 Stress intensity factor under different rk

表3 不同擴(kuò)充方式下的KI計(jì)算值Tab.3 Stress intensity factor under different extended modes

圖6 結(jié)點(diǎn)數(shù)對(duì)應(yīng)力強(qiáng)度因子計(jì)算值的影響

圖7 結(jié)點(diǎn)數(shù)對(duì)應(yīng)力強(qiáng)度因子計(jì)算值相對(duì)誤差的影響

從圖6可以看出，當(dāng)裂紋長(zhǎng)度a=0.3 m和 0.6 m 時(shí)，四種方法計(jì)算的KI隨結(jié)點(diǎn)數(shù)的增加均收斂于解析解，當(dāng)結(jié)點(diǎn)數(shù)為 70×140 時(shí)，所有方法的相對(duì)誤差都控制在5%左右，當(dāng)結(jié)點(diǎn)數(shù)為100×200時(shí)，XFEM-Q4方法的相對(duì)誤差在兩種裂紋長(zhǎng)度下分別可以達(dá)到0.309%和0.508%。通過(guò)比較可以看出，XFEM的計(jì)算精度和收斂速度稍?xún)?yōu)于 S -XFEM，相同算法下，Q4單元的計(jì)算精度和收斂速度稍?xún)?yōu)于T3單元。

3.1.5 不規(guī)則網(wǎng)格的影響研究

采用不規(guī)則網(wǎng)格檢驗(yàn)S -XFEM程序?qū)W(wǎng)格的依賴(lài)性，不規(guī)則網(wǎng)格采用式(17)隨機(jī)生成：

(17)

式中x和y為初始規(guī)則網(wǎng)格結(jié)點(diǎn)的坐標(biāo)值，Δx和Δy為單元直角邊在x和y方向的長(zhǎng)度，rd為在[-1,1]內(nèi)變化的隨機(jī)數(shù)，α為單元奇異系數(shù)。

通過(guò)改變?chǔ)林悼梢陨刹煌潭鹊牟灰?guī)則網(wǎng)格，如圖8所示，采用3.1.1的模型。逐漸增加不規(guī)則因子α以計(jì)算不同網(wǎng)格畸變程度下的KI，計(jì)算結(jié)果相對(duì)誤差如圖9所示?？梢钥闯?，S -XFEM的精確度雖然沒(méi)有XFEM高，但S -XFEM在不規(guī)則網(wǎng)格下的穩(wěn)定性要明顯優(yōu)于XFEM，其計(jì)算結(jié)果幾乎不受網(wǎng)格畸變程度的影響，而XFEM計(jì)算結(jié)果波動(dòng)較大。在實(shí)際的計(jì)算中，當(dāng)α=0.12時(shí)，XFEM-T3出現(xiàn)了三點(diǎn)幾乎共線(xiàn)的細(xì)長(zhǎng)單元，凸包計(jì)算無(wú)法完成，這直接導(dǎo)致了程序的崩潰。S -XFEM對(duì)畸變網(wǎng)格的特性將十分有助于分析大變形問(wèn)題。

3.2 剪切載荷下的邊界直裂紋

使用一個(gè)承受剪切載荷的邊界直裂紋算例對(duì)算法的精確性和收斂性進(jìn)行進(jìn)一步研究。材料參數(shù)為楊氏模量E=30 GPa，泊松比υ=0.3，假設(shè)為平面應(yīng)變問(wèn)題。平板的尺寸為16 mm×7 mm，邊界直裂紋長(zhǎng)度a=3.55 mm，在上邊界施加大小為τ=1.0 Pa的剪切應(yīng)力。平板底部邊界施加y方向上的位移約束，底部左端點(diǎn)同時(shí)施加x和y方向的位移約束。模型的幾何、載荷、邊界條件以及T3和Q4網(wǎng)格(結(jié)點(diǎn)數(shù)20×40)如圖10所示。文獻(xiàn)[18]給出了該問(wèn)題的應(yīng)力強(qiáng)度因子參考解：

(18)

裂尖單元采用拓?fù)鋽U(kuò)充，ngau=5，ns c=2，rk=5。圖11和圖12分別給出了隨著網(wǎng)格數(shù)的增加，四種方法計(jì)算的KI、KII和相對(duì)誤差的變化規(guī)律。可以看出，這四種方法計(jì)算的KI和KII隨著結(jié)構(gòu)結(jié)點(diǎn)數(shù)的增加，均能收斂于解析解。當(dāng)結(jié)點(diǎn)數(shù)為 60×120時(shí)，所有方法的相對(duì)誤差都控制在5%左右；當(dāng)結(jié)點(diǎn)數(shù)為100×200時(shí)，XFEM-Q4方法的相對(duì)誤差分別為0.494%和1.397%。通過(guò)比較可以看出，XFEM的計(jì)算精度和收斂速度稍好于S -XFEM，相同算法下，Q4單元的計(jì)算精度和收斂速度稍好于T3單元。對(duì)比3.1節(jié)的結(jié)果，可以發(fā)現(xiàn)S -XFEM方法在計(jì)算精度上與XFEM的差距有大幅減小，其數(shù)值計(jì)算精度幾乎與XFEM相當(dāng)。

圖8 不同α值下的三角形和四邊形網(wǎng)格

圖9 單元奇異系數(shù)α對(duì)相對(duì)誤差ek的影響

圖10 剪切載荷下邊界直裂紋模型及網(wǎng)格劃分

3.3 拉伸載荷下的邊界斜裂紋擴(kuò)展問(wèn)題

使用一個(gè)承受拉伸載荷作用的邊界斜裂紋擴(kuò)展算例驗(yàn)證本文S -XFEM程序在處理連續(xù)擴(kuò)展時(shí)的特性。材料參數(shù)為楊氏模量E=30 GPa，泊松比υ=0.3，假設(shè)為平面應(yīng)變問(wèn)題。平板的尺寸為1 m×2 m，含有與水平成30°夾角的邊界斜裂紋，裂紋初始坐標(biāo)為(0,0)和(0.3，0.173205)。在上邊界施加大小為σ=1.0 Pa的拉伸應(yīng)力。平板底部邊界施加y方向上的位移約束，底部左端點(diǎn)同時(shí)施加x和y方向的位移約束。模型的幾何、載荷和邊界條件如圖13(a)所示，設(shè)定擴(kuò)展步長(zhǎng)為0.03 m，擴(kuò)展步數(shù)為10。圖13(b，c)為該模型使用三角形網(wǎng)格和四邊形網(wǎng)格劃分示意圖(圖示布種為20×40)，實(shí)際計(jì)算使用布種數(shù)為100×200。

圖11 結(jié)點(diǎn)數(shù)對(duì)應(yīng)力強(qiáng)度因子計(jì)算值的影響

圖14給出四種方法計(jì)算裂紋擴(kuò)展過(guò)程中的應(yīng)力強(qiáng)度因子KI和KII的值?？梢钥闯?，隨著裂紋的向前擴(kuò)展，KI值均逐漸增大，KII值均逐漸收斂于0，裂紋保持水平擴(kuò)展，裂紋轉(zhuǎn)換為純II型。

圖12 結(jié)點(diǎn)數(shù)對(duì)應(yīng)力強(qiáng)度因子計(jì)算值相對(duì)誤差的影響

圖13 拉伸載荷下邊界斜裂紋擴(kuò)展模型及網(wǎng)格劃分

4 結(jié) 論

本文主要研究了求解二維平面問(wèn)題的S -XFEM，詳細(xì)闡述了其基本理論并推導(dǎo)數(shù)值積分公式，給出了光滑域內(nèi)高斯點(diǎn)的積分策略，設(shè)計(jì)并編制了算法的有限元計(jì)算程序，采用經(jīng)典算例驗(yàn)證并與XFEM對(duì)比。光滑應(yīng)變技術(shù)的運(yùn)用使得在計(jì)算單元?jiǎng)偠染仃嚨倪^(guò)程中無(wú)需對(duì)各點(diǎn)的形函數(shù)進(jìn)行求導(dǎo)，所以在S -XFEM中也無(wú)需對(duì)擴(kuò)充項(xiàng)求偏導(dǎo)，從而避免了對(duì)極坐標(biāo)函數(shù)求偏導(dǎo)數(shù)的復(fù)雜過(guò)程，有效地簡(jiǎn)化了單元?jiǎng)偠染仃嚨挠?jì)算過(guò)程。主要的研究?jī)?nèi)容和結(jié)論如下。

(1) 將子域光滑應(yīng)變技術(shù)和邊域光滑應(yīng)變技術(shù)擴(kuò)展到XFEM的框架中，建立了一種新的 S -XFEM，在單元及擴(kuò)充結(jié)點(diǎn)選取時(shí)采用ES -FEM的光滑域劃分方式，在數(shù)值積分計(jì)算剛度矩陣時(shí)采用基于三角形子域的CS -FEM積分思路。該方法通過(guò)光滑域?qū)?nèi)部積分轉(zhuǎn)化為邊界積分，消除了對(duì)受不連續(xù)結(jié)構(gòu)切割單元進(jìn)行細(xì)分的需要。而且，通過(guò)應(yīng)變平滑技術(shù)，計(jì)算剛度矩陣時(shí)不再需要對(duì)形函數(shù)求導(dǎo)，避免了裂尖奇異點(diǎn)的出現(xiàn)，極大簡(jiǎn)化了計(jì)算過(guò)程。

(2) 計(jì)算結(jié)果表明XFEM和S -XFEM均有很高的計(jì)算精度和穩(wěn)定的收斂性，XFEM略?xún)?yōu)于S -XFEM，而S -XFEM對(duì)于不規(guī)則網(wǎng)格的適應(yīng)性要明顯強(qiáng)于XFEM，可以很好地應(yīng)用于大變形問(wèn)題分析，甚至是動(dòng)態(tài)分析。本文提出的S -XFEM對(duì)于純II型裂紋、裂紋擴(kuò)展模擬皆有很好的計(jì)算結(jié)果，說(shuō)明本文編制的S -XFEM程序是有效的。