船用灰鑄鐵材料裂紋擴展的近場動力學研究

李 雙，呂海寧，黃小華，楊建民

（1.上海交通大學海洋工程國家重點實驗室，上海 200240；2.上海交通大學三亞崖州灣深海科技研究院，海南 三亞 572000；3.廣西大學防災與結構安全教育部重點實驗室，南寧 530004）

0 引 言

在船舶與海洋結構物的設計制造及使用過程中，不可避免地會出現裂紋的萌生和擴展問題，這對結構的承載能力和安全性能起著明顯的削弱作用[1]。綜合性能優(yōu)良的高強度灰鑄鐵作為廣泛應用于船舶工程各種復雜機械設備構件（內燃機的機架、汽缸蓋、燃油泵、船體水泵、高壓箱體、基座、推進器、軸承、高速齒輪、飛輪、防護外殼等）的主要金屬材料[2]，一旦產生裂紋擴展和脆性斷裂問題，將導致船體發(fā)動機無法正常運行、船舶動力裝置受損，帶來安全隱患并影響船舶的日常工作。因此，對船用灰鑄鐵構件進行裂紋擴展行為模擬，研究試件從裂紋萌生、擴展、交匯直至斷裂的整個服役過程，對掌握裂紋擴展規(guī)律、預防構件脆性斷裂、評估船體灰鑄鐵結構的完整性和安全性具有重要的指導意義。

目前，用來模擬裂紋擴展的數值方法中，有限元法、內聚力法和擴展有限元法等得到廣泛應用。余天堂等[3]采用改進的擴展有限元法對張開型和閉合型裂紋的裂紋擴展進行了模擬；唐志波等[4]基于船體表面微裂縫在應力和海水腐蝕耦合作用下的擴展機理，根據Charles-Hillig 模型的有限元方法模擬了裂縫的動態(tài)擴展過程；陳家旺等[5]基于有限元分析軟件ANSYS 研究了帶有初始表面裂紋的潛艇典型結構在不同下潛深度狀態(tài)下的裂紋擴展情況；錢怡等[6]提出一種基于累積損傷理論和裂紋尖端彈塑性應力場的數值方法來預測裂紋的擴展行為；王敏等[7]用擴展有限元方法來分析含裂紋板的裂紋擴展，采用求解Reissner 板問題的混合插值單元，避免了剛度矩陣無窮大導致的單元剛度矩陣變?yōu)槠娈惖膯栴}；苗婷等[8]采用內聚力方法對X65管線鋼進行了斷裂韌度評定及裂紋擴展行為研究；李亞政等[9]基于損傷力學和擴展有限元方法對艦船蒸汽輪機葉片進行了裂紋擴展研究；李卓等[10]采用內聚力和擴展有限元方法對含雙邊預制裂紋的灰鑄鐵金屬材料進行了裂紋擴展的試驗和模擬研究；Chang等[11]采用擴展有限元方法對含單邊裂紋的灰鑄鐵材料進行了裂紋擴展和斷裂分析。

上述基于連續(xù)介質力學的數值模擬方法，例如有限元雖然在結構變形分析中取得較大成功，但在解決破壞問題時，由于運動方程涉及位移對空間的微分，在不連續(xù)處微分型運動方程沒有定義，出現奇異性，必須人為地對不連續(xù)處進行專門的處理；而且上述有限元法、內聚力法和擴展有限元法等在處理不連續(xù)問題時，都不可避免地存在網格依賴性、裂紋的萌生和擴展需要引入額外的斷裂準則、涉及多裂紋相互作用時存在較大局限性等問題，這無疑增加了數值模擬的復雜性。目前，針對船用灰鑄鐵的裂紋擴展研究還較少，灰鑄鐵作為一種脆性金屬材料，具有良好的鑄造、減振、耐磨、切削加工性能和較低的缺口敏感性，被廣泛應用于結構較為復雜的船海工程中，因此，有必要研究新的數值計算方法來對船用灰鑄鐵的裂紋擴展行為做進一步分析。

