基于連續(xù)介質損傷力學的鈦合金耐壓球殼極限強度分析

張博文，萬正權，趙 青，徐 強，張愛鋒

（1.中國船舶科學研究中心，江蘇無錫 214082；2.深海技術科學太湖實驗室，江蘇無錫 214082）

0 引 言

鈦合金比強度高、耐腐蝕性好，是優(yōu)良的輕質高強“海洋金屬”，已廣泛應用于深海載人潛水器等深海裝備的耐壓結構中。但鈦合金材料塑韌性指標相對較低，導致其對微觀缺陷更加敏感。如圖1所示，相同強度等級的鈦合金耐壓球殼[1]的破壞形式與高強鋼[2]迥然不同，鈦合金結構呈現明顯的內潰式破碎。精準預測鈦合金材料與結構的損傷、失效以及結構極限強度是開展深海載人裝備結構技術研究的關鍵，對于指導深海裝備用鈦合金性能的提升與改進具有重要意義。

圖1 耐壓球殼失效模式[2，30]Fig.1 Failure modes of pressure hulls

自上個世紀60 年代以來，大量學者針對材料的本構關系與損傷模型展開了研究工作[3-7]，其中根據損傷是否影響材料塑性階段的力學行為，損傷模型大致可以分為非耦合類和耦合類。非耦合類模型結構形式簡單，模型參數易于確定，但由于沒有考慮損傷擴展對材料力學行為的影響，導致其使用范圍受限，代表模型有Johnson-Cook（J-C）模型[8]、Wierzbicki-Bao 模型[9-10]以及L-H 模型[11-13]等。耦合類損傷模型中最具影響力的代表是Gurson-Tvergaard-Needleman（GTN）模型[14-16]和連續(xù)介質力學（Continuum Damage Mechanics，CDM）模型[17]。GTN模型材料參數關系復雜，導致材料參數的確定較為困難，大多需通過人為調試，此外經典的GTN 模型中損傷的累積并不影響材料的彈性模量，沒有真正做到彈塑性與損傷的強耦合。CDM模型根據連續(xù)介質力學與熱力學的唯象學方法研究損傷隨變形發(fā)展最后導致破壞的力學過程，通過引入損傷變量來描述含缺陷材料的力學響應及材料自身的劣化過程，并從宏觀上研究損傷的發(fā)展對材料宏觀力學性能的影響以及材料損傷演化過程和規(guī)律，以便預測材料的變形、破壞和使用壽命。Bonora 等[18]在Lemaitre 研究的基礎上對損傷演化方程進行了改進，提出了考慮非線性的損傷勢函數，并通過了實驗驗證，預測精度相比Lemaitre 模型有所提升。王治磊等[19]利用有限元軟件研究了SUS304 不銹鋼的損傷和韌性斷裂行為，研究表明Lemaitre 模型比Johnson-Cook 模型具有更高的精度。Bonora等[20]與Testa等[21]深入研究了損傷起始應變（孔洞形核應變）與應力三軸度的關系，發(fā)現隨著應力三軸度的增加孔洞形核塑性應變越低，并修正了形核應變與應力狀態(tài)的關系，更新了非線性CDM 模型，使之物理意義更加清晰。Li等[22]考慮了應力三軸度、羅德參數以及剪切應力對損傷累積的影響，進一步修正了Lemaitre 模型，并通過有限元模擬與實驗證明了修正模型對鋁合金復雜成形過程損傷行為預測的準確性。張毅等[23-24]將CDM 模型應用在聚乙烯材料的斷裂行為預測上，模擬結果與實驗結果吻合良好。

