基于固定网格法的力–热多功能超材料拓扑优化设计及多级晶格填充

基金项目

国家自然科学基金面上项目(12272200);北京优解未来科技有限公司项目(20212002316)。

中图分类号:

V25TB3

文献标识码:

A

通信作者

杜建镔,副教授,工学博士,研究方向为结构拓扑优化和超材料设计。

编辑

责编 :晓月

引文格式

陈雪乾, 任宏塬, 张海西, 等. 基于固定网格法的力–热多功能超材料拓扑优化设计及多级晶格填充[J]. 航空制造技术, 2025, 68(15): 82–91.

Topology Optimization Design of Mechanical–Thermal Multifunctional Metamaterials Based on Fixed Grid Method and Multi Level Lattice Filling

Citations

CHEN Xueqian, REN Hongyuan, ZHANG Haixi, et al. Topology optimization design of mechanical–thermal multifunctional metamaterials based on fixed grid method and multi level lattice filling[J]. Aeronautical Manufacturing Technology, 2025, 68(15): 82–91.

航空制造技术    第68卷    第15期    82-91
Aeronautical Manufacturing Techinology    Vol.68    No.15 : 82-91
DOI: 10.16080/j.issn1671-833x.2025.15.082
论坛 >> 超材料(FORUM >> Metamaterials)

基于固定网格法的力–热多功能超材料拓扑优化设计及多级晶格填充

  • 陈雪乾 1
  • 任宏塬 1
  • 张海西 2
  • 周平章 2
  • 杜建镔 1,3
1.清华大学航天航空学院北京 100084
2.北京优解未来科技有限公司北京 100176
3.清华大学柔性电子技术国家级重点实验室北京 100084

通信作者

杜建镔,副教授,工学博士,研究方向为结构拓扑优化和超材料设计。

基金项目

国家自然科学基金面上项目(12272200);北京优解未来科技有限公司项目(20212002316)。

中图分类号:

V25TB3

文献标识码:

A

引文格式

陈雪乾, 任宏塬, 张海西, 等. 基于固定网格法的力–热多功能超材料拓扑优化设计及多级晶格填充[J]. 航空制造技术, 2025, 68(15): 82–91.

摘要

针对超材料的复杂几何形状和多层次等特性,提出了基于固定网格技术的超材料仿真分析和拓扑优化方法,无须人工划分贴合几何外形的有限元网格,即可大幅降低前处理时间。采用基于四叉树/八叉树的局部自适应加密技术,在不增加有限元分析计算量的前提下,既保证了结构分析精度,同时可获得更高分辨率的拓扑优化结果。基于该方法,建立了同时考虑传力、散热的多功能超材料拓扑优化模型,并以三周期极小曲面(Triply periodic minimal surface,TPMS)复合加筋结构为优化对象进行了数值实现和验证。本文方法已在国产全自主CAX工业软件OptFuture中实现。此外,面向超材料结构的多级减重需求,依托OptFuture软件实现了超材料多级晶格设计填充,进一步拓展了超材料的轻量化设计潜力。

关键词

固定网格法;拓扑优化;超材料;力热多功能设计;多级晶格填充;

Topology Optimization Design of Mechanical–Thermal Multifunctional Metamaterials Based on Fixed Grid Method and Multi Level Lattice Filling

  • CHEN Xueqian 1
  • REN Hongyuan 1
  • ZHANG Haixi 2
  • ZHOU Pingzhang 2
  • DU Jianbin 1, 3
1.School of Aerospace Engineering, Tsinghua University, Beijing 100084, China
2.Beijing OptFuture Technology Co., Ltd, Beijing 100176, China
3.State Key Laboratory of Flexible Electronics Technology, Tsinghua University, Beijing 100084, China

Citations

CHEN Xueqian, REN Hongyuan, ZHANG Haixi, et al. Topology optimization design of mechanical–thermal multifunctional metamaterials based on fixed grid method and multi level lattice filling[J]. Aeronautical Manufacturing Technology, 2025, 68(15): 82–91.

Abstract

Targeting the complex geometric shapes and multi-level characteristics of metamaterials, a simulation and topology optimization method based on fixed grid techniques is proposed. Users do not need to manually partition finite element meshes which conform to geometric shapes, thereby significantly reducing preprocessing time. By employing quadtree/octree based local adaptive refinement, the accuracy of structural analysis can be ensured without increasing finite element analysis computational costs. And higher resolution topology optimization results can be generated. Based on the proposed method, a multifunctional metamaterial topology optimization model considering both load bearing and heat conduction was established. Numerical implementation and verification were conducted using a triply periodic minimal surface (TPMS) composite rib-reinforced structure as the optimization object. The proposed method has been implemented in OptFuture, a domestically developed and fully autonomous CAX industrial software. In addition, to address the multi-level weight reduction requirements of metamaterial structures, OptFuture is used to achieve multi-level lattice design and filling of metamaterials, thereby further expanding the potential for lightweight design of metamaterials.

Keywords

Fixed grid method; Topology optimization; Metamaterials; Mechanical–thermal multifunction design; Multi level lattice filling;



