鎳鈦形狀記憶合金循環(huán)變形行為的晶體塑性有限元模擬

陸辰釗，周 廷，于 超，康國政

（西南交通大學(xué)力學(xué)與工程學(xué)院，成都 610031）

引言

鎳鈦形狀記憶合金（NiTi SMA）由于其固-固無擴(kuò)散熱彈性馬氏體相變而表現(xiàn)出超彈性、形狀記憶效應(yīng)和高阻尼特性，在血管內(nèi)支架、眼鏡架、阻尼器、微機(jī)電系統(tǒng)中的執(zhí)行器等領(lǐng)域得到了廣泛應(yīng)用［1-2］。由NiTi SMA制成的結(jié)構(gòu)部件和器件經(jīng)常承受著循環(huán)熱-機(jī)械載荷作用［3］。因此，需要對(duì)NiTi SMA的熱-機(jī)械循環(huán)變形行為進(jìn)行研究，以合理準(zhǔn)確地評(píng)估此類部件、器件的可靠性。

近年來，諸多學(xué)者對(duì)超彈性NiTi SMA的循環(huán)變形行為進(jìn)行了相關(guān)實(shí)驗(yàn)和理論研究。實(shí)驗(yàn)方面，Miyazaki等［4］的研究結(jié)果表明，NiTi SMA的超彈性在循環(huán)變形過程中會(huì)發(fā)生退化，主要體現(xiàn)在相變臨界應(yīng)力和耗散能降低、殘余應(yīng)變累積、相變硬化模量增加這幾個(gè)方面，且這些現(xiàn)象在一定循環(huán)周次后趨于飽和。進(jìn)一步的實(shí)驗(yàn)研究［5-11］揭示了這種超彈性退化與加載峰值、加載速率、加載路徑和環(huán)境溫度之間的相關(guān)性。此外，為了明確循環(huán)變形行為的微觀機(jī)理，Delville等［12］、Benafan等［13-14］、Hamilton等［15］、Norfleet等［16］、Hua等［17］、Xiao等［18］、Pfetzing-Micklich等［19］、Brinson等［20］通過透射電鏡等微觀觀測(cè)手段，解釋了殘余應(yīng)變累積的物理機(jī)制，即奧氏體相中的位錯(cuò)滑移、殘余馬氏體相的循環(huán)累積以及位錯(cuò)-相變交互作用。

理論模型方面，根據(jù)空間尺度的不同，描述形狀記憶合金熱-力變形行為的本構(gòu)模型大致可以分為宏觀唯象模型和晶體塑性模型兩大類［21］。宏觀唯象模型通過引入宏觀內(nèi)變量及相關(guān)演化方程來描述實(shí)驗(yàn)研究揭示的循環(huán)變形行為演化特征，此類模型不關(guān)注材料微結(jié)構(gòu)演化，具有計(jì)算量小，適用于工程應(yīng)用的優(yōu)點(diǎn)。晶體塑性模型則是建立在單晶尺度上，通過確定的晶體取向關(guān)系來度量非彈性應(yīng)變，可以方便地將不同的非彈性變形機(jī)制引入到模型中［22］，再通過尺度過渡準(zhǔn)則或有限元方法將單晶模型拓展到多晶情形。與宏觀唯象模型相比，晶體塑性模型充分地考慮了材料微結(jié)構(gòu)形態(tài)及其演化，具備堅(jiān)實(shí)的物理背景［22］，并且能夠?yàn)椴牧衔⒂^結(jié)構(gòu)設(shè)計(jì)與優(yōu)化提供理論指導(dǎo)。

然而，上述晶體塑性本構(gòu)模型僅僅考慮了馬氏體相變這一種非彈性變形機(jī)制，無法描述鎳鈦形狀記憶合金在循環(huán)變形過程中出現(xiàn)的超彈性退化現(xiàn)象。最近，Yu等［23］通過同時(shí)考慮馬氏體相變、相變誘發(fā)塑性以及二者交互作用，建立了基于晶體塑性的循環(huán)本構(gòu)模型。然而，該模型采用了顯示過渡準(zhǔn)則完成從單晶到多晶的過渡。該過渡準(zhǔn)則屬于平均場(chǎng)理論的范疇，即無法給出合金循環(huán)變形過程中應(yīng)力、應(yīng)變以及內(nèi)變量場(chǎng)在晶粒內(nèi)部的非均勻分布。因此，為了合理度量多晶體中不同晶粒之間的交互作用并反映非均勻變形特性，本文首先基于Yu等［22］建立的單晶循環(huán)本構(gòu)模型，通過編寫通用有限元軟件Abaqus中的用戶材料子程序（VUMAT）來完成對(duì)該模型的有限元移植。然后，建立多晶代表性體積單元，引入＜111＞型初始織構(gòu)。最后，通過模型計(jì)算討論了外部加載條件（應(yīng)力水平）和內(nèi)部微結(jié)構(gòu)特征（織構(gòu)強(qiáng)度）對(duì)鎳鈦形狀記憶合金超彈性循環(huán)退化的影響。