2000 年，新興的近場動力學（Peridynamics，PD）理論[12]的提出彌補了這些不足。它采用位移的空間積分形式完全重構了傳統連續(xù)介質力學理論，這種積分型方程在不連續(xù)處仍有定義，避免了奇異性問題，且不存在網格依賴性。PD 的斷裂準則通過一個記憶函數嵌入本構方程中，能夠模擬裂紋自發(fā)的萌生和擴展，實現了對連續(xù)介質和非連續(xù)介質力學行為的統一描述，是一種特別適合模擬材料破壞的力學理論[13]。Silling[14]采用PD 方法研究了沖擊荷載作用下Kalthoff-Winkler 韌性鋼試件的裂紋擴展和動態(tài)斷裂行為；黃小華等[13]采用雙參數PD 方法模擬了單軸拉伸荷載下含平行雙裂紋薄板的裂紋擴展行為；Ha等[15]采用PD方法研究了脆性材料PMMA在單軸拉伸荷載下裂紋動態(tài)擴展和分叉的力學行為；Ran 等[16]采用改進的PD 方法進行了金屬鋼材的的裂紋擴展分析和由氫原子導致的氫脆斷裂問題研究。

上述非局部PD 模型中，雖然可以應用黃小華等[13]和Silling[14]的PD 方法進行船用灰鑄鐵金屬脆性材料的裂紋擴展分析，但是其微觀參數鍵剛度系數都為常數，不能反映物質點鄰域內長程力隨距離增大而衰減的空間分布規(guī)律，影響計算精度。而且，上述PD 模型存在的“邊界效應”問題也需要進一步研究，PD“邊界效應”問題[17]，即物體內部的中心物質點具有完整的鄰域，而邊界處物質點的鄰域卻小于物體內部物質點的鄰域，造成邊界處物質點連接鍵的缺失，在計算時會產生較大誤差。

因此本文基于近場動力學理論，通過校正PD的“邊界效應”，考慮PD力非局部長程力的影響并推導出一種新的損傷斷裂準則來研究船用灰鑄鐵金屬脆性材料的裂紋擴展行為，預測灰鑄鐵構件的裂紋萌生、擴展和斷裂的整個完整過程，采用近場動力學數值方法研究船用灰鑄鐵材料的裂紋擴展行為。本研究可為預防船用構件脆性斷裂和評估船體灰鑄鐵結構的安全性提供參考。

1 鍵基PD概述

近場動力學理論[12]是使用節(jié)點力的積分來代替運動方程中的空間導數，物體被離散為一系列帶質量的物質點，各物質點在有限的范圍δ內彼此相互作用。假設空間域R任一物質點X與δ范圍內的其他物質點X'∈R:‖‖X'-X≤δ在某一時刻t之間存在相互作用力f，則根據牛頓第二定律可得PD 方法的運動方程為

式中，ρ為材料密度，u為物質點的位移，b為外荷載密度，HX為X的近場范圍，dVX'表示物質點X’處的體積微元。近場范圍外物質點與X的作用力f為零。

物質點間的本構力函數f包含了材料的本構信息，也稱為本構力函數，在微彈性材料[12]中，f=f(η,ξ)可用一個標量函數ω(η,ξ)表示,

式中，η為相對位移η=u'-u；ξ為相對位置ξ=X'-X；ω(η,ξ)為微觀應變能密度，類似于分子動力學方法中的對勢函數，用于描述近場范圍內兩物質點間相互作用的強弱[12]。對于線彈性材料，ω(η,ξ)可表達為

式中：c(ξ,δ)為微觀模量函數，在不同維度下的表達式如表1所示，其中E為彈性模量；s為物質點對的相對伸長率，s=(|ξ+η|- |ξ|)/ |ξ|。由式（2）～（3），對點力函數f可表示為

表1 鍵基PD微模量函數表達式Tab.1 Micro-modulus functions of bond-based PD

式中，μ(t,ξ)是一個歷史依賴的標量函數，表示物質點對的破壞情況，

式中，s0為臨界伸長率。在PD理論中，通過定義標量函數φ(X,t)來反映物質點X的損傷，

