周期性功能梯度结构稳态热传导拓扑优化设计

李信卿 赵清海， 张洪信 张铁柱 陈建良

1.青岛大学机电工程学院，青岛，2660712.青岛大学动力集成及储能系统工程技术中心，青岛，266071

0 引言

功能梯度结构作为一种新型非均匀结构，具有缓冲吸振、隔热降噪等优异的综合性能，已得到航空航天、轨道交通、通信电子等行业的广泛关注，近年来逐渐成为跨结构、材料领域的研究热点。自20世纪80年代BENDSØE和KIKUCHI[1]首次提出连续体结构拓扑优化概念以来，拓扑优化方法已发展出均匀化法、渐近结构优化法[2]、变密度法[3]、独立连续映射法[4]、水平集法[5]和移动变形组件法[6-7]等，各类方法相辅相成推动了宏/微观结构的协同发展。PAULINO等[8]率先通过准则法搭建了功能梯度材料的固体各向同性微结构惩罚(solid isotropic microstructure with penalization，SIMP)模型，探讨了功能梯度结构拓扑优化设计的有效性。此后，功能梯度宏观结构拓扑优化不断得到发展和完善。邱克鹏等[9]采用凸规划求解策略，获得了功能梯度夹层结构中夹芯的拓扑构型。杨旭静等[10]结合节点梯度位移函数，建立了功能梯度材料结构拓扑优化近似显式模型。CHENG等[11]提出了一种应力约束下基于结构拓扑优化设计的梯度点阵结构新方法，以生成具有可预测屈服性能的轻量化点阵结构。上述研究拓展了功能梯度结构拓扑优化方法，但未考虑微观材料的结构、布置对宏观结构的影响，极大限制了微观材料性能的充分利用。

鉴于此，部分学者在宏观结构拓扑优化中施加周期性约束，为发掘微观材料结构提供了新思路。SIGMUND[12-13]运用逆均匀化技术实现了具有极值力学性能(如零剪切性能、负泊松比)的周期性材料微结构设计。除均匀化法外，焦洪宇等[14]提出了一种基于变密度理论SIMP模型的周期性拓扑优化方法。杜义贤等[15]以周期性单胞为研究对象，构建出了基于宏观力学性能的细观点阵结构优化模型。付志方等[16]通过搭建线弹性结构周期性稳健拓扑优化模型，得到了具有良好稳定性的周期性结构。此外，周期性设置在热传导结构上的应用已逐渐成为新的研究热点。龙凯等[17]通过建立周期性结构拓扑优化模型，实现了传热微观复合材料优化设计。贾娇等[18]基于SIMP周期性结构模型，验证了具有宏观导热性能的周期性材料的设计可行性问题。赵清海等[19]考虑周期性约束下的多材料结构稳态热传导拓扑优化设计模型，探讨了一种周期性多材料结构拓扑设计方法。上述文献探讨了宏观结构中最优微观材料结构，但并未考虑宏/微观结构之间的最佳性能匹配设计，往往难以实现结构的综合性能最优。

在此背景下，周期性功能梯度结构拓扑优化设计应运而生，它可综合考虑宏/微观结构的性能匹配。ZHOU等[20]借助两相梯度微结构的逆均匀化方法，确定了周期性基底单元沿平行于梯度的变化方向。DENG等[21]搭建并行优化模型应用到轻量化结构拓扑和材料微观结构设计，探究了功能梯度结构的多目标设计问题。GUO等[22]研究了多尺度框架下有界载荷不确定性材料和结构的鲁棒并行优化问题。此类设计方法主要针对由单一材料微结构组成的宏观结构设计，研究具有多层微结构的功能梯度拓扑优化时，则根据材料性能需要在不同区域分布不同的微结构，这可增加结构的多样性。张卫红等[23]针对功能梯度材料，给出了与尺寸关联的结构和材料拓扑优化设计结果。COELHO等[24]基于层次计算方法探讨了特定局部微结构的等效材料特性和层次优化设计问题，并实现了三维功能梯度结构设计。XIA等[25]实现了功能梯度结构材料性能和拓扑布局的并行优化设计。HUANG等[26]采用双向进化结构优化方法对材料单元进行灵敏度分析，得到了多孔材料和复合材料的各向异性微结构。LI等[27]借助蜂窝复合材料和功能梯度材料的优势，提出了一种用于集成材料和结构设计的分层多尺度的拓扑优化方法。上述研究集中在结构场下力学性能对宏观结构和微观构型的影响，有关热传导功能梯度方面的研究报道较少。相对于力场问题，受热载荷影响热场结构的拓扑优化问题则更加复杂[28]。在散热结构设计中，应综合考虑周期性与功能梯度约束，实现装配简单、生产成本低、易于模块化设计等周期性结构制造加工的同时，并通过对结构进行梯度分层来保证结构具有良好的散热性能。