超材料一般是指通过人工微结构设计而非化学组分设计组成的复合材料或结构,通常能实现传统材料无法或很难具有的属性[  胡更开, 刘晓宁. 弹性超材料设计与波动控制[M]. 北京: 科学出版社, 2024.HU Gengkai, LIU Xiaoning. Design and fluctuation control of elastic metamaterials[M]. Beijing: Science Press, 2024.
1
]
(例如轻质高承载[  EVANS A G, HUTCHINSON J W, FLECK N A, et al. The topological design of multifunctional cellular metals[J]. Progress in Materials Science, 2001, 46(3–4): 309–327.
2
]
、吸能抗冲击[  程乾, 尹剑飞, 温激鸿, 等. 极小曲面力学超材料抗冲吸能特性分析[J]. 动力学与控制学报, 2023, 21(7): 43–50.CHENG Qian, YIN Jianfei, WEN Jihong, et al. Impact resistance and energy absorption of mechanical metamaterials with minimal surfaces[J]. Journal of Dynamics and Control, 2023, 21(7): 43–50.
 LI X W, YU X, ZHAI W. Additively manufactured deformation-recoverable and broadband sound-absorbing microlattice inspired by the concept of traditional perforated panels[J]. Advanced Materials, 2021, 33(44): 2104552.
3-4
]
、定制变形[  HAGHPANAH B, SALARI-SHARIF L, POURRAJAB P, et al. Multistable shape-reconfigurable architected materials[J]. Advanced Materials, 2016, 28(36): 7915–7920.
 FRENZEL T, KADIC M, WEGENER M. Three-dimensional mechanical metamaterials with a twist[J]. Science, 2017, 358(6366): 1072–1074.
 REN X, DAS R, TRAN P, et al. Auxetic metamaterials and structures: A review[J]. Smart Materials and Structures, 2018, 27(2): 023001.
5-7
]
、热性能可控[  PERALTA I, FACHINOTTI V D, ÁLVAREZ HOSTOS J C. A brief review on thermal metamaterials for cloaking and heat flux manipulation[J]. Advanced Engineering Materials, 2020, 22(2): 1901034.
 LI Y, LI W, HAN T C, et al. Transforming heat transfer with thermal metamaterials and devices[J]. Nature Reviews Materials, 2021, 6: 488–507.
8-9
]
、电磁隐身[  CHEN M J, PEI Y M, FANG D N. Design, fabrication, and characterization of lightweight and broadband microwave absorbing structure reinforced by two dimensional composite lattice[J]. Applied Physics A, 2012, 108(1): 75–80.
10
]
、声学降噪[  ZIGONEANU L, POPA B I, CUMMER S A. Three-dimensional broadband omnidirectional acoustic ground cloak[J]. Nature Materials, 2014, 13(4): 352–355.
11
]
等特性),在航空航天等复杂多物理场环境具有较高的应用潜力。

在超材料的设计方法方面,除了传统的经验式预设结构配置的调参设计方法、参数优化设计方法以及形状优化设计方法,近年来拓扑优化方法被越来越多地引入到超材料设计中。相较传统设计方法,拓扑优化因可以彻底改变结构的连接方式和构型从而提供了更大的灵活性和自由度。常见的拓扑优化方法主要有基于密度描述的方法、基于边界演化的方法,以及基于几何组件描述的方法。基于密度描述的方法是以每个单元或节点的密度为设计变量,通过材料的密度分布表征结构拓扑[  BENDSØE M P, SIGMUND O. Material interpolation schemes in topology optimization[J]. Archive of Applied Mechanics, 1999, 69(9): 635–654.
 SIGMUND O. A 99 line topology optimization code written in Matlab[J]. Structural and Multidisciplinary Optimization, 2001, 21(2): 120–127.
12-13
]
;基于边界演化的方法中,结构边界通常通过低维平面切割高维曲面而形成,结构拓扑随着更高维度函数的变化而演变[  SETHIAN J A, WIEGMANN A. Structural boundary design via level set and immersed interface methods[J]. Journal of Computational Physics, 2000, 163(2): 489–528.
 WANG M Y, WANG X M, GUO D M. A level set method for structural topology optimization[J]. Computer Methods in Applied Mechanics and Engineering, 2003, 192(1–2): 227–246.
14-15
]
;基于几何组件描述的方法是将结构分解成若干个离散构件,通过构件的移动、变形、消失和融合等形式实现结构拓扑的变化[  GUO X, ZHANG W S, ZHONG W L. Doing topology optimization explicitly and geometrically—A new moving morphable components based framework[J]. Journal of Applied Mechanics, 2014, 81(8): 081009.
 ZHANG W S, YUAN J, ZHANG J, et al. A new topology optimization approach based on moving morphable components (MMC) and the ersatz material model[J]. Structural and Multidisciplinary Optimization, 2016, 53(6): 1243–1260.
16-17
]