1 本構(gòu)模型及數(shù)值算法

1.1 本構(gòu)模型

本節(jié)首先對(duì)Yu等［23］提出的形狀記憶合金晶體塑性循環(huán)本構(gòu)模型進(jìn)行簡(jiǎn)要介紹。該模型建立在單晶尺度上，考慮了24個(gè)馬氏體變體和奧氏體-馬氏體界面錯(cuò)滑移系。進(jìn)一步建立了材料點(diǎn)的Helmholtz自由能，并在熱力學(xué)框架下，通過耗散能不等式推導(dǎo)出相變和位錯(cuò)滑移的熱力學(xué)驅(qū)動(dòng)力，進(jìn)而提出了內(nèi)變量（包括馬氏體體積分?jǐn)?shù)和塑性滑移量）演化方程，具體細(xì)節(jié)可參見原始文獻(xiàn)［23］。

1.1.1 應(yīng)變分解

基于小變形假設(shè)，將總應(yīng)變?chǔ)欧纸鉃閺椥詰?yīng)變?chǔ)舉、應(yīng)力誘發(fā)相變應(yīng)變?chǔ)舤r和相變誘發(fā)塑性應(yīng)變?chǔ)舤p：

1.1.2 相變驅(qū)動(dòng)力及內(nèi)變量演化

其中，T為溫度，為第α個(gè)馬氏體變體的平衡溫度，Btr為相變內(nèi)應(yīng)力張量，是第α個(gè)變體的相變阻力，可表示為：

其中d、c1是材料參數(shù)，是馬氏體相變相關(guān)的第 α個(gè)馬氏體變體導(dǎo)致的內(nèi)應(yīng)力的大小，是第α個(gè)馬氏體變體的相變模量。是位錯(cuò)密度為零狀態(tài)下的初始相變模量，Htr-p反映位錯(cuò)滑移對(duì)相變模量的影響。

其中，Yα是一個(gè)控制馬氏體相變滯回環(huán)寬度的正變量，Y0是位錯(cuò)密度為零狀態(tài)下的初始量，反映位錯(cuò)滑移對(duì)滯回環(huán)寬度的影響，m表征了馬氏體相變的粘性，由于鎳鈦形狀記憶合金中馬氏體相變粘性很弱，因此，這里m值取為50。

為反映位錯(cuò)滑移對(duì)相變模量和滯回環(huán)寬度的影響，Htr-p和的演化率可表示為：

其中，c2和c3是材料參數(shù)。

1.1.3 塑性驅(qū)動(dòng)力及內(nèi)變量演化

其中，ρα是第α個(gè)滑移系的位錯(cuò)密度，k1和k2是材料參數(shù)，k1項(xiàng)反映了位錯(cuò)的生成，項(xiàng)反映了位錯(cuò)的湮滅，b是鎳鈦形狀記憶合金的Burgers矢量的大小。

1.2 數(shù)值算法

本文通過編寫用戶材料子程序VUMAT將上述本構(gòu)關(guān)系移植到有限元軟件Abaqus中。假設(shè)第n步的應(yīng)力σn、應(yīng)變 εn以 及 內(nèi) 變 量 如 （Btr）n、（、（（Yα）n、（、（ξtr）n、、（ρα）n均為已知。在第n+1步，給出應(yīng)變?cè)隽?（Δε）n+1和時(shí)間增量 （Δt）n+1，可得應(yīng)變張量εn+1=εn+（Δε）n+1，需要計(jì)算的未知量為（Btr）n+1、 （、 （、 （Yα）n+1、 （（ξtr）n+1、（、（ρα）n+1、σn+1。在第 n+1個(gè)加載步中，使用以下的顯式數(shù)值算法：

（1）在第n+1步開始時(shí)，VUMAT子程序從用戶定義變量中讀取需要的變量。從Abaqus主程序得到應(yīng)力、應(yīng)變?cè)隽?、時(shí)間增量和溫度。

