三維負泊松比超材料拓撲優化設計

Topology optimization for 3D metamaterials with negative Poisson's ratio

HU Tiannan1， GUO Honghu2 （1. DepartmentofechanicalEngneeringandience，Kyoto UnversityKyoto6553o，Jan;2.Facultyofied Engineering， Waseda University， Kyoto l69-8555， Japan）

Abstract： Negative Poisson's ratio metamaterials have significant application potential in aerospace， biomedical engineering， and flexible electronics due to their anomalous mechanical behavior of expanding laterally under tension and contracting laterally under compression. However， existing research predominantly focuses on 2D or 3D isotropic designs， which struggle to meet practical engineering demands for direction-dependent material properties. This paper proposed a density-based topology optimization method for designing 3D orthotropic negative Poisson's ratio metamaterials. By constructing a novel multi-objective optimization function combined with homogenization theory， negative Poisson's ratio characteristics in three orthogonal directions for unit cell structures was achieved. First， based on the SIMP （solid isotropic material with penalization） material interpolation model， geometric constraints were introduced to ensure unit cell symmetry while eectively reducing computational scale. Second， a new objective function and optimization model were established through penalty functions. Finally， the equivalent mechanical properties were calculated using finite element homogenization under periodic boundary conditions， with design variables updated through sensitivity analysis. Numerical examples demonstrate that the optimized unit cells exhibit negative Poisson's ratio behavior in all three principal directions. This study provides a theoretical foundation for the controllable fabrication of 3D orthotropic auxetic metamaterials and expands the design and application scope of mechanical metamaterials.

Keywords： topology optimization; negative Poisson's ratio; metamaterial; homogenization method; density method

負泊松比材料作為一種力學超材料，因其在單軸載荷下呈現反常的橫向變形行為（拉伸時橫向膨脹，壓縮時橫向收縮），近年來在航空航天、生物醫學、防護裝備及智能傳感器等領域展現出廣闊的應用前景[1]。與傳統材料相比，負泊松比材料在抗剪切性、抗斷裂性和能量吸收率等方面更有優勢[2]。然而，現有研究多集中于二維結構或基于經驗設計的周期性胞元，三維負泊松比超材料的設計仍面臨幾何自由度受限、力學性能可調性不足及制造成本高昂等挑戰[3]。針對這些問題，提出結合拓撲優化方法與均勻化理論，系統研究三維負泊松比超材料的優化設計方法，旨在突破傳統設計框架，實現高性能、可定制化的負泊松比超材料結構設計。

