三維切口/裂紋結構的擴展邊界元法分析1)

李 聰 牛忠榮 胡宗軍 胡 斌

*(合肥工業大學土木與水利工程學院力學系，合肥 230009)

?(安徽建筑大學土木工程學院，合肥 230601)

引言

在土木、機械和航空等工程結構中經常出現三維V 形切口/裂紋問題，如地下室方形洞室角區、重力壩踵壩趾區等.由于幾何形狀或材料組成的突變會造成三維V 形切口/裂紋尖端的強應力集中，嚴重時會導致受力構件甚至結構的整體破壞.因而，準確計算三維V 形切口/裂紋結構的應力場具有重要的理論意義和工程實用意義.

目前國內外學者對于三維V 形切口/裂紋結構的研究大多局限于裂尖場的奇異項、應力強度因子、裂紋擴展路徑等.Chue 和Liu[1]給出了正交各向異性切口應力奇異階的一般解，并指出解的奇異性取決于切口角度、邊界條件及材料特性.Ungamornrat[2]用弱奇異對稱Galerkin 邊界元法(SGBEM)分析了三維多相材料線彈性裂紋結構的應力奇異性.Sator 和Becker[3]采用復變函數法獲得了各向同性多材料接頭端部應力奇異性指數的解析解.胡宗軍等[4]建立了三維位勢問題高階邊界單元幾乎強奇異和幾乎超奇異積分的半解析算法，其數值結果表明該方法可高效計算邊界元法中幾乎奇異面積分.Wu 等[5]基于離散的快速傅里葉變換(FFT)理論，推導出線彈性各向異性雙材料結構的界面應力場的半解析解.

此外，對應力強度因子和裂紋擴展的研究也比較深入.文龍飛等[6]基于改進型擴展有限元法，研究了動載荷作用下擴展裂紋尖端應力強度因子的求解方法.Xu 等[7]在變長徑比和雙材料參數下，采用體力法計算了三維矩形界面裂紋的應力強度因子.王振和余天堂[8]采用自適應多尺度擴展有限元法計算了三維I 型裂紋和I-II 復合型裂紋的應力強度因子.賈旭等[9]采用裂紋彈性應力場的有限元解計算了三維裂紋的應力強度因子，其結果符合三維裂紋體裂紋尖端實際的應力應變狀態.Fakovr 和Ghoreishi[10]研究了機械支柱、轉速、裂紋取向等因素對旋轉盤中應力強度因子的影響.Huang 等[11]采用有限塊法(FBM)計算了靜載荷和瞬態動態載荷下的應力強度因子和t應力.賈旭等[12]提出了三維矩形平板內承受復雜非線性載荷下，僅含有3個系數的通用權函數，并采用FEM 計算了三維含穿透裂紋的矩形板的3 組參考應力強度因子.Magnus[13]采用FEM 分析了三維各向異性低滲透巖石在水壓作用下的應力場.Dhanesh 等[14]采用迭代分析技術MMEKM 分析了三維彈性復合材料層合板在均勻軸向拉伸、彎曲、扭轉和熱載荷下的應力場，結果顯示了自由邊應力場的MMEKM 解在項數和迭代次數方面具有良好的收斂特性.Stasyuk[15]基于改進的邊界積分方程法，分析了三維彈性無限體和有限彈性體中兩個相互作用的圓盤形裂紋附近三維應力強度因子，并給出了裂紋沿法向均勻加載下的應力場.

