基于密度法的传热结构拓扑优化设计

崔天福， 丁晓红， 侯丽园

（上海理工大学 机械工程学院，上海 200093）

具有散热通道的传热结构的传热性能与散热通道的材料、形态分布和外部热环境密切相关.当材料和外部热环境一定时，传热结构中的散热通道就对整个结构的热量传递起到主要作用，因此，如何在传热结构中合理布置散热通道对提高结构的传热性能、延长其使用寿命至关重要.在发热量比较集中、空间有限的集成电路中，可采用由高导热材料构成的散热通道直接嵌入到元器件内部，将热量导出到外部环境中，再采用其它高效冷却方式进行冷却，可以有效解决空间限制问题和冷却效率问题.传统的传热结构设计一般基于传热学基本理论及实践经验，然而这种方法很难解决复杂传热结构的设计问题.针对这种复杂模型，可使用结构拓扑优化方法，计算出高自由度的散热通道拓扑形态，在此形态的基础上，通过经验改进以及进一步的尺寸优化获得最佳的结构尺寸.

结构拓扑优化以材料分布为优化对象，在设计空间内寻求最佳的分布方案，以得到结构某种性能最优.在拓扑优化方法中，理论相对成熟且有一定工程应用的方 法 主 要 有 密 度 法［1］、均 匀 化 方 法［2－3］及Level Set方法［4－5］等.近年来这些方法在传热结构拓扑优化设计中已有所应用.如Bendsøe等［6］将结构拓扑优化中的密度法直接拓展到了传热结构优化当中，并对简单的稳态热传导问题进行了研究.Iga［7］在均匀化理论的基础上建立了以结构总势能为目标函数的优化数学模型，确定了结构具有最佳传热效应时，在Dirichlet边界、Neumann边界、对流边界情况下目标 函 数 的 最 大 值、最 小 值 问 题.Yamada［8－9］以Lveleset方法为基础建立了热扩散最大为目标的传热拓扑优化数学模型，并对热力耦合情况下考虑热扩散效果的结构刚度最大化问题进行了研究.对于体－点热 传 导 问 题，Bejan［10］提 出 构 形 理 论（constructal theory）进行求解，该方法首先对初始单元进行优化，获得高导热通道的尺寸参数，然后对由初始单元组装而成的一级装配体进行优化，获得高导热通道的尺寸参数，逐级装配直至将高导热材料组成的散热通道覆盖至所需要的散热面积.Gersborg－Hansen［11］提出了在有限体积法基础上的传热结构拓扑优化的方法.丁晓红［12－13］通过研究植物根系的成长机理，提出自适应成长法来解决薄板结构的加强筋分布的设计问题，并将自适应成长法应用于体－点散热问题，该方法通过对局部分枝的灵敏度分析，使灵敏度大的区域分枝生长迅速，灵敏度小的区域分枝生长缓慢或者退化，从而在迭代过程中实现分枝结构的自适应功能，使得最终分枝的拓扑结构具有最优的形态布局.

本文以变密度法为理论基础，并以SIMP（solid isortopic material with penalization model）插值模型［14－15］为材料插值模型，建立结构散热弱度为目标函数的优化数学模型，推导出满足KKT（Karush－Kuhn－Tucker）条件的传热结构拓扑优化迭代公式.采用高阶单元方法来解决优化过程中存在的数值不稳定现象，然后通过不同边界组合下的数值算例，研究了拓扑形态分布上的一些特性.

1 理论分析与优化模型

1.1 传热结构的拓扑优化数学模型

如图1所示，在设计区域Ωd内有生热率载荷Q，存 在Dirichlet 边 界ΓT，边 界 温 度 为T0.在Neumann边界Γq，热流密度为q.传热结构的拓扑优化问题为：如何在区域Ωd中布置高导热材料，形成散热通道，将区域内的热量传送到边界.