（2）計(jì)算相變驅(qū)動(dòng)力和馬氏體相變體積分?jǐn)?shù)增量。

（3）為了使馬氏體體積分?jǐn)?shù)保持在合理范圍，防止算法發(fā)散，應(yīng)該對(duì)馬氏體體積分?jǐn)?shù)增量做一些微調(diào)，隨后更新馬氏體體積分?jǐn)?shù)。

（4）計(jì)算相變誘發(fā)塑性驅(qū)動(dòng)力和塑性滑移量。

（6）根據(jù)位錯(cuò)密度更新滑動(dòng)阻力和相變內(nèi)應(yīng)力。

（7）更新相變模量和相變阻力。

（8）更新滯回環(huán)控制參數(shù)。

（9）更新相變應(yīng)變、相變誘發(fā)塑性應(yīng)變及彈性應(yīng)變。

（11）將更新的變量寫入子程序的自定義變量中，將數(shù)據(jù)返回給Abaqus主程序。

1.3 材料參數(shù)

參考Yu等［23］的工作，計(jì)算中使用的材料參數(shù)見表1。

表1 模型參數(shù)

2 多晶代表性體積單元有限元模型

計(jì)算模型如圖1所示，使用一個(gè)5×5×5的立方體來構(gòu)建一個(gè)多晶代表性體積單元。模型采用8節(jié)點(diǎn)立方體單元C3D8，模型中每一個(gè)單元代表一個(gè)奧氏體相晶粒，模型共包含125個(gè)晶粒。通過在VUMAT程序中使用旋轉(zhuǎn)矩陣對(duì)24個(gè)馬氏體變體的晶體學(xué)參數(shù)矢量進(jìn)行旋轉(zhuǎn)，以得到整體坐標(biāo)系下每個(gè)晶粒的相變切邊方向、慣習(xí)面法線方向、位錯(cuò)滑移方向和位錯(cuò)滑移面法線方向。對(duì)于常用的NiTi形狀記憶合金板材、管材、絲材和棒材，其加工成型工藝通常會(huì)導(dǎo)致試樣具有強(qiáng)＜111＞型初始織構(gòu)。因此，本文將在有限元模型中考慮這一種類型的織構(gòu)，采用晶粒取向與＜111＞方向的擺動(dòng)角度大小來描述織構(gòu)強(qiáng)度，角度越小代表織構(gòu)越強(qiáng)，即0°角時(shí)退化到單晶、90°角時(shí)退化到無織構(gòu)隨機(jī)取向多晶。本文計(jì)算中，規(guī)定底面（y=0）的位移uy=0，并規(guī)定頂面（y=1）處位移uy=˙Eyt，其中，˙Ey為應(yīng)變率。

圖1 立方體模型網(wǎng)格劃分

3 模擬結(jié)果與討論

本節(jié)將利用Yu等［23］建立的晶體塑性循環(huán)本構(gòu)模型和本文建立的多晶代表性體積單元模型對(duì)NiTi形狀記憶合金在循環(huán)變形過程中出現(xiàn)的超彈性退化現(xiàn)象進(jìn)行模擬，并討論外部加載條件（應(yīng)力水平）和內(nèi)部微結(jié)構(gòu)特征（織構(gòu)強(qiáng)度）對(duì)循環(huán)變形行為的影響。

3.1 加載水平對(duì)循環(huán)變形的影響

首先固定織構(gòu)強(qiáng)度（15°擺動(dòng)角），討論加載水平對(duì)循環(huán)變形行為的影響。峰值應(yīng)變分別取4%，7%和10%，加載10圈。第1圈和第10圈應(yīng)力-應(yīng)變曲線以及殘余應(yīng)變隨循環(huán)演化曲線分別如圖2和圖3所示。

圖2 不同加載水平下應(yīng)力應(yīng)變曲線圖

圖3 不同加載水平下殘余應(yīng)變的循環(huán)演化曲線

