比例邊界有限元法求解裂紋面接觸問(wèn)題1)

章 鵬 杜成斌 江守燕

(河海大學(xué)工程力學(xué)系，南京211100)

比例邊界有限元法求解裂紋面接觸問(wèn)題1)

章 鵬 杜成斌2)江守燕

(河海大學(xué)工程力學(xué)系，南京211100)

比例邊界有限元側(cè)面上有任意荷載時(shí)，將側(cè)面載荷分解成關(guān)于徑向方向局部坐標(biāo)的多項(xiàng)式函數(shù)的和，推導(dǎo)給出了考慮側(cè)面載荷存在的新型形函數(shù)，并基于該形函數(shù)推導(dǎo)了剛度矩陣和等效節(jié)點(diǎn)載荷列陣.首次對(duì)比例邊界有限元法求解裂紋面接觸問(wèn)題進(jìn)行了研究，運(yùn)用Lagrange乘子引入接觸界面約束條件，推導(dǎo)給出了比例邊界有限元求解裂紋面接觸問(wèn)題的控制方程.將裂紋面單元分為非裂尖單元和含有側(cè)面的裂尖單元.在非裂尖單元中的裂紋面，裂紋面作為多邊形單元的邊界，邊界上的接觸力可等效到節(jié)點(diǎn)上，通過(guò)在節(jié)點(diǎn)上構(gòu)造Lagrange乘子，采用點(diǎn)對(duì)點(diǎn)接觸約束進(jìn)行處理.對(duì)于含有側(cè)面的裂尖單元，在整個(gè)側(cè)面上構(gòu)造Lagrange乘子的插值場(chǎng)，采用邊對(duì)邊接觸約束進(jìn)行處理.對(duì)三個(gè)不同的接觸約束狀態(tài)下的算例進(jìn)行了數(shù)值計(jì)算，通過(guò)與解析解及有限元軟件ABAQUS計(jì)算結(jié)果的對(duì)比，驗(yàn)證了本文提出的比例邊界有限元點(diǎn)對(duì)點(diǎn)和邊對(duì)邊接觸求解裂紋面接觸問(wèn)題的精確性與有效性.

比例邊界有限元,形函數(shù),裂紋,點(diǎn)對(duì)點(diǎn)接觸,邊對(duì)邊接觸

引言

在工程實(shí)際中，由于各種原因，工程結(jié)構(gòu)或多或少都存在著裂縫(紋)[1].在外力作用下，既存在受拉狀態(tài)的張開型裂紋，也有壓剪狀態(tài)的閉合型裂紋.用數(shù)值方法分析閉合型裂紋時(shí)，需要在裂紋面上施加合適的接觸條件，以反映裂紋面間的接觸狀態(tài)，否則裂紋面會(huì)發(fā)生與實(shí)際不符的相互嵌入現(xiàn)象[2-3].目前施加接觸條件的方法主要有Lagrange乘子法[4-5]，增廣的Lagrange乘子法[6]，罰函數(shù)法[7-8]，線性互補(bǔ)法[9]等.

比例邊界有限元法(scaled boundary fi nite element method，SBFEM)是由Wolf和Song在20世紀(jì)90年代末率先提出和發(fā)展起來(lái)的一種半解析的數(shù)值計(jì)算方法[10-11].該方法在計(jì)算過(guò)程中僅需離散結(jié)構(gòu)的邊界，降低了數(shù)值模擬維度，而且可以半解析地表征裂紋尖端的奇異性[12-14].基于上述優(yōu)點(diǎn)，SBFEM在工程問(wèn)題方面的應(yīng)用正逐漸成為學(xué)術(shù)研究的前沿和熱點(diǎn)[15-18].比例邊界有限單元可以為含有任意邊數(shù)的多邊形，該特點(diǎn)也使在求解裂紋擴(kuò)展時(shí)重劃分網(wǎng)格具有簡(jiǎn)單性和高效性.Ooi等[19]推導(dǎo)給出了不考慮側(cè)面載荷存在的多邊形SBFEM位移形函數(shù)，并利用該形函數(shù)，推導(dǎo)了采用SBFEM求解彈塑性問(wèn)題的有關(guān)公式，拓寬了SBFEM的應(yīng)用范圍，使得SBFEM可以求解一些非線性問(wèn)題.