在具体超材料设计实现方面,近年的工作主要集中于3类:第1类是基于经验的调参式设计和验证。Zhang等[  ZHANG P, Biligetu, QI D X, et al. Mechanical design and energy absorption of 3D novel hybrid lattice metamaterials[J]. Science China Technological Sciences, 2021, 64(10): 2220–2228.
18
]
构建了不同类型的基本立方单元混合的梁系点阵超材料,研究结果表明,新型混合点阵的变形机制有助于提升致密应变和能量吸收效率,并获得了应变硬化和双线性特性。Wei等[  WEI Y L, YANG Q S, LIU X, et al. Multi-bionic mechanical metamaterials: A composite of FCC lattice and bone structures[J]. International Journal of Mechanical Sciences, 2022, 213: 106857.
19
]
通过仿生设计,将具有轻质和高比强度的面心对称结构与具有高韧性的同心圆结构相结合,通过调整硬相和软相材料的比例,使得吸能效果大幅提升。Li等[  LI D M, QIN R X, XU J X, et al. Improving mechanical properties and energy absorption of additive manufacturing lattice structure by struts’ node strengthening[J]. Acta Mechanica Solida Sinica, 2022, 35(6): 1004–1020.
20
]
针对梁系点阵结构在外荷载作用下的应力集中问题,提出添加圆角、圆弧过渡等策略;结果表明,优化后结构的变形模式从剪切破坏变为逐层坍塌,结构的力学性能改善显著。第2类是基于简单的参数优化的设计。Messner[  MESSNER M C. Optimal lattice-structured materials[J]. Journal of the Mechanics and Physics of Solids, 2016, 96: 162–183.
21
]
提出一种用于描述梁系点阵材料周期结构的参数化模型,并通过参数调整和优化得到了一种弹性各向同性且刚度最大的微结构。Berger等[  BERGER J B, WADLEY H G, MCMEEKING R M. Mechanical metamaterials at the theoretical limit of isotropic elastic stiffness[J]. Nature, 2017, 543(7646): 533–537.
22
]
采用启发式优化策略,提出了一种板系构型微结构,研究结果显示,该构型具有各向同性性质,且材料刚度能达到理论上界。第3类是基于拓扑优化的设计。Long等[  LONG K, DU X R, XU S Q, et al. Maximizing the effective Young’s modulus of a composite material by exploiting the Poisson effect[J]. Composite Structures, 2016, 153: 593–600.
23
]
提出了一种用于最大化且具有多材料的复合微结构的有效杨氏模量的超材料设计方法。Wang等[  WANG Y Q, LUO Z, ZHANG N, et al. Topological shape optimization of microstructural metamaterials using a level set method[J]. Computational Materials Science, 2014, 87: 178–186.
24
]
将均匀化方法与参数化水平集方法相结合,用于负泊松比超材料设计。Liu等[  LIU Y, WANG Y Z, REN H Y, et al. Ultrastiff metamaterials generated through a multilayer strategy and topology optimization[J]. Nature Communications, 2024, 15(1): 2984.
25
]
采用边界惩罚法,提出了多层壳系超材料设计与优化策略,超材料单胞的优化结果可以自动形成梁、板、壳形式的综合体结构,并能获得逼近理论极限的材料刚度。

目前围绕超材料设计已开展了大量研究工作,然而现有研究还存在以下3项困难或不足。

(1)由于超材料的复杂几何形状和多层次等特性,传统的数值模拟方法通常会导致计算成本高且复杂。例如,采用传统贴体网格的有限元方法,首先需要划分贴合几何外形的网格,然后再对结构进行分析或优化设计。对于超材料这类复杂结构,高质量的贴体网格划分(尤其是六面体网格)通常需要花费大量时间成本。

(2)在超材料设计中,模型的离散化通常会影响优化的精度和效率。即为了获得性能更好、细节更丰富的优化结果,通常需要更精细的有限元网格,这将会大大增加分析和优化的成本。虽然引入多重网格技术(即分析和优化采用粗细两套不同的全局网格)可以在一定程度上缓解此问题,却无法从根本上解决问题,同时还会带来精度的下降。

(3)现有的基于优化理论的超材料设计,其建立的优化数学模型通常仅考虑单一功能的指标进行优化求解,对多功能需求的考虑方式通常是在单目标优化求解的基础上对结果进行仿真验证和微调,难以充分挖掘超材料设计潜力从而最大程度地满足多功能设计需求。

针对以上困难和不足,本文提出了基于固定网格技术的超材料仿真分析和拓扑优化方法。固定网格法是将经典的有限单元法和虚拟域方法结合,将复杂的物理域嵌入到一个规则的扩展域中,在扩展域上进行网格划分及结构分析[  DÜSTER A, PARVIZIAN J, YANG Z, et al. The finite cell method for three-dimensional problems of solid mechanics[J]. Computer Methods in Applied Mechanics and Engineering, 2008, 197(45–48): 3768–3782.
26
]
。由于扩展域是规则的矩形或立方体,因此网格划分几乎是无成本的,从而可实现复杂结构的高效分析。此外,为了高精度的计算边界单元积分,固定网格法中采用基于四叉树(对于二维问题)或八叉树(对于三维问题)的局部自适应加密技术,将边界单元细分为子域。该类细分策略并不会引入新的自由度,即不会导致有限元分析计算量的增加,可较好的平衡计算精度和效率。同时,细分后的子域可自然的作为拓扑优化的设计变量,从而实现在不增加有限元分析计算量的前提下获得更高分辨率的拓扑优化结果。基于所发展的固定网格技术,本文建立了同时考虑传力与散热的多功能超材料拓扑优化模型并进行了数值实现和验证,上述方法已在国产全自主版权的CAX工业软件OptFuture中实现落地。此外,面向超材料结构的多级减重需求,依托OptFuture软件,对拓扑优化后的结构进行了多级晶格设计填充,进一步提升了轻量化设计水平。

1     基于固定网格法的结构分析与拓扑优化理论

1.1     基于固定网格法的结构分析