定义 TTP 为结构的散热弱度［16］，其中，T 和P是离散结构中的节点温度矢量和等效节点载荷矢量，等效节点载荷矢量包括生热率载荷和热流边界引起的等效载荷.结构的散热弱度是结构对外部热载荷的响应，由于散热弱度与热势能具有一致性，当散热弱度最小时，在相同的外部热载荷作用下所能传递的热量最大，使结构具有最佳的温度分布［7］.最佳的温度分布体现在结构的平均温度梯度小、平均温度低等方面.因此，可以通过使结构的散热弱度最低来确定热传导问题中高导热材料的最优分布.

图1 二维热传导问题Fig.1 Two－dimensional heat conduction problem

在基于密度法的结构拓扑优化设计中，SIMP插值模型是广泛采用的模型，如图2所示（见下页）.在该插值模型中引入了一个伪密度变量的概念，每一个离散单元都对应一个伪密度变量.对于传热结构拓扑优化问题，可使伪密度变量关联的材料属性为离散单元的热传导系数.图2的SIMP 插值关系可以表示为

式中，xi为离散单元伪密度变量，即优化设计变量；λi表示伪密度变量xi的单元热传导系数；λmax表示传热结构中高导热材料的热传导系数；λmin表示低导热材料的热传导系数；p 为惩罚系数.

引入p 的目的是在优化过程中避免中间密度单元的形成，使优化后的传热结构中尽可能只存在λmax和λmin两种导热系数的材料.当p＜1时，对中间密度具有强化作用；当p＞1时，对中间密度具有惩罚作用；当p＝1时，不对中间密度进行惩罚或者强化.通常p 取2～5范围内的值.

图2 SIMP插值曲线Fig.2 SIMP interpolation curve

为了便于表达数学模型，引入一个高、低导热材料的热传导系数比h＝λmax／λmin.式（1）可以表示为

式中，mi为单元导热比例系数，表示伪密度变量xi的导热系数λi与低导热材料导热系数λmin之间的比值.

由单元热传导系数与单元热传导矩阵的线性关系可推出

式中，kmin是当热传导系数为λmin时的单元热传导矩阵；ki是导热系数为λi时的单元热传导矩阵.

以结构的散热弱度C＝TTP 为目标函数，在SIMP插值模型下，其拓扑优化模型可表示为

式中，ti为第i个单元的节点温度矢量；vi为第i个单元的体积；K 为优化模型的总热传导矩阵；V＊为体积约束上限值；V 为优化后的体积；V0为设计域体积；f 为体积约束百分率；xmin为伪密度下限值；n为离散单元总数.

1.2 基于最优准则法的迭代公式

由上述目标函数、体积约束、平衡方程约束、设计变量约束构造拉格朗日函数，当该构造函数满足KKT 条件时，目标函数C 具有最小值.综合考虑优化设计变量的边界约束以及热传导矩阵的对称性，对于离散单元应满足KKT 条件的表达式为

式中，λl为Lagrangian体积乘子.

由（6）式可得

其中

式中，w 为当前迭代次数；si为第i个单元的热耗散（heat dissipation）；为第w 次迭代过程中第i个单元的跌倒系数.

将式（8）和式（9）代入式（7），整理得

将式（10）作为优化设计准则，得出基于该准则下的迭代公式.

式中，η是为了保证数值计算的稳定性而引入的一个阻尼系数；z为密度变化上限值.

基于最优准则法实现的热传导体结构拓扑优化程序流程如图3所示.首先建立有限元模型，包括几何模型的建立以及材料、热载荷、边界约束的定义.然后初始化设计变量，将材料参数用设计变量来表示.计算优化模型的目标值以及敏度值.计算满足体积约束条件的Lagrangian体积乘子，并通过当前迭代步中的体积乘子和设计变量来更新下一个迭代步中的设计变量.当设计变量满足收敛条件时，优化过程结束.否则，进入下一个循环过程，直到满足收敛条件.对于迭代过程中的收敛性判定，采用相邻迭代的最大密度变化值为判定依据，表达式为

式中，Xw为第w 次迭代的伪密度矢量；ε是很小的密度参考值，一般取0.01.

图3 优化程序流程图Fig.3 Flow of topology optimization

与以目标变化值为收敛判定条件时相比，采用密度变化值判定条件时，所得到的拓扑形态密度分布更为集中，中间密度单元也相对较少，但迭代次数也相应增多.