接觸問(wèn)題屬于典型的非線性問(wèn)題，在有限元法[20-21]、擴(kuò)展有限元法[22-23]、邊界元法[24]中已經(jīng)得到了較長(zhǎng)時(shí)間的發(fā)展.但是在各種方法中或多或少存在一些問(wèn)題，有限元求解裂紋問(wèn)題時(shí)，因?yàn)榱鸭y的不連續(xù)性，需要將裂紋定義在單元邊界上，且在裂尖區(qū)域需要?jiǎng)澐址浅＜?xì)的網(wǎng)格來(lái)反映應(yīng)力奇異性，這將很大程度上降低求解裂紋接觸問(wèn)題的效率.擴(kuò)展有限元中裂紋位于單元內(nèi)部，采用點(diǎn)對(duì)點(diǎn)接觸方案不能正確的表示裂紋面接觸，只能采用砂漿法(mortar method)或邊對(duì)邊(segment-to-segment)接觸方法來(lái)求解裂紋面接觸應(yīng)力，對(duì)每一個(gè)裂紋面單元都需數(shù)值積分求解，且對(duì)于裂尖單元的接觸應(yīng)力只可假設(shè)為恒定值[22]，這將很大程度上降低接觸應(yīng)力的求解精度.邊界元法求解接觸較有優(yōu)勢(shì)，但是其本身方法需要求得嚴(yán)格滿足控制方程的基本解，不適于對(duì)非均質(zhì)、非各向同性介質(zhì)的求解，還需進(jìn)一步發(fā)展研究.比例邊界有限元作為一種新型的高精度的數(shù)值方法，尚未發(fā)現(xiàn)有文獻(xiàn)采用該方法求解裂紋面接觸問(wèn)題.采用SBFEM求解裂紋面接觸問(wèn)題時(shí)會(huì)遇到縫面存在載荷的情況，如縫面接觸應(yīng)力等.所以首先必須給出側(cè)面存在載荷的比例邊界有限元的有關(guān)列式.本文將任意側(cè)邊載荷表示為關(guān)于徑向坐標(biāo)ξ的多項(xiàng)式函數(shù)[25-26]，推導(dǎo)了含有側(cè)面載荷影響的SBFEM的形函數(shù)，并給出了相應(yīng)的剛度矩陣以及等效載荷列式.將裂紋面單元分為非裂尖單元和含有側(cè)面的裂尖單元.對(duì)于不含有裂尖的單元，裂紋面位于單元邊界上，在節(jié)點(diǎn)上構(gòu)造Lagrange乘子，采用點(diǎn)對(duì)點(diǎn)約束模擬裂紋面的接觸.對(duì)于含有裂尖的單元，其裂紋面位于裂尖單元的側(cè)面(side faces)上，由于側(cè)面為非邊界面，按常規(guī)思路施加點(diǎn)對(duì)約束不能正確的表示裂紋面接觸，因此本文提出了在裂尖單元的裂紋面上采用邊對(duì)邊約束方法處理接觸問(wèn)題，即在整個(gè)側(cè)面上構(gòu)造Lagrange乘子的插值場(chǎng)，Lagrange乘子呈線性變化，通過(guò)Lagrange乘子推導(dǎo)出了SBFEM接觸控制方程.通過(guò)若干不同的算例，驗(yàn)證了采用多邊形SBFEM點(diǎn)對(duì)點(diǎn)--邊對(duì)邊接觸可精確的模擬含有裂尖單元的接觸問(wèn)題.

1 含有側(cè)面載荷的比例邊界多邊形有限元支配方程

1.1 考慮側(cè)面載荷的比例邊界多邊形有限元的形函數(shù)

圖1為一個(gè)含有裂紋的任意多邊形，采用4個(gè)任意多邊形SBFEM單元[27]進(jìn)行離散，其中子域1包含裂紋信息(裂紋單元).根據(jù)比例邊界有限元的單元形態(tài)要求，在每一個(gè)單元子域內(nèi)，均需要選取一個(gè)比例中心，通過(guò)此中心該子域的全部邊界都可見(jiàn)[28].處理復(fù)雜結(jié)構(gòu)問(wèn)題時(shí)，可通過(guò)合理布置比例中心位置將模型離散成多個(gè)多邊形子域來(lái)滿足這一基本要求.

圖1 比例邊界有限元多邊形離散模型Fig.1 Scaled boundary fi nite element method polygon discrete model

圖2為具有裂紋面任意載荷的比例邊界有限元單元模型，O(裂尖)為比例中心，斷裂問(wèn)題比例中心通常選在裂尖處，定義ξ(0 6ξ6 1)為徑向坐標(biāo)，s(0 6s6 1)為環(huán)向坐標(biāo)，模型邊界(ξ=1)離散成一維線單元，而裂紋側(cè)面(side faces)不需要離散.ξ?s形成比例邊界有限元局部坐標(biāo)系.其中環(huán)向坐標(biāo)s正向沿著單元邊界逆時(shí)針變化[11]，徑向坐標(biāo)ξ正向沿著比例中心向外變化，其中ξ=0代表比例中心O，ξ=1代表邊界上的點(diǎn).

圖2 比例邊界有限元裂紋單元模型(?為比例中心，?為頂點(diǎn)，·為節(jié)點(diǎn))Fig.2 Scaled boundary fi nite element method crack element model(? scaling centres，? vertices，·nodes)

對(duì)于裂紋單元內(nèi)任意一點(diǎn)，相對(duì)于比例中心的局部坐標(biāo)(x,y)T可用ξ和s表示為

其中，(x0y0)為邊界上的節(jié)點(diǎn)坐標(biāo).N(s)表示為沿s方向的形函數(shù)

其中n為邊界節(jié)點(diǎn)個(gè)數(shù).

考慮側(cè)面載荷的比例邊界有限元的控制方程可通過(guò)加權(quán)平均法[7]或者虛功原理[8]得到

該方程為徑向上的平衡方程.其中E0,E1,E2為單元的系數(shù)矩陣，F(xiàn)t(ξ)為側(cè)面載荷，u(ξ)為節(jié)點(diǎn)位移函數(shù).

式(3)為含有2n個(gè)未知數(shù)u(ξ)的二階微分方程，可化為一階微分方程

其中q(ξ)為內(nèi)部節(jié)點(diǎn)力，Z為2n×2n階Hamiltonian矩陣，滿足

采用Schur分解[29]

得到單元的模態(tài)位移ψd，模態(tài)載荷ψq及對(duì)應(yīng)的特征值矩陣Λλ.

裂紋側(cè)面上的任意載荷可以分解為關(guān)于ξ的M階多項(xiàng)式函數(shù)[25]

式中Wi為單元側(cè)邊力分布矩陣，ai為第i階多項(xiàng)式系數(shù).

側(cè)面載荷的第i階位移模態(tài)[11]為

則徑向位移解為

式(10)中第1項(xiàng)為式(3)的齊次通解.

