基于拓扑优化的新型热沉翅片结构设计

季晓彤　解立垚　王坤　闵春华

文章编号：1007-2373（2022）01-0061-08

摘要 采用基于变密度法的拓扑优化方法对强制对流空气热沉结构进行优化设计，以最小压降作为优化目标，传热性能为约束条件，采用二维双层模型代替传统三维模型进行优化设计。采用在优化过程中改变插值参数方法，有效避免优化结果中阻塞结构的形成。将拓扑结构与直翅片结构进行对比，在入口速度为1.2 m/s时，拓扑结构热沉的平均温度比直翅片热沉降低了约5.4%，努塞尔数增加了约为94.1%，拓扑结构热沉有效地增强了换热效果。研究结果发现拓扑结构热沉翅片在入口和出口处分别存在局部高速区域，能够迅速带走热量，有利于翅片散热。

关 键 词 拓扑优化;二维双层模型;插值公式;传热性能;变密度法

中图分类号 TK124     文献标志码 A

Topology optimization of heat sink structure based on variable density method

JI Xiaotong， XIE Liyao， WANG Kun， MIN Chunhua

（School of Energy and Environmental Engineering， Hebei University of Technology， Tianjin 300401， China）

Abstract The topology optimization method based on the variable density method is used to optimize the structure of the forced convection air heat sink. The minimum pressure drop is used as the optimization goal， the heat transfer performance is the constraint condition， and the two-dimensional double-layer model is used to replace the traditional three-dimensional model for optimization design. This paper adopts the method of changing the interpolation parameters in the optimization process to avoid the blocking structure formation in the optimization results. The comparison of the topological structure with the straight fin structure shows as follow： when the inlet velocity is 1.2 m/s， the average temperature of the topological heat sink is about 5.4% lower than that of the straight fin heat sink， and the Nusselt number increases by about 94.1% ， The topological structure heat sink effectively enhances the heat transfer effect. It is found that the topological structure heat sink fins have local high-speed areas at the entrance and exit respectively， which can quickly take away the heat and facilitate the heat dissipation of the fins.

Key words topology optimization; two-dimensional two-layer model; interpolation formula; heat transfer performance; variable density method

0 引言

隨着微电子技术的发展，电子芯片的热流密度越来越大，从而散热问题变得越来越突出，给热管理带来了很大挑战。传统的热沉结构设计通常是通过热力学计算进行，可是传统方法优化周期过长、自由度较低且无法进行复杂传热结构的优化设计问题。目前新兴的设计方法为拓扑优化设计，此方法具有设计自由度大、灵活性高等特点，因此被广泛应用到热沉结构优化设计中。

结构拓扑优化是以材料分布为优化对象，寻求设计空间内材料分布的最佳方案，得到结构某种性能的最优解。其主要应用包括体积约束下的最小柔度问题、柔性机构拓扑优化设计问题和多物理场拓扑优化设计问题等。在拓扑优化方法中，较为成熟且具有一定工程应用的方法主要有变密度法、均匀化方法及Level Set方法等。近年来这些方法已经应用到传热结构优化设计中。如Yamada等[1]以Level Set方法为基础建立了以热扩散为目标的传热拓扑优化数学模型，在热力耦合情况下考虑热扩散效果，研究了结构刚度最大化问题。Gersborg-Hansen等[2]提出了在有限体积法基础上的传热结构拓扑优化方法。Yoon等[3]在强制对流的换热器结构设计中考虑了传热和压降2个物理场，并且得到流体热传导与对流之间的平衡是结构设计核心问题的结论。裴元帅等[4]利用变密度法，以流体输入功率恒定时的最大化散热为目标进行拓扑优化。发现雷诺数和无量纲化产热系数越大，流道越复杂，散热效果越好;如果出现雷诺数过小，无量纲化产热系数过大则会出现热量富集。游吟等[5]建立了二维相变储热系统模型，对储热系统模型的传热能力进行了研究，讨论了不同肋片设计对于传热能力的影响。李昊等[6]以换热量最大为目标进行不同进出口布置的液冷通道分布优化设计。并以最高温度为评价指标，对比各进口布置的拓扑优化通道与传统直通道的散热性能。安丽[7]进行了流固耦合传热结构拓扑优化设计。尤灏等[8]进行了二维低雷诺数对流强化换热拓扑优化研究。吴璇等[9]基于火积耗散极值原理进行了二维自然对流热沉拓扑优化设计。张晖等[10]根据热源随结构形式变化的热传导拓扑优化问题，提出以温度方差为目标函数的优化设计，并通过数值算例验证了方法的合理性与有效性。乔赫廷等[11]将几何平均温度作为近似最高温度引入到散热结构拓扑优化中，建立了新的散热结构拓扑优化模型，并通过对比最高温度、最大温度梯度等散热性能指标，研究了新优化模型与传统优化模型在热传导拓扑优化设计中的差异及适用范围。赵青海等[12-13]进行了热传导拓扑优化研究，对二维和三维条件下拓扑相关热载荷与周期性多材料结构的拓扑优化设计。魏啸等[14]讨论了不同目标函数对传热结构拓扑优化结果的影响，研究发现以平均温度梯度最小化为目标得出的优化结构具有较好的散热效果。但上述学者或未考虑热对流条件下结论是否依旧适用，又或仅进行了二维条件下拓扑优化研究，未拓展到三维条件下进行研究。

