王 雷 ， 屈 平 ， 黃進浩 ， 張愛鋒 ， 萬正權

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

0 引 言

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

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

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

1 蠕變數值計算方法

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

（1）創(chuàng)建目標結構的幾何模型；

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

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

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

（5）提取蠕變前的結構應力、應變特征；

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

（7）提取蠕變應變、蠕變變形以及蠕變后的結構應力、應變特征。

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

2 耐壓結構蠕變本構方程

應用有限元數值方法計算結構的蠕變特性具有求解速度快、計算成本低的特點，在工程領域有廣泛的應用。ANSYS、ABAQUS、ADINA、MARC等大型商用有限元軟件均有成熟的蠕變計算模塊。

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

表1 ANSYS隱式蠕變本構方程Tab.1 Implicit creep constitutive equations of ANSYS

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

表2 ABAQUS蠕變本構方程Tab.2 Creep constitutive equations of ABAQUS

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

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

廣義Graham蠕變本構方程基于Graham&Walles數學模型，以指數多項式的形式對時間強化本構方程的時間項進行細化[10]。通過增加與時間相關的冪指數多項式系數，可以提升模型對蠕變應變率變化的適應性，同時考慮有限元法對變量的迭代要求。

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

3 蠕變本構模型系數與試驗驗證

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

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

將試驗模型吊入壓力筒內，向壓力筒內注水至目標壓力，保持壓力不變，進行蠕變試驗。在穩(wěn)定的壓力環(huán)境下，應用應變儀以間隔測量的方式記錄測點應變隨時間變化。模型在蠕變試驗時間為1 960 h時，壓力卸載至0，應變儀記錄了該次模型測點的縱向和周向蠕變應變。

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

表3 不同應力水平的蠕變系數和蠕變應變Tab.3 Coefficients and strain of creep for various stress levels

（1）蠕變應力水平的影響

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

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

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

（2）蠕變本構模型的影響

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

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

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

表4 不同本構模型的蠕變系數和蠕變應變Tab.4 Coefficients and strain of creep for various constitutive models

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

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

（3）蠕變時間的影響

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

表5 不同蠕變時間的蠕變系數和蠕變應變Tab.5 Coefficients and strain of creep for various creep time

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

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

表6 基于1stOpt回歸的蠕變系數和蠕變應變Tab.6 Coefficients and strain of creep for different constitutive models

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

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

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

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

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

圖6 蠕變模型測點蠕變應變曲線Fig.6 Creep curve of creep model

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

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

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

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

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

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

5 結 論

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

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

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

（3）產生蠕變變形后鈦合金耐壓結構應力重新分配，高應力區(qū)范圍有所擴大，蠕變后彈性應變和總應變均減?。幌啾扔诳v向，蠕變對環(huán)肋圓柱殼結構的徑向變形影響更大。