式(10)的矩陣形式為

其中Λλ為特征值矩陣，Λt=dig(1,2,3···,M),a為ai組成的系數(shù)向量，c為積分常數(shù)，由邊界條件確定.

邊界上的節(jié)點(diǎn)位移為

其中積分常數(shù)c用邊界位移表達(dá)為

將式(13)代入式(11)可得

式(14)的矩陣形式可寫為

將式(15b)和式(15c)代入式(15a)可得

比例邊界有限元在環(huán)向s上采用與有限單元法中形函數(shù)類似方法，通過(guò)N(s)進(jìn)行插值，即

將式(16)代入式(17)可得

令

可得

由式(20)可知，Nt(ξ,s)為采用比例邊界有限元法考慮側(cè)面任意載荷時(shí)的形函數(shù).將式(19)各變量代入式(15b)和式(15c)得到

對(duì)于不考慮側(cè)面載荷的比例邊界元形函數(shù)為[16]

由式 (21)和式 (22)對(duì)比可看出Nt(ξ,s)左半部分N(s)ψdξ?Λλψ?1d為不考慮側(cè)面載荷的形函數(shù)，右半部分為多項(xiàng)式函數(shù)

該多項(xiàng)式函數(shù)可看作是由側(cè)面載荷對(duì)位移模式的影響，式(16)u0中的a作為特解位移的值，可看作為額外自由度.式(22)和式(23)兩部分組成了考慮側(cè)邊力任意載荷時(shí)的形函數(shù).從中可看出，無(wú)側(cè)面載荷的形函數(shù)為考慮側(cè)面載荷的形函數(shù)的一個(gè)特例，即a=0.

SBFEM的應(yīng)變矩陣[11]為

其中B1(s)和B2(s)是應(yīng)變位移矩陣.

將式(16)代入式(24)得

引入應(yīng)變位移矩陣

應(yīng)變矩陣可寫為

應(yīng)力矩陣為

其中D是材料的彈性矩陣.

1.2 比例邊界多邊形有限元支配方程

考慮面載荷和側(cè)邊載荷時(shí)，采用虛功原理的表達(dá)式為

式中t(s)為邊界上的面力[11]，其他符號(hào)同前.

將式(20)、式(28)、式(29)代入式(30)得

式(31)等價(jià)為

式(31)中左邊括號(hào)部分可視為剛度矩陣，即

將式(27)代入式(33)可得

將dV進(jìn)行積分得到

定義矩陣Y

式(36)可以采用高斯積分或者Guass-Lobatto-Legendre積分[30]方法進(jìn)行積分求解.

考慮式(36)，式(35)可改寫為

定義矩陣X

采用分部積分，式(38)可化簡(jiǎn)為L(zhǎng)yapunov方程[31]

式(39)可通過(guò)Matlab自帶函數(shù)求解Lyapunov方程，得出矩陣X.

將式(38)代入式(37)得出剛度矩陣K

式(32)中等式右側(cè)為等效節(jié)點(diǎn)載荷Fp

將式(33)、式(41)代入式(31)中可得到考慮側(cè)面載荷的比例邊界有限元支配方程

2 接觸問(wèn)題的有關(guān)公式

兩物體接觸受壓時(shí)，在物體表面需引入物體表面的接觸條件，以避免物體發(fā)生相互侵入[32].如圖3所示，考慮二維問(wèn)題，可能發(fā)生接觸的兩個(gè)表面記為SA和SB，為了系統(tǒng)地分析兩物體表面的接觸條件的施加方法，建立局部坐標(biāo)系，設(shè)xA為SA上任一指定點(diǎn)P的坐標(biāo)，則該點(diǎn)至SB面上最接近點(diǎn)Q的法向相對(duì)距離gN為

圖3 接觸點(diǎn)對(duì)與點(diǎn)對(duì)間的距離Fig.3 The contact points and distance between the points

同理定義切向相對(duì)距離gT

對(duì)于彈性接觸問(wèn)題，滿足Coulomb定律的接觸力和物體相對(duì)距離的關(guān)系可用下列等式和不等式[33]表示

接觸問(wèn)題可描述為求區(qū)域內(nèi)位移場(chǎng)，使得系統(tǒng)的勢(shì)能達(dá)到最小.Lagrange乘子法[4]是求解接觸約束最小化問(wèn)題的常用方法之一，通過(guò)引入Lagrange乘子將接觸問(wèn)題轉(zhuǎn)化為無(wú)約束問(wèn)題求解.接觸狀態(tài)有3種：張開、粘結(jié)、滑移.法向Lagrange乘子場(chǎng)λN表示當(dāng)接觸發(fā)生時(shí)迫使gN等于0(裂紋面接觸但無(wú)嵌入)的壓力.在裂紋面為粘結(jié)接觸時(shí)，切向Lagrange乘子場(chǎng)λT表示在粘結(jié)區(qū)迫使?jié)M足粘結(jié)條件的切向力，在裂紋面發(fā)生滑移接觸時(shí)，切向Lagrange乘子場(chǎng)λT表示為滑移的摩擦力.

應(yīng)用Lagrange乘子法在物體表面考慮摩擦接觸條件時(shí)可導(dǎo)出3個(gè)方程：

(1)接觸不嵌入條件

(2)在切向滿足粘結(jié)條件時(shí)，即

在切向滿足有摩擦的滑動(dòng)接觸狀態(tài)時(shí)：切向運(yùn)動(dòng)不再受約束，但是切向力滿足

(3)虛功原理(外力虛功包括接觸應(yīng)力)，即

