環形缺口小沖桿試樣結合內聚力模型提取斷裂韌性參數

張 宇， 劉海亭， 翁 琳， 沈 耀

(上海交通大學 材料科學與工程學院, 上海 201100)

在高溫高壓或者輻照的惡劣環境下，材料長期服役會導致其力學性能退化，容易引發事故[1].對于核電設備材料而言，輻照后組織和性能的均勻性顯著降低，如果韌性不足則會造成斷裂[2].所以從安全以及經濟效益方面考慮，評價在役設備和輻照材料的斷裂韌性至關重要.常見的斷裂韌性指標為臨界應力強度因子KIC或斷裂韌度JIC，可通過常規力學性能試驗獲得.然而，某些在役設備和輻照材料取樣條件有限，無法滿足常規試驗樣品的尺寸要求[3]，因此需要采用小樣品外推獲取.小沖桿測試作為小樣品技術的一種，最早由Baik等[4]提出，它有著近乎無損取樣、耗材少等優點[5]，可以在任何承壓設備的表面取樣，且無需對設備進行修補，彌補了常規力學性能試驗因材料不足而無法進行的缺陷.起初，小沖桿測試技術用于核反應堆殼體材料輻照后的脆化評定，大幅度節省了輻照樣品的取樣量.到20世紀90年代，小沖桿測試技術開始應用到各個工程領域，可以用來評估材料的屈服強度、抗拉強度以及斷裂韌性等力學性能，這對在役設備的剩余壽命預測有著重要意義.

基于小沖桿試驗提取材料斷裂韌性的主要方法主要有經驗關聯法[6]、應變能密度法[7]、人工神經網絡法[8]以及反向有限元法[9].經驗關聯法對一系列小沖桿試驗和標準斷裂韌性試驗的數據進行關聯，建立兩者之間的關系.但關聯公式因人而異，受試驗設備、試樣類型和尺寸的影響很大.應變能密度法最早由Petersen等[10]提出，他們將應變能密度定義為應力應變曲線下的面積，啟裂點處的值對應為臨界應變能密度.當裂尖附近最小應變能密度超過臨界值時，裂紋發生擴展[11]，但是該方法的啟裂點確定比較困難.人工神經網絡法需預先構造一個包含很多材料參數及其對應有限元模擬結果的數據庫，再利用人工神經網絡通過試驗結果從數據庫中找出相匹配的材料參數.該方法雖然準確度高，但前期準備工作量大，訓練過程復雜.反向有限元法基于某種回歸算法調整模型材料參數，使模擬曲線逐漸逼近試驗結果，當兩者之間的差別滿足設置的精度要求時，便可獲取材料參數.與前3種方法相比，反向有限元法有著更高的精度，且實現過程相對容易.Yang等[12]采用基于黃金分割搜索算法的反向有限元方法提取材料參數，很好地預測了材料的屈服應力和硬化指數.但是想要預測材料的斷裂性能，還需對算法進行改進和優化.

采用反向有限元法提取材料的斷裂韌性需要引入具體的損傷模型.目前比較常見的損傷模型有GTN模型和內聚力模型.GTN模型是最初由Gurson提出，而后由Tvergaard和Needleman發展和完善的細觀韌性斷裂模型[13].GTN模型能夠顯式考慮材料的塑性損傷，包含9個模型參數，其中需要標定的參數有3個，每個參數影響載荷位移曲線的不同階段[14].內聚力模型將材料分為無損傷連續基體和內聚力單元兩部分[15]，能夠顯式描述材料內部的裂紋擴展行為.相比于GTN模型，內聚力模型中與斷裂相關的控制性參數少，只有極限分離應力σc與斷裂能Γ0兩個參數，有利于提高反向有限元法的準確性.

常規的小沖桿樣品在加載過程中不會產生應力集中，適用于提取彈塑性本構關系.若要獲取材料的斷裂韌性，需要采用預制缺口的小沖桿試樣.常見的預制缺口有中心貫穿孔型、直線型和環形.本文采用環形缺口試樣，因其最接近于小沖桿試樣的軸對稱特征[16].為確保反向有限元法的精度，需對樣品的幾何尺寸和缺口尺寸進行優化，以提升模擬載荷位移曲線對模型參數的敏感程度.

為了從小沖桿測試中準確獲取被測材料的斷裂韌性，本文基于內聚力損傷模型，對缺口樣品的尺寸優化進行了研究，對反向有限元法的算法進行改進.首先，建立環形缺口小沖桿試驗的有限元模型，從樣品的徑厚比、裂紋深度和有無預制裂紋3個方面分析內聚力模型參數的敏感性.然后，采用基于遺傳算法的反向有限元法，驗證樣品優化設計方案的反算準確性.最后，結合隨機游走模型對收斂算法進行優化，進一步提高反向有限元法的效率和精度.