現有深潛器耐壓球殼結構極限強度分析是在板殼理論和彈塑性力學的理論框架下，隨著工作潛深與殼板厚度的增加，力學模型由薄殼結構力學模型漸變?yōu)楹駳ちW模型?，F有潛水器規(guī)范不適用于厚球殼結構的極限強度分析[25]，經大量學者廣泛研究[26-27，29]，引入了材料非線性和幾何缺陷的有限元技術可以預報厚球殼結構的極限強度。然而，為了簡化計算，眾多學者普遍選用雙線性的理想彈塑性材料本構關系模型，未充分考慮加載過程中材料力學性能硬化與軟化過程。本文將在耐壓球殼極限強度有限元分析技術的基礎上，研究Bonora連續(xù)介質損傷力學模型[28]，模擬典型球殼結構的失效模式與破壞壓力，為后續(xù)深海鈦合金耐壓結構精細化數值仿真技術研究提供基礎。

1 損傷模型理論

材料或結構失效是由于微損傷增長所致，屬于不可逆能量耗散過程。連續(xù)介質損傷力學采用損傷變量D從宏觀層面描述材料的損傷演化過程。損傷變量定義為材料中存在的不可逆缺陷的比率，D=0表示無損傷的完整材料，D=Dcr表示材料完全損傷失去承載能力，Dcr表示臨界損傷變量，一般認為等于1，但是學者普遍認為Dcr＜1。由于損傷變量D是一個難以直接進行實驗測量的量，一般只能間接測量，目前最常用的方法是利用楊氏模量的退化來測量損傷，如式（1）所示，但由于損傷對材料特性的影響是通過微結構的變化實現的，一般情況下組織不敏感量（如材料變形特性）對損傷不敏感，而組織敏感量（如強度特性）對損傷非常敏感。

式中，ED為材料損傷后的有效彈性模量，E0為無損材料的彈性模量。

由于損傷導致材料組織結構的有效承載面積有所減少，有效應力則相應會隨之增加，根據等效應變假設，作用在損傷材料上的名義應力所產生的應變等于作用在相應無損材料上有效應力所產生的應變。

式中，σeff為有效應力，σ為實際應力。根據Lemaitre實驗觀察，發(fā)現損傷和材料的彈性相關，于是把材料耗散勢函數分為不耦合的彈性部分和塑性部分，只把損傷引入到彈性勢函數（損傷耗散勢函數FD）中，塑性勢函數fp不反應損傷，如式（3）所示。

Bonora在Lemaitre準則的基礎上提出了考慮損傷演化的損傷耗散勢函數：

式中，S0是材料參數，n是材料硬化參數，α是決定損傷演化曲線形狀的損傷指數，ēp為累積等效塑性應變，Y表示應變能釋放率，

式中，σeq表示等效應力，η表示應力三軸度（靜水壓力與等效應力之比），μ為泊松比。

根據Lemaitre準則，損傷演化的動力學規(guī)律可以描述為

在比例加載的情況下，損傷變量可以表達為

式中，εˉth為單軸拉伸損傷起始時的等效塑性應變，εˉcr為單軸拉伸失效時的等效塑性應變，pth為材料損傷起始時的等效塑性應變，為了簡化計算，通常情況下認為材料損傷起始時的等效塑性應變近似等于單軸拉伸損傷起始時的等效塑性應變。從式（7）可以看出，Bonora模型通過應變能釋放率引入了應力三軸度的影響，模型中的主要參數分別為εˉth、εˉcr、Dcr和α。

2 材料試驗

試驗材料為TA31鈦合金，其主要成分如表1所示。

表1 鈦合金主要化學成分（%）Tab.1 Main chemical composition of titanium alloy（%）

試樣根據國家標準GB/T 228.1-2010 進行制備，形狀與尺寸如圖2所示。

圖2 鈦合金標準試樣Fig.2 Sample of titanium alloy

為測量損傷參數，進行兩類試驗：單向拉伸試驗和拉伸加載卸載試驗，每類試驗進行3 組重復性試驗。兩類試驗的試驗加載設備均采用UTM 5105 型電子萬能試驗機，試驗機最大拉力為100 kN，在每個試件中部懸掛電子引伸計，標距為25 mm。試驗加載均采用位移控制，加載速率設定為2 mm/min。

