周 博 薛世峰

(中國石油大學(xué)儲運(yùn)與建筑工程學(xué)院，山東青島266580)

基于擴(kuò)展有限元的應(yīng)力強(qiáng)度因子的位移外推法1)

周 博2)薛世峰

(中國石油大學(xué)儲運(yùn)與建筑工程學(xué)院，山東青島266580)

針對平面裂紋問題，闡述了擴(kuò)展有限元法的單元位移模式、推導(dǎo)了擴(kuò)展有限元法的控制方程、介紹了特殊單元的數(shù)值積分技術(shù).基于最小二乘法，建立了應(yīng)力強(qiáng)度因子位移外推法的計算公式.利用MATLAB編寫計算程序，對平面裂紋問題用擴(kuò)展有限元法進(jìn)行了計算.基于擴(kuò)展有限元法的計算結(jié)果，分別利用位移外推法和相互作用積分法，對平面裂紋的應(yīng)力強(qiáng)度因子進(jìn)行了計算.計算結(jié)果表明,位移外推法比相互作用積分法能更方便和準(zhǔn)確地計算平面裂紋的應(yīng)力強(qiáng)度因子.

擴(kuò)展有限元法，應(yīng)力強(qiáng)度因子，位移外推法，相互作用積分法

在工程實際中經(jīng)常發(fā)生斷裂引起的破壞事故，如地震所致地質(zhì)構(gòu)造和工程結(jié)構(gòu)的開裂、壓力管道的裂紋失穩(wěn)擴(kuò)展、機(jī)械構(gòu)件的疲勞斷裂等，這些事故均可以對人民生命和財產(chǎn)造成重大損失.裂紋的萌生是不可避免和難以量化研究的，確定裂紋產(chǎn)生后是否繼續(xù)擴(kuò)展或發(fā)生止裂的斷裂力學(xué)具有理論和實際意義[1].有限元法理論基礎(chǔ)成熟且通用性強(qiáng)，是處理斷裂力學(xué)問題的常用數(shù)值方法.在應(yīng)用有限元法處理裂紋時，需要將裂紋設(shè)置為單元的邊，將裂紋尖端設(shè)置為單元的節(jié)點(diǎn).為提高計算精度在裂紋尖端處還需使用奇異單元和進(jìn)行網(wǎng)格加密，在模擬裂紋擴(kuò)展時還需要重新劃分網(wǎng)格.因此有限元法處理斷裂力學(xué)問題的過程，相對繁瑣和困難，遇到復(fù)雜問題時就可能無能為力.

1999–2002年，Belytschko等[24]，在傳統(tǒng)有限元法的框架內(nèi)建立了擴(kuò)展有限元法，很好地解決了由于材料或幾何等因素引起的不連續(xù)問題，特別適合于處理斷裂力學(xué)問題，近年來正不斷得到成功應(yīng)用.2005年，余天堂[5]詳細(xì)介紹了擴(kuò)展有限元法求解裂紋體問題的基本過程，利用相互作用積分法計算了I型裂紋的應(yīng)力強(qiáng)度因子，取得了良好的計算結(jié)果.2007年，方修君等[67]基于擴(kuò)展有限元法基本原理，開發(fā)了ABAQUS調(diào)用的黏聚裂紋模型，對三點(diǎn)彎曲梁的裂紋擴(kuò)展過程進(jìn)行了模擬，結(jié)果表明擴(kuò)展有限元法對非連續(xù)位移場的求解不依賴于單元邊界，是模擬裂紋擴(kuò)展的有效方法.2008年，懂玉文等[8]提出了通過求解擴(kuò)展有限元法總體剛度方程中的附加自由度，直接計算應(yīng)力強(qiáng)度因子的計算方法，取得了良好數(shù)值計算結(jié)果.2009年，宋娜等[9]研究了單元類型、網(wǎng)格密度、J積分區(qū)域半徑等因素對擴(kuò)展有限元法計算應(yīng)力強(qiáng)度因子的精度的影響，結(jié)果表明利用擴(kuò)展有限元法計算應(yīng)力強(qiáng)度因子時，上述因素對計算結(jié)果精度的影響不顯著.2010年，余天堂[10]建立了閉合型摩擦裂紋問題的擴(kuò)展有限元法線性互補(bǔ)模型，將裂紋面非線性摩擦接觸轉(zhuǎn)化為一個無需迭代求解的線性互補(bǔ)問題，取得了良好的計算效果.

