本構模型對GH4169缺口件蠕變變形響應計算結果的影響

胡賢明, 嵇大偉, 莊書穎, 胡緒騰

（南京航空航天大學 能源與動力學院，南京 210016）

在持續載荷作用下，高溫構件經過一定的時間會發生不可恢復的蠕變變形，經過足夠長的時間后會發生斷裂失效，在恒定載荷下直至失效所經歷的時間稱為構件的持久壽命（或稱應力斷裂壽命）。對于航空發動機的高溫構件而言，在設計時需控制其應力水平，使其持久壽命滿足要求[1]。此外，還需對結構的蠕變變形進行限制，防止變形過大導致轉靜子碰摩和影響發動機結構裝拆[1-2]。

目前，對金屬材料的高溫蠕變變形以及持久性能方面內容已進行了大量研究，發展出了不同的預測方法[3-10]。葉文明等[3]和馬曉健[4]采用θ參數法蠕變本構模型，通過對TC11鈦合金在500 ℃下的蠕變/持久實驗進行分析，建立了基于大變形蠕變分析的持久壽命預測方法，較準確地預測缺口試件的高溫蠕變響應和持久壽命，預測結果與實驗結果的誤差在 ± 40%以內，該方法的預測精度優于一些國外文獻提出的基于關鍵點斷裂應變[5]、缺口凈截面平均有效應力[6-7]以及骨點應力[8-10]的小變形有限元分析的壽命預測方法。王延榮等[11-13]分析了不同蠕變模型的差別，發現一般的蠕變模型只能描述蠕變前兩個階段（初始階段和穩定階段）的變形。構件在其使用壽命內幾乎都處于蠕變前兩個階段，但同時也需要能夠描述第三階段（加速階段）蠕變變形的模型，以便判斷構件是否處于比較危險的蠕變加速階段；因此提出了一種歸一化參數的三階段蠕變模型，并對高溫材料渦輪盤和定向結晶材料渦輪葉片結構進行了蠕變變形及應力松弛效應計算分析，分析結果較為合理。總體來看，在復雜應力狀態下，準確預測高溫構件的蠕變變形和持久壽命仍是一個較為困難的問題。

由于實際構件在蠕變失效時已累積較大變形，因此在進行有限元分析時不能忽略因變形導致的結構尺寸和形狀的改變。GH4169合金具有良好的綜合性能，目前被廣泛用于航空發動機高溫部件中[14-15]。GH4169合金長時間工作最高溫度為650 ℃，開展該溫度下蠕變實驗和失效規律研究具有較為重要的意義。本工作基于大變形有限元方法，采用彈塑性耦合蠕變分析對GH4169光滑和缺口試件的高溫蠕變變形響應和持久壽命進行有限元計算，分析蠕變本構模型對光滑和缺口試件蠕變變形響應和持久壽命計算結果的影響。

1 材料與實驗方案

實驗材料為直接時效GH4169合金。

設計如圖1所示的光滑平板高溫拉伸試件和3種雙邊缺口平板高溫拉伸試件，缺口半徑R分別為1 mm、5 mm和20 mm，應力集中系數分別為2.8、1.5和1.1，拉伸實驗溫度為650 ℃。

光滑和缺口蠕變試件設計方案如圖2所示，蠕變實驗方案見表1，施加的載荷根據平板試件最小截面上所承受的名義應力確定，實驗溫度650 ℃。

圖 3 GH4169合金拉伸試件 （a）光滑平板；（b）R1雙邊缺口平板；（c）R5雙邊缺口平板；（d）R20雙邊缺口平板Fig. 3 Tensile specimens of GH4169 alloy （a）smooth plate；（b）R1 double notched plate；（c）R5 double notched plate；（d）R20 double notched plate

表 1 GH4169合金平板試件高溫蠕變實驗方案Table 1 High temperature creep test scheme of GH4169 alloy plate specimen

圖 1 GH4169合金平板拉伸試件 （a）光滑平板；（b）R1雙邊缺口平板；（c）R5雙邊缺口平板；（d）R20雙邊缺口平板Fig. 1 Tensile specimens of GH4169 alloy （a）smooth plate；（b）R1 double notched plate；（c）R5 double notched plate；（d）R20 double notched plate