拓撲優化作為一種基于數學模型的材料分布優化方法，能夠通過目標函數與約束條件的靈活配置，高效探索材料與結構的最優構型[4]。均勻化法通過表征微觀結構與宏觀等效性能的關系，為超材料的多尺度設計提供了理論框架。Sigmund 等[5]系統總結了拓撲優化在多物理場超材料設計中的應用，指出其在實現負泊松比、負熱膨脹等多功能耦合方面的潛力。Anaya-Jaimes等提出了基于雙向漸進結構優化的無約束設計框架，結合鄰近材料相間插值、靈敏度穩定化策略及動態添加比率調控，實現了二維正交各向異性超材料的熱膨脹系數精準定制。Murai等運用拓撲優化技術實現了可捕捉局部共振效應，并精準調控多頻電磁波行為的超材料設計。Paulino等構建了針對功能梯度材料的多尺度架構設計框架，并在一維、二維梯度及極端性能（近零剪切模量、負泊松比）材料中驗證了方法的普適性與有效性。Ha等[9]提出了一種基于手性負泊松比晶格（內凹六邊形或手性三角形結構）的材料設計方法，通過采用兩種具有不同熱膨脹系數的材料，構建雙材料梁單元，在宏觀尺度上實現了對材料整體熱膨脹系數的主動調控。Ai等[1基于參數化水平集方法和緊支徑向基函數無網格方法，提出了一種二維微結構機械超材料拓撲優化設計機制，實現了負泊松比拉脹超材料設計，并可擴展至設計具有負熱膨脹系數或頻率帶隙的超材料。Qin 等[1]以最大化初始機械阻抗水平為目標，提出了一種用于提升減振性能的超材料拓撲優化設計方法。Ye等[12]設計了一種自適應漸進變硬機械超材料，克服了傳統均勻負泊松比材料常存在的剛度突變或應力局限性，為開發具有定制化非線性力學響應的超材料提供了設計理論與方法。Jebellat等[13]研究提出了一種由智能材料構成的超結構設計，該結構不僅具有形狀記憶效應，還具備負泊松比性能，可應用于新型傳感器、執行器及生物醫學領域。Agrawal等[14]基于密度法，提出了一種魯棒拓撲優化（robust topologyoptimization，RTO）設計方法，以應對材料不確定性對負泊松比超材料性能的影響。

雖然運用拓撲優化實現負泊松比結構設計已十分成熟，三維負泊松比超材料的優化設計研究卻相對較少，限制了超材料的工程應用。Evans等[15]提出了三維正交內凹蜂窩負泊松比結構。基于拉一扭耦合效應，Duan 等[提出了一種新型三維負泊松比材料設計方法，通過將多層二維負泊松比材料與三維手性鏈結構耦合連接，構建出兼具負泊松比和正泊松比性能的三維材料體系。然而，這些三維超材料設計都并非基于拓撲優化設計，難以通過合理的材料分布實現最優的機械性能。Clausen等[17]運用幾何非線性拓撲優化技術，首次實現了泊松比在 -0.8～0.8 范圍內可編程調控的三維超材料設計，并通過參數化超級橢圓結構優化與多噴嘴連續打印策略，解決了傳統拓撲優化超材料結構難以規模化制造的難題。Zhang等[18]利用GPU并行計算優化設計了高分辨率的三維負泊松比微結構。Vogiatzis等[9]采用調和水平集方法實現了多材料負泊松比超材料的優化設計，顯著提升了結構的剛度與能量吸收性能。雖然基于拓撲優化的三維負泊松比超材料設計已逐漸獲得關注，如何便捷有效地實現三維正交異性負泊松比性能仍然需要進一步研究。

針對上述挑戰，本文提出一種基于拓撲優化與均勻化理論的三維負泊松比超材料設計方法。首先，通過密度法構建參數化設計空間，結合均勻化理論計算等效力學性能；其次，引入調節因子調控結構的負泊松比行為，實現拉壓載荷下力學響應的對稱性優化；最后，通過光固化（stereolithography，SLA）技術驗證設計方案的可行性。本文的創新點包括：a.提出了一種兼顧幾何自由度與計算效率的三維拓撲優化框架；b.構建了新的三維負泊松比拓撲優化目標函數；c.通過增材制造驗證了優化結構的可制造性。研究結果將為高性能負泊松比超材料的定制化設計提供理論支撐與技術路徑。

1 三維負泊松比超材料拓撲優化設計方法

1.1 均勻化法

利用均勻化法計算超材料單胞的等效彈性矩陣，通過施加周期性邊界條件對微結構單胞進行有限元分析，結構平衡方程為

KU=f

式中： K 為單胞的整體剛度矩陣； U 為單胞的位移場； f 為由均勻應變場引起的外力，可以表示為

式中： B 為應變-位移矩陣； Y 為單胞； Y_e 為單胞中第 e 個單元的體積； D_e 為單元 e 的彈性矩陣；m 為單元總個數。針對三維問題，應變 ε=[ε₁，ε₂，ε₃ ε₄， ε₅ ， ε₆] ， ε₁=[1，0，0，0，0，0]^T ， ε₂=[0，1，0，0，0，0]^T ，ε₃=[0，0，1，0，0，0]^T ， ε₄=[0，0，0，1，0，0]^T ， ε₅=[0，0，0，0，1，0]^T ，ε₆=[0，0，0，0，0，1]^T 。

得到位移場后，利用均勻化法，超材料的等效彈性矩陣 D^H 可以表示為

式中： |Y| 為單胞的總體積； u_e⁰ 為通過施加單位應變計算得到的節點位移場； U_e 為通過施加全局單位應變計算得到的位移場[20]； k_e 為單元剛度矩陣。

k_eu_e⁰=f_e

矩陣 D^H 中的彈性張量 D_rs 可以表示為

式中： r 對應應力分量， r=1，2，3，4，5，6 s 對應應變分量， s=1，2，3，4，5，6

D^H=

式中： c 為彈性系數。對于施加對稱性約束后的三維正交各向異性材料而言， C₁₁₂₂^H=C₂₂₁₁^H ，C₁₁₃₃^H=C₃₃₁₁^H ， C₂₂₃₃^H=C₃₃₂₂^H°

通過對彈性矩陣中相關項的計算，可得到3個方向的泊松比的計算公式為

1.2 優化數學模型

實現正交各向異性負泊松比超材料設計的關鍵是彈性矩陣中 為負數，理論上，通過最小化 C₁₁₂₂^H+C₁₁₃₃^H+C₂₂₃₃^H 可實現負泊松比設計。因此，負泊松比超材料設計的核心轉變為如何通過迭代，實現 C₁₁₂₂ ！ C₁₁₃₃ 和 C₂₂₃₃ 的最小化設計。為了改善迭代過程的收斂性，基于文獻[14]和[21]，引人松弛項 wβ^l（C₁₁₁₁^H+C₂₂₂₂^H+C₃₃₃₃^H） ，構建新的目標函數。最終，通過引入松弛函數構造的新目標函數為

式中： β 為常數，取值范圍為（0，1）；指數 l 為迭代次數； w 為調節參數。

為了調節3個正交方向上的泊松比數值，引入額外的3個調節參數 w₁，w₂ 及 w₃ ，將式（8）修改為

將固體材料體積設定為上限值，其作為不等式約束，具體上限值可根據實際應用需求進行靈活調整。因此，最終的優化數學模型為

式中： u 是實體材料的體積； u_f 是體積分數； u₀ 是微結構的總體積； x 是設計變量； x_min 是為了避免數值不穩定現象設置的1個小的正值。體積分數約束定義了微觀結構中材料的最大體積。基于SIMP法，單元剛度矩陣 k_e 可以表示為

k_e=x_e^p（k_s-k₀）+k₀

式中： k_s 為實體單元的剛度矩陣； k₀ 為空材料單元對應的剛度矩陣； p 為懲罰因子。

此外，為了保證最終設計結果的正交各向異性，在 xy 、yz、 zx 面對單胞結構施加幾何對稱約束，如圖1所示。由于結構是對稱的，因此，拓撲優化對應的獨立設計變量為圖1藍色部分所示的單胞微結構1/8部分。對稱約束的施加，不僅保證了設計結果的正交性，也減小了設計變量，提高了計算效率。

1.3 并行計算

三維結構拓撲優化相較于二維結構，其計算規模呈指數增加。為了實現三維負泊松比超材料的高分辨率設計，采用開源并行框架對上述問題進行求解。使用開源軟件FreeFEM ⁺⁺ 對均勻化問題進行離散，采用PETSc并行求解線性方程組。并行求解算法邏輯如下：首先，使用圖分割包，如METIS，對網格進行分解，并指派給CPU各核心；隨后，由FreeFEM執行偏微分方程離散，并將離散后的線性代數方程組傳遞給PETSc進行并行求解；最后，將分別求解得到的數據進行回傳處理，得到最終分析結果。本文所有數值算例均在DELL桌面級工作站Precision3660上完成，其主要配置如下：中央處理器為Inteli913900K，基頻為 3.0GHz ；內存為128GB；顯卡為NVIDIAGeForceRTX4080。

1.4 優化算法及流程

優化算法的流程如圖2所示。其中，有限元離散和線性方程求解通過并行計算處理，兩者均基于MPI（message passing interface）實現。

優化算法的描述如下：

a.建立設計域和參數初始化

b.對設計域進行離散，并網格進行分解。

c.有限元法求解控制方程。

d.計算目標函數，進行收斂判斷。如果目標

函數收斂，則優化結束，否則跳轉至步驟e進行迭代直至收斂。收斂條件：優化迭代的目標函數連續兩次的差值小于收斂容差，或者迭代次數到達設定的最大值 N_c

e.計算靈敏度。

f.基于移動漸近線算法（methodofmovingasymptotes，MMA）更新設計變量，并跳轉至步驟c。

2 數值算例

2.1 算例設置與優化過程

為了驗證所提設計方法的有效性，選取邊長為1的三維單胞微結構作為設計域，采用結構化的四面體網格進行離散，單邊網格數為30，網格總數為 32.4×10⁴ 。為了實現結果的對稱性，在單胞構型上施加三向對稱面 （xy，yz，zx） 約束并構建網格，如圖3所示。優化數學模型中的參數設定為：松弛因子 β=0.8 ，調節參數 w=1000 ， w₁= 3， w₂=1 ， w₃=4 ，體積分數 u_f=0.5 。材料屬性為：彈性模量 E₀=0.91 ，泊松比 μ=0.3 。

迭代歷程及優化結果如圖4所示。在初始階段，迭代步 llt;50 時，設計域內材料呈均勻分布，3 個正交方向的泊松比 μ_xy，μ_yz 和 μ_zx 相等，為材料本身的泊松比0.3。隨著迭代推進，迭代步 l 為50～200 時，目標函數通過懲罰項動態調整彈性矩陣項 C₁₁₁₁^H+C₂₂₂₂^H+C₃₃₃₃^H ，驅動泊松比向負值演變。

如圖5所示，至迭代步 l=200 時， μ_zx 率先下降到 -0.5 ，表明 zx 方向上率先形成穩定的拉脹結構。最終，迭代步 lgt;300 時，目標函數趨于穩定，3個方向上的泊松比分別達到 μ_xy=-0.558 ，μyz=-0.161，μzx=-0.801，實現三維正交各向異性負泊松比超材料結構逆向設計。可以看到，最終獲得的3個方向的負泊松比數值間的比值關系基本滿足設計要求 w₁=3 ， w₂=1 ， w₃=4 。

2.2 優化構型與力學機理

圖6為優化后的單胞結構及其與經典二維手性超材料的對比。圖 6（a）～（b） 所示的優化結果呈現出獨特的空間螺旋拓撲結構： z 軸方向通過交叉肋條形成類似二維手性結構的扭轉耦合效應，而x、 y 軸方向通過弧形空腔實現橫向膨脹抑制。這種三維耦合機制使得單軸拉伸時，交叉肋條的彈性變形引發多向收縮響應。值得注意的是，優化結構的材料分布與二維手性超材料具有幾何相似性，如圖6（c所示，驗證了目標函數對負泊松性能的捕捉能力。

由于優化是基于密度法進行的，從圖6（b）可以看到，優化后的單胞有較為明顯的鋸齒形邊界面，這一現象可以通過引入額外的技術來消除，如貼體網格技術。然而，考慮到單胞的幾何尺寸，圖6（b）所示的鋸齒形邊界面在實際制造中并不會變得明顯，且對結構的性能影響較小，因此，在本文中不作額外的處理

進一步將優化單胞進行 3×3 周期擴展，如圖7所示，宏觀等效泊松比呈現顯著方向依賴性：沿主對角線方向（ μ_zx=-0.801 ）的拉脹效應最強，而yz方向 μ_yz=-0.161 ）受對稱約束影響較弱。這種各向異性可通過調整體積分數和幾何約束進行定制，為多軸載荷場景提供設計靈活性。

3 超材料增材制造

為了驗證三維負泊松比超材料的可制造性，使用3D技術對設計結果進行增材制造。首先，將文中的優化結果導出為STL文件；然后，將STL文件導入切片軟件以生成G代碼；最后，將G代碼導入3D打印機打印結構。3D打印機的打印最高分辨率為 16μm ，打印精度為： 20～85μm （尺寸小于 50mm 零件）； 200μm （全尺寸模型， 260mm 零件），成型材料為光敏樹脂。3D打印的結果如圖8所示。從圖中可以看出，打印的結構具有很高的分辨率，與優化結果基本一致，驗證了優化結果的可制造性。

4結論

本文提出了一種基于密度法的三維正交各向異性負泊松比超材料并行拓撲優化框架，通過理論推導、數值仿真與增材制造驗證，得出以下結論：a.通過融合幾何對稱約束與動態懲罰函數，成功實現了三維單胞在3個正交方向的負泊松比各向異性調控 'μ_xy=-0.558 ， μ_yz=-0.161 ， μ_zx=-0.801 。優化結構呈現空間螺旋拓撲，其力學行為與二維手性超材料具有機理一致性，驗證了目標函數的物理合理性。

b.采用PETSc并行求解器與MPI任務分配策略，有效縮短了計算時間，相較傳統串行算法效率提升8倍，為高分辨率設計提供技術支撐。

c.增材制造樣件證實了優化結構的可制造性，其多向拉脹特性可應用于定制化緩沖材料、柔性傳感器及航空航天輕量化結構。通過調節體積分數與對稱約束，可進一步擴展至負熱膨脹、能量吸收等多功能耦合設計。

未來工作將聚焦于非線性大變形超材料拓撲優化算法開發，并探索金屬基超材料的SLA成形工藝，以提升結構的承載能力與耐久性。

參考文獻：

[1]LAKES R. Foam structures with a negative Poisson's ratio[J]. Science， 1987， 235（4792）： 1038-1040.

[2]RENX，DASR， TRANP， et al. Auxetic metamaterials and structures： a review[J]. Smart Materials and Structures， 2018，27（2）： 023001.

[3]楊智春，鄧慶田.負泊松比材料與結構的力學性能研究 及應用[J].力學進展，2011，41（3）：335-350.

[4]BENDSOE M P， KIKUCHI N. Generating optimal topologies in structural design using a homogenization method[J]. Computer Methods in Applied Mechanics and Engineering， 1988， 71（2）： 197-224.

[5]SIGMUND O， MAUTEK. Topologyoptimization approaches：a comparative review[J]. Structuraland Multidisciplinary Optimization， 2013， 48（6）： 1031-1055.

[6]ANAYA-JAIMES L M， VICENTE W M， PAVANELLO R.Metamaterials design with a desired thermal expansion using a multi-material BESO method[J]. Structural and Multidisciplinary Optimization， 2022， 65（12）： 355.

[7]MURAI N， NOGUCHI Y， MATSUSHIMA K， et al. Multiscale topology optimizationofelectromagnetic metamaterialsusingahigh-contrast homogenization method[J]. Computer Methods in Applied Mechanics and Engineering， 2023， 403： 115728.

[8]PAULINO GH， SILVA ECN，LE CH. Optimal design of periodic functionally graded composites with prescribed properties[J]. Structural and Multidisciplinary Optimization，2009， 38（5）： 469-489.

[9]HA C S， HESTEKINE，LI J H， et al. Controllable thermal expansion of large magnitude in chiral negative Poisson's ratio latices[J]. Physica Status Solidi （B）， 2015， 252（7）： 1431-1434.

[10] AI L， GAO X L. Topology optimization of 2-D mechanical metamaterials using a parametric level set method combined with a meshfree algorithm[J]. Composite Structures， 2019， 229： 111318.

[11] QIN H X， YANG D Q. Vibration reduction design method of metamaterials with negative Poisson's ratio[J]. Journal of Materials Science， 2019， 54（22）： 14038-14054.

[12] YE ML， GAO L， LI H. A design framework for gradually stiffer mechanical metamaterial induced by negative Poisson's ratio property[J].Materials amp; Design，2020，192： 108751.

[13] JEBELLAT E，BANIASSADI M， MOSHKI A， et al. Numerical investigation of smart auxetic three-dimensional meta-structures based on shape memory polymers via topology optimization[J]. Journal of Intelligent Material Systems and Structures， 2020，31（15）： 1838-1852.

[14] AGRAWAL G， GUPTA A， CHOWDHURY R，et al. Robust topology optimization of negative Poisson’ s ratio metamaterialsundermaterialuncertainty[J].Finite Elements in Analysis and Design， 2022， 198： 103649.

[15]EVANS K E， NKANSAH M A， HUTCHINSON I J. Auxetic foams： modelling negative Poisson’ s ratios[J]. Acta Metallurgica et Materialia， 1994， 42（4）： 1289-1294

[16] DUAN S Y，XI L，WEN W B，et al. A novel design method for 3D positive and negative Poisson's ratio material based on tension-twist couplingeffects[J]. Composite Structures， 2020，236： 111899.

[17] CLAUSEN A， WANG F W， JENSEN J S， et al. Topology optimized architectures with programmable Poisson's ratio over large deformations[J]. Advanced Materials， 2015， 27（37）： 5523-5527.

[18] ZHANG D， ZHAI X Y， LIU L G， et al. An optimized， easy-to-use， open-source GPU solver for large-scale inversehomogenizationproblems[J]. Structuraland Multidisciplinary Optimization， 2023， 66（9）： 207.

[19] VOGIATZIS P， CHEN S K， WANG X， et al. Topology optimization of multi-material negative Poisson’ sratio metamaterials using a reconciled level set method[J]. Computer-Aided Design，2017， 83： 15-32.

[20] ZHANG H， DING X H， WANG Q， et al. Topology optimization of composite material with high broadband damping[J]. Computers amp; Structures， 2020，239： 106331.

[21] ZHANG H K， LUO Y J， KANG Z. Bi-material microstructural design of chiral auxetic metamaterials using topology optimization[J]. Composite Structures， 2018，195： 232-248.

（編輯：董偉）