2011年，茹忠亮等[11]分析了應(yīng)力強(qiáng)度因子J積分方法中積分區(qū)域、網(wǎng)格密度等因素對擴(kuò)展有限元法計算精度的影響，給出了計算應(yīng)力強(qiáng)度因子的合適參數(shù)，有效提高了應(yīng)力強(qiáng)度因子的計算精度. 2012年，杜修力等[12]利用擴(kuò)展有限元法模擬了濕篩混凝土單軸拉伸作用下混凝土板的細(xì)觀斷裂破壞過程，分析了混凝土裂紋萌生、擴(kuò)展的過程及破壞形態(tài).2013年，茹忠亮等[13]建立了預(yù)置裂紋混凝土梁的三維擴(kuò)展有限元模型，對鋼筋混凝土梁的復(fù)合斷裂過程進(jìn)行了模擬分析，取得了良好的計算結(jié)果.2014年，傅玉珍等[14]利用擴(kuò)展有限元法模擬了沖擊載荷作用下混凝土巴西圓盤的裂紋擴(kuò)展過程，分析了應(yīng)力場分布規(guī)律和網(wǎng)格形態(tài)對計算精度的影響；師訪等[15]基于擴(kuò)展有限元法，研究了正交各向異性巖體的應(yīng)力強(qiáng)度因子和裂紋擴(kuò)展特性，模擬了四點(diǎn)彎曲正交各向異性巖體式樣的裂紋擴(kuò)展過程，詳細(xì)分析了彈性模量比對裂紋擴(kuò)展行為的影響；盛茂等[16]建立了水力壓裂問題的擴(kuò)展有限元法，將縫內(nèi)水壓轉(zhuǎn)化為單元等效節(jié)點(diǎn)力，利用相互作用積分法計算縫尖的應(yīng)力強(qiáng)度因子，采用最大能量釋放率準(zhǔn)則作為裂紋擴(kuò)展判據(jù)，模擬了水力壓裂裂紋擴(kuò)展路徑.2015年，江守燕等[17]基于擴(kuò)展有限元法，建立了求解雙材料界面裂紋應(yīng)力強(qiáng)度因子的相互作用積分法，對雙材料界面裂紋的應(yīng)力強(qiáng)度因子進(jìn)行了計算.

計算斷裂參數(shù)是斷裂力學(xué)的重要內(nèi)容，應(yīng)力強(qiáng)度因子、J積分、應(yīng)變能釋放率是斷裂力學(xué)中三個最基本的斷裂參數(shù).應(yīng)力強(qiáng)度因子反映了彈性裂紋尖端應(yīng)力場的強(qiáng)弱，J積分描述了由于裂紋存在所吸收的能量，應(yīng)變能釋放率描述的是產(chǎn)生裂紋面所需要的能量.對于線彈性材料，這三個斷裂參數(shù)可通過材料參數(shù)聯(lián)系起來，且J積分和應(yīng)變能釋放率是等價的.本文針對彈性平面裂紋問題，闡述了擴(kuò)展有限元法平面4節(jié)點(diǎn)等參單元的位移模式，推導(dǎo)了平面裂紋的擴(kuò)展有限元法控制方程，介紹了特殊單元的數(shù)值積分方法.基于最小二乘法和擴(kuò)展有限元法的計算結(jié)果，建立了應(yīng)力強(qiáng)度因子位移外推法的計算公式.利用MATLAB編寫相關(guān)計算程序，對平面裂紋問題用擴(kuò)展有限元法進(jìn)行了計算.基于擴(kuò)展有限元法的計算結(jié)果，分別利用位移外推法和相互作用積分法，對平面裂紋的應(yīng)力強(qiáng)度因子進(jìn)行了計算.計算結(jié)果表明，本文建立的基于擴(kuò)展有限元法的應(yīng)力強(qiáng)度因子的位移外推法，能方便和準(zhǔn)確地計算平面裂紋的應(yīng)力強(qiáng)度因子.

1 擴(kuò)展有限元法

1.1 位移模式

