形狀記憶合金結構的拓撲優化設計方法

中圖分類號： TH122 文獻標志碼：A

Topology optimization design method of shape memory alloy structure

CHEN Shitong， DING Xiaohong， XIONG Min （SchoolofMechanical Engineering， UniversityofShanghai for Scienceand TechnologyShanghai 2ooo93，China）

Abstract： Based on the superelastic properties exhibited by shape memory alloy， a cosine type constitutive model and a solid isotropic material penalty model were combined to describe the characteristics of shape memory alloy. A material interpolation function was established and a variable density topology optimization method was used to achieve structural topology optimization design considering the nonlinearity of shape memory alloy. The proposed method was applied to design the typical two-dimensional and three-dimensional structures. Static analysis and martensitic transformation evolution analysis were used to analysis the structure obtained by the proposed optimization method.

The results show that the structure optimized by the proposed method has beter mechanical properties compared to the traditional one. In static analysis， the maximum internal stress decreases by 19.57% and the maximum displacement decreases by 31.03% . In the analysis of martensitic transformation evolution， the maximum volume fraction of martensite at the end of loading decreases by 20.27% . The proposed optimization design method results in the even overal stress distribution within the structure， less prone to local stress concentration， and less deformation caused by martensitic transformation.

Keywords： shape memory aloy; topology optimization; structure design; material nonlinearity; Martensitic transformation evolution

