郝 林，陳子雄

(上海飛機設計研究院，上海 201210)

形狀記憶合金(shape memory alloys, SMAs)是一種強非線性功能材料，因其具有超彈性效應、形狀記憶效應、高阻尼能力和生物相容性等特性被廣泛應用于航空航天器結(jié)構(gòu)中[1-3]。其特殊的性質(zhì)給仿真分析帶來一定的困難，諸多學者曾嘗試對形狀記憶合金建立其本構(gòu)模型以求在仿真過程中能真實地再現(xiàn)其熱力響應力學特性[4]，如Lagoudas團隊提出的SMAs三維唯象本構(gòu)模型能較好地模擬三維形狀記憶合金結(jié)構(gòu)的熱力學特性[5-6]。然而，由于熱載荷引起相變時形狀記憶合金材料的橫向剪切剛度不是常數(shù)[7]，因此采用軟件自帶的殼體單元無法滿足描述形狀記憶合金薄殼結(jié)構(gòu)的熱力位移響應的精度要求，而大型薄殼結(jié)構(gòu)若采用三維單元進行仿真分析則代價非常大。

本文基于一階剪切變形理論(FSDT)建立了一種新的適用于形狀記憶合金薄殼結(jié)構(gòu)分析的4節(jié)點、20自由度殼單元，采用Lagoudas本構(gòu)模型模擬形狀記憶合金材料特性，運用Newton-Raphson方法迭代求解單元矩陣。最后針對給定初始相變應變條件的形狀記憶合金薄殼結(jié)構(gòu)進行分析，驗證新單元的有效性。

1 形狀記憶合金本構(gòu)方程

基于材料內(nèi)部和外部的狀態(tài)變量，用熱力學勢函數(shù)來表征材料一點的熱力學狀態(tài)?；贚agoudas本構(gòu)模型[8]，形狀記憶合金的Gibbs自由能與應力張量、溫度和馬氏體體積分數(shù)的關(guān)系如式(1)所示：

(1)

式中：G(σ,T,ξ,εt)為材料中一點的Gibbs自由能；σ,T,T0,εt,ξ分別為應力張量、溫度、參考溫度、相變應變張量和馬氏體體積分數(shù)；S,α,c,ρ,s0,u0分別為柔度張量、熱脹系數(shù)、比熱、密度、熵和參考溫度下的內(nèi)能；f(ξ)為反映相變應變演化的硬化函數(shù)。f(ξ)的具體形式如式(2)所示：

(2)

式中：Δs0,Ms,Mf,As,Af,Δu0分別為熵增、馬氏體相變起始溫度、馬氏體相變結(jié)束溫度、馬氏體逆相變起始溫度、馬氏體逆相變結(jié)束溫度以及內(nèi)能變化量。

根據(jù)Gibbs自由能方程可得到形狀記憶合金材料的應力-應變關(guān)系[5]，如式(3)所示：

(3)

式中：ε為應變。

類比于塑性流動法則，假設形狀記憶合金相變應變與馬氏體體積分數(shù)有如下關(guān)系：

(4)

式中：Λ為用以描述相變應變演化方向與幅值的張量。根據(jù)Lagoudas提出的觀點[5,9],形狀記憶合金材料相變演化的起始與結(jié)束時刻所遵循的路徑嚴格遵循式(5)所示的Clausius-Planck不等式：

(5)

式中:π為相變驅(qū)動力。

于是形狀記憶合金發(fā)生相變時的相變函數(shù)可假定為以下形式:

(6)

式中：Φ為形狀記憶合金相變函數(shù)；Y為相變熱力學力的臨界值。馬氏體體積分數(shù)演化的約束條件采用Kuhn-Tucker條件，具體表示為：

(7)

式中：Φ(σ,T,ξ)為相變函數(shù)。

所給出的形狀記憶合金本構(gòu)方程式(3)依然是一組非線性方程組，求解依舊存在很大的困難，為此采用位移增量法進行求解。

為能夠得到形狀記憶合金的切線剛度張量，對形狀記憶合金本構(gòu)關(guān)系關(guān)系式進行微分，寫作如下形式：

dσ=R∶dε+N ∶dT

(8)

式中：R為切線剛度張量；N為熱彈性張量。

以增量形式改寫式(3)并代入式(4)，得到：

dσ=S-1∶{dε-αdT-[ΔS∶σ+Δα(T-T0)+Λ]dξ}

(9)

式中：ΔS為柔度；Δα為熱脹系數(shù)。

考慮到相變函數(shù)具體形式，則有：

(10)

式中：?σΦ為相變函數(shù)對應力導數(shù)。

形狀記憶合金發(fā)生相變時有Φ=0。相變函數(shù)全微分形式為：

dΦ=?σΦ∶dσ+?TΦdT+?ξΦdξ=0

(11)

通過消除馬氏體體積分數(shù)dξ，并使用式(10)、(11)改寫式(9)即可得到切線剛度張量以及熱彈性張量具體形式，如式(12)所示：

(12)

式中：?代表馬氏體正/逆相變過程。

2 殼單元模型

形狀記憶合金殼單元采用4節(jié)點單元，每個節(jié)點有5個自由度，即4節(jié)點20自由度單元。對于每個自由度，使用線性形函數(shù)。為了使單元適應復雜邊界，采用如圖1所示的等參數(shù)變換。

圖1 殼單元等參變換示意圖

2.1 位移場假設