1 理論模型

1.1 內聚力模型

內聚力模型最早是由Matos等[17]提出，它們將裂紋問題視為非線性邊值問題，適用于大變形[18]過程，可以很好地描述韌性斷裂和裂紋擴展過程[19].在有限元中，內聚力模型通過基于牽引分離模型的內聚力單元實現.牽引分離模型有很多種，包括雙線性模型、指數模型、常數模型和梯形模型[20]等.本文采用雙線性模型，該模型只有σc和極限分離位移δ0兩個控制性參數，兩者都與樣品的幾何形狀及尺寸無關.考慮到雙線性模型中，Γ0(曲線下圍成的面積)可以用σc和δ0唯一確定，因此δ0可以和Γ0相互轉換.本文選取σc和Γ0作為有限元中內聚力模型參數的輸入，采用反向有限元法將其提取后，可用于標準斷裂試樣的模擬，進而獲得材料的斷裂韌性.

1.2 有限元模型

本文考慮的樣品為直徑3 mm的環形缺口小沖桿薄片，厚度d分別為0.1、0.2及0.3 mm.設計的缺口深度d′分別為各自厚度的0.4、0.5及0.6倍，缺口寬度為各自厚度的0.1倍.考慮到試樣為環形缺口的薄圓片，采用二維軸對稱的方式建模，有限元軟件為Abaqus.內聚力模型采用零厚度的內聚力單元實現.以T91鋼的J2本構模型作為有限元模型本構參數的輸入.

在裂紋尖端附近，應力梯度較大，網格要適當密集，而遠離裂紋尖端的部位，網格要適當稀疏，這樣才能同時保證計算效率和精度.因此，本文采用過渡網格策略[21]進行網格劃分.圖1為小沖桿試驗的有限元模型，其中紅線為0厚度的內聚力單元，單元類型為COHAX4，綠色網格為基體單元，單元類型為CAX4R.上夾具、下夾具和沖桿都設為解析剛體，RP-1、RP-2、RP-3分別為各自剛體的參考點.沖桿的直徑為1 mm，上卡具的內徑為1.02 mm，下卡具的內徑為1.5 mm.接觸方式為面對面接觸，滑移方式為有限滑移，摩擦因子為0.1.

圖1 小沖桿試驗二維軸對稱模型

1.3 算法模型

普通的回歸算法在處理最優化問題時，容易陷入局部最優解，為獲取全局最優解，本文采用兩種智能算法：遺傳算法和隨機游走算法.遺傳算法在搜索過程中引入了隨機因素，能夠消除對初始值的依賴性.初始生成的參數組集合稱為初始種群，種群中每個參數組被稱為個體，每個個體都有自身的適應度，本文將適應度定義為殘差的倒數.每一輪迭代根據適應度的大小以幾何概率來選擇遺傳到下一代的個體，再根據概率進行交叉和變異，產生下一代的種群.最終，末代種群產生最優個體，得到全局最優解.隨機游走模型需要一個參數組作為起點，通過隨機生成的n維單位向量，按照下式進行游走：

(1)

圖2 算法流程圖

2 內聚力模型的參數敏感性分析

反向有限元法的核心思想為基于某種回歸算法調整模型材料參數，使模擬曲線逐漸逼近試驗結果.若改變模型參數對載荷位移曲線的影響較小，則反算結果難以收斂，進而導致較大的誤差.因此，需要從樣品的幾何和缺口設計出發，對內聚力模型參數作敏感性分析.本文首先根據未開缺口的小沖桿試樣的模擬結果優選環形缺口的位置，再從試樣直徑與厚度的比例、缺口深度以及有無預制裂紋3個方面來分析內聚力模型參數的敏感性.

2.1 小沖桿試樣的缺口位置

對無缺口的試樣進行有限元模擬，研究試樣變形時最大主應力點的位置以及應力方向，以此優選缺口的位置.

3種厚度d的無缺口小沖桿試樣有限元模擬結果如圖3所示.在變形的過程中，將最大主應力Smax的位置標記出來，觀察其移動軌跡.結果表明：當加載位移到一定值后，3種厚度試樣最大主應力的位置都不再變化，始終停留在圖中黃色圓點處.同時，從圖中可以看出，最大主應力在徑向的分量較大，軸向的分量較小，因此斷裂以徑向拉開為主，符合環形缺口小沖桿試樣的斷裂規律.經過測量可知，不同厚度試樣的最大主應力的位置距離試樣中心都是0.3 mm左右，因此在距離試樣中心0.3 mm的圓周上開缺口較為合理.

圖3 無缺口小沖桿試樣最大主應力分布圖

2.2 小沖桿試樣直徑與厚度的比例