式中，0≤φ(X,t)≤1。φ(X,t)=0 表示材料未損傷；φ(X,t)=1 表示該點不再與其他點相互作用，即完全損傷。

2 改進鍵基PD

2.1 考慮PD非局部長程力修正

Huang 等[18]和章青等[19]指出傳統鍵基PD 模型中的微觀模量c(ξ,δ)為常數，不能反映長程力隨距離增大而衰減的空間分布規(guī)律，影響計算精度。本文引入長程力函數g(ξ,δ) =(1-(ξ/δ)4)，則PD 微觀模量可重新表達為

方程中的本構力函數f(η,ξ)可改進為

式中，c(0,δ)為空間維數相關的材料微觀模量；g(ξ,δ)為反映非局部長程力特性的空間分布函數，其強化了近距離物質點間的相互作用，弱化了遠距離物質點的相互作用[19]。PD應變能密度表達式為

通過PD應變能密度與連續(xù)介質力學應變能密度相等的方法，可推導得到在不同維度下材料微觀模量c(0,δ)的表達式，如表2所示。

表2 考慮非局部長程力的鍵基PD微模量函數表達式Tab.2 Micro-modulus functions of bond-based PD with long-range force correction

2.2“邊界效應”修正

在傳統鍵基PD 和2.1節(jié)考慮長程力特性的PD 模型中，定義PD 相互作用的參數是通過假設PD 點具有完整的鄰域來計算的。然而，物體內部的中心物質點具有完整的鄰域，而物體邊界處物質點的鄰域卻小于物體內部物質點的鄰域，邊界點或近邊界點比體內點有更少的鍵。因此，根據應變能密度等效的方法求得微觀模量c(ξ,δ)，在物體內部和邊界處均取相同的常數，弱化了邊界處的微觀模量值，導致計算結果不準確，即存在所謂的“邊界效應”問題。本文在2.1節(jié)的PD修正模型基礎上，借鑒Ganzenmuller等[20]的思想，重新推導物質點在體內和邊界處的PD鍵剛度系數，對邊界效應進行修正。

在計算邊界區(qū)域的應變能密度時，僅考慮邊界區(qū)域近場范圍內物質點的貢獻，則當發(fā)生空間或平面內各向均勻變形時，PD應變能密度為

式中，Hi為HX的離散形式，ηij為物質點Xi和Xj的相對位移，ξij為物質點Xi和Xj在參考構型中的相對位置。為了方便起見，物質點X和X'被替換為物質點Xi和Xj進行表示。因此，在2.1節(jié)考慮長程力修正的PD模型基礎上，PD應變能密度可進一步表示為

令PD 應變能密度與連續(xù)介質力學應變能密度相等，可得長程力效應修正和邊界效應修正后的PD微模量函數表達式，如表3所示。

表3 考慮邊界效應和長程力修正后的鍵基PD微模量函數表達式Tab.3 Micro-modulus functions of bond-based PD with long-range force and surface effect correction

在傳統PD 模型中，所有相互作用的粒子微模量常數都是相同的，因此滿足fij=-fji。在當前改進的PD模型中，fij=-fji必然成立，但采用離散求和方式計算時，Hi和Hj、ci(0,δ)和cj(0,δ)不一定相同，為保證動量守恒，取

2.3 損傷斷裂準則

PD 中材料的損傷是通過消除物質點之間的相互作用進行考慮的，假設當兩個物質點鍵的拉伸s超過其臨界值s0時，損傷發(fā)生。在基于鍵基PD 模擬斷裂的方法上，鍵的臨界拉伸值s0可以由臨界能量釋放率G0求得，臨界能量釋放率表示產生單位面積裂紋所需的能量，對于脆性材料是一個可采用試驗測量的值。

由于考慮了長程力效應修正和邊界效應修正，傳統基于鍵基PD 的臨界拉伸準則便不再適用，本文根據修正過的臨界微觀應變能密度，推導出一種新的求解物質點鍵臨界拉伸表達式。鍵基PD模型臨界微觀應變能密度ω0(ξ,η)可表示為