3 基于比例邊界有限元求解摩擦接觸問(wèn)題

3.1 接觸約束公式

如圖4所示為比例邊界有限元網(wǎng)格，假設(shè)滿足粘結(jié)接觸，則法向在粘結(jié)約束條件下，裂紋面上的P點(diǎn)和Q點(diǎn)法向接觸間距為0.根據(jù)式(48)可得

其中外法向向量nQ={?sinαcosα}T.

圖4 點(diǎn)對(duì)點(diǎn)與邊對(duì)邊約束示意圖Fig.4 Point-to-segment constraints diagram

切向方向上，在粘結(jié)接觸條件下，根據(jù)式(49)可得

正切向向量sQ={cosαsinα}T.

3.1.1 不含裂尖的裂紋面單元

在不含裂尖的裂紋面單元中，如圖4的第1,2個(gè)單元，裂紋面接觸分別在兩個(gè)單元的邊界上(ξ=1)，在該邊界上其形函數(shù)N(s,ξ=1)=N(s)只是關(guān)于s的函數(shù)，可以采用虛功原理將邊界上的接觸應(yīng)力等效到該裂紋面邊界節(jié)點(diǎn)上，邊界節(jié)點(diǎn)上的接觸載荷為等效節(jié)點(diǎn)載荷.所以只需點(diǎn)約束即可精確模擬接觸，在節(jié)點(diǎn)上構(gòu)造Lagrange乘子(λN,λT)，則在第k個(gè)點(diǎn)對(duì)P點(diǎn)和Q點(diǎn)間，式(52)可寫為

其中，CNk為矩陣CN第k行，Cindex為點(diǎn)對(duì)P和Q兩節(jié)點(diǎn)在整體節(jié)點(diǎn)中的坐標(biāo)變換.

由式(54)可得裂紋面點(diǎn)對(duì)點(diǎn)接觸法向方向約束方程

同樣在切向方向，令

可得裂紋面點(diǎn)對(duì)接觸切向方向約束方程

3.1.2 包含裂尖的裂紋面單元

將含有裂尖單元的比例中心放在裂尖處，如圖4的單元3所示.在該單元的裂紋面?zhèn)让?side faces)上，ξ從0變化到1，該單元所受的裂紋面接觸應(yīng)力為側(cè)面應(yīng)力.該側(cè)面的形函數(shù)Nt(ξ,s)為關(guān)于(ξ,s)的函數(shù)，因?yàn)榇藭r(shí)裂紋面為非邊界面(ξ(0 6ξ6 1))，按常規(guī)思路施加點(diǎn)對(duì)約束不能正確地表示裂紋面接觸.為此本文提出了在裂尖單元中采用邊對(duì)邊接觸方法，即在整個(gè)側(cè)面上構(gòu)造Lagrange乘子的插值場(chǎng)，Lagrange乘子呈線性變化，其插值可表示為

其中,λj1為側(cè)邊中心(裂尖中心)的接觸應(yīng)力，λj2為側(cè)邊端點(diǎn)的接觸應(yīng)力.

Mi為L(zhǎng)agrange乘子插值函數(shù)

則式(52)可寫為

由式(19)形函數(shù)可知

同理可求出uPy，uQx，uQy代入式(61)并進(jìn)行運(yùn)算積分可得

其中CNc為考慮式(61)和式(62)運(yùn)算積分所得.

最終得到裂尖單元邊對(duì)邊接觸法向約束方程

同理可得到裂紋側(cè)面邊對(duì)邊接觸切向方向約束方程

考慮插值得到的Lagrange乘子，式(51)的虛功方程用矩陣形式可表示為

其中，K和F為SBFEM的剛度矩陣和載荷列陣，u0為SBFEM位移列陣；λN和λT為法向和切向Lagrange乘子.

根據(jù)虛功方程(66)和前文所得的約束方程，假設(shè)裂紋面滿足粘結(jié)接觸條件，則系統(tǒng)的控制方程為

式(67)第1行為虛功方程推導(dǎo)出的平衡方程，第2行為法向接觸條件，第3行為切向粘結(jié)接觸條件.在滿足有摩擦的滑移條件的狀態(tài)下，其切向運(yùn)動(dòng)不受約束，即式(67)第3行將不再成立.而切向摩擦力滿足式(50).

故滿足滑動(dòng)接觸狀態(tài)的控制方程為

4 數(shù)值算例

4.1 貫穿裂紋板接觸

圖5為含貫穿裂紋的矩形板，板左端固定，右端受到q=100Pa的均勻壓應(yīng)力.板長(zhǎng)3m，寬為2m，裂紋面與x軸夾角α=50?，裂紋與上邊界的交點(diǎn)坐標(biāo)為(1.9m,2m).楊氏模量E=76kPa，泊松比ν=0.3，采用平面應(yīng)變假設(shè).假設(shè)斜裂紋摩擦系數(shù)f足夠大，裂紋面處于完全粘結(jié)接觸狀態(tài).采用多邊形單元[34-35]進(jìn)行離散.由于結(jié)構(gòu)比較簡(jiǎn)單，因此將模型劃分成2個(gè)子域，并采用2節(jié)點(diǎn)單元離散子域的邊界(節(jié)點(diǎn)用?表示)，如圖6所示.由于是完全貫穿的裂紋，只需在裂紋面上采用點(diǎn)對(duì)點(diǎn)約束即可模擬接觸問(wèn)題.

圖5 含貫穿裂紋的矩形板(m)Fig.5 A rectangular plate divided by through crack(m)

圖6 多邊形比例邊界有限元網(wǎng)格Fig.6 Polygon SBFEM mesh