圖 2 GH4169合金平板蠕變試件 （a）光滑平板；（b）R1雙邊缺口平板；（c）R5雙邊缺口平板；（d）R20雙邊缺口平板Fig. 2 Creep specimens of GH4169 alloy （a）smooth plate；（b）R1 double notched plate；（c）R5 double notched plate；（d）R20 double notched plate

2 GH4169合金彈塑性和蠕變本構模型

材料的彈塑性和蠕變本構模型均采用真應力、真應變數據來建立。彈性變形采用廣義Hooke定律來描述。塑性變形采用多線性各向同性硬化模型來描述。短時彈塑性變形的真應力、真應變與工程應力、工程應變之間的轉換關系如下：

式中：σ為真應力；ε為真應變；σe為 工程應力；εe為工程應變。

蠕變變形階段的真應力、真應變與工程應力和工程應變之間的轉換關系如下：

式中：σ′為 真應力；εc′為 真蠕變應變；ε0為初始加載彈塑性應變；εc為 工程蠕變應變。

本工作采取三種能描述蠕變三階段變形曲線的本構模型來建模，一種是Wilshire等[16-17]提出的θ參數模型，單軸應力下其公式為：

式中：εc為 工程蠕變應變；t為 時間；ai、bi（i= 1、2、3、4）為待定參數。第二種是Ye等[18]提出的修正蠕變本構模型，單軸應力下該模型的表達式為：

式中：k、l、m、n、p、q為待定參數。

3 基于大變形有限元分析的彈塑性和蠕變響應計算

3.1 基于大變形有限元分析的彈塑性響應計算

由實驗得到的GH4169合金平板件650 ℃拉伸數據如圖3黑色實線所示。根據式（1），對材料光滑拉伸數據進行真應力-真應變轉換，由于本構模型的擬合只能采用最大應力之前的數據，因此圖3 （a）中只給出了試件從加載到最大應力之間的真應力-真應變曲線，由彈性段可擬合得到材料的彈性模量為161.2 GPa，泊松比為0.3，并基于塑性段獲取多線性各向同性硬化塑性本構模型參數來對試件的彈塑性變形響應進行有限元計算。

根據試件的形狀在有限元軟件中進行三維建模計算，三種缺口試件的彈塑性響應大變形有限元計算結果如圖3（b）～（d）所示，由圖可知，實驗和有限元計算得到的載荷-位移曲線十分貼合，缺口試件的彈塑性變形均得到較好的預測。根據有限元計算的最大載荷可以得到缺口試件的極限強度，即名義抗拉強度，三種缺口試件計算得到的極限強度分別為：1399 MPa（R1缺口平板）、1342 MPa（R5缺口平板）以及1270 MPa（R20缺口平板），與實驗值之間的誤差分別為：+3.0%、+0.65%以及-2.8%，相對誤差均在 ± 3%以內。總之，基于大變形有限元分析方法，三種缺口平板試件的彈塑性變形和極限強度都得到較好的預測。

3.2 基于大變形有限元分析的蠕變響應計算

由實驗得到GH4169合金光滑試件650 ℃下的蠕變數據，根據式（2）對工程應力和工程蠕變應變進行轉換，轉換成真實應力和真實蠕變應變，并采用θ參數法、修正蠕變模型以及Batsoulas模型對轉換后的數據進行擬合，擬合結果如圖4所示，得到的本構模型的參數見表2～表4。將這三種蠕變本構模型分別和商用有限元軟件進行編譯鏈接，形成可執行的用戶子程序，在軟件中進行調用即可使用這三種模型完成所需的蠕變響應計算。

表 2 θ參數法模型擬合參數值Table 2 Fitting parameter values of θ-projection model

表 3 修正蠕變模型擬合參數值Table 3 Fitting parameter values of modified creep model

表 4 Batsoulas模型擬合參數值Table 4 Fitting parameter values of Batsoulas model