本文的研究内容是基于拓扑优化进行强制对流热沉结构优化设计，并与传统直翅片进行对比。为了降低计算成本，将三维耦合传热模型简化为二维模型。同时，对插值参数进行分段处理，减少可能出现的堵塞结构，使模型更加合理化。最后，基于三维数值模拟，分析了热沉结构的流动与传热特性。

1 模型构建

1.1 双层二维模型

三维耦合传热模型的计算成本较高，而拓扑优化计算过程中需要进行数百次迭代，每次迭代都需要求解耦合传热模型一次，文献[15]表明一般大型计算集群可用于三维耦合传热模型拓扑优化计算。因此，本文基于文献[16]构建了双层二维模型，以减少计算成本。如图1所示，热沉由顶部和底部的实心底座2部分组成，其中顶部再由翅片和通道组成，为了缩短计算所需时间，本文选择截取宽度为10 mm的长条区域进行优化模拟。在二维模型中，假設在高度方向上各部分温度均匀分布，则可以用顶层表示翅片和通道，用底层来表示实心底座。在顶层，存在翅片与流体的对流换热;在底层内部只有热传导，而应用于热沉结构底部表面的恒定热流边界条件则可以用整个底层均匀热源来表示。最后，采用对流换热边界条件将2层结合起来。简化后的双层二维模型如图2所示，其中，左侧为空气入口，右侧为空气出口，为了保证空气能在通道内充分发展，入口和出口处各设有1个30 cm的延长段。中间矩形区域为优化设计部分。整个区域的上下两边设定为对称边界条件。入口速度为1.2 m/s，出口设为压力出口，入口温度为20 ℃。

1.2 控制方程

顶层的控制方程包括连续性方程、动量方程和能量方程，如式（1）～式（3）所示：

式中：[u]为流体的速度，m/s;[ρ]为空气密度，kg/m3;[P]为压力项，Pa;[μ]为空气动力黏度，N·s/m2，[α]为摩擦力系数或称逆渗透率系数;[c1]为空气比热，J/（kg·K）;[∇]为梯度算符。摩擦力系数是通过插值公式与设计变量[γ]进行联系。在顶层中[γ=0]时表示固体区域，此时[α]取最大值，通过式（2）中的摩擦力项保证局部速度为0;当[γ=1]时表示流体区域，此时[α=0]。其中，[k（γ）]和[h（γ）]分别为关于设计变量[γ]的热导率和对流传热系数，[T1]和[T2]分别为顶层温度场和底层温度场，[dz1]为顶层高度。

底层固体域的能量方程为

式中：[ks]为固体热导率，W/（m×K）;[q]为底层施加的热源，W;[dz2]为底层基座高度，m。

2 拓扑优化设计

2.1 插值公式

上述公式中的逆渗透率、热导率和对流传热系数都是设计变量的插值函数。变密度法主要采用SIMP或RAMP材料插值模型，由于RAMP插值模型具有收敛速度快、优化效率高、迭代次数少的优点，因此本文采用RAMP类型函数进行插值[17]，如式（5）～式（7）所示：

式中：[αmin]为最小逆渗透率，本研究中取值为0;[αmax]为最大逆渗透率;[qa]、[qk]和[qh]为插值公式的控制参数，简称为插值参数;[k1]和[k2]分别为固体和流体的导热系数;[h1]为空气与底座之间的对流传热系数;[h2]为翅片与底座之间的对流传热系数。

2.2 密度过滤

为了消除拓扑优化结果中出现的棋盘格和网格依赖性问题，在设计域中采用霍兹海姆偏微分方程进行密度过滤处理，其表达式为