首先，對TA31鈦合金標準試樣進行單向拉伸至斷裂，獲得TA31鈦合金載荷位移曲線；隨后，將載荷位移曲線塑性段進行20 等分，分割點作為拉伸加載卸載試驗時的加載位置，開始時將試樣拉伸到第1 個分割點位置后進行卸載，當載荷卸載至0 后再拉伸至下個分割點位置，如此反復進行加載卸載實驗，根據每個卸載段楊氏模量衰減變化規(guī)律測定式（1）中的損傷因子。

3 損傷參數分析

在拉伸過程中材料分別經歷了線彈性階段、塑性硬化階段、頸縮階段以及失效斷裂。塑性硬化模型如Johnson-Cook（J-C）模型，本文僅考慮室溫準靜態(tài)過程，因此忽略應變速率和溫度的影響，簡化公式為

式中，A為材料的屈服極限，B和n分別為材料應變硬化模量和硬化指數。在材料進入塑性初期，材料內部孔洞缺陷形核，但所占比例較少，并不影響材料的宏觀力學性能，此時J-C 硬化模型可以較好地模擬材料硬化階段的應力應變響應。隨著材料內部孔洞缺陷擴張積累到一定程度，有效承載面積逐漸減少，有效應力逐漸增加，材料宏觀軟化行為開始凸顯，未考慮損傷因素的J-C 模型會高估材料的硬化行為。當仿真結果與試驗結果出現分離時，定義此刻為損傷起始位置，截面中心位置的等效塑性應變?yōu)棣拧h，如圖3 所示。當材料內部的微小單元出現孔洞貫穿形成宏觀裂紋時，結構承載能力喪失，在裂紋發(fā)生失穩(wěn)擴展后，載荷驟降至0。

圖3 TA31鈦合金光滑圓棒式樣載荷位移曲線Fig.3 The load-displacement curve of TA31 titanium alloy with smooth round bar sample

損傷指數α決定了損傷演化曲線的形狀，對于給定的εˉcr和Dcr，如圖4（a）所示，α越大損傷演化曲線越陡峭，斷裂時刻引伸計位移和臨界斷裂載荷越小，材料軟化速率越快，而當α減小時載荷位移曲線趨于平緩，α趨近于0 時，D也趨近于0，表明材料無軟化行為，此時，基于Bonora 模型預測的載荷位移曲線與J-C模型一致。上述規(guī)律表明，損傷指數α與材料內部孔洞缺陷的形核速率有關：當α較大，且塑性應變超過門檻值時，孔洞缺陷擴張速率逐漸增大，材料的有效承載面積加速減少，表明材料內部孔洞較多且塑韌性能較差；當α較小時，孔洞缺陷形核速率較慢，材料力學性能變化較為平緩，材料內部孔洞較少，韌塑性能較好。

圖4 載荷位移曲線Fig.4 Load-displacement curve

Dcr是臨界損傷因子，取值介于[0，1]，損傷因子的理論上限是1，表示材料完全喪失承載能力，式（2）中的有效應力趨近于無限大。Dcr=0 表示材料無損傷行為，載荷位移曲線與不考慮損傷演化行為的J-C模型一致，如圖4（b）所示，臨界損傷因子Dcr對載荷位移曲線的影響波動較損傷指數α小。

部分學者[20-21，28]認為Dcr=1 是理論情況，實際材料失效時的損傷臨界值小于1，損傷臨界值與材料內部孔洞密度有關，不同材料損傷臨界值差別較大。但隨著外部載荷的增加，材料內部孔洞缺陷動態(tài)演化，由于孔洞分布與生長的隨機性，導致很難準確地測量，不同學者利用超聲與射線技術測量得到的楊氏模量與單位體積孔洞數量之間的經驗公式會存在差別，針對相同材料得到的損傷臨界值也會不同。當Dcr＜1 時，有效應力不會出現奇異解，在加載歷程中積分點處的損傷值達到臨界值而單元被刪除時，積分點處的應力不會出現奇異，無法模擬因微缺陷演化而導致的應力集中現象。因此，本文認為在有限元模擬時Dcr=1 具有較明確的物理意義，一個單元為模型最小宏觀尺度，當單元內部等效塑性應變達到起始值εˉth時，單元內部缺陷演化加速，表現為力學性能軟化；當等效塑性應變達到臨界值εˉcr時，Dcr=1，單元失去承載能力，積分點位置有效應力發(fā)生奇異，采用單元刪除技術，將Dcr=1 的單元刪除，可以進一步模擬裂紋擴展的過程及應力應變演化規(guī)律。