圖5 給出了GH4169合金光滑平板和3種缺口平板高溫蠕變試件的有限元模型。由于試件為帶厚度的板件，需建立三維單元模型，因此可采用Solid 186單元進行分網建模。根據實驗條件，在所建立模型底端的面上施加固支約束，在另一端面的上方建立一個剛性節點與端面上的節點進行耦合，并在剛性節點上施加對應于實驗的載荷，目的是為了將平板一端的力集中于一點，使有限元計算過程更加接近真實實驗過程。

圖 4 GH4169合金光滑試件蠕變數據擬合結果 （a）θ參數法；（b）修正蠕變模型；（c）Batsoulas模型Fig. 4 Creep data fitting results of GH4169 alloy smooth specimen （a）θ-projection model；（b）modified creep model；（c）Batsoulas model

圖 5 GH4169合金平板蠕變試件有限元模型 （a）光滑平板；（b）R1缺口平板；c）R5缺口平板；（d）R20缺口平板Fig. 5 Plate tensile specimens of GH4169 alloy （a）smooth plate；（b）R1 double notched plate；（c）R5 double notched plate；（d）R20 double notched plate

基于大變形有限元分析的彈塑性耦合蠕變計算步驟為：（1）計算彈塑性加載過程，設置加載結束時間為10-8h，并打開大變形計算選項；（2）計算后繼恒定載荷下的蠕變過程，打開率相關的選項，分段設置時間增量步，以10倍時間間隔為基礎，可避免蠕變計算不收斂的問題，終止分析時間的最大值為該應力條件下的蠕變壽命，分析子步可根據計算模型的幾何尺寸進行調整，其取值在100～500 h左右。

由上述建立的有限元模型，在完成彈塑性加載和蠕變保載計算過程后，通過后處理所建模型上下凸臺對應節點位置的變形位移數據，減去初始彈塑性變形位移值后，即可得到所需的蠕變變形-時間曲線數據。GH4169合金光滑平板試件高溫蠕變有限元計算結果如圖6所示，給出了光滑試件在三種不同蠕變本構模型大變形有限元計算下的蠕變響應計算結果。從圖6可以看出，采用θ參數法模型計算出的光滑試件蠕變斷裂變形值和蠕變斷裂時間與實驗值最為接近，但蠕變各階段的變形曲線與實驗曲線存在一定差距，蠕變前兩個階段的變形計算值偏大；修正蠕變模型計算出的光滑試件蠕變斷裂變形值和蠕變斷裂時間與實驗值相差最大，但其描述蠕變前兩個階段的曲線效果最好；而對于Batsoulas模型而言，其對于光滑試件的描述效果介于前兩種模型之間，且也能較好描述前兩個階段的蠕變變形。

圖 6 光滑平板大變形蠕變響應計算結果 （a）θ參數法；（b）修正蠕變模型；（c）Batsoulas模型Fig. 6 Calculation results of large deformation creep response of smooth specimen （a）θ-projection model；（b）modified creep model；（c）Batsoulas model

圖7 ～圖9給出了缺口試件在三種不同蠕變本構模型大變形有限元計算下的蠕變響應計算結果。從圖7～圖9可以看出，在大變形計算下，θ參數法模型計算出的缺口試件蠕變持久壽命與實驗值最為接近，而蠕變斷裂變形值與實驗結果相差較大，且差值隨著缺口半徑變小而增大，蠕變前兩個階段的變形計算值與實驗值相比普遍偏大；同樣，修正蠕變模型計算出的缺口試件蠕變斷裂變形值和蠕變斷裂壽命與實驗值相差最大，但其描述蠕變前兩個階段的曲線效果仍是最好；而對于Batsoulas模型而言，其對于缺口試件的變形描述效果介于前兩種模型之間，且也能較好描述前兩個階段的蠕變變形。

圖 7 R1缺口平板大變形蠕變響應預測結果 （a）θ參數法；（b）修正蠕變模型；（c）Batsoulas模型Fig. 7 Prediction results of large deformation creep response of R1 notched specimen （a）θ-projection model；（b）modified creep model；（c）Batsoulas model