擴(kuò)展有限元法的基本原理是基于單位分解思想，在傳統(tǒng)有限元法的位移模式中加進(jìn)一些加強(qiáng)函數(shù)以反映特殊單元位移場的不連續(xù)性.擴(kuò)展有限元法在分析裂紋問題時，計算網(wǎng)格和裂紋是相互獨(dú)立的，比傳統(tǒng)有限元法處理裂紋問題方便很多.在擴(kuò)展有限元法的具體計算過程中，通過水平集函數(shù)來識別和追蹤裂紋及其擴(kuò)展路徑.

對于平面裂紋問題，擴(kuò)展有限元法的建模包括兩部分：一是不考慮裂紋的幾何形狀和位置，對計算區(qū)域進(jìn)行網(wǎng)格劃分；二是對和裂紋接觸的單元采用特殊的位移模式，改進(jìn)裂紋附近位移場的數(shù)值精度.

如圖1所示含裂紋平面結(jié)構(gòu)的計算網(wǎng)格，其中的單元類型有三類：(1)不接觸裂紋的常規(guī)單元，節(jié)點(diǎn)無特殊標(biāo)記；(2)裂紋貫穿單元，節(jié)點(diǎn)用方塊標(biāo)記；(3)包含裂紋尖端單元，節(jié)點(diǎn)用圓圈標(biāo)記.對于不接觸裂紋的常規(guī)單元，其位移模式的選擇和有限元法位移模式的選擇完全相同，即

其中,Ni為單元形函數(shù)，ui為節(jié)點(diǎn)位移分量.

圖1 含裂平面結(jié)構(gòu)的計算網(wǎng)格

對于裂紋貫穿單元，單元位移模式[4]為

式中，A+和A-分別代表如圖2所示的裂紋貫穿單元內(nèi)位于裂紋兩側(cè)的區(qū)域.

圖2 裂紋貫穿單元的區(qū)域劃分

對于包含裂紋尖端單元，單元位移模式[4]為

其中r和θ為如圖3所示裂紋尖端局部極坐標(biāo)系內(nèi)的極坐標(biāo).

圖3 包含裂紋尖端單元的裂紋尖端局部坐標(biāo)系

1.2 控制方程

和傳統(tǒng)有限元法控制方程的推導(dǎo)類似，擴(kuò)展有限元法的控制方程也可以根據(jù)虛功原理得到.設(shè)含裂紋體結(jié)構(gòu)的虛位移列陣為δu，由虛位移引起的虛應(yīng)變列陣為δε，根據(jù)虛位移原理得到

其中,ε為應(yīng)變列陣，F(xiàn)b為體力列陣，F(xiàn)s為面力列陣，D為彈性矩陣.對于平面應(yīng)力問題

其中E和μ分別為彈性模量和泊松比.若將式(7)中的E和μ分別用和替換，即可得到平面應(yīng)變問題的彈性矩陣.

將單元位移模式式(1)、式(2)、式(4)帶入式(6)，得到含裂紋體結(jié)構(gòu)的擴(kuò)展有限元法的整體控制方程為其中,K為結(jié)構(gòu)剛度矩陣，u′為節(jié)點(diǎn)未知量列陣，P為節(jié)點(diǎn)載荷列陣.和傳統(tǒng)有限元法不同的是，節(jié)點(diǎn)未知量列陣u′中不僅包括節(jié)點(diǎn)位移分量ui，還包括式(2)和式(4)中節(jié)點(diǎn)增強(qiáng)變量.

和傳統(tǒng)有限元法相似，擴(kuò)展有限元法的結(jié)構(gòu)剛度矩陣K，可其由單元剛度矩陣，按照對號入座的方法集合而成.不同的是在擴(kuò)展有限元法中，不同種類單元的單元剛度矩陣的規(guī)模大小是不同的.以平面四節(jié)點(diǎn)等參元為例，在擴(kuò)展有限元法中，所有單元剛度矩陣中的子矩陣可統(tǒng)一表示為

其中,non,cut和tip分別代表常規(guī)單元、裂紋貫穿單元和包含裂紋尖端單元，Bαi為幾何子矩陣.

1.3 特殊單元積分

在利用式(9)計算單元剛度矩陣時，涉及到數(shù)值積分運(yùn)算，對于平面四節(jié)點(diǎn)等參元，式(9)可以轉(zhuǎn)化為