εˉcr是材料失效臨界時刻的等效塑性應變，試驗中得到的塑性應變是引伸計標距范圍內的平均值，加載過程中試樣不同截面的缺陷孔洞分布不均勻，最大塑性應變發(fā)生在斷裂截面的中心位置大于0.4，最小塑性應變發(fā)生在引伸計夾持端0.02左右，根據試驗推算得到的平均斷裂塑性應變?yōu)?.14，如圖5 所示。中心位置等效塑性應變達到失效臨界值時，邊緣位置尚未達到損傷起始臨界值，而結構失效判據是敏感參量，因此，在Bonora 模型中εˉcr設定為最小截面中心位置的最大塑性應變。

圖5 應力應變曲線Fig.5 Curve of stress-strain

4 損傷參數校準

4.1 參數校準方法

根據上述分析，本文提出了一種Bonora損傷模型參數校準方法，如圖6所示。

圖6 模型參數校準流程圖Fig.6 Flow chart of model parameter calibration

首先，基于單向拉伸試驗獲得的載荷位移曲線校準J-C 硬化模型，根據仿真計算得到的載荷位移曲線與試驗值比對，確定損傷起始時刻等效塑性應變與損傷臨界時刻等效塑性應變初始值。

然后，根據加載卸載試驗獲得的真實應力應變關系換算得到損傷變量與等效塑性應變關系，認定損傷臨界值為1，根據損傷變量演化趨勢，換算得到損傷指數初始值。

最后，將Bonora 損傷模型通過VUMAT 子程序嵌入到ABAQUS 有限元軟件中，將上述初始參數代入模型，如果計算得到的載荷位移曲線和斷裂失效位置與試驗值相吻合，則認為上述參數為準確值，反之則修改參數重新計算，直到與試驗值吻合為止。

4.2 參數校準

本文材料硬化模型選用J-C 模型（如式（8）所示），基于最小二乘法并根據頸縮以前階段的單向拉伸試驗數據（平均值）進行擬合，得到模型參數如表2所示

表2 模型參數Tab.2 Model parameters

通過試驗數據處理，可以得到加載卸載狀態(tài)下真實應力應變曲線，如圖7 所示，在最后一次加載進入塑性后試件發(fā)生斷裂，通過最后一次卸載加載楊氏模量變化得到的損傷因子最大值為0.1，表明當損傷變量超過0.1時，孔洞擴張呈現加速狀態(tài)，大孔洞之間相互貫通最終導致整個截面喪失承載能力。但由于加載卸載試驗測得的損傷因子是標距范圍內的平均值，屬于鈍感參量，而在有限元計算時引入的損傷因子，表征微小單元積分點位置的材料軟化特性，當損傷因子達到損傷臨界值時，該單元有效應力發(fā)生奇異，單元喪失承載能力被刪除，屬于敏感參量。當一個單元被刪除時，周圍單元發(fā)生應力重分配，截面仍具有一定的承載能力。因此，基于加載卸載試驗得到的楊氏模量衰減規(guī)律滯后于基于有限元方法計算得到的結果，加載卸載試驗僅能反應損傷變量的演化趨勢，無法準確測定最終的損傷臨界值。

圖7 TA31鈦合金拉伸加載卸載真應力應變曲線Fig.7 True stress-strain curve of TA31 titanium alloy under tensile loading and unloading

試驗所用試樣為標準圓棒拉伸試樣，平均應力三軸度為1/3，代入式（7），經代數變換得到如下線性表達式：

加載卸載試驗測得的ε～D為平均統(tǒng)計值，但其相互間比例關系與有限元計算值一致，因此基于式（9），將試驗值代入可以推算得到損傷指數初始值α0，如圖8所示，其平均值為0.0359。