圖7為板受壓變形后的位移云圖.由于裂紋面接觸條件的施加，從位移云圖中可看出，位移具有連續(xù)性，在受壓裂紋面上沒(méi)有發(fā)生相互嵌入，從而表明接觸條件施加是正確的.

為了驗(yàn)證該方法精確性與收斂性，采用結(jié)構(gòu)應(yīng)變能相對(duì)誤差作為指標(biāo).應(yīng)變能誤差的表達(dá)式[36]可以寫為應(yīng)變能的相對(duì)誤差為

圖7 位移云圖Fig.7 Displacement fi eld

其中，σh是本文方法求得的應(yīng)力，σexact為應(yīng)力值解析解，Uexact為應(yīng)變能解析解.

為了比較SBFEM求解裂紋面接觸問(wèn)題的精度和效率，采用目前應(yīng)用較為廣泛的商業(yè)有限元軟件ABAQUS進(jìn)行模擬對(duì)比.ABAQUS在接觸程序處理中，采用Lagrange乘子接觸[37].圖8分別給出了兩種方法計(jì)算得到的應(yīng)變能相對(duì)誤差隨不同網(wǎng)格(粗網(wǎng)格、中等網(wǎng)格和細(xì)網(wǎng)格)的自由度變化，從圖中可看出，SBFEM方法求解裂紋面接觸問(wèn)題收斂極快，僅需16個(gè)自由度就無(wú)限趨近于解析解(應(yīng)變能相對(duì)誤差為(3×10?20)%)，隨著網(wǎng)格的變化，求解精度基本不變.而ABAQUS收斂性較慢(粗網(wǎng)格下應(yīng)變能相對(duì)誤差(3.6×10?3)%)，而且在網(wǎng)格較細(xì)，自由度較高的情況下，其計(jì)算精度仍遠(yuǎn)低于SBFEM較低的自由度計(jì)算精度，這說(shuō)明了SBFEM求解裂紋面接觸問(wèn)題的高精度與高效性.

圖8 應(yīng)變能相對(duì)誤差隨網(wǎng)格自由度變化Fig.8 The relative error in the strain energy varies with the grid degrees of freedom

4.2 含有圓孔的縫內(nèi)為非均勻壓力的裂紋板接觸

為了進(jìn)一步驗(yàn)證該方法在一般情況下的求解能力，本文對(duì)含有孔洞及裂紋的無(wú)限大板在壓力狀態(tài)下進(jìn)行了數(shù)值模擬.式(71)～式(73)給出了含有半圓孔的半無(wú)限大板承受壓力載荷時(shí)的解析解[38].該板在極坐標(biāo)下的解為

當(dāng)板的尺寸長(zhǎng)度L與圓孔半徑R的比例足夠大時(shí)，可以近似用來(lái)模擬含有半孔洞的無(wú)限大板.圖9為該模型的幾何尺寸.坐標(biāo)系原點(diǎn)設(shè)在圓孔的圓心處，板的長(zhǎng)度為30m，板的高度為15m，圓孔的半徑R為1m，在圓孔的正下方有一條長(zhǎng)度為a=1.68m的裂紋，板的兩端受到σ∞=100Pa的均勻壓力，板的上端數(shù)值方向位移為0.假設(shè)摩擦因子f足夠大，裂紋處于完全粘結(jié)接觸狀態(tài).

圖9 包含有裂縫及圓孔的矩形板(m)Fig.9 Model of a rectangular plate with a crack and a circular hole(m)

板在孔洞周圍存在應(yīng)力集中現(xiàn)象，因此該模型采用較密的多邊形網(wǎng)格進(jìn)行離散，如圖10所示(圖中*為裂紋尖端).采用點(diǎn)對(duì)點(diǎn)與邊對(duì)邊接觸方法求解該裂紋接觸問(wèn)題.

圖11和圖 12給出了結(jié)構(gòu)的變形圖和位移云圖，由于原結(jié)構(gòu)變形極小，圖11對(duì)原結(jié)構(gòu)變形進(jìn)行了一定比例的放大(放大了70倍).從結(jié)構(gòu)的變形圖和位移云圖可看出，采用該方法考慮接觸時(shí)，變形與沒(méi)有裂紋情況幾乎一樣，受壓裂紋面上位移具有連續(xù)性，沒(méi)有發(fā)生相互嵌入.從而表明接觸算法是正確的.

圖10 包含有裂縫及圓孔的矩形板多邊形網(wǎng)格Fig.10 Polygon SBFEM mesh of a rectangular plate with a crack and a circular hole

圖11 結(jié)構(gòu)變形圖(放大70倍，虛線為未變形圖)Fig.11 Structural deformation diagram(magni fi ed 70 times，dashed line is the original shape)

圖12 位移云圖Fig.12 Displacement fi eld

圖13為采用SBFEM和ABAQUS兩種方法計(jì)算的應(yīng)變能相對(duì)誤差隨不同網(wǎng)格自由度的變化.由圖中可看出，基于SBFEM求解的結(jié)果在自由度數(shù)較小情況下，誤差較小(相對(duì)誤差為0.0032%)，而基于ABAQUS的求解誤差在自由度數(shù)較小情況下所得誤差較大(相對(duì)誤差為1.7%)，隨著網(wǎng)格的加密，自由度數(shù)的增大，誤差逐漸減小，但ABAQUS求解精度遠(yuǎn)小于SBFEM求解精度.這說(shuō)明SBFEM點(diǎn)對(duì)點(diǎn)和邊對(duì)邊約束方法在非均布應(yīng)力下的裂紋接觸仍具有較高的精度和效率.

