裂紋面分布加載裂尖SIFs分析的廣義參數Williams單元*

徐 華， 曹 政， 鄒云鵬， 楊綠峰

（1.廣西大學 土木建筑工程學院 工程防災與結構安全教育部重點實驗室，南寧 530004；2.廣西大學 土木建筑工程學院 廣西防災減災與工程安全重點實驗室，南寧 530004；3.廣西桂能工程咨詢集團有限公司，南寧 530015）

引 言

工程結構通常是帶裂縫工作的，如水文環境中的混凝土壩、海洋工程、地下工程和壓力容器等，其結構表面的裂縫易遭受流體侵入，有可能加速裂縫擴展與連通，引起結構的局部強度和整體承載力降低，甚至引發結構失效破壞，這些重大工程一旦破壞將引起難以估量的損失.因此，對裂紋面受荷的工程結構斷裂行為研究具有十分重要的現實意義.

裂紋面受荷時，其裂尖附近應力場和位移場較復雜，難以定量描述，通常仍以SIFs作為裂尖附近應力-應變場強弱程度的度量，并可據此進一步判斷結構的安全性.基于此，國內外學者圍繞裂紋面加載的裂尖SIFs開展了大量研究工作.劉鈞玉等[1-2]、Zhong等[3]和陳白斌等[4]基于比例邊界有限元法，分別對裂紋面受拉伸荷載的單邊裂紋、多裂紋和界面裂紋的裂尖SIFs進行了分析，該方法相較于普通有限單元法，需引入多個系數矩陣，增加了計算工作量.李亞等[5]采用線場分析方法，對無限大板中心裂紋面受均布荷載的應力場進行修正，進而得到了有限寬板中心裂紋的裂尖SIFs，但不同的裂紋模型均需相對應的應力函數.Fett等[6]利用ABAQUS有限元軟件擬合半無限大板邊界斜裂紋面上受集中力的權函數，以此求解了裂尖復合型SIFs.文獻[5-6]針對不同的裂紋模型，均需找到相對應的函數才能求解SIFs，不具有通用性.Li等[7]采用柔度法求解有限寬板邊界裂紋面受分布荷載的裂尖SIFs，需通過疊加法對原始模型進行拆分處理，不能直接求解.Walters等[8]針對三維模型中裂紋面受荷的平面裂紋SIFs，基于相互作用積分法分離出獨立的裂紋面荷載積分項，提高了計算精度.Muthu等[9]提出了一種基于擴展無網格Galerkin法的裂紋閉合積分理論，通過求解應變能釋放率獲得裂尖SIFs的方法.在實際工程方面，賈金生等[10]基于斷裂力學理論估算了重力壩壩踵水平裂縫中含高壓水的裂尖SIFs，但基本假設較為粗略，與實際結果存在較大誤差.唐世斌等[11]建立了在土壓力、裂紋面和井筒內水壓耦合作用下的射孔裂尖SIFs計算方程.綜上所述，當前的裂紋面受荷問題仍以有限單元法為主，但該方法本身仍存在較大的改進空間.廣義參數有限元法[12]是有限單元法的一個分支，該方法通過整體控制方程求解出與SIFs直接相關的廣義參數，避免了線性外推擬合等繁瑣的后處理過程引起的二次誤差，從而具有很高的精度.但該方法建立的W單元需滿足奇異區內裂紋面自由的邊界條件，即 σφ=0，τρφ=0(φ=±π)，而裂紋面受荷時不滿足該邊界條件，仍需改進.

針對同樣不滿足邊界條件的曲線裂紋[13]，課題組在裂尖建立等效區，將等效區內的曲線轉換為折線以滿足W單元的邊界條件，受此思想啟發，對裂紋面直接受荷模型，本文提出了一種基于SIFs互等的裂紋面裂尖奇異區荷載等效處理方法，即將奇異區內的分布荷載等效為奇異區邊界的集中荷載，避免奇異區裂紋面受荷而無法使用W單元求解，充分發揮了W單元求解SIFs的高精高效優勢.同時，通過算例分析，給定荷載類型等效系數建議值，為裂紋面受荷求解裂尖SIFs提供了新的思路.

1 裂尖局部分布荷載的等效處理

以一含邊界裂紋的矩形板為例，其板高為2H，板寬為W，裂紋長度為a，在裂紋張口處建立整體坐標系XOY，取裂尖局部等效區邊長為2c，裂紋面受到分布壓力 σ(X)=σ0|1-X/a|λ， σ0為O點處壓力值，如圖1所示.

圖1 裂尖局部面荷載等效Fig.1 Crack tip local surface load equivalence

