树状分支传热结构层次生长优化设计技术

随着现代工业及科学技术的飞速发展，电子元器件日益向高集成度、高性能和小型化发展，导致其热流密度大幅增加，例如，大型计算机CPU芯片在运行过程中产生的热流密度已达到200 W/cm2以上[1]。随着热流密度和温度的上升，电子元器件使用可靠性和使用寿命都会急剧下降。电子元器件的正常工作温度一般低于70℃，据统计，芯片温度每升高2℃，其性能可靠性会降低10%[2-3]；高达55%的电子元器件失效是温度超过规定值引起的[4]。目前常用的散热技术种类繁多，传统的空气冷却方式的极限为50 W/cm2，显然已经不能满足高热流密度电子元件的散热需求，同时也受制于设计空间的约束；而液冷流道和相变热管等散热技术需要较大的空间才能实现。为了将小空间高热流密度电子元器件内部的高发热量快速有效地导出，需要通过嵌入高导热材料形成高效的散热通道进行导热散热，即通过将高导热材料形成的散热通道敷设于电子元器件表面或直接嵌入元器件内部，将热量快速导出至外界环境。

为了得到高效的散热通道，设计出最优的散热通道分布/拓扑形态，近年来众多的国内外学者和学术团体开展了广泛的研究。1997年BEJAN[5]提出构形理论(constructal theory)，可以用来设计散热通道。构形理论以最小热阻为目标进行设计，但是设计中须假设高导热材料的热导率远大于低导热材料的热导率，因此在低热导率比和填充率较高的情况下不适用[6]；同时构形理论存在最小热阻极限即装配极限，导致设计结果无法取得最优的传热效果[7]。此后随着拓扑优化技术的发展，许多学者采用结构拓扑优化技术设计散热通道，主要有SIMP(solid isotropic material with penalization)方法[8-10]、ESO(evolutionary structural optimization)方法[11-12]和水平集(level set)方法[13]。上述方法可突破构形理论的局限，但是仍存在灰度单元、棋盘格现象等问题，且最终的散热通道拓扑形态中会存在无效的细小分支，从而导致制造加工困难。另一方面，自然界中的动植物经过亿万年的进化发展，形成了各种优异高效的分支结构以满足自身生存需求。由此，DING等[6]基于自然界植物根系形态生成机理，提出了一种自适应成长法，构建了高效散热通道，然而，该方法是单纯从生物分支网形态构筑机理出发，没有考虑散热系统热设计的要求。

本文基于植物叶片有效的水(养)分输送、能量交换和散热冷却机制与高热流密度条件下树状分支散热通道的热传输功能具有本质的相似性这一特性，研究植物网状叶脉多层次结构特征，形成多层次树状分支散热通道生成方法，并通过不同热设计条件下二维和三维传热结构的设计结果，验证此设计技术的有效性。

1 树状分支散热通道层次生长策略

1.1 植物叶脉分层结构特征

叶脉分支拓扑构型是水分运输的重要通道，叶脉宏观的分支形态和微观的叶片组织使得植物叶片进化出了十分有效的散热冷却机制来应付高温[14]，避免自身被过度光照晒伤。植物叶脉的分支拓扑构型有多种，如分叉状脉、网状脉和平行脉等，其形成机理和脉序等级有关系，而高级脉序与低级脉序的形成影响因素及生物功能各不相同。细观网状叶脉(图1)，其叶脉系统主要由主脉和次脉组成，主脉是低级脉序，负责叶片内外部的能量和养分输送与交换；次脉为高级脉序，有序地和主脉连通，为主脉收集能量和养分。低级脉序和高级脉序由于其功能不同，具有分层次生长的特点，从而形成了不同的分布形态。

1.2 散热通道分层次生长方法

在小空间内，高热流密度条件下，通过布置一定体积的高导热材料形成散热通道将热量有效导出，使得设计结构的散热性能最优。如图2所示，设计域Ωd内存在均匀或非均匀生热热源，热量通过散热边界散出，其余边界绝热。设计问题中散热通道的生成控制策略将决定散热通道的分布形式，以及散热结构的散热性能。

散热结构中热的传输本质上讲是一种“流”，与自然界植物叶脉中水分输送一样都属于“物质流”的输送过程，因此散热通道的生长可依据叶脉脉序的功能不同导致形态不同的原理，采用类似的生成方式，主次通道依次分层生长。

如图3所示，定义散热通道由不同尺寸的矩形组成，采用起点、终点和通道宽度3个参数即可确定任一通道的位置和尺寸。