圖13 應(yīng)變能相對(duì)誤差隨網(wǎng)格自由度變化Fig.13 The relative error in the strain energy varies with the grid degrees of freedom

4.3 考慮摩擦的滑移接觸

在滿足粘結(jié)接觸下，裂紋接觸面在法向和切向方向都無(wú)相對(duì)位移.在有摩擦的滑移接觸下，接觸在法向方向仍然滿足不嵌入條件，在切向方向，發(fā)生相互移動(dòng)，系統(tǒng)方程滿足式(68).為了驗(yàn)證SBFEM滑移接觸的有效性，本文對(duì)一個(gè)含有裂紋尖端的矩形板進(jìn)行了模擬，圖14為該矩形板的多邊形網(wǎng)格，其中裂尖坐標(biāo)為(0.809m,0.7m)，摩擦系數(shù)f=0.2，裂紋在切向上發(fā)生了滑移.其他力學(xué)參數(shù)與算例4.1中參數(shù)相同.

圖14 多邊形比例邊界有限元網(wǎng)格Fig.14 Polygon SBFEM mesh

在裂紋接觸面上產(chǎn)生滑移，兩個(gè)接觸面產(chǎn)生不連續(xù)位移，故圖14采用了相對(duì)較多的多邊形SBFEM網(wǎng)格對(duì)矩形板進(jìn)行離散.圖15和圖16分別為矩形板變形圖和位移云圖，由圖15中可看到，切向上產(chǎn)生了滑移，在法向上裂紋面間沒(méi)有相互侵入.從而表明該滑移接觸算法是正確的.

圖15 結(jié)構(gòu)變形圖(放大30倍，虛線為未變形圖)Fig.15 Structural deformation diagram(magni fi ed 30 times，dashed line is the original shape)

圖16 位移云圖Fig.16 Displacement fi eld

本算例沒(méi)有解析解，仍采用商業(yè)有限元軟件ABAQUS模擬結(jié)果來(lái)對(duì)比本文提出的SBFEM求解裂紋面接觸方法的精度和效率.在使用ABAQUS模擬過(guò)程中，分別采用粗網(wǎng)格、中等網(wǎng)格和細(xì)網(wǎng)格求解該滑移接觸問(wèn)題，圖17分別給出了ABAQUS不同網(wǎng)格的模型.

圖17 ABAQUS不同網(wǎng)格模型Fig.17 Di ff erent mesh models of ABAQUS

圖18為采用不同模型模擬得到的裂紋面下表面σxx分布，可看出在ABAQUS粗網(wǎng)格下，其應(yīng)力分布與細(xì)網(wǎng)格相差較大,特別是靠近裂紋接觸的兩端，其應(yīng)力誤差較為明顯，同時(shí)可看出SBFEM模擬結(jié)果與細(xì)網(wǎng)格的模擬結(jié)果較為一致.從網(wǎng)格對(duì)比來(lái)看，采用SBFEM求解的網(wǎng)格自由度數(shù)(數(shù)量為1560)遠(yuǎn)小于ABAQUS細(xì)網(wǎng)格的自由度數(shù)(數(shù)量為16325)，SBFEM僅需較低的自由度即可達(dá)到較高精度的解.說(shuō)明了本文SBFEM模型計(jì)算考慮摩擦的滑移接觸的精確性與高效性.

圖18 裂紋面下表面σxx分布Fig.18 Stress σxxdistribution along right crack face

5 結(jié)論

本文基于比例邊界有限元法，推導(dǎo)給出了側(cè)面在任意載荷下的比例邊界有限元的新型形函數(shù)，以及相應(yīng)的SBFEM有關(guān)公式.運(yùn)用Lagrange乘子法，推導(dǎo)給出了點(diǎn)對(duì)點(diǎn)與邊對(duì)邊接觸的SBFEM求解裂紋面接觸問(wèn)題的支配方程，較為精確地模擬了裂紋面的接觸，解決了裂紋面在受壓時(shí)的粘接接觸和滑移接觸問(wèn)題.通過(guò)3個(gè)算例的數(shù)值分析，驗(yàn)證了本文方法的精確性與高效性.算例表明，本文推導(dǎo)的考慮側(cè)邊載荷存在的SBFEM形函數(shù)是合適的；在不含有裂尖單元的裂紋面接觸時(shí)，采用點(diǎn)對(duì)點(diǎn)接觸模型即可精確地模擬裂紋面的接觸條件；在含有裂尖單元的裂紋面接觸時(shí)，需要對(duì)裂尖單元的側(cè)面采用邊對(duì)邊約束接觸模型.

1 嵇醒.斷裂力學(xué)判據(jù)的評(píng)述.力學(xué)學(xué)報(bào),2016,48(4):741-753(Ji Xing.A critical review on criteria of fracture mechanics.Chinese Journal of Theoretical and Applied Mechanics,2016,48(4):741-753(in Chinese))

2 李臥東,陳勝宏.接觸摩擦問(wèn)題的數(shù)值模擬.巖土力學(xué),2003,24(3):385-388(Li Wodong，Chen Shenghong.Numerical model-ing for frictional contact problems.Rock and Soil Mechanics,2003,24(3):385-388(in Chinese))

3 高志強(qiáng),傅衛(wèi)平,王雯等.彈塑性微凸體側(cè)向接觸相互作用能耗.力學(xué)學(xué)報(bào),2017,49(4):858-869(Gao Zhiqiang,Fu Weiping,Wang Wen,et al.Study on the contact energy dissipation of the lateraland interactional between the elastic-plastic asperities.Chinese Journal of Theoretical and Applied Mechanics,2017,49(4):858-869(in Chinese))