式中：r为长度参数，控制滤波的大小;[γ]为原始设计变量;[γ]为密度过滤后的设计变量，当密度滤波器开启时，所有控制方程中的[γ]将被[γ]取代。

2.3 优化设置

在进行优化设计时，同时考虑了压降及传热性能这两个方面。以前的一些研究[18-20] 多为将两者结合成1个单一的目标函数。然而，2个目标值的量级会有很大差距，且加权比的引入也增加了拓扑优化的复杂性。因此，在本文中，以最小压降为唯一优化目标，传热性能则设定为约束条件。施加传热性能约束的方法为定义1个恒定的底层温度，并根据应用需求设定空气吸收热量的下限。设置底层温度[T2]为314 K，最小吸收热量为2 W。

综上所述，拓扑优化用公式表示为：

最小化：           [ΔP]，

设计变量范围：[0≤γ（r）≤1，]

体积约束：        [Ωγ（r）dΩ≤V，]

最小吸收热量：[Ωh（γ）（T2-T1）dΩ≥q，]

式中：[ΔP]为压力差;[r]为局部坐标;[V]为体积量上限，设为0.6。

2.4 优化计算流程

在前文中进行了几何结构、边界条件、目标和约束条件的定义，接下来分析实现拓扑优化的过程。首先对计算域进行三角单元离散化、有限元分析和求流场和温度场;其次对目标函数、约束条件和灵敏度进行计算，在计算基础上采用全局收敛移动渐近法（Globally Convergent Method of Moving Asymptotes，GCMMA）更新设计变量场。利用调整后的设计变量场重新求解控制方程，得到新的流场和温度场，使其更接近于约束条件或目标函数。再根据灵敏度分析对设计变量场进行二次更新，这个过程反复进行;最后，生成一个满足所有约束条件、最小目标函数的优化结构。本文中拓扑优化过程在有限元软件COMSOL Multiphysics上进行，插值参数的变化则通过MATLAB执行，具体的计算流程如图3所示。其中，theta_0表示设计变量的初始值，本文中设置初始值为0.5。

3 结果分析

3.1 插值参数的动态变化

拓扑优化在计算中将设计域视为多孔介质，优化过程中可能产生非物理结果。例如，强制流体渗透穿过整个设计变量为0的固体区域，这种结构在物理上是不可能的，但是在优化算法中被视为高性能结构。这个阻塞问题在文献[21]中有明确的说明和讨论，结果中也存在可疑的阻塞结构，即通道死角和出口部分之前的固体材料层。为了解决此问题，文献[11]提出了一种在优化过程中逐渐增大插值参数的方法，并且通过多种条件下的拓扑优化计算证明了方法的合理性和有效性。本文采用插值参数值随优化计算先增大后降低的方法，将[qk]和[qh]先增大后降低，以避免阻塞结构的形成;同时逐步增加[αmax]和[qa]的值，增强摩擦阻力项，使局部无限流体的速度趋近于0。

图4为不同插值参数情况下的优化结构图，其中[γ=1]的部分为流体区域，[γ=0]的部分为固体区域。可以看出当[qk]、[qh]增大到500时，基本的拓扑结构已经形成，但是在设计域出口部分开始形成明显的阻塞结构;当[qk]、[qh]减小为80时，能够看出阻塞区域已经明显缩小，但仍有中间值存在;当[qk]、[qh]继续减小至8时，阻塞结构已经基本消除，只有少量中间值存在。

以热导率的插值公式（6）为例对上述结果进行分析。如图5所示，当[qk]的逐渐增大时，导热系数的取值逐渐趋向于两端。可以看出，在[qk=500]的情况下，即便[γ]的值为0.5也能取得与较接近流体状态下的物性参数，但是表现在结构图中则是出现了阻塞结构。因此，逐渐减小[qk]的取值，使[γ]的中间值区域更偏向固体状态下的物性参數，并在此基础上进行优化计算，使优化结构更加合理化。

3.2 模型验证

为了证明双层二维模型可以较为准确的代替三维模型直接进行拓扑优化，即物理有效性。以设计变量等值线0.6提取出翅片的二维形状并进行拉伸和阵列，得到三维热沉结构模型，如图6所示，图7为相同体积比下直翅片结构模型。用COMSOL Multiphysics软件进行优化结构的三维模拟，并与二维优化结果进行对比分析。计算模型设定为稳态不可压缩流动，边界条件的设定与优化问题设置相同，入口速度为1.2 m/s，空气温度设为293.15 K，出口压力设为0 Pa。底部设为恒定热流密度1 666 W/m2，两侧为对称边界条件，其余为绝热面。将三维模拟得出的结果与优化结果比较，具体数据如表1所示。其中，三维热沉基座的平均温度[Tave]与双层模型的底层平均温度[Tave]相差1.5%，双层模型温度场的物理有效性得到了较好证明。同时，二者压差之间相差有1.8%，表明双层模型流场具有较好的物理有效性，可以应用到模拟中。