在三維狀態(tài)下，臨界能量釋放率G0的積分形式可通過ω0(ξ,η)表達為

二維狀態(tài)下，臨界能量釋放率為

其中令板厚h=1，邊界效應修正后的微模量cij替代常數微模量c進行求解，根據方程（14）和方程（15）可得不同維度下損傷斷裂準則，即鍵的拉伸臨界值sij0，如表4所示。

表4 改進鍵基PD的鍵拉伸臨界值Tab.4 Critical stretch of the improved PD model

需要注意的是，物體邊界或近邊界點不具有完整鄰域，這將導致邊界或近邊界點的臨界拉伸小于體內點的臨界拉伸值。在數值模擬中，邊界可能會首先發(fā)生破壞，為了避免這種情況而更好地模擬試件真實的裂紋擴展行為，在邊界處設置虛擬物質點來填充邊界處不完整的鄰域空間，而這些虛擬點并不參與相關的計算，使邊界處物質點和內部物質的鍵臨界伸長量保持一致。

3 算例驗證

3.1 微模量系數和應變能密度驗證

為了驗證本文提出的微觀模量計算方法的合理性，以及改進PD模型與原PD模型的應變能密度的計算精確性，考慮一邊長為0.5 m、厚度為0.01 m 的正方形薄板，其彈性模量E=1 GPa，泊松比ν=1/3，取物質點間距Δ=0.01 m 進行均勻離散，采用鄰域半徑δ=3.015Δ進行計算。采用改進的鍵基PD 模型計算物質點（y=0 和y=0.25）的鍵剛度系數，如圖1（a）所示；計算的薄板均勻變形時物質點（y=0）的應變能密度誤差如圖1（b）所示。

圖1 微模量系數和應變能密度相對誤差對比Fig.1 Comparison of micro-modulus coefficients and strain energy density relative errors

由圖1（a）可知，傳統鍵基PD模型和長程力修正的鍵基PD模型的微模量系數均為常數；而改進的鍵基PD模型計算的邊界處微模量系數值不再與體內物質點的微模量系數值相同，邊界處物質點的鍵剛度系數的數值得到強化，合理有效地削弱了“邊界效應”。由圖1（b）可知，傳統鍵基PD 模型在邊界處物質點的應變能密度與理論解（彈性應變能）有較大偏差；考慮長程力修正的鍵基PD模型減小了這一相對誤差；而本文改進的鍵基PD 模型在此基礎上進一步減小了誤差，邊界處物質點應變能密度的相對誤差控制在0.03%左右，內部誤差降至零。

3.2 非破壞誤差分析

為了驗證提出的改進鍵基PD模型的計算精度，通過數值解與解析解進行對比分析。采用長1 m、寬0.5 m的長方形薄板，彈性模量E=200 GPa，密度矢量ρ=7850 kg/m3，泊松比ν=1/3，幾何參數和材料參數與黃小華等[13]采用的誤差分析算例相同，薄板被離散為100×50 個物質點，物質點間距Δ=0.1 mm，鄰域半徑δ=3.015Δ，沿矩形短邊施加的單軸應力荷載為p=200 MPa，應力邊界條件的施加方式是從薄板邊界區(qū)域上以體力密度的形式施加在板內物質點上。

根據彈性力學基本理論，二維平面x和y位移解析解計算公式為ux(x,y=0)=px/E，uy(x=0,y)=-νpy/E。x和y位移的數值解與解析解ux、uy的相對誤差計算公式為ex=(-ux)/ux，ey=(-uy)/uy，其中，ex、ey分別表示x和y位移的相對誤差。

圖2（a）和2（b）分別展示了傳統鍵基PD 模型和改進PD 模型計算的位移相對誤差云圖。圖2（c）展示了改進鍵基PD模型相對誤差絕對值 |ex|≥5%的物質點數量和比例。

圖2 相對誤差位移云Fig.2 Contours of relative errors

如圖2 所示，傳統鍵基PD 模型和改進鍵基PD 模型的最大相對誤差分別為25.4%和6.26%，內部相對誤差分別為16.3%和0.41%。采用改進PD模型得到的具有相對誤差絕對值 |ex|≥5%的物質點個數和比例僅為24個和0.5%，表明了改進鍵基PD方法的有效性和準確性。