4 Tur M,Fuenmayor FJ,Wriggers P.A mortar-based frictional contact formulation for large deformations using Lagrange multipliers.Computer Methods in Applied Mechanics and Engineering,2009,198(37):2860-2873

5 MateiA.WeaksolvabilityviaLagrangemultipliersforcontactproblems involving multi-contact zones.Mathematics and Mechanics of Solids,2016,21(7):826-841

6 Hansbo P,Rashid A,Salomonsson K.Least-squares stabilized augmented Lagrangian multiplier method for elastic contact.Finite Elements in Analysis and Design,2016,116(1):32-37

7 Benkhira EH,Essou fiEH,Fakhar R.On convergence of the penalty method for a static unilateral contact problem with nonlocal friction in electro-elasticity.European Journal of Applied Mathematics,2016,27(1):1-22

8 Bourichi S,Essou fi E.Penalty method for unilateral contact problem with Coulomb’s friction for locking material.International Journal of Mathematical Modelling&Computations,2016,6(1):61-81

9 李建宇,潘少華,張洪武.解三維摩擦接觸問(wèn)題的一個(gè)二階錐線性互補(bǔ)法.力學(xué)學(xué)報(bào),2009,41(6):869-877(Li Jianyu,Pan Shaohua,Zhang Hongwu.A second-order cone linear complementarity approach for three-dimensional frictional contact problems.Chinese Journal of Theoretical and Applied Mechanics,2009,41(6):869-877(in Chinese))

10 Song C,Wolf JP.The scaled boundary fi nite-element method—alias consistent in fi nitesimal fi nite-element cell method—for elastodynamics.Computer Methods in Applied Mechanics and Engineering,1997,147(3):329-355

11 Deeks AJ,Wolf JP.A virtual work derivation of the scaled boundary finite-element method for elastostatics.Computational Mechanics,2002,28(6):489-504

12 Song C,Wolf JP.Semi-analytical representation of stress singularities as occurring in cracks in anisotropic multi-materials with the scaled boundary fi nite-element method.Computers&Structures,2002,80(2):183-197

13 Birk C,Prempramote S,Song C.An improved continued-fractionbased high-order transmitting boundary for time-domain analyses in unbounded domains.International Journal for Numerical Methods in Engineering,2012,89(3):269-298

14 Saputra A,Talebi H,Tran D,et al.Automatic image-based stress analysisbythescaledboundary fi niteelementmethod.International Journal for Numerical Methods in Engineering,2017,109(5):697-738

15 Behnke R,Mundil M,Birk C,et al.A physically and geometrically nonlinear scaled-boundary-based fi nite element formulation for fracture in elastomers.International Journal for Numerical Methods in Engineering,2014,99(13):966-999

16 吳澤艷,王立峰,武哲.比例邊界坐標(biāo)插值方法在譜元法中的應(yīng)用——無(wú)窮域Euler方程的數(shù)值模擬.力學(xué)學(xué)報(bào),2013,45(4):619-623(Wu Zeyan,Wang Lifeng,Wu Zhe.The scaled boundary coordinate interpolation method and its application to spectral element method:numerical simulation of the Euler equations over unbounded domains.Chinese Journal of Theoretical and Applied Mechanics,2013,45(4):619-623(in Chinese))

17 薛冰寒,林皋,胡志強(qiáng)等.摩擦接觸問(wèn)題的比例邊界等幾何B可微方程組方法.力學(xué)學(xué)報(bào),2016,48(3):615-623(Xue Binghan,Lin Gao,Hu Zhiqiang,et al.Analysis of frictional contact problems by SBIGA-BDE method.Chinese Journal of Theoretical and Applied Mechanics,2016,48(3):615-623(in Chinese))

18 陳燈紅,杜成斌.基于SBFE和改進(jìn)連分式的有限域動(dòng)力分析.力學(xué)學(xué)報(bào),2013,45(2):297-301(Chen Denghong,Du Chengbin.Dynamic analysis of bounded domains by SBFE and the improved continued-fraction expansion.Chinese Journal of Theoretical and Applied Mechanics,2013,45(2):297-301(in Chinese))

19 Ooi ET,Song C,Tin-Loi F.A scaled boundary polygon formulation for elasto-plastic analyses.Computer Methods in Applied Mechanics and Engineering,2014,268(1):905-937

20 Zheng Z,Zang M,Chen S,et al.An improved 3D DEM-FEM contact detection algorithm for the interaction simulations between particles and structures.Powder Technology,2017,305(1):308-322

21 王福軍,王利萍,程建鋼等.并行有限元計(jì)算中的接觸算法.力學(xué)學(xué)報(bào),2007,39(3):422-427(Wang Fujun,Wang Liping,Cheng Jianggang,et al.A contact algorithm for parallel computation of FEM.Chinese Journal of Theoretical and Applied Mechanics,2007,39(3):422-427(in Chinese))

22 Giner E,Tur M,Taranc′on JE,et al.Crack face contact in X-FEM usingasegment-to-segmentapproach.InternationalJournalforNumerical Methods in Engineering,2010,82(11):1424-1449

23 Zheng A,Luo X.A mathematical programming approach for frictional contact problems with the extended fi nite element method.Archive of Applied Mechanics,2016,86(4):599-616

24 劉永健,姚振漢.三維彈性體移動(dòng)接觸邊界元法的一類新方案.工程力學(xué),2005,22(1):6-11(Liu Yongjian,Yao Zhenhan.A new scheme of BEM for moving contact of 3D elastic solid.Engineering Mechanics,2005,22(1):6-11(in Chinese))