圖 8 R5缺口平板大變形蠕變響應預測結果 （a）θ參數法；（b）修正蠕變模型；（c）Batsoulas模型Fig. 8 Prediction results of large deformation creep response of R5 notched specimen （a）θ-projection model；（b）modified creep model；（c）Batsoulas model

圖 9 R20缺口平板大變形蠕變響應預測結果 （a）θ參數法；（b）修正蠕變模型；（c）Batsoulas模型Fig. 9 Prediction results of large deformation creep response of R20 notched specimen （a）θ-projection model；（b）modified creep model；（c）Batsoulas model

總體來看，三種蠕變模型對于缺口平板試件在大變形有限元分析下的蠕變響應預測效果不一，θ參數法對于缺口持久壽命的預測較為準確，而蠕變變形預測效果較差；修正蠕變模型對于缺口持久壽命的預測較差，而蠕變前兩個階段的變形預測效果較好；而Batsoulas模型在大變形有限元分析下的蠕變響應預測效果介于前兩種模型之間。同時，三種模型預測的蠕變斷裂變形值與實驗值相比均偏大，且差值隨著缺口半徑變小而增大。

3.3 基于大變形有限元分析的蠕變/持久壽命預測分析

上文已采用大變形有限元分析方法計算了光滑和缺口試件在高溫下的整體變形響應曲線，隨計算時間的增加，試件變形量將逐漸增大，其內部危險部位變形逐漸累積，當整體變形響應因其局部蠕變變形累積到第三階段時急劇增大，有限元計算過程將無法收斂而停止，此時可認為試件發生持久斷裂破壞[3]。因此，上文計算獲得的蠕變斷裂時間即為基于大變形有限元方法預測的蠕變/持久壽命。

為了更加直觀地對比三種蠕變本構模型對試件蠕變/持久壽命的預測效果，采用壽命分散帶圖形進行描述，如圖10所示。對于θ參數法模型而言，采用基于大變形有限元分析的蠕變/持久壽命預測方法可以較準確地預測光滑試件的蠕變壽命及三種缺口試件的持久壽命，預測誤差基本在 ±2倍分散帶以內，個別數據點在+2.5倍以內；修正蠕變模型對于蠕變/持久壽命預測效果較差，預測誤差基本在+7倍分散帶以內，壽命預測偏大；Batsoulas模型蠕變/持久壽命預測結果同樣偏大，預測誤差基本在+5倍分散帶以內。該結論進一步說明，在使用大變形有限元方法預測高溫構件蠕變響應時，采用θ參數法模型可較為準確地預測構件的持久壽命，而采用Batsoulas模型和修正蠕變模型可較為準確地預測蠕變前兩個階段的變形。

圖 10 大變形蠕變/持久壽命預測結果 （a）θ參數法；（b）修正蠕變模型；（c）Batsoulas模型Fig. 10 Prediction results of large deformation creep rupture life （a）θ-projection model；（b）modified creep model；（c）Batsoulas model

4 結論

（1）基于大變形有限元分析的方法能較好地預測彈塑性變形響應，三種缺口平板試件的彈塑性變形和極限強度都得到良好的預測，極限強度預測誤差均在 ± 3%以內。

（2）三種蠕變本構模型對于缺口平板試件在大變形有限元分析下的蠕變響應預測效果不一，θ參數法對于缺口持久壽命的預測較為準確，而蠕變變形預測效果較差；修正蠕變模型對于缺口持久壽命的預測較差，而蠕變前兩個階段的變形預測效果較好；而Batsoulas模型在大變形有限元分析下的蠕變響應預測效果介于前兩種模型之間。

（3）基于大變形有限元分析的蠕變/持久壽命預測方法，θ參數法模型可以較為準確地預測光滑試件的蠕變壽命及三種缺口試件的持久壽命，預測誤差基本在 ± 2倍分散帶以內，個別數據點在+2.5倍以內；修正蠕變模型對于蠕變/持久壽命預測效果較差，預測誤差基本在+7倍分散帶以內，壽命預測偏大；Batsoulas模型蠕變/持久壽命預測結果同樣偏大，預測誤差基本在+5倍分散帶以內。