內聚力模型的參數需要在合理區間內選取.參數值過高，內聚力模型難以收斂；參數值過低，則載荷位移曲線過早下降，與實際情況不符.經過嘗試，對于T91鋼，極限分離應力取值約在800～1 200 MPa，斷裂能取值約在5～8 MPa·mm時，能夠較好滿足以上兩點.因此，本文采用的5組內聚力模型參數分別為(σc=900 MPa，Γ0=5 MPa·mm)、(σc=800 MPa，Γ0=8 MPa·mm)、(σc=1 100 MPa，Γ0=7 MPa·mm)、(σc=1 000 MPa，Γ0=8 MPa·mm)及(σc=800 MPa，Γ0=5 MPa·mm).

圖4為3種厚度試樣在不同參數下模擬得到的結果，圖中D為加載位移，F為加載反力.可以看到，其他條件相同時，厚度越小的試樣，對內聚力模型參數的敏感性越強.對于d=0.3 mm的試樣，改變參數對載荷位移曲線的影響不大，5條曲線幾乎重合.對于d=0.2 mm的試樣，5條曲線在斷裂之后的區別較為明顯，斷裂之前區分度不高，且(σc=1 100 MPa，Γ0=7 MPa·mm)和(σc=1 000 MPa，Γ0=8 MPa·mm)對應的曲線較為接近.因為在斷裂過程中網格畸變嚴重，所以載荷位移曲線后期不太光滑.對于d=0.1 mm的試樣，5條曲線在斷裂前后的區別都很明顯，層次感分明，不同的參數組對應不同的最大載荷及斷裂時刻的位移.由此可知，相比之下，d=0.1 mm的試樣對內聚力模型的參數最為敏感.

圖4 厚度對內聚力模型參數的影響

2.3 小沖桿試樣的缺口深度

考慮到0.1 mm厚度的試樣對內聚力模型參數最為敏感， 所以使用0.1 mm厚度的試樣來分析缺口深度對參數敏感性的影響.缺口深度在厚度一半的附近選取，分別為0.04、0.05和0.06 mm.采用的內聚力模型參數與3.2節相同.

圖5為3種缺口深度的試樣在不同參數下模擬得到的結果.可以看到，在其他條件相同的情況下，缺口深度越大，對內聚力模型參數的敏感性越強，但是模型的收斂性會變差.對于缺口深度0.04 mm的試樣，(σc=1 100 MPa，Γ0=7 MPa·mm)和(σc=1 000 MPa，Γ0=8 MPa·mm)得到的兩條曲線幾乎一致，說明該樣品的參數敏感程度較低.對于缺口深度為0.05 mm和0.06 mm的試樣，5條曲線差別都比較明顯.然而，后者采用 (σc=800 MPa，Γ0=8 MPa·mm)和(σc=1 100 MPa，Γ0=7 MPa·mm)得到的曲線幾乎沒有下降的階段，說明缺口深度為0.06 mm試樣的內聚力模型收斂較為困難.基于敏感性和收斂性兩個方面考慮， 擇優選取的缺口深度為0.05 mm.

圖5 厚度0.1 mm試樣缺口深度對內聚力模型參數的影響

2.4 有無預制裂紋

由于實際加工的不穩定性，環形缺口小沖桿試樣的裂紋尖端可能不夠尖銳，應力集中也可能偏離裂尖的正上方.為了能夠按照預期的路徑發生斷裂，可以考慮在缺口的正上方預制裂紋.但是，對于厚度較小的小沖桿試樣，在開缺口的基礎上預制裂紋難度很大.因此，有必要研究預制裂紋樣品的參數敏感性，以確定引入預制裂紋的必要程度.

緊湊拉伸型(CT)試樣中，預制裂紋的長度與剩余韌帶長度的比例約為1∶10.基于相同比例，小沖桿樣品的預制裂紋長度設為5 μm，如圖6(a)所示，其中紅色部分為內聚力單元.以d=0.1 mm的試樣進行模擬，內聚力模型參數與3.3節相同，模擬結果如圖6(b)所示. (σc=1 100 MPa，Γ0=7 MPa·mm)和(σc=1 000 MPa，Γ0=8 MPa·mm)得到的曲線重疊，參數組(σc=900 MPa，Γ0=5 MPa·mm ) 和(σc=800 MPa，Γ0=5 MPa·mm)的結果也很接近，載荷位移曲線的區分度不高.與無預制裂紋的結果相比，開預制裂紋后的試樣對參數的敏感性有所降低.原因可能是有預制裂紋后，應力集中更明顯，斷裂更加容易，改變參數對斷裂過程的影響變小.因此，預制裂紋對參數的敏感性的提升較小，且實際加工也較為困難，所以從經濟效益和反算精度兩個方面來看，都無預制裂紋的必要.

圖6 厚度0.1 mm試樣開預制裂紋的有限元模型及模擬結果