筆者嘗試推導整個裂紋面加載時裂尖應力場和位移場的Williams級數表達式，即裂紋面為非自由面，此時，因邊界條件較難處理，其表達式異常復雜.根據《應力強度因子手冊》[14]，裂尖局部加載和裂紋面上作用集中力情況均有高精度的參考解.基于此，現提出利用SIFs互等，即將裂尖局部分布力等效為裂紋面上一對集中力Q，如圖1所示.以λ=0為例，裂尖局部均布荷載σ0和裂紋面集中荷載Q作用時，其裂尖Ⅰ型SIFs表達式為

式中，KI，分別表示僅裂尖局部受分布壓力和裂紋面受一對集中壓力的Ⅰ型SIFs；p，p′分別表示裂尖局部加載和裂紋面作用集中力時的荷載系數，詳見文獻[14].

假定KI=則集中力Q與裂尖局部荷載σ0有如下關系：

由此，可得到裂紋面局部受均布荷載的等效公式.同理，可以推導裂紋面裂尖局部長度c上受分布荷載σ (X)的等效公式，即

式中，P表示等效荷載系數，

2 基于W單元確定裂紋面受荷的SIFs

2.1 裂紋面奇異區受荷的W單元模型

如圖1所示，綠色正方形區域為裂尖奇異區，其外圍為常規區，因W單元對奇異區尺寸不敏感，假定本文選取的奇異區尺寸等于上節中的裂尖局部尺寸.奇異區內采用W單元，外圍常規區選用四邊形8節點等參單元；常規區網格可通過ANSYS有限元軟件進行自由離散，奇異區的網格劃分規則如下：將奇異區均分為8個三角形條元，各條元離散為n個基于徑向離散比例因子α自相似的梯形子單元和一個裂尖三角形微單元，因裂尖三角形微單元極小，可忽略其剛度貢獻，則將內部n個位移場由Williams級數控制的子單元集成的梯形條元定義為W單元.W單元最外層子單元因其一條邊（圖2中的綠色邊界）為奇異區與常規區的交界，保留邊界上3個實節點，內部5個節點轉化為虛節點，稱該單元為過渡單元.奇異區內分布荷載σ (X)可通過等效處理為奇異區邊界一對集中力Q，如圖2所示.

圖2 裂尖奇異區網格劃分及分布荷載等效Fig.2 Discretization and distributed load equivalence in the crack tip singular region

以裂尖為局部坐標原點o，建立局部直角坐標系xoy，如圖2所示.以Williams級數表示裂尖奇異區位移場，并取級數的前m+ 1項：

式中，u，v分別表示裂尖局部直角坐標系下x，y軸對應的位移分量；ρ，φ分別表示裂尖局部極坐標系下的極徑和極角分量，fl,11(φ)，fl,12(φ)，fl,21(φ)，fl,22(φ)為三角函數，具體取值可參考文獻[12].

2.2 奇異區剛度方程集成

以任一W單元為例，根據四邊形8節點等參單元理論，該條元內任意第k層子單元的剛度方程可表示為

結合式(4)，將任意第k層子單元位移列陣表示為矩陣形式：

式中，T(k)表示第k層子單元轉換矩陣；?表示廣義參數列陣.

將式(5)等號兩邊同時左乘T(k)T得

并可改寫為

式中，K(k)，p(k)分 別為第k層子單元的廣義剛度矩陣和廣義荷載列陣，且當k=2 ～n時，p(k)=0.

根據式(11)可知，各層子單元剛度矩陣大小相同，且存在相同乘數?，故可將條元內第2 ～n層子單元剛度方程進行疊加：

即Kss?=0，Kss表示第2～n層子單元廣義剛度矩陣之和.

奇異區最外層的過渡單元，即W單元中第1層子單元，有3個實節點位于常規區與奇異區交界處，5個虛節點位于奇異區內部.根據節點所在區域不同，將過渡單元的位移向量分塊表示為

將過渡單元的廣義剛度方程按照虛、實節點分塊表示，并與第2～n層子單元的剛度方程進行疊加，因此可得一個W單元的剛度方程：

根據式(14)，將8個W單元剛度方程進行集成，并將奇異區邊界的裂紋面等效荷載Q疊加至相應節點可得到奇異區的整體剛度方程：

2.3 整體剛度方程的集成

將常規區所有單元的剛度方程進行集成，并以奇異區外圍節點為界，分塊表示為

計算模型分為常規區和奇異區兩個部分，即將式(15)和(16)的剛度方程進行集成，得到整個模型的剛度方程：

2.4 SIFs求解