基于一階剪切變形假設，形狀記憶合金殼單元中任意一點位移模式具有如下形式：

(13)

式中：u(x,y,z)、v(x,y,z)、w(x,y,z)分別為該點在X,Y,Z方向的總位移；符號(?)0表示中性面的值；Ψx，Ψy分別為過點的法線繞X，Y軸轉(zhuǎn)角?；谝浑A剪切理論變形示意圖如圖2所示。

圖2 基于一階剪切理論變形示意圖

于是可通過幾何方程得到應變?yōu)椋?/p>

(14)

2.2 單元建立

單元的平衡方程可通過虛功原理建立：

(15)

式中：U為廣義位移；u為位移向量；ΨU為廣義力；f為單元體力向量；τ為面力向量；δ是變分符號；Γ為邊界；Ω為體積。

形狀記憶合金殼單元在每個節(jié)點上有5個自由度，采用形函數(shù)Ni(x,y), 單元的自由度表示為如下形式：

(16)

式中：nodes為具體點的編號。

3 非線性分析過程

由公式(15)可知，形狀記憶合金本構(gòu)方程求解為強非線性問題。綜合考慮公式(15)及(16)，應力平衡方程的弱形式可表述為：

Ψ4}TdA+F)=0

(17)

其中，

(18)

式中：σx，σxy，…為不同方向應力；elm表示單元信息；nelm為單元數(shù)；Ae為單元面積。

式(17)大括號中的表達式表示單元的內(nèi)力, 其中F包含了外面施加的面力與場力。

由于形狀記憶合金本構(gòu)方程求解是強非線性響應問題，因此本文采用Newton-Raphson迭代法以及向后差分積分方法求解。分析時，載荷以增量形式進行施加，每個增量步施加后均需求解式(19)：

-t+Δt[KU]t+Δt{ΔU}=t+Δt{ΨU}

(19)

式中：KU為單元剛度陣；Δt為時間步增量；ΔU為時間步增量Δt下廣義位移U的增量。

(20)

4 基于ABAQUS軟件UEL建立

為了在ABAQUS中實現(xiàn)所開發(fā)的形狀記憶合金殼單元，基于ABAQUS軟件用戶單元子程序(UEL)功能將開發(fā)單元數(shù)學模型按ABAQUS軟件UEL要求編寫成FORTRAN子程序。為了可視化有限元分析結(jié)果，使用了虛擬單元技術(shù)和子程序UVARM。虛擬單元和展開單元具有相同的幾何特征，可以如實反映UEL分析所得結(jié)果。

5 應用分析以及結(jié)果討論

為了檢驗所開發(fā)單元的實際應用能力，本文將其應用于一個形狀記憶合金薄殼結(jié)構(gòu)進行仿真分析。形狀記憶合金薄殼模型示意圖如圖3所示，其幾何尺寸為長140 mm、寬20 mm、厚0.55 mm，初始時室溫下的撓度約為5 mm。設薄殼沿著X方向存在大小為-0.5%的初始相變應變，形狀記憶合金初始處于完全馬氏體相，此時具有升溫鼓起的形狀記憶效應。表1、表2分別給出了材料參數(shù)以及部分模型信息。

圖3 形狀記憶合金薄殼模型

仿真過程中假定形狀記憶合金薄殼結(jié)構(gòu)與外界絕熱，結(jié)構(gòu)溫度從260 K線性上升至290 K，并保持恒定，溫度加載歷程如圖4所示。形狀記憶合金薄殼結(jié)構(gòu)隨著溫度的升高在Y軸方向上的撓度增大，如圖5所示，加載結(jié)束后，薄殼中心點變形撓度達到2.57 mm。圖6所示為形狀記憶合金薄殼中心點Y方向撓度隨溫度變化歷程，由圖可知，隨著薄殼溫度升高，盡管結(jié)構(gòu)因熱脹冷縮效應產(chǎn)生變形，但在Y方向上撓度變化并不明顯。當溫度高于262 K時，此時薄殼開始發(fā)生相變(圖7)，形狀記憶合金開始向高溫形狀轉(zhuǎn)變，Y方向上變形撓度迅速增加。當溫度達到280 K時，薄殼完全相變?yōu)閵W氏體相，隨后因熱脹效應其撓度有所增加。

表1 材料參數(shù)

圖4 加載歷程

表2 模型信息

圖5 溫度加載結(jié)束后位移分布

圖6 薄殼中心點Y方向撓度隨溫度變化歷程

圖7 薄殼馬氏體體積分數(shù)隨溫度變化歷程

圖8所示為加載過程中形狀記憶合金薄殼中性面形狀變化過程。從圖中可以看出，隨著溫度逐漸升高，薄殼在Y方向上產(chǎn)生的變形撓度也在增大，且撓度最大位置始終在薄殼中心區(qū)域。

圖8 薄殼隨溫度變形過程

6 結(jié)束語

本文提出了一種4節(jié)點20自由度的形狀記憶合金殼單元，并采用一階剪切理論描述其位移場，使該殼單元能夠有效地模擬形狀記憶合金的熱力行為?；赨EL功能將所開發(fā)的殼單元納入ABAQUS軟件實現(xiàn)仿真應用，用所開發(fā)的形狀記憶合金殼單元對懸臂薄殼結(jié)構(gòu)進行數(shù)值分析，結(jié)果表明，該單元能夠很好地預測SMA薄殼結(jié)構(gòu)的時間響應問題。