1.3 数值不稳定处理方法

在拓扑优化过程中，由于场函数不稳定或数值奇异解的存在导致优化后易出现数值不稳定现象，如棋盘格现象、网格依赖性及局部极值等［17］.棋盘格和网格依赖这两种现象一般同时出现在优化结果中，而能够有效去除棋盘格的方法通常也能有效克服网格依赖现象.棋盘格的产生与分析单元的选择有关.研究表明［18－19］，合理选择高阶单元或采用非协调元，可有效降低或消除棋盘格.如图4 所示，8节点的高阶单元比4节点的低阶单元在每个边中心点处多出了1个中间节点，使节点温度函数的一阶导数在单元节点更加圆滑，可有效降低棋盘格现象.

图4 8节点单元Fig.4 8node element

2 数值算例

2.1 平面问题的Dirichlet边界算例

首先，考虑如图5所示的设计问题.设计空间是0.1m×0.1m 的正方形.在内部设计域Ωd存在均匀 生 热 率Q1＝6×104W／m3，在 下 边 界 存 在Dirichlet边界ΓT，边界温度为T0＝0 ℃，边界长度L1＝0.01m，低导热系数λmin＝1 W／（m·K），热传导系数比h＝400，体积约束分别为10%，20%，30%，50%时的传热结构拓扑形态如图6所示（见下页）.设计域内黑色部分表示高导热材料，白色部分表示低导热材料.由拓扑形态可知，高导热材料主要分布在底部温度边界附近，并以枝状形态向四周分歧延伸.随着高导热材料体积约束值的提高，出现了枝的形态逐渐变粗、分歧细化、高导热材料覆盖区域扩大等现象.

图5 均匀生热率载荷Dirichlet边界设计模型Fig.5 Uniform heat generation model satisfying Dirichlet boundary condition

图6 不同体积约束下均匀生热率载荷Dirichlet边界的设计结果Fig.6 Topology forms of uniform heat generation models under different volume fraction

图7 非均匀生热率载荷Dirichlet边界设计模型Fig.7 Non－uniform heat generation design model satisfying Dirichlet boundary condition

表1 生热率载荷Tab.1 Uniform heat generation rate （W·m－3）

其次，考虑在Dirichlet边界下，非均匀生热率载荷的算例模型，如图7所示（见下页）.设计空间是由面积相等的4个子域Ω1，Ω2，Ω3，Ω4组成的正方形区域，每个子域中分别存在均匀生热率载荷Q1，Q2，Q3，Q4，如表1所示（见下页）.通过不同的生热率载荷组合来实现不同的非均匀生热率载荷工况.工况1和工况2生热率载荷为上下分布，工况1的上半部分Ω1和Ω2是高热载荷区域，下半部分Ω3和Ω4是低热载荷区域，工况2与工况1的载荷分布相反.工况3生热率载荷为对角分布，Ω2和Ω3是高载荷区域，Ω1和Ω4是低载荷区域.工况4生热率载荷为对角分布，Ω1和Ω3是高载荷区域，Ω2和Ω4是低载荷区域.体积约束均为30%，低导热系数λmin＝1 W／（m·K），热传导系数比均为400.在不同生热率组合下的传热结构拓扑形态如图8所示.在工况1中，高导热材料主要集中在上半部的高载荷区域和高载荷区与温度边界的最短路径上.在工况2中，高导热材料主要集中在下半部高载荷区域，并且有一部分延伸到了上半部的低载荷区域.在工况3中，高导热材料主要集中在对角形的高载荷区域，在载荷区域高导热材料未被填充.在工况4 中，大部分高导热材料集中在左侧的高载荷区，只有一少部分分支结构延伸到了低载荷区.在4个工况中都明显出现了高导热材料往高热载区域集中的趋势.

图8 相同体积约束下非均匀生热率载荷Dirichlet边界设计结果Fig.8 Topology forms of non－uniform heat generation models under same volume fraction

