王 雷 ， 屈 平 ， 黃進浩 ， 張愛鋒 ， 萬正權(quán)

（1.中國船舶科學(xué)研究中心，江蘇 無錫 214082；2.深海載人裝備國家重點實驗室，江蘇 無錫 214082）

0 引 言

蠕變是材料在恒定應(yīng)力作用下，應(yīng)變隨時間的延長而增長的流變現(xiàn)象。與傳統(tǒng)的塑性變形不同，蠕變在應(yīng)力小于材料屈服極限時也會出現(xiàn)[1-2]。鈦合金在深海環(huán)境下會產(chǎn)生不同程度的蠕變變形[3]。與航空領(lǐng)域鈦合金在高溫下的中低應(yīng)力拉伸蠕變現(xiàn)象不同[4]，深海環(huán)境的鈦合金耐壓結(jié)構(gòu)蠕變屬于常溫下的高應(yīng)力壓縮蠕變。深海高壓下耐壓結(jié)構(gòu)局部高應(yīng)力區(qū)域的應(yīng)力會接近材料屈服強度，會產(chǎn)生明顯的蠕變響應(yīng)，進而蠕變應(yīng)變對鈦合金耐壓結(jié)構(gòu)的強度和穩(wěn)定性會產(chǎn)生影響。由深海環(huán)境壓縮載荷引起的、高應(yīng)力下的蠕變應(yīng)變是鈦合金耐壓結(jié)構(gòu)安全性評估與計算的重要方面。

目前國內(nèi)外對工程結(jié)構(gòu)蠕變計算方法的研究已經(jīng)取得了一定成果。Ankem等[5]基于冪次法則的蠕變本構(gòu)方程對核廢料箱保護罩開展了蠕變計算，認(rèn)為Ti-6Al-4V合金在接近屈服極限的應(yīng)力水平、持續(xù)10 000 h蠕變時間的情況下易于產(chǎn)生蠕變失效。劉學(xué)等[6]應(yīng)用修正的Graham蠕變模型對P91鋼在625℃下的不同應(yīng)力水平進行了單軸蠕變過程的數(shù)值模擬，計算結(jié)果和實驗數(shù)據(jù)吻合較好。蔡志勤等[4]對鎳基合金的航空發(fā)動機渦輪葉片在980℃下單軸蠕變的全過程作了計算。苗克軍等[7]采用θ概念投影法根據(jù)Graham蠕變模型對拉伸蠕變試驗數(shù)據(jù)進行非線性擬合，對油膜軸承襯套巴氏合金進行了蠕變計算。現(xiàn)已有較多文獻對高溫結(jié)構(gòu)在拉應(yīng)力狀態(tài)的蠕變計算方法進行了研究，但是針對常溫高應(yīng)力的鈦合金耐壓結(jié)構(gòu)壓縮蠕變尚未建立成熟的計算方法，適合深海環(huán)境下鈦合金蠕變計算模型系數(shù)尚不明確。

本文對比分析多種典型蠕變本構(gòu)模型，考慮應(yīng)力水平、本構(gòu)模型、蠕變時間等三個因素的影響，確定鈦合金材料常溫壓縮蠕變本構(gòu)方程及其系數(shù)，初步建立鈦合金耐壓結(jié)構(gòu)蠕變數(shù)值計算方法，并以耐壓結(jié)構(gòu)蠕變試驗進行驗證。開展環(huán)肋圓柱殼結(jié)構(gòu)的蠕變數(shù)值計算，給出蠕變前后模型的應(yīng)力、應(yīng)變和位移的變化情況，為鈦合金耐壓結(jié)構(gòu)蠕變特性研究提供參考。

1 蠕變數(shù)值計算方法

基于ANSYS有限元軟件的蠕變數(shù)值計算程序：

（1）創(chuàng)建目標(biāo)結(jié)構(gòu)的幾何模型；

（2）對幾何模型劃分網(wǎng)格，定義單元屬性；

（3）施加載荷與邊界條件；