Karsten 和Gunther[16]采用對偶邊界元法模擬了復雜三維線彈性裂紋的裂紋擴展路徑，其結果和實驗結果吻合較好.Li 和Guo[17]采用三維FEM 研究了有限厚度板在單軸拉應力作用下的V 形切口尖端彈性應力場，研究結果表明了板厚、切口角度對應力集中系數有明顯的影響等.Park 和Nikishkov[18]提出了SGBEM-FEM 交替方法，該方法可有效分析三維平面和非平面裂紋及其擴展.Lan 等[19]采用三維FEM 模擬了鋁合金板在混合型加載下的裂紋擴展過程，并研究了裂紋尖端有限元尺寸和裂紋擴展增量尺寸對模擬預測結果的影響.Dias 等[20]采用有限元方法模擬了準脆性材料的三維裂紋擴展問題，該方法可高效模擬位移不連續擴展的性能.Ding 和Wang[21]采用三維邊界層(小范圍屈服)公式對三維小應變進行了彈塑性模擬，模擬結果表明了尖端區域非奇異T 應力對塑性區的尺寸和形狀有顯著的影響，且應力的大小、方向也受到T 應力的影響.Wang 等[22]通過建立三維復合材料的2 種幾何結構模型：微觀結構模型和多單元模型，分析了三維碳/環氧復合材料在-100°C 至-140°C 不同溫度場下的熱膨脹行為和界面熱應力.

研究三維線彈性V 形切口/裂紋的文獻很多，但大多局限在研究第1 階應力奇性指數及其在此基礎上模擬裂紋擴展路徑，而關于三維切口/裂紋結構完整位移和應力場的分析缺乏有效方法.本文基于二維V 形切口/裂紋結構的全域彈性應力場[23-24]的求解思想，聯合尖端扇形柱的漸近分析和外圍結構的邊界元法分析求解三維V 形切口/裂紋結構完整的位移和應力場，稱之為擴展邊界元法(XBEM).XBEM 計及了離尖端徑向距離r的漸近級數展開式中若干項的特征指數和特征函數.文中將通過算例定量展現XBEM 求解三維V 形切口/裂紋結構完整應力場的準確性和適用性.

1 線彈性三維V 形切口/裂紋尖端附近應力奇異性分析

考慮三維柱狀V 形切口結構，如圖1 所示，其柱形切口張角為α=2π-θ1-θ2.當α=0°時，圖1 所示結構為三維裂紋結構.將三維V 形切口結構分為一個半徑r=ρ 的扇形柱體(圖2(b))和外圍結構(圖2(a))兩部分.外圍結構新形成的邊界為，在切口尖端處同時定義笛卡爾坐標系ox1x2x3和柱坐標系orθz(坐標原點o在切口尖端).扇形柱體是應力奇異性區域，外圍結構是無應力奇異性的普通彈性區域.對于一個三維V 形切口結構，圖2(b)所示切口尖端附近應力奇異性區域的位移場可表達成關于徑向r的漸近級數展開式[25-26]

式中，Ak為位移幅值系數，λk為切口尖端的位移特征指數，N表示截取的級數項數，(θ,z,λk)，(θ,z,λk)，(θ,z,λk) 為位移特征角函數.由于圖2(b)結構形狀沿z方向是相同的，并不涉及外載荷，因此分析切口尖端扇形柱區域的奇異性時，式(1)中的特征角函數(θ,z,λk)，(θ,z,λk)，(θ,z,λk)與z無關，于是式(1)可簡化寫成

圖1 三維V形切口結構Fig.1 A 3-D V-notched structure

圖2 三維V形切口結構挖去小扇形柱體Fig.2 A sectoral column is separated from the 3-D V-notched structure