图9（a）的算例模型是四边Dirichlet边界ΓT，并在内部存在均匀的生热率载荷 Qd＝3×105W／m3.由图9（b）拓扑形态可知，高导热材料从四边向中心区域扩散，并且具有四边对称分布.图10（a）的算例模型是四角Dirichlet边界ΓT，并在内部存在均匀的生热率载荷Qd＝3×105W／m3.由图10（b）拓扑形态可知，高导热材料从四角向中心区域扩散，且在模型中出现了十字型的空白区域.图11（a）的算例模型是四角Dirichlet边界ΓT，并在中心点存在集中生热率载荷Q1＝3×106W／m3.由图11（b）拓扑形态可知，高导热材料由四角和中心点形成了一个X 型的材料分布.图12（a）的算例模型是四边Dirichlet边界ΓT，并在中心点存在集中生热率载荷Q1＝3×106W／m3.由图12（b）拓扑形态可知，高导热材料由四角和中心点形成了一个十字型的材料分布.综上所述，图9（b）和图10（b）均明显出现了有分歧状的拓扑结构，而图11（b）和图12（b）出现了清晰的X和十型的拓扑结构，这是由于集中生热时，热量流动不易受到分支结构的引导，所以，拓扑结构中高导热材料都集中在生热点与热量出口的最短路径上.

图10 四角温度边界－均匀生热Fig.10 Four corners temperature boundary with uniform heat generation

图11 四角温度边界－中心生热Fig.11 Four corners temperature boundary with centre point heat generation

2.2 三维问题的Dirichlet边界算例

图12 四边温度边界－中心生热Fig.12 Four sides temperature boundary with centre point heat generation

对于薄壁管状体，由于管壁相对于管径而言小很多，所以，此类问题可以采用壳单元进行计算.图13是薄壁管状体算例模型，在管的表面有尺寸0.04m×0.04 m 方 形Dirichlet 边 界，边 界 温 度T0＝30 ℃，管体表面有生热率载荷Q1＝10 W／m3，高导热系数λmax＝80 W／（m·K），低导热系数λmin＝0.2 W／（m·K），体积约束为20%.图14（见下页）为薄壁管状体传热结构拓扑形态的不同视角显示，图14（a）是未旋转时的轴测图，图14（b）是以管体中心为旋转点绕X 轴顺时针旋转60°时的视角，图14（c）是旋转180°的视角，可以看出，高导热材料并未出现在方形边界内部，这是因为在边界内部高导热材料的有无并不能改变边界上的温度，所以，边界上的目标敏度在优化过程中始终保持为零，从而高导热材料并未填充到Dirichlet边界内.在拓扑形态上，与平面算例的分布相似，但由于是管体结构，所以，在轴向上并没有实现形态闭合.

图13 三维薄壁管状体结构模型Fig.13 Thin walled tube model

2.3 平面问题的混合边界算例

图14 薄壁管状体结构拓扑结果（不同视角）Fig.14 Different views of the topology form

图15 混合边界模型Fig.15 Model satisfying both Dirichlet and neumann boundary conditions

热传导问题中存在Dirichlet边界和Neumann边界的混合边界的算例模型如图15所示（见下页）.设计空间规格为0.1m×0.1m 的正方形，左下角有宽度为L1的流入热量的Dirichlet边界ΓT，右上角有宽度为L2的流出热量的Neumann边界Γq，4种工况参数如表2所示.体积约束均为30%，低导热系数为λmin＝1 W／（m·K），热传导系数比为400.图16是4种工况参数下的拓扑优化结果.从图16（a）中可以看出，高导热材料集中在两种边界的最短路径上，随着Dirichlet边界长度L1的增加，中间材料向两边分开，如图16（b）所示.当L2的边界长度也由0.01m 增加到0.10m 时，高导热材料越来越多地集中在热量的入口边界上，如图16（c）所示.由此可知，当模型中只有Dirichlet边界和Neumann边界时，高导热材料主要分布在Dirichlet边界和Neumann边界的最短路径上.

表2 工况参数Table.2 Parameters of working condition

图16 混合边界下的优化结果Fig.16 Results of different length of boundary

3 结 论