（4）設(shè)置求解控制選項，進行靜力求解；

（5）提取蠕變前的結(jié)構(gòu)應(yīng)力、應(yīng)變特征；

（6）選擇蠕變本構(gòu)方程，輸入材料的蠕變系數(shù)和蠕變時間，設(shè)置加載方式和步長調(diào)節(jié)，進行蠕變求解；

（7）提取蠕變應(yīng)變、蠕變變形以及蠕變后的結(jié)構(gòu)應(yīng)力、應(yīng)變特征。

在上述的蠕變數(shù)值計算步驟中，蠕變本構(gòu)方程的選取和蠕變系數(shù)的確定是蠕變分析的關(guān)鍵要素，決定了蠕變計算的方法適用性和結(jié)果準(zhǔn)確度。本文對蠕變計算的本構(gòu)模型進行梳理，給出每個模型特征和適用階段。結(jié)合深海鈦合金耐壓結(jié)構(gòu)的蠕變特點，初步選取幾種典型的本構(gòu)模型進行比較分析，對鈦合金材料壓縮蠕變試驗數(shù)據(jù)進行非線性回歸分析，將基于不同本構(gòu)模型的蠕變數(shù)值計算結(jié)果與試驗值進行對比，最終確定適用于深海鈦合金耐壓結(jié)構(gòu)的蠕變本構(gòu)模型及其系數(shù)。

2 耐壓結(jié)構(gòu)蠕變本構(gòu)方程

應(yīng)用有限元數(shù)值方法計算結(jié)構(gòu)的蠕變特性具有求解速度快、計算成本低的特點，在工程領(lǐng)域有廣泛的應(yīng)用。ANSYS、ABAQUS、ADINA、MARC等大型商用有限元軟件均有成熟的蠕變計算模塊。

商用軟件ANSYS的蠕變分析包括隱式蠕變模型和顯式蠕變模型。顯示蠕變分析的求解采用歐拉朝前法，每個子步內(nèi)蠕變應(yīng)變率為常數(shù)，不執(zhí)行平衡迭代，計算精度較低。隱式蠕變分析計算時間較長，數(shù)據(jù)存儲量大，收斂要求較高，計算結(jié)果誤差小。對于需要很少時間步的情況適用于顯示蠕變分析，隱式蠕變分析更精確，ANSYS提供的隱式蠕變本構(gòu)模型有13個，如表1所示。

表1 ANSYS隱式蠕變本構(gòu)方程Tab.1 Implicit creep constitutive equations of ANSYS

與ANSYS類似，有限元軟件ABAQUS的蠕變計算主要包括兩個部分：獲得該結(jié)構(gòu)材料的蠕變模型參數(shù)和建立蠕變分析步。ABAQUS提供的蠕變模型有三種，分別為時間強化模型、應(yīng)變強化模型以及雙曲正弦模型，其蠕變本構(gòu)方程及特點如表2所示。

表2 ABAQUS蠕變本構(gòu)方程Tab.2 Creep constitutive equations of ABAQUS

綜合ANSYS和ABAQUS兩種有限元軟件的蠕變方程，基于表1和表2可以發(fā)現(xiàn)，雙曲正弦模型（廣義Garofalo模型）、指數(shù)模型和Norton模型認(rèn)為蠕變應(yīng)變率與時間無關(guān)，僅適用于蠕變應(yīng)變呈線性變化的穩(wěn)態(tài)蠕變階段，不能完全描述深海鈦合金耐壓結(jié)構(gòu)兩個階段的蠕變特點[8-9]。除此之外，其它形式的蠕變本構(gòu)模型根據(jù)強化準(zhǔn)則可分為兩類：時間強化和應(yīng)變強化，這兩種本構(gòu)方程均能描述初始蠕變階段的蠕變特征。相比于以應(yīng)變?yōu)橹饕卣髯兞康膽?yīng)變強化準(zhǔn)則本構(gòu)模型，時間強化準(zhǔn)則更能凸顯材料的蠕變應(yīng)變隨時間變化的特點，直觀性更好。