考慮各向同性均質材料，將式(2)代入三維彈性力學幾何方程和本構方程，可用位移場函數表示應力場函數(為簡明起見，以下分別簡寫為

式中，(···)′=d(···)/dθ，余下同.G=E/[2(1+v)]，E為材料楊氏模量，v是泊松比.

注意到式(1)的位移特征指數λk和位移特征角函數僅取決于切口的材料性質和楔形邊界條件，不依賴載荷，因此應力特征分析不用考慮體力.取式(2)中的典型項Akrλk+1(θ,λk)，Akrλk+1(θ,λk)，Akrλk+1(θ,λk)代入到忽略體力項的三維彈性力學平衡方程

通過一系列推導，可得關于切口尖端扇形柱區域應力奇異性特征值問題的常微分方程組(為方便起見，分別簡寫為)

式中，(···)′′表示對坐標θ 的二階導數，余類同.X=Eν/[(1+ν)(1-2ν)]=2Gν/(1-2ν).

假定切口兩楔邊自由，則其上面力為零，邊界條件可由式(6)表示

因此，三維V 形切口尖端扇形柱附近的位移特征指數λk的計算變成了求解常微分方程組(5)和相應的邊界條件(6)的特征值問題.采用插值矩陣法[27-28]求解式(5)和式(6)，可得切口尖端扇形柱附近的位移特征指數λk、相應的位移特征角函數及其導函數.

式中，下標R 表示復數的實部，下標I 表示復數的虛部.若λk是實根，則其相應虛部項為零，系數ζ=1；如果λk是共軛復根，則相應虛部項不為零，系數ζ=2，因需計及復數共軛項對位移場和應力場的貢獻.

注意：式(7)和式(8)中位移幅值系數AkR和AkI(k=1,2,···,N)待求.基于切口尖端應力奇異性的特征分析結果式(7)和式(8)，下面建立求解三維切口結構位移和應力漸近場表達式中幅值系數AkR和AkI的擴展邊界元法.

2 分析三維V形切口/裂紋結構完整應力場的擴展邊界元法

在圖1 切口結構中挖去尖端扇形柱后剩余的外圍結構區域Ω′無應力奇異性，見圖2(a).其邊界位移um(x)和面力分量tm(x)可由如下常規邊界積分方程描述

外圍結構Ω′內任意點的位移和應力可分別由位移和應力邊界積分方程表達

3 三維V 形切口/裂紋問題算例

3.1 含V 形切口的長方形柱體

含V 形切口的長方形柱體，見圖3.結構尺寸L=15 mm，H=5 mm，B=10 mm，彈性模量E=200 GPa，泊松比ν=0.3.沿長方形柱體正表面對稱線上有一個等深度的V 形切口，切口深度a=2.5 mm，切口張角2α=60°.長方形柱體結構左側面(x2=L)固定，右側面(x2=-L)施加均勻的外部拉力qx=0 MPa，qy=30 MPa.

首先從三維V形切口結構(如圖3 所示)中挖取一個ρ=0.5 mm 的扇形柱體(如圖4(b)所示)，采用插值矩陣法對扇形柱體進行特征分析，獲得三維V形切口尖端扇形柱區域的前若干階位移特征指數λk和其相對應的位移特征角函數,,，如表1 和圖5 和圖6 所示.

圖3 三維V形切口結構Fig.3 A 3-D V-notched structure

圖4 三維V形切口挖去小扇形柱體Fig.4 A sectoral column is separated from the 3-D V-notched structure

表1 三維V 型切口2α=60° 尖端的位移特征指數λkTable 1 The displacement indexes λk of the 3-D V-notched structure with 2α=60°

圖5 三維V 形切口2α=60° 對應的位移特征角函數Fig.5 The displacement eigenvectors for the 3-D V-notched structure when 2α=60°

圖5 三維V 形切口2α=60° 對應的位移特征角函數(續)Fig.5 The displacement eigenvectors for the 3-D V-notched structure when 2α=60° (continued)

圖6 三維V 形切口2α=60° 對應的位移特征角函數Fig.6 The displacement eigenvectors for the 3-D V-notched structure when 2α=60°

圖6 三維V 形切口2α=60° 對應的位移特征角函數(續)Fig.6 The displacement eigenvectors for the 3-D V-notched structure when 2α=60° (continued)

表1 中λk前4 項是剛體位移項，相應的特征角函數見圖5.圖5(a)～圖5(c)呈現的是平面問題的3 個剛體位移項，圖5(d)呈現的是反平面問題的剛體位移項.這些剛體位移項對應力場的計算沒有影響，但對于位移場的計算不可省略.

注意到圖6(b)、圖6(d)、圖6(h)中?ur和?uθ均為0，故λ6,λ8和λ10為反平面V 形切口對應的位移特征指數.圖6(a)、圖6(c)、圖6(e)、圖6(f)中?uz均為0，故λ5,λ7和λ9為平面V 形切口對應的位移特征指數；其中，λ5和λ9為平面內對稱位移特征指數，λ7為平面內反對稱位移特征指數.

擴展邊界元法(XBEM)對三維V 形切口結構進行計算時，初始選取漸近級數展開式(1)中截取級數項數N=7(4 個剛體位移項，1 個平面內對稱項，1 個平面內反對稱項和1 個反平面項)作為漸近級數展開式(1)的基準項數.隨后在N=7 的基礎上增大N，討論N的取值對三維V 形切口結構完整位移和應力場精度的影響.

XBEM 分析三維V 形切口結構的實施過程：首先，采用插值矩陣法求解常微分方程組(5)獲得切口尖端的位移特征指數λk及其相對應的特征角函數，然后對外圍結構采用邊界元法進行分析，邊界采用8 節點四邊形單元，單元數424，節點數1274.聯立求解兩組方程組獲得了位移幅值系數AkR,AkI以及外圍結構邊界未知的位移和面力.在XBEM 執行中，不可避免遇到有些源點到邊界單元距離很小的情形—三維BEM 中幾乎奇異積分的難題，本文采用半解析法[4,30]高效解決了這一難題.最后，將求得的邊界位移和面力代入式(10)和式(11)得外圍結構任意內點的位移和應力.對于扇形柱體區域的位移和應力場，采用漸近級數展開式(1)和式(3)計算，外圍區域(圖4(a))內點的位移和應力場采用內點積分方程(10)和式(11)計算，其臨近扇形柱附近的外圍區域內點的位移和應力場也可采用漸近級數展開式(1)和式(3)計算.

表2 和表3 以及圖7 分別給出了切口尖端扇形柱一定范圍內應力和位移的XBEM 計算結果，其中，XBEM-ASY 表示采用漸近級數展開式(1)和式(3)計算的結果(扇形柱及其臨近區域)，XBEMBEM 表示采用積分方程式(10)和式(11)計算的結果(外圍區域)，FEM 表示采用有限元分析軟件計算的結果.表2 和表3 中XBEM-BEM 計算值和圖7(a)曲線顯示：當截取級數項數N＞7 時，N的增加對于外圍結構區域的XBEM-BEM 計算結果影響很小.說明選取了足夠的N之后，XBEM 的計算精度對N的取值不敏感，穩定性很好，這對于XBEM 求解三維V 形切口結構的完整應力場非常便利.

表2 顯示，對外圍區域點(1.438 mm,0 mm,5 mm)處，當N=7 時，XBEM-ASY 計算的應力值(σx2=32.56 MPa)和XBEM-BEM 計算相同點的應力值(σx2=40.04 MPa)的相對差為-18.7%.當N=8時，XBEM-ASY 計算的應力值(σx2=34.89 MPa)和XBEM-BEM 計算相同點的應力值(σx2=40.04 MPa)的相對差為-12.9%.當N=9 時，XBEM-ASY 計算的應力值(σx2=36.57 MPa)和XBEM-BEM 計算相同點的應力值(σx2=40.04 MPa)的相對差為-8.7%.說明相對XBEM-BEM計算值而言，XBEM-ASY 計算同一點的相對差隨N的增大逐漸減小，即XBEMASY 計算切口尖端應力場的有效范圍隨N的增大逐漸增大.

通常情況，外圍結構的位移和應力由XBEMBEM 解勝任，無需采用XBEM-ASY 計算，這里關于兩者計算值的比較，意在說明XBEM-ASY 解的有效范圍可為扇形柱半徑的選取提供參考.例如，若XBEM-ASY 計算的切口尖端應力場在r＜1.0 mm 的范圍內是足夠準確的，說明XBEM 使用時，尖端扇形區域半徑ρ 的選取范圍較為寬松，ρ ＜1.0 mm 之間均可行.表3 中切口結構位移解顯示了同樣的規律.

表2 XBEM和FEM 計算沿x2=0 mm, x3=5 mm 附近內點的應力分量σx2(MPa)Table 2 Variation of stress component σx2(MPa)along x2=0 mm, x3=5 mm

表3 XBEM和FEM 計算沿x2=0 mm, x3=5 mm 附近內點的位移分量ux2(mm)Table 3 Variation of displacement component ux2(mm)along x2=0 mm, x3=5 mm

圖7 XBEM-BEM 和FEM 計算3-D V 形切口沿x1=2 mm, x3=5 mm 線段的σx2 和ux2 值Fig.7 σx2 and ux2 along x1=2 mm, x3=5 mm calculated by XBEM-BEM and FEM

此外，本文還采用有限元法(FEM)對圖3 結構進行計算.這里FEM 采用20 節點六邊形單元，劃分單元數111 000.表2 顯示，當x2=0 mm，x3=5 mm，x1＞ 1.262 mm 的遠場區域，FEM 的結果和XBEM-BEM 的結果吻合較好；而當x2=0 mm，x3=5 mm，x1＜ 1.262 mm 時(切口尖端近場區域)，FEM 計算的應力場與XBEM-BEM 的結果有較大的差異.對于三維V 形切口尖端的強應力集中，盡管FEM 在尖端使用了細密的網格劃分，其采用的單元數為XBEM 單元數的261 倍，極大地增加了計算量和計算機內存需求，但仍無法準確描述其尖端近場區域的應力場.隨著距尖端距離的增大，應力集中的現象逐漸減弱，FEM 的計算結果方為可信.

3.2 含裂紋的長方形柱體單向受拉和受剪

含裂紋的各向同性均質材料長方形柱體，結構尺寸L=10.0 mm，H=3 mm，B=4 mm，彈性模量E=200 GPa，泊松比ν=0.3.沿長方形柱體正表面對稱線上有一個等深度的裂紋，裂紋深度a=0.8 mm，邊界條件為長方形柱體結構左側面固定，右側面施加均布載荷qx=10 MPa，qy=10 MPa，為復合型裂紋，如圖8 所示.

圖8 三維裂紋結構Fig.8 A 3-D cracked structure

圖9 三維裂紋挖去小扇形柱體Fig.9 Asectoral column is separated from the 3-D cracked structure

首先從圖8 三維復合型裂紋結構中挖取一個ρ=0.2 mm 的扇形柱體(如圖9(b)所示)，采用插值矩陣法對扇形柱體進行特征分析，獲得裂紋尖端扇形柱區域的前若干階位移特征指數λ1,λ2,λ3,···,λ10分別為-1,-1,-1,0,-0.5,-0.5,-0.5,0,0.000 01,0.5.其相對應的位移特征角函數?ur,?uθ,?uz如圖10 所示.

圖10 三維裂紋對應的位移特征角函數Fig.10 The displacement eigenvectors for the 3-D cracked structure

λ1～ λ4前4 項是剛體位移項，λ5～ λ10為主導的各階位移特征指數，其對應的位移特征角函數見圖10.注意到圖10(b)、圖10(d)、圖10(h)為反平面裂紋對應的位移特征角函數；圖10(a)、圖10(c)、圖10(e)、圖10(f)為平面裂紋對應的位移特征角函數.

圖11 三維裂紋結構邊界元網格Fig.11 The BE meshes of 3-D cracked structure

表4 和表5 以及圖12 結果顯示：選取不同的截取級數項數N，XBEM-BEM 計算相同點的應力和位移值基本相等，說明XBEM 計算三維復合型裂紋結構應力場和位移場的數值穩定，XBEM 的計算精度對截取級數項數N的選取不敏感，便于XBEM 的應用.

表4 還顯示，在外圍區域點(1.70 mm,0 mm,3 mm)處，當N=7 時，XBEM-ASY 計算的應力值(σx1=4.665 MPa)和XBEM-BEM 計算相同點的應力值(σx1=6.407 MPa)的相對差為-27.2%.當N=8時，XBEM-ASY 計算的應力值(σx1=5.314 MPa)和XBEM-BEM 計算相同點的應力值(σx1=6.417 MPa)的相對差為-17.2%.當N=9 時，XBEM-ASY 計算的應力值(σx1=5.544 MPa)和XBEM-BEM 計算相同點的應力值(σx1=6.527 MPa) 的相對差為-15.1%.說明相對XBEM-BEM 計算值而言，XBEMASY 計算同一點應力的相對誤差隨N的增大逐漸減小，即XBEM-ASY 計算裂紋尖端應力場的有效范圍隨N的增大逐漸增大.同理分析表5 可得，XBEM-ASY 計算裂紋尖端位移場的有效范圍隨N的增大逐漸增大.

表4 XBEM和FEM 計算沿x2=0 mm, x3=3 mm 附近內點的應力分量σx1(MPa)Table 4 Variation of σx1(MPa)along x2=0 mm, x3=3 mm calculated by XBEM and FEM

表5 XBEM和FEM 計算沿x2=0 mm, x3=3 mm 附近內點的位移分量ux1(mm)Table 5 Variation of ux1(mm)along x2=0 mm, x3=3 mm calculated by XBEM and FEM

圖12 XBEM-BEM 和FEM計算沿x1=2.88 mm, x3=3 mm 附近內點的應力分量Fig.12 The stress component along x1=2.88 mm, x3=3 mm calculated by XBEM-BEM and FEM

通常情況，外圍結構的位移和應力由XBEMBEM 解勝任，無需用XBEM-ASY 計算.這里關于兩者計算值比較，意在說明XBEM-ASY 解的有效范圍內，均可作為扇形柱半徑的選擇范圍.例如，若XBEM-ASY 計算的裂紋尖端應力場在r＜0.2 mm 的范圍內是足夠準確的，說明XBEM使用時，尖端扇形區域半徑ρ 的選取范圍較為寬松，ρ ＜0.2 mm 之間均可行.

表5 顯示，當x2=0 mm，x3=3 mm，x1＞1.38 mm 的遠場區域，FEM 的結果和XBEM-BEM 的應力結果吻合較好，而當x2=0 mm，x3=5 mm，x1≤1.38 mm 時，FEM計算的應力場(裂紋尖端近場區域)與XBEM-BEM 的結果有較大的差異.對于三維裂紋尖端的強應力集中，盡管FEM 在尖端使用了細密單元，結構采用20 節點六邊形單元，劃分單元數141 000，XBEM 使用單元數808，FEM 采用的單元數為XBEM 的163 倍，這極大地增加了計算量和計算機內存需求，但仍無法準確描述其尖端近場區域的應力場.隨著距尖端距離的增大，應力集中的現象逐漸減弱，FEM的計算結果變得可信.而本文方法可采用較少的單元數來準確獲取三維裂紋結構完整的應力場，進一步表明了XBEM 處理V 形切口應力集中問題的優越性.

4 結論

本文建立了三維線彈性V 形切口/裂紋結構完整應力場分析的擴展邊界元法，并對三維V 形切口/裂紋結構在典型載荷工況下進行了計算，結論如下.

(1) 將三維柱狀V 形切口/裂紋結構分離為尖端扇形柱區域和外圍結構區域，擴展邊界元法發揮了切口尖端扇形柱區域漸近分析和外圍結構邊界元法分析特長.XBEM 類似于半解析方法，在切口尖端扇形柱區域位移和應力由漸近級數展開式表達，遠端區域的力學場由邊界元法分析.文中算例結果顯示了XBEM 可準確獲取三維V 形切口/裂紋結構完整的位移和應力場，有利用于判定三維裂紋結構的起裂擴展.

(2)尖端扇形柱區域的漸近級數展開式選取的截取級數項數N越大，則漸近級數展開式計算尖端應力和位移場的有效范圍越大，說明尖端扇形柱半徑的選取范圍越大.

(3) 按照線彈性理論分析，三維V 形切口/裂紋尖端附近的應力值隨r→0 趨近無限大，尖端區域發生了塑性變形.后續將研究考慮三維V 形切口/裂紋結構尖端附近塑性變形的擴展邊界元法分析.