4 新型结构与传统平直翅片热沉性能比较

以体积占比相同的平直翅片结构热沉为基准，比较了拓扑结构热沉的传热和压降性能。拓扑结构热沉和直翅片热沉的三维模型如图6、图7所示，其中直翅片模型底座部分尺寸为100 mm[×]120 mm[×]17 mm，翅片部分为10 mm[×]108 mm[×]8 mm。运用ANSYS-Fluent软件进行拓扑结构热沉和平直翅片结构热沉三维数值模拟，除入口速度不同外，外其他边界条件的设定与3.2节相同。

图8为不同入口速度下，拓扑结构热沉和直翅片热沉基座平均温度对比图。由图可知，随着入口速度逐渐增加，流体迅速流过翅片带走更多热量，使得两种热沉的基座平均温度逐渐降低。并且随着速度逐渐升高，温度的降低趋势逐渐变缓。同时，2种结构平均温度之间的差值随着入口速度的增加逐渐缩小，当进口速度为0.5 m/s时，拓扑结构热沉的平均温度为344.2 K，直翅片热沉的平均温度为359.9 K，两者之间相差约为4.3%;当进口速度为1.2 m/s时，拓扑结构热沉的平均温度为321.2K，直翅片热沉的平均温度为339.8 K，两者之间相差约为5.4%;当进口速度为2 m/s时，拓扑结构热沉的平均温度为312.2 K，直翅片的平均温度为326.5 K，两者之间相差降低至4.3%。

图9是入口速度为1.2 m/s条件下，通道高度为4 mm部分的拓扑结构热沉和直翅片热沉的温度场与速度场。直翅片热沉的温度场如图9a）所示，翅片温度分布均匀，平均温度约为330 K。流体温度边界层逐渐变厚，使得对流换热系数减小，不利于直翅片热沉的散热。图9b）为优化结构热沉温度云图，流体温度边界层在入口微通道处较小，说明此处换热系数较大，能够进行有效对流换热，中间通道部分流体温度逐渐升高。拓扑结构热沉的3个翅片温度均大约为318 K，远低于直翅片结构的平均温度。

图9c）为直翅片热沉的速度云图，其翅片区域流体速度分布均匀，平均速度约为1.8 m/s。图9d）为拓扑结构热沉的速度云图，可以看出前段部分翅片的不规则结构有效阻碍了速度边界层的生长，使得前段流体速度保持较高的状态，同时微通道处流速进一步提高能够快速带走热量;后段翅片使空气进入2个分支通道，并且在下侧通道处形成了局部高速区域，这有助于提高局部传热系数。

在入口温度为293 K的工况下，不同入口速度对拓扑结构热沉和直翅片热沉的影响，如图10所示。图10a）中横坐标为热沉底座平均温度，纵坐标为进出口压力差。从图中可以看出，进出口压差随着平均温度逐渐增加而迅速降低，同时当平均温度低于300 K或平均温度高于320 K时，两者压差之间的差值较小;在平均温度为300 K和320 K时，两者之间的差值开始逐渐增大，在平均温度为310 K时，两者之间的差值达到最大约为56%，随后逐渐降低。图10b）为拓扑结构热沉与直翅片热沉努塞尔数随着入口速度的变化图。从图中可以看出，两种热沉的努塞尔数均随着入口速度增大逐渐增大，这是由于入口速度的增加导致空气带走更多的热量，换热效率增强。同时拓扑结构热沉的努塞尔数始终高于直翅片热沉，当进口速度分别为0.5 m/s、1.2 m/s和2 m/s时，拓扑结构的努塞尔数较直翅片热沉的分别提高了75.5%、94.3%和59.6%。在拓扑结构热沉中，流体入口段和出口段均有局部高速区域，提高了局部换热系数，增强了流体和翅片间的对流换热。

5 结语

1）为了减少拓扑优化的计算成本且让拓扑模型更加贴合实际情况，本文采用了双层二维模型进行优化设计，顶层代表流体域和翅片结构，底层代表了热沉结构的基座。并且，运用模拟软件对得出的拓扑结构进行模拟验证，证明了双层二维模型的可行性。

2）通过运用动态改变插值参数的方法，有效防止了阻塞结构的产生。使得到的优化结构更加符合真实物理结构，便于加工。

