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

章 鵬 杜成斌 江守燕

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

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

章 鵬 杜成斌2)江守燕

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

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

比例邊界有限元,形函數,裂紋,點對點接觸,邊對邊接觸

引言

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

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

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

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

1.1 考慮側面載荷的比例邊界多邊形有限元的形函數

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

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

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

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

對于裂紋單元內任意一點，相對于比例中心的局部坐標(x,y)T可用ξ和s表示為

其中，(x0y0)為邊界上的節點坐標.N(s)表示為沿s方向的形函數

其中n為邊界節點個數.

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

該方程為徑向上的平衡方程.其中E0,E1,E2為單元的系數矩陣，Ft(ξ)為側面載荷，u(ξ)為節點位移函數.

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

其中q(ξ)為內部節點力，Z為2n×2n階Hamiltonian矩陣，滿足

采用Schur分解[29]

得到單元的模態位移ψd，模態載荷ψq及對應的特征值矩陣Λλ.

裂紋側面上的任意載荷可以分解為關于ξ的M階多項式函數[25]

式中Wi為單元側邊力分布矩陣，ai為第i階多項式系數.

側面載荷的第i階位移模態[11]為

則徑向位移解為

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

式(10)的矩陣形式為

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

邊界上的節點位移為

其中積分常數c用邊界位移表達為

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

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

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

比例邊界有限元在環向s上采用與有限單元法中形函數類似方法，通過N(s)進行插值，即

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

令

可得

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

對于不考慮側面載荷的比例邊界元形函數為[16]

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

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

SBFEM的應變矩陣[11]為

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

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

引入應變位移矩陣

應變矩陣可寫為

應力矩陣為

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

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

考慮面載荷和側邊載荷時，采用虛功原理的表達式為

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

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

式(31)等價為

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

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

將dV進行積分得到

定義矩陣Y

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

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

定義矩陣X

采用分部積分，式(38)可化簡為Lyapunov方程[31]

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

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

式(32)中等式右側為等效節點載荷Fp

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

2 接觸問題的有關公式

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

圖3 接觸點對與點對間的距離Fig.3 The contact points and distance between the points

同理定義切向相對距離gT

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

接觸問題可描述為求區域內位移場，使得系統的勢能達到最小.Lagrange乘子法[4]是求解接觸約束最小化問題的常用方法之一，通過引入Lagrange乘子將接觸問題轉化為無約束問題求解.接觸狀態有3種：張開、粘結、滑移.法向Lagrange乘子場λN表示當接觸發生時迫使gN等于0(裂紋面接觸但無嵌入)的壓力.在裂紋面為粘結接觸時，切向Lagrange乘子場λT表示在粘結區迫使滿足粘結條件的切向力，在裂紋面發生滑移接觸時，切向Lagrange乘子場λT表示為滑移的摩擦力.

應用Lagrange乘子法在物體表面考慮摩擦接觸條件時可導出3個方程：

(1)接觸不嵌入條件

(2)在切向滿足粘結條件時，即

在切向滿足有摩擦的滑動接觸狀態時：切向運動不再受約束，但是切向力滿足

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

3 基于比例邊界有限元求解摩擦接觸問題

3.1 接觸約束公式

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

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

圖4 點對點與邊對邊約束示意圖Fig.4 Point-to-segment constraints diagram

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

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

3.1.1 不含裂尖的裂紋面單元

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

其中，CNk為矩陣CN第k行，Cindex為點對P和Q兩節點在整體節點中的坐標變換.

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

同樣在切向方向，令

可得裂紋面點對接觸切向方向約束方程

3.1.2 包含裂尖的裂紋面單元

將含有裂尖單元的比例中心放在裂尖處，如圖4的單元3所示.在該單元的裂紋面側面(side faces)上，ξ從0變化到1，該單元所受的裂紋面接觸應力為側面應力.該側面的形函數Nt(ξ,s)為關于(ξ,s)的函數，因為此時裂紋面為非邊界面(ξ(0 6ξ6 1))，按常規思路施加點對約束不能正確地表示裂紋面接觸.為此本文提出了在裂尖單元中采用邊對邊接觸方法，即在整個側面上構造Lagrange乘子的插值場，Lagrange乘子呈線性變化，其插值可表示為