形狀記憶合金作為一種新型的智能材料，在航空航天、生物醫療等領域都取得了廣泛的應用。與不銹鋼、鑄鐵等傳統的線彈性金屬材料相比，形狀記憶合金的顯著特性是在外部溫度和力學刺激下發生奧氏體相和馬氏體相的無擴散固態相變。這種相變往往伴隨著高達 6%～8% 的可恢復彈性形變以及材料屬性的變化，這類非線性材料特性通常被稱為形狀記憶效應或超彈性[1]。鑒于這一獨特的材料性質，形狀記憶合金被廣泛應用于多個領域。王奇等[2將形狀記憶合金絲組合布置在機翼后緣上下關節鉸鏈處，通過加熱觸發形狀記憶效應，帶動鉸鏈關節部分實現機翼后緣結構的偏轉，設計了機翼后緣的智能可控變形結構。戴剋戎等[3使用鎳鈦基形狀記憶合金制作了骨折環抱器，這種環抱器利用形狀記憶效應，使環抱臂變形并抱緊骨斷端兩側，實現了對斷骨的固定。Boroschek等[4利用形狀記憶合金受力學刺激發生超彈性變形時會消耗大量能量的滯回特性，制作了形狀記憶合金阻尼器，將形狀記憶合金絲嵌人到剛性框架結構中，并引入振動臺進行測試。結果表明，相比原始框架結構，帶阻尼結構的峰值相對位移減少了 30%～60% ，峰值加速度減少了近 58% 。

結構優化設計指通過改變結構的幾何尺寸、形狀或排列以及材料分布等方式，改善結構在載荷工況下的負載表現，從而提升結構性能的一種設計方法[5]。常見的結構優化設計方法有尺寸優化、形狀優化、拓撲優化等。結構優化設計通過結合實際工況，可以顯著改善結構在工作狀態下的承載效能，進一步提升零件及材料在實際應用狀態下的性能。目前，該類優化設計方法已經廣泛應用于航空航天、生物醫療、機械制造領域。Al-Tamimi等基于變密度法，對長骨骨板在3種不同受力情況下的固定結構進行拓撲優化設計。王小華等基于模態分析結果對離心式壓縮機底座進行拓撲優化設計，在保證強度、剛度要求的前提下使底座減重 30% ，取得了良好的輕量化設計效果。肖密等[提出一種基于各向異性材料插值的點陣散熱結構跨尺度拓撲優化方法，通過考慮最小化平均溫度、最小化溫差、最小化最高溫度3類典型熱學性能指標，開展跨尺度拓撲優化設計，結果表明，優化結構具有良好的熱學性能。但由于形狀記憶合金的復雜非線性特性，目前，針對這種材料的拓撲優化方法研究仍處于初級階段。

在考慮形狀記憶合金超彈性材料特性的結構優化設計中，需要選取合適的本構模型來描述材料特性。相關學者對形狀記憶合金的材料特性進行了大量研究，結合材料實驗及數學理論計算等多種方法，建立了多種描述其在不同工況下的宏觀、細觀以及介觀本構模型。Thamburaja等[]基于晶體塑性滑移理論建立了單晶模型，分析了形狀記憶合金相變過程中馬氏體重定向現象。張揚[]通過實驗觀察，發現鎳鈦合金在馬氏體相變時具有層狀微結構，并基于此現象構建了馬氏體相變的鎳鈦合金單晶體和多晶體的細觀力學本構模型。Auricchio 等[1]基于廣義塑性原理，引人Drucker-Prager加載函數，得到馬氏體體積分數變化規律及增量型應力應變關系，進而建立了馬氏體體積分數與外力、溫度之間的本構關系，用于計算超彈性和單程記憶效應。Kan等[12]綜合考慮形狀記憶合金相變過程中的馬氏體取向過程和殘余應變的累積，以及兩相中的塑性變形過程和相變誘發的塑性變形現象，構建了一個較為全面的熱力學耦合本構模型，并使用該模型模擬了形狀記憶合金彈簧在循環固態制冷過程中的彈熱性能演化特征。在眾多本構模型中，宏觀唯象模型因其計算過程簡單、參數較少、與實際應用環境下的載荷工況結合較緊密的特點，在實際工程設計中取得了較多的應用[13]。

本文基于描述鎳鈦基形狀記憶合金非線性材料特性的余弦型本構模型，結合變密度拓撲優化設計方法，通過改進優化迭代過程中的材料插值模型，實現考慮形狀記憶合金材料非線性的結構拓撲優化設計，并將其應用于典型的二維和三維結構設計中。

1 設計方法

1.1 形狀記憶合金本構模型

為了研究外部力學刺激和溫度變化對形狀記憶合金材料特性的影響，王振清等[14]通過示差掃描熱量（differential scanningcalorimetry，DSC）實驗得到了形狀記憶合金馬氏體體積分數與相變自由能之間的微分關系，并通過對Tanaka[15增量型本構方程的推導，得到了以余弦函數描述的全量型本構方程。在馬氏體正相變過程中，馬氏體體積分數 ξ 與環境溫度 T_? 材料應力 σ 之間的函數關系為

在馬氏體逆相變過程中， ξ 與 之間的函數關系為

式中： ξ_AM0 為初始馬氏體體積分數； T_ms，T_mp，T_mf 分別為馬氏體相變的起始溫度、峰值溫度和結束溫度； T_as 、 T_ap 、 T_af 分別為奧氏體相變的起始溫度、峰值溫度和結束溫度； C_m 1 C_a 為形狀記憶合金材料的特性常數； m₁ 、 m₂ 、 a₁ 、 a₂ 為與材料相變溫度有關的常數。

在一定溫度下，結合材料加載和卸載時單元內部應力大小，采用式（1）和式（2）即可獲得材料當前狀態下的馬氏體體積分數，進而可求得材料在當前狀態下的彈性模量：

E=[（1-ξ）E_a^-1+ξE_m^-1]^-1

式中： E 為材料當前彈性模量； E_a 為材料奧氏體相的彈性模量； E_m 為材料馬氏體相的彈性模量。

1.2 材料插值模型

固體各向同性材料懲罰（solid isotropicmaterialwithpenalization，SIMP）模型是一種被廣泛應用于結構拓撲優化設計的數學模型[1d]。該模型基于有限元理論，將連續體結構轉化為大量離散單元的集合，設定每個單元內部的材料密度為常數，通過指數函數構建單元的材料密度與彈性模量之間的插值關系。SIMP模型主要與變密度優化設計方法結合使用來實現對結構的拓撲優化設計。該方法把離散單元的材料密度 x_i 作為設計變量，通過SIMP模型計算得到 x_i ，并以此來判斷單元材料的有無：當 x_i=1 時，表示該處單元為實體單元；當x_i=0 時，表示該處單元無材料。之后，通過優化迭代算法完成對結構的拓撲優化設計[17]。

假設環境溫度 T 保持在奧氏體轉變結束溫度T_af 以上，材料各轉變起始與結束溫度不參與插值，則基于SIMP模型的插值函數可以表示為

式中： e_ai?σ_asi 、 σ_afi 分別為實體單元奧氏體的彈性模量、轉變起始應力和轉變結束應力；emi、σmsi＼`σ_mfi 分別為實體單元馬氏體的彈性模量、轉變起始應力和轉變結束應力； E_min 、 σ_min 為2個值極小但不為0的參量，可以防止有限元分析計算過程中出現材料彈性模量為0時所導致的結構剛度矩陣奇異； P 為懲罰因子，引入 P 的目的是為了使x_i 的最終取值更接近0或1，以便得到清晰的優化結構。

1.3 拓撲優化設計方法

1.3.1優化數學模型

結構拓撲優化的數學模型為

min c

式中： c 為優化結構的柔度，是優化設計的目標函數； 分別為結構整體的位移矩陣和剛度矩陣； u_i 人 k_i 分別為單元的位移矩陣和剛度矩陣；V₀ 表示結構初始體積； f_v 表示優化結構的體積分數約束， 0vlt;1 ； F 為結構整體受到的外載荷力矩陣。該優化設計流程以結構的最小柔度作為優化設計目標。

1.3.2靈敏度分析

在每個優化迭代步對結構進行有限元分析，計算結構柔度 Ψ_c ，再通過求導得到靈敏度。在單個單元上，目標函數 Ψ_c 對于設計變量 x_i 的靈敏度可以表示為

根據方程 F=KU ，兩邊同時對 x_i 進行求導，移項得：

將式（7）代入式（6）中，得到：

體積約束 V 關于 x_i 的靈敏度為

將單個單元計算得到的靈敏度加權疊加，進行平均計算，得到目標函數對 x_i 的靈敏度過濾公式：

式中， 為卷積算子，其計算公式為

式中： dist（i，j） 表示單元 i，j 中心的距離； N 為離

散單元數目； r_min 為過濾半徑，本文最小過濾半徑取單元長度的2倍。

1.3.3優化設計流程

在拓撲優化迭代過程中，考慮到形狀記憶合金的性質會發生非線性變化，引入如式（4）所示的插值函數，對單元關鍵力學性能指標進行插值，并借助有限元分析方法得到單元內部應力。再以單元內部應力作為參量，借助余弦型本構模型方程，對結構內部各單元彈性模量進行同步迭代，拓撲優化迭代流程如圖1所示，形狀記憶合金的材料參數如表1所示[18]。在不考慮形狀記憶合金材料非線性的傳統線彈性插值拓撲優化方法中，假設彈性模量為 55000MPa ，泊松比為 v=0.3 。

2 優化設計算例

為了驗證本文提出的優化設計方法的有效性，以二維懸臂梁和三維支撐結構為例，進行結構拓撲優化設計，并與傳統的線彈性材料優化設計方法進行比較。

2.1 二維懸臂梁

優化設計力學模型為左端固定、右端加載100N 豎直向下集中力 F 的懸臂梁結構，將結構離散為 80×40 個四邊形單元，體積分數 f_v=0.5 ，優化過程經280步迭代收斂。對2種方法所得的結構進行靜力分析，優化結構及應力分布情況如圖2所示。可以發現，本文方法得到的結構低應力區（深藍色區域）面積較少，低應力區和高應力區之間的應力差值較小，應力在結構內部分布得更均勻。

有關應力云圖的關鍵指標對比如表2所示。結果表明，本文方法相比傳統方法，所得結構整體剛度有所下降，但是結構的最大應力較小、最小應力較大，結構內部整體應力數值較大。這是由于形狀記憶合金的材料特性導致拓撲優化的迭代初期階段結構內部出現了大量呈現材料非線性的中間密度單元。這些中間密度單元大部分會因為剛度較低被優化算法刪除，從而導致結構整體剛度偏低。在后續的迭代步中，為了滿足目標函數及保證整體結構的剛度，會在結構右端生長出支撐桿以補強結構剛度。如果忽略形狀記憶合金的材料非線性，以針對傳統線彈性材料的優化方法來對形狀記憶合金結構進行拓撲優化，則在優化過程中的非線性中間密度單元會被忽略，從而導致優化結構剛度偏高，同時，在結構剛度較大區域出現應力集中現象。

優化迭代歷程如圖3所示：在迭代初期，目標函數變化存在波動；在迭代后期，迭代過程趨于平穩。

2.2 三維支撐結構

將本文方法進一步應用于三維支撐結構，設計模型如圖4所示。結構上表面受到 20MPa 的均布面載荷，下部兩腳處底面完全固定，上部受載面及以下 1mm 厚度為非設計域，其余部分為設計域。使用Abaqus軟件中的UMAT材料子程序，實現對形狀記憶合金非線性材料特性本構方程在有限元分析和拓撲優化迭代過程中的彈性模量插值過程的函數描述；使用Abaqus軟件中的Tosca拓撲優化設計模塊完成優化設計迭代流程。設置體積分數 f_v=0.2 ，采用兩種方法得到的設計結果及迭代曲線如圖5所示。可以發現，兩種優化方法所得到的橋梁幾何構型類似。本文方法獲得的結構在橋梁拱形支撐的中心部位有較密集的分支支撐，受載面邊緣部位通過一對較粗的支撐桿與主體拱形支撐連接；而傳統線彈性材料優化方法獲得的結構支撐桿分布較為平均，受載面中心部位的支撐桿分布較為稀疏，且在受載面邊緣部位分出一對額外支撐桿，其延伸到底部約束面以加固支撐。比較2種優化方法的優化迭代歷程，可以發現，在本文方法的優化迭代過程中，結構的體積在迭代初期出現了快速上升、快速下降、短暫波動、平穩上升的過程。這是由于結構內部在受載時，非線性中間密度單元剛度較低，導致局部結構不穩定。隨著迭代流程的進行，低剛度的非線性中間密度單元經過優化算法的過濾被消除，結構的成型過程趨于穩定。在傳統方法迭代過程中，由于材料為線彈性，迭代過程不存在剛度下降的非線性中間密度單元，因此，整體迭代流程較為平穩，體積平穩增加至約束范圍。導出片體模型，并將其導入Solidworks軟件，再使用曲面拼接、實體填充以及拉伸、切除等功能，對設計結果進行工程化建模，結果如圖6所示。

2.3 優化結果分析與對比

2.3.1靜力學性能分析

以圖6所示的三維結構為例，比較2種優化結構在外部載荷下的力學響應差異。使用Abaqus軟件，頂面施加 20MPa 均布載荷，底部兩腳完全固定，對分析模型進行網格劃分。在相同計算資源消耗下，六面體網格的有限元分析精度和準確度均高于四面體網格。為了保證靜力分析的準確度和后續相變演化過程仿真的計算效率，將優化設計結果的工程化模型用 1mm×lmm×lmm 的六面體C3D8單元進行網格劃分，如圖7所示。其中，本文方法、傳統方法所得結構的離散六面體實體單元個數分別為7542個和7950個。2種結構在 20MPa 外載下的應力云圖如圖8所示，其中，本文方法所得結構應力分布更加均勻，整個拱形支撐下表面的應力在 250～300MPa ，整體平均應力水平較高，未出現局部高應力區；而傳統方法所得結構的應力集中分布在拱形支撐下半部分，拱形支撐和支撐桿連接處出現局部應力集中，應力達到 400MPa 以上，而拱形支撐上半部分應力較低，低于 150MPa 。位移云圖如圖9所示，2種結構的主要變形區域都位于受載面中心區域，本文方法所得結構在受載面中段的平均位移較高，但未出現局部大位移區域；而傳統方法所得結構雖然大部分區域位移數值較低，但是在受載面中心區域出現了 0.2mm 以上的局部大變形（云圖中的紅色部分）。提取各單元的應力和位移數據進行靜力學分析，結果如表3所示。

通過靜力學分析可知，相比傳統方法，本文方法所得結構在受載時的最大應力下降 19.57% 最大位移下降 31.03% ，平均應力和平均位移則有小幅上升。這表明考慮形狀記憶合金的非線性后，優化得到的結構在應對載荷的整體剛度上有所下降，但是結構內部應力集中和局部大變形的情況有所減少，結構整體應力分布更加均勻。

2.3.2形狀記憶合金相變演化模式分析

為了進一步研究2種優化方法對形狀記憶合金的力學響應，基于Python語言的Abaqus軟件二次開發子程序，對形狀記憶合金結構在漸變外載荷下的相變過程進行仿真，對比2種優化結構在 0～ 20MPa 頂部均布載荷下結構內部的馬氏體相變情況

采用漸進迭代法，共100個迭代步，每個迭代步都在上一迭代步的基礎上對結構受載區域額外增加 0.2MPa 的載荷。每次加載結束后，通過靜力分析讀取每個單元的應力值，再利用前文所述余弦型本構模型對每個單元當前迭代步下的馬氏體體積分數 ξ 進行計算。具體分析流程如圖10所示。

將2種優化方法所得結構進行有限元網格劃分，將其導人Abaqus軟件中，并按照預先設置好的載荷施加方式添加載荷及約束條件，然后開始運行子程序。運行完畢后，得到2種結構馬氏體體積分數演化情況，如圖11和圖12所示。2種結構的平均馬氏體體積分數隨載荷增加的變化如圖13所示。在 20MPa 載荷下，2種結構的平均馬氏體體積分數如表4所示。

觀察馬氏體體積分數的演化情況可以發現：本文方法所得結構在載荷逐漸增大的過程中，馬氏體體積分數增長速度顯著變慢；在同迭代步載荷下，馬氏體體積分數相對較小。這表明在相同載荷下，本文方法所得結構內部達到馬氏體相變臨界應力的單元數量較少，發生馬氏體相變的單元數量也較少。在最大載荷（ 20MPa， 時，該結構平均馬氏體體積分數比傳統方法所得結構下降20.27% 。且在加載迭代的全過程內，本文方法所得結構內部馬氏體體積分數增長曲線較為平穩，表明內部應力變化較為均勻；而傳統方法所得結構在馬氏體體積分數增長過程中會出現階躍現象，這表明在該載荷下該結構內部部分區域出現應力集中現象。這些部位附近的應力水平較高，馬氏體相變快速發生。

從相變發生分布情況來看，本文方法所得結構發生相變的單元先集中出現在支撐結構上表面加載平面中間部位，應力傳遞遵循自上而下的規律。到第100個迭代步，即載荷為 20MPa 時，該結構發生馬氏體相變的單元均勻分布在橋梁拱形支撐區域，且發生相變的單元馬氏體體積分數為0.5～0.8 。傳統方法所得結構發生相變的單元集中出現在橋梁拱形支撐腳部區域，應力傳遞自下而上，且單元相變演化速度較快。這表明在載荷較大時，該結構拱形支撐腳部位開始出現應力集中現象，局部區域快速發生相變，并由于形狀記憶合金的超彈性特性伴隨著較大的相變變形，導致結構不穩定。

3結論

提出了一種基于SIMP模型和形狀記憶合金余弦型本構模型的非線性材料插值拓撲優化設計方法。首先，通過余弦模型對形狀記憶合金材料非線性特性進行描述，對馬氏體和奧氏體彈性模量以及4個轉變應力進行插值，建立密度設計變量與材料屬性之間的對應關系。之后，通過二維和三維結構優化算例，與忽略材料非線性條件下獲得的優化結構進行靜力力學分析對比以及基于形狀記憶合金超彈性的相變演化過程對比。研究表明：考慮了形狀記憶合金非線性優化方法所得結構對應力誘發的馬氏體相變有更優的抵抗性能。由于馬氏體彈性模量低于奧氏體，相變過程對于結構而言是一種“軟化過程”。將形狀記憶合金的材料非線性納入剛度最大化為目標函數的拓撲優化流程后，優化算法會通過改進結構內部應力分布的方式抵消“軟化過程”的影響，以提高結構剛度，滿足目標函數。“優化應力分布”體現在分析結果上，即減少了會觸發相變的高應力區域，弱化應力集中情況，將高應力區域的應力分流到低應力區域，形成大量應力水平持平或低于馬氏體相變臨界點的中間應力區域，這種優化結構整體剛度低于線彈性方法優化結構，但是減少了局部高應力區和應力集中現象，避免了應力誘發馬氏體相變帶來的局部大變形，使結構更加穩定，不易發生相變。這種優化設計方法可以應用在用形狀記憶合金作為阻尼器的元件上，提升其阻尼性能，如建筑物阻尼器和減震彈簧等，也可以用于形狀記憶合金材料的承載部件上，例如心腦血管支架，以增加其在血流沖擊和血管搏動循環載荷下的使用壽命。

參考文獻：

[1]OTSUKA K， WAYMAN C M. Shape memory materials[M]. New York： Cambridge University Press， 1999.

[2]王奇，徐志偉.SMA驅動變體機翼后緣結構力學分析 [J].南京航空航天大學學報，2015，47（6）：904-910.

[3]戴尅戎，倪誠，吳小濤，等.形狀記憶鋸齒臂環抱內固定 器的實驗研究與臨床應用[J].中華外科雜志，1994， 32（10）： 629-632.

[4]BOROSCHEKRL，FARIASG，MORONIO，etal.Effect of SMA braces in a steel frame building[J]. Journal of Earthquake Engineering， 2007，11（3）： 326-342.

[5]朱繼宏，周涵，王創，等.面向增材制造的拓撲優化技術 發展現狀與木米[J].航空制造技木，2020，63（10）：24-38.

[6]AL-TAMIMI A A， FERNANDES PR A， PEACH C， et al. Metallic bone fixation implants： a novel design approach for reducing the stress shielding phenomenon[J]. Virtual and Physical Prototyping， 2017， 12（2）： 141-151.

[7]王小華，王偉.基于拓撲優化的離心式壓縮機底座輕量 化設計[J].石油化工設備，2024，53（3）：48-51.

[8]肖蜜，李奇石，高亮，等.基于各向異性材料插值的點陣 散熱結構跨尺度拓撲優化[J].機械工程學報，2024， 60（13）： 71-80.

[9]THAMBURAJA P， ANAND L. Polycrystalline shapememory materials： effect of crystallographic texture[J]. Journal of the Mechanics and Physics of Solids，2001， 49（4）： 709-737.

[10]張楊.鎳鈦形狀記憶合金超彈性本構關系[D].大連：大 連理工大學，2013.

[11]AURICCHIO F， TAYLOR R L. Shape-memory alloys： modelling and numerical simulations of the finite-strain superelastic behavior[J]. Computer Methods in Applied Mechanics and Engineering，1997， 143（1/2）： 175-194.

[12]KAN Q H， ZHANG Y， SHI W X， et al. Functional fatigue ofsuperelasticityand elastocaloric effect for NiTi springs[J]. International Journal of Mechanical Sciences， 2024，265：108889

[13] ZHU J H， ZHANG W H， XIA L. Topology optimization in aircraft and aerospace structures design[J]. Archives of Computational Methods in Engineering， 2016，23（4）： 595-622.

[14]王振清，周博，梁文彥.形狀記憶合金的本構關系[J].金 屬學報，2007，43（11）： 1211-1216.

[15] TANAKA K. A thermomechanical sketch of shape memoryeffect：one-dimensional tensile behavior[J]. Resmech Mechanica， 1986， 2（3）： 59-72.

[16] SIGMUND O. A 99 line topology optimization code writen in Matlab[J]. Structural andMultidisciplinary Optimization， 2001， 21（2）： 120-127，

[17] BENDSOE M P， SIGMUND O. Topology optimization： theory，method， and applications[M]. 2nd ed. Berlin， Heidelberg： Springer， 2004.

[18]王振清，梁文彥，周博.形狀記憶材料的本構模型[M].北 京：科學出版社，2017.

（編輯：董偉）