应用密度法设计传热结构中散热通道的分布，建立了在稳态热传导下的拓扑优化数学模型，推出了满足KKT 条件的传热结构拓扑优化的最优准则法迭代公式，通过对Dirichlet边界和Neumann边界的不同组合下的数值算列，验证了算法的有效性.在拓扑结果的形态上，Dirichlet边界和非集中生热率载荷时，高导热材料分布形态为类似植物分枝结构的拓扑形态；Dirichlet边界和集中生热率载荷时，高导热材料分布在集中载荷点与Dirichlet边界的最短路径上；Dirichlet边界和Neumann边界的混合边界时，高导热材料分布在两种边界的最短路径上，且集中在Neumann边界附近.

［1］Bendsøe M P.Optimal shape design as a material distribution problem［J］.Structural Optimization，1989，4（1）：193－202.

［2］Bendsøe M P，Kikuchi N.Generating optimal topologies instructuraldesign using a homogenizationmethod［J］.Computer Methods in Applied Mechanics and Engineering，1988，71（2）：197－224.

［3］Fujii D，Chen B C，Kikuchi N.Composite material design of two－dimensional structures using the homogenization design method［J］.International Journal for Numerical Methods in Engineering，2001，50（9）：2031－2051.

［4］Wang M，Wang X，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.

［5］Allaire G，Jouve F，Toader A.Structural optimization using sensitivity analysis and a level－set method［J］.Journal of Computational Physics，2004，194（1）：363－393.

［6］Bendsøe M P，Sigmund O.Topology optimization theory，methods and applications［M］.Berlin：Springer Verlag，2003.

［7］Iga A，Nishiwaki S，Izui K，et al.Topology optimization for thermal problems based on assumed continuous approximation of material distributions［J］.Transacions of the Japan Society of Mechanical Engineers，Series C，2007，73（733）：2426－2433.

［8］Yamata T，Nishiwaki S.Level set－based topolog optimization method for thermal problems［J］.Transacions of the Japan Society of Mechanical Engineers，Series C，2009，75（759）：2868－2876.

［9］Iga A，Yamada T，Nishiwaki S，et al.Topogy optimization for coupled thermal and srtuctural problems using the level set method［J］.Transactions of the Japan Society of Mechanical Engineers，Series C，2010，76（761）：36－43.

［10］Bejan A.Constructal－theory network of conducting paths for cooling a heat generating volume［J］.International Journal of Heat and Mass Transfer，1997，40（4）：799－816.

［11］Gersborg－Hansen A，Bendsøe M P，Sigmund O.Topology optimization of heat conduction problems using the finite volume method［J］.Structural and Multidisciplinary Optimization，2004，31（4）：251－259.

［12］Ding X H，Yamazaki K.Stiffener layout design for plate structures by growing and branching tree model（application to vibration－proof design）［J］.Structural and Multidisciplinary Optimization，2004，26（2）：99－110.

［13］Ding X H，Yamazaki K.Constructal design of cooling channel in heat transfer system by utilizing optimality of branch systems in nature［J］.Journal of Heat Transfer，2007，129（3）：245－255.

［14］Bendsøe M，Sigmund O.Material interpolation schemes in topology optimization［J］.Archive of Applied Mechanics，1999，69（9／10）：635－654.

［15］Rozvany G I N.Aims，scope，methods，history and unified terminology of computer－aided topology optimization in structural mechanics［J］.Structural and Multidisciplinary Optimization，2001，21（2）：90－108.

［16］左孔天，陈立平，张云清，等.用拓扑优化方法进行热传导散热体的结构优化设计［J］.机械工程学报，2005，41（4）：13－16.

［17］Sigmund O，Petersson J.Numerical instabilities in topology optimization：a survey on procedures dealing with checkerboards，mesh－dependencies and local minima［J］.Structural Optimization，1998，16（1）：68－75.

［18］Diaz A，Sigmund O.Checkerboard patterns in layout optimization［J］.Structural Optimization，1995，10（1）：40－45.

［19］袁振，吴长春.采用非协调元的连续体拓扑优化设计［J］.力学学报，2003，35（2）：176－180.