時間強化和修正的時間強化本構(gòu)方程均以時間和應(yīng)力為表征蠕變的主要變量，時間強化模型可以看作修正的時間強化模型與Norton模型的疊加，不僅可以表征蠕變應(yīng)變率不斷減小的初始蠕變階段，而且能夠描述應(yīng)變呈線性增加的穩(wěn)態(tài)蠕變階段。

廣義Graham蠕變本構(gòu)方程基于Graham&Walles數(shù)學(xué)模型，以指數(shù)多項式的形式對時間強化本構(gòu)方程的時間項進行細(xì)化[10]。通過增加與時間相關(guān)的冪指數(shù)多項式系數(shù)，可以提升模型對蠕變應(yīng)變率變化的適應(yīng)性，同時考慮有限元法對變量的迭代要求。

由此可見，修正的時間強化模型、時間強化模型和廣義Graham模型可以很好地描述初始蠕變階段和穩(wěn)態(tài)蠕變階段的蠕變特征，本文初步選擇這三種模型作為典型的鈦合金壓縮蠕變本構(gòu)方程。

3 蠕變本構(gòu)模型系數(shù)與試驗驗證

為了驗證經(jīng)過非線性回歸分析的蠕變本構(gòu)方程系數(shù)的正確性，開展環(huán)肋圓柱殼模型的蠕變試驗，檢測典型位置的蠕變應(yīng)變，揭示鈦合金耐壓結(jié)構(gòu)蠕變應(yīng)變的分布特征和變化規(guī)律。蠕變試驗系統(tǒng)由壓力筒、加卸載系統(tǒng)、應(yīng)變測量系統(tǒng)、數(shù)據(jù)采集系統(tǒng)和試驗?zāi)Ｐ偷冉M成。

蠕變試驗?zāi)Ｐ筒捎肨C4鈦合金材料，在環(huán)材基礎(chǔ)上精車加工而成。其結(jié)構(gòu)形式為兩端帶有封板的環(huán)肋圓柱殼結(jié)構(gòu)，由圓柱殼、肋骨、封板、加強筋等結(jié)構(gòu)組成，其中環(huán)肋圓柱殼包括試驗段和過渡段，結(jié)構(gòu)如圖1所示。模型試驗段圓柱殼的長徑比L/R=3.1，徑厚比R/t=58.3。

將試驗?zāi)Ｐ偷跞雺毫ν矁?nèi)，向壓力筒內(nèi)注水至目標(biāo)壓力，保持壓力不變，進行蠕變試驗。在穩(wěn)定的壓力環(huán)境下，應(yīng)用應(yīng)變儀以間隔測量的方式記錄測點應(yīng)變隨時間變化。模型在蠕變試驗時間為1 960 h時，壓力卸載至0，應(yīng)變儀記錄了該次模型測點的縱向和周向蠕變應(yīng)變。

選擇文獻[8]中TC4ELI材料的壓縮蠕變試驗數(shù)據(jù)進行回歸分析。蠕變本構(gòu)模型系數(shù)的影響因素主要包括應(yīng)力水平、本構(gòu)模型、蠕變時間等三個方面，下面針對這三個主要影響因素分別進行比較分析。

表3 不同應(yīng)力水平的蠕變系數(shù)和蠕變應(yīng)變Tab.3 Coefficients and strain of creep for various stress levels

（1）蠕變應(yīng)力水平的影響