其中,J為描述單元整體坐標(biāo) (x,y)和局部坐標(biāo)(ξ,η)間導(dǎo)數(shù)關(guān)系的雅可比矩陣的行列式.對于常規(guī)單元，積分式(10)計算和有限元法相同，采用高斯積分，取4個積分點(diǎn)就可以滿足精度要求.

對于非常規(guī)單元，積分式(10)是不連續(xù)函數(shù)的積分，可以將含裂紋的單元沿著裂紋所在位置分割成若干個四邊形子域，在每個子域進(jìn)行4個積分點(diǎn)的高斯積分，再將各子域內(nèi)的積分相加，得到整個單元的積分，具體分割方案如圖4所示.

圖4 非常規(guī)單元的分割方案

2 應(yīng)力強(qiáng)度因子的位移外推法

2.1 應(yīng)力強(qiáng)度因子

圖5為含裂紋的無限大板及其坐標(biāo)系.I型和II型裂紋應(yīng)力強(qiáng)度因子[19]分別為

和

裂紋尖端附近的位移場[19]為

其中

圖5 無限大板內(nèi)裂紋及坐標(biāo)系

2.2 位移外推法

根據(jù)式(14b)，可知當(dāng)θ=π且r→0時

因此，I型裂紋的應(yīng)力強(qiáng)度因子還可以定義為

根據(jù)式(16)，利用

根據(jù)式(16)的定義，可知式(18)中等號右端的系數(shù)C2即I型裂紋的應(yīng)力強(qiáng)度因子KI的近似值，因此可以將式(18)改寫為

根據(jù)式(17)和式(19)，可得到所有數(shù)據(jù)點(diǎn)的偏差平方和為

式(21b)即為I型裂紋應(yīng)力強(qiáng)度因子位移外推法的計算公式.

根據(jù)式(14a)，可知當(dāng)θ=π且r→0時，

因此，II型裂紋的應(yīng)力強(qiáng)度因子可以定義為

根據(jù)II型裂紋的應(yīng)力強(qiáng)度因子的定義式(31)，利用

與式(19)分析過程相似，可以設(shè)

并利用最小二乘法可以得到

式(26b)即為II型裂紋的應(yīng)力強(qiáng)度因子位移外推法的計算公式.

3 數(shù)值算例

利用MATLAB編寫計算程序，對平面裂紋問題用擴(kuò)展有限元法進(jìn)行計算.基于擴(kuò)展有限元法的計算結(jié)果，分別利用本文介紹位移外推法和文獻(xiàn)[18]介紹的相互作用積分法，對平面裂紋的應(yīng)力強(qiáng)度因子進(jìn)行計算.圖6為有中心裂紋的有限寬板，其寬和高分別為2w和2h，裂紋長度為2a，在上下兩端承受均勻拉伸應(yīng)力σ，根據(jù)對稱性，在擴(kuò)展有限元模擬時可取右半部為計算對象，計算中用到幾何尺寸、材料參數(shù)和外載荷列于表1中，其中E和μ為彈性模量和泊松比.與有限元法相比較，用擴(kuò)展有限元法計算裂紋問題的結(jié)果與網(wǎng)格的關(guān)系不大，不必在裂紋尖端處加密網(wǎng)格.為保證數(shù)值精度，在數(shù)據(jù)擬合時，極半徑r的取值范圍不能過大，通常取單元邊長的5～15倍即可.

圖6 有中心裂紋的有限寬板

表1 計算模型的幾何尺寸、材料參數(shù)和外載荷

圖7為擴(kuò)展有限元法的計算模型，其中長度單位為m，中間的粗實線代表裂紋所在位置，水平方向和豎直方向的單元數(shù)量分別為50和101，單元類型為平面應(yīng)力四節(jié)點(diǎn)等參元.邊界條件為：x=0的左側(cè)邊的水平位移u=0；在裂紋下側(cè)離裂紋垂直距離最短的、右側(cè)4個節(jié)點(diǎn)的豎直位移v=0.外力邊界條件為：在y=1和y=-1的上下兩邊，施加σ=10MPa拉伸靜載荷.