3）将优化后得到的拓扑结构与相同体积比的直翅片结构进行传热和压降性能对比，发现拓扑结构热沉翅片在入口和出口处分别存在局部高速区域，能够迅速带走热量，有利于翅片散热。分析了不同入口速度下2种热沉的基座平均温度和努赛尔数。在入口速度为1.2 m/s时，拓扑结构热沉的平均温度比直翅片热沉降低了约5.4%，努塞尔数增加了约为94.3%，拓扑结构热沉有效地增强了换热效果。

参考文献：

[1]    YAMADA T，NISHIWAKI S，IGA A，et al. Level set-based topology optimization method for thermal problems（mechanical systems）[J]. Transactions of the Japan Society of Mechanical Engineers Series C，2009，75（759）：2868-2876.

[2]    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，2006，31（4）：251-259.

[3]    YOON G H. Topological design of heat dissipating structure with forced convective heat transfer[J]. Journal of Mechanical Science and Technology，2010，24（6）：1225-1233.

[4]    裴元帅，王定标，王光辉，等. 恒定压降下电子元件热沉的拓扑优化设计[J]. 低温与超导，2020，48（1）：56-61.

[5]    游吟，赵耀，赵长颖，等. 相变储热单元内肋片结构的拓扑优化[J]. 科学通报，2019，64（11）：1191-1199.

[6]    李昊，丁晓红，景大雷. 液冷通道分布优化设计的仿真和试验研究[J]. 机械工程学报，2019，55（10）：198-206.

[7]    安丽. 流固共轭传热结构拓扑优化[D]. 大连：大连理工大学，2020.

[8]    尤灏，陈熹挺，王文. 风冷散热翅片结构的拓扑优化[J]. 制冷技术，2020，40（3）：6.

[9]    吴璇，陈群. 基于不动点迭代的自然对流热沉拓扑优化[J]. 工程热物理学报，2020，41（9）：6.

[10]  张晖，刘书田，张雄. 拓扑相关热载荷作用下稳态热传导结构拓扑优化[J]. 中国机械工程，2009，20（11）：1339-1343.

[11]  乔赫廷，张永存，刘书田. 散热结构拓扑优化目标函数的讨论[J]. 中国机械工程，2011，22（9）：1112-1117，1122.

[12]  赵清海，张洪信，蒋荣超，等. 考虑拓扑相关热载荷的散热结构多相材料拓扑优化设计[J]. 中国机械工程，2020，31（20）：9.

[13]  赵清海，张洪信，华青松，等. 周期性多材料结构稳态热传导拓扑优化设计[J]. 工程力学，2019，36（3）：10.

[14]  魏啸，丁晓红. 不同目標函数的传热结构拓扑优化研究[J]. 电子科技，2017，30（2）：5.

[15]  ALEXANDERSEN J，SIGMUND O，AAGE N. Large scale three-dimensional topology optimisation of heat sinks cooled by natural convection[J]. International Journal of Heat and Mass Transfer，2016，100：876-891.

[16]  ZENG S，KANARGI B，LEE P S. Experimental and numerical investigation of a mini channel forced air heat sink designed by topology optimization[J]. International Journal of Heat and Mass Transfer，2018，121：663-679.

[17]  陈祥，刘辛军. 基于RAMP插值模型结合导重法求解拓扑优化问题[J]. 机械工程学报，2012，48（1）：135-140.

[18]  KOGA A A，LOPES E C C，VILLA NOVA H F，et al. Development of heat sink device by using topology optimization[J]. International Journal of Heat and Mass Transfer，2013，64：759-772.

[19]  QIAN X P，DEDE E M. Topology optimization of a coupled thermal-fluid system under a tangential thermal gradient constraint[J]. Structural and Multidisciplinary Optimization，2016，54（3）：531-551.

[20]  Ozguc S，Liang P，Weibel J A. Topology optimization of microchannel heat sinks using a homogenization approach[J]. International Journal of Heat and Mass Transfer，2021，169（1/2）：120896.

[21]  VAN OEVELEN T，BAELMANS M. Numerical topology optimization of heat sinks[C]// Proceedings of the 15th International Heat Transfer Conference. August 10-15，2014. Kyoto，Japan. Connecticut：Begellhouse，2014.

收稿日期：2021-03-15

基金项目：河北省青年科学基金（E2019202255）

第一作者：季晓彤（1995—），女，硕士研究生。通信作者：解立垚 （1986—），男，讲师，xieliyao@hebut.edu.cn。