通過引入模型的應力、位移邊界條件，對整體剛度方程式（17）進行求解，即可獲得裂尖奇異區的廣義參數列陣?.將列陣中的廣義參數a1，b1分別代入式(18)，計算獲得裂尖應力強度因子：

3 算例分析

3.1 裂紋面受均布壓力

例1如圖1所示的邊界裂紋矩形板，其板高為2H，寬為W，且有 W=2H=100cm，裂紋長度為a，板厚t取單位厚度，奇異區范圍取邊長為2c的正方形.裂紋面施加均布壓力（ λ=0）， σ0=1kN/cm2， 彈性模量E=3×104MPa ，Poisson比ν =0.3.因篇幅所限，本文先研究此方法在Ⅰ型裂紋中的適用性，Ⅱ型裂紋另做研究.

取 a /W=0.3的 模型為例，根據前期研究成果建議奇異區尺寸取 c =a/20，使用ANSYS有限元軟件對模型常規區進行網格自由離散，此時離散為146個單元，501個節點，如圖3所示.將常規區的單元、節點信息以及約束等導入自編程的FORTRAN程序，通過求解模型整體剛度方程以獲得裂尖SIFs.P=1.94+0.57(a/W)-1.51(a/W)2.本模型參數的具體取值如表1所示.

表1 模型參數（λ=0）Table 1 Model parameters (λ=0)

圖3 常規區網格離散Fig.3 Discretization in the regular region

根據《應力強度因子手冊》（P303、P304）的裂尖局部面荷載作用下荷載系數p，p′隨a/W變化曲線圖[14]，通過GetData軟件提取荷載系數p，p′代表性點的坐標，利用Excel軟件擬合即可獲得P與a/W的關系式：

由表1可以看出，隨著a/W的變化，等效荷載系數P在1.60～2.00范圍取值.鑒于此，本算例適當放大P的取值范圍，研究在P取不同值時，a/W的變化對SIFs的影響.

W單元的3個重要參數的建議值[12]為 α =0.9， m =20， n =300.將W單元計算結果與文獻[14]的解對比，為直觀體現各計算結果的區別，將SIFs無量綱化如圖4所示.

圖4 邊界裂紋面受均布壓力計算結果Fig.4 Calculation results for uniformly distributed pressure on the crack surface

均布壓力作用于邊界裂紋面時，若忽略奇異區裂紋面的荷載（P=0），計算誤差為4%～17%，且隨著a/W的增大而減小，表明隨著裂紋長度的增加，裂尖局部荷載對SIFs的影響逐漸變小.對奇異區裂紋面荷載等效處理后，現選取1%為誤差限，根據疊加法可知，當P的取值范圍為[2.0，2.1]時，任意a/W均滿足該誤差限.整體上，隨著a/W的增加，無量綱SIFs呈增大的趨勢，且增大的趨勢愈加顯著.

例2當均布壓力作用于中心裂紋上，如圖5所示，其板高為2H，寬為2W，且有W=H=100cm，矩形板厚度t為單位厚度，裂紋長度為2a，奇異區邊長2c，材料參數同本算例的邊界裂紋板.

圖5 中心裂紋面受均布荷載Fig.5 The uniformly distributed load on the central

根據對稱性，取半模型進行分析，以a/W=0.3的 模型為例，根據前期研究成果建議奇異區尺寸取c=a/20，使用ANSYS有限元軟件對模型常規區進行網格自由離散，此時離散為181個單元，602個節點，如圖6所示.將常規區的單元、節點信息以及約束等導入自編程的FORTRAN程序，通過求解模型整體剛度方程以獲得裂尖SIFs.

圖6 常規區網格離散Fig.6 Discretization in the regular region

根據《應力強度因子手冊》[14]擬合出中心裂紋的等效荷載公式P=2.03+0.11(a/W)+0.11(a/W)2.當a/W取值在[0.01，0.8]區間時，P的取值范圍為[2.03，2.18].適當放大P的取值范圍，研究P取不同值的情況下，a/W的變化對SIFs的影響.

W單元的3個重要參數的建議值[12]為 α =0.9，m=20，n=300.將W單元計算結果與文獻[14]的解對比，如圖7所示.

圖7 中心裂紋面受均布壓力計算結果Fig.7 Calculation results for uniformly distributed pressure on the central crack surface

均布壓力作用于中心裂紋面時，若忽略奇異區裂紋面的荷載（P=0），計算誤差均在11%以上，最大達到20.2%，且該誤差隨著a/W的增大而逐漸減小.同樣以1%為誤差限，根據疊加原理可以推斷出，對任意a/W，均滿足該誤差限要求的P取值范圍為[1.95, 2.05].隨著中心裂紋板的裂紋長度增加，無量綱SIFs呈遞增的趨勢，且遞增幅度逐漸變大.