圖8為根據(jù)節(jié)點(diǎn)位移分量，繪制的擴(kuò)展有限元法計算模型的變形圖，為便于觀察與分析計算模型的變形規(guī)律，圖中的節(jié)點(diǎn)位移分量被放大了2000倍.不難發(fā)現(xiàn)計算模型的變形圖所反映的變形規(guī)律，符合計算模型中的位移邊界條件和外載荷條件，因此根據(jù)擴(kuò)展有限元法計算出的節(jié)點(diǎn)位移符合實際情況.

圖9為根據(jù)節(jié)點(diǎn)位移幅值，繪制的計算模型位移幅值云圖，其中的位移幅值單位為m.可以看出計算模型的位移幅值，關(guān)于y=0描述的水平線具有對稱性，這符合計算模型的位移邊界條件和外力邊界條件所決定的變形規(guī)律.

圖7 擴(kuò)展有限元法的計算模型

圖8 計算模型的變形圖

圖10為根據(jù)節(jié)點(diǎn)豎向位移分量，繪制的計算模型豎向位移云圖，其中豎向位移單位為m.可以看出計算模型的豎向位移，關(guān)于y=0描述的水平線具有反對稱性，這符合計算模型的位移邊界條件和外力邊界條件所決定的變形規(guī)律.

圖9 計算模型的位移幅值云圖

圖10 計算模型的豎向位移云圖

圖 11為利用擴(kuò)展有限元法計算出裂紋面附近節(jié)點(diǎn)豎向位移，通過MATLAB編寫位移外推法計算程序，得到位移外推法數(shù)據(jù)擬合曲線，其中的星號標(biāo)出的點(diǎn)為利用式(24)構(gòu)造的數(shù)據(jù)點(diǎn)虛線為利用式(26)和式(29)擬合的曲線，曲線縱坐標(biāo)截距即為利用式(29b)計算出的I型裂紋的應(yīng)力強(qiáng)度因子，具體結(jié)果列于表2中.

圖11 位移外推法的擬合曲線

表2 不同方法的計算結(jié)果

結(jié)合前面擴(kuò)展有限元法的結(jié)算結(jié)果，利用MATLAB編寫相互作用積分法的計算程序，計算I型裂紋的應(yīng)力強(qiáng)度因子，結(jié)果也列于表2中.為便于比較兩種數(shù)值算法的精度，利用文獻(xiàn)[19]中的解析公式

計算I型裂紋的應(yīng)力強(qiáng)度因子，并計算兩種數(shù)值方法與解析法的相對誤差，結(jié)果均列于表2中.比較兩種數(shù)值算法與解析法計算結(jié)果的相對誤差，可以發(fā)現(xiàn)位移外推法的計算精度高于相互作用積分法的計算精度.而位移外推法的求解過程和程序設(shè)計，比相互作用積分法的求解過程和程序設(shè)計簡單很多.因此本文建立的基于擴(kuò)展有限元的位移外推法，能夠更高效和準(zhǔn)確地計算平面裂紋的應(yīng)力強(qiáng)度因子.

4 結(jié)論

介紹了擴(kuò)展有限元法基本原理，闡述了平面裂紋問題的四節(jié)點(diǎn)等參元單元的位移模式，推導(dǎo)了擴(kuò)展有限元法的控制方程，介紹了特殊單元的數(shù)值積分方法.基于最小二乘法和擴(kuò)展有限元法的計算結(jié)果，建立了應(yīng)力強(qiáng)度因子位移外推法的計算公式.利用MATLAB編寫了相關(guān)計算程序，對平面裂紋問題用擴(kuò)展有限元法進(jìn)行了計算.結(jié)合擴(kuò)展有限元法的計算結(jié)果，分別利用位移外推法和相互作用積分法，對平面裂紋的應(yīng)力強(qiáng)度因子進(jìn)行了計算.數(shù)值計算與解析結(jié)果的對比表明,位移外推法比相互作用積分法能更方便和準(zhǔn)確地計算平面裂紋的應(yīng)力強(qiáng)度因子.

1范天佑.斷裂理論基礎(chǔ).北京:科學(xué)出版社,2003:73-112

2 Belytschko T,Black T.Elastic crack growth in f i nite elements with minimal remeshing.International Journal for Numerical Method in Engineering,1999,45(5):601-620

3 Moes N,Dolbow J,Belytschko T.A f i nite element method for crack growth without without remeshing.International Journal for Numerical Method in Engineering,1999,46(1): 131-150