1.2.1 主通道的生长

主通道是传导热量的主要结构，其主要作用是负责将热量由结构内部直接有效地导出至热沉点，基于这一思想，主通道以热沉作为起点，采用考虑温度场的空间殖民算法计算通道终点。空间殖民算法是一种分割空间的算法，RUNIONS等[15-16]基于空间竞争的思想提出空间殖民算法，并将其用于树形结构的可视化建模以及叶脉的可视化建模。LOHAN等[17]利用空间殖民算法对热传导问题中传热结构进行拓扑优化设计。本文在此基础上，考虑结构温度场，结合高温点的温度和位置信息，采用考虑温度信息的空间殖民算法(space colonization algorithm with temperature information，SCATI)搜寻主通道终点。

如图4a所示，首先找出结构内最高温度点s1，并找出与该点距离小于或等于d的点s2，s3，…，sn，得到集合S(v)；以热沉点位置p为起点，集合S(v)内点的位置和温度值形成各点矢量，如图4b虚线箭头所示；由下式可得矢量和n(如图4c实线箭头所示)：

式中，Ti为S(v)中各个点si的温度值；

为S(v)中各个点si的温度平均值。

按照下式确定主通道终点(即图4d所示p′位置点)：

根据叶脉层序特征的研究，叶脉不同层序体积分数与叶脉面积之间具有一定的统计规律[18]，在此，主通道体积分数按照叶脉中主脉体积分数统计规律确定。

1.2.2 次通道的生长

次通道的功能是收集热量，将设计域内的热量快速有效地传导至主通道，因此选择温度最高点作为次通道的终点，而起点的选择必须能使得热传导效率尽可能高。从能量学角度出发考虑，遵循类似于自然分支系统中普遍存在的最小能耗原则，故采用最小热阻原理来进行次通道的生长。

一根分支通道的热阻r可以定义为

其中，L，h和D分别为分支通道的长度、高度和宽度；κ为材料的热导率。从热沉点到通道末端，包含n根分支通道的总热阻R可定义为

其中，ri为第i个分支通道的热阻；n为导热通道中分支通道的个数。如图5所示，新生长的通道为子通道，与其相连已存在的通道为其母通道。沿母通道建立局部坐标系X′Y′如图5所示，从而总热阻R采用局部坐标系可表示为

其中，xs为子通道起点在X′轴上的坐标值；xe和ye分别为子通道终点在局部坐标系中的X′和Y′坐标值；Dn和Dn-1分别为子通道和与其相连的母通道的宽度；R′为传热路径中除新子通道之外的通道热阻。为求最小热阻，对热阻R表达式求导数，获得极值点：

由此可确定次通道的起点和终点。给定子通道初始枝宽D0，完成一根次通道的生长。由图5可知，原母通道被子通道起点分割成两段：

式中,λ为分歧指数，根据能量损耗最小原理，λ取值为3[19]。

1.2.3 分层次生长流程图

图6为传热结构分层次生长流程图。首先设定设计域以及热设计条件和参数，包括设计域尺寸、热沉条件、热源条件、高导热材料和低导热材料的热导率、主通道体积分数上限、通道初始宽度、所有通道体积分数上限等。然后利用有限元热分析，得到设计域温度场，采用SCATI算法实现主通道生长，直至满足主通道体积分数要求，主通道生长完成。主通道生长完成后，开始次通道的生长。次通道的生长依据最小热阻原理确定生长分歧点，按照Murray法则更新枝宽。最后，当高导热材料的体积分数达到设计要求时生长迭代停止，即传热通道生长完成。

2 典型数值算例

以热通道分层生长的方法，对典型的二维均匀生热及非均匀生热和三维热传导算例进行散热通道优化设计。

2.1 单点热沉均匀生热算例

图7所示为典型均匀生热单点热沉问题。设计域为0.1 m×0.1 m的正方形，其内部存在均匀生热，生热率Q=3 000 W/m2。在设计域的底边中部有热沉作为散热边界，其长度b=0.01 m，热沉温度t0=0 ℃，其余边界绝热。高导热材料与低导热材料的热导率比K=400。根据相应生物学叶脉中主脉体积的统计规律，主通道体积分数上限为16.93%[18]，高导热材料体积分数上限为20%。