固定网格法是将待分析的物理域置于一个更大的、几何规则的区域内,如图1所示。类似于有限单元法,扩展区域可以离散成独立于初始物理域的网格单元。由于这些嵌入区域的单元可能会被初始物理域的边界分割,因此它们不一定满足初始物理域的一般几何性质。为将固定网格单元与经典的有限单元区分开来,将其命名为“胞元”,所有的胞元组成了完整的扩展域。按照胞元与物理域边界的位置关系与分割情况,可以将胞元分为3类,即物理胞元、边界胞元、虚拟胞元,如图1所示。

图1     固定网格法示意图
Fig.1     Schematic diagram of fixedgrid method

本文将以线弹性静力学问题为例推导固定网格下的有限元方程。定义物理域为,规则的扩展域为′。扩展域上线弹性体静力学微分方程弱形式为

ΩδεijσijhdΩ=ΓtδuitidΓ+ΩδuibihdΩ
(1)

式中,δεij为虚应变;σij为应力;δui为虚位移;ti为给定的面力载荷;Γt为面力边界(Γ为面力边界上的1个微元);bi为给定的体力载荷;h为示性阶跃函数。

h(x)={1  (xΩ)0  (other)
(2)

式中,x为任意点的坐标。对应的边界条件为

σijnj=ti(on Γt)ui=u¯i(on Γu)
(3)

式中,nj为边界外法线单位矢量;ui为位移矢量;u¯i为给定的位移矢量;Γu为给定的位移边界。面力边界条件已经自然地蕴含在虚功原理中,无须特殊处理。位移边界条件可通过罚函数法施加,此时弱形式可改写为

ΩδεijσijhdΩ+βΓuδui(uiu¯i)dΓ=ΓtδuitidΓ+ΩδuibihdΩ
(4)

式中,β为惩罚系数。将扩展域离散为有限胞元,胞元内任意点的位移向量ue及其变分δue可由该胞元对应的节点向量qeδqe插值得到。

ue=Nqeδue=Nδqe
(5)

式中,N为形函数,其形式与传统有限元方法完全相同。将分别列出,可得到如下形式的刚度方程。

KU=FK=e=1neleKeF=e=1neleFe
(6)

式中,K为总体刚度矩阵;UF为定义在胞元节点上的全局位移列阵和全局载荷列阵;KeFe为第e个胞元的刚度矩阵和载荷向量,其具体形式为

Ke=ΩeBTDBhdΩ+βΓuNTNdΓFe=ΓteNTtdΓ+ΩeNTbhdΩ
(7)

式中,h为示性阶跃函数。对于物理胞元,取h=1;对于虚拟胞元,取h=0;对于边界胞元,由于物理域的边界不规则,无法直接应用高斯积分等方式进行体积分计算,因此采用四叉树/八叉树方法将其划分为多个子域进行计算,使子域尽量逼近物理域边界。图2(a)为二维四节点边界胞元,划分一层四叉树后如图2(b)所示。将新的积分域Ωii=1,…,4)称为子域。此时子域中任意点的位移仍然是通过胞元节点位移插值得到,即仍然是通过公式进行计算,而不是通过子域节点位移插值得到,因此划分子域后的胞元自由度仍然是8个(即子域划分不会导致自由度的增加),此时刚度矩阵可通过式(8)计算。

ΩeBTDBhdΩ=i=13ΩiBTDBdΩ+Ω4BTDBhdΩ
(8)

图2     采用二层四叉树将胞元划分为子域进行刚度矩阵计算
Fig.2     Partitioning cell into subdomains for stiffness matrix calculation based on the quadtree algorithm with a recursion-depth of 2

式(8)中等号右端第1项,由于积分域为规则的矩形,并且函数h取值均为1,因此可以直接采用高斯积分精确计算。而等号右端第2项,由于函数h的间断性及结构边界不规则性,仍然无法精确计算,因此对4继续采用四叉树细分,如图2(c)所示,此时刚度矩阵可通过式(9)计算。

ΩeBTDBhdΩ=i=13ΩiBTDBdΩ+i=14Ω4iBTDBhdΩ
(9)

以此类推,理论上经过无限次四叉树细分,可精确计算刚度矩阵,实际编程中则需给定一个细分上限,通常经过6次以内的细分即可得到较高精度的结果。

稳态热传导问题的求解方法与线弹性静力问题类似,因此略去推导过程,最终有限元方程为

KtT=Ft
(10)

式中,Kt为总体传热矩阵;TFt分别为在胞元节点上的全局温度列阵和全局热载荷列阵。

为了验证固定网格法的求解精度及效率,本文针对图3(a)所示的4×4阵列三周期极小曲面(Triply periodic minimal surface,TPMS)结构,分别采用固定网格法和Abaqus 2023软件进行线性静力分析。静力工况如图3(b)所示,结构底部位移完全固定,顶部施加沿Z轴负方向的1 MPa均布压力,材料弹性模量为71000 MPa,泊松比为0.3。其中,固定网格法采用尺寸为0.5 mm的二阶六面体单元,Abaqus仿真分别采用平均尺寸为0.5 mm和0.2 mm的二阶四面体单元。计算域划分规则六面体固定背景网格与不规则四面体贴体网格示意图分别如图3(c)和(d)所示,所有模型均采用16核心CPU并行计算。

图3     极小曲面阵列结构、边界条件和网格划分示意图
Fig.3     Schematic diagram of minimal surface array structure, boundary conditions, and mesh division