4 Mose N,Belytschko T.Extended f i nite element method for cohesive crack growth.Engineering Fracture Mechanics, 2002,69(7):813-833

5余天堂.含裂紋體的數(shù)值模擬.巖石力學(xué)與工程學(xué)報,2005,24: 4432-4439

6方修君,金峰,王進(jìn)廷.基于擴(kuò)展有限元法的粘聚裂紋模型.清華大學(xué)學(xué)報(自然科學(xué)版),2007,47(3):344-347

7方修君,金峰.基于ABAQUS平臺的擴(kuò)展有限元法.工程力學(xué),2007,24(7):6-10

8懂玉文,余天堂,任青文.直接計算應(yīng)力強(qiáng)度因子的擴(kuò)展有限元法.計算力學(xué)學(xué)報,2008,25(1):72-77

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

10余天堂.摩擦接觸裂紋問題的擴(kuò)展有限元法.工程力學(xué),2010, 27(4):84-89

11茹忠亮,朱傳銳,張友良等.裂紋問題的擴(kuò)展有限元法研究.巖土力學(xué),2011,32(7):2171-2176

12杜修力,金瀏,黃景琦.基于擴(kuò)展有限元法的混凝土細(xì)觀斷裂破壞過程模擬.計算力學(xué)學(xué)報,2012,29(6):940-947

13茹忠亮，申葳，趙洪波.基于擴(kuò)展有限元法的鋼筋混凝土梁復(fù)合斷裂過程模擬研究.工程力學(xué)，2013,30(5):215-220

14傅玉珍,張華,談?wù)?基于擴(kuò)展有限元的巴西圓盤動態(tài)裂縫擴(kuò)展分析.四川建筑科學(xué)研究,2014,40(2):47-54

15師訪,高峰,楊玉貴.正交各向異性巖體裂紋擴(kuò)展的擴(kuò)展有限元法研究.巖土力學(xué),2014,35(4):1203-1210

16盛茂，李根生.水力壓裂過程的擴(kuò)展有限元數(shù)值模擬方法.工程力學(xué)，2014,31(10):123-128

17江守燕，杜成斌，顧沖時等.求解雙材料界面裂紋應(yīng)力強(qiáng)度因子的擴(kuò)展有限元法.工程力學(xué)，2015,32(3):22-27

18莊茁，柳占立，成斌斌等.擴(kuò)展有限單元法.北京：清華大學(xué)出版社，2012

19薛世峰，侯密山.工程斷裂力學(xué).東營：中國石油大學(xué)出版社，2011

(責(zé)任編輯:劉希國)

DISPLACEMENT EXTRAPOLATION METHOD FOR STRESS INTENSITY FACTOR BASED ON EXTENDED FINITE ELEMENT METHOD1)

ZHOU Bo2)XUE Shifeng
(College of Pipeline and Civil Engineering,China University of Petroleum,Qingdao 266580,Shandong,China)

The element displacement modes of the extended f i nite element method for plane crack problems are discussed in this paper.The governing equations for the extended f i nite element method are derived based on the principle of virtual work.The numerical integration technique is adopted for special elements.The displacement extrapolation method is formulated to calculate the stress intensity factor based on the least square method.The numerical computation programs are developed by using the MATLAB,and the numerical calculations are carried out by the extended f i nite element method for the problem of plane crack.Based on the results of the extended f i nite element method,the stress intensity factors of plane crack are calculated by the displacement extrapolation method and the interaction integral method,respectively.Numerical results show that with the displacement extrapolation method,the stress intensity factors can be obtained more easily and precisely than the interaction integral method.

extended f i nite element method,stress intensity factor,displacement extrapolation method,interaction integral method

O346,TB115

：Adoi：10.6052/1000-0879-16-413

2016–12–20收到第1稿，2017–02–05收到修改稿.

1)國家自然科學(xué)基金(11472309)和中石油重大專項(2012-ZG-002)資助項目.

2)周博，教授.E-mail:zhoubo@upc.edu.cn

周博,薛世峰.基于擴(kuò)展有限元的應(yīng)力強(qiáng)度因子的位移外推法.力學(xué)與實踐,2017,39(4):371-378

Zhou Bo,Xue Shifeng.Displacement extrapolation method for stress intensity factor based on extended f i nite element method.Mechanics in Engineering,2017,39(4):371-378