從圖2和圖3中可以看出，模型能夠很好的捕捉到NiTi形狀記憶合金在循環(huán)變形中出現(xiàn)的超彈性退化現(xiàn)象，即殘余應(yīng)變累積、相變臨界應(yīng)力降低、滯回環(huán)面積減小和相變模量升高。同時(shí)，隨著加載水平的增加，相同循環(huán)圈數(shù)下的殘余應(yīng)變的累積與相變開始應(yīng)力的下降更加明顯。這是由于應(yīng)變幅值越大，發(fā)生相變的馬氏體體積分?jǐn)?shù)越高，使得相變誘發(fā)塑性量越多，位錯(cuò)密度也就越大。且隨著循環(huán)的進(jìn)行，高加載峰值下較高的位錯(cuò)密度對(duì)超彈性退化的影響越顯著，從而導(dǎo)致殘余應(yīng)變的累積和相變臨界應(yīng)力的降低越明顯。也就是說，超彈性循環(huán)退化行為隨著加載水平的提升而加劇。模擬結(jié)果與現(xiàn)有的實(shí)驗(yàn)結(jié)論相吻合［8］。

3.2 織構(gòu)強(qiáng)度對(duì)循環(huán)變形的影響

接下來固定峰值應(yīng)力（800 MPa），討論織構(gòu)強(qiáng)度對(duì)循環(huán)變形的影響。擺動(dòng)角分別取為 5°、10°、15°、20°和25°，第1圈和第10圈應(yīng)力-應(yīng)變曲線以及殘余應(yīng)變隨循環(huán)演化曲線分別如圖4和圖5所示。

圖4 不同織構(gòu)強(qiáng)度應(yīng)力應(yīng)變曲線圖

圖5 不同織構(gòu)強(qiáng)度殘余應(yīng)變的循環(huán)演化曲線

圖4和圖5的結(jié)果表明，在加載的第1圈，織構(gòu)強(qiáng)度越強(qiáng)，相變模量越低，導(dǎo)致相變束應(yīng)力和逆相變開始應(yīng)力越低；在加載的第10圈，織構(gòu)強(qiáng)度越強(qiáng)，累積的塑性應(yīng)變?cè)缴伲嘧冮_始應(yīng)力隨循環(huán)減少的幅度也就越少。也就是說，超彈性循環(huán)退化行為隨著織構(gòu)強(qiáng)度的提升而減弱。

圖6 擺動(dòng)角為15度模型循環(huán)變形第10圈應(yīng)變分布云圖

圖6給出了多晶代表性體積單元在循環(huán)變形第10圈變形最大點(diǎn)的等效應(yīng)力（Mises應(yīng)力）和等效應(yīng)變分布云圖?？梢钥吹?，由于晶粒取向的不同，各個(gè)晶粒存在著較強(qiáng)的交互作用，這使得多晶代表性體積單元在循環(huán)變形過程中應(yīng)力-應(yīng)變場(chǎng)出現(xiàn)明顯的非均勻性。

圖7 擺動(dòng)角為15度模型循環(huán)變形10圈后分布云圖

圖7給出了多晶代表性體積單元在循環(huán)變形10圈后的等效殘余應(yīng)力（Mises應(yīng)力）和殘余應(yīng)變分布云圖。圖中可見，由于塑性變形的發(fā)生，多晶代表性體積單元在循環(huán)變形后存在殘余應(yīng)變，并且由于殘余應(yīng)變的非均勻性，誘發(fā)了較高的殘余應(yīng)力。

4 結(jié)論

本文基于Yu等［23］建立的考慮馬氏體相變和相變誘發(fā)塑性的鎳鈦合金晶體塑性循環(huán)本構(gòu)模型，通過對(duì)本構(gòu)方程進(jìn)行數(shù)值離散化編寫了Abaqus用戶材料子程序（VUMAT）。進(jìn)一步建立了有多晶代表性體積單元的有限元模型，考慮合金在成型過程中最常見的＜111＞織構(gòu)，討論了加載水平和織構(gòu)強(qiáng)度對(duì)循環(huán)應(yīng)力-應(yīng)變響應(yīng)的影響。可以得到以下幾點(diǎn)結(jié)論：

（1）循環(huán)應(yīng)力-應(yīng)變響應(yīng)顯著依賴于應(yīng)變加載水平，超彈性退化現(xiàn)象隨著加載水平的提高而更加顯著；

（2）循環(huán)應(yīng)力-應(yīng)變響應(yīng)顯著依賴于織構(gòu)強(qiáng)度，超彈性退化現(xiàn)象隨著織構(gòu)強(qiáng)度的提高而減弱；

（3）由于多晶體中各個(gè)晶粒晶向的差異性，多晶代表性體積單元應(yīng)力-應(yīng)變場(chǎng)呈現(xiàn)非均勻性。同時(shí)，由于塑性變形的逐漸累積和其非均勻性，合金在循環(huán)變形后存在較高的殘余應(yīng)力。