選取4組不同應(yīng)力水平的蠕變數(shù)據(jù)進行擬合：第1組，0.8Rpc0.2、0.85Rpc0.2和0.9Rpc0.2三個應(yīng)力水平；第 2 組，0.8Rpc0.2、0.85Rpc0.2、0.9Rpc0.2和 1.1Rpc0.2四個應(yīng)力水平；第 3 組，0.7Rpc0.2、0.8Rpc0.2、0.85Rpc0.2和 0.9Rpc0.2四個應(yīng)力水平；第 4 組，0.7Rpc0.2、0.8Rpc0.2、0.85Rpc0.2、0.9Rpc0.2和 1.1Rpc0.2五個應(yīng)力水平。應(yīng)用Origin數(shù)據(jù)分析軟件將上述4組數(shù)據(jù)均以修正的時間強化模型為蠕變本構(gòu)方程進行蠕變系數(shù)的非線性回歸分析。以擬合得到的蠕變系數(shù)作為輸入，對環(huán)肋圓柱殼模型進行蠕變數(shù)值計算，并與跨中位置縱向蠕變應(yīng)變的試驗值進行對比，誤差如表3所示。

圖2 不同應(yīng)力水平的蠕變曲線擬合Fig.2 Fitting of creep curves for various stresses

將由0.8Rpc0.2、0.85Rpc0.2和0.9Rpc0.2三個應(yīng)力水平組成的第1組作為參照組。第2組加入了1.1Rpc0.2這一應(yīng)力水平的蠕變數(shù)據(jù)，使系數(shù)C1減小了一個數(shù)量級，該組的校正決定系數(shù)最大，擬合結(jié)果最好，同時蠕變應(yīng)變的計算值也最接近試驗值。第3組加入了未進入穩(wěn)態(tài)階段的0.7Rpc0.2應(yīng)力水平蠕變數(shù)據(jù)，相比于第1組使系數(shù)C1增加了3個數(shù)量級，系數(shù)C2減小至2，得到了較大的蠕變應(yīng)變計算值。第4組綜合了5個應(yīng)力水平的蠕變數(shù)據(jù)進行回歸分析，系數(shù)C1、C2與第1組相似，計算值與試驗值的誤差較大。由此可見，從校正決定系數(shù)與試驗值誤差兩個方面考慮，第2組的蠕變數(shù)據(jù)擬合結(jié)果最優(yōu)。

（2）蠕變本構(gòu)模型的影響

基于0.8Rpc0.2、0.85Rpc0.2、0.9Rpc0.2和1.1Rpc0.2四個應(yīng)力水平的蠕變應(yīng)變數(shù)據(jù)，應(yīng)用Origin數(shù)據(jù)分析軟件以修正的時間強化模型、時間強化模型和廣義Graham模型為蠕變本構(gòu)方程進行蠕變系數(shù)的非線性回歸分析。

時間強化模型采用兩種方法開展回歸分析：第一、基于模型直接擬合；第二、分段擬合。直接擬合方法即根據(jù)本構(gòu)方程對試驗數(shù)據(jù)直接擬合，與上一節(jié)應(yīng)力水平影響分析中修正的時間強化模型擬合方法相同。分段擬合方法是將初始蠕變階段和穩(wěn)態(tài)蠕變階段的試驗數(shù)據(jù)分開，分別以時間強化模型的第一部分和第二部分進行回歸分析，最后將擬合的系數(shù)組合。

廣義Graham模型具有7個蠕變系數(shù)，應(yīng)用Origin數(shù)據(jù)分析軟件對蠕變應(yīng)變數(shù)據(jù)進行回歸分析時，初始值對擬合結(jié)果影響較大，當(dāng)初始值設(shè)定距離理論值偏差較大時將無法得到參數(shù)結(jié)果，校正決定系數(shù)呈現(xiàn)負(fù)值，擬合失敗。本文以多組不同的初始值進行回歸分析后，得到兩組擬合結(jié)果；將兩組結(jié)果應(yīng)用ANSYS軟件分別進行數(shù)值計算后，蠕變積分不收斂，未得到蠕變計算結(jié)果。三個模型的回歸分析結(jié)果如表4所示。

表4 不同本構(gòu)模型的蠕變系數(shù)和蠕變應(yīng)變Tab.4 Coefficients and strain of creep for various constitutive models

圖3 不同本構(gòu)模型的蠕變曲線擬合Fig.3 Fitting of creep curves for various constitutive models