4 船用灰鑄鐵裂紋擴展模擬

4.1 含單邊缺口灰鑄鐵裂紋擴展

為了證明改進的鍵基PD方法對船用灰鑄鐵的損傷預測和裂紋擴展模擬的能力，本文采用與李卓等[10]相同的單邊缺口灰鑄鐵試件進行模擬計算。試件長寬厚度分別為0.2 m、0.03 m 和0.007 m，試件的左右兩端0.05 m 范圍內是機械夾持部位，如圖3 所示，分別在單邊預制長度為a=0.005 和0.01m、寬度為b=0.002 m 的裂紋。試件的材料參數分別為：彈性模量E=210 Gpa，泊松比ν=1/3，單軸抗拉強度ft=100 MPa，質量密度ρ=7150 kg/m3，斷裂伸長率φ=0.3%。在PD 模擬中，物質點均勻離散尺寸為Δ=0.5 mm，物質點鄰域半徑δ=3.015Δ。在試件的兩端均勻地施加位移控制荷載，每一步的荷載增量為Δu=1.0×10-8m，采用動態(tài)松弛法進行求解，時間步為Δt=1.0。

圖3 不同初始長度單邊裂紋試件的裂紋擴展示意圖Fig.3 Crack propagation of single-edge specimens with different pre-crack lengths

圖4 展示了采用改進鍵基PD 模型模擬的含5 mm 單邊預制裂紋灰鑄鐵試件在t=0、1720、1900 和2090 時刻的裂紋擴展示意圖。圖5 展示了采用改進鍵基PD 模型模擬的含10 mm 單邊預制裂紋灰鑄鐵試件在t=0、1350、1600和2010時刻的裂紋擴展示意圖。由圖4和圖5可知，兩種試件的破壞模式相似，裂紋首先從預切口尖端開始擴展，然后在單軸拉伸荷載下沿垂直于加載方向擴展，最后到達自由邊界，兩試件發(fā)生破壞。但是，含0.01 m 裂紋灰鑄鐵試件的起裂時刻和貫通時刻早于0.005 m 灰鑄鐵試件。可知，在一定條件下，預制裂紋的長度對裂紋起裂和貫通時間有較大影響，預制裂紋越大，裂紋萌生越早，試件越容易斷裂。

圖4 5 mm單邊預制裂紋灰鑄鐵試件的裂紋擴展示意圖Fig.4 Crack propagation of gray cast iron specimens with 5 mm unilateral notch

圖5 10 mm單邊預制裂紋灰鑄鐵試件的裂紋擴展示意圖Fig.5 Crack propagation of gray cast iron specimens with 10 mm unilateral notch

圖6 展示了兩種不同預制單裂紋灰鑄鐵試件的試驗結果和模擬結果對比情況，模擬結果和試驗結果基本保持一致。由試驗結果可知，裂紋雖然也沿垂直于加載方向擴展，但擴展呈鋸齒狀，而數值模擬的裂紋路徑更光滑。可能的原因是，試驗所用的灰鑄鐵材料并非均質材料，從拉斷之后的斷面可看出試件內部存在很多細微的孔洞[10]，而PD 模型將灰鑄鐵試件假設為理想均質材料，不考慮內部微觀孔洞對裂紋擴展路徑的影響，因此，裂紋形狀對比呈現一定差異性。

圖6 兩種不同預制單裂紋灰鑄鐵試件的試驗[10]和模擬破壞形態(tài)對比示意圖（擴大十倍）Fig.6 Numerical results of gray cast iron and experiments of two different single-edge specimens(magnified ten times)

4.2 含雙邊缺口灰鑄鐵裂紋擴展