设计域内黑色部分表示由高导热材料敷设的散热通道，白色部分表示低导热材料。如图8a所示，主通道较为均匀地分布在整个设计域内，将设计域分成多个区域；主通道的长度较长，从热沉出发延展至设计域的边界处。次通道主要出现在设计域内主通道分布稀疏处，如图8b所示，可以发现，随着次通道的生长，大部分主通道及个别次通道的宽度会随之更新、增宽，其中有分支的通道宽度与无分支的通道宽度相比更宽。比较次通道与主通道的长度可以发现，次通道的长度要远小于主通道。可以发现，分支现象主要发生在主通道上。当高导热材料的体积分数到达设计上限时生长过程停止，获得最终拓扑形态。计算最终结果，主通道体积分数为16.98%，与计算要求的主通道体积分数16.93%十分接近。设计域温度场如图8c所示。

为了验证所设计散热结构的传热性能，同条件下采用SIMP方法设计散热结构拓扑形态进行对比，如图8d所示。两种方法所设计传热通道形态类似，其中SIMP方法设计结果与天然树形结果更相似，但其热通道细小、不规则，存在灰度单元，不利于传统加工方式制造。分层生长方式设计的热通道边界清晰、无灰度单元，可采用线切割等低成本加工方式进行制造。两种设计结果的热性能对比见表1，与SIMP方法所设计结构相比，分层生长方法所设计结构的最高温度、平均温度和温度方差分别下降了28.9%、6.8%和55.4%。

2.2 两点热沉非均匀生热算例

图9所示为两点热沉非均匀生热模型，热沉的温度和尺寸保持不变，数量增加1，相应的设置方式如图8所示；非均匀生热率Q1=6 000 W/m2，Q2=600 W/m2，布置方式如图9所示，其余参数与边界条件均保持不变。

图10a所示为主通道生长结束后散热通道拓扑形态，主通道主要分布在高热流密度区，与热沉点连接。图10b所示次通道主要分布在高热流密度且主通道分布较为稀疏的区域。散热通道最终的拓扑形态及其温度场分布图见图10c。

同条件下采用SIMP方法设计散热结构拓扑形态，如图10d所示，其热性能对比见表2。在非均匀热源两点热沉条件下，与SIMP方法所设计结构相比，分层生长方法所设计结构的最高温度、平均温度和温度方差分别下降了42.1%、8.2%和73.9%。分层生长方法所设计的散热结构结构简单且拓扑边界清晰，不会出现SIMP方法中存在的灰度单元或边界模糊的现象。

2.3 四点热沉非均匀生热算例

图11所示为四点热沉非均匀生热模型，热沉的温度和尺寸保持不变，数量增加为4，分别位于设计域的4个角点；非均匀生热率Q1=6 000 W/m2，Q2=600 W/m2，布置方式如图11所示，其余参数与边界条件均保持不变。

图12a所示为主通道生长结束后散热通道拓扑形态，主通道主要分布在高热流密度区，与热沉点连接。次通道主要分布在高热流密度区且主通道分布较为稀疏的区域。图12b所示为散热通道最终的拓扑形态，高热流密度区出现较多细小次通道分支。散热结构温度场分布如图12c所示。

同条件下采用SIMP方法设计散热结构拓扑形态，如图12d所示，其热性能对比见表3。在非均匀热源四点热沉条件下，与SIMP方法所设计结构相比，分层生长方法所设计结构的最高温度、平均温度和温度方差分别下降了49.0%、19.7%和64.7%。同样地，分层生长方法所设计的散热结构结构简单且拓扑边界清晰，不会出现SIMP方法中存在的灰度单元或边界模糊的现象。

2.4 三维传热结构算例

图13所示为三维传热结构算例。设计域为半圆柱曲面，圆柱半径r=0.1 m，圆柱长度为0.2 m，曲面均匀生热，生热率Q=3×103 W/m2。如图13所示，在设计域的底边中部有一个热沉作为散热边界，其长度b=0.01 m，热沉温度t0=0 ℃，其余边界绝热。

图14a和图14b所示为分层生长方法所设计三维半圆柱设计域内散热通道拓扑形态，其分布与单点热沉平面结构类似。主通道较为均匀地分布在整个设计域内,次通道主要分布在主通道分布较为稀疏的区域；散热结构温度场分布如图14c所示。设计域内最高温度、平均温度和温度方差分别为1.026 ℃、0.790 ℃和0.050 0。

3 结论

本文针对高热流密度电子元件所需小空间散热方式，提出了一种分层生长散热通道拓扑优化设计方法。与SIMP的算例比较，所提出的设计方法设计结果轮廓清晰、无灰度单元，易于采用线切割等传统低成本加工方法进行加工，并且具有较好传热性能。通过单点热沉、两点热沉、四角热沉、均匀生热、非均匀生热以及三维设计域等多种设计算例，验证了所提方法的适用性。