由表4可知，基于修正的時間強化模型回歸分析的蠕變應(yīng)變計算值更接近試驗值，基于時間強化模型直接擬合與分段擬合的計算結(jié)果均偏大。7個參數(shù)的Graham模型給蠕變參數(shù)回歸帶來了復(fù)雜的分析過程，時間項的冪指數(shù)多項式細(xì)分與系數(shù)增加雖然提升了Graham模型對蠕變應(yīng)變率變化的適應(yīng)性，卻給蠕變積分的收斂性增加了阻礙。修正的時間強化模型蠕變系數(shù)少，收斂計算快，與試驗值吻合度更好，可以較準(zhǔn)確地表征鈦合金材料的壓縮蠕變特性規(guī)律。

（3）蠕變時間的影響

選取4組不同蠕變時間的蠕變數(shù)據(jù)進行擬合：第1組，全部試驗數(shù)據(jù)；第2組，1 600 h的試驗數(shù)據(jù)；第3組，1 300 h的試驗數(shù)據(jù)；第4組，1 000 h的試驗數(shù)據(jù)。應(yīng)用Origin數(shù)據(jù)分析軟件將上述4組數(shù)據(jù)均以修正的時間強化模型為蠕變本構(gòu)方程進行蠕變系數(shù)的非線性回歸分析。以擬合得到的蠕變系數(shù)作為輸入，對環(huán)肋圓柱殼模型進行蠕變數(shù)值計算，并與跨中位置縱向蠕變應(yīng)變的試驗值進行對比，誤差如表5所示。

表5 不同蠕變時間的蠕變系數(shù)和蠕變應(yīng)變Tab.5 Coefficients and strain of creep for various creep time

由表5可知，選取全部試驗數(shù)據(jù)進行回歸分析的擬合度最高，蠕變應(yīng)變計算值與試驗值的誤差最小，擬合結(jié)果最優(yōu)。隨著蠕變時間的減少，試驗曲線與擬合曲線的擬合度逐漸降低；同時，計算值與試驗值的誤差逐漸增大。

應(yīng)用1stOpt軟件對0.8Rpc0.2、0.85Rpc0.2、0.9Rpc0.2和1.1Rpc0.2四個應(yīng)力水平的蠕變應(yīng)變數(shù)據(jù)，基于修正的時間強化模型、時間強化模型、廣義Graham模型進行蠕變系數(shù)的非線性回歸分析。將三種模型回歸得到的蠕變系數(shù)作為材料特征輸入，應(yīng)用ANSYS軟件分別進行數(shù)值計算，并與試驗值進行對比。回歸分析后計算結(jié)果如表6所示，應(yīng)用1stOpt軟件的蠕變參數(shù)回歸分析結(jié)果與Origin數(shù)據(jù)分析軟件基本相同，修正的時間強化模型計算值最優(yōu)，時間強化模型的結(jié)果偏差較大，廣義Graham模型的蠕變計算仍未收斂?；阝伜辖饓嚎s蠕變試驗數(shù)據(jù)，1stOpt軟件的非線性回歸分析進一步驗證了修正的時間強化模型的適用性。

表6 基于1stOpt回歸的蠕變系數(shù)和蠕變應(yīng)變Tab.6 Coefficients and strain of creep for different constitutive models

4 鈦合金環(huán)肋圓柱殼模型的蠕變計算

建立鈦合金環(huán)肋圓柱殼結(jié)構(gòu)的幾何模型。在此基礎(chǔ)上，采用20節(jié)點的SOLID186單元建立有限元模型。選擇映射網(wǎng)格依次對環(huán)肋圓柱殼、封板、加筋板等三部分結(jié)構(gòu)進行網(wǎng)格劃分，共劃分為87 248個單元、487 842個節(jié)點，結(jié)構(gòu)有限元模型如圖5所示。

圖5 環(huán)肋圓柱殼結(jié)構(gòu)的有限元模型Fig.5 FEM model of ring-stiffened cylindrical shell