上述3方面的參數敏感性分析表明，通過對小沖桿樣品進行幾何尺寸的優化能夠顯著提升模擬結果對內聚力參數的敏感性.優化樣品具有如下關鍵幾何特征：厚度為0.1 mm，缺口深度為0.05 mm，無需預制疲勞裂紋.

3 反向有限元法提取模型參數

3.1 反向有限元法的準確性驗證

理論上講，可以基于被測材料的試驗載荷位移曲線獲取內聚力模型參數，通過對比模型參數的理論值，驗證反向有限元法的準確性.然而，每種材料對應的內聚力模型參數具體數值未知，也無標準可查，難以判斷反算結果是否準確.此外，某些在役設備和輻照材料的取樣條件有限，進一步增加了試驗驗證的困難程度.因此，基于上述優化樣品設計，本文以給定的 (σc=800 MPa，Γ0=5 MPa·mm，記為參數組1)和(σc=1 000 MPa，Γ0=6 MPa·mm，記為參數組2)進行有限元模擬，再將模擬所得載荷位移曲線作為逼近目標，采用反向有限元進行參數提取.通過比較反算參數與給定參數之間的差別，對反向有限元法的準確性進行驗證.

殘差函數定義為目標曲線和模擬曲線之間的面積差.基于遺傳算法的反向有限元計算結果如圖7和表1所示.可以看出，經過 1 500 次左右的迭代后，模擬和試驗的載荷位移曲線差別已經很小，斷裂之前的部分幾乎重合，斷裂后的部分稍有差異.迭代過程中，殘差函數值整體呈下降趨勢，最終反算獲得的參數與給定參數較為接近，誤差在6%以內.可以看出，基于遺傳算法的反向有限元方法能夠較為準確地提取內聚力模型參數.

圖7 遺傳算法的反算結果

表1 初始參數值與遺傳算法的反算結果

3.2 搜索算法的改進

以上結果若是繼續以遺傳算法迭代，精度會有所提高，但是耗時過長，計算成本太高.所以將遺傳算法和隨機游走算法相結合，先將迭代解逼近到全局最優解附近，再采用收斂速度較快的隨機游走算法，以達到縮短反算時間并提高精度的目的.采用種群內最多的重復個體數量來判斷迭代解是否到達全局最優解附近.遺傳算法是根據幾何概率進行選擇，個體適應度越大，被選擇的次數也就越多，遺傳到下一代就會出現重復的個體.經過多輪遺傳，個體適應度會越來越大，重復個體的數量也隨之增加.因此，當種群內最多的重復個體數量超過給定閾值后，便可認為迭代解逼近已到達全局最優解的附近.改進算法后的結果如表2和圖8所示.可以看出，改進算法的迭代次數和迭代時間均得到了大幅降低，降幅約為2/3；另一方面，模擬和試驗曲線幾乎完全重合，反算所得參數與給定參數間的誤差減小至1%以內，說明收斂精度得到了進一步提升.可以看出，改進搜索算法后的反向有限元法具有更高的計算效率和精度.

圖8 遺傳算法和隨機游走算法相結合的反算結果

表2 初始參數值與遺傳算法和隨機游走算法相結合的反算結果

4 結論

本文以環形缺口小沖桿試樣為研究對象，從試樣的徑厚比、缺口深度以及有無預制裂紋3個方面，系統地研究了內聚力參數對模擬結果的敏感性，得到了優化的缺口樣品幾何尺寸.在此基礎上，采用預先選取的內聚力模型參數進行模擬，并將載荷位移曲線作為目標結果；利用改進算法后的反向有限元法來提取模型參數，驗證了反向有限元法的準確性.主要結論如下：

(1) 對于直徑3 mm的無缺口小沖桿試樣，在加載后期，最大主應力始終停留在試樣下表面的確定位置，距離圓心的距離為0.3 mm.將此處設為環形缺口位置較為合理.

(2) 試樣直徑與厚度的比例越大，缺口深度越大，對內聚力參數的敏感性越高，反向有限元法的結果就越精確.但是缺口深度過大，內聚力模型的收斂性會下降.從敏感性和收斂性兩方面綜合考慮，缺口深度選為厚度的一半.

(3) 預制裂紋對參數敏感性的提升較小，且實際加工困難，在樣品設計時不作考慮.

(4) 基于遺傳算法的反向有限元計算表明，反算所得參數與預先選取的參數之間誤差在6%以內，模擬和目標試驗的曲線基本重合，驗證了反向有限元法的準確性.將遺傳算法和隨機游走算法相結合，對反向有限元法進行優化，誤差則減小到1%以內，模擬和目標試驗的曲線幾乎完全重合.這說明，改進搜索算法后反向有限元法的效率與精度都得到了大幅提升，也進一步驗證了樣品設計的有效性.