綜上，對奇異區荷載進行等效處理是有必要的，在1%的誤差限下，現給出P建議取值為2.0，不僅對裂紋面受均布壓力的邊界裂紋和中心裂紋都有良好的適用性，且保證了計算精度.均布荷載和線性荷載是工程中常見的荷載類型，故下文將以裂紋面受線性壓力為例，驗證P建議取值的通用性.

3.2 裂紋面受線性壓力

算例各參數具體取值與模型網格劃分同本文的例1和例2，僅荷載形式不一樣，此時施加于裂紋面的分布壓力為σ (X)=σ0|1-X/a|λ，σ0=1kN/cm2，λ =1.結合本文的荷載等效處理方法和W單元求解算例中的SIFs.

將常規區的單元、節點信息以及約束等導入自編程的FORTRAN程序，通過求解模型整體剛度方程以獲得裂尖SIFs.本算例中，W單元重要參數建議值為α =0.9，m=20，n=300.

① 裂紋面受到線性壓力時，由于文獻[14]無荷載系數p，p′的具體取值，故無法擬合出等效荷載系數P的關系式.現嘗試P=2.0的建議取值時，研究不同a/W模型下的邊界裂紋和中心裂紋SIFs，并與文獻[14]的參考解進行比較.無量綱SIFs的計算結果如圖8所示.

由圖8可知，當P=2.0時，計算結果誤差均滿足1%的誤差限，證明在分析裂紋面受線性壓力的中心裂紋或邊界裂紋問題時，荷載等效的處理方法均有很好的通用性和高精性.當裂紋長度較小時，邊界裂紋與中心裂紋的無量綱SIFs接近，而隨著a/W的增大，無量綱SIFs呈增大的趨勢，且邊界裂紋的增長幅度大于中心裂紋.

圖8 裂紋面受線性壓力計算結果Fig.8 Calculation results for linear pressure on the crack surface

② 當P=2.0 時 ，現取邊界裂紋和中心裂紋的矩形板寬為W=100cm，研究不同的H/W情況下，a/W的變化對SIFs的影響，并與文獻[14]的參考解進行比較.無量綱SIFs結果如圖9所示.

由圖9可看出，H/W和a/W取不同值時，W單元計算結果均滿足1%的誤差限，這說明在奇異區尺寸取建議值的情況下，等效荷載系數P對裂紋位置或模型尺寸都不敏感，且適用于常見的裂紋面分布荷載類型，證明了本文的荷載等效處理方法具有很強的適用性，但由于篇幅所限，算例中只對均布荷載與線性荷載進行了論證.當H/W增大時，無量綱SIFs呈減小的趨勢；當a/W增大時，無量綱SIFs呈增大的趨勢，這表明裂紋模型是存在尺寸效應的.

圖9 無量綱SIFs隨尺寸變化的計算結果Fig.9 Calculation results for dimensionless SIFs varying with sizes

4 結 論

本文對裂紋面奇異區荷載進行等效處理，使得等效奇異區內裂紋面滿足W單元裂紋面自由的邊界條件，通過裂紋面受不同荷載的邊界裂紋和中心裂紋奇異區荷載等效算例分析，驗證了本文等效處理的正確性與通用性，得到了以下結論：

1) 確定等效奇異區尺寸后，結合算例簡化了等效荷載系數P的關系式，并建議取值為2.0.算例結果表明，本文方法與《應力強度因子手冊》對比，其SIFs誤差均小于1%；而且《應力強度因子手冊》所采用的權函數法，當模型、荷載變化時，均需重新確定權函數，相較于W單元缺乏靈活性.

2) 等效荷載系數P的取值不再是由模型尺寸、裂紋長度決定的關系式，無需通過《應力強度因子手冊》擬合，且適用于裂紋面受均布或線性壓力的邊界裂紋和中心裂紋，證明了本文對奇異區荷載進行等效處理方法具有很強的簡便性與通用性.

3) 含裂紋的矩形板，裂紋面受荷時，隨著a/W的增大，Ⅰ型無量綱SIFs均呈增大的趨勢，且增大幅度愈加明顯；隨著H/W的增大，無量綱SIFs呈減小的趨勢，這些都表明結構內部缺陷的增多會加速其破壞進程.相同模型尺寸下，邊界裂紋的無量綱SIFs整體大于中心裂紋，說明了裂紋的約束條件對SIFs的影響較大.

參考文獻（ References ）：