本文提出了一种周期性约束下的功能梯度结构稳态热传导拓扑优化方法。采用基于SIMP模型的材料插值方法，利用热传导结构提取最优构型中各预设梯度层的体积分数，通过重新分配散热弱度来实现周期性约束设置，同时采用移动渐近线法推导出设计变量的迭代公式，借助2D和3D数值算例，探讨周期性约束下功能梯度稳态热传导结构的材料分布、各区域的宏微观结构以及优化结果。

1 热传导周期性功能梯度结构描述

对于热传导周期性功能梯度结构，宏观设计域被划分为多个梯度层(L1，…，Lj，…，Ln)，其中j(j=1，2，…，n)为梯度层序号，n为划分的总梯度层数，每个梯度层由周期性材料微结构组成。首先给定热传导结构的Dirichlet边界ГD、Neumann边界ГN和内热源等初始条件，对宏观结构进行拓扑优化设计，获得各梯度层材料密度最优结果；进而对梯度层宏观结构和微观材料协同优化设计，获得周期性功能梯度材料微结构最佳分布，实现功能梯度拓扑优化设计。热传导周期性功能梯度结构如图1所示。

图1 热传导周期性功能梯度结构Fig.1 Heat conduction periodic functional gradient structure

2 周期性功能梯度结构的建模方法

2.1 基于SIMP的传统热传导结构建模

针对热传导结构设计，基于传统SIMP法建立的相对密度插值模型为

(1)

式中，xi为单元i的相对密度(即设计变量)，xi∈[0，1]；λ为插值后的热导率；λ0、λv分别为材料初始热导率和孔洞区域热导率；p为惩罚因子。

建立稳态热传导的结构拓扑优化数学模型：

(2)

式中，C为结构散热弱度；Q、T和Kλ分别为热载荷向量、节点温度向量和整体热传导刚度矩阵；Ne为设计域单元总数；vi为单元i的体积；V、f、V0分别为优化后体积、设计域的体积比和总体积；xmin为单元相对密度的下限值。

2.2 功能梯度结构建模和求解

为计算体积约束下功能梯度结构的总体布局，借助热传导结构拓扑优化模型，可提取最优拓扑构型中各预设梯度层的体积分数。将二维矩形设计域划分为n个梯度层(L1，L2，…，Lj，…，Ln)，分层梯度设置如图2所示。

图2 分层梯度设置示意图Fig 2 Layered gradient setting diagram

针对功能梯度结构设计，建立SIMP插值模型，可表示为

(3)

式中，xi，j为第j梯度层单元i的相对密度，xi，j∈[0，1]。

其拓扑优化数学模型可表示为

(4)

式中，N为梯度层内的单元数；ti，j、kλ分别为第j梯度层单元i的节点温度向量和单元热传导刚度矩阵；vi，j为第j梯度层单元i的体积；fj、Vj分别为第j梯度层的体积比和总体积。

热传导结构拓扑优化设计灵敏度求解过程中，利用伴随变量法推导出目标函数C相对于单元密度xi，j的灵敏度，即

(5)

不考虑温度载荷的相关性，则∂QT/∂xi，j=0，根据平衡方程KλT=Q对设计变量xi，j求导，可得

(6)

将式(6)代入式(5)可得

(7)

体积约束V相对于单元密度xi，j的导数为

(8)

2.3 设计变量更新和灵敏度过滤