求解结果如图4表1所示。由图4可知,固定网格法与Abaqus求解所得位移云图分布及趋势基本相同。由表1可知,Abaqus模型在网格细化的情况下,位移最大值几乎不再变化,即网格已收敛,因此可将0.2 mm网格的Abaqus计算结果作为标准解用于精度对标。与标准解相比,固定网格法的最大位移误差仅相差2.38%,因此固定网格法求解结果在位移数值和云图分布方面均满足精度要求。前处理效率方面,如表1图5所示,固定网格法无须人工手动划分贴体网格,程序可自动对完整阵列模型划分规则的背景网格,因此网格划分用时仅1.97 s。Abaqus则需要人工手动进行网格划分,仅以简单的四面体单元为例,首先划分单胞网格,然后经过网格复制、对称、连接性检查等操作完成整体结构划分。该部分用时与人员操作熟练度相关,图5展示了本文作者实测采用平均尺寸为0.5 mm的四面体网格,在熟练操作的前提下,完成超材料阵列模型网格划分的完整流程,用时共计9 min。如果采用0.2 mm的网格,由于网格数量大幅增加,Abaqus进行网格清理、复制、对称、连接性检查等操作所花费时间同样大幅增加,其完整流程用时会增加至30 min左右。而众所周知,对于复杂几何模型进行高质量的六面体网格划分则需要更长时间。综上所述,在满足计算精度的前提下,固定网格法可大幅降低用户前处理所需时间。

图4     固定网格法与Abaqus求解所得位移对比
Fig.4     Comparison of displacement solved by fixed grid method and Abaqus
表1     固定网格模型与Abaqus贴体网格模型的网格生成及分析效率对比
Table 1     Comparison of mesh generation and analysis efficiency between fixed grid model and Abaqus body fitted mesh model
特征 单元类型 单元平均尺寸/mm 单元数量 节点数量 合位移最大值/mm 网格划分用时/s 有限元求解用时/s 网格划分与有限元求解总用时/s
固定网格法 二阶六面体 0.5 111926 621885 1.419e-3 1.97 129 130.97
Abaqus 二阶四面体 0.5 655735 1052415 1.383e-3 540 105 645
二阶四面体 0.2 4665056 7058888 1.386e-3 1800 2195 3995

图5     贴体网格与固定网格划分流程和用时对比
Fig.5     Comparison of the process and time consumption for Abaqus body fitted and fixedgrid mesh generation

1.2     基于固定网格法的结构传力–散热多功能一体化拓扑优化

本文考虑体积约束下的结构静力与传热多工况优化问题,基于固定网格法的拓扑优化数学列式为

findx={x1,x2,,xn}
(11)

式中,xi表示第i个子域的设计变量(即伪密度),而不是第i个胞元的设计变量,这是与传统变密度拓扑优化不同之处,即在结构分析计算量不变的情况下,可大幅增加设计变量的数量,从而提高拓扑优化的分辨率;n为子域的数量;w1w2分别为给定的静力工况和传热工况的权系数;UTKUTTKtT分别为结构静柔度与热柔度;U0TK0U0T0TKt0T分别为结构优化前的静柔度和热柔度;目标函数定义为无量纲化后的静柔度和热柔度的加权组合;V为结构初始总体积;Vi为第i个子域的体积;vf为给定的体积分数。应用固体各向同性材料惩罚模型(Solid isotropic material with penalization,SIMP),第i个子域的弹性模量及导热系数与设计变量的插值模型为

Ei=xi3E0+(1xi3)Eminλi=xi3λ0+(1xi3)λminEmin=103E0λmin=103λ0
(12)

式中,E0λ0分别为实体材料的弹性模量和导热系数;E minλ min分别为空材料的弹性模量和导热系数。为了避免数值奇异,E minλ min均为较小的正数而非精确等于0。

为了避免拓扑优化过程中出现棋盘格现象,本文应使用密度过滤方法,对设计变量进行过滤。第i个子域经过过滤后的设计变量为

x˜i=j=1nHi(Xj)Vjxjj=1nHi(Xj)Vj
(13)

式中,Vj为第j个子域的体积;Xj为第j个子域的中心坐标;Hi为权函数,其表达式为

Hi(Xj)=max(rmin|XjXi|,0)
(14)

式中,r min为用户自定义的过滤半径。

此外,为了得到0–1二值化的拓扑优化结果,避免出现灰度单元,使用平滑的Heaviside函数对过滤后的设计变量进行投影。经过投影后的设计变量可通过式(15)计算。

x˜¯i=tanh(θη)+tanh(θ(x˜iη))tanh(θη)+tanh(θ(1η))
(15)

式中,η为投影阈值,本文取值0.5,即高于0.5的设计变量将被投影至趋近于1,低于0.5的设计变量将被投影至趋近于0;θ为用于控制Heaviside函数在η附近的斜率。

设计变量经过过滤和投影操作后,SIMP材料插值模型将被修改为

Ei=x˜¯i3E0+(1x˜¯i3)Eminλi=x˜¯i3λ0+(1x˜¯i3)λmin
(16)

相应的,目标函数对第k个设计变量的灵敏度可通过链式法则计算。

fxk=i=1nfx˜¯ix˜¯ix˜ix˜ixk
(17)

式(17)中等号右端第1项可通过伴随变量法求解,后两项可由公式和推导得到。