25 Yang ZJ,Deeks AJ.Fully-automatic modelling of cohesive crack growth using a fi nite element–scaled boundary fi nite element coupled method.Engineering Fracture Mechanics,2007,74(16):2547-2573

26 鐘紅,宋平平.任意裂紋面載荷作用下界面斷裂分析.哈爾濱工業(yè)大學(xué)學(xué)報(bào),2016,48(2):152-157(Zhong Hong,Song Pingping.Analysis of interface crack with arbitrary crack tractions.Journal of Harbin Institute of Technology,2016,48(2):152-157(in Chinese))

27 Song C.A super-element for crack analysis in the time domain.International Journal for Numerical Methods in Engineering,2004,61(8):1332-1357

28 Deeks AJ,Wolf JP.Semi-analytical elastostatic analysis of unbounded two-dimensional domains.International Journal for Numerical and Analytical Methods in Geomechanics,2002,26(11):1031-1057

29 Song C.A matrix function solution for the scaled boundary fi niteelement equation in statics.Computer Methods in Applied Mechanics and Engineering,2004,193(23):2325-2356

30 Vu TH,Deeks AJ.Use of higher-order shape functions in the scaled boundary fi nite element method.International Journal for Numerical Methods in Engineering,2006,65(10):1714-1739

31 Ooi ET,Song C,Natarajan S.Construction of high-order complete scaled boundary shape functions over arbitrary polygons with bubble functions.International Journal for Numerical Methods in Engineering,2016,108(9):1086-1120

32 宣兆成,李興斯.接觸分析的光滑模型及迭代算法.力學(xué)學(xué)報(bào),2001,33(3):340-348(Xuan Zhaocheng,Li Xingsi.Smooth model and iteration algorithm for contact analysis.Acta Mechanica Sinica,2001,33(3):340-348(in Chinese))

33 王勖成.有限單元法.北京：清華大學(xué)出版社，2013(Wang Xucheng.Finite Element Method.Beijing:Tsinghua University Press，2013(in Chinese))

34 Dai S,Augarde C,Du C,et al.A fully automatic polygon scaled boundary fi nite element method for modelling crack propagation.Engineering Fracture Mechanics,2015,133(1):163-178

35 Ooi ET,Song C,Tin-Loi F,et al.Polygon scaled boundary fi nite elements for crack propagation modelling.International Journal for Numerical Methods in Engineering,2012,91(3):319-342

36 Deeks AJ,Augarde CE.A meshless local Petrov-Galerkin scaled boundary method.Computational Mechanics,2005,36(3):159-170

37 莊茁.基于ABAQUS的有限元分析和應(yīng)用.北京：清華大學(xué)出版社,2008(Zhuang Zhuo.Based on ABAQUS Finite Element Analysis and Application.Beijing:Tsinghua University Press，2008(in Chinese))

38 徐芝綸.彈性力學(xué).北京：高等教育出版社，2006(XuZhilun.Elasticity.Beijing:Higher Education Press,2006(in Chinese))

CRACK FACE CONTACT PROBLEM ANALYSIS USING THE SCALED BOUNDARY FINITE ELEMENT METHOD1)

Zhang Peng Du Chengbin2)Jiang Shouyan
(Department of Engineering Mechanics,Hohai University,Nanjing211100,China)

In the case of arbitrary tractions on the side faces of the crack,a polynomial function of the radial coordinate can be employed to describe the side face loads in the scaled boundary fi nite element method(SBFEM).The SBFEM new shape function considering the side face loads is presented.The corresponding sti ff ness matrix and equivalent node load is derived together based on the SBFEM new shape function.The model of crack face contact using the SBFEM is proposed in this paper for the fi rst time.Lagrange’s multiplier method is used to establish the contact constraints of contact model between crack faces.The governing equations for the nonlinear surface contact problems in SBFEM is derived,including adhesion contact problems and sliding friction problems.The elements where the crack faces lie are divided into non crack tip elements and the crack tip element.For the former,the crack faces act as the boundary of the SBFEM element,the contact tractions on the boundary can be assigned to the nodes equivalently and the Lagrange’s multiplier is applied for the point constraints.For the latter,the interpolation fi eld of Lagrange’s multiplier is constructed on the whole side faces.The Lagrange’s multiplier is assumed to be linear along the side faces，the segment constraint approach is proposed to optimize the ful fi lment of the contact constraints along the crack faces.By comparing the results of the calculation of analytical solution and software ABAQUS on three di ff erent numerical contact problems of crack faces,the accuracy and e ff ectiveness of the proposed point-to-point and segment-to-segment contact model for fracture surfaces contact problems is veri fi ed in this paper.

scaled boundary fi nite element method，shape functions，crack，point constraints，segment constraints

O346.1

A doi：10.6052/0459-1879-17-195

2017–05–23 收稿，2017–08–07 錄用,2017–08–11 網(wǎng)絡(luò)版發(fā)表.

1)國(guó)家自然科學(xué)基金資助項(xiàng)目(11372098,51579084,51309088).

2)杜成斌，教授，主要研究方向：水工結(jié)構(gòu)工程中的力學(xué)問(wèn)題.E-mail:cbdu@hhu.edu.cn.

章鵬，杜成斌，江守燕.比例邊界有限元法求解裂紋面接觸問(wèn)題.力學(xué)學(xué)報(bào),2017,49(6):1335-1347

Zhang Peng,Du Chengbin,Jiang Shouyan.Crack face contact problem analysis using the scaled boundary fi nite element method.Chinese Journal of Theoretical and Applied Mechanics,2017,49(6):1335-1347