圖8 塑性應變與損傷因子關系Fig.8 Relationship between plastic strain and damage factor

將損傷指數初始值代入到Bonora模型中，進行迭代運算，確定最終模型參數，如表2 所示?；贐onora 模型仿真計算得到的載荷位移曲線與試驗值對比，如圖9所示，起裂點在最小截面中心位置，隨著載荷增加逐漸向邊緣擴展，在截面邊緣形成近似45°斜切面。從圖中可以看出，在J-C 塑性模型的基礎上，Bonora 模型可以延伸模擬材料軟化行為，有助于精準模擬結構的極限強度。

圖9 仿真計算與試驗值對比Fig.9 Comparison between simulation and test values

5 耐壓球殼極限強度分析

剛度破壞（失穩(wěn)）與強度破壞是耐壓球殼結構失效的主要形式，薄球殼結構彈性失穩(wěn)壓力與結構厚徑比（厚度t/半徑R）的平方成正比[29]。厚徑比較低時，球殼結構失效模式以剛度破壞為主；而當厚徑比較大時，球殼結構彈性失穩(wěn)壓力遠高于實際破壞壓力，材料非線性與初始幾何缺陷的影響不能忽視，強度破壞起主導作用。1996 版中國船級社《潛水系統(tǒng)和潛水器入級建造規(guī)范》給出了外壓球殼承載能力的計算方法。在此基礎上，潘彬彬[1，30]開展了鈦合金球殼結構極限承載能力模型試驗與仿真計算分析，基于大量有限元仿真計算，提出了以材料拉伸極限σb為材料輸入的鈦合金球殼結構極限承載能力計算公式，該公式經調整被接受為2018版中國船級社《潛水系統(tǒng)和潛水器入級規(guī)范》[31]（簡稱“潛規(guī)2018”）。石佳睿等[32]基于強度計算公式，提出了以材料屈服極限σs為材料輸入的考慮球殼開孔和初始缺陷的極限強度計算公式。熊志鑫等[33-34]基于切線模量理論提出了一種耐壓球殼極限強度的簡化算法，并深入研究了初始缺陷對耐壓球殼結構極限強度的影響規(guī)律。上述研究工作側重研究球殼結構幾何初撓度偏差對結構極限承載能力的影響，未在數值仿真計算模型中考慮材料強化及損傷演化帶來的材料非線性影響，導致屈服強度或抗拉強度相同，但屈強比不同的材料，通過彈塑性模型計算得到的結構極限強度是相同的。盡管基于上述公式的預報誤差在可接受范圍內，但由于缺少考慮材料的非線性行為，無法準確指導深海裝備用鈦合金材料性能的優(yōu)化與改進。因此，本文將Bonora模型引入到球殼結構的極限承載能力分析中，充分考慮材料強化及軟化行為對結構極限承載能力的影響。

本文依據Pan 等[30]研究工作中TA31 鈦合金（又名Ti80）耐壓球殼縮比模型壓力筒考核試驗（具體模型參數詳見文獻[30]），開展Bonora 模型擴展應用研究，利用有限元軟件ABAQUS 對典型球殼結構進行數值仿真計算，首先開展球殼結構屈曲模態(tài)分析，隨后根據實測缺陷幅值將1 階屈曲模態(tài)轉換為幾何缺陷引入到ABAQUS/Explicit 模塊中，編寫考慮Bonora 損傷模型的用戶自定義子程序進行非線性有限元分析。計算結果如圖10 所示，當靜水外壓載荷達到57.39 MPa 時，在模態(tài)缺陷凹坑邊緣位置損傷因子升至臨界值Dcr=1，單元被刪除，凹坑中心位置單元Mises 應力達到最大值，隨著載荷繼續(xù)增加，凹坑中心位置單元Mises應力驟降，裂紋貫穿并快速擴展致模型破壞。

圖10 基于Bonora模型的球殼結構典型位置應力應變Fig.10 Stress-strain of spherical shell based on Bonora model