采用和4.1 節(jié)相同的灰鑄鐵材料屬性，預制六種雙邊裂紋形式：雙邊裂紋長度各為5 mm，沿垂直于加載方向的對稱軸進行對稱分布，雙裂紋尖端距離分別為5 mm、10 mm和15 mm，為方便描述，分別記為5-5 試件、5-10 試件和5-15 試件；雙邊裂紋長度各為10 mm，沿垂直于加載方向的對稱軸進行對稱分布，雙裂紋尖端距離分別為5 mm、10 mm 和15 mm，分別記為10-5 試件、10-10 試件和10-15 試件。采用改進的PD 方法模擬的六種灰鑄鐵試件的最終破壞形態(tài)如圖7所示。如圖7（a）～（c）所示，5-5試件在荷載作用下，初始裂紋張開并向前擴展之后，在兩條裂紋之間的區(qū)域出現了交匯串接現象，5-10試件和5-15試件的雙裂紋各自獨立擴展，擴展速率相同，且并無串接趨勢；三種試件的破壞模式分別和10-5 試件、10-10 試件和10-15 試件的破壞模式相似，當雙裂紋尖端距離為5 mm 時，裂紋串接，裂紋尖端距離是10 mm 和15 mm 時，裂紋獨自沿水平方向擴展，直至破壞。這也證實了影響裂紋串接的并不是初始裂紋長度，而是裂尖縱向距離的大小，裂紋初始長度改變的是裂紋萌發(fā)和試件斷裂的時間。

圖7 不同預制雙裂紋灰鑄鐵試件的數值模擬裂紋擴展示意圖Fig.7 Numerical simulation of crack propagation of different double-edge gray cast iron specimens

圖8 給出了六種不同預制雙裂紋灰鑄鐵試件的裂紋擴展試驗結果，與圖7 數值模擬結果對比可知，試驗結果和模擬結果具有很好的一致性。由于材料的非均勻性問題，導致兩條裂紋的擴展速率并不相同，在同一時刻的裂紋擴展長度不同，但最終破壞模式一致，都是沿垂直于加載方向擴展，直至裂紋尖端到達試件邊界處，試件破壞。

圖8 不同預制雙裂紋灰鑄鐵試件的裂紋擴展試驗結果[10]Fig.8 Experimental results of crack propagation of different double-edge gray cast iron specimens

實例分析結果證明了改進的鍵基PD 模型能夠很好地模擬灰鑄鐵試件從連續(xù)體到非連續(xù)體的整個破壞過程，驗證了該方法用于船用灰鑄鐵類材料的脆性裂紋擴展和破壞分析是可行的。此外，由于采用空間積分方程求解來模擬裂紋擴展，在近場動力學模擬過程中不存在常規(guī)連續(xù)介質方法模擬時的網格二次剖分或裂尖應力奇異性等問題，隨著荷載的施加，初始裂紋自然起裂、擴展和貫通。

5 結 論

本文基于近場動力學理論，提出一種改進鍵基PD 方法來模擬船用灰鑄鐵材料的裂紋擴展問題。該模型考慮了PD非局部長程力的影響，修正了PD邊界效應，推導了對應的鍵剛度系數和斷裂準則表達式，研究了含單邊裂紋和雙邊裂紋的船用灰鑄鐵構件在單軸拉伸荷載下的裂紋擴展和斷裂失效問題，得出以下結論：

（1）通過薄板的變形計算了改進PD 模型的鍵剛度值、應變能密度和物質點位移相對誤差，結果顯示邊界處物質點的鍵剛度系數的數值得到強化，合理有效地削弱了PD 邊界效應；邊界處物質點應變能密度的相對誤差控制在0.03%左右，內部誤差降為零；物質點位移相對誤差大大減小，有效地提高了PD模型的計算精度。

（2）改進的PD 數值模擬的船用灰鑄鐵構件裂紋開裂位置、擴展方向和最終破壞路徑與試驗結果吻合較好，證明了本文提出模型的適用性和準確性。由于試驗試件內部存在微孔洞的原因，模擬和試驗的部分裂紋形狀雖存在一定差異，但對宏觀裂紋擴展路徑和最終破壞形態(tài)影響較小。

（3）預制裂紋的裂尖縱向間距是影響裂紋串接的重要因素，預制裂紋的初始長度對裂紋萌發(fā)和試件斷裂的時間有較大影響。