目前常用的优化求解算法包括数学规划法和优化准则法，由于移动渐近线法(method of moving asymptotes，MMA)[29]对复杂目标和多约束拓扑优化问题具有更好的适定性，因此，本文选取MMA法进行功能梯度结构拓扑优化设计求解。

对于拓扑优化结果中存在的数值不稳定现象，目前常见的解决方法包括周长约束法、过滤法等。基于Helmholts方程的偏微分方程(partial differential equation，PDE)过滤方法具有高效的计算效率，适合用于大规模的优化求解问题。本文借助PDE过滤方法进行灵敏度过滤，其过滤方程组为

(9)

(10)

(11)

2.4 周期性结构设置

为获得周期性功能梯度结构拓扑优化形式，可将每个梯度层划分为Mxj×Myj个相同的子区域(其中Mxj、Myj分别为第j梯度层在x与y轴方向的子区域数)，使得各梯度层不同子结构在相同位置处的单元具有相同的材料属性。如图3所示，梯度层内周期性结构的数学模型可描述为

(12)

图3 周期性分层设置示意图Fig.3 Schematic diagram of periodic layered setup

(13)

图4 周期性功能梯度拓扑优化流程图Fig.4 Flowchart of periodic functional gradient structure

3 算例

3.1 2D算例

算例1 以二维热传导问题为例，宏观设计域如图5a所示，选取几何尺寸为180 mm×180 mm 的正方形结构，将其离散为180×180个四节点四边形单元。在整个设计域施加强度为Q=0.01 W/mm2的均布热源，左边界中间位置给定温度为零，其余边界绝热，宏观设计域的体积分数为0.55，惩罚因子p=3，过滤半径rmin=1.2 mm。分层设计域如图5b所示，沿x轴方向分别对设计域进行梯度预设，给定梯度层数n=3，5，10，由宏观拓扑优化可获得三种预设梯度层的体积分数如表1所示。

(a)宏观设计域 (b)分层设计域图5 二维设计域Fig.5 Two-dimensional design domain

表1 各梯度层的体积分数Tab.1 Volume fraction of gradient layers

全局周期性约束下的拓扑构型如图6所示。周期性约束为My=3，5，10；Mx=1，以此进行周期性功能梯度拓扑优化研究。功能梯度结构的拓扑构型以及温度分布分别如图7和图8所示，并将梯度层子结构的优化结果列于表2。

由图6可知，对于全局周期性约束下的拓扑构型，随子区域数量的增加，其微观子区域的细节特征趋于简洁，宏观结构则呈现类似“框架”的形式。由图7可知，拓扑构型呈现一种由绝热边界向零温点边界聚拢的“树枝”状分布，且在周期性约束相同的情况下，随梯度层数的增加，结构的拓扑构型差异愈发明显。由表2可以发现，每个梯度层的子区域均呈现出材料最优化分布。图8中，由蓝色到黄色的温度分布表示温度由低到高。由图8可知，零温处的温度最低，离零温处越远其温度呈梯度升高，具有一定的周期性特征，且温度分布较为均匀。

(a)My=Mx =1 (b)My=Mx =3

(a)My=3；Mx=1

(a)My=3；Mx=1

图9为全局周期(图6c)和周期性分层梯度(图7a～图7c中梯度层数划分为n=5时)设置下散热弱度随迭代次数的变化曲线，可以发现分层梯度会对结构的散热效果产生影响，在周期性约束相同的情况下，结构的散热性能随梯度层数的增加而降低；而对比全局周期以及周期性分层梯度设置下结构的散热性能发现，后者的散热性能明显优于前者的散热性能，说明周期性功能梯度结构具有优越性。由图9还可知，分层梯度设置下目标函数曲线经过迭代更新可收敛到最优解，且收敛速度较快，验证了所提方法的有效性。

表2 2D拓扑构型中梯度层内子结构Tab.2 Gradient layer sub-structure in 2D topological configurations

图9 目标性能随迭代次数的变化曲线Fig.9 Variation curve of target performance with number of iterations