為了對比分析，在ABAQUS/Riks模塊中不考慮材料非線性本構關系模型，設定材料雙線性各向同性硬化參數，材料極限強度分別設定為σs和σb，切線模量設定為0 MPa[35]。當耐壓球殼凹坑位置大范圍進入塑性時，結構喪失承載能力，達到結構極限強度，節(jié)點位移隨載荷下降而持續(xù)增加。計算結果匯總于表3，結果表明，當材料極限強度增加時，結構極限承載能力提升近7.2%，而耐壓球殼模型實際破壞壓力介于圖11所示計算區(qū)間，誤差不超過5%。強度破壞是在大潛深環(huán)境下厚球殼結構的主要失效模式，材料非線性行為的影響權重隨著厚徑比的增加而逐漸增加。計算結果表明，基于Bonora 模型的耐壓球殼結構極限承載能力計算值與文獻[30]中TA31鈦合金球殼模型試驗值較好吻合，預報誤差不到1%，預報精度有一定提升，相比于雙線性彈塑性簡化模型，Bonora 材料本構關系模型充分考慮了材料屈服強度、抗拉強度以及最大塑性應變的影響，結構破壞的物理意義更加明確。

表3 計算結果匯總Tab.3 Calculation results

圖11 基于理想彈塑性模型的球殼結構載荷位移曲線Fig.11 The load-displacement curve of spherical shell based on ideal elastic-plastic model

6 結 論

本文基于連續(xù)介質損傷力學理論框架，針對TA31 鈦合金開展了Bonora 模型適用性研究，分析了典型耐壓球殼結構的失效模式和極限承載能力，得出如下結論：

（1）本文設計并開展了鈦合金標準試樣的加載卸載試驗，獲得了材料損傷因子隨塑性應變的變化規(guī)律，但試驗測得的塑性應變與損傷因子均是標距范圍內的平均值，不同截面的損傷缺陷分布是不均勻的，最小截面中心位置的塑性應變是試驗平均值的2.8 倍，表明基于楊氏模量衰減的試驗方法測得的損傷因子演化規(guī)律滯后于有限元方法計算得到的結果，加載卸載試驗僅能反應損傷變量的演化趨勢，無法準確測定最終的損傷臨界值。

（2）未考慮損傷因素的J-C 塑性硬化模型在模擬標準試樣單向拉伸時，會高估材料的硬化行為，無法準確模擬材料頸縮后的應力應變響應。Bonora損傷模型可以較準確地模擬材料軟化行為和失效模式，模型中的主要參數分別為損傷起始塑性應變εˉth、臨界失效塑性應變εˉcr、臨界損傷因子Dcr和損傷指數α，基于本文提出的校準方法，在加載卸載試驗的基礎上，可以準確地獲得上述模型參數。損傷指數α與材料內部孔洞缺陷的形核速率有關，α越大，材料內部孔洞擴張速率越快，材料塑韌性能越差；臨界損傷因子Dcr=1 在數值仿真分析中具有明確物理意義，單元損傷因子趨近1 時，單元積分點有效應力趨于無窮大，能夠模擬因缺陷演化而導致的應力集中現象。

（3）本文首次在耐壓球殼極限強度分析中，引入了考慮損傷演化的Bonora 模型，可以準確模擬耐壓球殼結構的失效模式和破壞壓力，預報誤差不足1%。該模型在材料本構關系中考慮了屈服強度、抗拉強度以及最大塑性應變的影響，對指導材料性能的優(yōu)化與提升具有借鑒意義。但由于現階段在公開發(fā)表的論文中，基于TA31 鈦合金的球殼模型試驗結果較少，單一模型不足以定量說明Bonora 模型的預報精度，后續(xù)工作會在此基礎上針對性地開展重復試驗，進一步消除試驗分散度，提升Bonora模型的可靠性。

盡管Bonora 模型在應變能釋放率中考慮了應力三軸度的影響，但本文尚未深入研究應力狀態(tài)對損傷指數、損傷起始應變以及損傷臨界應變的影響規(guī)律，因此下一步研究重點是考慮裂紋閉合效應和應力狀態(tài)的Bonora模型適用性研究。