應(yīng)用Origin數(shù)據(jù)分析軟件將0.8Rpc0.2、0.85Rpc0.2、0.9Rpc0.2和1.1Rpc0.2四個應(yīng)力水平的蠕變數(shù)據(jù)以修正的時間強化模型為蠕變本構(gòu)方程進行蠕變系數(shù)的非線性回歸分析。以擬合得到的蠕變系數(shù)作為輸入，按照第1節(jié)所示的計算步驟對環(huán)肋圓柱殼模型進行蠕變數(shù)值計算。

提取蠕變模型跨中位置的等效蠕變應(yīng)變，將其隨時間變化情況以最大值為1的無量綱化表示，如圖6所示。蠕變模型跨中位置的等效彈性應(yīng)變隨時間變化情況以最大值為1的無量綱化如圖7所示。

圖6 蠕變模型測點蠕變應(yīng)變曲線Fig.6 Creep curve of creep model

圖7 蠕變模型測點彈性應(yīng)變變化曲線Fig.7 Elastic strain curve of creep model

提取蠕變模型跨中位置的等效應(yīng)力，將其隨時間變化情況以最大值為1的無量綱化表示，如圖8所示。環(huán)肋圓柱殼模型蠕變前后應(yīng)力和應(yīng)變的極值變化情況匯總于表7。對比蠕變前后應(yīng)力變化可知，模型的等效應(yīng)力極值下降了14.2%，蠕變后耐壓結(jié)構(gòu)應(yīng)力重新分配，高應(yīng)力區(qū)范圍有所擴大。同時，蠕變后模型的彈性應(yīng)變和總應(yīng)變均減小，總應(yīng)變減小7.6%。

產(chǎn)生蠕變后環(huán)肋圓柱殼模型的最大位移變化情況如表7所示。模型總位移與縱向位移增加7.1%，位于跨中殼板處的最大徑向位移增加24.4%。由此可見，相比于縱向，蠕變對環(huán)肋圓柱殼結(jié)構(gòu)的徑向變形影響更大。

圖8 蠕變模型測點等效應(yīng)力變化曲線Fig.8 Equivalent stress curve of creep model

表7 模型蠕變前后的應(yīng)力和應(yīng)變變化Tab.7 Variations of stress and strain of creep model

表8 模型蠕變前后的最大位移變化Tab.8 Variations of the maximum displacement of creep model

5 結(jié) 論

本文對比分析多種典型蠕變本構(gòu)模型，確定鈦合金材料常溫壓縮蠕變本構(gòu)方程和系數(shù)，初步建立鈦合金耐壓結(jié)構(gòu)蠕變數(shù)值計算方法。對鈦合金環(huán)肋圓柱殼模型開展蠕變計算，給出蠕變前后模型的應(yīng)力、應(yīng)變和位移的變化情況。本文得到以下結(jié)論：

（1）通過比較分析應(yīng)力水平、本構(gòu)模型、蠕變時間等三個主要影響因素，修正的時間強化模型可以表征鈦合金耐壓結(jié)構(gòu)初始蠕變階段和穩(wěn)態(tài)階段的蠕變特性，能夠適用于深海鈦合金耐壓結(jié)構(gòu)的蠕變計算；

（2）針對應(yīng)力變化范圍較大的蠕變數(shù)據(jù)非線性回歸分析，基于時間強化模型的蠕變計算值與蠕變模型的試驗值相比結(jié)果偏大；廣義Graham模型細(xì)分的7個參數(shù)給蠕變積分計算的收斂性增加了阻礙，未形成蠕變計算結(jié)果；

（3）產(chǎn)生蠕變變形后鈦合金耐壓結(jié)構(gòu)應(yīng)力重新分配，高應(yīng)力區(qū)范圍有所擴大，蠕變后彈性應(yīng)變和總應(yīng)變均減??；相比于縱向，蠕變對環(huán)肋圓柱殼結(jié)構(gòu)的徑向變形影響更大。