其中,λj1為側邊中心(裂尖中心)的接觸應力，λj2為側邊端點的接觸應力.

Mi為Lagrange乘子插值函數

則式(52)可寫為

由式(19)形函數可知

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

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

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

同理可得到裂紋側面邊對邊接觸切向方向約束方程

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

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

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

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

故滿足滑動接觸狀態的控制方程為

4 數值算例

4.1 貫穿裂紋板接觸

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

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

圖6 多邊形比例邊界有限元網格Fig.6 Polygon SBFEM mesh

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

為了驗證該方法精確性與收斂性，采用結構應變能相對誤差作為指標.應變能誤差的表達式[36]可以寫為應變能的相對誤差為

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

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

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

圖8 應變能相對誤差隨網格自由度變化Fig.8 The relative error in the strain energy varies with the grid degrees of freedom

4.2 含有圓孔的縫內為非均勻壓力的裂紋板接觸

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

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

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

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

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

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

圖11 結構變形圖(放大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兩種方法計算的應變能相對誤差隨不同網格自由度的變化.由圖中可看出，基于SBFEM求解的結果在自由度數較小情況下，誤差較小(相對誤差為0.0032%)，而基于ABAQUS的求解誤差在自由度數較小情況下所得誤差較大(相對誤差為1.7%)，隨著網格的加密，自由度數的增大，誤差逐漸減小，但ABAQUS求解精度遠小于SBFEM求解精度.這說明SBFEM點對點和邊對邊約束方法在非均布應力下的裂紋接觸仍具有較高的精度和效率.

圖13 應變能相對誤差隨網格自由度變化Fig.13 The relative error in the strain energy varies with the grid degrees of freedom

4.3 考慮摩擦的滑移接觸

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

圖14 多邊形比例邊界有限元網格Fig.14 Polygon SBFEM mesh

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

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

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

本算例沒有解析解，仍采用商業有限元軟件ABAQUS模擬結果來對比本文提出的SBFEM求解裂紋面接觸方法的精度和效率.在使用ABAQUS模擬過程中，分別采用粗網格、中等網格和細網格求解該滑移接觸問題，圖17分別給出了ABAQUS不同網格的模型.

圖17 ABAQUS不同網格模型Fig.17 Di ff erent mesh models of ABAQUS

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

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

5 結論

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

1 嵇醒.斷裂力學判據的評述.力學學報,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 李臥東,陳勝宏.接觸摩擦問題的數值模擬.巖土力學,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 高志強,傅衛平,王雯等.彈塑性微凸體側向接觸相互作用能耗.力學學報,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 李建宇,潘少華,張洪武.解三維摩擦接觸問題的一個二階錐線性互補法.力學學報,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 吳澤艷,王立峰,武哲.比例邊界坐標插值方法在譜元法中的應用——無窮域Euler方程的數值模擬.力學學報,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 薛冰寒,林皋,胡志強等.摩擦接觸問題的比例邊界等幾何B可微方程組方法.力學學報,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和改進連分式的有限域動力分析.力學學報,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 王福軍,王利萍,程建鋼等.并行有限元計算中的接觸算法.力學學報,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 劉永健,姚振漢.三維彈性體移動接觸邊界元法的一類新方案.工程力學,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 鐘紅,宋平平.任意裂紋面載荷作用下界面斷裂分析.哈爾濱工業大學學報,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 宣兆成,李興斯.接觸分析的光滑模型及迭代算法.力學學報,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 王勖成.有限單元法.北京：清華大學出版社，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的有限元分析和應用.北京：清華大學出版社,2008(Zhuang Zhuo.Based on ABAQUS Finite Element Analysis and Application.Beijing:Tsinghua University Press，2008(in Chinese))

38 徐芝綸.彈性力學.北京：高等教育出版社，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 網絡版發表.

1)國家自然科學基金資助項目(11372098,51579084,51309088).

2)杜成斌，教授，主要研究方向：水工結構工程中的力學問題.E-mail:cbdu@hhu.edu.cn.

章鵬，杜成斌，江守燕.比例邊界有限元法求解裂紋面接觸問題.力學學報,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