[1] 郝俊娇，潘日，周刚，等. 高热流密度电子元件中热管散热技术的进展[J]. 化工进展，2015，34(5)：1220-1231.

HAO Junjiao，PAN Ri，ZHOU Gang，et al. Development of Heat Pipe Cooling Technology in High Heat Flux Electronic Components[J]. Chemical Industry and Engineering Progress，2015，34(5)：1220-1231.

[2] CHOI J，JEONG M，YOO J，et al. A New CPU Cooler Design Based on an Active Cooling Heatsink Combined with Heat Pipes[J]. Applied Thermal Engineering，2012，44(2)：50-56.

[3] 诸凯，王彬，王雅博，等. 用于芯片冷却的重力式热管散热器实验研究[J]. 工程热物理学报，2017，38(8)：1748-1752.

ZHU Kai，WANG Bin，WANG Yabo，et al. Experimental Study of Gravity Heat Pipe Heat Sink Used for Chip Cooling[J]. Journal of Engineering Thermophysics，2017，38(8)：1748-1752.

[4] KANDA R，JAIN K. Thermal Management of Low Volume Complex Electronic Systems[C]∥ Proceedings of the International Conference on Research and Innovations in Mechanical Engineering. Ludhiana，2013：519-533.

[5] 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.

[6] 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.

[7] ZHANG Y，LIU S. Design of Conducting Paths Based on Topology Optimization[J]. Heat and Mass Transfer，2008，44(10)：1217-1227.

[8] ZHAO X，ZHOU M D，SIGMUND O，et al. A “Poor Man’s Approach” to Topology Optimization of Cooling Channels Based on a Darcy Flow Model[J]. International Journal of Heat and Mass Transfer，2018，116：1108-1123.

[9] QIAN S H，WANG W，GE C L，et al. Topology Optimization of Fluid Flow Channel in Cold Plate for Active Phased Array Antenna[J]. Structural and Multidisciplinary Optimization，2018，57：2223-2232.

[10] 崔天福，丁晓红，侯丽园. 基于密度法的传热结构拓扑优化[J]. 上海理工大学学报，2014，36(6)：548-555.

CUI Tianfu，DING Xiaohong，HOU Liyuan. Structural Topology Optimization Design of Heat Transfer Based on Density Method[J]. Journal of University of Shanghai for Science and Technology，2014，36(6)：548-555.

[11] XIE Y M，STEVEN G P. A Simple Evolutionary Procedure for Structural Optimization[J]. Computers and Structures，1993，49(5)：885-896.

[12] 张林森，丁晓红，熊敏. 基于渐进结构法的散热通道构建技术[J]. 通信电源技术，2015，32(2)：13-16.

ZHANG Linsen，DING Xiaohong，XIONG Min. Construction Technology of Cooling Channel Based on Evolutionary Structural Optimization Method[J]. Telecom Power Technology，2015，32(2)：13-16.

[13] SETHIAN J A，WIEGMANN A. Structural Boundary Design via Level Set and Immersed Interface Methods[J]. Journal of Computational Physics，2010，163(2)：489-528.

[14] PENG Y , LIU W , CHEN W , et al. A Conceptual Structure for Heat Transfer Imitating the Transporting Principle of Plant Leaf[J]. International Journal of Heat and Mass Transfer, 2014, 71:79-90.

[15] ROUNIONS A，LANE B，PRUSINKIEWICZ P. Modeling Trees with a Space Colonization Algorithm Modeling Trees with a Space Colonization Algorithm[R]. Prague : Eurographics Workshop on Natural Phenomena，2007.

[16] ROUNIONS A，FUHRER M，LANE B，et al. Modeling and Visualization of Leaf Venation Patterns[J]. ACM Transactions on Graphics，2005，24 (3)：702-711.

[17] LOHAN D J，DEDE E M，ALLISON J T. Topology Optimization for Heat Conduction Using Generative Design Algorithms[J]. Structural and Multidisciplinary Optimization，2017，55(3)：1063-1077.

[18] SACK L，SCOFFONI C，MCKOWN A D，et al. Developmentally Based Scaling of Leaf Venation Architecture Explains Global Ecological Patterns[J]. Nature Communications，2012，3(1)：810-837.

[19] MURRAY C D. The Physiological Principle of Minimum Work. I. the Vascular System and the Cost of Blood Volume[J]. Proceedings of the National Academy of Sciences of the United States of America，1926，12(3)：207-214.