x˜ixk=Hi(Xk)Vkj=1nHi(Xj)Vj
(18)

x˜¯ix˜i=θ(1tanh2(θ(x˜iη)))tanh(θη)+tanh(θ(1η)) 
(19)

综上所述,基于固定网格法的结构拓扑优化流程如图6所示,具体步骤为:

图6     基于固定网格法的结构拓扑优化流程
Fig.6     Structural topology optimization process based on fixed grid method

(1)将结构置于更大的、规则的扩展域,并将扩展域离散为规则网格;

(2)利用四叉树/八叉树对物理胞元和边界胞元进行细分;

(3)初始化子域设计变量。仅物理胞元和边界胞元子域有设计变量,虚拟胞元不参与有限元分析和优化过程,所有子域密度初始化为vf

(4)对设计变量进行过滤和投影操作;

(5)采用SIMP模型求解有限元方程得到位移场和温度场;

(6)更新目标函数值和约束值;

(7)灵敏度分析;

(8)判断结果是否收敛。如果不收敛,利用移动渐近线法[  SVANBERG K. The method of moving asymptotes—A new method for structural optimization[J]. International Journal for Numerical Methods in Engineering, 1987, 24(2): 359–373.
27
]
更新设计变量,然后转到步骤(4),进行新一轮的优化设计。如果收敛,则结束流程,将优化后的结构边界光滑化[  REN H Y, XIA B, WANG W R, et al. AMRTO: Automated CAD model reconstruction of topology optimization result[J]. Computer Methods in Applied Mechanics and Engineering, 2025, 435: 117673.
28
]
并输出目标函数值。

2     极小曲面复合加筋超材料模型构造与优化工况说明

应用上述方法,面向传力与散热一体化多功能超材料设计,本文对基于TPMS的超材料单胞进行拓扑优化。TPMS是极小曲面中的一类可周期性阵列的曲面,该类曲面满足以下方程。

i=12xi(fxi1+(fx1)2+(fx2)2)=0
(20)

上式可解出不同几何构型的TPMS解析式。本文选择具有体心对称的Schwarz’s Primitive曲面,其表达式为

cos(2πx15)+cos(2πy15)+cos(2πz15)=0
(21)

为了进一步扩大设计空间,本文在原始Primitive单胞的内部添加十字加筋结构,如图7(a)所示。为了降低计算量且保持结构对称性,选取单胞的八分之一作为拓扑优化的初始设计域。八分之一单胞在固定网格下的计算模型如图7(b)所示(为了便于读者查看,示意图胞元尺寸采用1 mm,实际计算模型胞元尺寸为0.2 mm)。静力与传热工况分别如图7(c)和(d)所示。静力工况中,在顶部施加沿Z轴负方向的均布1 MPa压力,在结构对称面上施加位移对称边界条件,即XY平面上Z方向位移为0;YZ平面上X方向位移为0;ZX平面上Y方向位移为0。传热工况中,对称面上施加绝热边界,结构内部施加105 W/m3的均匀体热源,顶部施加23 ℃的固定温度热沉,用以模拟结构内部生热,两端散热工况。单胞材料为铝合金,弹性模量为71000 MPa,泊松比为0.3,导热系数为173 W/(m·K)。胞元尺寸为0.2 mm,采用一阶六面体单元及2层八叉树细分,物理胞元与边界胞元总数为50653,结构总自由度数为72729,子域/设计变量总数1036658。设计变量数量远高于胞元数量,验证了固定网格法可在结构分析自由度较少的情况下,实现较高分辨率的拓扑优化。体积分数vf设置为0.5,即优化后的结构体积不超过初始设计域的0.5倍。

图7     TPMS超材料单胞初始设计域及工况设置
Fig.7     Initial design domain and load conditions of TPMS metamaterial unit cell

3     优化结果及讨论

不同权系数下的拓扑优化结果如图8表2所示。其中,权系数组合(1,0)和(0,1)分别为静力单工况和传热单工况。静力单工况下,优化构型呈现出显著的载荷定向聚集效应,使得主承载方向(Z方向)的十字加筋结构得以完整保留。并且为了有效抑制X/Y方向的位移分量,底部保留了呈对称分布的斜交加强肋,整体呈现梁、板、壳形式的综合体结构,与文献[  LIU Y, WANG Y Z, REN H Y, et al. Ultrastiff metamaterials generated through a multilayer strategy and topology optimization[J]. Nature Communications, 2024, 15(1): 2984.
25
]所得结论一致。传热单工况下,结构呈现出大量不规则的分叉形式,通过构建大量导热通路,可实现更强的导热性能。由于不考虑静力性能,因此破坏了构型的十字加筋结构,导致其静柔度值大幅增加。多工况优化结果既保留了十字加筋结构以保证结构刚度,同时也具有大量导热通路以保证导热性能。并且随着静力权系数w1的降低,传热权系数w2的增加,静柔度逐渐增大,热柔度逐渐降低。此外,由图8可观察到,单工况优化结果虽然在设计工况下可以获得最佳性能,但在非设计工况下性能呈现显著劣化(静柔度增幅1350%,热柔度增幅370%),难以适应复杂应用场景下的多功能设计要求。相较而言,多工况协同优化方法通过有限度的性能折损(静柔度与热柔度增幅均在60%以内),实现了多工况条件下结构响应的均衡性。