算例2 选取几何尺寸为180 mm×180 mm的正方形结构，分别离散为180×180个四节点四边形单元，将强度为Q=324 W/mm2的内热源施加在设计区域(45 mm，45 mm)、(1 mm，90 mm)和(90 mm，90 mm)位置点，四个角点的温度恒为零，四边界绝热宏观设计域的体积分数为0.55，惩罚因子p=3，过滤半径rmin=1.2 mm。

宏观结构设计域如图10所示，沿x轴方向分别对设计域进行梯度预设，给定梯度层数n=1，3，6，由宏观拓扑优化可获得三种预设梯度层的体积分数。周期性约束为My=3；Mx=1，不同热源位置下功能梯度结构优化构型结果及其温度分布分别如图11和图12所示，功能梯度结构的散热性能对比见表3。

(a)(45，45)mm (b)(1，90)mm (c)(90，90)mm图10 结构设计域Fig.10 Structural design domain

(a)(45 mm，45 mm)位置

(a)(45 mm，45 mm)位置

由图11可知，不同受热位置下热传导结构的拓扑构型有明显差异。由表3可以发现，同一分层梯度设置下(即梯度层数相同时)，不同热源位置结构的散热性能会有较大差异，这说明热源位置可能会影响周期性功能梯度结构的最短散热路径。由图12可知，温度最高点出现在热源位置，整个设计区域的温度分布较为均匀，整体温度梯度较小，内热源附近的温度梯度较大。

3.2 3D算例

以三维热传导问题为例，如图13所示，选取几何尺寸为120 mm×6 mm×120 mm的长方体结构，将其离散为120×6×120个八节点六面体单元，三维结构下端面中间位置边界温度设为零，其余边界表面绝热。将强度为Q=0.01 W/mm2的均布热源施加在整个三维区域，体积分数约束为0.55，惩罚因子p=3，过滤半径rmin=1.4 mm，沿z轴方向分别对设计域进行梯度预设，给定梯度层数n=3，6，由宏观拓扑优化可获得两种预设梯度层的体积分数。周期性约束为Mx=3，6；My=Mz=1，功能梯度结构拓扑构型如图14所示，图14b中梯度层数划分为n=3时梯度层内的子结构如图15所示。

表3 功能梯度结构散热性能对比Tab.3 Comparison of heat dissipation performance of functional gradient structures

图13 三维宏观设计域和结构拓扑构型Fig.13 Three-dimensional macro design domain and macro-structural topology configurations

(a)My=Mz=1；Mx=3

图15 3D拓扑构型中梯度层内子结构Fig.15 Gradient layer sub-structure in 3D topological configurations

由图14可以看出，周期性功能梯度结构子区域构型呈现出了与图13相似的“树枝状”拓扑构型。在周期性约束相同的情况下，随着梯度层数的增加，结构拓扑构型越发明显且边界清晰合理，散热性能逐渐降低，从而证明了所提方法可实现热传导结构周期性功能梯度拓扑优化设计，保证了该方法在工程上的适用性。

以图14a为例，三维热传导结构的散热曲线随迭代次数变化情况见图16，可以看出，当迭代到60步后，变化曲线趋于平稳。不同周期性约束下功能梯度热传导结构的散热弱度变化不大，与二维设计相对应，进一步验证了功能梯度结构具有良好的散热性能。

图16 目标性能随迭代次数的变化曲线Fig.16 Variation curve of target performance with number of iterations

4 结论

(1)提出了一种考虑周期性约束的功能梯度结构稳态热传导拓扑优化设计方法，对比探讨了全局周期以及周期性功能梯度拓扑优化，通过2D与3D数值算例验证了所提方法的有效性。

(2)通过对宏观设计区域进行梯度分层设置，并结合固体各向同性微结构惩罚(SIMP)插值方法和重新分配散热弱度基值，可获得功能梯度结构周期性热传导构型。

(3)数值算例结果表明，各梯度层不同设计变量设置下，结构在宏微观构型方面具有差异性，表明了宏观性能约束下，梯度层微观结构的改变会对宏观拓扑构型产生影响。

(4)采用所提方法获得的周期性功能梯度结构具有构型材料分布清晰、方法设置简便的优点，且有广阔的工程适用前景。