［1］劉鈞玉, 林皋, 范書立, 等.裂紋面受荷載作用的應力強度因子的計算[J].計算力學學報, 2008, 25（5）: 621-626.（LIU Junyu, LIN Gao, FAN Shuli, et al.The calculation of stress intensity factor including the effects of surface tractions[J].Chinese Journal of Computational Mechanics, 2008, 25（5）: 621-626.（in Chinese））

［2］LIU J Y, LIN G, LI X C, et al.Evaluation of stress intensity factors for multiple cracked circular disks under crack surface tractions with SBFEM[J].China Ocean Engineering, 2013, 27（3）: 417-426.

［3］ZHONG H, LI C L, LI H J, et al.Stress intensity factors of interfacial crack with arbitrary crack tractions[J].IOP Conference Series:Earth and Environmental Science, 2019, 304（5）: 052111.

［4］陳白斌, 李建波, 林皋.無需裂尖增強函數的擴展比例邊界有限元法[J].水利學報, 2015, 46（4）: 489-496, 504.（CHEN Baibin, LI Jianbo, LIN Gao.An extended scaled boundary finite element method without asymptotic enrichment of the crack tip[J].Journal of Hydraulic Engineering, 2015, 46（4）: 489-496, 504.（in Chinese））

［5］李亞, 易志堅, 王敏, 等.裂紋面局部均布荷載下Ⅰ型裂紋有限寬板應力強度因子[J].應用數學和力學, 2020, 41（10）:1083-1091.（LI Ya, YI Zhijian, WANG Min, et al.The stress intensity factor of a finite-width plate with a mode-Ⅰcenter crack subjected to uniform stress on the crack surface near the crack tip[J].Applied Mathematics and Mechanics, 2020, 41（10）: 1083-1091.（in Chinese））

［6］FETT T, RIZZI G.Weight functions for stress intensity factors and T-stress for oblique cracks in a halfspace[J].International Journal of Fracture, 2005, 132（1）: L9-L16.

［7］LI J, WANG X, TAN C L.Weight functions for the determination of stress intensity factor and T-stress for edge-cracked plates with built-in ends[J].International Journal of Pressure Vessels and Piping, 2004, 81（3）: 285-296.

［8］WALTERS M C, PAULINO G H, DODDS JR R H.Interaction integral procedures for 3-D curved cracks including surface tractions[J].Engineering Fracture Mechanics, 2005, 72（11）: 1635-1663.

［9］MUTHU N, MAITI S K, FALZON B G, et al.A comparison of stress intensity factors obtained through crack closure integral and other approaches using extended element-free Galerkin method[J].Computational Mechanics, 2013, 52（3）: 587-605.

［10］賈金生, 汪洋, 馮煒, 等.重力壩高壓水劈裂模擬方法與特高重力壩設計準則初步探討[J].水利學報, 2013, 44（2）:127-133.（JIA Jinsheng, WANG Yang, FENG Wei, et al.Simulation method of hydraulic fracturing and discussions on design criteria for super high gravity dams[J].Journal of Hydraulic Engineering, 2013, 44（2）: 127-133.（in Chinese））

［11］唐世斌, 劉向君, 羅江, 等.水壓誘發裂縫拉伸與剪切破裂的理論模型研究[J].巖石力學與工程學報, 2017, 36（9）:2124-2135.（TANG Shibin, LIU Xiangjun, LUO Jiang, et al.Theoretical model for tensile and shear crack initiation at the crack tip in rock subjected to hydraulic pressure[J].Chinese Journal of Rock Mechanics and Engineering, 2017, 36（9）: 2124-2135.（in Chinese））

［12］楊綠峰, 徐華, 李冉, 等.廣義參數有限元法計算應力強度因子[J].工程力學, 2009, 26（3）: 48-54.（YANG Lüfeng,XU Hua, LI Ran, et al.The finite element with generalized coefficients for stress intensity factor[J].Engineering Mechanics, 2009, 26（3）: 48-54.（in Chinese））

［13］徐華, 鄧鵬, 藍淞耀, 等.曲線裂紋裂尖SIFs等效分析的廣義參數Williams單元確定方法[J].工程力學, 2020, 37（6）:34-41.（XU Hua, DENG Peng, LAN Songyao, et al.The determination method of Williams element with generalized degrees of freedom for equivalent analysis of SIFs at the curved crack tip[J].Engineering Mechanics, 2020,37（6）: 34-41.（in Chinese））

［14］中國航空研究院.應力強度因子手冊[M].增訂版.北京: 科學出版社, 1993.（Chinese Aeronautical Establishment.Handbook of Stress Intensity Factors[M].Revised ed.Beijing: Science Press, 1993.（in Chinese））