图8     不同权系数下优化所得结构静柔度与热柔度
Fig.8     Static compliance and thermal compliance of optimized structures fordifferent weight coefficients
表2     不同权系数下的拓扑优化结果
Table 2     Topology optimization results of different weight coefficients
权系数组合 (w1w2 优化构型 静柔度/10–3 mJ 热柔度/(mW·K)
等轴测视图 右视图 底视图
(1,0) 1.31 2.54
(0.9,0.1) 1.46 0.86
(0.8,0.2) 1.49 0.78
(0.6,0.4) 1.59 0.74
(0.5,0.5) 1.67 0.71
(0.4,0.6) 1.70 0.69
(0.2,0.8) 1.96 0.58
(0.1,0.9) 2.10 0.56
(0,1) 19.0 0.54

以静力单工况优化结果为例(对应表2中权系数(1,0)),经过3次对称操作,再进行4×4阵列优化后,可得到如图9所示的全尺寸超材料阵列模型,该超材料可作为抗压层合夹心结构的夹心层使用。

图9     优化后的全尺寸超材料阵列模型
Fig.9     Optimized full-size metamaterial array model

为了验证优化结果的正确性,对权系数组合(1,0)所得优化结果重新划分了足够细密的二阶四面体贴体网格(图10(a)),并分别采用OptiStruct 2021软件和Abaqus 2023软件进行了仿真验证。同时采用固定网格法对优化结果进行仿真分析,进一步验证固定网格法分析的准确性,结果如图10表3所示。两款软件采用了相同的网格划分及边界条件,计算结果几乎完全相同(仅第5位有效数字开始不同)。以Abaqus计算结果为基准,固定网格法优化所得结果以及重新仿真验证的静柔度和热柔度误差分别低于1.11%和3.10%,且在工程误差许可范围内,从而验证了固定网格法优化与分析的正确性。

图10     权系数组合(1,0)对应优化结果的重新仿真验证
Fig.10     Re-simulation verification corresponding to the optimization result of weight coefficient combination(1,0)
表3     权系数组合(1,0)对应优化结果的模型重建参数与仿真验证
Table 3     Model reconstruction parameters and simulation verification of the optimization result corresponding to theweight coefficient combination (1,0)
参数 OptiStruct Abaqus 固定网格法优化结果 固定网格法重新验证
单元类型 二阶四面体 二阶四面体 一阶六面体 二阶六面体
单元平均尺寸/mm 0.05 0.05 0.2 0.2
单元数量 1044320 1044320 50653 11061
节点数量 1544374 1544374 24243 57254
静柔度/10–3 mJ 1.295758 1.295672 1.31 1.306253
热柔度/(mW·K) 2.620584 2.620721 2.54 2.577459

本文所研发的基于固定网格技术的结构分析和拓扑优化方法已在国产全自主CAX工业软件OptFuture中实现了落地。此外,面向工程中的多级减重设计需求,OptFuture软件还实现了多级晶格设计填充,以图9中的静力工况优化结果为例,将其作为一级晶格设计构型,在其基础上可以进一步进行二级晶格设计填充,次级填充晶格可自动适应上一级优化设计结果的非规则外形。如图11(a)和(b)所示,作为示例,本文给出了Gyroid型TPMS薄壁结构填充以及体心立方型杆类晶格填充的结果(为了便于查看,仅展示八分之一模型的填充结果)。

图11     优化结果的二级晶格填充
Fig.11     Secondary-level lattice filling of optimized result

4     结论

(1)实现了基于固定网格法的结构分析及拓扑优化方法。无须手动划分贴体的有限元网格,从而大幅降低前处理所需时间。结合基于四叉树/八叉树的局部自适应加密技术,可实现高效率、高精度的复杂结构有限元分析。以固定网格的子域密度为设计变量,在不增加有限元分析计算量的前提下提高设计变量的数量,从而在兼顾结构分析精度和效率的前提下,实现更高分辨率的拓扑优化。

(2)面向承重与散热一体化需求,对TPMS复合加筋单胞进行了静力工况、传热工况、力热多工况拓扑优化。单工况优化中,优化结果虽然能在目标工况下获得最佳性能,但在非设计工况下性能呈现显著劣化,难以适应复杂应用场景下的多功能设计要求。应用多工况优化,通过有限度的性能折损(相较单工况优化获得的最佳性能,静柔度与热柔度增幅均在60%以内),提高了结构响应的均衡性,拓展了超材料的多功能应用潜力。

(3)针对超材料结构多级轻量化的工程需求,依托国产全自主CAX工业软件OptFuture,实现了多级晶格设计填充,填充结构可自动适应拓扑优化结果的非规则外形,进一步拓展了超材料的轻量化设计潜力。

作者介绍



陈雪乾 博士研究生,研究方向为结构拓扑优化和超材料设计。

参考文献

[1]

胡更开, 刘晓宁. 弹性超材料设计与波动控制[M]. 北京: 科学出版社, 2024.
HU Gengkai, LIU Xiaoning. Design and fluctuation control of elastic metamaterials[M]. Beijing: Science Press, 2024.

[2]

EVANS A G, HUTCHINSON J W, FLECK N A, et al. The topological design of multifunctional cellular metals[J]. Progress in Materials Science, 2001, 46(3–4): 309327.

[3]

程乾, 尹剑飞, 温激鸿, . 极小曲面力学超材料抗冲吸能特性分析[J]. 动力学与控制学报, 2023, 21(7): 4350.
CHENG Qian, YIN Jianfei, WEN Jihong, et al. Impact resistance and energy absorption of mechanical metamaterials with minimal surfaces[J]. Journal of Dynamics and Control, 2023, 21(7): 4350.

[4]

LI X W, YU X, ZHAI W. Additively manufactured deformation-recoverable and broadband sound-absorbing microlattice inspired by the concept of traditional perforated panels[J]. Advanced Materials, 2021, 33(44): 2104552.

[5]

HAGHPANAH B, SALARI-SHARIF L, POURRAJAB P, et al. Multistable shape-reconfigurable architected materials[J]. Advanced Materials, 2016, 28(36): 79157920.

[6]

FRENZEL T, KADIC M, WEGENER M. Three-dimensional mechanical metamaterials with a twist[J]. Science, 2017, 358(6366): 10721074.

[7]

REN X, DAS R, TRAN P, et al. Auxetic metamaterials and structures: A review[J]. Smart Materials and Structures, 2018, 27(2): 023001.

[8]

PERALTA I, FACHINOTTI V D, ÁLVAREZ HOSTOS J C. A brief review on thermal metamaterials for cloaking and heat flux manipulation[J]. Advanced Engineering Materials, 2020, 22(2): 1901034.

[9]

LI Y, LI W, HAN T C, et al. Transforming heat transfer with thermal metamaterials and devices[J]. Nature Reviews Materials, 2021, 6: 488507.

[10]

CHEN M J, PEI Y M, FANG D N. Design, fabrication, and characterization of lightweight and broadband microwave absorbing structure reinforced by two dimensional composite lattice[J]. Applied Physics A, 2012, 108(1): 7580.

[11]

ZIGONEANU L, POPA B I, CUMMER S A. Three-dimensional broadband omnidirectional acoustic ground cloak[J]. Nature Materials, 2014, 13(4): 352355.

[12]

BENDSØE M P, SIGMUND O. Material interpolation schemes in topology optimization[J]. Archive of Applied Mechanics, 1999, 69(9): 635654.

[13]

SIGMUND O. A 99 line topology optimization code written in Matlab[J]. Structural and Multidisciplinary Optimization, 2001, 21(2): 120127.

[14]

SETHIAN J A, WIEGMANN A. Structural boundary design via level set and immersed interface methods[J]. Journal of Computational Physics, 2000, 163(2): 489528.

[15]

WANG M Y, WANG X M, GUO D M. A level set method for structural topology optimization[J]. Computer Methods in Applied Mechanics and Engineering, 2003, 192(1–2): 227246.

[16]

GUO X, ZHANG W S, ZHONG W L. Doing topology optimization explicitly and geometrically—A new moving morphable components based framework[J]. Journal of Applied Mechanics, 2014, 81(8): 081009.

[17]

ZHANG W S, YUAN J, ZHANG J, et al. A new topology optimization approach based on moving morphable components (MMC) and the ersatz material model[J]. Structural and Multidisciplinary Optimization, 2016, 53(6): 12431260.

[18]

ZHANG P, Biligetu, QI D X, et al. Mechanical design and energy absorption of 3D novel hybrid lattice metamaterials[J]. Science China Technological Sciences, 2021, 64(10): 22202228.

[19]

WEI Y L, YANG Q S, LIU X, et al. Multi-bionic mechanical metamaterials: A composite of FCC lattice and bone structures[J]. International Journal of Mechanical Sciences, 2022, 213: 106857.

[20]

LI D M, QIN R X, XU J X, et al. Improving mechanical properties and energy absorption of additive manufacturing lattice structure by struts’ node strengthening[J]. Acta Mechanica Solida Sinica, 2022, 35(6): 10041020.

[21]

MESSNER M C. Optimal lattice-structured materials[J]. Journal of the Mechanics and Physics of Solids, 2016, 96: 162183.

[22]

BERGER J B, WADLEY H G, MCMEEKING R M. Mechanical metamaterials at the theoretical limit of isotropic elastic stiffness[J]. Nature, 2017, 543(7646): 533537.

[23]

LONG K, DU X R, XU S Q, et al. Maximizing the effective Young’s modulus of a composite material by exploiting the Poisson effect[J]. Composite Structures, 2016, 153: 593600.

[24]

WANG Y Q, LUO Z, ZHANG N, et al. Topological shape optimization of microstructural metamaterials using a level set method[J]. Computational Materials Science, 2014, 87: 178186.

[25]

LIU Y, WANG Y Z, REN H Y, et al. Ultrastiff metamaterials generated through a multilayer strategy and topology optimization[J]. Nature Communications, 2024, 15(1): 2984.

[26]

DÜSTER A, PARVIZIAN J, YANG Z, et al. The finite cell method for three-dimensional problems of solid mechanics[J]. Computer Methods in Applied Mechanics and Engineering, 2008, 197(45–48): 37683782.

[27]

SVANBERG K. The method of moving asymptotes—A new method for structural optimization[J]. International Journal for Numerical Methods in Engineering, 1987, 24(2): 359373.

[28]

REN H Y, XIA B, WANG W R, et al. AMRTO: Automated CAD model reconstruction of topology optimization result[J]. Computer Methods in Applied Mechanics and Engineering, 2025, 435: 117673.

目录